From 839fe7711eefb868a670f882c5d8d8209dd2eebf Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane@gmail.com>
Date: Tue, 12 Dec 2023 15:10:24 +0100
Subject: [PATCH 01/36] Maintain edge tag consistency in group_split

First attempt to maintain edge tags consistency in grp_split: for now, edges of xtetra of new mesh are updated if they belong to a boundary face.
---
 src/grpsplit_pmmg.c      | 77 ++++++++++++++++++++++++++++++++++++----
 src/loadbalancing_pmmg.c |  2 +-
 2 files changed, 72 insertions(+), 7 deletions(-)

diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index f737b819..0fe219bb 100644
--- a/src/grpsplit_pmmg.c
+++ b/src/grpsplit_pmmg.c
@@ -800,6 +800,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 +818,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];
@@ -1088,6 +1090,39 @@ 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<4; j++ ) {
+      /* We recover edge tag infos from boundary faces */
+      if ( !(pxt->ftag[j] & MG_BDY) ) continue;
+
+      int i;
+      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]];
+
+        int16_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[ia] |= tag;
+      }
+    }
+
   }
 
   return 1;
@@ -1254,6 +1289,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 +1301,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 +1388,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 +1497,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 +1510,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 +1549,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 +1573,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
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index dfcbcd47..a90e3868 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -137,7 +137,7 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
     ier = PMMG_split_n2mGrps(parmesh,PMMG_GRPSPL_MMG_TARGET,0,partitioning_mode);
     if ( ier<=0 )
       fprintf(stderr,"\n  ## Problem when splitting into a lower number of groups.\n");
-    }
+  }
 
   // Algiane: Optim: is this reduce needed?
   MPI_Allreduce( &ier, &ier_glob, 1, MPI_INT, MPI_MIN, parmesh->comm);

From 451599b365a0c50cde03a81a728942c84e19d073 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane@gmail.com>
Date: Tue, 12 Dec 2023 17:14:44 +0100
Subject: [PATCH 02/36] Modif of edge tag update in updateTag.

---
 src/libparmmg1.c |  2 +-
 src/tag_pmmg.c   | 21 +++++++++++++++++++--
 2 files changed, 20 insertions(+), 3 deletions(-)

diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 16d4b22f..6c463f98 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -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;
diff --git a/src/tag_pmmg.c b/src/tag_pmmg.c
index a49fed2b..41b7c278 100644
--- a/src/tag_pmmg.c
+++ b/src/tag_pmmg.c
@@ -433,10 +433,27 @@ int PMMG_updateTag(PMMG_pParMesh parmesh) {
          * so remove the MG_NOSURF tag if the edge is truly required */
         if( pxt->tag[j] & MG_REQ )
           gettag &= ~MG_NOSURF;
+
         /* set edge tag (without NOSURF tag if the edge is required by the
          * user): here we preserve the initial MG_REQ tag of each tetra, thus,
-         * potential inconsistencies will not be solved. */
-        pxt->tag[j] |= gettag;
+         * potential inconsistencies will not be solved.
+         *
+         * A xtetra may have an edge that is boundary but doesn't belong to any
+         * boundary face:
+         *  - if this edge is marked as MG_BDY, the edge tag should be
+         *    consistent with edge tag stored from a boundary face and we have
+         *    to maintain this consistency;
+         *
+         * - if this edge is not marked as MG_BDY (tag == 0), we are not able to
+         *    know if the edge is ref or required or if it has any other tag so
+         *    we are not able to maintain the tag consistency and we have to
+         *    preserve the fact that the edge is not MG_BDY.
+         *
+         */
+        if ( (pxt->tag[j] & MG_BDY) ||
+             ( (pxt->ftag[MMG5_ifar[j][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[j][1]] & MG_BDY) ) ) {
+               pxt->tag[j] |= gettag;
+        }
       }
     }
     PMMG_DEL_MEM( mesh, hash.geom, MMG5_hgeom, "Edge hash table" );

From 8ee20c9d3398faecc8ea0f774d019e80951831aa Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane@gmail.com>
Date: Tue, 12 Dec 2023 17:19:52 +0100
Subject: [PATCH 03/36] Modif of edge tag update after group splitting.

---
 src/grpsplit_pmmg.c | 33 +++++++++++++++------------------
 1 file changed, 15 insertions(+), 18 deletions(-)

diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index 0fe219bb..f1b1d7d4 100644
--- a/src/grpsplit_pmmg.c
+++ b/src/grpsplit_pmmg.c
@@ -1100,27 +1100,24 @@ PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int g
     if ( !tetraCur->xt ) continue;
 
     pxt = &mesh->xtetra[tetraCur->xt];
-    for ( j=0; j<4; j++ ) {
-      /* We recover edge tag infos from boundary faces */
-      if ( !(pxt->ftag[j] & MG_BDY) ) continue;
+    for ( j=0; j<6; j++ ) {
 
-      int i;
-      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]];
+      /* Tag infos have to be consistent for all edges marked as boundary */
+      if ( !(pxt->tag[j] & MG_BDY) ) continue;
 
-        int16_t tag;
-        int     ref;
+      int ip0 = pt->v[MMG5_iare[j][0]];
+      int ip1 = pt->v[MMG5_iare[j][1]];
 
-        /* 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[ia] |= tag;
-      }
+      int16_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;
     }
 
   }

From 4a5e6870dbd6979ebe38d8d1addf14d609956b9f Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane@gmail.com>
Date: Tue, 12 Dec 2023 17:38:09 +0100
Subject: [PATCH 04/36] TO REVERT : Add  checks to validate the previous
 modifications.

---
 src/grpsplit_pmmg.c      | 22 ++++++++++++++++++++++
 src/libparmmg1.c         | 24 ++++++++++++++++++++++--
 src/loadbalancing_pmmg.c | 16 ++++++++++++++++
 3 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index f1b1d7d4..b7f6418d 100644
--- a/src/grpsplit_pmmg.c
+++ b/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
@@ -1843,6 +1844,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]));
@@ -1907,9 +1915,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);
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 6c463f98..0f6ae98d 100644
--- a/src/libparmmg1.c
+++ b/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;
@@ -624,6 +624,12 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
     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,6 +671,12 @@ 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++ )
@@ -939,6 +951,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);
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index a90e3868..2397eea2 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -32,6 +32,8 @@
  *
  */
 #include "parmmg.h"
+#include "mmgexterns_private.h"
+
 
 /**
  * \param parmesh pointer toward a parmesh structure
@@ -89,6 +91,14 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
     ier = PMMG_split_n2mGrps(parmesh,PMMG_GRPSPL_DISTR_TARGET,1,partitioning_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;
+    }
+  }
+
+
   /* There is mpi comms in distribute_grps thus we don't want that one proc
    * enters the function and not the other proc */
   MPI_Allreduce( &ier, &ier_glob, 1, MPI_INT, MPI_MIN, parmesh->comm);
@@ -152,6 +162,12 @@ 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) ) {
+      fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
+      return 0;
+    }
+  }
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));

From 9a8ae17f8b9d8fbc1c5e3a4b683657c07b2fc9e3 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 23 Aug 2024 16:10:26 +0200
Subject: [PATCH 05/36] Whitespace cleanup + improve code comments.

---
 src/analys_pmmg.c        | 16 +++++++++++++---
 src/loadbalancing_pmmg.c |  2 +-
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 5e084eb0..3c4adcbb 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -1551,7 +1551,7 @@ int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
   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;
@@ -1582,7 +1582,7 @@ int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
     p0->tag |= MG_REQ;
     p0->tag &= ~MG_NOSURF;
   }*/
-  
+
   return 1;
 }
 
@@ -1731,6 +1731,7 @@ int PMMG_loopr(PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
  * For all ather calls, comm has to be the communicator to use for computations.
  *
  * Check for singularities.
+ *
  * \remark Modeled after the MMG5_singul function.
  */
 int PMMG_singul(PMMG_pParMesh parmesh,MMG5_pMesh mesh,PMMG_hn_loopvar *var,MPI_Comm comm) {
@@ -2729,6 +2730,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   }
 
   /* check for ridges: check dihedral angle using adjacent triangle normals */
+  /* 1. Call \ref MMG5_setdhd to analyze edges at interfaces of 2 "true"
+   * boundary faces. Skip pure parallel faces */
   if ( mesh->info.dhd > MMG5_ANGLIM && !MMG5_setdhd(mesh) ) {
     fprintf(stderr,"\n  ## Geometry problem. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
@@ -2740,6 +2743,10 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
+  /* 2. Call PMMG_setdhd to analyze edges splitted by the parallel interface:
+   *   2.1. adds possibly missing MG_NOM tags;
+   *   2.2. computes dihedral angle and set MG_GEO (ridge) tag if needed. */
+#warning setdhd analysis the NOM edges: what if -nr is triggered (wrong NOM setting)?
   if ( mesh->info.dhd > MMG5_ANGLIM && !PMMG_setdhd( parmesh,mesh,&hpar,comm ) ) {
     fprintf(stderr,"\n  ## Geometry problem on parallel edges. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
@@ -2794,7 +2801,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* Tag parallel faces on material interfaces as boundary */
+  /* Tag parallel faces on material interfaces as boundary (ie, add \ref MG_PARBDYBDY tag) */
   if( !PMMG_parbdySet( parmesh ) ) {
     fprintf(stderr,"\n  ## Unable to recognize parallel faces on material interfaces. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
@@ -2811,6 +2818,9 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   if ( abs(mesh->info.imprim) > 5  || mesh->info.ddebug )
     fprintf(stdout,"  ** UPDATING TOPOLOGY AT NON-MANIFOLD POINTS\n");
 
+  /* 1. set non-manifold edges sharing non-intersecting multidomains as required
+     2. travel points lying along non-manifold edges and set tags depending on
+     the number of feature edges passing through the point */
   if ( !MMG5_setNmTag(mesh,&hash) ) {
     fprintf(stderr,"\n  ## Non-manifold topology problem. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index 2397eea2..e27eafca 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -143,7 +143,7 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   }
 
   if ( ier ) {
-    /** Redistribute the ngrp groups of listgrp into a higher number of groups */
+    /** Redistribute the ngrp groups of listgrp into a lower number of groups */
     ier = PMMG_split_n2mGrps(parmesh,PMMG_GRPSPL_MMG_TARGET,0,partitioning_mode);
     if ( ier<=0 )
       fprintf(stderr,"\n  ## Problem when splitting into a lower number of groups.\n");

From 7edb110d4105daf15b452ef653e15657548ea525 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 23 Aug 2024 16:11:08 +0200
Subject: [PATCH 06/36] Solve inconsistencies in // analysis of non-manifold
 edges.

Solves missing MG_NOM and MG_REF tags when a rank has a nom edge but without any physical boudary on it.
See the following test:
```
parmmg_debug TEST_OUTPUTS/sphere_nom_0-2.mesh -v 6 -centralized-output -metis-ratio 82 -noinsert -noswap -nomove -rn 0
```
---
 src/analys_pmmg.c | 41 +++++++++++++++++++++++++++++++----------
 1 file changed, 31 insertions(+), 10 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 3c4adcbb..0ae22597 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2199,14 +2199,19 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
 
     /* Get parallel edge touched by a boundary face and store normal vectors */
     for (i=0; i<3; i++) {
-      /* Skip non-manifold edges */
-      if ( (ptr->tag[i] & MG_NOM) ) continue;
 
       i1 = MMG5_inxt2[i];
       i2 = MMG5_inxt2[i1];
       if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
       idx = edg-1;
 
+      if ( (ptr->tag[i] & MG_NOM) ) {
+        /* We only need to store that the edge is non-manifold to share the
+         * information with the other MPI processes */
+        intvalues[2*idx] = 3;
+        continue;
+      }
+
       /* Count how many times the edge is seen locally */
       intvalues[2*idx]++;
       /* Do not store anything else for non-manifold */
@@ -2342,19 +2347,30 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
 
     /* Get parallel edge touched by a MG_BDY face and store normal vectors */
     for (i=0; i<3; i++) {
-      /* Skip non-manifold edges */
-      if ( (ptr->tag[i] & MG_NOM) ) continue;
 
       i1 = MMG5_inxt2[i];
       i2 = MMG5_inxt2[i1];
       if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
       idx = edg-1;
+
+      /* Skip non-manifold edges */
+#ifndef NDEBUG
+      if ( (ptr->tag[i] & MG_NOM) ) {
+        assert ( intvalues[2*idx] != 2 );
+      }
+#endif
+
       if( intvalues[2*idx] == 1 ) { /* no adjacent */
-        ptr->tag[i] |= MG_GEO + MG_NOM;
+#warning remove MG_GEO fpr Mmg consistency ?
+        /* MG_REF info is not analyzed in parallel for non-manifold edges (only
+           serially by Mmy).  As we need to ensure the tag consistency across
+           the processes and for sake of simplicity, we simply mark all the
+           MG_NOM edges as MG_REF */
+        ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
         i1 = MMG5_inxt2[i];
         i2 = MMG5_inxt2[i1];
-        mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM;
-        mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM;
+        mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM + MG_REF;
+        mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM + MG_REF;
         nr++;
       } else {
         if( (intvalues[2*idx+1] ==   PMMG_UNSET) ||
@@ -2376,11 +2392,16 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
           nr++;
         }
         if( intvalues[2*idx] > 2 ) { /* non-manifold edge */
-          ptr->tag[i] |= MG_GEO + MG_NOM;
+#warning remove MG_GEO fpr Mmg consistency ?
+          /* MG_REF info is not analyzed in parallel for non-manifold edges (only
+           serially by Mmy).  As we need to ensure the tag consistency across
+           the processes and for sake of simplicity, we simply mark all the
+           MG_NOM edges as MG_REF */
+          ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
           i1 = MMG5_inxt2[i];
           i2 = MMG5_inxt2[i1];
-          mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM;
-          mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM;
+          mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM + MG_REF;
+          mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM + MG_REF;
           nm++;
         }
       }

From f57d19425c567048a0c438bdb3a6da708331a47a Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 23 Aug 2024 16:21:57 +0200
Subject: [PATCH 07/36] Solve missing // analysis of non-manifold edges in -nr
 mode (no ridge detection).

    Along // edges, non-manifold analysis and information sharing is done inside the PMMG_setdhd function, thus it was skipped in -nr mode.

    See the following test:
    ```
    parmmg_debug TEST_OUTPUTS/sphere_nom_0-2.mesh -v 6 -centralized-output -metis-ratio 82 -noinsert -noswap -nomove -rn 0 -nr
    ```
---
 src/analys_pmmg.c | 40 +++++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 19 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 0ae22597..3e9f51e4 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2310,27 +2310,29 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
     }
   }
 
-  /** Second pass: Check dihedral angle and mark geometric edge with
-   *  2*PMMG_UNSET (3x if it is already a reference edge) */
-  for ( k = 0; k < parmesh->next_edge_comm; ++k ) {
-    ext_edge_comm = &parmesh->ext_edge_comm[k];
+  if ( mesh->info.dhd > MMG5_ANGLIM ) {
+    /** Second pass: Check dihedral angle and mark geometric edge with
+     *  2*PMMG_UNSET (3x if it is already a reference edge) */
+    for ( k = 0; k < parmesh->next_edge_comm; ++k ) {
+      ext_edge_comm = &parmesh->ext_edge_comm[k];
 
-    for ( i=0; i<ext_edge_comm->nitem; ++i ) {
-      idx  = ext_edge_comm->int_comm_index[i];
+      for ( i=0; i<ext_edge_comm->nitem; ++i ) {
+        idx  = ext_edge_comm->int_comm_index[i];
 
-      nt1 = intvalues[2*idx];
+        nt1 = intvalues[2*idx];
 
-      if( nt1 == 2 ) {
-        for( d = 0; d < 3; d++ ) {
-          n1[d] = doublevalues[6*idx+d];
-          n2[d] = doublevalues[6*idx+3+d];
-        }
-        dhd = n1[0]*n2[0] + n1[1]*n2[1] + n1[2]*n2[2];
-        if ( dhd <= mesh->info.dhd ) {
-          if( intvalues[2*idx+1] != PMMG_UNSET )
-            intvalues[2*idx+1] = 2*PMMG_UNSET;
-          else
-            intvalues[2*idx+1] = 3*PMMG_UNSET;
+        if( nt1 == 2 ) {
+          for( d = 0; d < 3; d++ ) {
+            n1[d] = doublevalues[6*idx+d];
+            n2[d] = doublevalues[6*idx+3+d];
+          }
+          dhd = n1[0]*n2[0] + n1[1]*n2[1] + n1[2]*n2[2];
+          if ( dhd <= mesh->info.dhd ) {
+            if( intvalues[2*idx+1] != PMMG_UNSET )
+              intvalues[2*idx+1] = 2*PMMG_UNSET;
+            else
+              intvalues[2*idx+1] = 3*PMMG_UNSET;
+          }
         }
       }
     }
@@ -2768,7 +2770,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
    *   2.1. adds possibly missing MG_NOM tags;
    *   2.2. computes dihedral angle and set MG_GEO (ridge) tag if needed. */
 #warning setdhd analysis the NOM edges: what if -nr is triggered (wrong NOM setting)?
-  if ( mesh->info.dhd > MMG5_ANGLIM && !PMMG_setdhd( parmesh,mesh,&hpar,comm ) ) {
+  if ( !PMMG_setdhd( parmesh,mesh,&hpar,comm ) ) {
     fprintf(stderr,"\n  ## Geometry problem on parallel edges. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,mesh->htab.geom);

From 00a564e63fe712fd8c4c3588748633f9e0e4fdee Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 23 Aug 2024 16:31:42 +0200
Subject: [PATCH 08/36] Rename PMMG_setdhd function into PMMG_setfeatures.

---
 src/analys_pmmg.c | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 3e9f51e4..4d316391 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2044,7 +2044,10 @@ int PMMG_singul(PMMG_pParMesh parmesh,MMG5_pMesh mesh,PMMG_hn_loopvar *var,MPI_C
  *
  * \return 1 if success, 0 if failure.
  *
- * Check dihedral angle to detect ridges on parallel edges.
+ * Check features along parallel edges:
+ *   - check for non-manifold edges
+ *   - check for reference edges (the edge shares 2 surfaces with different refs)
+ *   - if needed, check dihedral angle to detect ridges on parallel edges.
  *
  * The integer communicator is dimensioned to store the number of triangles seen
  * by a parallel edge on each partition, and a "flag" to check the references of
@@ -2058,7 +2061,7 @@ int PMMG_singul(PMMG_pParMesh parmesh,MMG5_pMesh mesh,PMMG_hn_loopvar *var,MPI_C
  * MG_PARBDYBDY only on the process who has them with the right orientation
  * (by PMMG_parbdyTria), so they will be processed only once.
  */
-int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm comm ) {
+int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm comm ) {
   PMMG_pGrp      grp;
   PMMG_pInt_comm int_edge_comm;
   PMMG_pExt_comm ext_edge_comm;
@@ -2766,11 +2769,13 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* 2. Call PMMG_setdhd to analyze edges splitted by the parallel interface:
-   *   2.1. adds possibly missing MG_NOM tags;
-   *   2.2. computes dihedral angle and set MG_GEO (ridge) tag if needed. */
-#warning setdhd analysis the NOM edges: what if -nr is triggered (wrong NOM setting)?
-  if ( !PMMG_setdhd( parmesh,mesh,&hpar,comm ) ) {
+  /* 2. Call PMMG_setfeatures to analyze edges splitted by the parallel interface:
+   *   2.1. adds possibly missing MG_NOM tags (add MG_REF tags to MG_NOM edges);
+   *   2.2. computes dihedral angle and set MG_GEO (ridge) tag if needed;
+   *   2.3. adds MG_REF tag for edges separating surfaces with different refs;
+   *   2.4. transfer edges tags to edge vertices.
+   */
+  if ( !PMMG_setfeatures( parmesh,mesh,&hpar,comm ) ) {
     fprintf(stderr,"\n  ## Geometry problem on parallel edges. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,mesh->htab.geom);

From f119be5b048ed247db296baf591789f931ea9d64 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 23 Aug 2024 17:09:39 +0200
Subject: [PATCH 09/36] Point to Mmg commit in which MMG5_OLD_PARBDY tags are
 ignored when checking mesh edge consistency.

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c93e373d..27ce238f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 2263f92c
+    GIT_TAG 7c75282c
     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}

From a17a3e953db9b8f7be22aede0cd3be27cfdc55e5 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 26 Aug 2024 18:19:09 +0200
Subject: [PATCH 10/36] Update Mmg tag after edge tag fix in swap23.

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 27ce238f..beac893b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 7c75282c
+    GIT_TAG a99d555c02
     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}

From 0b98e30232ba0a2852cf5dd7a7fb06b97addde68 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 2 Sep 2024 10:17:21 +0200
Subject: [PATCH 11/36] Update Mmg tag after fix of wrong edge tag (due to
 missing tag after collapse) in swap23.

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index beac893b..cd9d0a46 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG a99d555c02
+    GIT_TAG 75579a7305
     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}

From b22d42bce23e01c3f87f88687a94b1ad160fa71e Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 3 Oct 2024 09:21:14 +0200
Subject: [PATCH 12/36] Update Mmg tag after swap23 fix.

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index cd9d0a46..e61958cb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 75579a7305
+    GIT_TAG 8a844888bfeaa9aff1136a7782ff561dbd08f40c
     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}

From 9a554ceabb98a79c2a80e9037d3f5ca13de64488 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 3 Oct 2024 10:25:28 +0200
Subject: [PATCH 13/36] Update tag types (int16_t -> unint16_t).

---
 src/analys_pmmg.c        | 32 ++++++++++++++++----------------
 src/communicators_pmmg.c |  4 ++--
 src/grpsplit_pmmg.c      |  4 ++--
 src/hash_pmmg.c          |  2 +-
 src/quality_pmmg.c       |  2 +-
 src/tag_pmmg.c           |  6 +++---
 6 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 4d316391..1ae9564d 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -109,8 +109,8 @@ typedef struct {
  * stored for each of the two edge orientations.
  */
 static inline
-int PMMG_hGetOri( MMG5_HGeom *hash,int ip0,int ip1,int *ref,int16_t *color ) {
-  int16_t tag;
+int PMMG_hGetOri( MMG5_HGeom *hash,int ip0,int ip1,int *ref,uint16_t *color ) {
+  uint16_t tag;
 
   /* Get edge from hash table */
   if( !MMG5_hGet( hash,
@@ -136,8 +136,8 @@ int PMMG_hGetOri( MMG5_HGeom *hash,int ip0,int ip1,int *ref,int16_t *color ) {
  * each of the two edge orientations.
  */
 static inline
-int PMMG_hTagOri( MMG5_HGeom *hash,int ip0,int ip1,int ref,int16_t color ) {
-  int16_t tag;
+int PMMG_hTagOri( MMG5_HGeom *hash,int ip0,int ip1,int ref,uint16_t color ) {
+  uint16_t tag;
 
   /* Set bitwise tag from color */
   if( color ) {
@@ -309,11 +309,11 @@ int PMMG_hash_nearParEdges( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
 static inline
 int PMMG_hashNorver_edges( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
   MMG5_pPoint ppt[2];
-  double *doublevalues;
-  int    ia[2],ip[2],gip;
-  int    *intvalues,idx,d,edg,j,pos;
-  int16_t tag;
-  int8_t  found;
+  double   *doublevalues;
+  int      ia[2],ip[2],gip;
+  int      *intvalues,idx,d,edg,j,pos;
+  uint16_t tag;
+  int8_t   found;
 
   doublevalues = parmesh->int_node_comm->doublevalues;
   intvalues    = parmesh->int_node_comm->intvalues;
@@ -441,7 +441,7 @@ int PMMG_hashNorver_switch( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
 static inline
 int PMMG_hashNorver_sweep( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
   int        edg;
-  int16_t    color_old,color_new;
+  uint16_t   color_old,color_new;
 
   /* If non-manifold, only process exterior points */
   if( (var->ppt->tag & MG_NOM) && var->iadj ) return 1;
@@ -495,7 +495,7 @@ int PMMG_hashNorver_edge2paredge( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var,
                                   int idx ) {
   MMG5_pEdge pa;
   int       *intvalues,edg,j,i[2],ip,ip1;
-  int16_t    color_old,color_new;
+  uint16_t   color_old,color_new;
 
   /* Get internal communicator */
   intvalues = parmesh->int_edge_comm->intvalues;
@@ -543,7 +543,7 @@ int PMMG_hashNorver_paredge2edge( PMMG_pParMesh parmesh,MMG5_HGeom *hash,
   MMG5_pMesh mesh = parmesh->listgrp[0].mesh;
   MMG5_pEdge pa;
   int       *intvalues,edg,j,i[2],ip,ip1;
-  int16_t    color_old,color_new;
+  uint16_t   color_old,color_new;
 
   assert( parmesh->ngrp == 1 );
 
@@ -1207,7 +1207,7 @@ int PMMG_set_edge_owners( PMMG_pParMesh parmesh,MMG5_HGeom *hpar,MPI_Comm comm )
   MMG5_pEdge     pa;
   int            *intvalues,*itosend,*itorecv;
   int            idx,k,nitem,color,edg,ia,ie,ifac,ip[2],i;
-  int16_t        tag;
+  uint16_t       tag;
   MPI_Status     status;
 
   assert( parmesh->ngrp == 1 );
@@ -1587,7 +1587,7 @@ int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
 }
 
 static inline
-int MMG5_skip_nonOldParBdy ( int8_t tag ) {
+uint16_t MMG5_skip_nonOldParBdy ( uint16_t tag ) {
   return !(tag & MG_OLDPARBDY);
 }
 
@@ -1670,7 +1670,7 @@ int PMMG_loopr(PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
   MMG5_pPoint ppt[2];
   double      *doublevalues;
   int         *intvalues,ip[2],k,j,idx,ns0,edg,d;
-  int16_t     tag;
+  uint16_t    tag;
   int8_t      isEdg;
 
   /* Get node communicator */
@@ -2073,7 +2073,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
   int            k,ne,nr,nm,j;
   int            i,i1,i2;
   int            idx,edg,d;
-  int16_t        tag;
+  uint16_t       tag;
   MPI_Status     status;
 
   assert( parmesh->ngrp == 1 );
diff --git a/src/communicators_pmmg.c b/src/communicators_pmmg.c
index ad18fb13..955f6b21 100644
--- a/src/communicators_pmmg.c
+++ b/src/communicators_pmmg.c
@@ -262,7 +262,7 @@ int PMMG_fillExtEdgeComm_fromFace( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HG
                                    PMMG_pExt_comm ext_edge_comm,MMG5_pTetra pt,int ifac,int iloc,int j,int color,int *item ) {
   MMG5_pEdge pa;
   int        edg;
-  int16_t    tag;
+  uint16_t   tag;
   int8_t     i1,i2;
 
   /* Take the edge opposite to vertex iloc+j on face ifac */
@@ -661,7 +661,7 @@ int PMMG_build_edgeComm( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *hpar,
   MMG5_hgeom     *ph;
   int            *nitems_ext_comm,color,k,i,idx,ie,ifac,iloc,j,item;
   int            edg;
-  int16_t        tag;
+  uint16_t       tag;
   int8_t         ia,i1,i2;
 
   assert( parmesh->ngrp == 1 );
diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index b7f6418d..baf3cb9c 100644
--- a/src/grpsplit_pmmg.c
+++ b/src/grpsplit_pmmg.c
@@ -1109,8 +1109,8 @@ PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int g
       int ip0 = pt->v[MMG5_iare[j][0]];
       int ip1 = pt->v[MMG5_iare[j][1]];
 
-      int16_t tag;
-      int     ref;
+      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
diff --git a/src/hash_pmmg.c b/src/hash_pmmg.c
index fc30f658..63ed9424 100644
--- a/src/hash_pmmg.c
+++ b/src/hash_pmmg.c
@@ -239,7 +239,7 @@ int PMMG_bdryUpdate( MMG5_pMesh mesh )
   MMG5_pxTetra pxt;
   MMG5_HGeom   hash;
   int          k,edg;
-  int16_t      tag;
+  uint16_t     tag;
   int8_t       i,i1,i2;
 
 
diff --git a/src/quality_pmmg.c b/src/quality_pmmg.c
index 737eb1e3..8bf4373a 100644
--- a/src/quality_pmmg.c
+++ b/src/quality_pmmg.c
@@ -393,7 +393,7 @@ int PMMG_computePrilen( PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_pSol met, do
   double          len;
   int             i,k,ia,np,nq,n;
   int             ref;
-  int16_t         tag;
+  uint16_t        tag;
   int8_t          i0,i1,ier;
   static double   bd[9]= {0.0, 0.3, 0.6, 0.7071, 0.9, 1.3, 1.4142, 2.0, 5.0};
 
diff --git a/src/tag_pmmg.c b/src/tag_pmmg.c
index 41b7c278..a086e549 100644
--- a/src/tag_pmmg.c
+++ b/src/tag_pmmg.c
@@ -68,8 +68,8 @@ void PMMG_tag_par_edge(MMG5_pxTetra pxt,int j){
  * Tag an edge as parallel.
  */
 int PMMG_tag_par_edge_hash(MMG5_pTetra pt,MMG5_HGeom hash,int ia){
-  int     ip0,ip1,getref;
-  int16_t gettag;
+  int      ip0,ip1,getref;
+  uint16_t gettag;
 
   ip0 = pt->v[MMG5_iare[ia][0]];
   ip1 = pt->v[MMG5_iare[ia][1]];
@@ -275,7 +275,7 @@ int PMMG_updateTag(PMMG_pParMesh parmesh) {
   MMG5_HGeom      hash;
   int             *node2int_node_comm0_index1,*face2int_face_comm0_index1;
   int             grpid,iel,ifac,ia,ip0,ip1,k,j,i,getref;
-  int16_t         gettag;
+  uint16_t        gettag;
   int8_t          isbdy;
 
   /* Loop on groups */

From d09630a42ca67fef441e41743eeb573c3dbe7644 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 7 Oct 2024 15:05:36 +0200
Subject: [PATCH 14/36] Update comments + remove useless preprocessor test.

---
 src/analys_pmmg.c | 8 ++++++--
 src/ls_pmmg.c     | 8 +++-----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 1ae9564d..5110b356 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2107,7 +2107,9 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
       idx = edg-1;
 
-      /* Store edge tag in the internal communicator */
+      /* Mark required in the internal communicator (note that we need intvalue
+       * to be a signed integer array while Mmg tags are unsigned ints so we
+       * must be careful to not directly store the edge tag).*/
       if( (ptr->tag[i] & MG_REQ) && !(ptr->tag[i] & MG_NOSURF) ) {
         intvalues[idx] |= MG_REQ;
       }
@@ -2140,7 +2142,6 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
     ext_edge_comm = &parmesh->ext_edge_comm[k];
 
     itorecv = ext_edge_comm->itorecv;
-    rtorecv = ext_edge_comm->rtorecv;
 
     for ( i=0; i<ext_edge_comm->nitem; ++i ) {
       idx  = ext_edge_comm->int_comm_index[i];
@@ -2774,6 +2775,9 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
    *   2.2. computes dihedral angle and set MG_GEO (ridge) tag if needed;
    *   2.3. adds MG_REF tag for edges separating surfaces with different refs;
    *   2.4. transfer edges tags to edge vertices.
+   *
+   * Analysis uses the tria array so if it is not allocated, this part of the
+   * analysis will not be complete.
    */
   if ( !PMMG_setfeatures( parmesh,mesh,&hpar,comm ) ) {
     fprintf(stderr,"\n  ## Geometry problem on parallel edges. Exit program.\n");
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 38c1ab8e..63946db8 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -1541,6 +1541,7 @@ int PMMG_snpval_ls(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_pSol sol) {
 int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met) {
   char str[16]="";
   MMG5_HGeom hpar;
+  MMG5_int k;
 
   /* Set function pointers */
   /** \todo TODO :: Surface ls and alias functions */
@@ -1560,7 +1561,6 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
   }
 
   /* Modify the value of the level-set to work with the 0 level-set  */
-  MMG5_int k;
   for (k=1; k<= sol->np; k++)
     sol->m[k] -= mesh->info.ls;
 
@@ -1652,7 +1652,7 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
     return 0;
   }
 
-  /* Not sure which function to be used to deallocate memory */
+  /* Delete outdated arrays */
   MMG5_DEL_MEM(mesh,mesh->adja);
   MMG5_DEL_MEM(mesh,mesh->adjt);
   MMG5_DEL_MEM(mesh,mesh->tria);
@@ -1678,7 +1678,7 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
   }
 
   /* Clean old bdy analysis */
-  for ( MMG5_int k=1; k<=mesh->np; ++k ) {
+  for ( k=1; k<=mesh->np; ++k ) {
     if ( mesh->point[k].tag & MG_BDY ) {
       mesh->point[k].tag &= ~MG_BDY;
     }
@@ -1690,13 +1690,11 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
   /* Clean memory */
   MMG5_DEL_MEM(mesh,sol->m);
 
-#ifndef NDEBUG
   /* Check communicators */
   assert ( PMMG_check_extFaceComm ( parmesh,parmesh->info.read_comm ) );
   assert ( PMMG_check_intFaceComm ( parmesh ) );
   assert ( PMMG_check_extNodeComm ( parmesh,parmesh->info.read_comm ) );
   assert ( PMMG_check_intNodeComm ( parmesh ) );
-#endif
 
   /* Dealloc edge comm  as it is not up-to-date */
   MMG5_DEL_MEM(mesh,hpar.geom);

From d2ccead996ecb347af85b1ae778c1a806811870b Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 10 Oct 2024 17:29:50 +0200
Subject: [PATCH 15/36] Refactoring in PMMG_setfeatures function.

---
 src/analys_pmmg.c | 32 ++++++++++++++------------------
 1 file changed, 14 insertions(+), 18 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 5110b356..70448dbb 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2356,7 +2356,11 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 
       i1 = MMG5_inxt2[i];
       i2 = MMG5_inxt2[i1];
-      if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
+
+      MMG5_int ip1 = ptr->v[i1];
+      MMG5_int ip2 = ptr->v[i2];
+
+      if ( !MMG5_hGet( pHash, ip1, ip2, &edg, &tag ) ) continue;
       idx = edg-1;
 
       /* Skip non-manifold edges */
@@ -2373,41 +2377,33 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
         ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
-        i1 = MMG5_inxt2[i];
-        i2 = MMG5_inxt2[i1];
-        mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM + MG_REF;
-        mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM + MG_REF;
+        mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
+        mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
         nr++;
       } else {
         if( (intvalues[2*idx+1] ==   PMMG_UNSET) ||
             (intvalues[2*idx+1] == 3*PMMG_UNSET) ) { /* reference edge */
           ptr->tag[i]   |= MG_REF;
-          i1 = MMG5_inxt2[i];
-          i2 = MMG5_inxt2[i1];
-          mesh->point[ptr->v[i1]].tag |= MG_REF;
-          mesh->point[ptr->v[i2]].tag |= MG_REF;
+          mesh->point[ip1].tag |= MG_REF;
+          mesh->point[ip2].tag |= MG_REF;
           ne++;
         }
         if( (intvalues[2*idx+1] == 2*PMMG_UNSET) ||
             (intvalues[2*idx+1] == 3*PMMG_UNSET) ) { /* geometric edge */
           ptr->tag[i]   |= MG_GEO;
-          i1 = MMG5_inxt2[i];
-          i2 = MMG5_inxt2[i1];
-          mesh->point[ptr->v[i1]].tag |= MG_GEO;
-          mesh->point[ptr->v[i2]].tag |= MG_GEO;
+          mesh->point[ip1].tag |= MG_GEO;
+          mesh->point[ip2].tag |= MG_GEO;
           nr++;
         }
         if( intvalues[2*idx] > 2 ) { /* non-manifold edge */
-#warning remove MG_GEO fpr Mmg consistency ?
+#warning remove MG_GEO for consistency with Mmg ?
           /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
           ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
-          i1 = MMG5_inxt2[i];
-          i2 = MMG5_inxt2[i1];
-          mesh->point[ptr->v[i1]].tag |= MG_GEO + MG_NOM + MG_REF;
-          mesh->point[ptr->v[i2]].tag |= MG_GEO + MG_NOM + MG_REF;
+          mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
+          mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
           nm++;
         }
       }

From 0b1afdc37e66e2d93ae594c93f9752a7ab0311c3 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 10 Oct 2024 17:31:53 +0200
Subject: [PATCH 16/36] Fix erroneous point tag after point insertion along REF
 edge in level-set splitting.

---
 src/analys_pmmg.c | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 70448dbb..f2b13520 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2406,6 +2406,25 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
           mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
           nm++;
         }
+
+        /* If a feature edge has been provided it is possible that the edge tag
+         * has not been transferred to the edge extremities (for example a
+         * MG_REF edge along boudary triangles of same references may be split
+         * by the level-set. In this case, the new point has tag 0 and its tag
+         * has not yet been updated).
+         */
+        tag = mesh->point[ip1].tag;
+        mesh->point[ip1].tag |= ptr->tag[i];
+        // Remove the MG_NOSURF tag if the vertex is really required.
+        if ( (tag & MG_REQ) && !(tag & MG_NOSURF) ) {
+          mesh->point[ip1].tag &= ~MG_NOSURF;
+        }
+        tag = mesh->point[ip2].tag;
+        mesh->point[ip2].tag |= ptr->tag[i];
+        // Remove the MG_NOSURF tag if the vertex is really required.
+        if ( (tag & MG_REQ) && !(tag & MG_NOSURF) ) {
+          mesh->point[ip2].tag &= ~MG_NOSURF;
+        }
       }
     }
   }

From 30df811565fad6f1458b2140611f07a5fe61e42b Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 10 Oct 2024 17:35:44 +0200
Subject: [PATCH 17/36] Add an assert on the edge tag.

---
 src/analys_pmmg.c | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index f2b13520..555f21f3 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2363,6 +2363,10 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       if ( !MMG5_hGet( pHash, ip1, ip2, &edg, &tag ) ) continue;
       idx = edg-1;
 
+      assert ( (((ptr->tag[i] & MG_PARBDY) && !(ptr->tag[i] & MG_PARBDYBDY))
+                || (ptr->tag[i] & MG_BDY))
+               && "edge at intersection of BDY and PARBDY tria or along // tria" );
+
       /* Skip non-manifold edges */
 #ifndef NDEBUG
       if ( (ptr->tag[i] & MG_NOM) ) {
@@ -2371,7 +2375,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 #endif
 
       if( intvalues[2*idx] == 1 ) { /* no adjacent */
-#warning remove MG_GEO fpr Mmg consistency ?
+#warning remove MG_GEO for consistency with Mmg  ?
         /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the

From 18a3d8ac49a9d5a5f06400b56980f6c4f5599b60 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 11 Oct 2024 10:41:03 +0200
Subject: [PATCH 18/36] Remove erroneous assert: edges may lie along PARBDYBDY
 faces too.

---
 src/analys_pmmg.c | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 555f21f3..4bd1372c 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2363,10 +2363,6 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       if ( !MMG5_hGet( pHash, ip1, ip2, &edg, &tag ) ) continue;
       idx = edg-1;
 
-      assert ( (((ptr->tag[i] & MG_PARBDY) && !(ptr->tag[i] & MG_PARBDYBDY))
-                || (ptr->tag[i] & MG_BDY))
-               && "edge at intersection of BDY and PARBDY tria or along // tria" );
-
       /* Skip non-manifold edges */
 #ifndef NDEBUG
       if ( (ptr->tag[i] & MG_NOM) ) {

From fbb8f0c287bad8dbf1f44ae1d400868a4454ad78 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 11 Oct 2024 10:42:19 +0200
Subject: [PATCH 19/36] Remove obsolete warning: with all the performed tests,
 it seems that edge tags are consistents.

---
 src/analys_pmmg.c | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 4bd1372c..88a703f4 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2429,8 +2429,6 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
     }
   }
 
-#warning Luca: should we exchange tags on the edge communicator to update MG_PARBDYBDY?
-
   if ( abs(mesh->info.imprim) > 3 && nr > 0 )
     fprintf(stdout,"     %d ridges, %d edges updated\n",nr,ne);
 

From cc2c79755ef813ebd5c87d6ffd607dd326c3a4ad Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 17 Oct 2024 11:30:27 +0200
Subject: [PATCH 20/36] Update Mmg tag after second fix in swap23 (should solve
 lsCenIn-2 test case).

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e61958cb..f31a211b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 8a844888bfeaa9aff1136a7782ff561dbd08f40c
+    GIT_TAG  668f6752cb69d5f5256d52093735e8ff84188e7a
     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}

From b94434b51f32df2b45864fcc2e2d66b7f82f070a Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 22 Oct 2024 10:51:58 +0200
Subject: [PATCH 21/36] Refactoring: rename functions + update comments.

---
 src/analys_pmmg.c        | 20 +++++++++++---------
 src/communicators_pmmg.c |  2 +-
 src/hash_pmmg.c          | 14 +++++++++-----
 src/ls_pmmg.c            |  2 +-
 src/parmmg.h             |  4 ++--
 src/quality_pmmg.c       |  2 +-
 src/tag_pmmg.c           |  2 +-
 7 files changed, 26 insertions(+), 20 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 88a703f4..e29b9f1d 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2091,7 +2091,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
   memset(doublevalues,0x00,6*parmesh->int_edge_comm->nitem*sizeof(double));
 
 
-  /** Loop on boundary triangles and store a MG_REQ tag in the edge internal
+  /** Step 1: Loop on boundary triangles and store a MG_REQ tag in the edge internal
    *  communicator where the triangle touches a parallel edge.
    *  (Loop on all triangles, as the tags on corresponding edges are not
    *  required to match yet) */
@@ -2104,6 +2104,8 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 
       i1 = MMG5_inxt2[i];
       i2 = MMG5_inxt2[i1];
+      /* At this stage, the `tag` field store only the MG_PARBDY tag if pHash
+       * has been created by the `PMMG_hashPar_fromFaceComm` function */
       if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
       idx = edg-1;
 
@@ -2186,7 +2188,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
   }
 
 
-  /** Loop on true boundary triangles and store the normal in the edge internal
+  /** Step 2: Loop on true boundary triangles and store the normal in the edge internal
    *  communicator where the triangle touches a parallel edge. */
   for( k = 1; k <= mesh->nt; k++ ) {
     ptr = &mesh->tria[k];
@@ -2238,7 +2240,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
     }
   }
 
-  /** Exchange values on the interfaces among procs */
+  /** Step 3: Exchange values on the interfaces among procs */
   for ( k = 0; k < parmesh->next_edge_comm; ++k ) {
     ext_edge_comm = &parmesh->ext_edge_comm[k];
     nitem         = ext_edge_comm->nitem;
@@ -2271,7 +2273,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
                    comm,&status),return 0 );
   }
 
-  /** First pass: Increment the number of seen triangles, check for reference
+  /** Step 4 - First pass: Increment the number of seen triangles, check for reference
    *  edges and mark them with PMMG_UNSET, and store new triangles normals if
    *  there is room for them. */
   for ( k = 0; k < parmesh->next_edge_comm; ++k ) {
@@ -2315,7 +2317,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
   }
 
   if ( mesh->info.dhd > MMG5_ANGLIM ) {
-    /** Second pass: Check dihedral angle and mark geometric edge with
+    /** Step 4 - Second pass: Check dihedral angle and mark geometric edge with
      *  2*PMMG_UNSET (3x if it is already a reference edge) */
     for ( k = 0; k < parmesh->next_edge_comm; ++k ) {
       ext_edge_comm = &parmesh->ext_edge_comm[k];
@@ -2342,7 +2344,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
     }
   }
 
-  /** Third pass: Loop on triangles to tag edges and points.
+  /** Step 4 - Third pass: Loop on triangles to tag edges and points.
    *  Now we loop on all triangles because there could be parallel boundary
    *  edges not touched by triangles on the local process, but we want to add
    *  tags on them. */
@@ -2702,7 +2704,6 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-
   /* Set surface triangles to required in nosurf mode or for parallel boundaries */
   MMG3D_set_reqBoundaries(mesh);
 
@@ -2733,8 +2734,9 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* Hash parallel edges */
-  if( PMMG_hashPar_pmmg( parmesh,&hpar ) != PMMG_SUCCESS ) {
+  /* Hash parallel edges from tetra and face communicator: store edges and the
+     MG_PARBDY tag (other edge tags are not stored).*/
+  if( PMMG_hashPar_fromFaceComm( parmesh,&hpar ) != PMMG_SUCCESS ) {
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,hpar.geom);
diff --git a/src/communicators_pmmg.c b/src/communicators_pmmg.c
index 955f6b21..25bfa715 100644
--- a/src/communicators_pmmg.c
+++ b/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).
diff --git a/src/hash_pmmg.c b/src/hash_pmmg.c
index 63ed9424..8aabf04d 100644
--- a/src/hash_pmmg.c
+++ b/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) {
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 63946db8..3e4a8ee6 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -1630,7 +1630,7 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
 
   /* Hash parallel edges
      This step is needed to compute the edge communicator */
-  if( PMMG_hashPar_pmmg( parmesh,&hpar ) != PMMG_SUCCESS ) {
+  if( PMMG_hashPar_fromFaceComm( parmesh,&hpar ) != PMMG_SUCCESS ) {
     fprintf(stderr,"\n\n\n  -- WARNING: Impossible to compute the hash parallel edge \n\n\n");
     return 0;
   }
diff --git a/src/parmmg.h b/src/parmmg.h
index 57c6dd39..453939a0 100644
--- a/src/parmmg.h
+++ b/src/parmmg.h
@@ -450,8 +450,8 @@ int PMMG_boulen(PMMG_pParMesh parmesh,MMG5_pMesh mesh,int start,int ip,int iface
 int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh);
 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 );
-int PMMG_hashPar_pmmg( PMMG_pParMesh parmesh,MMG5_HGeom *pHash );
+int PMMG_hashParTag_fromXtet( MMG5_pMesh mesh,MMG5_HGeom *pHash );
+int PMMG_hashPar_fromFaceComm( PMMG_pParMesh parmesh,MMG5_HGeom *pHash );
 int PMMG_hashOldPar_pmmg( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_Hash *hash );
 
 /* Isovalue discretization functions */
diff --git a/src/quality_pmmg.c b/src/quality_pmmg.c
index 8bf4373a..f0a16b30 100644
--- a/src/quality_pmmg.c
+++ b/src/quality_pmmg.c
@@ -413,7 +413,7 @@ int PMMG_computePrilen( PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_pSol met, do
   }
 
   /* Hash parallel edges in the mesh */
-  if ( PMMG_hashPar(mesh,&hpar) != PMMG_SUCCESS ) return 0;
+  if ( PMMG_hashParTag_fromXtet(mesh,&hpar) != PMMG_SUCCESS ) return 0;
 
   /* Build parallel edge communicator */
   if( !PMMG_build_edgeComm( parmesh,mesh,&hpar,comm ) ) return 0;
diff --git a/src/tag_pmmg.c b/src/tag_pmmg.c
index a086e549..8980e9db 100644
--- a/src/tag_pmmg.c
+++ b/src/tag_pmmg.c
@@ -512,6 +512,7 @@ void PMMG_updateTagRef_node(PMMG_pParMesh parmesh, MMG5_pMesh mesh) {
  *
  * Check if faces on a parallel communicator connect elements with different
  * references, and tag them as a "true" boundary (thus PARBDYBDY).
+ *
  */
 int PMMG_parbdySet( PMMG_pParMesh parmesh ) {
   PMMG_pGrp      grp;
@@ -660,7 +661,6 @@ int PMMG_parbdySet( PMMG_pParMesh parmesh ) {
     }
   }
 
-
   /* Tag parallel points touched by simple MG_BDY faces as MG_PARBDYBDY
    * (a parallel surface can pinch a regular surface in just one point).
    * The same problem on edges is handled by MMG5_mmgHashTria. */

From c97fe457f787cd220ccc5b6a11a4777b5afb8195 Mon Sep 17 00:00:00 2001
From: AlgianeFroehly <algiane.froehly@inria.fr>
Date: Tue, 22 Oct 2024 11:25:34 +0200
Subject: [PATCH 22/36] Fix possible inconsistency between MG_REF edge tags.

  - If a mesh is subdivided into 3 partitions with 3 MPI processes, an
edge provided as ref edge at input by the user will miss the MG_REF tag
if it blongs to only purely parallel faces on a partition. The tag is
suitably set to REF on the other partitions, leading to non-consistency.

  - fixes ls ls-DisIn-toygeom-faces-5 test case
---
 src/analys_pmmg.c | 30 +++++++++++++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index e29b9f1d..e1999bb8 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2105,7 +2105,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       i1 = MMG5_inxt2[i];
       i2 = MMG5_inxt2[i1];
       /* At this stage, the `tag` field store only the MG_PARBDY tag if pHash
-       * has been created by the `PMMG_hashPar_fromFaceComm` function */
+       * has been created by the `PMMG_hashPar_fromFaceComm` function. */
       if ( !MMG5_hGet( pHash, ptr->v[i1], ptr->v[i2], &edg, &tag ) ) continue;
       idx = edg-1;
 
@@ -2115,6 +2115,24 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       if( (ptr->tag[i] & MG_REQ) && !(ptr->tag[i] & MG_NOSURF) ) {
         intvalues[idx] |= MG_REQ;
       }
+
+      /* Store other relevant edge tag that we want to synchronize through the
+       * parallel interfaces:
+         - MG_NOM, MG_OPNBDY and MG_GEO will be analyzed after (this analysis
+           should be consistent through the procs);
+         - MG_REQ and MG_NOSURF tags are already dealed too;
+         - OLDPARBDY, PARBDY, PARBDYBDY and OVERLAP are related to ParMmg and
+           should be consistent.
+
+         It left us with:
+           - the MG_REF tag that may be not consistent (if, on a
+             partition, the edge belongs to only PARBDY faces (non PARBDYBDY), it is
+             not marked as REF, while it may be marked as ref if it belongs to a
+             true boundary and is provided as a user ref edge between triangles with
+             same references on another partition.)
+      */
+      tag = ptr->tag[i] & MG_REF;
+      intvalues[idx] |= tag;
     }
   }
 
@@ -2152,6 +2170,11 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
         intvalues[idx] |= MG_REQ;
         intvalues[idx] &= ~MG_NOSURF;
       }
+
+      if ( itorecv[i] & MG_REF ) {
+        /* Sync ref tags */
+        intvalues[idx] |= MG_REF;
+      }
     }
   }
 
@@ -2175,6 +2198,11 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
         ptr->tag[i] |= MG_REQ;
         ptr->tag[i] &= ~MG_NOSURF;
       }
+
+      if ( intvalues[idx] & MG_REF ) {
+        ptr->tag[i] |= MG_REF;
+      }
+
     }
   }
 

From ddfc4da323db4dc72e768a19a8fdb7f1b83e19e6 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 22 Oct 2024 16:44:20 +0200
Subject: [PATCH 23/36] Update Mmg tag to not check PARBDYBDY consistency +
 update doc + add a consistency check at end of analysis.

---
 CMakeLists.txt    |  2 +-
 src/analys_pmmg.c | 38 +++++++++++++++++++++++++++++++-------
 src/tag_pmmg.c    | 13 +++++++++++--
 3 files changed, 43 insertions(+), 10 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f31a211b..e12f5b23 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG  668f6752cb69d5f5256d52093735e8ff84188e7a
+    GIT_TAG 86de8073a91e3d1b03ea3bc6a839db2b9fd1d6e3
     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}
diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index e1999bb8..69e6143f 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -34,6 +34,7 @@
 
 #include "parmmg.h"
 #include "libmmg3d.h"
+#include "mmgexterns_private.h"
 
 /**
  * \param ppt pointer toward the point structure
@@ -2121,15 +2122,27 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
          - MG_NOM, MG_OPNBDY and MG_GEO will be analyzed after (this analysis
            should be consistent through the procs);
          - MG_REQ and MG_NOSURF tags are already dealed too;
-         - OLDPARBDY, PARBDY, PARBDYBDY and OVERLAP are related to ParMmg and
+         - OLDPARBDY, PARBDY and OVERLAP are related to ParMmg and
            should be consistent.
 
          It left us with:
            - the MG_REF tag that may be not consistent (if, on a
-             partition, the edge belongs to only PARBDY faces (non PARBDYBDY), it is
+             partition, the edge belongs to only PARBDY faces (non PARBDYBDY),
+             it is
              not marked as REF, while it may be marked as ref if it belongs to a
-             true boundary and is provided as a user ref edge between triangles with
-             same references on another partition.)
+             true boundary and is provided as a user ref edge between triangles
+             with same references on another partition.)
+           - the MG_PARBDYBDY tag that may be inconsistent between trias due to
+             the call of PMMG_parbdyTria: for a physical  boundary triangle
+             at partition interface. On the domain with lower ref, the
+             PARBDYPARBDY tag is removed from edges (to ensure the tria
+             orientation during analysis, the tria will be looked from one rank
+             only, the rank that owned the domain with higer ref). If the edge
+             also belongs to another boundary triangle, it will still have
+             the PARBDYBDY tag on this triangle.
+             We don't want to synchronize this specific tag as it is used for
+             the parallel analysis (for example in hashNorver_loop to loop
+             on well oriented true boundary faces).
       */
       tag = ptr->tag[i] & MG_REF;
       intvalues[idx] |= tag;
@@ -2202,7 +2215,6 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
       if ( intvalues[idx] & MG_REF ) {
         ptr->tag[i] |= MG_REF;
       }
-
     }
   }
 
@@ -2999,7 +3011,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 #warning Luca: check that parbdy are skipped
   MMG5_chkVertexConnectedDomains(mesh);
 
-  /* build hash table for geometric edges */
+  /* build hash table for geometric edges: gather tag infos from edges and
+   * triangles and store these infos in tria. Skip non PARBDYBDY // edges. */
   if ( !mesh->na && !MMG5_hGeom(mesh) ) {
     fprintf(stderr,"\n  ## Hashing problem (0). Exit program.\n");
     PMMG_edge_comm_free( parmesh );
@@ -3013,7 +3026,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* Update edges tags and references for xtetras */
+  /* Update MG_REQ and MG_NOSURF edges tags as well as references for xtetras */
   if ( !MMG5_bdryUpdate(mesh) ) {
     fprintf(stderr,"\n  ## Boundary problem. Exit program.\n");
     PMMG_edge_comm_free( parmesh );
@@ -3063,5 +3076,16 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 
   if ( mesh->nprism ) MMG5_DEL_MEM(mesh,mesh->adjapr);
 
+#ifndef NDEBUG
+  MMG5_int 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;
+    }
+  }
+#endif
+
   return 1;
 }
diff --git a/src/tag_pmmg.c b/src/tag_pmmg.c
index 8980e9db..8938b603 100644
--- a/src/tag_pmmg.c
+++ b/src/tag_pmmg.c
@@ -513,6 +513,11 @@ void PMMG_updateTagRef_node(PMMG_pParMesh parmesh, MMG5_pMesh mesh) {
  * Check if faces on a parallel communicator connect elements with different
  * references, and tag them as a "true" boundary (thus PARBDYBDY).
  *
+ * \remark: Edge tags are not updated along faces with the PARBDYBDY tag as
+ * it is not sufficient to maintain the consistency of PARBDYBDY tags through the mesh.
+ * Morover, even if we manage to have consistent tags inside one mesh, we will
+ * still have to synchronize the edge tags through the partition interfaces. In consequence,
+ * the PARBDYBDY tags may be not consistent throught the entire remeshing process.
  */
 int PMMG_parbdySet( PMMG_pParMesh parmesh ) {
   PMMG_pGrp      grp;
@@ -635,8 +640,8 @@ int PMMG_parbdySet( PMMG_pParMesh parmesh ) {
       if( seenFace[idx] != 1 ) continue;
 
       /* Tag face as "true" boundary if its ref is different (or if triangle has
-       * a non nul ref in openbdy mode), delete reference if it is only a
-       * parallel boundary */
+       * a non nul ref in opnbdy mode), delete reference if it is only a
+       * parallel boundary. */
       if ( (mesh->info.opnbdy && pxt->ref[ifac]>0) ) {
         pxt->ftag[ifac] |= MG_PARBDYBDY;
         for( i = 0; i < 3; i++)
@@ -835,6 +840,10 @@ int PMMG_parbdyTria( PMMG_pParMesh parmesh ) {
         ptt->v[1] = ib;
         ptt->v[2] = ic;
       } else {
+        /* The boundary will belong to one partition only, for the other
+         * partition, the triangle will be considered as simply parallel.  This
+         * means that the PARBDYBDY triangle tags may be not consistent from
+         * here.  */
         ptt->tag[0] &= ~MG_PARBDYBDY;
         ptt->tag[1] &= ~MG_PARBDYBDY;
         ptt->tag[2] &= ~MG_PARBDYBDY;

From 35ddd83901ba8bbf1259c3dc6f71be276f3d38de Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Thu, 24 Oct 2024 18:33:56 +0200
Subject: [PATCH 24/36] Fix division by 0 when computing the tangent at EDG
 point (issue with EDG tag at point along non EDG edge due to spurious input
 surface triangles matching the // interface).

---
 src/analys_pmmg.c | 22 +++++++++++++++++++++-
 src/ls_pmmg.c     |  6 +++++-
 2 files changed, 26 insertions(+), 2 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 69e6143f..53db5d70 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -953,6 +953,26 @@ int PMMG_hashNorver_normals( PMMG_pParMesh parmesh, PMMG_hn_loopvar *var,MPI_Com
 #warning Luca: why not like in MMG5_boulec?
       if( MG_EDG(var->ppt->tag) ) {
 
+        if ( (!intvalues[2*idx]) || (!intvalues[2*idx+1]) ) {
+           /* Issue with input edges along parallel interfaces: an input REF or
+            * GEO edge may be stored inside a boundary triangle at interface
+            * between tetra with same reference. In this case, the LS split,
+            * creates a points with REF or GEO tag to match the tag of the
+            * triangle edge (in setfeatures) but the triangle and associated
+            * edge tag are not stored in the xtetra, ended with a division by 0
+            * when computing the tangent at point. */
+
+          if ( parmesh->ddebug  ) {
+            printf("  ## Warning: %s:%d: rank %d: tag inconsistency: ppt tag %u"
+                   " - edge extremities %d %d\n              Point tag is removed.\n",
+                   __func__,__LINE__,parmesh->myrank,
+                   var->ppt->tag,intvalues[2*idx],intvalues[2*idx+1]);
+          }
+
+          var->ppt->tag &= ( (~MG_REF) && (~MG_GEO) );
+          continue;
+        }
+
         c[0] = &doublevalues[6*idx];
         c[1] = &doublevalues[6*idx+3];
 
@@ -2451,7 +2471,7 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 
         /* If a feature edge has been provided it is possible that the edge tag
          * has not been transferred to the edge extremities (for example a
-         * MG_REF edge along boudary triangles of same references may be split
+         * MG_REF edge along boundary triangles of same references may be split
          * by the level-set. In this case, the new point has tag 0 and its tag
          * has not yet been updated).
          */
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 3e4a8ee6..cb46ea33 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -1583,7 +1583,11 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
   }
 
   /* Identify surface mesh. Clean triangle array: remove useless or double
-     triangles and add the missing ones */
+     triangles and add the missing ones.  Remark: spurious boundary triangles
+     across parallel interface cannot be removed by the serial function but will
+     not be stored inthe xtetra by the MMG5_bdrySet function during analysis.
+     This may create inconsistencies between edge and point tags.
+  */
   if ( !MMG5_chkBdryTria(mesh) ) {
     fprintf(stderr,"\n  ## Boundary problem. Exit program.\n");
     return 0;

From ba27447c79486cb72eda9436e71952269834455e Mon Sep 17 00:00:00 2001
From: AlgianeFroehly <algiane.froehly@inria.fr>
Date: Fri, 25 Oct 2024 10:23:33 +0200
Subject: [PATCH 25/36] Fix previous commit: intvalue array is erased at
 function beginning so cannot be used.

---
 src/analys_pmmg.c | 61 +++++++++++++++++++++++++++++++----------------
 1 file changed, 41 insertions(+), 20 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 53db5d70..bd209ca5 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -356,7 +356,11 @@ int PMMG_hashNorver_edges( PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
           pos++;
         }
         assert(found < 2);
+
+#ifndef NDEBUG
         if( pos == 2 ) assert(found);
+#endif
+
         if( !found ) {
           assert( pos < 2 );
           intvalues[2*idx+pos] = gip;
@@ -932,6 +936,43 @@ int PMMG_hashNorver_normals( PMMG_pParMesh parmesh, PMMG_hn_loopvar *var,MPI_Com
   intvalues    = parmesh->int_node_comm->intvalues;
   doublevalues = parmesh->int_node_comm->doublevalues;
 
+  /* Check consistency between point tag and stored infos before erasing
+   * intvalues: if MG_EDG edges have been found on both side of a non singular
+   * MG_EDG point, the intvalues array should be filled with the global id of
+   * the 2 neigbouring points along the feature edge.
+   * It aims to solve the
+   * following issue with input edges along parallel interfaces: an input REF or
+   * GEO edge may be stored inside a boundary triangle at interface * between
+   * tetra with same reference. In this case, the LS split, * creates a points
+   * with REF or GEO tag to match the tag of the * triangle edge (in
+   * setfeatures) but the triangle and associated * edge tag are not stored in
+   * the xtetra, ended with a division by 0 * when computing the tangent at
+   * point. */
+  for( var->ip = 1; var->ip <= var->mesh->np; var->ip++ ) {
+    var->ppt = &var->mesh->point[var->ip];
+
+    /* Loop on parallel, non-singular points (they have been flagged in
+     * PMMG_hashNorver_xp_init()) */
+    if( var->ppt->flag && MG_EDG(var->ppt->tag) ) {
+
+      idx = PMMG_point2int_comm_index_get( var->ppt );
+
+      if ( (!intvalues[2*idx]) || (!intvalues[2*idx+1]) ) {
+        /* We will miss infos to compute the tangent... Erase point tag */
+
+        if ( parmesh->ddebug  ) {
+          printf("  ## Warning: %s:%d: rank %d: tag inconsistency: ppt %d tag %u"
+                 " - edge extremities %d %d\n              Point tag is removed.\n",
+                 __func__,__LINE__,parmesh->myrank,var->ip,
+                 var->ppt->tag,intvalues[2*idx],intvalues[2*idx+1]);
+        }
+
+        var->ppt->tag &= ( (~MG_REF) & (~MG_GEO) );
+      }
+    }
+  }
+
+
   memset(intvalues,0,parmesh->int_node_comm->nitem*sizeof(int));
 
   /* Accumulate normal vector contributions */
@@ -953,26 +994,6 @@ int PMMG_hashNorver_normals( PMMG_pParMesh parmesh, PMMG_hn_loopvar *var,MPI_Com
 #warning Luca: why not like in MMG5_boulec?
       if( MG_EDG(var->ppt->tag) ) {
 
-        if ( (!intvalues[2*idx]) || (!intvalues[2*idx+1]) ) {
-           /* Issue with input edges along parallel interfaces: an input REF or
-            * GEO edge may be stored inside a boundary triangle at interface
-            * between tetra with same reference. In this case, the LS split,
-            * creates a points with REF or GEO tag to match the tag of the
-            * triangle edge (in setfeatures) but the triangle and associated
-            * edge tag are not stored in the xtetra, ended with a division by 0
-            * when computing the tangent at point. */
-
-          if ( parmesh->ddebug  ) {
-            printf("  ## Warning: %s:%d: rank %d: tag inconsistency: ppt tag %u"
-                   " - edge extremities %d %d\n              Point tag is removed.\n",
-                   __func__,__LINE__,parmesh->myrank,
-                   var->ppt->tag,intvalues[2*idx],intvalues[2*idx+1]);
-          }
-
-          var->ppt->tag &= ( (~MG_REF) && (~MG_GEO) );
-          continue;
-        }
-
         c[0] = &doublevalues[6*idx];
         c[1] = &doublevalues[6*idx+3];
 

From cdf19dce660170d054f4be3ffd94c93a2af9f3a7 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 25 Oct 2024 15:21:35 +0200
Subject: [PATCH 26/36] Reset ls references before analysing parallel triangles
 in order to not ad PARBDYBDY tags on temporary triangles that are finally
 deleted when storing the tags inside the xtetra.

---
 src/analys_pmmg.c |  8 ++++++--
 src/ls_pmmg.c     | 22 +++++++++++++++-------
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index bd209ca5..52dfb4ca 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2795,7 +2795,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* build hash table for geometric edges */
+  /* build hash table for geometric edges: gather tag infos from edges and
+   * triangles and store these infos in tria. Skip non PARBDYBDY // edges. */
   if ( !MMG5_hGeom(mesh) ) {
     fprintf(stderr,"\n  ## Hashing problem (0). Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
@@ -2929,7 +2930,10 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
     return 0;
   }
 
-  /* Tag parallel faces on material interfaces as boundary (ie, add \ref MG_PARBDYBDY tag) */
+  /* Tag parallel faces on material interfaces as boundary (ie, add \ref
+   * MG_PARBDYBDY tag) and remove spurious boundary MG_PARBDYBDY tags coming
+   * from spurious internal triangles (between tetras of same references) along
+   * partition interfaces.  */
   if( !PMMG_parbdySet( parmesh ) ) {
     fprintf(stderr,"\n  ## Unable to recognize parallel faces on material interfaces. Exit program.\n");
     MMG5_DEL_MEM(mesh,hash.item);
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index cb46ea33..13fa3655 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -1576,6 +1576,19 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
     return 0;
   }
 
+  /* Reset the mesh->info.isoref field everywhere */
+  if ( !MMG3D_resetRef_ls(mesh) ) {
+    fprintf(stderr,"\n  ## Problem in resetting references. Exit program.\n");
+    return 0;
+  }
+
+  /* Tag parallel triangles on material interfaces as boundary. */
+  if( !PMMG_parbdyTria( parmesh ) ) {
+    fprintf(stderr,"\n  ## Unable to recognize parallel triangles on material interfaces."
+            " Exit program.\n");
+    return 0;
+  }
+
   /* Check the compatibility of triangle orientation with tetra faces */
   if ( !MMG5_bdryPerm(mesh) ) {
     fprintf(stderr,"\n  ## Boundary orientation problem. Exit program.\n");
@@ -1593,7 +1606,8 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
     return 0;
   }
 
-  /* Build hash table for initial edges */
+  /* Build hash table for initial edges: gather tag infos from edges and
+   * triangles and store these infos in tria. Skip non PARBDYBDY // edges. */
   if ( !MMG5_hGeom(mesh) ) {
     fprintf(stderr,"\n  ## Hashing problem (0). Exit program.\n");
     return 0;
@@ -1605,12 +1619,6 @@ int PMMG_ls(PMMG_pParMesh parmesh, MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met)
     return 0;
   }
 
-  /* Reset the mesh->info.isoref field everywhere */
-  if ( !MMG3D_resetRef_ls(mesh) ) {
-    fprintf(stderr,"\n  ## Problem in resetting references. Exit program.\n");
-    return 0;
-  }
-
   /** \todo TODO :: Removal of small parasitic components */
   if ( mesh->info.rmc > 0 ) {
     PMMG_rmc(parmesh,mesh,sol);

From 0e467d35bbf4aa92e420c560e9bc61702d2a8fee Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 25 Oct 2024 15:22:02 +0200
Subject: [PATCH 27/36] Update Mmg tag.

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e12f5b23..0831f45f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 86de8073a91e3d1b03ea3bc6a839db2b9fd1d6e3
+    GIT_TAG 04aad1fdf2d0d960321a983301744e94276a1a61
     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}

From cd157e1e43bc90ed5c863936a659794f4286272e Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 25 Oct 2024 15:27:56 +0200
Subject: [PATCH 28/36] Add support for empty group in parbdyTria.

---
 src/tag_pmmg.c | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/tag_pmmg.c b/src/tag_pmmg.c
index 8938b603..f36131b0 100644
--- a/src/tag_pmmg.c
+++ b/src/tag_pmmg.c
@@ -721,7 +721,13 @@ int PMMG_parbdyTria( PMMG_pParMesh parmesh ) {
 
   comm   = parmesh->comm;
   myrank = parmesh->myrank;
-  assert( parmesh->ngrp == 1 );
+
+  if ( !parmesh->ngrp ) {
+    return 1;
+  }
+
+  assert( parmesh->ngrp == 1 && "Not implemented for multiple groups");
+
   mesh = grp->mesh;
 
   /* intvalues will be used to store tetra ref */

From fe69e6071c9f8cf60c36368a9b4ff93b1eaaa9ec Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 25 Oct 2024 15:29:20 +0200
Subject: [PATCH 29/36] Add checks for mesh consistency in CMake Debug mode.

---
 src/libparmmg1.c         | 10 ++++++++++
 src/loadbalancing_pmmg.c | 34 ++++++++++++++++++++++++++++++++--
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 0f6ae98d..44478505 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -567,6 +567,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 ) {
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index e27eafca..28094f1e 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -86,18 +86,29 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   for ( igrp=0; igrp < parmesh->ngrp; igrp++ )
     ne += parmesh->listgrp[igrp].mesh->ne;
 
+  int k;
+#ifndef NDEBUG
+  for (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 ) {
     /** Split the ngrp groups of listgrp into a higher number of groups */
     ier = PMMG_split_n2mGrps(parmesh,PMMG_GRPSPL_DISTR_TARGET,1,partitioning_mode);
   }
 
-  for (int k=0; k<parmesh->ngrp; ++k ) {
+#ifndef NDEBUG
+  for ( k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
-
+#endif
 
   /* There is mpi comms in distribute_grps thus we don't want that one proc
    * enters the function and not the other proc */
@@ -119,6 +130,16 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
     chrono(ON,&(ctim[tim]));
   }
 
+#ifndef NDEBUG
+  int 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;
+    }
+  }
+#endif
+
   ier = PMMG_distribute_grps(parmesh,partitioning_mode);
 
   if ( ier <= 0 ) {
@@ -142,6 +163,15 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
     chrono(ON,&(ctim[tim]));
   }
 
+#ifndef NDEBUG
+  for ( int 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 ( ier ) {
     /** Redistribute the ngrp groups of listgrp into a lower number of groups */
     ier = PMMG_split_n2mGrps(parmesh,PMMG_GRPSPL_MMG_TARGET,0,partitioning_mode);

From 12b0390319e75bc9c3a293fb14d1221e64dc6f07 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 16:58:08 +0100
Subject: [PATCH 30/36] Update git tag to avoid issue with wrong normal along
 nom edge with 2 singular extremities.

---
 CMakeLists.txt           |  2 +-
 src/libparmmg.c          | 13 +++++++++++++
 src/libparmmg1.c         | 28 +++++++++++++++++++++++++++-
 src/loadbalancing_pmmg.c |  3 +++
 src/mergemesh_pmmg.c     | 22 ++++++++++++++++++++++
 5 files changed, 66 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0831f45f..f5e0f6da 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 04aad1fdf2d0d960321a983301744e94276a1a61
+    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}
diff --git a/src/libparmmg.c b/src/libparmmg.c
index 9fd6201a..c38e34a1 100644
--- a/src/libparmmg.c
+++ b/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
  *
@@ -408,6 +410,17 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
       return PMMG_STRONGFAILURE;
     }
 
+     const char* s = getenv("SAVE");
+
+     printf("ngrps ? %d\n",parmesh->ngrp);
+
+     for ( MMG5_int i=0; i<parmesh->ngrp; ++i ) {
+
+       if ( s ) {
+         PMMG_grp_to_saveMesh( parmesh, i, "AfterLs-br" );
+       }
+     }
+
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);
     if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 44478505..e9f4e672 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -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
@@ -568,6 +569,16 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   /** Set inputMet flag */
   parmesh->info.inputMet = 0;
 
+  const char* s = getenv("SAVE");
+
+  for ( i=0; i<parmesh->ngrp; ++i ) {
+
+     if ( s ) {
+      PMMG_grp_to_saveMesh( parmesh, i, "BeforeSplit" );
+    }
+  }
+
+
 #ifndef NDEBUG
   for ( i=0; i<parmesh->ngrp; ++i ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {
@@ -632,6 +643,8 @@ 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) ) {
@@ -687,10 +700,14 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
         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 */
@@ -761,6 +778,10 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           }
         }
 
+        if ( s ) {
+          PMMG_grp_to_saveMesh( parmesh, i, "BeforeAdp" );
+        }
+
 #ifdef PATTERN
         ier = MMG5_mmg3d1_pattern( mesh, met, permNodGlob );
 #else
@@ -769,6 +790,11 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
         mesh->npi = mesh->np;
         mesh->nei = mesh->ne;
 
+
+        if ( s ) {
+          PMMG_grp_to_saveMesh( parmesh, i, "AfterAdp" );
+        }
+
         if ( !ier ) {
           fprintf(stderr,"\n  ## MMG remeshing problem. Exit program.\n");
         }
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index 28094f1e..25ec6190 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -44,6 +44,8 @@
  * 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;
@@ -83,6 +85,7 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   }
 
   ne = 0;
+  printf("ngrps before loadbal = %d\n", parmesh->ngrp);
   for ( igrp=0; igrp < parmesh->ngrp; igrp++ )
     ne += parmesh->listgrp[igrp].mesh->ne;
 
diff --git a/src/mergemesh_pmmg.c b/src/mergemesh_pmmg.c
index ff450f20..2b3e3e8a 100644
--- a/src/mergemesh_pmmg.c
+++ b/src/mergemesh_pmmg.c
@@ -456,6 +456,9 @@ int PMMG_mergeGrps_interfacePoints( PMMG_pParMesh parmesh ) {
            && ( poi_id_glo < parmesh->int_node_comm->nitem )
            && "check intvalues indices" );
     intvalues[ poi_id_glo ] = poi_id_int;
+    if ( poi_id_int == 1 ) {
+      printf("point id glob is %d\n",poi_id_glo);
+    }
   }
 
   /** Step 2: add points referenced by the rest of the groups meshes' internal
@@ -676,6 +679,15 @@ int PMMG_mergeGrpJinI_interfaceTetra(PMMG_pParMesh parmesh,PMMG_pGrp grpI,
       }
       else break;
     }
+/*
+    if ( ie==112310 ) {
+      printf("merge interface tetra %d: %d %d %d %d - %d\n",ie,ptI->v[0],ptI->v[1],ptI->v[2],ptI->v[3],ptI->xt);
+      printf("xt %d %d %d %d - %d %d %d %d %d %d\n",pxtI->ftag[0],pxtI->ftag[1],pxtI->ftag[2],
+             pxtI->ftag[3],pxtI->tag[0],pxtI->tag[1],pxtI->tag[2],
+             pxtI->tag[3],pxtI->tag[4],pxtI->tag[5]);
+      printf("from tetra: %d: %d %d %d %d - ref: %d - xt: %d\n",iel,ptJ->v[0],ptJ->v[1],ptJ->v[2],ptJ->v[3],ptJ->ref,ptJ->xt);
+    }
+*/
   }
 
   return ier;
@@ -745,6 +757,16 @@ int PMMG_mergeGrpJinI_internalTetra( PMMG_pGrp grpI, PMMG_pGrp grpJ ) {
       }
       else break;
     }
+
+/*
+    if ( ie==112310 ) {
+      printf("merge tetra %d: %d %d %d %d - %d\n",ie,ptI->v[0],ptI->v[1],ptI->v[2],ptI->v[3],ptI->xt);
+      printf("xt %d %d %d %d - %d %d %d %d %d %d\n",pxtI->ftag[0],pxtI->ftag[1],pxtI->ftag[2],
+             pxtI->ftag[3],pxtI->tag[0],pxtI->tag[1],pxtI->tag[2],
+             pxtI->tag[3],pxtI->tag[4],pxtI->tag[5]);
+      printf("from tetra: %d: %d %d %d %d - ref: %d - xt: %d\n",k,ptJ->v[0],ptJ->v[1],ptJ->v[2],ptJ->v[3],ptJ->ref,ptJ->xt);
+    }
+*/
   }
 
   return ier;

From b9a3cc809cb9a6d132d80e8206c4174894c62432 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 16:58:42 +0100
Subject: [PATCH 31/36] Check for truely required edges when ensuring edge
 consistency after group splitting.

---
 src/grpsplit_pmmg.c | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index baf3cb9c..9362402c 100644
--- a/src/grpsplit_pmmg.c
+++ b/src/grpsplit_pmmg.c
@@ -1044,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);
@@ -1119,6 +1119,11 @@ PMMG_splitGrps_fillGroup( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int ngrp,int g
        * 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;
+      }
     }
 
   }

From fdabace9312ad7de289f081b1babf4cc979b4169 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 17:38:48 +0100
Subject: [PATCH 32/36] Remove MG_GEO marker from MG_NOM edges in // analysis
 (in Mmg, nom edges are not marked as ridges anymore).

---
 src/analys_pmmg.c | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 55585b70..cfd7dc71 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2345,14 +2345,13 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 #endif
 
       if( intvalues[2*idx] == 1 ) { /* no adjacent */
-#warning remove MG_GEO for consistency with Mmg  ?
         /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
-        ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
-        mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
-        mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
+        ptr->tag[i] |=  MG_NOM + MG_REF;
+        mesh->point[ip1].tag |= MG_NOM + MG_REF;
+        mesh->point[ip2].tag |= MG_NOM + MG_REF;
         nr++;
       } else {
         if( (intvalues[2*idx+1] ==   PMMG_UNSET) ||
@@ -2370,14 +2369,13 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
           nr++;
         }
         if( intvalues[2*idx] > 2 ) { /* non-manifold edge */
-#warning remove MG_GEO for consistency with Mmg ?
           /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
-          ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
-          mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
-          mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
+          ptr->tag[i] |= MG_NOM + MG_REF;
+          mesh->point[ip1].tag |= MG_NOM + MG_REF;
+          mesh->point[ip2].tag |= MG_NOM + MG_REF;
           nm++;
         }
 

From 5a780dbfaaf3fa230b8f65774d06a70f633de0ac Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 18:43:22 +0100
Subject: [PATCH 33/36] Revert "Remove MG_GEO marker from MG_NOM edges in //
 analysis (in Mmg, nom edges are not marked as ridges anymore)."

This reverts commit fdabace9312ad7de289f081b1babf4cc979b4169: for now, we tag non-consistency if we don't add the MG_GEO tag at the same time than the MG_NOM one. I have no time to look at this now... postponed
---
 src/analys_pmmg.c | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index cfd7dc71..de50bd60 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -1993,6 +1993,9 @@ int PMMG_singul(PMMG_pParMesh parmesh,MMG5_pMesh mesh,PMMG_hn_loopvar *var,MPI_C
  * Boundary triangles shared between two processes have been tagged as
  * MG_PARBDYBDY only on the process who has them with the right orientation
  * (by PMMG_parbdyTria), so they will be processed only once.
+ *
+ * \todo Do not add MG_GEO tag to MG_NOM edges and fix tag non-consistencies in this case.
+ *
  */
 int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm comm ) {
   PMMG_pGrp      grp;
@@ -2345,13 +2348,14 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
 #endif
 
       if( intvalues[2*idx] == 1 ) { /* no adjacent */
+#warning remove MG_GEO for consistency with Mmg  ?
         /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
-        ptr->tag[i] |=  MG_NOM + MG_REF;
-        mesh->point[ip1].tag |= MG_NOM + MG_REF;
-        mesh->point[ip2].tag |= MG_NOM + MG_REF;
+        ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
+        mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
+        mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
         nr++;
       } else {
         if( (intvalues[2*idx+1] ==   PMMG_UNSET) ||
@@ -2369,13 +2373,14 @@ int PMMG_setfeatures(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI
           nr++;
         }
         if( intvalues[2*idx] > 2 ) { /* non-manifold edge */
+#warning remove MG_GEO for consistency with Mmg ?
           /* MG_REF info is not analyzed in parallel for non-manifold edges (only
            serially by Mmy).  As we need to ensure the tag consistency across
            the processes and for sake of simplicity, we simply mark all the
            MG_NOM edges as MG_REF */
-          ptr->tag[i] |= MG_NOM + MG_REF;
-          mesh->point[ip1].tag |= MG_NOM + MG_REF;
-          mesh->point[ip2].tag |= MG_NOM + MG_REF;
+          ptr->tag[i] |= MG_GEO + MG_NOM + MG_REF;
+          mesh->point[ip1].tag |= MG_GEO + MG_NOM + MG_REF;
+          mesh->point[ip2].tag |= MG_GEO + MG_NOM + MG_REF;
           nm++;
         }
 

From d7b2439a5cb2d2e4eab74271eee94c1a138df796 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 22:10:37 +0100
Subject: [PATCH 34/36] Remove debug prints.

---
 src/loadbalancing_pmmg.c | 1 -
 src/mergemesh_pmmg.c     | 3 ---
 2 files changed, 4 deletions(-)

diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index 25ec6190..8ba8b3c3 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -85,7 +85,6 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   }
 
   ne = 0;
-  printf("ngrps before loadbal = %d\n", parmesh->ngrp);
   for ( igrp=0; igrp < parmesh->ngrp; igrp++ )
     ne += parmesh->listgrp[igrp].mesh->ne;
 
diff --git a/src/mergemesh_pmmg.c b/src/mergemesh_pmmg.c
index 2b3e3e8a..010fdf36 100644
--- a/src/mergemesh_pmmg.c
+++ b/src/mergemesh_pmmg.c
@@ -456,9 +456,6 @@ int PMMG_mergeGrps_interfacePoints( PMMG_pParMesh parmesh ) {
            && ( poi_id_glo < parmesh->int_node_comm->nitem )
            && "check intvalues indices" );
     intvalues[ poi_id_glo ] = poi_id_int;
-    if ( poi_id_int == 1 ) {
-      printf("point id glob is %d\n",poi_id_glo);
-    }
   }
 
   /** Step 2: add points referenced by the rest of the groups meshes' internal

From fc412d0d9c84f8f0454274a2e99c5704f8114a98 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 22:19:38 +0100
Subject: [PATCH 35/36] Remove debug prints.

---
 src/libparmmg.c  | 11 -----------
 src/libparmmg1.c | 10 ----------
 2 files changed, 21 deletions(-)

diff --git a/src/libparmmg.c b/src/libparmmg.c
index c38e34a1..401fdce6 100644
--- a/src/libparmmg.c
+++ b/src/libparmmg.c
@@ -410,17 +410,6 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
       return PMMG_STRONGFAILURE;
     }
 
-     const char* s = getenv("SAVE");
-
-     printf("ngrps ? %d\n",parmesh->ngrp);
-
-     for ( MMG5_int i=0; i<parmesh->ngrp; ++i ) {
-
-       if ( s ) {
-         PMMG_grp_to_saveMesh( parmesh, i, "AfterLs-br" );
-       }
-     }
-
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);
     if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index e9f4e672..c34b2be9 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -569,16 +569,6 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   /** Set inputMet flag */
   parmesh->info.inputMet = 0;
 
-  const char* s = getenv("SAVE");
-
-  for ( i=0; i<parmesh->ngrp; ++i ) {
-
-     if ( s ) {
-      PMMG_grp_to_saveMesh( parmesh, i, "BeforeSplit" );
-    }
-  }
-
-
 #ifndef NDEBUG
   for ( i=0; i<parmesh->ngrp; ++i ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {

From 6063bc0cbc2003b47c62b4692bf73637f60f3bdd Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Mon, 28 Oct 2024 22:22:31 +0100
Subject: [PATCH 36/36] Remove debug prints.

---
 src/libparmmg1.c | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index c34b2be9..d4700387 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -768,10 +768,6 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           }
         }
 
-        if ( s ) {
-          PMMG_grp_to_saveMesh( parmesh, i, "BeforeAdp" );
-        }
-
 #ifdef PATTERN
         ier = MMG5_mmg3d1_pattern( mesh, met, permNodGlob );
 #else
@@ -780,11 +776,6 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
         mesh->npi = mesh->np;
         mesh->nei = mesh->ne;
 
-
-        if ( s ) {
-          PMMG_grp_to_saveMesh( parmesh, i, "AfterAdp" );
-        }
-
         if ( !ier ) {
           fprintf(stderr,"\n  ## MMG remeshing problem. Exit program.\n");
         }