Add comment to highlight the inconsistencies in edge lengths computat…

…ions.
MmgTools · May 31, 2024 · 5a73683 · 5a73683
1 parent f09cee8
commit 5a73683
Show file tree

Hide file tree

Showing 6 changed files with 90 additions and 13 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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;
     }