From 694439c49ae2071f34545f9d4bf12a2b226ecb73 Mon Sep 17 00:00:00 2001 From: Corentin Prigent Date: Mon, 15 Apr 2024 16:45:40 +0200 Subject: [PATCH 1/5] Creation of skeleton of refactored MMG5_chkBdryTria; for draft PR purposes --- src/mmg3d/hash_3d.c | 313 ++++++++++++++++++++--------------- src/mmg3d/libmmg3d_private.h | 1 + 2 files changed, 182 insertions(+), 132 deletions(-) diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c index 0f15bbf7a..23a817021 100644 --- a/src/mmg3d/hash_3d.c +++ b/src/mmg3d/hash_3d.c @@ -1550,11 +1550,49 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { MMG5_pTria ptt,pttnew; MMG5_int *adja,adj,k,kk,i,j,ntmesh; MMG5_int ia,ib,ic, nbl,nt,ntpres; - int iface; + int iface, ier; MMG5_Hash hashElt, hashTri; /** Step 1: scan the mesh and count the boundaries */ - ntmesh = ntpres = 0; + ier = MMG5_chkBdryTria_boundaryCount(mesh,&ntmesh,&ntpres); + + /** Step 2: detect the extra boundaries (that will be ignored) provided by the + * user */ + if ( mesh->nt ) { + ier = MMG5_chkBdryTria_extraBoundary(mesh,ntmesh,&hashElt); + // Travel through the tria, delete those that are not in the hash tab or + // that are stored more that once. + ier = MMG5_chkbdryTria_deleteExtraBoundaries(); + } + ntmesh +=ntpres; + + /** Step 3: add the missing boundary triangles or, if the mesh contains + * prisms, set to required the triangles at interface betwen prisms and tet */ + ier = MMG5_chkBdryTria_step3(); + + /* Fill missing bdy triangles */ + return MMG5_bdryTria(mesh,ntmesh); +} + +/** + * \param mesh pointer to the mesh structure. + * \param ntmesh number of boundary triangles in the mesh. + * \param ntpres number of preserved boundaries in the mesh. + * \return 0 if failed, 1 if success. + * + * Step 1 of MMG5_chkBdryTria : scan the mesh and count the boundaries + * + */ +int MMG5_chkBdryTria_boundaryCount(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres) { + + MMG5_pTetra pt, pt1; + MMG5_pPrism pp, pp1; + MMG5_int *adja,adj,k; + MMG5_int ia,ib,ic,j; + MMG5_Hash hashTri; + int i; + + *ntmesh = *ntpres = 0; for (k=1; k<=mesh->ne; k++) { pt = &mesh->tetra[k]; if ( !MG_EOK(pt) ) continue; @@ -1563,13 +1601,13 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { adj = adja[i]; if ( !adj ) { - ++ntmesh; + ++(*ntmesh); continue; } adj /= 4; pt1 = &mesh->tetra[adj]; if ( pt->ref > pt1->ref ) - ++ntmesh; + ++(*ntmesh); } } @@ -1591,7 +1629,7 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { if ( pt->ref != pt1->ref ) continue; if ( (mesh->xtetra[pt->xt].ftag[i] & MG_BDY) && - (MG_GET(mesh->xtetra[pt->xt].ori,i) ) ) ++ntpres; + (MG_GET(mesh->xtetra[pt->xt].ori,i) ) ) ++(*ntpres); } } } @@ -1606,7 +1644,7 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { adj = adja[i]; if ( !adj ) { - ++ntmesh; + ++(*ntmesh); continue; } else if ( adj<0 ) { @@ -1616,14 +1654,14 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { adj /= 5; pp1 = &mesh->prism[adj]; if ( pp->ref > pp1->ref) { - ++ntmesh; + ++(*ntmesh); } } } /* Detect the triangles at the interface of the prisms and tetra (they have been * counted twice) */ - if ( ! MMG5_hashNew(mesh,&hashTri,0.51*ntmesh,1.51*ntmesh) ) return 0; + if ( ! MMG5_hashNew(mesh,&hashTri,0.51*(*ntmesh),1.51*(*ntmesh)) ) return 0; for (k=1; k<=mesh->ne; k++) { pt = &mesh->tetra[k]; if ( !MG_EOK(pt) ) continue; @@ -1661,151 +1699,101 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { j = MMG5_hashGetFace(&hashTri,ia,ib,ic); if ( !j ) continue; - --ntmesh; + --(*ntmesh); adja[i] = -j; } } MMG5_DEL_MEM(mesh,hashTri.item); } + return 1; +} - /** Step 2: detect the extra boundaries (that will be ignored) provided by the - * user */ - if ( mesh->nt ) { - if ( ! MMG5_hashNew(mesh,&hashElt,0.51*ntmesh,1.51*ntmesh) ) return 0; - // Hash the boundaries found in the mesh - if ( mesh->info.opnbdy) { - /* We want to keep the internal triangles: we must hash all the tetra faces */ - for (k=1; k<=mesh->ne; k++) { - pt = &mesh->tetra[k]; - if ( !MG_EOK(pt) ) continue; - - for (i=0; i<4; i++) { - ia = pt->v[MMG5_idir[i][0]]; - ib = pt->v[MMG5_idir[i][1]]; - ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; - } - } - } else { - for (k=1; k<=mesh->ne; k++) { - pt = &mesh->tetra[k]; - if ( !MG_EOK(pt) ) continue; - adja = &mesh->adja[4*(k-1)+1]; - for (i=0; i<4; i++) { - adj = adja[i]; - if ( !adj ) { - ia = pt->v[MMG5_idir[i][0]]; - ib = pt->v[MMG5_idir[i][1]]; - ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; - } - adj /= 4; +/** + * \param mesh pointer to the mesh structure. + * \param ntmesh number of boundary triangles in the mesh. + * \return 0 if failed, 1 if success. + * + * Step 2 of MMG5_chkBdryTria : detect extra boundaries + * + */ +int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt) { + + MMG5_pTetra pt, pt1; + MMG5_pPrism pp, pp1; + MMG5_int *adja,adj,k; + MMG5_int ia,ib,ic; + int i; + + if ( ! MMG5_hashNew(mesh,&hashElt,0.51*ntmesh,1.51*ntmesh) ) return 0; + // Hash the boundaries found in the mesh + if ( mesh->info.opnbdy) { + /* We want to keep the internal triangles: we must hash all the tetra faces */ + for (k=1; k<=mesh->ne; k++) { + pt = &mesh->tetra[k]; + if ( !MG_EOK(pt) ) continue; - pt1 = &mesh->tetra[adj]; - if ( pt->ref > pt1->ref ) { - ia = pt->v[MMG5_idir[i][0]]; - ib = pt->v[MMG5_idir[i][1]]; - ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; - } - } + for (i=0; i<4; i++) { + ia = pt->v[MMG5_idir[i][0]]; + ib = pt->v[MMG5_idir[i][1]]; + ic = pt->v[MMG5_idir[i][2]]; + if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; } } - for (k=1; k<=mesh->nprism; k++) { - pp = &mesh->prism[k]; - if ( !MG_EOK(pp) ) continue; - adja = &mesh->adjapr[5*(k-1)+1]; - for (i=0; i<2; i++) { + } else { + for (k=1; k<=mesh->ne; k++) { + pt = &mesh->tetra[k]; + if ( !MG_EOK(pt) ) continue; + adja = &mesh->adja[4*(k-1)+1]; + for (i=0; i<4; i++) { adj = adja[i]; if ( !adj ) { - ia = pp->v[MMG5_idir_pr[i][0]]; - ib = pp->v[MMG5_idir_pr[i][1]]; - ic = pp->v[MMG5_idir_pr[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; + ia = pt->v[MMG5_idir[i][0]]; + ib = pt->v[MMG5_idir[i][1]]; + ic = pt->v[MMG5_idir[i][2]]; + if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; } - else if ( adj<0 ) continue; - - adj /= 5; + adj /= 4; - pp1 = &mesh->prism[MMG5_abs(adj)]; - if ( pp->ref > pp1->ref ) { - ia = pp->v[MMG5_idir_pr[i][0]]; - ib = pp->v[MMG5_idir_pr[i][1]]; - ic = pp->v[MMG5_idir_pr[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; + pt1 = &mesh->tetra[adj]; + if ( pt->ref > pt1->ref ) { + ia = pt->v[MMG5_idir[i][0]]; + ib = pt->v[MMG5_idir[i][1]]; + ic = pt->v[MMG5_idir[i][2]]; + if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; } } } - - - // Travel through the tria, delete those that are not in the hash tab or - // that are stored more that once. - nt=0; nbl=1; - - if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0; - - for (k=1; k<=mesh->nt; k++) { - ptt = &mesh->tria[k]; - - ia = ptt->v[0]; - ib = ptt->v[1]; - ic = ptt->v[2]; - - i = MMG5_hashGetFace(&hashElt,ia,ib,ic); - j = MMG5_hashFace(mesh,&hashTri,ia,ib,ic,k); - - ptt->cc = i; - - if ( !j ) { - MMG5_DEL_MEM(mesh,hashElt.item); - MMG5_DEL_MEM(mesh,hashTri.item); - return 0; - } - else if ( j > 0 ) { - /* the face already exists in the tria table */ - continue; - } - - if ( !i ) { - /* the triangle is not a boundary tri or a tri at the interface of two - * subdomains with different references and the user don't ask to keep - * it. */ - continue; + } + for (k=1; k<=mesh->nprism; k++) { + pp = &mesh->prism[k]; + if ( !MG_EOK(pp) ) continue; + adja = &mesh->adjapr[5*(k-1)+1]; + for (i=0; i<2; i++) { + adj = adja[i]; + if ( !adj ) { + ia = pp->v[MMG5_idir_pr[i][0]]; + ib = pp->v[MMG5_idir_pr[i][1]]; + ic = pp->v[MMG5_idir_pr[i][2]]; + if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; } + else if ( adj<0 ) continue; - if ( mesh->info.opnbdy ) { - kk = i/4; - iface = i%4; - adj = mesh->adja[4*(kk-1)+1+iface]; - /* Check if we have found a triangle at the interface of 2 doms of same - * ref */ - if ( adj && mesh->tetra[kk].ref == mesh->tetra[adj/4].ref ) { - ++ntpres; - } - } + adj /= 5; - ++nt; - if ( k!=nbl ) { - pttnew = &mesh->tria[nbl]; - memcpy(pttnew,ptt,sizeof(MMG5_Tria)); + pp1 = &mesh->prism[MMG5_abs(adj)]; + if ( pp->ref > pp1->ref ) { + ia = pp->v[MMG5_idir_pr[i][0]]; + ib = pp->v[MMG5_idir_pr[i][1]]; + ic = pp->v[MMG5_idir_pr[i][2]]; + if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; } - ++nbl; } - nbl = mesh->nt-nt; - if ( nbl ) { - fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries provided." - " Ignored\n",__func__,nbl); - MMG5_ADD_MEM(mesh,(-nbl)*sizeof(MMG5_Tria),"triangles",return 0); - MMG5_SAFE_REALLOC(mesh->tria,mesh->nt+1,nt+1,MMG5_Tria,"triangles",return 0); - mesh->nt = nt; - } - MMG5_DEL_MEM(mesh,hashElt.item); - MMG5_DEL_MEM(mesh,hashTri.item); } - ntmesh +=ntpres; + return 1; +} + +int MMG5_chkBdryTria_step3() { - /** Step 3: add the missing boundary triangles or, if the mesh contains - * prisms, set to required the triangles at interface betwen prisms and tet */ if ( ntpres && (mesh->info.imprim > 5 || mesh->info.ddebug) ) printf(" %" MMG5_PRId " triangles between 2 tetrahdra with same" " references\n",ntpres); @@ -1837,10 +1825,71 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries founded\n", __func__,nbl); - /* Fill missing bdy triangles */ - return MMG5_bdryTria(mesh,ntmesh); } +int MMG5_chkbdryTria_deleteExtraBoundaries() { + nt=0; nbl=1; + + if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0; + + for (k=1; k<=mesh->nt; k++) { + ptt = &mesh->tria[k]; + + ia = ptt->v[0]; + ib = ptt->v[1]; + ic = ptt->v[2]; + + i = MMG5_hashGetFace(&hashElt,ia,ib,ic); + j = MMG5_hashFace(mesh,&hashTri,ia,ib,ic,k); + + ptt->cc = i; + + if ( !j ) { + MMG5_DEL_MEM(mesh,hashElt.item); + MMG5_DEL_MEM(mesh,hashTri.item); + return 0; + } + else if ( j > 0 ) { + /* the face already exists in the tria table */ + continue; + } + + if ( !i ) { + /* the triangle is not a boundary tri or a tri at the interface of two + * subdomains with different references and the user don't ask to keep + * it. */ + continue; + } + + if ( mesh->info.opnbdy ) { + kk = i/4; + iface = i%4; + adj = mesh->adja[4*(kk-1)+1+iface]; + /* Check if we have found a triangle at the interface of 2 doms of same + * ref */ + if ( adj && mesh->tetra[kk].ref == mesh->tetra[adj/4].ref ) { + ++ntpres; + } + } + + ++nt; + if ( k!=nbl ) { + pttnew = &mesh->tria[nbl]; + memcpy(pttnew,ptt,sizeof(MMG5_Tria)); + } + ++nbl; + } + nbl = mesh->nt-nt; + if ( nbl ) { + fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries provided." + " Ignored\n",__func__,nbl); + MMG5_ADD_MEM(mesh,(-nbl)*sizeof(MMG5_Tria),"triangles",return 0); + MMG5_SAFE_REALLOC(mesh->tria,mesh->nt+1,nt+1,MMG5_Tria,"triangles",return 0); + mesh->nt = nt; + } + MMG5_DEL_MEM(mesh,hashElt.item); + MMG5_DEL_MEM(mesh,hashTri.item); +} /** * \param mesh pointer to the mesh structure. diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h index cf80a2e78..3857fa24d 100644 --- a/src/mmg3d/libmmg3d_private.h +++ b/src/mmg3d/libmmg3d_private.h @@ -233,6 +233,7 @@ void MMG5_freeXTets(MMG5_pMesh mesh); void MMG5_freeXPrisms(MMG5_pMesh mesh); void MMG3D_Free_topoTables(MMG5_pMesh mesh); int MMG5_chkBdryTria(MMG5_pMesh mesh); +int MMG5_chkBdryTria_boundaryCount(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres); int MMG5_mmg3dBezierCP(MMG5_pMesh mesh,MMG5_Tria *pt,MMG5_pBezier pb,int8_t ori); extern int MMG5_BezierTgt(double c1[3],double c2[3],double n1[3],double n2[3],double t1[3],double t2[3]); extern double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3]); From c0d254d6f729f64e4bb1f1ae162246620afd8a11 Mon Sep 17 00:00:00 2001 From: Corentin Prigent Date: Tue, 16 Apr 2024 14:49:34 +0200 Subject: [PATCH 2/5] Proper factorization of function MMG5_chkBdryTria --- src/mmg3d/hash_3d.c | 146 +++++++++++++++++++++-------------- src/mmg3d/libmmg3d_private.h | 5 +- 2 files changed, 90 insertions(+), 61 deletions(-) diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c index 23a817021..5c4f1ba4c 100644 --- a/src/mmg3d/hash_3d.c +++ b/src/mmg3d/hash_3d.c @@ -1545,33 +1545,28 @@ int MMG5_bdryTria(MMG5_pMesh mesh, MMG5_int ntmesh) { * */ int MMG5_chkBdryTria(MMG5_pMesh mesh) { - MMG5_pTetra pt,pt1; - MMG5_pPrism pp,pp1; - MMG5_pTria ptt,pttnew; - MMG5_int *adja,adj,k,kk,i,j,ntmesh; - MMG5_int ia,ib,ic, nbl,nt,ntpres; - int iface, ier; - MMG5_Hash hashElt, hashTri; + MMG5_int ntmesh,ntpres; + int ier; + MMG5_Hash hashElt; /** Step 1: scan the mesh and count the boundaries */ - ier = MMG5_chkBdryTria_boundaryCount(mesh,&ntmesh,&ntpres); + ier = MMG5_chkBdryTria_countBoundaries(mesh,&ntmesh,&ntpres); /** Step 2: detect the extra boundaries (that will be ignored) provided by the * user */ if ( mesh->nt ) { - ier = MMG5_chkBdryTria_extraBoundary(mesh,ntmesh,&hashElt); + ier = MMG5_chkBdryTria_hashBoundaries(mesh,ntmesh,&hashElt); // Travel through the tria, delete those that are not in the hash tab or // that are stored more that once. - ier = MMG5_chkbdryTria_deleteExtraBoundaries(); + ier = MMG5_chkBdryTria_deleteExtraBoundaries(mesh,&ntpres,&hashElt); } ntmesh +=ntpres; /** Step 3: add the missing boundary triangles or, if the mesh contains * prisms, set to required the triangles at interface betwen prisms and tet */ - ier = MMG5_chkBdryTria_step3(); + ier = MMG5_chkBdryTria_addMissingTriangles(mesh,ntmesh,ntpres); - /* Fill missing bdy triangles */ - return MMG5_bdryTria(mesh,ntmesh); + return 1; } /** @@ -1583,7 +1578,7 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { * Step 1 of MMG5_chkBdryTria : scan the mesh and count the boundaries * */ -int MMG5_chkBdryTria_boundaryCount(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres) { +int MMG5_chkBdryTria_countBoundaries(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres) { MMG5_pTetra pt, pt1; MMG5_pPrism pp, pp1; @@ -1711,12 +1706,13 @@ int MMG5_chkBdryTria_boundaryCount(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int * /** * \param mesh pointer to the mesh structure. * \param ntmesh number of boundary triangles in the mesh. + * \param hashElt pointer towards face hash table. * \return 0 if failed, 1 if success. * - * Step 2 of MMG5_chkBdryTria : detect extra boundaries + * Step 2 of MMG5_chkBdryTria : create hash tables of boundaries in the mesh * */ -int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt) { +int MMG5_chkBdryTria_hashBoundaries(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt) { MMG5_pTetra pt, pt1; MMG5_pPrism pp, pp1; @@ -1724,7 +1720,7 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * MMG5_int ia,ib,ic; int i; - if ( ! MMG5_hashNew(mesh,&hashElt,0.51*ntmesh,1.51*ntmesh) ) return 0; + if ( ! MMG5_hashNew(mesh,hashElt,0.51*ntmesh,1.51*ntmesh) ) return 0; // Hash the boundaries found in the mesh if ( mesh->info.opnbdy) { /* We want to keep the internal triangles: we must hash all the tetra faces */ @@ -1736,7 +1732,7 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * ia = pt->v[MMG5_idir[i][0]]; ib = pt->v[MMG5_idir[i][1]]; ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; + if ( !MMG5_hashFace(mesh,hashElt,ia,ib,ic,4*k+i) ) return 0; } } } else { @@ -1750,7 +1746,7 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * ia = pt->v[MMG5_idir[i][0]]; ib = pt->v[MMG5_idir[i][1]]; ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; + if ( !MMG5_hashFace(mesh,hashElt,ia,ib,ic,4*k+i) ) return 0; } adj /= 4; @@ -1759,7 +1755,7 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * ia = pt->v[MMG5_idir[i][0]]; ib = pt->v[MMG5_idir[i][1]]; ic = pt->v[MMG5_idir[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0; + if ( !MMG5_hashFace(mesh,hashElt,ia,ib,ic,4*k+i) ) return 0; } } } @@ -1774,7 +1770,7 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * ia = pp->v[MMG5_idir_pr[i][0]]; ib = pp->v[MMG5_idir_pr[i][1]]; ic = pp->v[MMG5_idir_pr[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; + if ( !MMG5_hashFace(mesh,hashElt,ia,ib,ic,5*k+i) ) return 0; } else if ( adj<0 ) continue; @@ -1785,49 +1781,29 @@ int MMG5_chkBdryTria_extraBoundary(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * ia = pp->v[MMG5_idir_pr[i][0]]; ib = pp->v[MMG5_idir_pr[i][1]]; ic = pp->v[MMG5_idir_pr[i][2]]; - if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0; + if ( !MMG5_hashFace(mesh,hashElt,ia,ib,ic,5*k+i) ) return 0; } } } return 1; } -int MMG5_chkBdryTria_step3() { - - if ( ntpres && (mesh->info.imprim > 5 || mesh->info.ddebug) ) - printf(" %" MMG5_PRId " triangles between 2 tetrahdra with same" - " references\n",ntpres); - - if ( mesh->nt==ntmesh && !mesh->nprism ) { - for (k=1; k<=mesh->nt; k++) { - ptt = &mesh->tria[k]; - for (i=0; i<3; i++) { - mesh->point[ptt->v[i]].tag |= MG_BDY; - } - } - return 1; - } - - nbl = 0; - if ( !mesh->nt ) { - MMG5_ADD_MEM(mesh,(ntmesh+1)*sizeof(MMG5_Tria),"triangles",return 0); - MMG5_SAFE_CALLOC(mesh->tria,ntmesh+1,MMG5_Tria,return 0); - } - else { - assert((!mesh->nprism && ntmesh>mesh->nt)||(mesh->nprism && ntmesh>=mesh->nt)); - if ( ntmesh > mesh->nt ) { - MMG5_ADD_MEM(mesh,(ntmesh-mesh->nt)*sizeof(MMG5_Tria),"triangles",return 0); - MMG5_SAFE_RECALLOC(mesh->tria,mesh->nt+1,ntmesh+1,MMG5_Tria,"triangles",return 0); - nbl = ntmesh-mesh->nt; - } - } - if ( nbl && (mesh->info.imprim > 5 || mesh->info.ddebug) ) - fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries founded\n", - __func__,nbl); +/** + * \param mesh pointer to the mesh structure. + * \return 0 if failed, 1 if success. + * + * Step 2.5 of MMG5_chkBdryTria : Travel through the tria, delete those that are not in the hash tab or + * that are stored more that once. + * + */ +int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash *hashElt) { -} + MMG5_pTria ptt, pttnew; + MMG5_Hash hashTri; + MMG5_int nt, nbl, k, kk, i, j; + MMG5_int ia, ib, ic, adj; + int iface; -int MMG5_chkbdryTria_deleteExtraBoundaries() { nt=0; nbl=1; if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0; @@ -1839,13 +1815,13 @@ int MMG5_chkbdryTria_deleteExtraBoundaries() { ib = ptt->v[1]; ic = ptt->v[2]; - i = MMG5_hashGetFace(&hashElt,ia,ib,ic); + i = MMG5_hashGetFace(hashElt,ia,ib,ic); j = MMG5_hashFace(mesh,&hashTri,ia,ib,ic,k); ptt->cc = i; if ( !j ) { - MMG5_DEL_MEM(mesh,hashElt.item); + MMG5_DEL_MEM(mesh,hashElt->item); MMG5_DEL_MEM(mesh,hashTri.item); return 0; } @@ -1868,7 +1844,7 @@ int MMG5_chkbdryTria_deleteExtraBoundaries() { /* Check if we have found a triangle at the interface of 2 doms of same * ref */ if ( adj && mesh->tetra[kk].ref == mesh->tetra[adj/4].ref ) { - ++ntpres; + ++(*ntpres); } } @@ -1887,8 +1863,58 @@ int MMG5_chkbdryTria_deleteExtraBoundaries() { MMG5_SAFE_REALLOC(mesh->tria,mesh->nt+1,nt+1,MMG5_Tria,"triangles",return 0); mesh->nt = nt; } - MMG5_DEL_MEM(mesh,hashElt.item); + MMG5_DEL_MEM(mesh,hashElt->item); MMG5_DEL_MEM(mesh,hashTri.item); + return 1; +} + +/** + * \param mesh pointer to the mesh structure. + * \return 0 if failed, 1 if success. + * + * Step 3 of MMG5_chkBdryTria : add the missing boundary triangles or, if the mesh contains + * prisms, set to required the triangles at interface betwen prisms and tet + * + */ +int MMG5_chkBdryTria_addMissingTriangles(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_int ntpres) { + + MMG5_pTria ptt; + MMG5_int k, nbl; + int i; + + if ( ntpres && (mesh->info.imprim > 5 || mesh->info.ddebug) ) + printf(" %" MMG5_PRId " triangles between 2 tetrahdra with same" + " references\n",ntpres); + + if ( mesh->nt==ntmesh && !mesh->nprism ) { + for (k=1; k<=mesh->nt; k++) { + ptt = &mesh->tria[k]; + for (i=0; i<3; i++) { + mesh->point[ptt->v[i]].tag |= MG_BDY; + } + } + return 1; + } + + nbl = 0; + if ( !mesh->nt ) { + MMG5_ADD_MEM(mesh,(ntmesh+1)*sizeof(MMG5_Tria),"triangles",return 0); + MMG5_SAFE_CALLOC(mesh->tria,ntmesh+1,MMG5_Tria,return 0); + } + else { + assert((!mesh->nprism && ntmesh>mesh->nt)||(mesh->nprism && ntmesh>=mesh->nt)); + if ( ntmesh > mesh->nt ) { + MMG5_ADD_MEM(mesh,(ntmesh-mesh->nt)*sizeof(MMG5_Tria),"triangles",return 0); + MMG5_SAFE_RECALLOC(mesh->tria,mesh->nt+1,ntmesh+1,MMG5_Tria,"triangles",return 0); + nbl = ntmesh-mesh->nt; + } + } + if ( nbl && (mesh->info.imprim > 5 || mesh->info.ddebug) ) + fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries founded\n", + __func__,nbl); + + /* Fill missing bdy triangles */ + return MMG5_bdryTria(mesh,ntmesh); } /** diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h index 3857fa24d..9eb4ffd22 100644 --- a/src/mmg3d/libmmg3d_private.h +++ b/src/mmg3d/libmmg3d_private.h @@ -233,7 +233,10 @@ void MMG5_freeXTets(MMG5_pMesh mesh); void MMG5_freeXPrisms(MMG5_pMesh mesh); void MMG3D_Free_topoTables(MMG5_pMesh mesh); int MMG5_chkBdryTria(MMG5_pMesh mesh); -int MMG5_chkBdryTria_boundaryCount(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres); +int MMG5_chkBdryTria_countBoundaries(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres); +int MMG5_chkBdryTria_hashBoundaries(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt); +int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash* hashElt); +int MMG5_chkBdryTria_addMissingTriangles(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_int ntpres); int MMG5_mmg3dBezierCP(MMG5_pMesh mesh,MMG5_Tria *pt,MMG5_pBezier pb,int8_t ori); extern int MMG5_BezierTgt(double c1[3],double c2[3],double n1[3],double n2[3],double t1[3],double t2[3]); extern double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3]); From 9c45b08ac87e3ec768afe834ac1b673ae6b226eb Mon Sep 17 00:00:00 2001 From: Corentin Prigent Date: Tue, 16 Apr 2024 15:47:44 +0200 Subject: [PATCH 3/5] Extracted triangle deletion phase outside of triangle flagging loop --- src/mmg3d/hash_3d.c | 32 ++++++++++++++++++++++++-------- src/mmg3d/libmmg3d_private.h | 3 ++- 2 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c index 5c4f1ba4c..a55554c78 100644 --- a/src/mmg3d/hash_3d.c +++ b/src/mmg3d/hash_3d.c @@ -1556,9 +1556,11 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { * user */ if ( mesh->nt ) { ier = MMG5_chkBdryTria_hashBoundaries(mesh,ntmesh,&hashElt); - // Travel through the tria, delete those that are not in the hash tab or + // Travel through the tria, flag those that are not in the hash tab or // that are stored more that once. - ier = MMG5_chkBdryTria_deleteExtraBoundaries(mesh,&ntpres,&hashElt); + ier = MMG5_chkBdryTria_flagExtraTriangles(mesh,&ntpres,&hashElt); + // Delete flagged triangles + ier = MMG5_chkBdryTria_deleteExtraTriangles(mesh); } ntmesh +=ntpres; @@ -1796,16 +1798,14 @@ int MMG5_chkBdryTria_hashBoundaries(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash * that are stored more that once. * */ -int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash *hashElt) { +int MMG5_chkBdryTria_flagExtraTriangles(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash *hashElt) { MMG5_pTria ptt, pttnew; MMG5_Hash hashTri; - MMG5_int nt, nbl, k, kk, i, j; + MMG5_int k, kk, i, j; MMG5_int ia, ib, ic, adj; int iface; - nt=0; nbl=1; - if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0; for (k=1; k<=mesh->nt; k++) { @@ -1827,6 +1827,7 @@ int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MM } else if ( j > 0 ) { /* the face already exists in the tria table */ + ptt->v[0] = 0; continue; } @@ -1834,6 +1835,7 @@ int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MM /* the triangle is not a boundary tri or a tri at the interface of two * subdomains with different references and the user don't ask to keep * it. */ + ptt->v[0] = 0; continue; } @@ -1847,6 +1849,22 @@ int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MM ++(*ntpres); } } + } + MMG5_DEL_MEM(mesh,hashElt->item); + MMG5_DEL_MEM(mesh,hashTri.item); + return 1; +} + +int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh) { + + MMG5_pTria ptt, pttnew; + MMG5_int nt, nbl, k; + + nt = 0; nbl = 1; + for (k=1; k<=mesh->nt; k++) { + ptt = &mesh->tria[k]; + + if ( !MG_EOK(ptt) ) continue; ++nt; if ( k!=nbl ) { @@ -1863,8 +1881,6 @@ int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MM MMG5_SAFE_REALLOC(mesh->tria,mesh->nt+1,nt+1,MMG5_Tria,"triangles",return 0); mesh->nt = nt; } - MMG5_DEL_MEM(mesh,hashElt->item); - MMG5_DEL_MEM(mesh,hashTri.item); return 1; } diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h index 9eb4ffd22..949d4252b 100644 --- a/src/mmg3d/libmmg3d_private.h +++ b/src/mmg3d/libmmg3d_private.h @@ -235,8 +235,9 @@ void MMG3D_Free_topoTables(MMG5_pMesh mesh); int MMG5_chkBdryTria(MMG5_pMesh mesh); int MMG5_chkBdryTria_countBoundaries(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_int *ntpres); int MMG5_chkBdryTria_hashBoundaries(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt); -int MMG5_chkBdryTria_deleteExtraBoundaries(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash* hashElt); +int MMG5_chkBdryTria_flagExtraTriangles(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash* hashElt); int MMG5_chkBdryTria_addMissingTriangles(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_int ntpres); +int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh); int MMG5_mmg3dBezierCP(MMG5_pMesh mesh,MMG5_Tria *pt,MMG5_pBezier pb,int8_t ori); extern int MMG5_BezierTgt(double c1[3],double c2[3],double n1[3],double n2[3],double t1[3],double t2[3]); extern double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3]); From a82f2c328cb0f0200cf47eb6a0ad9925d342c51c Mon Sep 17 00:00:00 2001 From: Corentin Prigent Date: Tue, 16 Apr 2024 16:06:52 +0200 Subject: [PATCH 4/5] added permtria array to mmg --- src/mmg3d/hash_3d.c | 7 +++++-- src/mmg3d/libmmg3d_private.h | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c index a55554c78..50fa9080f 100644 --- a/src/mmg3d/hash_3d.c +++ b/src/mmg3d/hash_3d.c @@ -1560,7 +1560,7 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { // that are stored more that once. ier = MMG5_chkBdryTria_flagExtraTriangles(mesh,&ntpres,&hashElt); // Delete flagged triangles - ier = MMG5_chkBdryTria_deleteExtraTriangles(mesh); + ier = MMG5_chkBdryTria_deleteExtraTriangles(mesh,NULL); } ntmesh +=ntpres; @@ -1855,7 +1855,7 @@ int MMG5_chkBdryTria_flagExtraTriangles(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_ return 1; } -int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh) { +int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh, MMG5_int* permtria) { MMG5_pTria ptt, pttnew; MMG5_int nt, nbl, k; @@ -1869,6 +1869,9 @@ int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh) { ++nt; if ( k!=nbl ) { pttnew = &mesh->tria[nbl]; + if ( permtria ) { + permtria[k] = nbl; + } memcpy(pttnew,ptt,sizeof(MMG5_Tria)); } ++nbl; diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h index 949d4252b..6b841411f 100644 --- a/src/mmg3d/libmmg3d_private.h +++ b/src/mmg3d/libmmg3d_private.h @@ -237,7 +237,7 @@ int MMG5_chkBdryTria_countBoundaries(MMG5_pMesh mesh, MMG5_int *ntmesh, MMG5_in int MMG5_chkBdryTria_hashBoundaries(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_Hash *hashElt); int MMG5_chkBdryTria_flagExtraTriangles(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash* hashElt); int MMG5_chkBdryTria_addMissingTriangles(MMG5_pMesh mesh, MMG5_int ntmesh, MMG5_int ntpres); -int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh); +int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh, MMG5_int* permtria); int MMG5_mmg3dBezierCP(MMG5_pMesh mesh,MMG5_Tria *pt,MMG5_pBezier pb,int8_t ori); extern int MMG5_BezierTgt(double c1[3],double c2[3],double n1[3],double n2[3],double t1[3],double t2[3]); extern double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3]); From cf436e4297247e7d09502940be5302447f437511 Mon Sep 17 00:00:00 2001 From: Corentin Prigent Date: Thu, 18 Apr 2024 16:37:49 +0200 Subject: [PATCH 5/5] Added test coverage of duplicate triangles --- cmake/testing/mmg3d_tests.cmake | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cmake/testing/mmg3d_tests.cmake b/cmake/testing/mmg3d_tests.cmake index fa5baf52a..0c64adeff 100644 --- a/cmake/testing/mmg3d_tests.cmake +++ b/cmake/testing/mmg3d_tests.cmake @@ -539,6 +539,11 @@ ADD_TEST(NAME mmg3d_opnbdy_ref_island -in ${MMG3D_CI_TESTS}/OpnBdy_island/island -out ${CTEST_OUTPUT_DIR}/mmg3d_OpnBdy_island.o.meshb) +ADD_TEST(NAME mmg3d_duplicate_triangle +COMMAND ${EXECUT_MMG3D} -v 5 +-in ${MMG3D_CI_TESTS}/DuplicateTriangle/duplicate_triangle +-out ${CTEST_OUTPUT_DIR}/duplicate-triangle.o.meshb) + ############################################################################### ##### ##### Check Lagrangian motion option