Merge pull request #110 from coprigent/feature/permutation-array

Fix problems with extra triangles in input mesh

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 0b1fce5e92fd42514fd75e8542ff67adeb6d77f8
+        COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} checkout a8b02c1196945bab072664f4f30f053dabb2639c
         TIMEOUT 20
         WORKING_DIRECTORY ${CI_DIR}
         #COMMAND_ECHO STDOUT
@@ -434,6 +434,13 @@ IF( BUILD_TESTING )
       ${CI_DIR}/LevelSet/2p_toygeom/cube-distributed-faces-nomat-1edge.mesh -v 10 -hsiz 0.1
       -out ${CI_DIR_RESULTS}/update-ref-tag.o.mesh)
 
+    # Test to check that when not using -opnbdy option, internal triangles are correctly removed.
+    # See ParMmg PR#110
+    add_test( NAME extrainternaltriangles
+      COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} 3 $<TARGET_FILE:${PROJECT_NAME}>
+      ${CI_DIR}/Cube/internaltriangles-P3.mesh -v 10
+      -out ${CI_DIR_RESULTS}/internaltriangles-P3.o.mesh)
+
     ###############################################################################
     #####
     #####        Test isovalue mode - ls discretization

src/analys_pmmg.c

@@ -2403,7 +2403,7 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
  *
  * Check all boundary triangles.
  */
-int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh) {
+int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_int *permtria) {
   int       ier;
 
   /**--- stage 1: data structures for surface */
@@ -2438,7 +2438,7 @@ int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh) {
   }
 
   /* identify surface mesh */
-  if ( !MMG5_chkBdryTria(mesh) ) {
+  if ( !PMMG_chkBdryTria(mesh,permtria) ) {
     fprintf(stderr,"\n  ## Boundary problem. Exit program.\n");
     return 0;
   }

src/communicators_pmmg.c

@@ -866,7 +866,7 @@ void PMMG_tria2elmFace_coords( PMMG_pParMesh parmesh ) {
   /* Process tria stored in index1 */
   for( i=0; i<grp->nitem_int_face_comm; i++ ) {
     kt    = grp->face2int_face_comm_index1[i];
-    ptt = &mesh->tria[kt];
+    ptt    = &mesh->tria[kt];
     ie     = ptt->cc/4;
     ifac   = ptt->cc%4;
 
@@ -1008,7 +1008,7 @@ int PMMG_build_nodeCommIndex( PMMG_pParMesh parmesh ) {
  * stored in the external face communicator.
  *
  */
-int PMMG_build_faceCommIndex( PMMG_pParMesh parmesh ) {
+int PMMG_build_faceCommIndex( PMMG_pParMesh parmesh, MMG5_int* permtria ) {
   PMMG_pGrp      grp;
   PMMG_pInt_comm int_face_comm;
   PMMG_pExt_comm ext_face_comm;
@@ -1038,7 +1038,12 @@ int PMMG_build_faceCommIndex( PMMG_pParMesh parmesh ) {
   for( iext_comm = 0; iext_comm < parmesh->next_face_comm; iext_comm++ ) {
     ext_face_comm = &parmesh->ext_face_comm[iext_comm];
     for( iext = 0; iext < ext_face_comm->nitem; iext++ ) {
-      grp->face2int_face_comm_index1[iint] = ext_face_comm->int_comm_index[iext];
+      if (permtria) {
+        grp->face2int_face_comm_index1[iint] = permtria[ext_face_comm->int_comm_index[iext]];
+      }
+      else {
+        grp->face2int_face_comm_index1[iint] = ext_face_comm->int_comm_index[iext];
+      }
       grp->face2int_face_comm_index2[iint] = iint;
       ext_face_comm->int_comm_index[iext] = iint++;
     }
@@ -1214,7 +1219,7 @@ int PMMG_build_faceCommFromNodes( PMMG_pParMesh parmesh,MPI_Comm comm ) {
   }
 
   /** 6) Set communicators indexing, convert tria index into iel face index */
-  ier = PMMG_build_faceCommIndex( parmesh );
+  ier = PMMG_build_faceCommIndex( parmesh, NULL );
   PMMG_tria2elmFace_flags( parmesh );
 
 

src/hash_pmmg.c

@@ -275,6 +275,55 @@ int PMMG_bdryUpdate( MMG5_pMesh mesh )
   return PMMG_SUCCESS;
 }
 
+/**
+ * \param mesh pointer to the mesh structure.
+ * \return 1 if success, 0 otherwise.
+ *
+ * - Remove double triangles from tria array.
+ *
+ * - Remove triangles that do not belong to a boundary (non opnbdy mode) from
+ *   tria array.
+ *
+ * - Check the matching between actual and given number of faces in the mesh:
+ * Count the number of faces in mesh and compare this number to the number of
+ *   given triangles.
+ *
+ * - If the founded number exceed the given one, add the missing
+ *   boundary triangles (call to MMG5_bdryTria). Do nothing otherwise.
+ *
+ * - Fill the adjacency relationship between prisms and tetra (fill adjapr with
+ *   a negative value to mark this special faces).
+ *
+ * - Set to required the triangles at interface betwen prisms and tet.
+ *
+ */
+int PMMG_chkBdryTria(MMG5_pMesh mesh, MMG5_int* permtria) {
+  MMG5_int       ntmesh,ntpres;
+  int            ier;
+  MMG5_Hash      hashElt;
+
+  /** Step 1: scan the mesh and count the boundaries */
+  ier = MMG5_chkBdryTria_countBoundaries(mesh,&ntmesh,&ntpres);
+
+  /** Step 2: detect the extra boundaries (that will be ignored) provided by the
+   * user */
+  if ( mesh->nt ) {
+    ier = MMG5_chkBdryTria_hashBoundaries(mesh,ntmesh,&hashElt);
+    // Travel through the tria, flag those that are not in the hash tab or
+    // that are stored more that once.
+    ier = MMG5_chkBdryTria_flagExtraTriangles(mesh,&ntpres,&hashElt);
+    // Delete flagged triangles
+    ier = MMG5_chkBdryTria_deleteExtraTriangles(mesh, permtria);
+  }
+  ntmesh +=ntpres;
+
+  /** Step 3: add the missing boundary triangles or, if the mesh contains
+   * prisms, set to required the triangles at interface betwen prisms and tet */
+  ier = MMG5_chkBdryTria_addMissingTriangles(mesh,ntmesh,ntpres);
+
+  return 1;
+}
+
 /**
  * \param mesh pointer toward the mesh structure.
  * \param hash pointer toward the hash table of edges.

src/libparmmg.c

@@ -243,6 +243,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
   int8_t     tim;
   char       stim[32];
   mytime     ctim[TIMEMAX];
+  MMG5_int   *permtria;
   int ier = PMMG_SUCCESS;
 
   /* Chrono initialization */
@@ -342,17 +343,21 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
 
   /** Mesh analysis I: Needed to create communicators
    *  Check triangles, create xtetras */
+  PMMG_CALLOC(parmesh,permtria,mesh->nt+1,MMG5_int,"permtria",return 0);
+  MMG5_int k;
+  for (k=0;k<=mesh->nt;k++) {
+    permtria[k] = k;
+  }
   if ( parmesh->myrank < parmesh->info.npartin ) {
-    if ( !PMMG_analys_tria(parmesh,mesh) ) {
+    if ( !PMMG_analys_tria(parmesh,mesh,permtria) ) {
       return PMMG_STRONGFAILURE;
     }
   }
-
   /* For both API modes, build communicators indices and set xtetra as PARBDY */
   switch( parmesh->info.API_mode ) {
     case PMMG_APIDISTRIB_faces :
       /* 1) Set face communicators indexing */
-      if( !PMMG_build_faceCommIndex( parmesh ) ) return 0;
+      if( !PMMG_build_faceCommIndex( parmesh, permtria ) ) return 0;
 
       /* Convert tria index into iel face index (it needs a valid cc field in
        * each tria), and tag xtetra face as PARBDY before the tag is transmitted
@@ -383,6 +388,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
       if ( !PMMG_build_faceCommFromNodes(parmesh,parmesh->info.read_comm) ) return PMMG_STRONGFAILURE;
       break;
   }
+  MMG5_SAFE_FREE( permtria );
 
   /** Discretization of the isovalue  */
   if (mesh->info.iso) {
@@ -409,7 +415,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
     /** Mesh analysis Ib : After LS discretization
      * Check triangles, create xtetras */
     if ( parmesh->myrank < parmesh->info.npartin ) {
-      if ( !PMMG_analys_tria(parmesh,mesh) ) {
+      if ( !PMMG_analys_tria(parmesh,mesh,permtria) ) {
         return PMMG_STRONGFAILURE;
       }
     }

src/parmmg.h

@@ -447,7 +447,8 @@ void PMMG_Analys_Init_SurfNormIndex( MMG5_pTetra pt );
 int PMMG_Analys_Get_SurfNormalIndex( MMG5_pTetra pt,int ifac,int i );
 int PMMG_boulernm(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_Hash *hash,int start,int ip,int *ng,int *nr);
 int PMMG_boulen(PMMG_pParMesh parmesh,MMG5_pMesh mesh,int start,int ip,int iface,double t[3]);
-int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh);
+int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_int *permtria);
+int PMMG_chkBdryTria(MMG5_pMesh mesh, MMG5_int* permtria);
 int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm);
 int PMMG_update_analys(PMMG_pParMesh parmesh);
 int PMMG_hashPar( MMG5_pMesh mesh,MMG5_HGeom *pHash );
@@ -516,7 +517,7 @@ void PMMG_tria2elmFace_flags( PMMG_pParMesh parmesh );
 void PMMG_tria2elmFace_coords( PMMG_pParMesh parmesh );
 int PMMG_tria_highestcoord( MMG5_pMesh mesh, MMG5_int *v_t);
 int PMMG_build_nodeCommIndex( PMMG_pParMesh parmesh );
-int PMMG_build_faceCommIndex( PMMG_pParMesh parmesh );
+int PMMG_build_faceCommIndex( PMMG_pParMesh parmesh, MMG5_int* permtria );
 int PMMG_build_nodeCommFromFaces( PMMG_pParMesh parmesh, MPI_Comm comm );
 int PMMG_build_faceCommFromNodes( PMMG_pParMesh parmesh, MPI_Comm comm );
 int PMMG_build_simpleExtNodeComm( PMMG_pParMesh parmesh );