Merge pull request #129 from MmgTools/feature/edge-tag-consistency

Edge tag consistency

CMakeLists.txt

@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 88d9d1ed6746a1e68713623859ca070224a5abbc
+    GIT_TAG 1e303f8cea3b9fcbcaf9f1b8625de7e9b434b1aa
     INSTALL_COMMAND ${CMAKE_MAKE_PROGRAM} install
     CMAKE_ARGS ${MMG_ARGS} -DUSE_ELAS=OFF ${COMPILER_CFG} ${FLAGS_CFG}
     ${SCOTCH_CFG} ${VTK_CFG} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}

src/communicators_pmmg.c

@@ -644,7 +644,7 @@ int PMMG_build_completeExtEdgeComm( PMMG_pParMesh parmesh, MPI_Comm comm ) {
 /**
  * \param parmesh pointer to parmesh structure
  * \param mesh pointer to the mesh structure
- * \param hpar hash table of parallel edges
+ * \param hpar hash table of parallel edges. Read only array.
  * \param comm pointer toward the MPI communicator to use: when called before
  * the first mesh balancing (at preprocessing stage) we have to use the
  * read_comm communicator (i.e. the communicator used to provide the inputs).

src/grpsplit_pmmg.c

@@ -33,6 +33,7 @@
  */
 #include "parmmg.h"
 #include "metis_pmmg.h"
+#include "mmgexterns_private.h"
 
 /**
  * \param nelem number of elements in the initial group
@@ -800,6 +801,7 @@ static int PMMG_splitGrps_updateFaceCommOld( PMMG_pParMesh parmesh,
  * \param ngrp nb. of new groups
  * \param grpIdOld index of the group that is splitted in the old list of groups
  * \param grpId index of the group that we create in the list of groups
+ * \param hash table storing tags of boundary edges of original mesh.
  * \param ne number of elements in the new group mesh
  * \param np pointer toward number of points in the new group mesh
  * \param f2ifc_max maximum number of elements in the face2int_face_comm arrays
@@ -817,7 +819,8 @@ static int PMMG_splitGrps_updateFaceCommOld( PMMG_pParMesh parmesh,
  *
  */
 static int
-PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int grpIdOld,int grpId,int ne,
+PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int grpIdOld,int grpId,
+                          MMG5_HGeom hash,int ne,
                           int *np,int *f2ifc_max,int *n2inc_max,idx_t *part,
                           int* posInIntFaceComm,int* iplocFaceComm ) {
   PMMG_pGrp  const grp    = &listgrp[grpId];
@@ -1041,7 +1044,7 @@ PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int g
               pos ) ) return 0;
 
         for ( j=0; j<3; ++j ) {
-          /* Update the face and face vertices tags */
+          /** Update the face and face vertices tags */
           PMMG_tag_par_edge(pxt,MMG5_iarf[fac][j]);
           ppt = &mesh->point[tetraCur->v[MMG5_idir[fac][j]]];
           PMMG_tag_par_node(ppt);
@@ -1088,6 +1091,41 @@ PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int g
 
     if( !PMMG_splitGrps_updateNodeCommNew( parmesh,grp,mesh,meshOld,tetraCur,pt,
           tet,grpId,n2inc_max,part ) ) return 0;
+
+
+    /* Xtetra have been added by the previous loop, take advantage of the
+     * current one to update possible inconsistencies in edge tags (if a
+     * boundary face has been added along an edge that was previously boundary
+     * but not belonging to a boundary face, some of the edge tags may be
+     * missing). */
+    if ( !tetraCur->xt ) continue;
+
+    pxt = &mesh->xtetra[tetraCur->xt];
+    for ( j=0; j<6; j++ ) {
+
+      /* Tag infos have to be consistent for all edges marked as boundary */
+      if ( !(pxt->tag[j] & MG_BDY) ) continue;
+
+      int ip0 = pt->v[MMG5_iare[j][0]];
+      int ip1 = pt->v[MMG5_iare[j][1]];
+
+      uint16_t tag;
+      int      ref;
+
+      /* get the tag stored in the hash table (old mesh) and set it the xtetra
+       * edge (new mesh): hGet may return 0 as edges of the old mesh are not
+       * hashed if they were not belonging to a boundary face (but due to the
+       * new partitionning, it is possible that they are now belonging to a bdy
+       * face). */
+      MMG5_hGet( &hash, ip0, ip1, &ref, &tag );
+      pxt->tag[j] |= tag;
+
+      /* Remove spurious NOSURF tag for user required edges */
+      if ( (tag & MG_REQ) && !(tag & MG_NOSURF) ) {
+        pxt->tag[j] &= ~MG_NOSURF;
+      }
+    }
+
   }
 
   return 1;
@@ -1254,6 +1292,7 @@ int PMMG_update_oldGrps( PMMG_pParMesh parmesh ) {
  * \param ngrp number of groups to split the mesh into.
  * \param countPerGrp number of tetras in each new groups.
  * \param part array of the mesh element partitioning.
+ * \param hash table storing tags of boundary edges of original mesh.
  *
  * \return -1 : no possibility to save the mesh
  *         0  : failed but the mesh is correct
@@ -1265,7 +1304,7 @@ int PMMG_update_oldGrps( PMMG_pParMesh parmesh ) {
  *
  */
 int PMMG_split_eachGrp( PMMG_pParMesh parmesh,int grpIdOld,PMMG_pGrp grpsNew,
-                        idx_t ngrp,int *countPerGrp,idx_t *part ) {
+                        idx_t ngrp,int *countPerGrp,idx_t *part,MMG5_HGeom hash ) {
   PMMG_pGrp grpOld,grpCur;
   MMG5_pMesh meshOld,meshCur;
   /** size of allocated node2int_node_comm_idx. when comm is ready trim to
@@ -1352,7 +1391,7 @@ int PMMG_split_eachGrp( PMMG_pParMesh parmesh,int grpIdOld,PMMG_pGrp grpsNew,
     grpCur  = &grpsNew[grpId];
     meshCur = grpCur->mesh;
 
-    if ( !PMMG_splitGrps_fillGroup(parmesh,grpsNew,ngrp,grpIdOld,grpId,
+    if ( !PMMG_splitGrps_fillGroup(parmesh,grpsNew,ngrp,grpIdOld,grpId,hash,
                                    countPerGrp[grpId],&poiPerGrp[grpId],
                                    &f2ifc_max[grpId],
                                    &n2inc_max[grpId],part,posInIntFaceComm,
@@ -1461,7 +1500,7 @@ int PMMG_split_eachGrp( PMMG_pParMesh parmesh,int grpIdOld,PMMG_pGrp grpsNew,
  * \warning tetra must be packed.
  *
  */
-int PMMG_split_grps( PMMG_pParMesh parmesh,int grpIdOld,int ngrp,idx_t *part,int fitMesh ) {
+int PMMG_split_grps( PMMG_pParMesh parmesh,int grpIdOld,int ngrp,idx_t *part,int fitMesh) {
   PMMG_pGrp grpOld;
   PMMG_pGrp grpsNew = NULL;
   MMG5_pMesh meshOld;
@@ -1474,6 +1513,34 @@ int PMMG_split_grps( PMMG_pParMesh parmesh,int grpIdOld,int ngrp,idx_t *part,int
   grpOld = &parmesh->listgrp[grpIdOld];
   meshOld = parmesh->listgrp[grpIdOld].mesh;
 
+
+  /* Create hash table to store edge tags (to keep tag consistency along
+   * boundary edges that doesn't belong to a boundary face in meshOld (and
+   * doesn't has valid tags) but that will belongs to a PARBDY face after group
+   * splitting). Such kind of inconsistencies may be detected by calling the
+   * MMG3D_chkmesh function. */
+  MMG5_HGeom hash;
+  if ( !MMG5_hNew(meshOld, &hash, 6*meshOld->xt, 8*meshOld->xt) ) return 0;
+
+  int k,j,i;
+  for ( k=1; k<=meshOld->ne; k++ ) {
+    MMG5_pTetra pt = &meshOld->tetra[k];
+    if ( !pt->xt ) continue;
+
+    MMG5_pxTetra pxt = &meshOld->xtetra[pt->xt];
+    for ( j=0; j<4; j++ ) {
+      /* We recover edge tag infos from boundary faces */
+      if ( !(pxt->ftag[j] & MG_BDY) ) continue;
+
+      for ( i=0; i<3; ++i ) {
+        int ia = MMG5_iarf[j][i];
+        int ip0 = pt->v[MMG5_iare[ia][0]];
+        int ip1 = pt->v[MMG5_iare[ia][1]];
+        if( !MMG5_hEdge( meshOld, &hash, ip0, ip1, 0, pxt->tag[ia] ) ) return 0;
+      }
+    }
+  }
+
   /* count_per_grp: count new tetra per group, and store new ID in the old
    * tetra flag */
   PMMG_CALLOC(parmesh,countPerGrp,ngrp,int,"counter buffer ",return 0);
@@ -1485,7 +1552,9 @@ int PMMG_split_grps( PMMG_pParMesh parmesh,int grpIdOld,int ngrp,idx_t *part,int
 
 
   /** Perform group splitting */
-  ret_val = PMMG_split_eachGrp( parmesh,grpIdOld,grpsNew,ngrp,countPerGrp,part );
+  ret_val = PMMG_split_eachGrp( parmesh,grpIdOld,grpsNew,ngrp,countPerGrp,part,hash );
+  PMMG_DEL_MEM( meshOld, hash.geom, MMG5_hgeom, "Edge hash table" );
+
   if( ret_val != 1) goto fail_counters;
 
   PMMG_listgrp_free(parmesh, &parmesh->listgrp, parmesh->ngrp);
@@ -1507,7 +1576,6 @@ int PMMG_split_grps( PMMG_pParMesh parmesh,int grpIdOld,int ngrp,idx_t *part,int
   PMMG_DEL_MEM(parmesh,countPerGrp,int,"counter buffer ");
 
 #ifndef NDEBUG
-  int i;
   for( i = 0; i < parmesh->ngrp; i++ )
     PMMG_MEM_CHECK(parmesh,parmesh->listgrp[i].mesh,return 0);
 #endif
@@ -1781,6 +1849,13 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   if ( !ier ) {
     fprintf(stderr,"\n  ## Merge groups problem.\n");
   }
+  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;
+    }
+  }
+
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));
@@ -1845,9 +1920,23 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   }
 
   /** Split the group into the suitable number of groups */
+  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;
+    }
+  }
+
   if ( ier )
     ier = PMMG_splitPart_grps(parmesh,target,fitMesh,repartitioning_mode);
 
+  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;
+    }
+  }
+
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);

src/hash_pmmg.c

@@ -144,7 +144,7 @@ int PMMG_hashOldPar_pmmg( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_Hash *hash
  * Hash the parallel edges. Only use face communicators to this purpose.
  *
  */
-int PMMG_hashPar_pmmg( PMMG_pParMesh parmesh,MMG5_HGeom *pHash ) {
+int PMMG_hashPar_fromFaceComm( PMMG_pParMesh parmesh,MMG5_HGeom *pHash ) {
   PMMG_pGrp    grp = &parmesh->listgrp[0];
   MMG5_pMesh   mesh = grp->mesh;
   MMG5_pTetra  pt;
@@ -181,14 +181,18 @@ int PMMG_hashPar_pmmg( PMMG_pParMesh parmesh,MMG5_HGeom *pHash ) {
 /**
  * \param mesh pointer toward a MMG5 mesh structure.
  * \param pHash pointer to the edge hash table.
+ *
  * \return PMMG_FAILURE
  *         PMMG_SUCCESS
  *
- * Hash the edges. Use the assumption that all paralle edges are seen by a
- * MG_PARBDY face on an xtetra.
+ * Hash the edges belonging to parallel faces and store their tags with the
+ * additionnal MG_PARBDY tag.
+ *
+ * \remark Use the assumption that all paralle edges are
+ * seen by a MG_PARBDY face on an xtetra.
  *
  */
-int PMMG_hashPar( MMG5_pMesh mesh,MMG5_HGeom *pHash ) {
+int PMMG_hashParTag_fromXtet( MMG5_pMesh mesh,MMG5_HGeom *pHash ) {
   MMG5_pTetra  pt;
   MMG5_pxTetra pxt;
   int          k,na;
@@ -246,7 +250,7 @@ int PMMG_bdryUpdate( MMG5_pMesh mesh )
   assert ( !mesh->htab.geom );
 
   /* Hash the MG_PARBDY edges */
-  if( PMMG_hashPar(mesh,&hash) != PMMG_SUCCESS ) return PMMG_FAILURE;
+  if( PMMG_hashParTag_fromXtet(mesh,&hash) != PMMG_SUCCESS ) return PMMG_FAILURE;
 
   /** Update xtetra edge tag if needed */
   for (k=1; k<=mesh->ne; ++k) {

src/libparmmg.c

@@ -226,6 +226,8 @@ 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
  *

src/libparmmg1.c

@@ -269,7 +269,7 @@ 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__);
@@ -279,7 +279,7 @@ int PMMG_packParMesh( PMMG_pParMesh parmesh )
         fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
         return 0;
       }
-    }
+//    }
   }
 
   return 1;
@@ -401,7 +401,7 @@ int  PMMG_update_face2intInterfaceTetra( PMMG_pParMesh parmesh, int igrp,
     }
   }
 
-  /** Step 2: Travel through the \a facesData array, get int the hash table the
+  /** Step 2: Travel through the \a facesData array, get in the hash table the
    * index of the element to which belong the face and update the face
    * communicator */
   face2int_face_comm_index1 = grp->face2int_face_comm_index1;
@@ -535,6 +535,7 @@ 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
@@ -567,6 +568,16 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
   /** Set inputMet flag */
   parmesh->info.inputMet = 0;
+
+#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
+
   for ( i=0; i<parmesh->ngrp; ++i ) {
     met         = parmesh->listgrp[i].met;
     if ( met && met->m ) {
@@ -622,8 +633,16 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   for ( i=0; i<parmesh->ngrp; ++i ) {
     mesh         = parmesh->listgrp[i].mesh;
 
+    PMMG_grp_to_saveMesh( parmesh, i, "AfterSplit" );
+
     if ( !mesh ) continue;
 
+    if ( !MMG5_chkmsh(mesh,1,1) ) {
+      fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
+      return 0;
+    }
+
+
     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));
 
@@ -665,10 +684,20 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
       met          = parmesh->listgrp[i].met;
       field        = parmesh->listgrp[i].field;
 
+
+      if ( !MMG5_chkmsh(mesh,1,1) ) {
+        fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
+        return 0;
+      }
+
+
 #warning Luca: until analysis is not ready
 #ifdef USE_POINTMAP
-      for( k = 1; k <= mesh->np; k++ )
+      for( k = 1; k <= mesh->np; k++ ) {
+#warning Algiane todo
+        //assert ( mesh->point[k].src == k);
         mesh->point[k].src = k;
+      }
 #endif
 
       /* Reset the value of the fem mode */
@@ -939,6 +968,14 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   ier = PMMG_merge_grps(parmesh,0);
   MPI_Allreduce( &ier, &ieresult, 1, MPI_INT, MPI_MIN, parmesh->comm );
 
+  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;
+    }
+  }
+
+
   if ( parmesh->info.imprim > PMMG_VERB_STEPS ) {
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);