From 94649905b711d36fcf8eb70dad4578a9a890ca09 Mon Sep 17 00:00:00 2001
From: Corentin Prigent <corentin.prigent@inria.fr>
Date: Wed, 28 May 2025 11:04:35 +0200
Subject: [PATCH 1/2] remove usage of rand() and srand() functions to have
 deterministic output

---
 src/mmg2d/enforcement_2d.c | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/mmg2d/enforcement_2d.c b/src/mmg2d/enforcement_2d.c
index 490c7929f..24b0f63ce 100644
--- a/src/mmg2d/enforcement_2d.c
+++ b/src/mmg2d/enforcement_2d.c
@@ -38,7 +38,7 @@ int MMG2D_bdryenforcement(MMG5_pMesh mesh,MMG5_pSol sol) {
   MMG5_int        list[MMG5_TRIA_LMAX],list2[3];
   MMG5_int        k,l,kk,nex,kdep,lon,iel;
   int             iare,ied;
-  MMG5_int        ia,ib,ilon,rnd,idep,*adja,ir,adj;
+  MMG5_int        ia,ib,ilon,ii,idep,*adja,ir,adj;
   int8_t          i,i1,i2,j;
 //  int       iadr2,*adja2,ndel,iadr,ped0,ped1;
   static int8_t   mmgWarn0=0,mmgWarn1=0,mmgWarn2=0,mmgWarn3=0;
@@ -183,11 +183,10 @@ int MMG2D_bdryenforcement(MMG5_pMesh mesh,MMG5_pSol sol) {
 
       /* Randomly swap edges in the list, while... */
       ilon = lon;
-      srand(time(NULL));
+      ii = 0;
 
       while ( ilon > 0 ) {
-        rnd = ( rand() % lon );
-        k   = list[rnd] / 3;
+        k   = list[ii] / 3;
 
         /* Check for triangle k in the pipe until one triangle with base+1 is met (indicating that it is
            crossed by the considered edge) */
@@ -195,11 +194,11 @@ 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[(++rnd)%lon] / 3;
+          k = list[(++ii)%lon] / 3;
         }
 
         assert ( l <= lon );
-        idep = list[rnd] % 3;
+        idep = list[ii] % 3;
         // if(mesh->info.ddebug) printf("i= %" MMG5_PRId " < %" MMG5_PRId " ? on demarre avec %" MMG5_PRId "\n",i,lon+1,k);
         adja = &mesh->adja[3*(k-1)+1];
 
@@ -248,7 +247,8 @@ int MMG2D_bdryenforcement(MMG5_pMesh mesh,MMG5_pSol sol) {
           }
           break;
         }
-     }
+        ii++;
+      }
     }
   }
 

From 6e32ef53820e9260bf38f3a58f1e01299872fc38 Mon Sep 17 00:00:00 2001
From: Corentin Prigent <corentin.prigent@inria.fr>
Date: Fri, 18 Jul 2025 19:02:12 +0200
Subject: [PATCH 2/2] Do not allow creation of flat triangles in edge
 enforcement. Modify threshold for float comparison that could create false
 positive. Fix index bug

---
 src/mmg2d/enforcement_2d.c | 5 +++--
 src/mmg2d/locate_2d.c      | 8 ++++----
 src/mmg2d/swapar_2d.c      | 5 +++++
 3 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/src/mmg2d/enforcement_2d.c b/src/mmg2d/enforcement_2d.c
index 24b0f63ce..93b2ce1ec 100644
--- a/src/mmg2d/enforcement_2d.c
+++ b/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;
       }
     }
   }
diff --git a/src/mmg2d/locate_2d.c b/src/mmg2d/locate_2d.c
index 0a734b1e2..5dda9a3de 100644
--- a/src/mmg2d/locate_2d.c
+++ b/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;
         }
diff --git a/src/mmg2d/swapar_2d.c b/src/mmg2d/swapar_2d.c
index 1522d63b6..f10d2d0fe 100644
--- a/src/mmg2d/swapar_2d.c
+++ b/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;