From 5a73683f84fe422031921bef4ced8905d8b9eb7e Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 31 May 2024 16:49:42 +0200 Subject: [PATCH] Add comment to highlight the inconsistencies in edge lengths computations. --- src/common/inlined_functions_private.h | 25 ++++++++++-------- src/mmg3d/colver_3d.c | 4 +++ src/mmg3d/inlined_functions_3d_private.h | 15 +++++++++-- src/mmg3d/libmmg3d_tools.c | 32 ++++++++++++++++++++++++ src/mmg3d/quality_3d.c | 6 +++++ src/mmg3d/split_3d.c | 21 ++++++++++++++++ 6 files changed, 90 insertions(+), 13 deletions(-) diff --git a/src/common/inlined_functions_private.h b/src/common/inlined_functions_private.h index 451c6a4ee..16bb1a202 100644 --- a/src/common/inlined_functions_private.h +++ b/src/common/inlined_functions_private.h @@ -43,11 +43,12 @@ * \param m0 metric at point np0. * \param m1 metric at point np1. * \param isedg 1 if the edge is a ridge, 0 otherwise. - * \return length of edge according to the prescribed metric, 0 if fail. + * \return length of a curve edge according to the prescribed metric, 0 if fail. * - * Compute length of surface edge \f$[np0;np1]\f$ according to the prescribed - * aniso metrics \a m0 and \a m1. + * Compute the curve length of surface edge \f$[np0;np1]\f$ according to the + * prescribed aniso metrics \a m0 and \a m1. * + * \remark the edge has to be a boundary edge */ static inline double MMG5_lenEdg(MMG5_pMesh mesh,MMG5_int np0,MMG5_int np1, @@ -83,7 +84,6 @@ double MMG5_lenEdg(MMG5_pMesh mesh,MMG5_int np0,MMG5_int np1, n2 = &mesh->xpoint[p0->xp].n2[0]; ps1 = ux*n1[0] + uy*n1[1] + uz*n1[2]; ps2 = ux*n2[0] + uy*n2[1] + uz*n2[2]; - if ( fabs(ps2) < fabs(ps1) ) { n1 = &mesh->xpoint[p0->xp].n2[0]; ps1 = ps2; @@ -188,10 +188,12 @@ double MMG5_lenEdg(MMG5_pMesh mesh,MMG5_int np0,MMG5_int np1, * \param isedg 1 if the edge is a ridge, 0 otherwise. * \return length of edge according to the prescribed metric, 0 if fail. * - * Compute length of surface edge \f$[i0;i1]\f$ according to the prescribed - * aniso metric (for special storage of metrics at ridges points). Here the - * length is computed taking into account the curve nature of the surface edge. + * Compute the curve length of surface edge \f$[i0;i1]\f$ according to the + * prescribed aniso metric (for special storage of metrics at ridges + * points). Here the length is computed taking into account the curve nature of + * the surface edge. * + * \remark the edge has to be a boundary edge */ static inline double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np0,MMG5_int np1,int8_t isedg) { @@ -257,9 +259,10 @@ double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np0,MMG5_int n * \param isedg 1 if the edge is a ridge, 0 otherwise. * \return length of edge according to the prescribed metric. * - * Compute length of surface edge \f$[i0;i1]\f$ according to the prescribed - * aniso metric (for classic storage of metrics at ridges points). + * Compute the curve length of surface edge \f$[i0;i1]\f$ according to the + * prescribed aniso metric (for classic storage of metrics at ridges points). * + * \remark the edge has to be a boundary edge */ static inline double MMG5_lenSurfEdg33_ani(MMG5_pMesh mesh,MMG5_pSol met, @@ -282,8 +285,8 @@ double MMG5_lenSurfEdg33_ani(MMG5_pMesh mesh,MMG5_pSol met, * compatibility with \a lenedg_ani). * \return length of edge according to the prescribed metric. * - * Compute length of surface edge \f$[i0;i1]\f$ according to the prescribed iso - * metric. + * Compute the "straight" length of edge \f$[i0;i1]\f$ according to the + * prescribed iso metric. * */ static diff --git a/src/mmg3d/colver_3d.c b/src/mmg3d/colver_3d.c index c7eefa722..f31f46606 100644 --- a/src/mmg3d/colver_3d.c +++ b/src/mmg3d/colver_3d.c @@ -161,6 +161,10 @@ int MMG5_chkcol_int(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t iface, /* Rough evaluation of edge length (doesn't take into account if some of * the modified edges of pt0 are boundaries): for a more precise * computation, we need to update the edge tags of pt0. */ + + // Algiane 06/24: to check and fix: If the edge is marked as MG_BDY + // (keeping in mind that this tag is not always consistent), the curve + // length is computed, I don't know if it is wanted. ll = MMG5_lenedgspl(mesh,met,jj,pt0); if ( (!ll) || (ll > lon) )//LOPTL too small, we need to put greater than 1.41 return 0; diff --git a/src/mmg3d/inlined_functions_3d_private.h b/src/mmg3d/inlined_functions_3d_private.h index 8e0fe5168..9ebd22940 100644 --- a/src/mmg3d/inlined_functions_3d_private.h +++ b/src/mmg3d/inlined_functions_3d_private.h @@ -89,6 +89,9 @@ inline double MMG5_lenedgCoor_ani(double *ca,double *cb,double *sa,double *sb) { * Compute length of edge \f$[i0;i1]\f$ according to the prescribed aniso * metric (for classic storage of metrics at ridges points). * + * \warning in this function we may erroneously approximate the length of a + * curve boundary edge by the length of the straight edge if the "MG_BDY" tag is + * missing along the edge. */ static inline double MMG5_lenedg33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, @@ -102,8 +105,10 @@ inline double MMG5_lenedg33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_BDY)) { isedg = ( mesh->xtetra[pt->xt].tag[ia] & MG_GEO); + // Computation of the length of a curve edge with 33 aniso metric return MMG5_lenSurfEdg33_ani(mesh, met, ip1, ip2, isedg); } else { + // Computation for an internal edge with 33 aniso metric return MMG5_lenedgCoor_ani(mesh->point[ip1].c,mesh->point[ip2].c, &met->m[6*ip1],&met->m[6*ip2]); } @@ -150,8 +155,8 @@ inline double MMG5_lenedgspl33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, * \param pt pointer to the tetra from which we come. * \return length of edge according to the prescribed metric, 0 if fail. * - * Compute length of edge \f$[i0;i1]\f$ according to the prescribed aniso - * metric (for special storage of metrics at ridges points). + * Compute length of a straight edge \f$[i0;i1]\f$ according to the prescribed + * aniso metric (for special storage of metrics at ridges points). * */ static @@ -196,6 +201,10 @@ inline double MMG5_lenedgspl_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, * Compute length of edge \f$[i0;i1]\f$ according to the prescribed aniso * metric (for special storage of metrics at ridges points). * + * \warning in this function we may erroneously approximate the length of a + * curve boundary edge by the length of the straight edge if the "MG_BDY" tag is + * missing along the edge. + * */ static inline double MMG5_lenedg_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, @@ -209,8 +218,10 @@ inline double MMG5_lenedg_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia, if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_BDY)) { isedg = ( mesh->xtetra[pt->xt].tag[ia] & MG_GEO); + // Computation of the length of a curve edge with ridge metric return MMG5_lenSurfEdg_ani(mesh, met, ip1, ip2, isedg); } else { + // Computation for an internal edge with ridge metric return MMG5_lenedgspl_ani(mesh ,met, ia, pt); } return 0.0; diff --git a/src/mmg3d/libmmg3d_tools.c b/src/mmg3d/libmmg3d_tools.c index 977a495ca..43fa2b7d0 100644 --- a/src/mmg3d/libmmg3d_tools.c +++ b/src/mmg3d/libmmg3d_tools.c @@ -57,19 +57,39 @@ void MMG3D_setfunc(MMG5_pMesh mesh,MMG5_pSol met) { MMG5_caltet = MMG5_caltet_iso; MMG5_caltri = MMG5_caltri_iso; MMG3D_doSol = MMG3D_doSol_iso; + // same as MMG5_lenSurfEdg_iso for iso metric. The edge can be boundary or + // intenal but the test relies on the MG_BDY tag that may be missing along + // boundary edges (it doesn't matter in iso mode as we always compute the + // "straight" edge length). It starts from tetra pointer and edge index. MMG5_lenedg = MMG5_lenedg_iso; + // has to be used to compute "straight" edge length from coordinates MMG3D_lenedgCoor = MMG5_lenedgCoor_iso; + // Straight edge length (edge is guessed to be a surface edge) from point indices MMG5_lenSurfEdg = MMG5_lenSurfEdg_iso; } else { MMG5_caltet = MMG5_caltet_ani; MMG5_caltri = MMG5_caltri_ani; MMG3D_doSol = MMG3D_doSol_ani; + // lenedg is meant to compute the curve length along boundary edges + // and the straight length for internal edges (from + // tetra pointer and an edge index) but it relies + // on the MG_BDY tag which is not always consistent to detect boundary + // edges. Moreover, it seems that the "straight" length is computed in iso + // mode - the origin and effect of computing curve lengths along boudary + // edges should be investigated... MMG5_lenedg = MMG5_lenedg_ani; + // lenedgCoor has to be used to compute "straight" edge length from + // coordinates MMG3D_lenedgCoor = MMG5_lenedgCoor_ani; + // lenSurfEdg can be called only from a boundary edge: curve length for + // aniso metric from point indices MMG5_lenSurfEdg = MMG5_lenSurfEdg_ani; } MMG5_intmet = MMG5_intmet_ani; + // warning the lenedg_ani function we may erroneously approximate the length + // of a curve boundary edge by the length of the straight edge if the + // "MG_BDY" tag is missing along the edge. MMG5_lenedgspl = MMG5_lenedg_ani; MMG5_movintpt = MMG5_movintpt_ani; MMG5_movbdyregpt = MMG5_movbdyregpt_ani; @@ -97,8 +117,14 @@ void MMG3D_setfunc(MMG5_pMesh mesh,MMG5_pSol met) { } MMG5_caltri = MMG5_caltri_iso; MMG3D_doSol = MMG3D_doSol_iso; + // same as MMG5_lenSurfEdg_iso for iso metric. The edge can be boundary or + // intenal but the test relies on the MG_BDY tag that may be missing along + // boundary edges (it doesn't matter in iso mode as we always compute the + // "straight" edge length). It starts from tetra pointer and edge index. MMG5_lenedg = MMG5_lenedg_iso; + // lenedgCoor has to be used to compute "straight" edge length from coordinates MMG3D_lenedgCoor = MMG5_lenedgCoor_iso; + // Straight edge length (edge is guessed to be a surface edge) from point indices MMG5_lenSurfEdg = MMG5_lenSurfEdg_iso; MMG5_intmet = MMG5_intmet_iso; MMG5_lenedgspl = MMG5_lenedg_iso; @@ -1402,9 +1428,15 @@ int MMG3D_searchlen(MMG5_pMesh mesh, MMG5_pSol met, double lmin, ier = MMG5_hashPop(&hash,np,nq); if( ier ) { if ( (!metRidTyp) && met->m && met->size>1 ) { + // Warning: for aniso metrics, we may erroneously approximate the + // length of a curve boundary edge by the length of the straight edge + // if the "MG_BDY" tag is missing along the edge. len = MMG5_lenedg(mesh,met,ia,pt); } else { + // Warning: we may erroneously approximate the length of a curve + // boundary edge by the length of the straight edge if the "MG_BDY" + // tag is missing along the edge. len = MMG5_lenedg33_ani(mesh,met,ia,pt); } diff --git a/src/mmg3d/quality_3d.c b/src/mmg3d/quality_3d.c index 40aa7732e..4f610839a 100644 --- a/src/mmg3d/quality_3d.c +++ b/src/mmg3d/quality_3d.c @@ -295,9 +295,15 @@ int MMG3D_computePrilen( MMG5_pMesh mesh, MMG5_pSol met, double* avlen, ier = MMG5_hashPop(&hash,np,nq); if( ier ) { if ( (!metRidTyp) && met->size==6 && met->m ) { + // Warning: we may erroneously approximate the length of a curve + // boundary edge by the length of the straight edge if the "MG_BDY" + // tag is missing along the edge. len = MMG5_lenedg33_ani(mesh,met,ia,pt); } else + // Warning: we may erroneously approximate the length of a curve + // boundary edge by the length of the straight edge if the "MG_BDY" + // tag is missing along the edge. len = MMG5_lenedg(mesh,met,ia,pt); diff --git a/src/mmg3d/split_3d.c b/src/mmg3d/split_3d.c index 20d6bb4f0..fbbe0a6ac 100644 --- a/src/mmg3d/split_3d.c +++ b/src/mmg3d/split_3d.c @@ -648,8 +648,14 @@ int MMG5_split1b(MMG5_pMesh mesh, MMG5_pSol met,int64_t *list, int ret, MMG5_int for (i=0; i<6; i++) { pt = &mesh->tetra[list[j]/6]; if ( (!metRidTyp) && met->m && met->size>1 ) + // Warning: we may erroneously approximate the length of a curve + // boundary edge by the length of the straight edge if the "MG_BDY" + // tag is missing along the edge. len = MMG5_lenedg33_ani(mesh,met,i,pt); else + // Warning: for aniso metrics we may erroneously approximate the + // length of a curve boundary edge by the length of the straight edge + // if the "MG_BDY" tag is missing along the edge. len = MMG5_lenedg(mesh,met,i,pt); if ( len < lmin) { lmin = len; @@ -684,9 +690,17 @@ int MMG5_split1b(MMG5_pMesh mesh, MMG5_pSol met,int64_t *list, int ret, MMG5_int } } if ( (!metRidTyp) && met->m && met->size>1 ) + /* Computation of straight edge length */ len = MMG5_lenedgspl33_ani(mesh,met,taued[5],pt0); else + /* if edge is marked MG_BDY, curve length is computed, if edge is not + * MG_BDY, straight length is computed (even if the edge is indeed along + * the boundary but misses the tag). // Algiane 06/24: to check and fix + * or comment (I don't know if it is useful to compute the accurate + * curve length here) + */ len = MMG5_lenedgspl(mesh,met,taued[5],pt0); + if ( len < lmin ) break; memcpy(pt0,pt,sizeof(MMG5_Tetra)); @@ -698,8 +712,15 @@ int MMG5_split1b(MMG5_pMesh mesh, MMG5_pSol met,int64_t *list, int ret, MMG5_int } if ( (!metRidTyp) && met->m && met->size>1 ) + /* Computation of straight edge length */ len = MMG5_lenedgspl33_ani(mesh,met,taued[5],pt0); else + /* if edge is marked MG_BDY, curve length is computed, if edge is not + * MG_BDY, straight length is computed (even if the edge is indeed along + * the boundary but misses the tag). // Algiane 06/24: to check and fix + * or comment (I don't know if it is useful to compute the accurate + * curve length here) + */ len = MMG5_lenedgspl(mesh,met,taued[5],pt0); if ( len < lmin ) break; }