Merge branch 'develop' into feature/fix-cuttet

MmgTools · Sep 13, 2024 · 118f4c1 · 118f4c1
2 parents 4a4b5fd + 0f82307
commit 118f4c1
Show file tree

Hide file tree

Showing 6 changed files with 32 additions and 121 deletions.
diff --git a/cmake/testing/pmmg_tests.cmake b/cmake/testing/pmmg_tests.cmake
@@ -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

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
@@ -1477,115 +1477,6 @@ int PMMG_hashNorver( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *hash,
   return 1;
 }
 
-/**
- * \param parmesh pointer to the parmesh structure
- * \param mesh pointer to the mesh structure
- *
- * \return 1 if success, 0 if failure.
- *
- * Compute continuous geometric support (normal and tangent vectors) on
- * non-manifold MG_OLDPARBDY points.
- *
- * \remark Analogous to the MMG3D_nmgeom function, but it only travels on
- * old parallel points.
- * \remark Normal and tangent vectors on these points are overwritten.
- *
- */
-int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
-  MMG5_pTetra     pt;
-  MMG5_pPoint     p0;
-  MMG5_pxPoint    pxp;
-  int             k,base;
-  int             *adja;
-  double          n[3],t[3];
-  int             ip;
-  int8_t          i,j,ier;
-
-  for( ip = 1; ip <= mesh->np; ip++ ) {
-    mesh->point[ip].flag = mesh->base;
-  }
-
-  base = ++mesh->base;
-  for (k=1; k<=mesh->ne; k++) {
-    pt   = &mesh->tetra[k];
-    if( !MG_EOK(pt) ) continue;
-    adja = &mesh->adja[4*(k-1)+1];
-    for (i=0; i<4; i++) {
-      if ( adja[i] ) continue;
-      for (j=0; j<3; j++) {
-        ip = MMG5_idir[i][j];
-        p0 = &mesh->point[pt->v[ip]];
-        if ( p0->flag == base )  continue;
-        else if ( !(p0->tag & MG_OLDPARBDY) ) continue;
-        else if ( !(p0->tag & MG_NOM) )  continue;
-
-        p0->flag = base;
-        ier = MMG5_boulenm(mesh,k,ip,i,n,t);
-
-        if ( ier < 0 )
-          return 0;
-        else if ( !ier ) {
-          p0->tag |= MG_REQ;
-          p0->tag &= ~MG_NOSURF;
-        }
-        else {
-          if ( !p0->xp ) {
-            ++mesh->xp;
-            if(mesh->xp > mesh->xpmax){
-              MMG5_TAB_RECALLOC(mesh,mesh->xpoint,mesh->xpmax,MMG5_GAP,MMG5_xPoint,
-                                 "larger xpoint table",
-                                 mesh->xp--;
-                                 fprintf(stderr,"  Exit program.\n");return 0;);
-            }
-            p0->xp = mesh->xp;
-          }
-          pxp = &mesh->xpoint[p0->xp];
-          memcpy(pxp->n1,n,3*sizeof(double));
-          memcpy(p0->n,t,3*sizeof(double));
-        }
-      }
-    }
-  }
-  /* Deal with the non-manifold points that do not belong to a surface
-   * tetra (a tetra that has a face without adjacent)*/
-  for (k=1; k<=mesh->ne; k++) {
-    pt   = &mesh->tetra[k];
-    if( !MG_EOK(pt) ) continue;
-
-    for (i=0; i<4; i++) {
-      p0 = &mesh->point[pt->v[i]];
-      if ( !(p0->tag & MG_OLDPARBDY) ) continue;
-      else if ( p0->tag & MG_PARBDY || p0->tag & MG_REQ || !(p0->tag & MG_NOM) || p0->xp ) continue;
-      ier = MMG5_boulenmInt(mesh,k,i,t);
-      if ( ier ) {
-        ++mesh->xp;
-        if(mesh->xp > mesh->xpmax){
-          MMG5_TAB_RECALLOC(mesh,mesh->xpoint,mesh->xpmax,MMG5_GAP,MMG5_xPoint,
-                            "larger xpoint table",
-                            mesh->xp--;
-                            fprintf(stderr,"  Exit program.\n");return 0;);
-        }
-        p0->xp = mesh->xp;
-        pxp = &mesh->xpoint[p0->xp];
-        memcpy(p0->n,t,3*sizeof(double));
-      }
-      else {
-        p0->tag |= MG_REQ;
-        p0->tag &= ~MG_NOSURF;
-      }
-    }
-  }
-
-  /*for (k=1; k<=mesh->np; k++) {
-    p0 = &mesh->point[k];
-    if ( !(p0->tag & MG_NOM) || p0->xp ) continue;
-    p0->tag |= MG_REQ;
-    p0->tag &= ~MG_NOSURF;
-  }*/
-
-  return 1;
-}
-
 static inline
 int MMG5_skip_nonOldParBdy ( int8_t tag ) {
   return !(tag & MG_OLDPARBDY);
@@ -2403,7 +2294,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 +2329,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;
   }

diff --git a/src/communicators_pmmg.c b/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 );
 
 

diff --git a/src/hash_pmmg.c b/src/hash_pmmg.c
@@ -273,4 +273,5 @@ int PMMG_bdryUpdate( MMG5_pMesh mesh )
     PMMG_DEL_MEM(mesh,mesh->edge,MMG5_Edge,"deallocating edges");
 
   return PMMG_SUCCESS;
+
 }
diff --git a/src/libparmmg.c b/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;
       }
     }

diff --git a/src/parmmg.h b/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 );
@@ -513,7 +514,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 );