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/boulep_3d.c b/src/mmg3d/boulep_3d.c
index 94eb81d9a..90234fff9 100644
--- a/src/mmg3d/boulep_3d.c
+++ b/src/mmg3d/boulep_3d.c
@@ -843,13 +843,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 ) {
@@ -909,7 +909,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)++;
@@ -921,7 +921,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++;
@@ -1726,15 +1726,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",
@@ -1746,16 +1747,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) );
 
@@ -1834,7 +1830,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.
  */
@@ -1848,6 +1848,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] ];
 
@@ -1886,7 +1890,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/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 */
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/mmg3d/libmmg3d.c b/src/mmg3d/libmmg3d.c
index aab14438e..be9c6643a 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];
@@ -1755,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/libmmg3d.h b/src/mmg3d/libmmg3d.h
index d80c3d282..7e923cc04 100644
--- a/src/mmg3d/libmmg3d.h
+++ b/src/mmg3d/libmmg3d.h
@@ -3441,4 +3441,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 863ea5444..b09c2ba0a 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,6 @@ int main(int argc,char *argv[]) {
   /* free mem */
   MMG5_RETURN_AND_FREE(mesh,met,ls,disp,ier);
 }
+
+
+
diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index 9cb8aef26..c99ce3d27 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++;
             }
           }
@@ -950,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;
 
@@ -963,15 +987,21 @@ 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;
+                }
               }
             }
             else {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
-                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
-                return -1;
+                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
+                MMG5_show_tet_location(mesh, pt, k);
+                return -2;
+              }
             }
           }
           else {
@@ -1094,15 +1124,30 @@ 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 )
-                  return -1;
+                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);
+                  // 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;
+                  }
+                }
               }
             }
             else {
               if (MMG5_boulesurfvolp(mesh,k,ip,i,
-                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
-                return -1;
+                                     list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 ){
+                MMG5_show_tet_location(mesh, pt, k);
+                return -4;
+              }
             }
           }
           else {
@@ -1134,13 +1179,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 -5;
+          }
           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 -6;
+        }
       }
       if ( ier ) {
         p1->flag = base;
@@ -2593,6 +2644,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 +3213,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 +3224,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 +3240,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 +3250,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/mmg3d2.c b/src/mmg3d/mmg3d2.c
index 107763935..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);
 
@@ -2147,9 +2156,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;
         }
diff --git a/src/mmg3d/quality_3d.c b/src/mmg3d/quality_3d.c
index 40aa7732e..1bc9ac926 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,
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"))) {