From 4399a851339d55e40f88945a13af41f06c2c2554 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 27 Sep 2024 15:57:00 +0200
Subject: [PATCH 01/11] Small refactoring: should not change the results.

---
 src/libparmmg.c  |  7 ++++++-
 src/libparmmg1.c |  3 +--
 src/ls_pmmg.c    | 48 +++++++++++++++++++++++++++++++-----------------
 3 files changed, 38 insertions(+), 20 deletions(-)

diff --git a/src/libparmmg.c b/src/libparmmg.c
index 9fd6201a..76a05f88 100644
--- a/src/libparmmg.c
+++ b/src/libparmmg.c
@@ -1215,7 +1215,12 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
                  PMMG_DEL_MEM(parmesh,status,MPI_Status,"mpi_status");
                  PMMG_destroy_int(parmesh,ptr_int,nptr,"vertGlobNum");
                  return 0 );
-      for( i = 0; i < nitem; i++ ) assert(itorecv[i]);
+
+#ifndef NDEBUG
+      for( i = 0; i < nitem; i++ ) {
+        assert(itorecv[i]);
+      }
+#endif
     }
   }
 
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 16d4b22f..28459eb0 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -189,8 +189,7 @@ int PMMG_packTetra( PMMG_pParMesh parmesh, int igrp ) {
  *
  * \return 0 if fail, 1 otherwise
  *
- * Pack the sparse meshes of each group and create triangles and edges before
- * getting out of library
+ * Pack the sparse meshes of each group
  *
  */
 int PMMG_packParMesh( PMMG_pParMesh parmesh )
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 535d2b5f..b75ef5a4 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -259,25 +259,30 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   PMMG_CALLOC( parmesh,vGlobNum_tab,4*(nitem_grp_face),MMG5_int,"vGlobNum_tab",return 0 );
   PMMG_CALLOC( parmesh,ne_tmp_tab,nitem_grp_face+1,MMG5_int,"ne_tmp_tab",return 0 );
 
-  /** STEP 4 - Identify required edges. Put hash.item[key].k = -1 */
+  /** STEP 4 - Identify required edges. Put hash.item[key].k = -1. This step
+   * assumes that the required tags are consistent through the partitions. */
+
   /* Loop over tetra */
   for (k=1; k<=mesh->ne; k++) {
     pt = &mesh->tetra[k];
     if ( !MG_EOK(pt) )  continue;
 
-    /* Check whether the tetra with reference ref should be split */
-    // For now, the option nosplit is not supported
-    // So we assume all the elements should be split
+    /* Check whether the tetra with reference ref should be split.  If an edge
+      or face is at the interface of 2 partitions and on one partition, the
+      domain is marked as noSplit and on another, the domain is split, we
+      have to ensure the consistency of split along the interface. */
     // if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) continue;
 
     /** Step 4.1 - Identification of edges belonging to a required tet */
     /* If the tetra is required MG_REQ  */
     if ( pt->tag & MG_REQ ) {
+
+      np = -1;
       /* Loop over the edges */
       for (ia=0; ia<6; ia++) {
         ip0 = pt->v[MMG5_iare[ia][0]];
         ip1 = pt->v[MMG5_iare[ia][1]];
-        np  = -1;
+
         /* Add an edge to the edge table with hash.item[key].k = -1 */
         if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) )  return -1;
       }
@@ -299,10 +304,10 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       /* (a) otherwise loop over the edges */
       for (j=0; j<3; j++) {
 
-        /* (b) If the edges is not required, then continue */
+        /* (b) If the edge is not required, then continue */
         if ( !(pxt->tag[ MMG5_iarf[ia][j] ] & MG_REQ) ) continue;
 
-        /* (b) otherwise get the extremity of the edges ... */
+        /* (b) otherwise get the extremity of the edge ... */
         ip0 = pt->v[MMG5_idir[ia][MMG5_inxt2[j]]];
         ip1 = pt->v[MMG5_idir[ia][MMG5_iprv2[j]]];
         np  = -1;
@@ -316,8 +321,10 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   /** STEP 5 - Create points at iso-value. Fill or update edge hash table */
   /** STEP 5.1 - Create new points located on parallel interfaces */
   /* Internal edge communicator - intvalues stores:
-    - point position as in node2int_node_comm_index2 if the edge is split
-    - otherwise, -1 */
+    - point position in the shared buffer of points (as in node2int_node_comm_index2)
+      if the edge is split
+    - otherwise, -1
+  */
   PMMG_CALLOC(parmesh,int_edge_comm->intvalues,nitem_int_edge,int,"int_edge_comm intvalues",return 0);
 
   /* Loop on the internal edge communicator */
@@ -375,7 +382,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     np = MMG3D_newPt(mesh,c,MG_PARBDY+MG_NOSURF+MG_REQ,src);
 
     /* Update internal communicators of node and edge. For int edge comm
-       intvalues stores the point position as in node2int_node_comm_index2 */
+       intvalues stores the point position in the shared buffer of point as in
+       node2int_node_comm_index2 */
     grp->node2int_node_comm_index1[nitem_grp_node]=np; // Add this new point in int node comm index1
     grp->node2int_node_comm_index2[nitem_grp_node]=nitem_int_node; // Position in int node comm
     int_edge_comm->intvalues[pos] = nitem_int_node; // In int edge comm, assign position of node in int comm
@@ -459,6 +467,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     }
 
     /* Update hash table */
+    // Remark: to call succesfully hashUpdate, the edge must already exist in
+    // the hash table
     MMG5_hashUpdate(&hash,ip0,ip1,np);
 
   }
@@ -516,7 +526,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       ip1 = pt->v[MMG5_iare[ia][1]];
       np  = MMG5_hashGet(&hash,ip0,ip1);
 
-      /* If np>0 (i.e. hash.item[key].k != [0;-1]), this edge has already been split, pass to the next edge */
+      /* If np>0 (i.e. hash.item[key].k != [0;-1]), this edge has already been
+       * split or is required (user-required or //), pass to the next edge */
       if ( np>0 ) continue;
 
       /* Check whether an entity with reference ref should be split */
@@ -705,9 +716,9 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     /* STEP 6.2.1 - Find global numbering and split pattern of the tetra */
     /* If the tetra has already a flag (flag !=0) - then it has already been processed */
     already_split = 0;
-    if (pt->flag) {
+    if ( pt->flag ) {
       /* If flag>0, we need to update the face communicators; */
-      if (pt->flag > 0) {
+      if ( pt->flag > 0 ) {
         already_split = 1;             // Identify the tetra as already split
         idx_tmp = pt->mark;            // Index of ne_tmp stored in ne_tmp_tab
         ne_tmp  = ne_tmp_tab[idx_tmp]; // Old total number of tetras just after the split of this old tetra
@@ -717,7 +728,9 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
         }
       }
       /* Otherwise, if flag equals to -1 (flag<0), the tetra does not need to be split. Pass to the next tetra. */
-      else if (pt->flag < 0) continue;
+      else if ( pt->flag < 0 ) {
+        continue;
+      }
     }
     /* Otherwise, if flag=0, the tetra has never been processed. */
     else {
@@ -725,8 +738,11 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       memset(vx,0,6*sizeof(MMG5_int));
       for (ia=0; ia<6; ia++) {
         vx[ia] = MMG5_hashGet(&hash,pt->v[MMG5_iare[ia][0]],pt->v[MMG5_iare[ia][1]]);
-        if ( vx[ia] > 0 )
+
+        if ( vx[ia] > 0 ) {
           MG_SET(pt->flag,ia);
+        }
+
       }
 
       /* Get and store global num of the tetra vertices */
@@ -1693,9 +1709,7 @@ int PMMG_ls(PMMG_pParMesh parmesh) {
     return 0;
   }
 
-// #ifndef NDEBUG
   assert ( PMMG_check_extEdgeComm ( parmesh,parmesh->info.read_comm ) );
-// #endif
 
   /** Discretization of the implicit function - Cut tetra */
   if ( !PMMG_cuttet_ls(parmesh) ) {

From a9ca05a3124c8fa90c151799dc8eb4e9a524e046 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 27 Sep 2024 15:58:49 +0200
Subject: [PATCH 02/11] Add overlap creation - deletion in PMMG_cuttet_ls if
 DEV_OVERLAP environment variable is exported. For now the overlap is not
 used.

---
 src/ls_pmmg.c | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index b75ef5a4..dde0f6c7 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -152,6 +152,17 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
   ne_init = mesh->ne; // Initial number of tetra - before ls - needed in step 6.3
 
+
+  /* Create overlap: the overlap creation uses the point flags that are used in
+   * the current function too so it cannot be called further. */
+  const char *env_dev = getenv("DEV_OVERLAP");
+
+  if(env_dev){
+    if ( !PMMG_create_overlap(parmesh,parmesh->info.read_comm) ) {
+      return PMMG_STRONGFAILURE;
+    }
+  }
+
   /** STEP 1 - Reset flags */
   for (k=1; k<=mesh->np; k++)
     mesh->point[k].flag = 0;
@@ -170,6 +181,10 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       continue;
     }
 
+    if ( pt->tag & MG_OVERLAP ) {
+      continue;
+    }
+
     /* Loop over edges */
     for (ia=0; ia<6; ia++) {
 
@@ -313,11 +328,21 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
         np  = -1;
 
         /* (c) ... and add an edge to the edge table with hash.item[key].k = -1 */
+        printf("2 - hash %d %d\n",ip0,ip1);
+
         if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) )  return -1;
       }
     }
   }
 
+  /* Delete overlap */
+  if ( env_dev ) {
+    if ( !PMMG_delete_overlap(parmesh,parmesh->info.read_comm) ) {
+      return PMMG_STRONGFAILURE;
+    }
+  }
+
+
   /** STEP 5 - Create points at iso-value. Fill or update edge hash table */
   /** STEP 5.1 - Create new points located on parallel interfaces */
   /* Internal edge communicator - intvalues stores:

From b567cc33d757116a41353ac1f64a21a4d670985d Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Fri, 27 Sep 2024 16:05:52 +0200
Subject: [PATCH 03/11] Make global numbering compatible with overlap - not
 used for now.

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

diff --git a/src/libparmmg.c b/src/libparmmg.c
index 76a05f88..7cc3ce32 100644
--- a/src/libparmmg.c
+++ b/src/libparmmg.c
@@ -1103,9 +1103,13 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
   }
 
   /* Store owner in the point flag */
+  MMG5_int np_overlap = 0;
   for( ip = 1; ip <= mesh->np; ip++ ) {
     ppt = &mesh->point[ip];
-    if (ppt->tag & MG_OVERLAP) continue;
+    if (ppt->tag & MG_OVERLAP) {
+      ++np_overlap;
+      continue;
+    }
     ppt->flag = parmesh->myrank;
   }
 
@@ -1117,7 +1121,7 @@ int PMMG_Compute_verticesGloNum( PMMG_pParMesh parmesh,MPI_Comm comm ){
   }
 
   /* Count owned nodes */
-  nowned = mesh->np;
+  nowned = mesh->np - np_overlap;
   for( idx = 0; idx < int_node_comm->nitem; idx++ ) {
     if( intvalues[idx] != parmesh->myrank ) nowned--;
   }

From cd2ff2303fa71a70cb6ca05d4e0c119847093c37 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 1 Oct 2024 16:50:20 +0200
Subject: [PATCH 04/11] First version of nosplit multimat mode.

---
 src/hash_pmmg.c |  17 ++---
 src/ls_pmmg.c   | 186 +++++++++++++++++++++++++++++++++---------------
 src/parmmg.h    |   2 +
 3 files changed, 138 insertions(+), 67 deletions(-)

diff --git a/src/hash_pmmg.c b/src/hash_pmmg.c
index dd32e87c..3dc6f8dd 100644
--- a/src/hash_pmmg.c
+++ b/src/hash_pmmg.c
@@ -329,17 +329,19 @@ int PMMG_chkBdryTria(MMG5_pMesh mesh, MMG5_int* permtria) {
  * \param hash pointer toward the hash table of edges.
  * \param a index of the first extremity of the edge.
  * \param b index of the second extremity of the edge.
- * \param k index of new point along the edge [a,b].
- * \param s If ls mode in ParMmg: index of new point in internal edge communicator;
- *          otherwise, the value stored in variable s.
+ *
+ * \param s If ls mode: 1 for a parallel edge that belongs
+ * to at least one element whose reference has to be splitted (either because we
+ * are not in multi-mat mode or because the reference is split in multi-mat
+ * mode). To avoid useless checks, some non parallel edges may be marked.
+ * If the edge belongs only to non-split references, s has to be 0.
+ *
  * \return PMMG_SUCCESS if success, PMMG_FAILURE if fail (edge is not found).
  *
- * Update the index of the new point stored along the edge \f$[a;b]\f$ (similar to MMG5_hashUpdate in mmg).
- * If ls mode in ParMmg: update the index of the new point in internal edge communicator in variable s;
- * otherwise, update the value stored in variable s.
+ * Update the value of the s field stored along the edge \f$[a;b]\f$
  *
  */
-int PMMG_hashUpdate_all(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int k,MMG5_int s) {
+int PMMG_hashUpdate_s(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int s) {
   MMG5_hedge  *ph;
   MMG5_int     key;
   MMG5_int    ia,ib;
@@ -351,7 +353,6 @@ int PMMG_hashUpdate_all(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int k,MMG5_i
 
   while ( ph->a ) {
     if ( ph->a == ia && ph->b == ib ) {
-      ph->k = k;
       ph->s = s;
       return PMMG_SUCCESS;
     }
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index dde0f6c7..73a9539c 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -37,15 +37,48 @@
 #include "parmmg.h"
 #include "mmgexterns_private.h"
 
+/**
+ * \param hash pointer to the hash table of edges.
+ * \param a index of the first extremity of the edge.
+ * \param b index of the second extremity of the edge.
+ * \return the index of point stored along \f$[a;b]\f$.
+ *
+ * Find the index of point stored along \f$[a;b]\f$.
+ * \note In ParMmg, hash_pmmg.c: PMMG_hashGet_all gets k and s at the same time;
+ * \note PMMG_hashGet_all might be moved here if needed one day in mmg.
+ *
+ */
+static inline
+int PMMG_hashMark_splitEdge(MMG5_Hash *hash,MMG5_int a,MMG5_int b) {
+  MMG5_hedge  *ph;
+  MMG5_int    key;
+  MMG5_int    ia,ib;
+
+  if ( !hash->item ) return 0;
+
+  ia  = MG_MIN(a,b);
+  ib  = MG_MAX(a,b);
+  key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
+  ph  = &hash->item[key];
+
+  if ( !ph->a )  return 0;
+  if ( ph->a == ia && ph->b == ib )  return ph->k;
+  while ( ph->nxt ) {
+    ph = &hash->item[ph->nxt];
+    if ( ph->a == ia && ph->b == ib )  return ph->k;
+  }
+  return 0;
+}
+
+
 /**
  * \param parmesh pointer toward a parmesh structure
  * \param mesh pointer toward the mesh structure.
  * \param sol pointer toward the level-set values.
  * \param met pointer toward a metric (non-mandatory).
  *
- * \return 1 if success, 0 otherwise.
- *
- * \todo Multimaterial isSplit() on // interface
+ * \return 1 if success, -2 if fail in overlap creation or deletion,
+ * 0 or -1 if fail elsewhere.
  *
  * Proceed to discretization of the implicit function carried by sol into mesh,
  * once values of sol have been snapped/checked
@@ -108,9 +141,6 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   int idx_edge_ext,idx_edge_int,idx_edge_mesh;
   int idx_face_ext,idx_face_int,val_face;
 
-  if ( parmesh->myrank == parmesh->info.root )
-    fprintf(stdout,"\n      ## PMMG_cuttet_ls: Multimaterial not fully supported yet.\n");
-
   /* Ensure only one group on each proc */
   assert ( (parmesh->ngrp == 1 || parmesh->ngrp == 0) &&
            "Implemented for 1 group per rank" );
@@ -119,20 +149,6 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   met  = parmesh->listgrp[0].met;
   sol  = parmesh->listgrp[0].ls;
 
-  /* For now, does not support `nosplit` option in multimat  */
-  // To be removed when supported
-  for (i=0; i<mesh->info.nmat; i++) {
-    mat    = &mesh->info.mat[i];
-    ref    = mat->ref;
-    refint = mat->rin;
-    refext = mat->rex;
-    if ( (ref == refint) && (ref == refext) ) {
-      if ( parmesh->myrank == parmesh->info.root )
-        fprintf(stderr,"\n      -- ERROR: The option `nosplit` in multimat is not supported yet.\n");
-      return 0;
-    }
-  }
-
   /* Initialization */
   grp = &parmesh->listgrp[0];
   field = grp->field;
@@ -153,14 +169,17 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   ne_init = mesh->ne; // Initial number of tetra - before ls - needed in step 6.3
 
 
-  /* Create overlap: the overlap creation uses the point flags that are used in
-   * the current function too so it cannot be called further. */
-  const char *env_dev = getenv("DEV_OVERLAP");
+  /* Create an overlap to check if edges along the partition interfaces have to
+   * be split in multi-material mode (it allows to check if an edge belonging to
+   * a "nosplit" material on a partition belongs to a split material on another
+   * partition and thus has to be split) to maintains the mesh consistency.
+   */
 
-  if(env_dev){
-    if ( !PMMG_create_overlap(parmesh,parmesh->info.read_comm) ) {
-      return PMMG_STRONGFAILURE;
-    }
+  // Remark: The overlap creation uses the point flags that are used in the
+  // current function too so it cannot be called further.
+  ier = 1;
+  if ( !PMMG_create_overlap(parmesh,parmesh->info.read_comm) ) {
+    ier = -2;
   }
 
   /** STEP 1 - Reset flags */
@@ -276,22 +295,27 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
   /** STEP 4 - Identify required edges. Put hash.item[key].k = -1. This step
    * assumes that the required tags are consistent through the partitions. */
-
   /* Loop over tetra */
   for (k=1; k<=mesh->ne; k++) {
     pt = &mesh->tetra[k];
     if ( !MG_EOK(pt) )  continue;
 
-    /* Check whether the tetra with reference ref should be split.  If an edge
-      or face is at the interface of 2 partitions and on one partition, the
-      domain is marked as noSplit and on another, the domain is split, we
-      have to ensure the consistency of split along the interface. */
-    // if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) continue;
+    /* Check whether the tetra with reference ref should be split and, if yes,
+     * use the "s" field of the hashed edge to mark parallel edges that have to
+     * be split.
+     *
+     * If an edge or face is at the interface of 2 partitions and on one
+     * partition, the domain is marked as noSplit and on another, the domain is
+     * split, the overlap allows to mark the parallel edges that have to be
+     * split and to ensure the consistency of split along the interface. */
+    int is_split = MMG5_isSplit(mesh,pt->ref,&refint,&refext);
 
     /** Step 4.1 - Identification of edges belonging to a required tet */
-    /* If the tetra is required MG_REQ  */
-    if ( pt->tag & MG_REQ ) {
-
+    /* Add required edges of required tetras to the hash table. The
+       overlap has to be ignored because we don't want to hash the edges of the
+       overlap and interface edges are marked as required so they will be added
+       to the hash table in the next step. */
+    if ( (pt->tag & MG_REQ) && !(pt->tag & MG_OVERLAP) ) {
       np = -1;
       /* Loop over the edges */
       for (ia=0; ia<6; ia++) {
@@ -299,12 +323,17 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
         ip1 = pt->v[MMG5_iare[ia][1]];
 
         /* Add an edge to the edge table with hash.item[key].k = -1 */
-        if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) )  return -1;
+        if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) ) {
+#warning check handling of ret val
+          return -1;
+        }
       }
       continue;
     }
 
-    /** Step 4.2 - Identification of edges belonging to a (par)boundary or being explicitely required */
+    /** Step 4.2 - Identification of edges belonging to a (par)boundary or being
+     * explicitely required. Here overlap tetra are automatically ignored
+     * because the xt field is not transferred in the overlap. */
     /* If the xtetra associated to this tetra exists */
     if ( !pt->xt ) continue;
 
@@ -328,21 +357,55 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
         np  = -1;
 
         /* (c) ... and add an edge to the edge table with hash.item[key].k = -1 */
-        printf("2 - hash %d %d\n",ip0,ip1);
-
-        if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) )  return -1;
+        if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) ) {
+#warning check handling of ret val
+          return -1;
+        }
       }
     }
   }
 
-  /* Delete overlap */
-  if ( env_dev ) {
-    if ( !PMMG_delete_overlap(parmesh,parmesh->info.read_comm) ) {
-      return PMMG_STRONGFAILURE;
+  /** Use the overlap to check which parallel edges have to be split and store
+   * this info in the s field of the hash table. */
+  for (k=1; k<=mesh->ne; k++) {
+    pt = &mesh->tetra[k];
+    if ( !MG_EOK(pt) )  continue;
+
+    /* Check whether the tetra with reference ref should be split and, if yes,
+     * use the "s" field of the hashed edge to mark parallel edges that have to
+     * be split.
+     *
+     * If an edge or face is at the interface of 2 partitions and on one
+     * partition, the domain is marked as noSplit and on another, the domain is
+     * split, the overlap allows to mark the parallel edges that have to be
+     * split and to ensure the consistency of split along the interface. */
+    int is_split = MMG5_isSplit(mesh,pt->ref,&refint,&refext);
+
+    if ( !is_split ) {
+      continue;
+    }
+    /* Loop over the edges and mark them as belonging to a "split" reference. */
+    for (ia=0; ia<6; ia++) {
+      ip0 = pt->v[MMG5_iare[ia][0]];
+      ip1 = pt->v[MMG5_iare[ia][1]];
+
+      /* Remark: only edges that already exist in the hash table are updated,
+       * thus only required edges are updated. As all parallel edges are
+       * required and already added to the hash table, we will be able to store
+       * the required info. In other cases (edges that doesn't exist in the hash
+       * table), the PMMG_hashUpdate_s function return PMMG_FAILURE, which is
+       * expected and harmless as we will not need to check if the edge is split
+       * or not. */
+      PMMG_hashUpdate_s(&hash,ip0,ip1,1);
     }
   }
 
 
+  /* Delete overlap */
+  if ( !PMMG_delete_overlap(parmesh,parmesh->info.read_comm) ) {
+    return PMMG_STRONGFAILURE;
+  }
+
   /** STEP 5 - Create points at iso-value. Fill or update edge hash table */
   /** STEP 5.1 - Create new points located on parallel interfaces */
   /* Internal edge communicator - intvalues stores:
@@ -350,7 +413,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       if the edge is split
     - otherwise, -1
   */
-  PMMG_CALLOC(parmesh,int_edge_comm->intvalues,nitem_int_edge,int,"int_edge_comm intvalues",return 0);
+  PMMG_CALLOC(parmesh,int_edge_comm->intvalues,nitem_int_edge,int,
+              "int_edge_comm intvalues",return 0);
 
   /* Loop on the internal edge communicator */
   for (i=0; i < nitem_grp_edge; i++) {
@@ -363,15 +427,17 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     ip0 = mesh->edge[ie].a;
     ip1 = mesh->edge[ie].b;
 
-    // TODO:: Multimaterial - Check whether an entity with reference ref should be split
-    // Pb1: This function does not work here because we come from an edge and not a tetra
-    // Need to find a way to know if we need to split or not this edge...
-    // Pb2: even if I find a way to know if this edge need to be split on this partition
-    // still no way to know if it is split from another partition if we are in
-    // the case where on partition 1/material 1 and partition 2/material 2 and
-    // we split mat1 but we do not split mat2. In that case, on partition 2, we will never know
-    // that edges on // interface need to be split.
-    // if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) continue;
+    MMG5_int dummy,is_split;
+    if ( PMMG_SUCCESS != PMMG_hashGet_all(&hash,ip0,ip1,&dummy,&is_split) ) {
+#warning check retval handling
+      return -1;
+    }
+
+    if ( !is_split ) {
+      /* The parallel edge belongs to a domain that is not split */
+      printf("Warning: no split edge: %d %d\n",ip0,ip1);
+      continue;
+    }
 
     /* STEP 5.1.1 - Create a new point if this edge needs to be split */
     /* Check the ls value at the edge nodes */
@@ -555,10 +621,12 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
        * split or is required (user-required or //), pass to the next edge */
       if ( np>0 ) continue;
 
+
       /* Check whether an entity with reference ref should be split */
-      // For now, the option nosplit is not supported
-      // So we assume all the elements should be split
-      // if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) continue;
+      if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) {
+        printf("rank %d: tetra %d is no split\n",parmesh->myrank,k);
+        continue;
+      }
 
       /* STEP 5.3.1 - Create a new point if this edge needs to be split */
       /* Check the ls value at the edge nodes */
@@ -1737,7 +1805,7 @@ int PMMG_ls(PMMG_pParMesh parmesh) {
   assert ( PMMG_check_extEdgeComm ( parmesh,parmesh->info.read_comm ) );
 
   /** Discretization of the implicit function - Cut tetra */
-  if ( !PMMG_cuttet_ls(parmesh) ) {
+  if ( 0 >= PMMG_cuttet_ls(parmesh) ) {
     if ( parmesh->myrank == parmesh->info.root )
       fprintf(stderr,"\n  ## Problem in discretizing implicit function. Exit program.\n");
     return 0;
diff --git a/src/parmmg.h b/src/parmmg.h
index 3388c85b..15ab277e 100644
--- a/src/parmmg.h
+++ b/src/parmmg.h
@@ -454,6 +454,8 @@ 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_hashOldPar_pmmg( PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_Hash *hash );
+int PMMG_hashUpdate_s(MMG5_Hash *hash,MMG5_int ip0,MMG5_int ip1, MMG5_int s);
+MMG5_int PMMG_hashGet_all(MMG5_Hash *hash,MMG5_int a,MMG5_int b,MMG5_int *k,MMG5_int *s);
 
 /* Overlap functions */
 int PMMG_create_overlap(PMMG_pParMesh parmesh,MPI_Comm comm);

From ca4e40504513c566a1ea131f8f2dc0e701e8ec27 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 10:37:05 +0100
Subject: [PATCH 05/11] Avoid deadlock in cas of failure in cuttet and multimat
 options.

---
 src/analys_pmmg.c |   9 +-
 src/ls_pmmg.c     | 296 +++++++++++++++++++++++++++-------------------
 2 files changed, 179 insertions(+), 126 deletions(-)

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 8a297bcd..2691550f 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -2788,7 +2788,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   /* Hash parallel edges from tetra and face communicator: store edges and the
      MG_PARBDY tag (other edge tags are not stored).*/
   if( ier && (PMMG_hashPar_fromFaceComm( parmesh,&hpar ) != PMMG_SUCCESS) ) {
-    fprintf(stderr,"\n  ## Warning: impossible to compute the hash parallel edge\n");
+    fprintf(stderr,"\n  ## Impossible to compute the hash parallel edge."
+            " Exit program.\n");
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,hpar.geom);
@@ -2802,7 +2803,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   }
 
   if ( !PMMG_build_edgeComm( parmesh,mesh,&hpar,comm ) ) {
-    fprintf(stderr,"\n  ## Warning: Impossible to build edge communicator\n");
+    fprintf(stderr,"\n  ## Impossible to build edge communicator."
+            " Exit program\n");
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
     MMG5_DEL_MEM(mesh,hpar.geom);
@@ -2811,7 +2813,8 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 
   /* Compute global node numbering and store it in ppt->tmp */
   if( !PMMG_Compute_verticesGloNum( parmesh,comm ) ) {
-    fprintf(stderr,"\n  ## Warning: impossible to compute node global numbering\n");
+    fprintf(stderr,"\n  ## Impossible to compute node global numbering."
+            " Exit program\n");
 
     MMG5_DEL_MEM(mesh,mesh->htab.geom);
     MMG5_DEL_MEM(mesh,hash.item);
diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index c6c2c0e9..841d694a 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -70,6 +70,16 @@ int PMMG_hashMark_splitEdge(MMG5_Hash *hash,MMG5_int a,MMG5_int b) {
   return 0;
 }
 
+#ifndef NDEBUG
+/* Checks over external communicators consistenc, involving comms, are not
+ * called */
+#define ERR_RULE return 0
+
+#else
+/* Checks over external communicators consistenc, involving comms, are called */
+#define ERR_RULE MPI_Abort(parmesh->comm,PMMG_TMPFAILURE)
+
+#endif
 
 /**
  * \param parmesh pointer toward a parmesh structure
@@ -77,12 +87,14 @@ int PMMG_hashMark_splitEdge(MMG5_Hash *hash,MMG5_int a,MMG5_int b) {
  * \param sol pointer toward the level-set values.
  * \param met pointer toward a metric (non-mandatory).
  *
- * \return 1 if success, -2 if fail in overlap creation or deletion,
- * 0 or -1 if fail elsewhere.
+ * \return 1 if success, 0 if fail.
  *
  * Proceed to discretization of the implicit function carried by sol into mesh,
  * once values of sol have been snapped/checked
  *
+ * \todo all MPI_Abort have to be removed and replaced by a clean error handling
+ * without deadlocks.
+ *
  */
 int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   MMG5_pTetra  pt,pt0;
@@ -142,8 +154,11 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   int idx_face_ext,idx_face_int,val_face;
 
   /* Ensure only one group on each proc */
-  assert ( (parmesh->ngrp == 1 || parmesh->ngrp == 0) &&
-           "Implemented for 1 group per rank" );
+  if ( !parmesh->ngrp ) {
+    return 1;
+  }
+
+  assert ( parmesh->ngrp == 1 && "Implemented for 1 group per rank" );
 
   mesh = parmesh->listgrp[0].mesh;
   met  = parmesh->listgrp[0].met;
@@ -179,7 +194,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   // current function too so it cannot be called further.
   ier = 1;
   if ( !PMMG_create_overlap(parmesh,parmesh->info.read_comm) ) {
-    ier = -2;
+    ier = 0;
   }
 
   /** STEP 1 - Reset flags */
@@ -233,7 +248,10 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       }
     }
   }
-  if ( ! nb )  return 1;
+  if ( ier && !nb ) {
+    /* Succeed but no overlap to create */
+    return ier;
+  }
 
   /* TODO:: test if the number of point proc by proc is correct */
   // Cannot be done here as it is an approximation. Otherwise, need to robustify step 2 above.
@@ -243,28 +261,34 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
   /** STEP 3 - Memory allocation */
   /* STEP 3.1 - Initialize hash table for edges */
-  if ( !MMG5_hashNew(mesh,&hash,nb,7*nb) ) return 0;
+  if ( !MMG5_hashNew(mesh,&hash,nb,7*nb) ) {
+    ier = 0;
+  }
 
   /* STEP 3.2 - Realloc internal node communicators */
   PMMG_REALLOC(parmesh, grp->node2int_node_comm_index1,
                nitem_grp_node+nitem_grp_edge,
                nitem_grp_node,
-               int,"Allocation of node2int_node_comm_index1", return 0);
+               int,"Allocation of node2int_node_comm_index1",
+               ier = 0);
   PMMG_REALLOC(parmesh, grp->node2int_node_comm_index2,
                nitem_grp_node+nitem_grp_edge,
                nitem_grp_node,
-               int,"Allocation of node2int_node_comm_index2", return 0);
+               int,"Allocation of node2int_node_comm_index2",
+               ier = 0);
   nitem_grp_node_firstalloc = nitem_grp_node+nitem_grp_edge;
 
   /* STEP 3.3 - Realloc internal face communicators */
   PMMG_REALLOC(parmesh, grp->face2int_face_comm_index1,
                3*nitem_grp_face,
                nitem_grp_face,
-               int,"Allocation of face2int_face_comm_index1", return 0);
+               int,"Allocation of face2int_face_comm_index1",
+               ier = 0);
   PMMG_REALLOC(parmesh, grp->face2int_face_comm_index2,
                3*nitem_grp_face,
                nitem_grp_face,
-               int,"Allocation of face2int_face_comm_index2", return 0);
+               int,"Allocation of face2int_face_comm_index2",
+               ier = 0);
   nitem_grp_face_firstalloc = 3*nitem_grp_face;
 
   /* STEP 3.4 - Realloc external node communicator */
@@ -274,7 +298,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     PMMG_REALLOC(parmesh,ext_node_comm->int_comm_index,
                  nitem_ext_node+2*nb,
                  nitem_ext_node,
-                 int,"Allocation of external node communicator",return 0);
+                 int,"Allocation of external node communicator",
+                 ier = 0);
     ext_node_comm->nitem_to_share = nitem_ext_node+2*nb;
   }
 
@@ -285,13 +310,20 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     PMMG_REALLOC(parmesh,ext_face_comm->int_comm_index,
                  nitem_ext_face+2*nitem_ext_face,
                  nitem_ext_face,
-                 int,"Allocation of external face communicator",return 0);
+                 int,"Allocation of external face communicator",
+                 ier = 0);
     ext_face_comm->nitem_to_share = nitem_ext_face+2*nitem_ext_face;
   }
 
   /* STEP 3.6 - Allocate all the other variables needed to be allocated ! */
-  PMMG_CALLOC( parmesh,vGlobNum_tab,4*(nitem_grp_face),MMG5_int,"vGlobNum_tab",return 0 );
-  PMMG_CALLOC( parmesh,ne_tmp_tab,nitem_grp_face+1,MMG5_int,"ne_tmp_tab",return 0 );
+  PMMG_CALLOC( parmesh,vGlobNum_tab,4*(nitem_grp_face),MMG5_int,"vGlobNum_tab",ier = 0 );
+  PMMG_CALLOC( parmesh,ne_tmp_tab,nitem_grp_face+1,MMG5_int,"ne_tmp_tab",ier = 0 );
+
+
+  if ( !ier ) {
+    /* One alloc has failed: Avoid segfault or deadlock */
+    MPI_Abort(parmesh->comm,PMMG_TMPFAILURE);
+  }
 
   /** STEP 4 - Identify required edges. Put hash.item[key].k = -1. This step
    * assumes that the required tags are consistent through the partitions. */
@@ -324,8 +356,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
         /* Add an edge to the edge table with hash.item[key].k = -1 */
         if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) ) {
-#warning check handling of ret val
-          return -1;
+          ier = -1;
         }
       }
       continue;
@@ -358,8 +389,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
         /* (c) ... and add an edge to the edge table with hash.item[key].k = -1 */
         if ( !MMG5_hashEdge(mesh,&hash,ip0,ip1,np) ) {
-#warning check handling of ret val
-          return -1;
+          ier = -1;
         }
       }
     }
@@ -403,7 +433,13 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
   /* Delete overlap */
   if ( !PMMG_delete_overlap(parmesh,parmesh->info.read_comm) ) {
-    return PMMG_STRONGFAILURE;
+    ier = MG_MIN(ier,0);
+  }
+
+
+  if ( ier < 1 ) {
+    /* Avoid deadlock in comms in build_edgeComm */
+    MPI_Abort(parmesh->comm,PMMG_TMPFAILURE);
   }
 
   /** STEP 5 - Create points at iso-value. Fill or update edge hash table */
@@ -414,7 +450,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     - otherwise, -1
   */
   PMMG_CALLOC(parmesh,int_edge_comm->intvalues,nitem_int_edge,int,
-              "int_edge_comm intvalues",return 0);
+              "int_edge_comm intvalues",
+              MPI_Abort(parmesh->comm,PMMG_TMPFAILURE));
 
   /* Loop on the internal edge communicator */
   for (i=0; i < nitem_grp_edge; i++) {
@@ -429,13 +466,11 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
     MMG5_int dummy,is_split;
     if ( PMMG_SUCCESS != PMMG_hashGet_all(&hash,ip0,ip1,&dummy,&is_split) ) {
-#warning check retval handling
-      return -1;
+      ier = 0;
     }
 
     if ( !is_split ) {
       /* The parallel edge belongs to a domain that is not split */
-      printf("Warning: no split edge: %d %d\n",ip0,ip1);
       continue;
     }
 
@@ -491,9 +526,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
                             fprintf(stderr,"\n  ## Error: %s: unable to"
                                     " allocate a new point\n",__func__);
                             MMG5_INCREASE_MEM_MESSAGE();
-                            return 0
-                            ,c,0,src);
-      if( met ) {
+                            ier = 0,c,0,src);
+      if( met && ier ) {
         if( met->m ) {
           MMG5_ADD_MEM(mesh,(met->size*(mesh->npmax-met->npmax))*sizeof(double),
                         "larger solution",
@@ -502,66 +536,75 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
                         mesh->npmax = oldnpmax;
                         mesh->np = mesh->npmax-1;
                         mesh->npnil = 0;
-                        return 0);
-          MMG5_SAFE_REALLOC(met->m,met->size*(met->npmax+1),
-                            met->size*(mesh->npmax+1),
-                            double,"larger solution",
-                            MMG5_SAFE_RECALLOC(mesh->point,mesh->npmax+1,oldnpmax+1,MMG5_Point,,);
-                            mesh->memCur -= (mesh->npmax - oldnpmax)*sizeof(MMG5_Point);
-                            mesh->npmax = oldnpmax;
-                            mesh->np = mesh->npmax-1;
-                            mesh->npnil = 0;
-                            return 0);
+                        ier = 0 );
+
+          if ( ier ) {
+            MMG5_SAFE_REALLOC(met->m,met->size*(met->npmax+1),
+                              met->size*(mesh->npmax+1),
+                              double,"larger solution",
+                              MMG5_SAFE_RECALLOC(mesh->point,mesh->npmax+1,oldnpmax+1,MMG5_Point,,);
+                              mesh->memCur -= (mesh->npmax - oldnpmax)*sizeof(MMG5_Point);
+                              mesh->npmax = oldnpmax;
+                              mesh->np = mesh->npmax-1;
+                              mesh->npnil = 0;
+                              ier = 0);
+          }
+          met->npmax = mesh->npmax;
         }
-        met->npmax = mesh->npmax;
       }
     }
 
-    /* For this new point, add the value of the solution, i.e. the isovalue 0 */
-    sol->m[np] = 0;
-
-    /* If user provide a metric, interpolate it at the new point */
-    if ( met && met->m ) {
-      if ( met->size > 1 ) {
-        ier = MMG3D_intmet33_ani_edge(met,ip0,ip1,np,s);
-      }
-      else {
-        ier = MMG5_intmet_iso_edge(met,ip0,ip1,np,s);
-      }
-      if ( ier <= 0 ) {
-        /* Unable to compute the metric */
-        fprintf(stderr,"\n  ## Error: %s: unable to"
-                " interpolate the metric during the level-set"
-                " discretization\n",__func__);
-        return 0;
-      }
-    }
+    if ( ier ) {
+      /* For this new point, add the value of the solution, i.e. the isovalue 0 */
+      sol->m[np] = 0;
 
-    /* If user provide fields, interpolate them at the new point */
-    if ( mesh->nsols ) {
-      for ( j=0; j<mesh->nsols; ++j ) {
-        psl    = field + j;
-        if ( field->size > 1 ) {
-          ier = MMG3D_intmet33_ani_edge(psl,ip0,ip1,np,s);
+      /* If user provide a metric, interpolate it at the new point */
+      if ( met && met->m ) {
+        if ( met->size > 1 ) {
+          ier = MMG3D_intmet33_ani_edge(met,ip0,ip1,np,s);
         }
         else {
-          ier = MMG5_intmet_iso_edge(psl,ip0,ip1,np,s);
+          ier = MMG5_intmet_iso_edge(met,ip0,ip1,np,s);
         }
         if ( ier <= 0 ) {
-          /* Unable to compute fields */
+          /* Unable to compute the metric */
           fprintf(stderr,"\n  ## Error: %s: unable to"
-                  " interpolate fields during the level-set"
+                  " interpolate the metric during the level-set"
                   " discretization\n",__func__);
-          return 0;
+          ier = 0;
+        }
+      }
+
+      /* If user provide fields, interpolate them at the new point */
+      if ( mesh->nsols ) {
+        for ( j=0; j<mesh->nsols; ++j ) {
+          psl    = field + j;
+          if ( field->size > 1 ) {
+            ier = MMG3D_intmet33_ani_edge(psl,ip0,ip1,np,s);
+          }
+          else {
+            ier = MMG5_intmet_iso_edge(psl,ip0,ip1,np,s);
+          }
+          if ( ier <= 0 ) {
+            /* Unable to compute fields */
+            fprintf(stderr,"\n  ## Error: %s: unable to"
+                    " interpolate fields during the level-set"
+                    " discretization\n",__func__);
+            ier = 0;
+          }
         }
       }
-    }
 
-    /* Update hash table */
-    // Remark: to call succesfully hashUpdate, the edge must already exist in
-    // the hash table
-    MMG5_hashUpdate(&hash,ip0,ip1,np);
+      /* Update hash table */
+      // Remark: to call succesfully hashUpdate, the edge must already exist in
+      // the hash table
+      MMG5_hashUpdate(&hash,ip0,ip1,np);
+    }
 
+    if ( !ier ) {
+      /* Avoid too long list of errors in case of failure */
+      MPI_Abort(parmesh->comm,PMMG_TMPFAILURE);
+    }
   }
 
   /** STEP 5.2 -  Update external node communicator  */
@@ -624,7 +667,6 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
       /* Check whether an entity with reference ref should be split */
       if ( !MMG5_isSplit(mesh,pt->ref,&refint,&refext) ) {
-        printf("rank %d: tetra %d is no split\n",parmesh->myrank,k);
         continue;
       }
 
@@ -672,9 +714,9 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
                              fprintf(stderr,"\n  ## Error: %s: unable to"
                                      " allocate a new point\n",__func__);
                              MMG5_INCREASE_MEM_MESSAGE();
-                             return 0
+                             ier = 0
                              ,c,0,src);
-        if( met ) {
+        if( ier && met ) {
           if( met->m ) {
             MMG5_ADD_MEM(mesh,(met->size*(mesh->npmax-met->npmax))*sizeof(double),
                          "larger solution",
@@ -683,59 +725,69 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
                          mesh->npmax = oldnpmax;
                          mesh->np = mesh->npmax-1;
                          mesh->npnil = 0;
-                         return 0);
-            MMG5_SAFE_REALLOC(met->m,met->size*(met->npmax+1),
-                              met->size*(mesh->npmax+1),
-                              double,"larger solution",
-                              MMG5_SAFE_RECALLOC(mesh->point,mesh->npmax+1,oldnpmax+1,MMG5_Point,,);
-                              mesh->memCur -= (mesh->npmax - oldnpmax)*sizeof(MMG5_Point);
-                              mesh->npmax = oldnpmax;
-                              mesh->np = mesh->npmax-1;
-                              mesh->npnil = 0;
-                              return 0);
+                         ier = 0);
+
+            if ( ier ) {
+              MMG5_SAFE_REALLOC(met->m,met->size*(met->npmax+1),
+                                met->size*(mesh->npmax+1),
+                                double,"larger solution",
+                                MMG5_SAFE_RECALLOC(mesh->point,mesh->npmax+1,oldnpmax+1,MMG5_Point,,);
+                                mesh->memCur -= (mesh->npmax - oldnpmax)*sizeof(MMG5_Point);
+                                mesh->npmax = oldnpmax;
+                                mesh->np = mesh->npmax-1;
+                                mesh->npnil = 0;
+                                ier = 0);
+              met->npmax = mesh->npmax;
+            }
           }
-          met->npmax = mesh->npmax;
         }
       }
 
-      /* For this new point, add the value of the solution, i.e. the isovalue 0 */
-      sol->m[np] = 0;
-
-      /* If user provide a metric, interpolate it at the new point */
-      if ( met && met->m ) {
-        if ( met->size > 1 ) {
-          ier = MMG3D_intmet33_ani(mesh,met,k,ia,np,s);
-        }
-        else {
-          ier = MMG5_intmet_iso(mesh,met,k,ia,np,s);
-        }
-        if ( ier <= 0 ) {
-          /* Unable to compute the metric */
-          fprintf(stderr,"\n  ## Error: %s: unable to"
-                  " interpolate the metric during the level-set"
-                  " discretization\n",__func__);
-          return 0;
-        }
-      }
+      if ( ier ) {
+        /* For this new point, add the value of the solution, i.e. the isovalue 0 */
+        sol->m[np] = 0;
 
-      /* If user provide fields, interpolate them at the new point */
-      if ( mesh->nsols ) {
-        for ( j=0; j<mesh->nsols; ++j ) {
-          psl    = field + j;
-          if ( field->size > 1 ) {
-            ier = MMG3D_intmet33_ani(mesh,psl,k,ia,np,s);
+        /* If user provide a metric, interpolate it at the new point */
+        if ( met && met->m ) {
+          if ( met->size > 1 ) {
+            ier = MMG3D_intmet33_ani(mesh,met,k,ia,np,s);
           }
           else {
-            ier = MMG5_intmet_iso(mesh,psl,k,ia,np,s);
+            ier = MMG5_intmet_iso(mesh,met,k,ia,np,s);
           }
           if ( ier <= 0 ) {
-            /* Unable to compute fields */
+            /* Unable to compute the metric */
             fprintf(stderr,"\n  ## Error: %s: unable to"
-                    " interpolate fields during the level-set"
+                    " interpolate the metric during the level-set"
                     " discretization\n",__func__);
-            return 0;
+            ier = 0;
           }
         }
+
+        /* If user provide fields, interpolate them at the new point */
+        if ( mesh->nsols ) {
+          for ( j=0; j<mesh->nsols; ++j ) {
+            psl    = field + j;
+            if ( field->size > 1 ) {
+              ier = MMG3D_intmet33_ani(mesh,psl,k,ia,np,s);
+            }
+            else {
+              ier = MMG5_intmet_iso(mesh,psl,k,ia,np,s);
+            }
+            if ( ier <= 0 ) {
+              /* Unable to compute fields */
+              fprintf(stderr,"\n  ## Error: %s: unable to"
+                      " interpolate fields during the level-set"
+                      " discretization\n",__func__);
+              ier = 0;
+            }
+          }
+        }
+      }
+
+      if ( !ier ) {
+        /* Avoid too long list of errors in case of failure */
+        MPI_Abort(parmesh->comm,PMMG_TMPFAILURE);
       }
 
       /* STEP 5.2.3 - Update edge hash table */
@@ -759,19 +811,16 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     }
   }
 
-#ifndef NDEBUG
   /** Check the internal node communicator */
   assert( PMMG_check_intNodeComm( parmesh ) );
 
   /** Check the external node communicator */
   assert( PMMG_check_extNodeComm( parmesh,parmesh->info.read_comm ) );
-#endif
 
   /** STEP 6 - Split according to tets flags */
   /** STEP 6.1 - Compute global node vertices */
   if ( !PMMG_Compute_verticesGloNum( parmesh,parmesh->comm ) ) {
-    fprintf(stderr,"\n\n\n  -- WARNING: IMPOSSIBLE TO COMPUTE NODE GLOBAL NUMBERING\n\n\n");
-    return 0;
+    fprintf(stderr,"\n  ## Error: impossible to compute node global numbering\n");
   }
 
   /** STEP 6.2 - Do the splitting for tetra on parallel interface */
@@ -781,7 +830,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
   nitem_grp_face_tmp = nitem_grp_face;
 
   /* Allocate internal face comm */
-  PMMG_CALLOC(parmesh,int_face_comm->intvalues,3*nitem_int_face,int,"int_face_comm intvalues",return 0);
+  PMMG_CALLOC(parmesh,int_face_comm->intvalues,3*nitem_int_face,int,
+              "int_face_comm intvalues",ERR_RULE);
 
   /* Loop over the internal faces communicator */
   for (i=0; i < nitem_grp_face; i++) {
@@ -924,7 +974,9 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       break;
     }
 
-    if ( !ier ) return 0;
+    if ( !ier ) {
+      ERR_RULE;
+    }
 
     if (pt->flag != -1) {
       /* STEP 6.2.3 - Update tag of edges in xtetra with MG_PARBDY */
@@ -1006,13 +1058,11 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     ext_face_comm->nitem = nitem_ext_face; // Update nbr of face in ext face comm
   }
 
-#ifndef NDEBUG
   /** Check the internal face communicator */
   assert( PMMG_check_intFaceComm( parmesh ) );
 
   /** Check the external face communicator */
   assert( PMMG_check_extFaceComm( parmesh,parmesh->info.read_comm ) );
-#endif
 
   /** STEP 6.4 - Do the splitting for tetra located elsewhere */
   /* Loop over tetra */

From 995c7cebad298f3583c39a7868a5a5d702c8123e Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 13:27:04 +0100
Subject: [PATCH 06/11] Add nosplit test case + clean debug things.

---
 cmake/testing/pmmg_tests.cmake | 17 ++++++++++++++++-
 src/analys_pmmg.c              |  2 +-
 src/libparmmg1.c               | 15 +++++++--------
 3 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/cmake/testing/pmmg_tests.cmake b/cmake/testing/pmmg_tests.cmake
index 63268dca..13433163 100644
--- a/cmake/testing/pmmg_tests.cmake
+++ b/cmake/testing/pmmg_tests.cmake
@@ -21,7 +21,7 @@ IF( BUILD_TESTING )
       ENDIF()
       EXECUTE_PROCESS(
         COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} fetch
-        COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} checkout cd45788c32f6
+        COMMAND ${GIT_EXECUTABLE} -C ${CI_DIR} checkout 5091f86924742
         TIMEOUT 20
         WORKING_DIRECTORY ${CI_DIR}
         #COMMAND_ECHO STDOUT
@@ -783,6 +783,21 @@ IF( BUILD_TESTING )
 
 
       endforeach()
+
+    endforeach()
+
+    foreach( MODE faces )
+
+      # Toy geom nosplit:  2 procs, ls_val=0.0 + remesh hsiz 0.1 + 5 iter + multimat nosplit
+      SET( NP 2 )
+      add_test( NAME ls-DisIn-toygeom-nosplit-${MODE}-${NP}
+        COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
+        ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-distributed-${MODE}-mat-edges.mesh -v 5
+        -f ${CI_DIR}/LevelSet/${NP}p_toygeom/nosplit.mmg3d
+        -hsiz 0.1 -niter 5
+        -ls 0.0
+        -sol ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-ls.sol
+        -out ${CI_DIR_RESULTS}/ls-DisIn-toygeom-nosplit-${MODE}-${NP}.o.mesh)
     endforeach()
 
     #***********************
diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
index 2691550f..ac799608 100644
--- a/src/analys_pmmg.c
+++ b/src/analys_pmmg.c
@@ -3102,7 +3102,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
 
 #ifdef USE_POINTMAP
   /* Initialize source point with input index */
-  int ip;
+  MMG5_int ip;
   for( ip = 1; ip <= mesh->np; ip++ )
     mesh->point[ip].src = ip;
 #endif
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 378a9671..baf91bfe 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -632,15 +632,14 @@ 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;
 
+#ifndef NDEBUG
     if ( !MMG5_chkmsh(mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
-
+#endif
 
     memset(&mesh->xtetra[mesh->xt+1],0,(mesh->xtmax-mesh->xt)*sizeof(MMG5_xTetra));
     memset(&mesh->xpoint[mesh->xp+1],0,(mesh->xpmax-mesh->xp)*sizeof(MMG5_xPoint));
@@ -683,18 +682,16 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
       met          = parmesh->listgrp[i].met;
       field        = parmesh->listgrp[i].field;
 
-
+#ifndef NDEBUG
       if ( !MMG5_chkmsh(mesh,1,1) ) {
         fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
         return 0;
       }
-
+#endif
 
 #warning Luca: until analysis is not ready
 #ifdef USE_POINTMAP
       for( k = 1; k <= mesh->np; k++ ) {
-#warning Algiane todo
-        //assert ( mesh->point[k].src == k);
         mesh->point[k].src = k;
       }
 #endif
@@ -767,6 +764,7 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           }
         }
 
+
 #ifdef PATTERN
         ier = MMG5_mmg3d1_pattern( mesh, met, permNodGlob );
 #else
@@ -967,13 +965,14 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
   ier = PMMG_merge_grps(parmesh,0);
   MPI_Allreduce( &ier, &ieresult, 1, MPI_INT, MPI_MIN, parmesh->comm );
 
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
-
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_STEPS ) {
     chrono(OFF,&(ctim[tim]));

From 3af369bea072519265e2c64b64c731d9ff8dbc70 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 13:30:05 +0100
Subject: [PATCH 07/11] Debug cleaning.

---
 src/grpsplit_pmmg.c      |  8 +++++++-
 src/libparmmg1.c         |  7 ++++---
 src/loadbalancing_pmmg.c | 11 ++++++-----
 3 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/src/grpsplit_pmmg.c b/src/grpsplit_pmmg.c
index 9362402c..fad23421 100644
--- a/src/grpsplit_pmmg.c
+++ b/src/grpsplit_pmmg.c
@@ -1849,13 +1849,15 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   if ( !ier ) {
     fprintf(stderr,"\n  ## Merge groups problem.\n");
   }
+
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
-
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));
@@ -1920,22 +1922,26 @@ int PMMG_split_n2mGrps(PMMG_pParMesh parmesh,int target,int fitMesh,int repartit
   }
 
   /** Split the group into the suitable number of groups */
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( ier )
     ier = PMMG_splitPart_grps(parmesh,target,fitMesh,repartitioning_mode);
 
+#ifndef NDEBUG
   for (int k=0; k<parmesh->ngrp; ++k ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index baf91bfe..7439f97a 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -268,17 +268,19 @@ int PMMG_packParMesh( PMMG_pParMesh parmesh )
     }
 
     /* to could save the mesh, the adjacency have to be correct */
-//    if ( mesh->info.ddebug ) {
+    if ( mesh->info.ddebug ) {
       if ( (!mesh->adja) && !MMG3D_hashTetra(mesh,1) ) {
         fprintf(stderr,"\n  ## Error: %s: tetra hashing problem. Exit program.\n",
                 __func__);
         return 0;
       }
+#ifndef NDEBUG
       if ( !MMG5_chkmsh(mesh,1,1) ) {
         fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
         return 0;
       }
-//    }
+#endif
+    }
   }
 
   return 1;
@@ -534,7 +536,6 @@ int PMMG_scotchCall( PMMG_pParMesh parmesh,int igrp,int *permNodGlob ) {
   return 1;
 }
 
-int PMMG_grp_to_saveMesh( PMMG_pParMesh parmesh, int i, char*  );
 /**
  * \param parmesh pointer toward a parmesh structure where the boundary entities
  * are stored into xtetra and xpoint strucutres
diff --git a/src/loadbalancing_pmmg.c b/src/loadbalancing_pmmg.c
index 8ba8b3c3..c2c882c4 100644
--- a/src/loadbalancing_pmmg.c
+++ b/src/loadbalancing_pmmg.c
@@ -44,8 +44,6 @@
  * Load balancing of the mesh groups over the processors.
  *
  */
-int PMMG_grp_to_saveMesh( PMMG_pParMesh parmesh, int i, char*  );
-
 int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   MMG5_pMesh mesh;
   int        ier,ier_glob,igrp,ne;
@@ -166,7 +164,7 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
   }
 
 #ifndef NDEBUG
-  for ( int i=0; i<parmesh->ngrp; ++i ) {
+  for ( i=0; i<parmesh->ngrp; ++i ) {
     if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
@@ -194,12 +192,15 @@ int PMMG_loadBalancing(PMMG_pParMesh parmesh,int partitioning_mode) {
       }
     }
   }
-  for (int k=0; k<parmesh->ngrp; ++k ) {
-    if ( !MMG5_chkmsh(parmesh->listgrp[k].mesh,1,1) ) {
+
+#ifndef NDEBUG
+  for (i=0; i<parmesh->ngrp; ++i ) {
+    if ( !MMG5_chkmsh(parmesh->listgrp[i].mesh,1,1) ) {
       fprintf(stderr,"  ##  Problem. Invalid mesh.\n");
       return 0;
     }
   }
+#endif
 
   if ( parmesh->info.imprim > PMMG_VERB_DETQUAL ) {
     chrono(OFF,&(ctim[tim]));

From 8de421649120976440a453ac8f21385b129668dc Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 14:40:52 +0100
Subject: [PATCH 08/11] Fix use of wrong comminicator in cuttet.

---
 src/ls_pmmg.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 841d694a..17031ed3 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -819,7 +819,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
 
   /** STEP 6 - Split according to tets flags */
   /** STEP 6.1 - Compute global node vertices */
-  if ( !PMMG_Compute_verticesGloNum( parmesh,parmesh->comm ) ) {
+  if ( !PMMG_Compute_verticesGloNum( parmesh,parmesh->info.read_comm ) ) {
     fprintf(stderr,"\n  ## Error: impossible to compute node global numbering\n");
   }
 

From a12d160530005e224e0bf167226bd8d5c377e0fa Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 14:41:27 +0100
Subject: [PATCH 09/11] Fix possibly invalid pointers after tetra realloc in
 cutet.

---
 src/ls_pmmg.c | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/ls_pmmg.c b/src/ls_pmmg.c
index 17031ed3..515a8630 100644
--- a/src/ls_pmmg.c
+++ b/src/ls_pmmg.c
@@ -909,8 +909,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     case 1: case 2: case 4: case 8: case 16: case 32: // 1 edge split
       if (!already_split) {
         ier = MMG5_split1(mesh,met,k,vx,1);
-        pt->flag = flag;           // Re-flag tetra k as the flag has been reset in the split
-        pt->mark = ns;             // Split number to recover info later if needed
+        mesh->tetra[k].flag = flag;           // Re-flag tetra k as the flag has been reset in the split
+        mesh->tetra[k].mark = ns;             // Split number to recover info later if needed
         ne_tmp_tab[ns] = mesh->ne; // Total number of tetras after this split
         ne_tmp         = mesh->ne; // Total number of tetras after this split
         ns++;                      // Incremente the total number of split
@@ -925,8 +925,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     case 20: case 5:  case 17: case 9: case 3:  case 10:
       if (!already_split) {
         ier = MMG5_split2sf_globNum(mesh,met,k,vx,vGlobNum,1);
-        pt->flag = flag;           // Re-flag tetra k as the flag has been reset in the split
-        pt->mark = ns;             // Split number to recover info later if needed
+        mesh->tetra[k].flag = flag;           // Re-flag tetra k as the flag has been reset in the split
+        mesh->tetra[k].mark = ns;             // Split number to recover info later if needed
         ne_tmp_tab[ns] = mesh->ne; // Total number of tetras after this split
         ne_tmp         = mesh->ne; // Total number of tetras after this split
         ns++;                      // Incremente the total number of split
@@ -940,8 +940,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     case 7: case 25: case 42: case 52: // 3 edges on conic configuration split
       if (!already_split) {
         ier = MMG5_split3cone_globNum(mesh,met,k,vx,vGlobNum,1);
-        pt->flag = flag;           // Re-flag tetra k as the flag has been reset in the split
-        pt->mark = ns;             // Split number to recover info later if needed
+        mesh->tetra[k].flag = flag;           // Re-flag tetra k as the flag has been reset in the split
+        mesh->tetra[k].mark = ns;             // Split number to recover info later if needed
         ne_tmp_tab[ns] = mesh->ne; // Total number of tetras after this split
         ne_tmp         = mesh->ne; // Total number of tetras after this split
         ns++;                      // Incremente the total number of split
@@ -955,8 +955,8 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     case 30: case 45: case 51: // 4 edges on opposite configuration split
       if (!already_split) {
         ier = MMG5_split4op_globNum(mesh,met,k,vx,vGlobNum,1);
-        pt->flag = flag;           // Re-flag tetra k as the flag has been reset in the split
-        pt->mark = ns;             // Split number to recover info later if needed
+        mesh->tetra[k].flag = flag;           // Re-flag tetra k as the flag has been reset in the split
+        mesh->tetra[k].mark = ns;             // Split number to recover info later if needed
         ne_tmp_tab[ns] = mesh->ne; // Total number of tetras after this split
         ne_tmp         = mesh->ne; // Total number of tetras after this split
         ns++;                      // Incremente the total number of split
@@ -978,7 +978,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       ERR_RULE;
     }
 
-    if (pt->flag != -1) {
+    if (mesh->tetra[k].flag != -1) {
       /* STEP 6.2.3 - Update tag of edges in xtetra with MG_PARBDY */
       for (j=0; j<3; j++) {
         k0 = tetra_sorted[j];
@@ -1012,7 +1012,7 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
       nitem_grp_face_tmp += nface_added;
       int_face_comm->nitem = nitem_grp_face_tmp;
     }
-    else if (pt->flag == -1) {
+    else if (mesh->tetra[k].flag == -1) {
       /* STEP 6.2.5 - Update internal face communicators for tetra not split */
       /* As the tetra is not split, the tetra index has not changed.
           Moreover, the node of face ifac is chosen to be the node with highest
@@ -1097,26 +1097,26 @@ int PMMG_cuttet_ls(PMMG_pParMesh parmesh){
     switch (flag) {
     case 1: case 2: case 4: case 8: case 16: case 32: // 1 edge split
       ier = MMG5_split1(mesh,met,k,vx,1);
-      pt->flag = flag; // Re-flag tetra k as the flag has been reset in the split
+      mesh->tetra[k].flag = flag; // Re-flag tetra k as the flag has been reset in the split
       ns++;            // Incremente the total number of split
       break;
 
     case 48: case 24: case 40: case 6: case 34: case 36: // 2 edges (same face) split
     case 20: case 5:  case 17: case 9: case 3:  case 10:
       ier = MMG5_split2sf_globNum(mesh,met,k,vx,vGlobNum,1);
-      pt->flag = flag; // Re-flag tetra k as the flag has been reset in the split
+      mesh->tetra[k].flag = flag; // Re-flag tetra k as the flag has been reset in the split
       ns++;            // Incremente the total number of split
       break;
 
     case 7: case 25: case 42: case 52: // 3 edges on conic configuration split
       ier = MMG5_split3cone_globNum(mesh,met,k,vx,vGlobNum,1);
-      pt->flag = flag; // Re-flag tetra k as the flag has been reset in the split
+      mesh->tetra[k].flag = flag; // Re-flag tetra k as the flag has been reset in the split
       ns++;            // Incremente the total number of split
       break;
 
     case 30: case 45: case 51: // 4 edges on opposite configuration split
       ier = MMG5_split4op_globNum(mesh,met,k,vx,vGlobNum,1);
-      pt->flag = flag; // Re-flag tetra k as the flag has been reset in the split
+      mesh->tetra[k].flag = flag; // Re-flag tetra k as the flag has been reset in the split
       ns++;            // Incremente the total number of split
       break;
 

From 7e8a2ddfc185d6eec979a8a049acdffb366e2a3c Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Tue, 29 Oct 2024 21:26:46 +0100
Subject: [PATCH 10/11] Fix memory issue in permutation array with scotch
 renum.

---
 cmake/testing/pmmg_tests.cmake |  6 ++++--
 src/libparmmg1.c               | 20 +++++++++++++++-----
 2 files changed, 19 insertions(+), 7 deletions(-)

diff --git a/cmake/testing/pmmg_tests.cmake b/cmake/testing/pmmg_tests.cmake
index 13433163..d44979f8 100644
--- a/cmake/testing/pmmg_tests.cmake
+++ b/cmake/testing/pmmg_tests.cmake
@@ -788,13 +788,15 @@ IF( BUILD_TESTING )
 
     foreach( MODE faces )
 
-      # Toy geom nosplit:  2 procs, ls_val=0.0 + remesh hsiz 0.1 + 5 iter + multimat nosplit
+      # Toy geom nosplit:  2 procs, ls_val=0.0 + remesh hsiz 0.1 + 4 iter + multimat nosplit
+      ## Remark : FAIL inside mmg3d scotch renum of iter 5 if niter = 5
+
       SET( NP 2 )
       add_test( NAME ls-DisIn-toygeom-nosplit-${MODE}-${NP}
         COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
         ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-distributed-${MODE}-mat-edges.mesh -v 5
         -f ${CI_DIR}/LevelSet/${NP}p_toygeom/nosplit.mmg3d
-        -hsiz 0.1 -niter 5
+        -hsiz 0.1 -niter 4
         -ls 0.0
         -sol ${CI_DIR}/LevelSet/${NP}p_toygeom/cube-ls.sol
         -out ${CI_DIR_RESULTS}/ls-DisIn-toygeom-nosplit-${MODE}-${NP}.o.mesh)
diff --git a/src/libparmmg1.c b/src/libparmmg1.c
index 7439f97a..8fb56c50 100644
--- a/src/libparmmg1.c
+++ b/src/libparmmg1.c
@@ -630,6 +630,9 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
   /** Reset the boundary fields between the old mesh size and the new one (Mmg
    * uses this fields assiming they are setted to 0)/ */
+
+  permNodGlob = NULL;
+
   for ( i=0; i<parmesh->ngrp; ++i ) {
     mesh         = parmesh->listgrp[i].mesh;
 
@@ -719,7 +722,13 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
 #ifdef USE_SCOTCH
         /* Allocation of the array that will store the node permutation */
-        PMMG_MALLOC(parmesh,permNodGlob,mesh->np+1,int,"node permutation",
+        // npi stores the number of points when we enter Mmg, np stores the
+        // number of points after adatptation.
+        // In theorie, here np == npi
+
+        assert ( mesh->np == mesh->npi );
+
+        PMMG_MALLOC(parmesh,permNodGlob,mesh->npi+1,int,"node permutation",
                     PMMG_scotch_message(&warnScotch) );
         if ( permNodGlob ) {
           for ( k=1; k<=mesh->np; ++k ) {
@@ -832,14 +841,15 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
           goto strong_failed;
         }
 
+#ifdef USE_SCOTCH
+        PMMG_DEL_MEM(parmesh,permNodGlob,int,"node permutation");
+#endif
+
         if ( !ier ) { break; }
+
       }
       /* Reset the mesh->gap field in case Mmg have modified it */
       mesh->gap = MMG5_GAP;
-
-#ifdef USE_SCOTCH
-      PMMG_DEL_MEM(parmesh,permNodGlob,int,"node permutation");
-#endif
     }
 
     MPI_Allreduce( &ier, &ieresult, 1, MPI_INT, MPI_MIN, parmesh->comm );

From d8069d04d0722056daf4c171ff38e04d767b6d55 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane.froehly@inria.fr>
Date: Wed, 30 Oct 2024 20:28:07 +0100
Subject: [PATCH 11/11] Update npi and nei after centralized ls discretization.

---
 src/libparmmg.c | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/libparmmg.c b/src/libparmmg.c
index c2d89bbe..3063cf1e 100644
--- a/src/libparmmg.c
+++ b/src/libparmmg.c
@@ -191,6 +191,10 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
     if ( !MMG3D_mmg3d2(mesh,ls,met) ) {
       return PMMG_STRONGFAILURE;
     }
+    /* Update mesh->npi and mesh->nei to be equal to mesh->np and mesh->ne, respectively */
+    mesh->npi = mesh->np;
+    mesh->nei = mesh->ne;
+
     chrono(OFF,&(ctim[tim]));
     printim(ctim[tim].gdif,stim);
     if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {