Merge pull request #131 from MmgTools/feature/split_ls_multimat

Add multimat mode to ls split

cmake/testing/pmmg_tests.cmake

@@ -21,7 +21,7 @@ IF( BUILD_TESTING )
       ENDIF()
       EXECUTE_PROCESS(
         COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} fetch
-        COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} checkout cd45788c32f6
+        COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} checkout 5091f86924742
         TIMEOUT 20
         WORKING_DIRECTORY ${CI_DIR}
         #COMMAND_ECHO STDOUT
@@ -783,6 +783,23 @@ IF( BUILD_TESTING )
 
 
       endforeach()
+
+    endforeach()
+
+    foreach( MODE faces )
+
+      # Toy geom nosplit:  2 procs, ls_val=0.0 + remesh hsiz 0.1 + 4 iter + multimat nosplit
+      ## Remark : FAIL inside mmg3d scotch renum of iter 5 if niter = 5
+
+      SET( NP 2 )
+      add_test( NAME ls-DisIn-toygeom-nosplit-${MODE}-${NP}
+        COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
+        ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-distributed-${MODE}-mat-edges.mesh -v 5
+        -f ${CI_DIR}/LevelSet/${NP}p_toygeom/nosplit.mmg3d
+        -hsiz 0.1 -niter 4
+        -ls 0.0
+        -sol ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-ls.sol
+        -out ${CI_DIR_RESULTS}/ls-DisIn-toygeom-nosplit-${MODE}-${NP}.o.mesh)
     endforeach()
 
     #***********************

src/analys_pmmg.c

@@ -2788,7 +2788,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   /* Hash parallel edges from tetra and face communicator: store edges and the
      MG_PARBDY tag (other edge tags are not stored).*/
   if( ier && (PMMG_hashPar_fromFaceComm( parmesh,&hpar ) != PMMG_SUCCESS) ) {
-    fprintf(stderr,"\n  ## Warning: impossible to compute the hash parallel edge\n");
+    fprintf(stderr,"\n  ## Impossible to compute the hash parallel edge."
+            " Exit program.\n");
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,hpar.geom);
@@ -2802,7 +2803,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   }
 
   if ( !PMMG_build_edgeComm( parmesh,mesh,&hpar,comm ) ) {
-    fprintf(stderr,"\n  ## Warning: Impossible to build edge communicator\n");
+    fprintf(stderr,"\n  ## Impossible to build edge communicator."
+            " Exit program\n");
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,hpar.geom);
@@ -2811,7 +2813,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 
   /* Compute global node numbering and store it in ppt->tmp */
   if( !PMMG_Compute_verticesGloNum( parmesh,comm ) ) {
-    fprintf(stderr,"\n  ## Warning: impossible to compute node global numbering\n");
+    fprintf(stderr,"\n  ## Impossible to compute node global numbering."
+            " Exit program\n");
 
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
@@ -3099,7 +3102,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 
 #ifdef USE_POINTMAP
   /* Initialize source point with input index */
-  int ip;
+  MMG5_int ip;
   for( ip = 1; ip <= mesh->np; ip++ )
     mesh->point[ip].src = ip;
 #endif

src/grpsplit_pmmg.c

@@ -1849,13 +1849,15 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   if ( !ier ) {
     fprintf(stderr,"\n  ## Merge groups problem.\n");
   }
+
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
-
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));
@@ -1920,22 +1922,26 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   }
 
   /** Split the group into the suitable number of groups */
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( ier )
     ier = PMMG_splitPart_grps(parmesh,target,fitMesh,repartitioning_mode);
 
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));

src/hash_pmmg.c

@@ -333,17 +333,19 @@ int PMMG_chkBdryTria(MMG5_pMesh mesh, MMG5_int* permtria) {
  * \param hash pointer toward the hash table of edges.
  * \param a index of the first extremity of the edge.
  * \param b index of the second extremity of the edge.
- * \param k index of new point along the edge [a,b].
- * \param s If ls mode in ParMmg: index of new point in internal edge communicator;
- *          otherwise, the value stored in variable s.
+ *
+ * \param s If ls mode: 1 for a parallel edge that belongs
+ * to at least one element whose reference has to be splitted (either because we
+ * are not in multi-mat mode or because the reference is split in multi-mat
+ * mode). To avoid useless checks, some non parallel edges may be marked.
+ * If the edge belongs only to non-split references, s has to be 0.
+ *
  * \return PMMG_SUCCESS if success, PMMG_FAILURE if fail (edge is not found).
  *
- * Update the index of the new point stored along the edge \f$[a;b]\f$ (similar to MMG5_hashUpdate in mmg).
- * If ls mode in ParMmg: update the index of the new point in internal edge communicator in variable s;
- * otherwise, update the value stored in variable s.
+ * Update the value of the s field stored along the edge \f$[a;b]\f$
  *
  */
-int PMMG_hashUpdate_all(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int k,MMG5_int s) {
+int PMMG_hashUpdate_s(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int s) {
   MMG5_hedge  *ph;
   MMG5_int     key;
   MMG5_int    ia,ib;
@@ -355,7 +357,6 @@ int PMMG_hashUpdate_all(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int k,MMG5_i
 
   while ( ph->a ) {
     if ( ph->a == ia && ph->b == ib ) {
-      ph->k = k;
       ph->s = s;
       return PMMG_SUCCESS;
     }

src/libparmmg.c

@@ -191,6 +191,10 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
     if ( !MMG3D_mmg3d2(mesh,ls,met) ) {
       return PMMG_STRONGFAILURE;
     }
+    /* Update mesh->npi and mesh->nei to be equal to mesh->np and mesh->ne, respectively */
+    mesh->npi = mesh->np;
+    mesh->nei = mesh->ne;
+
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);
     if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
@@ -226,8 +230,6 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
   return PMMG_SUCCESS;
 }
 
-int PMMG_grp_to_saveMesh( PMMG_pParMesh parmesh, int i, char*  );
-
 /**
  * \param  parmesh pointer to parmesh structure
  *
@@ -1116,9 +1118,13 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
   }
 
   /* Store owner in the point flag */
+  MMG5_int np_overlap = 0;
   for( ip = 1; ip <= mesh->np; ip++ ) {
     ppt = &mesh->point[ip];
-    if (ppt->tag & MG_OVERLAP) continue;
+    if (ppt->tag & MG_OVERLAP) {
+      ++np_overlap;
+      continue;
+    }
     ppt->flag = parmesh->myrank;
   }
 
@@ -1130,7 +1136,7 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
   }
 
   /* Count owned nodes */
-  nowned = mesh->np;
+  nowned = mesh->np - np_overlap;
   for( idx = 0; idx < int_node_comm->nitem; idx++ ) {
     if( intvalues[idx] != parmesh->myrank ) nowned--;
   }
@@ -1228,7 +1234,12 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
                  PMMG_DEL_MEM(parmesh,status,MPI_Status,"mpi_status");
                  PMMG_destroy_int(parmesh,ptr_int,nptr,"vertGlobNum");
                  MPI_Abort(parmesh->comm,PMMG_TMPFAILURE) );
-      for( i = 0; i < nitem; i++ ) assert(itorecv[i]);
+
+#ifndef NDEBUG
+      for( i = 0; i < nitem; i++ ) {
+        assert(itorecv[i]);
+      }
+#endif
     }
   }
 

src/libparmmg1.c

@@ -189,8 +189,7 @@ int PMMG_packTetra( PMMG_pParMesh parmesh, int igrp ) {
  *
  * \return 0 if fail, 1 otherwise
  *
- * Pack the sparse meshes of each group and create triangles and edges before
- * getting out of library
+ * Pack the sparse meshes of each group
  *
  */
 int PMMG_packParMesh( PMMG_pParMesh parmesh )
@@ -269,17 +268,19 @@ int PMMG_packParMesh( PMMG_pParMesh parmesh )
     }
 
     /* to could save the mesh, the adjacency have to be correct */
-//    if ( mesh->info.ddebug ) {
+    if ( mesh->info.ddebug ) {
       if ( (!mesh->adja) && !MMG3D_hashTetra(mesh,1) ) {
         fprintf(stderr,"\n  ## Error: %s: tetra hashing problem. Exit program.\n",
                 __func__);
         return 0;
       }
+#ifndef NDEBUG
       if ( !MMG5_chkmsh(mesh,1,1) ) {
         fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
         return 0;
       }
-//    }
+#endif
+    }
   }
 
   return 1;
@@ -535,7 +536,6 @@ int PMMG_scotchCall( PMMG_pParMesh parmesh,int igrp,int *permNodGlob ) {
   return 1;
 }
 
-int PMMG_grp_to_saveMesh( PMMG_pParMesh parmesh, int i, char*  );
 /**
  * \param parmesh pointer toward a parmesh structure where the boundary entities
  * are stored into xtetra and xpoint strucutres
@@ -630,18 +630,20 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
   /** Reset the boundary fields between the old mesh size and the new one (Mmg
    * uses this fields assiming they are setted to 0)/ */
+
+  permNodGlob = NULL;
+
   for ( i=0; i<parmesh->ngrp; ++i ) {
     mesh         = parmesh->listgrp[i].mesh;
 
-    PMMG_grp_to_saveMesh( parmesh, i, "AfterSplit" );
-
     if ( !mesh ) continue;
 
+#ifndef NDEBUG
     if ( !MMG5_chkmsh(mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
-
+#endif
 
     memset(&mesh->xtetra[mesh->xt+1],0,(mesh->xtmax-mesh->xt)*sizeof(MMG5_xTetra));
     memset(&mesh->xpoint[mesh->xp+1],0,(mesh->xpmax-mesh->xp)*sizeof(MMG5_xPoint));
@@ -684,18 +686,16 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
       met          = parmesh->listgrp[i].met;
       field        = parmesh->listgrp[i].field;
 
-
+#ifndef NDEBUG
       if ( !MMG5_chkmsh(mesh,1,1) ) {
         fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
         return 0;
       }
-
+#endif
 
 #warning Luca: until analysis is not ready
 #ifdef USE_POINTMAP
       for( k = 1; k <= mesh->np; k++ ) {
-#warning Algiane todo
-        //assert ( mesh->point[k].src == k);
         mesh->point[k].src = k;
       }
 #endif
@@ -722,7 +722,13 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
 #ifdef USE_SCOTCH
         /* Allocation of the array that will store the node permutation */
-        PMMG_MALLOC(parmesh,permNodGlob,mesh->np+1,int,"node permutation",
+        // npi stores the number of points when we enter Mmg, np stores the
+        // number of points after adatptation.
+        // In theorie, here np == npi
+
+        assert ( mesh->np == mesh->npi );
+
+        PMMG_MALLOC(parmesh,permNodGlob,mesh->npi+1,int,"node permutation",
                     PMMG_scotch_message(&warnScotch) );
         if ( permNodGlob ) {
           for ( k=1; k<=mesh->np; ++k ) {
@@ -768,6 +774,7 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           }
         }
 
+
 #ifdef PATTERN
         ier = MMG5_mmg3d1_pattern( mesh, met, permNodGlob );
 #else
@@ -834,14 +841,15 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           goto strong_failed;
         }
 
+#ifdef USE_SCOTCH
+        PMMG_DEL_MEM(parmesh,permNodGlob,int,"node permutation");
+#endif
+
         if ( !ier ) { break; }
+
       }
       /* Reset the mesh->gap field in case Mmg have modified it */
       mesh->gap = MMG5_GAP;
-
-#ifdef USE_SCOTCH
-      PMMG_DEL_MEM(parmesh,permNodGlob,int,"node permutation");
-#endif
     }
 
     MPI_Allreduce( &ier, &ieresult, 1, MPI_INT, MPI_MIN, parmesh->comm );
@@ -968,13 +976,14 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   ier = PMMG_merge_grps(parmesh,0);
   MPI_Allreduce( &ier, &ieresult, 1, MPI_INT, MPI_MIN, parmesh->comm );
 
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
-
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_STEPS ) {
     chrono(OFF,&(ctim[tim]));

src/loadbalancing_pmmg.c

@@ -44,8 +44,6 @@
  * Load balancing of the mesh groups over the processors.
  *
  */
-int PMMG_grp_to_saveMesh( PMMG_pParMesh parmesh, int i, char*  );
-
 int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   MMG5_pMesh mesh;
   int        ier,ier_glob,igrp,ne;
@@ -166,7 +164,7 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   }
 
 #ifndef NDEBUG
-  for ( int i=0; i<parmesh->ngrp; ++i ) {
+  for ( i=0; i<parmesh->ngrp; ++i ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
@@ -194,12 +192,15 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
       }
     }
   }
-  for (int k=0; k<parmesh->ngrp; ++k ) {
-    if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
+
+#ifndef NDEBUG
+  for (i=0; i<parmesh->ngrp; ++i ) {
+    if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));