From 9adff4072c7321687dbbf744d29f2bc274d63be9 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Sun, 31 Dec 2023 15:23:18 +0100
Subject: [PATCH 1/8] applied the same changes to 5.7.1/develop as to 5.7.0:
 print all coquilface warnings, print problem location, min mesh quality in %g
 format, a few lines of documentation

---
 src/mmg3d/boulep_3d.c  | 33 ++++++++++-------
 src/mmg3d/libmmg3d.c   |  2 +
 src/mmg3d/libmmg3d.h   |  6 +++
 src/mmg3d/mmg3d.c      | 22 +++++++++++
 src/mmg3d/mmg3d1.c     | 83 +++++++++++++++++++++++++++++++-----------
 src/mmg3d/quality_3d.c |  6 ++-
 6 files changed, 114 insertions(+), 38 deletions(-)

diff --git a/src/mmg3d/boulep_3d.c b/src/mmg3d/boulep_3d.c
index bd1c16989..c47f86d95 100644
--- a/src/mmg3d/boulep_3d.c
+++ b/src/mmg3d/boulep_3d.c
@@ -1723,15 +1723,16 @@ int MMG3D_coquilFaceFirstLoop(MMG5_pMesh mesh,MMG5_int start,MMG5_int na,MMG5_in
     pradj = (*adj);
     pri    = i;
 
-    /* travel through new tetra */
+    // travel through new tetra.
+    // ier=face tag if boundary, 0 if not, -1 if error
+    // Will update *adj, *piv, iface.
     ier = MMG5_coquilTravel(mesh,na,nb,adj,piv,&iface,&i);
 
     /* fill the shell */
     list[(*ilist)] = 6*(int64_t)pradj +pri;
     (*ilist)++;
 
-    /* overflow */
-    if ( (*ilist) > MMG3D_LMAX-2 ) {
+    if ( (*ilist) > MMG3D_LMAX-2 ) {     // overflow
       if ( !mmgErr0 ) {
         fprintf(stderr,"\n  ## Warning: %s: problem in remesh process."
                 " Coquil of edge %" MMG5_PRId "-%" MMG5_PRId " contains too many elts.\n",
@@ -1743,16 +1744,11 @@ int MMG3D_coquilFaceFirstLoop(MMG5_pMesh mesh,MMG5_int start,MMG5_int na,MMG5_in
       return -1;
     }
 
-    if ( ier<0 ) return -1;
-    else if ( !ier ) continue;
+    if ( ier<0 ) return -1;      // eror
+    else if ( !ier ) continue;   // not a boundary
 
-    if ( !(*it2) ) {
-      *it2 = 4*pradj+iface;
-      (*nbdy)++;
-    }
-    else {
-      (*nbdy)++;
-    }
+    if ( !(*it2) ) *it2 = 4*pradj+iface;
+    (*nbdy)++;
 
   } while ( (*adj) && ((*adj) != start) );
 
@@ -1831,7 +1827,11 @@ void MMG3D_coquilFaceSecondLoopInit(MMG5_pMesh mesh,MMG5_int piv,int8_t *iface,
  *
  * Find all tets sharing edge \a ia of tetra \a start, and stores boundary faces
  * when met. \f$ it1 \f$ and \f$ it2 = 6*iel + iface\f$, \a iel = index of
- * tetra, \a iface = index of face in tetra.
+ * tetra, \a iface = index of face in tetra. This function can print a warning
+ * or error message when it finds that the edge has more than one boundary
+ * face. This is an error condition if the edge is supposed to be a manifold
+ * edge. Indeed this function is supposed not to be called for non-manifold
+ * edges, i.e. edges where multiple boundaries join.
  *
  * \warning Don't work if \a ia has only one boundary face in its shell.
  */
@@ -1845,6 +1845,10 @@ int MMG5_coquilface(MMG5_pMesh mesh,MMG5_int start,int8_t iface,int ia,int64_t *
 
   pt = &mesh->tetra[start];
 
+  /* MMG5_coquilface is called only on edges marked as manifold, check this */
+  assert ( pt->xt );
+  assert ( !(mesh->xtetra[pt->xt].tag[ia] & MG_NOM) );
+
   na   = pt->v[ MMG5_iare[ia][0] ];
   nb   = pt->v[ MMG5_iare[ia][1] ];
 
@@ -1883,7 +1887,8 @@ int MMG5_coquilface(MMG5_pMesh mesh,MMG5_int start,int8_t iface,int ia,int64_t *
           printf("  ## Warning: %s: you have %d boundary triangles in the closed shell"
                  " of a manifold edge.\n",__func__,nbdy);
           printf("  Problem may occur during remesh process.\n");
-          mmgWarn0 = 1;
+          MMG5_show_tet_location(mesh, pt, start);
+          if(0) mmgWarn0 = 1;     // disabled, I want to see how many there are!
 
           /* MMG5_coquilface is called only on edges marked as manifold, check this */
           assert ( pt->xt );
diff --git a/src/mmg3d/libmmg3d.c b/src/mmg3d/libmmg3d.c
index 60fa581ec..b392ccc32 100644
--- a/src/mmg3d/libmmg3d.c
+++ b/src/mmg3d/libmmg3d.c
@@ -972,6 +972,8 @@ int MMG3D_packMesh(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met) {
   return 1;
 }
 
+// Do all the work in a remeshing run (apart from input etc)
+//
 int MMG3D_mmg3dlib(MMG5_pMesh mesh,MMG5_pSol met) {
   MMG5_pSol sol=NULL; // unused
   mytime    ctim[TIMEMAX];
diff --git a/src/mmg3d/libmmg3d.h b/src/mmg3d/libmmg3d.h
index 1f3a072b9..2a3774b96 100644
--- a/src/mmg3d/libmmg3d.h
+++ b/src/mmg3d/libmmg3d.h
@@ -3057,4 +3057,10 @@ LIBMMG3D_EXPORT int MMG3D_loadVtuMesh_and_allData(MMG5_pMesh mesh,MMG5_pSol *sol
 }
 #endif
 
+
+/* --------------------------------------  Mark's hacks ---------------------------------------------  */
+
+void MMG5_show_tet_location(MMG5_pMesh mesh, MMG5_pTetra pt, int iel);
+
+
 #endif
diff --git a/src/mmg3d/mmg3d.c b/src/mmg3d/mmg3d.c
index 392dd6ae8..a2312acfc 100644
--- a/src/mmg3d/mmg3d.c
+++ b/src/mmg3d/mmg3d.c
@@ -21,6 +21,7 @@
 ** =============================================================================
 */
 
+
 /**
  * \file mmg3d/mmg3d.c
  * \brief Main file for MMG3D executable: perform 3d mesh adaptation.
@@ -538,3 +539,24 @@ int main(int argc,char *argv[]) {
   /* free mem */
   MMG5_RETURN_AND_FREE(mesh,met,ls,disp,ier);
 }
+
+
+
+/* --------------------------------------  Mark's hacks ---------------------------------------------  */
+
+// Print the indices and coordinates of the vertices of one tet.  This is to
+// inform the user about the location of a problem when the tet index is not
+// helpful, for example because it does not correpond to the input or output
+// mesh.
+//
+void MMG5_show_tet_location(MMG5_pMesh mesh, MMG5_pTetra pt, int iel)
+{
+  fprintf(stderr, "  ## tet index %d\n", iel);
+  for(int j=0; j<4; j++){
+    double U[3], *S = mesh->point[pt->v[j]].c;        // unscaled and scaled coords
+    for(int i=0; i<3; i++) U[i] = S[i]*mesh->info.delta + mesh->info.min[i];
+    fprintf(stderr, "  ## vertex %d at (%f,%f,%f)\n", pt->v[j], U[0], U[1], U[2]);
+  }
+}
+
+
diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index be8afb1e3..7d445a821 100644
--- a/src/mmg3d/mmg3d1.c
+++ b/src/mmg3d/mmg3d1.c
@@ -625,16 +625,24 @@ MMG5_int MMG5_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree, int
 
           ret = MMG5_coquilface(mesh,k,i,ia,list,&it1,&it2,0);
           ilist = ret / 2;
-          if ( ret < 0 )  return -1;
+          if ( ret < 0 ){
+            MMG5_show_tet_location(mesh, pt, k);
+            return -1;
+          }
           /* CAUTION: trigger collapse with 2 elements */
           if ( ilist <= 1 )  continue;
           ier = MMG5_chkswpbdy(mesh,met,list,ilist,it1,it2,typchk);
-          if ( ier <  0 )
+          if ( ier <  0 ){
+            MMG5_show_tet_location(mesh, pt, k);
             return -1;
+          }
           else if ( ier ) {
             ier = MMG5_swpbdy(mesh,met,list,ret,it1,PROctree,typchk);
             if ( ier > 0 )  ns++;
-            else if ( ier < 0 )  return -1;
+            else if ( ier < 0 ){
+              MMG5_show_tet_location(mesh, pt, k);
+              return -1;
+            }
             break;
           }
         }
@@ -694,11 +702,17 @@ MMG5_int MMG5_swptet(MMG5_pMesh mesh,MMG5_pSol met,double crit,double declic,
         }
 
         nconf = MMG5_chkswpgen(mesh,met,k,i,&ilist,list,crit,typchk);
-        if ( nconf<0 ) return -1;
+        if ( nconf<0 ){
+          MMG5_show_tet_location(mesh, pt, k);
+          return -1;
+        }
         else if ( nconf ) {
           ier = MMG5_swpgen(mesh,met,nconf,ilist,list,PROctree,typchk);
           if ( ier > 0 )  ns++;
-          else if ( ier < 0 ) return -1;
+          else if ( ier < 0 ){
+            MMG5_show_tet_location(mesh, pt, k);
+            return -1;
+          }
           break;
         }
       }
@@ -807,8 +821,10 @@ MMG5_int MMG5_movtet(MMG5_pMesh mesh,MMG5_pSol met, MMG3D_pPROctree PROctree,
                 if( !ier )  continue;
                 else if ( ier>0 )
                   ier = MMG5_movbdynompt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
-                else
+                else{
+                  MMG5_show_tet_location(mesh, pt, k);
                   return -1;
+                }
               }
             }
             else if ( ppt->tag & MG_GEO ) {
@@ -816,8 +832,10 @@ MMG5_int MMG5_movtet(MMG5_pMesh mesh,MMG5_pSol met, MMG3D_pPROctree PROctree,
               if ( !ier )  continue;
               else if ( ier>0 )
                 ier = MMG5_movbdyridpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
-              else
+              else{
+                MMG5_show_tet_location(mesh, pt, k);
                 return -1;
+              }
             }
             else if ( ppt->tag & MG_REF ) {
               ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
@@ -825,16 +843,19 @@ MMG5_int MMG5_movtet(MMG5_pMesh mesh,MMG5_pSol met, MMG3D_pPROctree PROctree,
                 continue;
               else if ( ier>0 )
                 ier = MMG5_movbdyrefpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
-              else
+              else{
+                MMG5_show_tet_location(mesh, pt, k);
                 return -1;
+              }
             }
             else {
               ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
               if ( !ier )
                 continue;
-              else if ( ier<0 )
+              else if ( ier<0 ){
+                MMG5_show_tet_location(mesh, pt, k);
                 return -1;
-
+              }
               n = &(mesh->xpoint[ppt->xp].n1[0]);
 
               /* If the orientation of the tetra face is
@@ -846,7 +867,10 @@ MMG5_int MMG5_movtet(MMG5_pMesh mesh,MMG5_pSol met, MMG3D_pPROctree PROctree,
               }
               ier = MMG5_movbdyregpt(mesh,met,PROctree,listv,ilistv,
                                      lists,ilists,improveSurf,improveVolSurf);
-              if (ier < 0 ) return -1;
+              if (ier < 0 ){
+                MMG5_show_tet_location(mesh, pt, k);
+                return -1;
+              }
               else if ( ier )  ns++;
             }
           }
@@ -964,14 +988,18 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               else {
                 if ( mesh->adja[4*(k-1)+1+i] )  continue;
                 if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
-                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 )
+                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 ){
+                  MMG5_show_tet_location(mesh, pt, k);
                   return -1;
+                }
               }
             }
             else {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
-                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
+                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
+                MMG5_show_tet_location(mesh, pt, k);
                 return -1;
+              }
             }
           }
           else {
@@ -1095,14 +1123,18 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               else {
                 if ( mesh->adja[4*(k-1)+1+i] )  continue;
                 if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
-                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 )
+                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 ){
+                  MMG5_show_tet_location(mesh, pt, k);
                   return -1;
+                }
               }
             }
             else {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
-                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
+                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
+                MMG5_show_tet_location(mesh, pt, k);
                 return -1;
+              }
             }
           }
           else {
@@ -1134,13 +1166,19 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
 
         if ( ilist > 0 ) {
           ier = MMG5_colver(mesh,met,list,ilist,iq,typchk);
-          if ( ier < 0 ) return -1;
+          if ( ier < 0 ){
+            MMG5_show_tet_location(mesh, pt, k);
+            return -1;
+          }
           else if ( ier ) {
             MMG3D_delPt(mesh,ier);
             break;
           }
         }
-        else if (ilist < 0 ) return -1;
+        else if (ilist < 0 ){
+          MMG5_show_tet_location(mesh, pt, k);
+          return -1;
+        }
       }
       if ( ier ) {
         p1->flag = base;
@@ -2593,6 +2631,7 @@ MMG3D_anatets_iso(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               fprintf(stderr,"\n  ## Warning: %s: surfacic pattern: unable to find"
                       " a valid split for at least 1 point. Point(s) deletion.\n",
                       __func__ );
+              MMG5_show_tet_location(mesh, pt, k);
               mmgWarn2 = 1;
             }
 
@@ -3161,7 +3200,7 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
       ier = MMG3D_anatets(mesh,met,typchk);
 
       if ( ier < 0 ) {
-        fprintf(stderr,"\n  ## Unable to complete surface mesh. Exit program.\n");
+        fprintf(stderr,"\n  ## Unable to complete surface mesh in MMG5_anatet. Exit program.\n");
         return 0;
       }
       ns += ier;
@@ -3172,7 +3211,7 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
         /* analyze internal tetras */
         ier = MMG5_anatetv(mesh,met,typchk);
         if ( ier < 0 ) {
-          fprintf(stderr,"\n  ## Unable to complete volume mesh. Exit program.\n");
+          fprintf(stderr,"\n  ## Unable to complete volume mesh in MMG5_anatet. Exit program.\n");
           return 0;
         }
         ns += ier;
@@ -3188,7 +3227,7 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
     if ( !mesh->info.noinsert ) {
       nc = MMG5_coltet(mesh,met,typchk);
       if ( nc < 0 ) {
-        fprintf(stderr,"\n  ## Unable to collapse mesh. Exiting.\n");
+        fprintf(stderr,"\n  ## Unable to collapse mesh in MMG5_anatet: nc=%d. Exiting.\n", nc);
         return 0;
       }
     }
@@ -3198,14 +3237,14 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
     if ( !mesh->info.noswap ) {
       ier = MMG5_swpmsh(mesh,met,NULL,typchk);
       if ( ier < 0 ) {
-        fprintf(stderr,"\n  ## Unable to improve mesh. Exiting.\n");
+        fprintf(stderr,"\n  ## Unable to improve mesh in MMG5_anatet. Exiting.\n");
         return 0;
       }
       nf  += ier;
 
       ier = MMG5_swptet(mesh,met,MMG3D_LSWAPIMPROVE,MMG3D_SWAP06,NULL,typchk,mesh->mark-2);
       if ( ier < 0 ) {
-        fprintf(stderr,"\n  ## Unable to improve mesh. Exiting.\n");
+        fprintf(stderr,"\n  ## Unable to improve mesh in MMG5_anatet. Exiting.\n");
         return 0;
       }
       nf += ier;
diff --git a/src/mmg3d/quality_3d.c b/src/mmg3d/quality_3d.c
index ba82544d6..7ce9679de 100644
--- a/src/mmg3d/quality_3d.c
+++ b/src/mmg3d/quality_3d.c
@@ -91,7 +91,9 @@ int MMG3D_tetraQual(MMG5_pMesh mesh, MMG5_pSol met,int8_t metRidTyp) {
 
   /* Here the quality is not normalized by alpha, thus we need to
    * normalized it */
-  return  MMG5_minQualCheck(iel,minqual,MMG3D_ALPHAD);
+  int r = MMG5_minQualCheck(iel,minqual,MMG3D_ALPHAD);
+  if(r==0) MMG5_show_tet_location(mesh, pt, iel);
+  return  r;
 }
 
 /**
@@ -476,7 +478,7 @@ int MMG3D_displayQualHisto(MMG5_int ne,double max,double avg,double min,MMG5_int
     fprintf(stdout," (LES)");
   fprintf(stdout,"  %" MMG5_PRId "\n",ne);
 
-  fprintf(stdout,"     BEST   %8.6f  AVRG.   %8.6f  WRST.   %8.6f (%" MMG5_PRId ")\n",
+  fprintf(stdout,"     BEST   %8.6f  AVRG.   %8.6f  WRST.   %.6g (%" MMG5_PRId ")\n",
           max,avg / ne,min,iel);
 
   return ( MMG3D_displayQualHisto_internal(ne,max,avg,min,iel,good,med,his,

From 219af1fa31530e4156e2d43f3efe481e36980af4 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Thu, 25 Jan 2024 18:24:33 +0100
Subject: [PATCH 2/8] In MMG5_coltet(), if MMG5_boulesurfvolpNom() fails, then
 jump to the next tet instead of returning to the caller. Also changed more
 return values so I can see where errors come from.

---
 src/mmg3d/boulep_3d.c |  8 ++++----
 src/mmg3d/mmg3d1.c    | 30 +++++++++++++++++++-----------
 2 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/src/mmg3d/boulep_3d.c b/src/mmg3d/boulep_3d.c
index c47f86d95..a8522f650 100644
--- a/src/mmg3d/boulep_3d.c
+++ b/src/mmg3d/boulep_3d.c
@@ -840,13 +840,13 @@ int MMG5_boulesurfvolpNom(MMG5_pMesh mesh,MMG5_int start,int ip,int iface,
           if ( *refplus == -1 ) {
             if ( pt->ref != *refmin  ) *refplus = pt->ref;
           }
-          else if ( pt->ref != *refmin && pt->ref != *refplus ) return -1;
+          else if ( pt->ref != *refmin && pt->ref != *refplus ) return -2;
         }
         pt->flag = base;
       }
       
       /* identification of edge number in tetra k */
-      if ( !MMG3D_findEdge(mesh,pt,k,na,nb,0,&mmgErr2,&i) ) return -1;
+      if ( !MMG3D_findEdge(mesh,pt,k,na,nb,0,&mmgErr2,&i) ) return -3;
       
       /* set sense of travel */
       if ( pt->v[ MMG5_ifar[i][0] ] == piv ) {
@@ -906,7 +906,7 @@ int MMG5_boulesurfvolpNom(MMG5_pMesh mesh,MMG5_int start,int ip,int iface,
                   " or/and the maximum mesh.\n");
           mmgErr1 = 1;
         }
-        return -1;
+        return -4;
       }
       listv[(*ilistv)] = 4*k1+j;
       (*ilistv)++;
@@ -918,7 +918,7 @@ int MMG5_boulesurfvolpNom(MMG5_pMesh mesh,MMG5_int start,int ip,int iface,
         if ( *refplus == -1 ) {
           if ( pt1->ref != *refmin  ) *refplus = pt1->ref;
         }
-        else if ( pt1->ref != *refmin && pt1->ref != *refplus ) return -1;
+        else if ( pt1->ref != *refmin && pt1->ref != *refplus ) return -5;
       }
     }
     cur++;
diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index 7d445a821..bab88e1d6 100644
--- a/src/mmg3d/mmg3d1.c
+++ b/src/mmg3d/mmg3d1.c
@@ -974,8 +974,8 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
         if ( mesh->info.npar ) {
           if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
             tag = pxt->tag[MMG5_iarf[i][j]];
-            isnm = (tag & MG_NOM);
-            isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor );
+            isnm = (tag & MG_NOM);             // is non-manifold
+            isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor ); // is non-manifold interior
 
             if ( p0->tag > tag ) continue;
 
@@ -987,8 +987,10 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               }
               else {
                 if ( mesh->adja[4*(k-1)+1+i] )  continue;
-                if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
-                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 ){
+                int bsret = MMG5_boulesurfvolpNom(mesh,k,ip,i, list,&ilist,lists,&ilists,
+                                                  &refmin,&refplus,p0->tag & MG_NOM);
+                if(bsret < 0 ){
+                  printf("  ## (1) MMG5_boulesurfvolpNom returned %d\n", bsret);
                   MMG5_show_tet_location(mesh, pt, k);
                   return -1;
                 }
@@ -998,7 +1000,7 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
                                      list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
                 MMG5_show_tet_location(mesh, pt, k);
-                return -1;
+                return -2;
               }
             }
           }
@@ -1122,10 +1124,15 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               }
               else {
                 if ( mesh->adja[4*(k-1)+1+i] )  continue;
-                if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
-                                          list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 ){
+                int bsret = MMG5_boulesurfvolpNom(mesh,k,ip,i,list,&ilist,lists,&ilists,
+                                                  &refmin,&refplus,p0->tag & MG_NOM);
+                if(bsret < 0 ){
+                  printf("  ## (2) MMG5_boulesurfvolpNom returned %d\n", bsret);
                   MMG5_show_tet_location(mesh, pt, k);
-                  return -1;
+                  // HACK: jump to the next tet. Maybe a "continue" would be
+                  // good enough (to go to the next edge of the face).
+                  // The original code did "return -1".
+                  goto next;    // return -3;
                 }
               }
             }
@@ -1133,7 +1140,7 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
                                      list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
                 MMG5_show_tet_location(mesh, pt, k);
-                return -1;
+                return -4;
               }
             }
           }
@@ -1168,7 +1175,7 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
           ier = MMG5_colver(mesh,met,list,ilist,iq,typchk);
           if ( ier < 0 ){
             MMG5_show_tet_location(mesh, pt, k);
-            return -1;
+            return -5;
           }
           else if ( ier ) {
             MMG3D_delPt(mesh,ier);
@@ -1177,7 +1184,7 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
         }
         else if (ilist < 0 ){
           MMG5_show_tet_location(mesh, pt, k);
-          return -1;
+          return -6;
         }
       }
       if ( ier ) {
@@ -1187,6 +1194,7 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
         break;
       }
     }
+  next: ;        // jump here if giving up on the current edge
   }
   if ( nc > 0 && (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) )
     fprintf(stdout,"     %8" MMG5_PRId " vertices removed, %8" MMG5_PRId " non manifold,\n",nc,nnm);

From b1444c56139c77e3d6d977da0c0eb9cc063c90f9 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Mon, 29 Jan 2024 13:38:28 +0100
Subject: [PATCH 3/8] neater solution for MMG5_coltet(): we can continue
 instead of goto, and

distinguish harmful from harmless returns from  MMG5_boulesurfvolpNom()
---
 src/mmg3d/mmg3d1.c | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index bab88e1d6..acce86ea6 100644
--- a/src/mmg3d/mmg3d1.c
+++ b/src/mmg3d/mmg3d1.c
@@ -1129,10 +1129,16 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
                 if(bsret < 0 ){
                   printf("  ## (2) MMG5_boulesurfvolpNom returned %d\n", bsret);
                   MMG5_show_tet_location(mesh, pt, k);
-                  // HACK: jump to the next tet. Maybe a "continue" would be
-                  // good enough (to go to the next edge of the face).
-                  // The original code did "return -1".
-                  goto next;    // return -3;
+                  // Mark's modification. The original code in mmg-5.7.1
+                  // returned -1 here. Algiane confirmed that the cases where
+                  // MMG5_boulesurfvolpNom returns anything but -1 or -4 (in my
+                  // version) are harmless; a ball could not be computed but it
+                  // does not mean that the mesh is wrong.
+                  if(bsret==-1 || bsret==-4){
+                    return -3;
+                  }else{
+                    continue;
+                  }
                 }
               }
             }
@@ -1194,7 +1200,6 @@ static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
         break;
       }
     }
-  next: ;        // jump here if giving up on the current edge
   }
   if ( nc > 0 && (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) )
     fprintf(stdout,"     %8" MMG5_PRId " vertices removed, %8" MMG5_PRId " non manifold,\n",nc,nnm);

From 67d98701abba45cbedbd38db0e2cb251c37712d5 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Mon, 29 Jan 2024 13:44:31 +0100
Subject: [PATCH 4/8] comment on the "we should rarely passed here" warning

---
 src/mmg3d/mmg3d2.c | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/src/mmg3d/mmg3d2.c b/src/mmg3d/mmg3d2.c
index b10da60cb..161b85047 100644
--- a/src/mmg3d/mmg3d2.c
+++ b/src/mmg3d/mmg3d2.c
@@ -2147,9 +2147,21 @@ int MMG3D_chkmanicoll(MMG5_pMesh mesh,MMG5_int k,int iface,int iedg,MMG5_int nde
          * non manifold */
         if ( indq == -1 ) {
           if ( mesh->info.ddebug ) {
-          fprintf(stderr,"\n  ## Warning: %s: we should rarely passed here. "
-                  "tetra %" MMG5_PRId " =  %" MMG5_PRId " %" MMG5_PRId " %" MMG5_PRId " %" MMG5_PRId ", ref = %" MMG5_PRId ".",__func__,
-                  jel,pt1->v[0],pt1->v[1],pt1->v[2],pt1->v[3],pt1->ref);
+            //
+            // Algiane about these (wrt the "david" test cases in MICROCARD/WP7:
+            // "ça n’est pas grave du tout : ça arrive quand on prévoit que si
+            // on collapse un point on créera une situation non-manifold. Du
+            // coup on ne collapse pas. On a un warning « historique » car si ça
+            // arrive vraiment beaucoup, on peut imaginer que le maillage final
+            // ne sera pas très bon et pour des cas d’optimisation de forme on
+            // peut éviter ces situations. Sur tes géométries ça me semble
+            // évident qu’on a pas mal de cas où collapser créerait des
+            // situations non-manifold donc ça n’est pas surprenant de voir ce
+            // warning."
+            //
+            fprintf(stderr,"\n  ## Warning: %s: we should rarely passed here. "
+                    "tetra %" MMG5_PRId " =  %" MMG5_PRId " %" MMG5_PRId " %" MMG5_PRId " %" MMG5_PRId ", ref = %" MMG5_PRId ".",__func__,
+                    jel,pt1->v[0],pt1->v[1],pt1->v[2],pt1->v[3],pt1->ref);
           }
           return 0;
         }

From 0c583b2ecfbdbf76e4b98b82458e5f1cffbcfebf Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Tue, 6 Feb 2024 16:59:39 +0100
Subject: [PATCH 5/8] moved implementation of my  MMG5_show_tet_location() to
 libmmg3d.c so that it is also available in the library version

---
 src/mmg3d/libmmg3d.c | 22 ++++++++++++++++++++++
 src/mmg3d/mmg3d.c    | 18 ------------------
 2 files changed, 22 insertions(+), 18 deletions(-)

diff --git a/src/mmg3d/libmmg3d.c b/src/mmg3d/libmmg3d.c
index b392ccc32..f68226bf3 100644
--- a/src/mmg3d/libmmg3d.c
+++ b/src/mmg3d/libmmg3d.c
@@ -1757,3 +1757,25 @@ void MMG3D_Set_commonFunc(void) {
     MMG5_renumbering = MMG5_mmg3dRenumbering;
 #endif
 }
+
+
+
+
+/* --------------------------------  Mark's hacks --------------------------------------  */
+
+// Print the indices and coordinates of the vertices of one tet.  This is to
+// inform the user about the location of a problem when the tet index is not
+// helpful, for example because it does not correpond to the input or output
+// mesh.
+//
+void MMG5_show_tet_location(MMG5_pMesh mesh, MMG5_pTetra pt, int iel)
+{
+  if ( mesh->info.imprim > 0 ){
+    fprintf(stderr, "  ## tet index %d\n", iel);
+    for(int j=0; j<4; j++){
+      double U[3], *S = mesh->point[pt->v[j]].c;        // unscaled and scaled coords
+      for(int i=0; i<3; i++) U[i] = S[i]*mesh->info.delta + mesh->info.min[i];
+      fprintf(stderr, "  ## vertex %d at (%f,%f,%f)\n", pt->v[j], U[0], U[1], U[2]);
+    }
+  }
+}
diff --git a/src/mmg3d/mmg3d.c b/src/mmg3d/mmg3d.c
index a2312acfc..4e186f607 100644
--- a/src/mmg3d/mmg3d.c
+++ b/src/mmg3d/mmg3d.c
@@ -542,21 +542,3 @@ int main(int argc,char *argv[]) {
 
 
 
-/* --------------------------------------  Mark's hacks ---------------------------------------------  */
-
-// Print the indices and coordinates of the vertices of one tet.  This is to
-// inform the user about the location of a problem when the tet index is not
-// helpful, for example because it does not correpond to the input or output
-// mesh.
-//
-void MMG5_show_tet_location(MMG5_pMesh mesh, MMG5_pTetra pt, int iel)
-{
-  fprintf(stderr, "  ## tet index %d\n", iel);
-  for(int j=0; j<4; j++){
-    double U[3], *S = mesh->point[pt->v[j]].c;        // unscaled and scaled coords
-    for(int i=0; i<3; i++) U[i] = S[i]*mesh->info.delta + mesh->info.min[i];
-    fprintf(stderr, "  ## vertex %d at (%f,%f,%f)\n", pt->v[j], U[0], U[1], U[2]);
-  }
-}
-
-

From d270498022d8f6b740a38de2ef6349cdac3d0467 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Fri, 10 May 2024 09:06:30 +0200
Subject: [PATCH 6/8] Don't erase MG_REQ in MMG5_bdryTria()   This reverses one
 line of commit ba87623ffc4bf758708686cef76814d47878b381   where it was argued
 that this is a "parallel tag" which is to be removed   from faces that are
 not parallel faces. It may have been necessary for   ParMmg but it erases all
 RequiredTriangles from mmg3d output, and   I rely on them in the "mmg3d -ls
 ..." calls in WP7/synthetic/mason.sh

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

diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c
index 0f15bbf7a..2b8a0141c 100644
--- a/src/mmg3d/hash_3d.c
+++ b/src/mmg3d/hash_3d.c
@@ -1382,7 +1382,7 @@ int MMG5_bdryTria(MMG5_pMesh mesh, MMG5_int ntmesh) {
               if ( !(pxt->ftag[i] & MG_PARBDY)) {
                 ptt->tag[j] &= ~MG_PARBDY;
                 ptt->tag[j] &= ~MG_NOSURF;
-                ptt->tag[j] &= ~MG_REQ;
+                // ptt->tag[j] &= ~MG_REQ;           // FIXME
               }
             }
             /* Assign ref to tria from xtetra->edg */

From a90333bfcfced1c63f8879bd72cfc367ec030261 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Sat, 11 May 2024 17:47:45 +0200
Subject: [PATCH 7/8] added some code in mmg3d2.c to investigate LS snapping
 behaviour

---
 src/mmg3d/mmg3d2.c | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/src/mmg3d/mmg3d2.c b/src/mmg3d/mmg3d2.c
index 41c374a3b..4b7883487 100644
--- a/src/mmg3d/mmg3d2.c
+++ b/src/mmg3d/mmg3d2.c
@@ -701,7 +701,7 @@ int MMG3D_snpval_ls(MMG5_pMesh mesh,MMG5_pSol sol) {
       if ( i < 4 ) {
         for (i=0; i<4; i++) {
           ip = pt->v[i];
-          sol->m[ip] = -1000.0*MMG5_EPS;
+          sol->m[ip] = -1001.0*MMG5_EPS;
         }
       }
     }
@@ -718,7 +718,7 @@ int MMG3D_snpval_ls(MMG5_pMesh mesh,MMG5_pSol sol) {
                 "previous value: %E.\n",__func__,k,fabs(sol->m[k]));
 
       tmp[k] = ( fabs(sol->m[k]) < MMG5_EPSD ) ?
-        (-100.0*MMG5_EPS) : sol->m[k];
+        (-111.0*MMG5_EPS) : sol->m[k];
       p0->flag = 1;
       sol->m[k] = 0;
       ns++;
@@ -738,9 +738,9 @@ int MMG3D_snpval_ls(MMG5_pMesh mesh,MMG5_pSol sol) {
         if ( p0->flag == 1 ) {
           if ( !MMG3D_ismaniball(mesh,sol,k,i) ) {
             if ( tmp[ip] < 0.0 )
-              sol->m[ip] = -100.0*MMG5_EPS;
+              sol->m[ip] = -123.0*MMG5_EPS;
             else
-              sol->m[ip] = +100.0*MMG5_EPS;
+              sol->m[ip] = +123.0*MMG5_EPS;
 
             p0->flag = 0;
             nc++;
@@ -880,7 +880,7 @@ int MMG3D_rmc(MMG5_pMesh mesh, MMG5_pSol sol){
         for (i=0; i<4; i++) {
           ip0 = pt1->v[i];
           if ( sol->m[ip0] > 0.0 ) {
-            sol->m[ip0] = - 100*MMG5_EPS;
+            sol->m[ip0] = - 101*MMG5_EPS;
           }
         }
       }
@@ -960,7 +960,7 @@ int MMG3D_rmc(MMG5_pMesh mesh, MMG5_pSol sol){
         pt1 = &mesh->tetra[pile[l]];
         for (i=0; i<4; i++) {
           ip0 = pt1->v[i];
-          if ( sol->m[ip0] < 0.0 ) sol->m[ip0] = 100*MMG5_EPS;
+          if ( sol->m[ip0] < 0.0 ) sol->m[ip0] = 101*MMG5_EPS;
         }
       }
       ncm++;
@@ -993,7 +993,7 @@ int MMG3D_rmc(MMG5_pMesh mesh, MMG5_pSol sol){
           pt1 = &mesh->tetra[pile[l]];
           for (i=0; i<4; i++) {
             ip0 = pt1->v[i];
-            if ( sol->m[ip0] < 0.0 ) sol->m[ip0] = 100*MMG5_EPS;
+            if ( sol->m[ip0] < 0.0 ) sol->m[ip0] = 102*MMG5_EPS;
           }
         }
         ncm++;
@@ -1121,6 +1121,15 @@ int MMG3D_cuttet_ls(MMG5_pMesh mesh, MMG5_pSol sol,MMG5_pSol met){
       else if ( !p0->flag || !p1->flag )  continue;
 
       npneg = (np<0);
+      if(npneg){
+        double X[3];
+        printf("splitting a required edge:\n");
+        for(int i=0; i<3; i++) X[i] = p0->c[i]*mesh->info.delta + mesh->info.min[i];
+        printf("p0 index %6d at (%f, %f, %f) sol=%8g\n", ip0, X[0], X[1], X[2], v0);
+        for(int i=0; i<3; i++) X[i] = p1->c[i]*mesh->info.delta + mesh->info.min[i];
+        printf("p1 index %6d at (%f, %f, %f) sol=%8g\n", ip1, X[0], X[1], p1->c[2], v1);
+        printf("\n");
+      }
 
       s = v0 / (v0-v1);
 

From 56af1a190aafb054da29767e60396a501e8266b7 Mon Sep 17 00:00:00 2001
From: Mark Potse <mark@potse.nl>
Date: Mon, 13 May 2024 12:30:43 +0200
Subject: [PATCH 8/8] get rid of some compiler warnings on fgets (and skip
 comment lines longer than 128 chars correctly)

---
 src/mmg2d/inout_2d.c | 6 +++++-
 src/mmg3d/inout_3d.c | 6 +++++-
 src/mmgs/inout_s.c   | 6 +++++-
 3 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/src/mmg2d/inout_2d.c b/src/mmg2d/inout_2d.c
index 1d311d3a0..0ab90038b 100644
--- a/src/mmg2d/inout_2d.c
+++ b/src/mmg2d/inout_2d.c
@@ -84,7 +84,11 @@ int MMG2D_loadMesh(MMG5_pMesh mesh,const char *filename) {
     strcpy(chaine,"D");
     while(fscanf(inm,"%127s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if ( chaine[0] == '#' ) {
-        fgets(strskip,MMG5_FILESTR_LGTH,inm);
+        while(1){           // skip until end of line or file
+          char *s = fgets(strskip,MMG5_FILESTR_LGTH,inm);
+          if(!s) break;     // nothing could be read
+          if(s[strlen(s)-1]=='\n') break;   // end of line
+        }
         continue;
       }
 
diff --git a/src/mmg3d/inout_3d.c b/src/mmg3d/inout_3d.c
index d30096be6..a436afa17 100644
--- a/src/mmg3d/inout_3d.c
+++ b/src/mmg3d/inout_3d.c
@@ -145,7 +145,11 @@ int MMG3D_loadMesh_opened(MMG5_pMesh mesh,FILE *inm,int bin) {
     strcpy(chaine,"D");
     while(fscanf(inm,"%127s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if ( chaine[0] == '#' ) {
-        fgets(strskip,MMG5_FILESTR_LGTH,inm);
+        while(1){           // skip until end of line or file
+          char *s = fgets(strskip,MMG5_FILESTR_LGTH,inm);
+          if(!s) break;     // nothing could be read
+          if(s[strlen(s)-1]=='\n') break;   // end of line
+        }
         continue;
       }
 
diff --git a/src/mmgs/inout_s.c b/src/mmgs/inout_s.c
index 41c6db6d2..7068eac7b 100644
--- a/src/mmgs/inout_s.c
+++ b/src/mmgs/inout_s.c
@@ -101,7 +101,11 @@ int MMGS_loadMesh(MMG5_pMesh mesh, const char *filename) {
     strcpy(chaine,"D");
     while(fscanf(inm,"%127s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if ( chaine[0] == '#' ) {
-        fgets(strskip,MMG5_FILESTR_LGTH,inm);
+        while(1){           // skip until end of line or file
+          char *s = fgets(strskip,MMG5_FILESTR_LGTH,inm);
+          if(!s) break;     // nothing could be read
+          if(s[strlen(s)-1]=='\n') break;   // end of line
+        }
         continue;
       }
       if(!strncmp(chaine,"MeshVersionFormatted",strlen("MeshVersionFormatted"))) {