From ffe6fc25aa464c0f46042160669d43e85a0c652c Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Tue, 19 Dec 2023 15:32:38 +0100 Subject: [PATCH 1/6] Check aniso metric interpolation Add an assertion to ensure that no boundary edge is seen as internal during anisotropic metric interpolation. Note that some ci tests are failing because indeed, some boundary edges are not detected as boundary. --- src/mmg3d/intmet_3d.c | 41 +++++++++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/src/mmg3d/intmet_3d.c b/src/mmg3d/intmet_3d.c index 4db3c516e..f3f0376d2 100644 --- a/src/mmg3d/intmet_3d.c +++ b/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); + } /** From d2bf05b201787fd88eeb274fc542376f30ea8970 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Tue, 19 Dec 2023 17:19:05 +0100 Subject: [PATCH 2/6] 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; From 4d7e06ed787e1f800f73784b050a4f8e47a44fc1 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Tue, 19 Dec 2023 17:23:08 +0100 Subject: [PATCH 3/6] TO REVERT after merge of PR#234: Add check over bdy faces in get_shellEdgeTag. --- src/mmg3d/colver_3d.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mmg3d/colver_3d.c b/src/mmg3d/colver_3d.c index 5400b4a72..3c2b34dfa 100644 --- a/src/mmg3d/colver_3d.c +++ b/src/mmg3d/colver_3d.c @@ -393,8 +393,8 @@ int MMG3D_get_shellEdgeTag_oneDir(MMG5_pMesh mesh,MMG5_int start, MMG5_int na, if ( pt->xt ) { pxt = &mesh->xtetra[pt->xt]; *ref = pxt->edg[i]; - if ( pxt->tag[i] & MG_BDY ) { - *tag |= pxt->tag[i]; + if ((pxt->ftag[MMG5_ifar[i][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[i][1]] & MG_BDY)) { + *tag |= (pxt->tag[i] | MG_BDY); *filled = 1; return adj; } @@ -444,8 +444,8 @@ int MMG3D_get_shellEdgeTag(MMG5_pMesh mesh,MMG5_int start, int8_t ia,int16_t *t if ( pt->xt ) { pxt = &mesh->xtetra[pt->xt]; - if ( pxt->tag[ia] & MG_BDY ) { - *tag |= pxt->tag[ia]; + if ((pxt->ftag[MMG5_ifar[ia][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[ia][1]] & MG_BDY)) { + *tag |= (pxt->tag[ia] | MG_BDY); *ref = pxt->edg[ia]; return 1; } From b86b88f7264e013e5375620b93c1c9d8a3d1caca Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 5 Apr 2024 11:44:25 +0200 Subject: [PATCH 4/6] Revert "TO REVERT after merge of PR#234: Add check over bdy faces in get_shellEdgeTag." This reverts commit 4d7e06ed787e1f800f73784b050a4f8e47a44fc1. --- src/mmg3d/colver_3d.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mmg3d/colver_3d.c b/src/mmg3d/colver_3d.c index 3c2b34dfa..5400b4a72 100644 --- a/src/mmg3d/colver_3d.c +++ b/src/mmg3d/colver_3d.c @@ -393,8 +393,8 @@ int MMG3D_get_shellEdgeTag_oneDir(MMG5_pMesh mesh,MMG5_int start, MMG5_int na, if ( pt->xt ) { pxt = &mesh->xtetra[pt->xt]; *ref = pxt->edg[i]; - if ((pxt->ftag[MMG5_ifar[i][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[i][1]] & MG_BDY)) { - *tag |= (pxt->tag[i] | MG_BDY); + if ( pxt->tag[i] & MG_BDY ) { + *tag |= pxt->tag[i]; *filled = 1; return adj; } @@ -444,8 +444,8 @@ int MMG3D_get_shellEdgeTag(MMG5_pMesh mesh,MMG5_int start, int8_t ia,int16_t *t if ( pt->xt ) { pxt = &mesh->xtetra[pt->xt]; - if ((pxt->ftag[MMG5_ifar[ia][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[ia][1]] & MG_BDY)) { - *tag |= (pxt->tag[ia] | MG_BDY); + if ( pxt->tag[ia] & MG_BDY ) { + *tag |= pxt->tag[ia]; *ref = pxt->edg[ia]; return 1; } From 2c3b08691d6b64696e5b0a9784d678f048779337 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 5 Apr 2024 15:29:56 +0200 Subject: [PATCH 5/6] Add MG_BDY tag to the edge before the call to chkswpbdy in optbdry. --- src/mmg3d/mmg3d1.c | 4 ++++ src/mmg3d/optbdry_3d.c | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c index 67e70b59a..6d681d916 100644 --- a/src/mmg3d/mmg3d1.c +++ b/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; diff --git a/src/mmg3d/optbdry_3d.c b/src/mmg3d/optbdry_3d.c index 7d446eb60..916423062 100644 --- a/src/mmg3d/optbdry_3d.c +++ b/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; From cf23ad20afa89fb2bac27ddb611bc103f3dc01d4 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 5 Apr 2024 15:33:05 +0200 Subject: [PATCH 6/6] Remove the assumption that we can enter the chkswpbdy function from a non boundary tetra and add the associated check (a similar check was already added by commit ffe6fc25aa4). --- src/mmg3d/swap_3d.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/mmg3d/swap_3d.c b/src/mmg3d/swap_3d.c index 4888aff4e..3102816ae 100644 --- a/src/mmg3d/swap_3d.c +++ b/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