Increased diagnostic output and added a few lines of documentation

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;
       }
 

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 );

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 */

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;
       }
 

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]);
+    }
+  }
+}

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

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);
 }
+
+
+