From d2bf05b201787fd88eeb274fc542376f30ea8970 Mon Sep 17 00:00:00 2001
From: Algiane Froehly <algiane@gmail.com>
Date: Tue, 19 Dec 2023 17:19:05 +0100
Subject: [PATCH] Add MG_BDY tag setting before interpolation along boundary
 edges.

---
 src/mmg3d/mmg3d1.c         | 15 +++++++++++++--
 src/mmg3d/mmg3d1_delone.c  |  4 ++++
 src/mmg3d/mmg3d1_pattern.c |  2 +-
 src/mmg3d/swap_3d.c        | 19 +++++++++++++++++++
 4 files changed, 37 insertions(+), 3 deletions(-)

diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index be8afb1e3..8ca3b3b87 100644
--- a/src/mmg3d/mmg3d1.c
+++ b/src/mmg3d/mmg3d1.c
@@ -1936,7 +1936,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);
@@ -2348,6 +2355,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);
@@ -2418,7 +2429,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);
diff --git a/src/mmg3d/mmg3d1_delone.c b/src/mmg3d/mmg3d1_delone.c
index eb397a448..b3b0988cd 100644
--- a/src/mmg3d/mmg3d1_delone.c
+++ b/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);
     }
diff --git a/src/mmg3d/mmg3d1_pattern.c b/src/mmg3d/mmg3d1_pattern.c
index 46998fc24..56dcde424 100644
--- a/src/mmg3d/mmg3d1_pattern.c
+++ b/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);
diff --git a/src/mmg3d/swap_3d.c b/src/mmg3d/swap_3d.c
index 6bec02880..8e95dce0d 100644
--- a/src/mmg3d/swap_3d.c
+++ b/src/mmg3d/swap_3d.c
@@ -343,7 +343,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;