Do not allow creation of flat triangles in edge enforcement. Modify threshold for float comparison that could create false positive. Fix index bug

Author: corentin-prigent

src/mmg2d/enforcement_2d.c

@@ -194,7 +194,8 @@ int MMG2D_bdryenforcement(MMG5_pMesh mesh,MMG5_pSol sol) {
         while ( l++ < lon ) {
           pt = &mesh->tria[k];
           if ( pt->base == mesh->base+1 ) break;
-          k = list[(++ii)%lon] / 3;
+          ii = (ii+1)%lon;
+          k = list[ii] / 3;
         }
 
         assert ( l <= lon );
@@ -247,7 +248,7 @@ int MMG2D_bdryenforcement(MMG5_pMesh mesh,MMG5_pSol sol) {
           }
           break;
         }
-        ii++;
+        ii = (ii+1)%lon;
       }
     }
   }

src/mmg2d/locate_2d.c

@@ -167,21 +167,21 @@ int MMG2D_cutEdgeTriangle(MMG5_pMesh mesh,MMG5_int k,MMG5_int ia,MMG5_int ib) {
   prod3 = area3*area1;
 
   /* Both edges p2p3 and p1p3 corresponding to prod2 and prod3 are cut by edge ia,ib */
-  if ( prod1 > 0 && ((prod2 < 0 || prod3 < 0))) {
+  if ( prod1 > MMG2D_EPSD && ((prod2 < -MMG2D_EPSD || prod3 < -MMG2D_EPSD ))) {
     if ( (iare = MMG2D_cutEdge(mesh,pt,ppa,ppb)) ) {
       return iare;
     }
   }
 
   /* Both edges corresponding to prod1 and prod3 are cut by edge ia,ib */
-  if ( prod2 > 0 && ((prod1 < 0 || prod3 < 0))) {
+  if ( prod2 > MMG2D_EPSD && ((prod1 < -MMG2D_EPSD || prod3 < -MMG2D_EPSD ))) {
     if ( (iare = MMG2D_cutEdge(mesh,pt,ppa,ppb)) ) {
       return iare;
     }
   }
 
   /* Both edges corresponding to prod1 and prod2 are cut by edge ia,ib */
-  if ( prod3 > 0 && ((prod2 < 0 || prod1 < 0))) {
+  if ( prod3 > MMG2D_EPSD && ((prod2 < -MMG2D_EPSD || prod1 < -MMG2D_EPSD ))) {
     if ( (iare = MMG2D_cutEdge(mesh,pt,ppa,ppb)) ) {
       return iare;
     }
@@ -191,7 +191,7 @@ int MMG2D_cutEdgeTriangle(MMG5_pMesh mesh,MMG5_int k,MMG5_int ia,MMG5_int ib) {
   for(i=0; i<3; i++){
     if ( pt->v[i] == ia || ibreak ) {
       /* One vertex is ia, and the opposite edge is frankly crossed */
-      if ( (prod1 < 0) || (prod2 < 0) || (prod3 < 0) ) {
+      if ( (prod1 < -MMG2D_EPSD) || (prod2 < -MMG2D_EPSD) || (prod3 < -MMG2D_EPSD) ) {
         if ( (iare = MMG2D_cutEdge(mesh,pt,ppa,ppb)) ) {
           return iare;
         }

src/mmg2d/swapar_2d.c

@@ -78,6 +78,11 @@ int MMG2D_swapdelone(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_int k,int8_t i,double cr
   arean2 = MMG2D_quickarea(mesh->point[pt0->v[0]].c,mesh->point[pt0->v[1]].c,mesh->point[pt0->v[2]].c);
   if ( cal2 > crit )  return 0;
 
+  /* Avoid creation of a flat triangle*/
+  if (fabs(arean1) < MMG2D_EPSD || fabs(arean2) < MMG2D_EPSD) {
+    return 0;
+  }
+
   if ( arean1 < 0.0 || arean2 < 0.0 || fabs((area1+area2)-(arean1+arean2)) > MMG2D_EPSD ) {
     if(mesh->info.ddebug) printf("  ## Warning: non convex configuration\n");
     return 0;