fix various bug in reload

hxt/hxt_linear_system_petsc.c

@@ -68,7 +68,6 @@ static HXTStatus assembleIfNeeded(HXTLinearSystemPETSc *system)
     system->assemblyNeeded = 0;
   }
   if (system->nIdentRows > 0){
-    printf("%i identity rows\n", system->nIdentRows);
     HXT_PETSC_CHECK(MatZeroRows(system->a,
           system->nIdentRows, system->identRows,1., NULL, NULL));
   }

src/fluid_problem.c

@@ -782,6 +782,8 @@ void fluid_problem_after_import(FluidProblem *problem)
 {
   mesh_tree_free(problem->mesh_tree);
   problem->mesh_tree = mesh_tree_create(problem->mesh);
+  free(problem->old_porosity);
+  problem->old_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
   for (int i = 0; i < problem->n_particles; ++i)
     problem->particle_element_id[i] = -1;
   mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);

src/fluid_problem_io.c

@@ -54,11 +54,11 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
       mb_vtk_data_array_destroy(&a);
     }
     else if (!strcasecmp(e.name,"Cells")) {
-      MBVtkDataArray connectivity;
       int ncon = -1;
       while(mb_read_vtk_data_array(r,&a) != 0){
         if(!strcasecmp(a.name,"Connectivity")) {
-          connectivity = a;
+          m->elements = mb_vtk_data_array_to_int(&a,m->n_elements,DIMENSION+1);
+          mb_vtk_data_array_destroy(&a);
         }
         else if(!strcasecmp(a.name,"Offsets")){
           int *offsets = mb_vtk_data_array_to_int(&a,m->n_elements,1);
@@ -81,18 +81,20 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
       problem->solution = malloc(sizeof(double)*(DIMENSION+1)*m->n_nodes);
       while(mb_read_vtk_data_array(r,&a) != 0){
         if(!strcasecmp(a.name,"Porosity")) {
-          if(problem->porosity)
-            free(problem->porosity);
+          free(problem->porosity);
           problem->porosity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
         }
         else if(!strcasecmp(a.name,"Pressure")) {
           double *p = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
-          memcpy(problem->solution+DIMENSION*m->n_nodes,p,sizeof(double)*m->n_nodes);
+          for (int i = 0; i < m->n_nodes; ++i)
+            problem->solution[(DIMENSION+1)*i+DIMENSION] = p[i];
           free(p);
         }
         else if(!strcasecmp(a.name,"Velocity")) {
           double *v = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
-          memcpy(problem->solution,v,sizeof(double)*m->n_nodes*DIMENSION);
+          for (int i = 0; i < m->n_nodes; ++i)
+            for (int d = 0; d < DIMENSION; ++d)
+              problem->solution[(DIMENSION+1)*d+DIMENSION] = v[i*DIMENSION+d];
           free(v);
         }
         mb_vtk_data_array_destroy(&a);
@@ -109,13 +111,14 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
   mb_xml_read_element(r,&fielddataelement);
   int n_phys=0;
   int n_phys_alloc = 10;
-  MBVtkDataArray *phys = malloc(sizeof(MBVtkDataArray*)*10);
+  MBVtkDataArray *phys = malloc(sizeof(MBVtkDataArray)*10);
   while(mb_read_vtk_data_array(r,&a) != 0){
     if(!strncasecmp(a.name,"Physical ",9)) {
       if (n_phys == n_phys_alloc) {
         n_phys_alloc *= 2;
-        phys = realloc(phys,n_phys_alloc*sizeof(MBVtkDataArray*));
+        phys = realloc(phys,n_phys_alloc*sizeof(MBVtkDataArray));
       }
+      phys[n_phys++] = a;
     }
     else {
       printf("unknown data array \"%s\" in vtk piece\n", a.name);
@@ -130,12 +133,14 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
   m->phys_dimension=realloc(m->phys_dimension,n_phys*sizeof(int));
   m->phys_n_nodes=realloc(m->phys_n_nodes,n_phys*sizeof(int));
   m->phys_name = realloc(m->phys_name,n_phys*sizeof(char*));
-  m->phys_nodes = realloc(m->phys_name,n_phys*sizeof(int*));
+  m->phys_nodes = realloc(m->phys_nodes,n_phys*sizeof(int*));
   for(int i = 0; i < n_phys; ++i) {
-    int p;
+    MBVtkDataArray a = phys[i];
+    int p, n;
     sscanf(a.name+9,"%i %i%n",&m->phys_dimension[i],&m->phys_id[i],&p);
     m->phys_nodes[i] = mb_vtk_data_array_to_int(&phys[i],-1,-1);
     m->phys_n_nodes[i] = phys[i].number_of_tuples;
+    m->phys_name[i] = strdup(a.name+p+10);
     mb_vtk_data_array_destroy(&phys[i]);
   }
   free(phys);
@@ -215,7 +220,7 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
   fprintf(f,"</Piece>\n");
   fprintf(f,"<FieldData>\n");
   for(int i = 0; i < mesh->n_phys; ++i) {
-    fprintf(f,"<DataArray Name=\"Physical %i %i %s\" NumberOfTuples=\"%i\" NomberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n",mesh->phys_dimension[i],mesh->phys_id[i],mesh->phys_name[i],mesh->phys_n_nodes[i]);
+    fprintf(f,"<DataArray Name=\"Physical %i %i %s\" NumberOfTuples=\"%i\" NumberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n",mesh->phys_dimension[i],mesh->phys_id[i],mesh->phys_name[i],mesh->phys_n_nodes[i]);
     for (int j = 0; j < mesh->phys_n_nodes[i]; ++j) {
       fprintf(f, "%i ", mesh->phys_nodes[i][j]);
     }

src/mbxml.c

@@ -145,9 +145,11 @@ int mb_read_vtk_data_array(MbXmlReader *reader, MBVtkDataArray *a) {
       a->name = strdup(e.attr[i].value);
     else if(!strcasecmp(e.attr[i].name,"NumberOfTuples"))
       sscanf(e.attr[i].value,"%i",&a->number_of_tuples);
-    else if(!strcasecmp(e.attr[i].name,"NumberOfComponents")){
+    else if(!strcasecmp(e.attr[i].name,"NumberOfComponents"))
+      sscanf(e.attr[i].value,"%i",&a->number_of_components);
+    //temporary fix
+    else if(!strcasecmp(e.attr[i].name,"NomberOfComponents"))
       sscanf(e.attr[i].value,"%i",&a->number_of_components);
-    }
   }
   a->content = strdup(mb_xml_read_content(reader));
   mb_xml_end_element(reader,&e);
@@ -188,6 +190,7 @@ static void validate_array(MBVtkDataArray *a, int number_of_tuples, int number_o
   }
   if (a->number_of_tuples == -1) {
     printf("number of tuples not specified for vtk array \"%s\"\n", a->name);
+    exit(1);
   }
   if (number_of_components !=-1) {
     if(number_of_components != a->number_of_components && a->number_of_components != -1){
@@ -197,7 +200,8 @@ static void validate_array(MBVtkDataArray *a, int number_of_tuples, int number_o
     a->number_of_components = number_of_components;
   }
   if (a->number_of_components == -1) {
-    printf("number of components not specified for vtk array \"%s\"", a->name);
+    printf("number of components not specified for vtk array \"%s\"\n", a->name);
+    exit(1);
   }
 }
 
@@ -225,7 +229,7 @@ int *mb_vtk_data_array_to_int(MBVtkDataArray *a, int number_of_tuples, int numbe
   FILE *f = fmemopen(a->content, strlen(a->content),"r");
   for(size_t i = 0; i < size; ++i){
     if(fscanf(f,"%d",&v[i]) != 1) {
-      printf("cannot parse double in vtk array \"%s\"", a->name);
+      printf("cannot parse int in vtk array \"%s\"\n", a->name);
       exit(1);
     }
   }

testcases/drop-3d/Metzger/dropReload.py

@@ -9,7 +9,7 @@ import shutil
 import random
 
 
-outputdir = "MetzgerB1"
+outputdir = "MetzgerB1Reload"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
@@ -46,11 +46,11 @@ mu = nu*rho
 
 
 p = scontact3.ParticleProblem()
-p.read_vtk(outputdir,2)
+p.read_vtk("MetzgerB1",2)
 
 print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
 print("RHOP = %g" % rhop)
-outf = 10
+outf = 1
 outf1 = 100000
 
 #strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
@@ -58,7 +58,7 @@ outf1 = 100000
 strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
 fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
 print("Import start!!!")
-fluid.import_vtk(outputdir+"/fluid_00002.vtu")
+fluid.import_vtk("MetzgerB1"+"/fluid_00002.vtu")
 print("Import done!!!")
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())