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 diff --git a/src/mmg3d/hash_3d.c b/src/mmg3d/hash_3d.c index 0f15bbf7a..50fa9080f 100644 --- a/src/mmg3d/hash_3d.c +++ b/src/mmg3d/hash_3d.c @@ -1545,16 +1545,51 @@ 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; - MMG5_Hash hashElt, hashTri; + MMG5_int ntmesh,ntpres; + int ier; + MMG5_Hash hashElt; /** Step 1: scan the mesh and count the boundaries */ - ntmesh = ntpres = 0; + 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_hashBoundaries(mesh,ntmesh,&hashElt); + // Travel through the tria, flag those that are not in the hash tab or + // that are stored more that once. + ier = MMG5_chkBdryTria_flagExtraTriangles(mesh,&ntpres,&hashElt); + // Delete flagged triangles + ier = MMG5_chkBdryTria_deleteExtraTriangles(mesh,NULL); + } + 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_addMissingTriangles(mesh,ntmesh,ntpres); + + return 1; +} + +/** + * \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_countBoundaries(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 +1598,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 +1626,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 +1641,7 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { adj = adja[i]; if ( !adj ) { - ++ntmesh; + ++(*ntmesh); continue; } else if ( adj<0 ) { @@ -1616,14 +1651,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 +1696,211 @@ 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++) { +/** + * \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 : create hash tables of boundaries in the mesh + * + */ +int MMG5_chkBdryTria_hashBoundaries(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; + + 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; + 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; + adj /= 4; - 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; - } + 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 (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; + } + 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; - adj /= 5; + adj /= 5; - 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; - } + 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; } } + } + return 1; +} +/** + * \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_flagExtraTriangles(MMG5_pMesh mesh, MMG5_int* ntpres, MMG5_Hash *hashElt) { - // Travel through the tria, delete those that are not in the hash tab or - // that are stored more that once. - nt=0; nbl=1; + MMG5_pTria ptt, pttnew; + MMG5_Hash hashTri; + MMG5_int k, kk, i, j; + MMG5_int ia, ib, ic, adj; + int iface; - if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0; + 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]; + for (k=1; k<=mesh->nt; k++) { + ptt = &mesh->tria[k]; - ia = ptt->v[0]; - ib = ptt->v[1]; - ic = ptt->v[2]; + 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); + i = MMG5_hashGetFace(hashElt,ia,ib,ic); + j = MMG5_hashFace(mesh,&hashTri,ia,ib,ic,k); - ptt->cc = i; + 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 ( !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 */ + ptt->v[0] = 0; + 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 ( !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. */ + ptt->v[0] = 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; - } + 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); } + } + } + MMG5_DEL_MEM(mesh,hashElt->item); + MMG5_DEL_MEM(mesh,hashTri.item); + return 1; +} + +int MMG5_chkBdryTria_deleteExtraTriangles(MMG5_pMesh mesh, MMG5_int* permtria) { + + 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 ) { - pttnew = &mesh->tria[nbl]; - memcpy(pttnew,ptt,sizeof(MMG5_Tria)); + ++nt; + if ( k!=nbl ) { + pttnew = &mesh->tria[nbl]; + if ( permtria ) { + permtria[k] = nbl; } - ++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; + memcpy(pttnew,ptt,sizeof(MMG5_Tria)); } - MMG5_DEL_MEM(mesh,hashElt.item); - MMG5_DEL_MEM(mesh,hashTri.item); + ++nbl; } - ntmesh +=ntpres; + 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; + } + 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; - /** 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); @@ -1841,7 +1936,6 @@ int MMG5_chkBdryTria(MMG5_pMesh mesh) { return MMG5_bdryTria(mesh,ntmesh); } - /** * \param mesh pointer to the mesh structure. * \return 0 if failed, 1 if success. diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h index cf80a2e78..6b841411f 100644 --- a/src/mmg3d/libmmg3d_private.h +++ b/src/mmg3d/libmmg3d_private.h @@ -233,6 +233,11 @@ 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_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_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, 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]);