correction mesh.c

fluid/mesh.c

@@ -182,6 +182,7 @@ Mesh *mesh_load_msh(const char *filename)
   }
   for (int i = 0; i < n_enodes; ++i) {
     m->elements[i] = get_node_id(m->elements[i], nodes_map, m->n_nodes);
+    printf("m->elements[%d] = %d\n",i,m->elements[i]);
   }
 
 
@@ -226,10 +227,11 @@ Mesh *mesh_load_msh(const char *filename)
     for (int itag = 0; itag < e->n_physicals; ++itag) {
       for (int ip = 0; ip < gmsh_m->n_physical_names; ++ip) {
         if (gmsh_m->physical_name_dim[ip] == e->dim && gmsh_m->physical_name_tag[ip] == e->physicals[itag]) {
-          for (int iel = 0; iel< e->n_elements; ++iel) {
+          for (int iel = 0; iel < e->n_elements; ++iel) {
             int keep = 1;
             for (int inode = 0; inode < DIMENSION; ++inode) {
-              keep &= (get_node_id(m->elements[iel*DIMENSION+inode], nodes_map, m->n_nodes) >= 0);
+              int nid = get_node_id(e->elements[iel*DIMENSION+inode], nodes_map, m->n_nodes);
+              keep &= (nid >= 0);
             }
             if (keep)
               n_bnd_el +=1;

testcases/periodicTest/grains.py

@@ -91,7 +91,7 @@ outf = 1                                       # number of iterations between ou
 
 #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
 fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
-fluid.load_msh("mesh.msh")
+fluid.load_msh("mesh.msh") #=>invalid pointer
 #fluid.set_open_boundary("Right",pressure = 0)
 #fluid.set_open_boundary("Left", pressure = 0)
 fluid.set_wall_boundary("Bottom",velocity=[0,0])

testcases/periodicTest/mesh.geo

 L = 1;
 H = 1;
 y = 0;
-lc = 0.05;
+lc = 1;//0.05;
 
 Point(1) = {0,0,0,lc};
 Point(2) = {L,0,0,lc};