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;