Merge pull request #237 from MmgTools/feature/bdy-edge-interp

Feature/bdy edge interp

src/mmg3d/intmet_3d.c

@@ -47,6 +47,12 @@
  * Interpolation of anisotropic sizemap at parameter \a s along edge \a i of elt
  * \a k for a special storage of ridges metric (after defsiz call).
  *
+ * \remark for boundary edges this function has to be called from a boundary
+ * face to ensure that the boundary edge will not be interpreted as an internal
+ * edge by error. It is due to the fact that regular boundary edges are not
+ * always marked as boundary inside a xtetra and tetra with boundary edges but
+ * no boundary faces are not always associated to an xtetra (moreover, for edges
+ * not belonging to a face inside the xtetra, the tag is not always up to date).
  */
 int MMG5_intmet_ani(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t i,MMG5_int ip,
                       double s) {
@@ -73,16 +79,35 @@ int MMG5_intmet_ani(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t i,MMG5_int i
     else if ( pxt->tag[i] & MG_BDY ) {
      return MMG5_intregmet(mesh,met,k,i,s,m);
     }
-    else {
-      /* The edge is an internal edge. */
-      return MMG5_intvolmet(mesh,met,k,i,s,m);
-    }
   }
-  else {
-    /* The edge is an internal edge. */
-    return MMG5_intvolmet(mesh,met,k,i,s,m);
+
+  /* If we pass here, either we don't have an xtetra, or the edge is not marked
+   * as MG_BDY */
+  assert ( (!pt->xt) || !(pxt->tag[i] & MG_BDY) );
+
+#ifndef NDEBUG
+  /* Assert that we are not dealing by error with a boundary edge (we may have
+   * regular edges without MG_BDY tag or edges ) */
+
+  // If we come from anatets/v, mesh->adja is not allocated and we cannot called
+  // the get_shellEdgeTag function: in this case, the MG_BDY tag has been
+  // explicitely added if we arrive from a boundary face and a boundary edge
+  // should not be splitted from a non boundary face
+  if ( mesh->adja ) {
+    int16_t  tag = 0;
+    MMG5_int ref = 0;
+    if ( !MMG3D_get_shellEdgeTag(mesh,k,i,&tag,&ref) ) {
+      fprintf(stderr,"\n  ## Warning: %s: 0. unable to get edge info"
+              " (tetra %d).\n",__func__,MMG3D_indElt(mesh,k));
+      return 0;
+    }
+    assert ( (!(tag & MG_BDY)) && "Boundary edge is seen as internal");
   }
-  return 0;
+#endif
+
+  /* The edge is an internal edge. */
+  return MMG5_intvolmet(mesh,met,k,i,s,m);
+
 }
 
 /**

src/mmg3d/mmg3d1.c

@@ -618,6 +618,8 @@ MMG5_int MMG5_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree, int
         if ( !(pxt->ftag[i] & MG_BDY) ) continue;
         for (j=0; j<3; j++) {
           ia  = MMG5_iarf[i][j];
+          /* Mark the edge as boundary in case of missing tag */
+          pxt->tag[ia] |= MG_BDY;
 
           /* No swap of geometric edge */
           if ( MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) )
@@ -628,6 +630,8 @@ MMG5_int MMG5_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree, int
           if ( ret < 0 )  return -1;
           /* CAUTION: trigger collapse with 2 elements */
           if ( ilist <= 1 )  continue;
+
+          /* Here, we work on a boundary edge lying along a boundary face */
           ier = MMG5_chkswpbdy(mesh,met,list,ilist,it1,it2,typchk);
           if ( ier <  0 )
             return -1;
@@ -1947,7 +1951,14 @@ int MMG3D_splsurfedge( MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,
   int16_t      tag;
   int8_t       j,i,i1,i2;
 
-  assert ( pxt == &mesh->xtetra[pt->xt] );
+  assert ( pxt == &mesh->xtetra[pt->xt] && "suitable xtetra assignation" );
+
+  assert ( ((pxt->ftag[MMG5_ifar[imax][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[imax][1]] & MG_BDY) )
+           && "Boundary edge has to be splitted from a boundary face" );
+
+  /* Mark the edge as MG_BDY to avoid wrong evaluation as an internal edge by
+   * the interpolation function (see intmet_ani) */
+  pxt->tag[imax] |= MG_BDY;
 
   /* proceed edges according to lengths */
   MMG3D_find_bdyface_from_edge(mesh,pt,imax,&i,&j,&i1,&i2,&ip1,&ip2,&p0,&p1);
@@ -2359,6 +2370,10 @@ MMG3D_anatets_iso(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
           ppt = &mesh->point[ip];
 
           assert ( met );
+
+          // Add MG_BDY tag before interpolation
+          pxt->tag[ia] |= MG_BDY;
+
           if ( met->m ) {
             if ( typchk == 1 && (met->size>1) )
               ier = MMG3D_intmet33_ani(mesh,met,k,ia,ip,0.5);
@@ -2429,7 +2444,7 @@ MMG3D_anatets_iso(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
         }
         else if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
           /* Point at the interface of 2 boundary faces belonging to different
-           * tetra : Point has alredy been created from another tetra so we have
+           * tetra : Point has already been created from another tetra so we have
            * to store the tangent and the second normal at edge */
           ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
           assert(ier);

src/mmg3d/mmg3d1_delone.c

@@ -146,6 +146,10 @@ int MMG3D_mmg3d1_delone_split(MMG5_pMesh mesh, MMG5_pSol met,
     }
 
     ier = 1;
+
+    /* Mark edge as bdy to avoid issue in intmet */
+    pxt->tag[imax] |= MG_BDY;
+
     if ( met && met->m ) {
       ier = MMG5_intmet(mesh,met,k,imax,ip,0.5);
     }

src/mmg3d/mmg3d1_pattern.c

@@ -116,7 +116,7 @@ static MMG5_int MMG5_adpspl(MMG5_pMesh mesh,MMG5_pSol met, int* warn) {
     else {
       /* Case of an internal face */
 
-      /* Skip only boundary edges but try to trat internal edges connecting bdy
+      /* Skip only boundary edges but try to treat internal edges connecting bdy
        * points */
       int8_t isbdy;
       ilist = MMG5_coquil(mesh,k,imax,list,&isbdy);

src/mmg3d/optbdry_3d.c

@@ -321,6 +321,9 @@ int MMG3D_optbdry(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree,MMG5_in
     for (j=0; j<3; j++) {
       ia  = MMG5_iarf[i][j];
 
+      /* Mark the edge as boundary in case that the tag is missing */
+      pxt->tag[ia] |= MG_BDY;
+
       /* No swap of geometric edge */
       if ( MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) ) {
         continue;
@@ -331,6 +334,8 @@ int MMG3D_optbdry(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree,MMG5_in
       if ( ret < 0 )  return -1;
       /* CAUTION: trigger collapse with 2 elements */
       if ( ilist <= 1 )  continue;
+
+      /* Here, we work on a boundary edge lying along a boundary face */
       ier = MMG5_chkswpbdy(mesh,met,list,ilist,it1,it2,2);
       if ( ier <  0 )
         return -1;

src/mmg3d/swap_3d.c

@@ -77,12 +77,15 @@ int MMG5_chkswpbdy(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list,int ilist,
   np = pt->v[MMG5_iare[ia][0]];
   nq = pt->v[MMG5_iare[ia][1]];
 
+  // Algiane 05/04/24: I think that the assumption that was previously made that
+  // we can arrive from a tetrahedra without a boundary face (i.e. without an
+  // xtetra) never happens
+  assert ( pt->xt && "Boundary edges have to be swapped from a boundary face" );
+
   /* No swap of geometric edge */
-  if ( pt->xt ) {
-    pxt = &mesh->xtetra[pt->xt];
-    if ( (pxt->edg[ia]>0) || MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) ) {
-      return 0;
-    }
+  pxt = &mesh->xtetra[pt->xt];
+  if ( (pxt->edg[ia]>0) || MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) ) {
+    return 0;
   }
 
   /* No swap when either internal or external component has only 1 element (as
@@ -343,7 +346,26 @@ int MMG5_chkswpbdy(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list,int ilist,
   ppt0->c[1] = 0.5*(p0->c[1] + p1->c[1]);
   ppt0->c[2] = 0.5*(p0->c[2] + p1->c[2]);
 
+#ifndef NDEBUG
+  /* Security check: ensure that the edge is boundary */
+  int16_t  tag = 0;
+  MMG5_int ref = 0;
+  if ( !MMG3D_get_shellEdgeTag(mesh,list[0]/6,list[0]%6,&tag,&ref) ) {
+    fprintf(stderr,"\n  ## Warning: %s: 0. unable to get edge info"
+            " (tetra %d).\n",__func__,MMG3D_indElt(mesh,list[0]/6));
+    return 0;
+  }
+  assert ( (tag & MG_BDY)  && "Edge should be boundary but is not");
+#endif
+
   if ( met->m ) {
+    pt  = &mesh->tetra[list[0]/6];
+    assert ( pt->xt && "Boundary edge interpolated from non-boundary face");
+
+    /* Mark edge as boundary to ensure suitable detection of bdy edge during
+     * interpolation */
+    mesh->xtetra[pt->xt].tag[list[0]%6] |= MG_BDY;
+
     if ( typchk == 1 && (met->size>1) ) {
       if ( MMG3D_intmet33_ani(mesh,met,list[0]/6,list[0]%6,0,0.5) <= 0 )
         return 0;