Additional option to remesh only the elements of bad quality

src/common/libmmgtypes.h

@@ -524,9 +524,13 @@ typedef struct {
   MMG5_pPar     par;
   double        dhd,hmin,hmax,hsiz,hgrad,hgradreq,hausd;
   double        min[3],max[3],delta,ls,lxreg,rmc;
+  double        limit_angle; /*< Angle threshold for modifying triangles or not */
+  double        relative_min_tolerance, relative_max_tolerance; /* < Tolerance on the minimal and maximal element sizes */
   MMG5_int      *br; /*!< list of based references to which an implicit surface can be attached */
   MMG5_int      isoref; /*!< isovalue reference in ls mode */
   MMG5_int      nsd; /*!< index of subdomain to save (0 by default == all subdomains are saved) */
+  int           isotropic_pt_relocation; /*!< force the use of the isotropic point relocation */
+  int           bdy_adaptation; /*!< extend the remeshing close to boundaries when limit_angle is activated */
   int           mem,npar,npari;
   int           nbr,nbri; /*!< number of based references for level-set (BC to which a material can be attached) */
   int           opnbdy; /*!< floating surfaces */
@@ -646,6 +650,7 @@ typedef struct {
                     \f$k^th\f$ quadrilaterals are adjacent and share their
                     edges \a j and \a l (resp.) */
   int       *ipar;  /*!< Store indices of the local parameters */
+  double    *lagrangian_velocity; /*!< Velocity of the vertices when Lagrangian resolution */
   MMG5_pPoint    point; /*!< Pointer toward the \ref MMG5_Point structure */
   MMG5_pxPoint   xpoint; /*!< Pointer toward the \ref MMG5_xPoint structure */
   MMG5_pTetra    tetra; /*!< Pointer toward the \ref MMG5_Tetra structure */

src/mmg2d/API_functions_2d.c

@@ -101,11 +101,19 @@ void MMG2D_Init_parameters(MMG5_pMesh mesh) {
   /* default values for doubles */
   /* level set value */
   mesh->info.ls       = MMG5_LS;
-/* xreg relaxation parameter value */
+  /* xreg relaxation parameter value */
   mesh->info.lxreg    = MMG5_XREG;
-
+  /* [0/1]   ,avoid/enforce isotropic smoothing even with anisotropic metric */
+  mesh->info.isotropic_pt_relocation = MMG5_OFF;
+  /* limit angle to avoid remeshing some good triangles */
+  mesh->info.limit_angle = -1.; // Deactivated when negative or > PI/3
+  /* Tolerances on the element sizes when the limit_angle is activated */
+  mesh->info.relative_min_tolerance = 1.e10;
+  mesh->info.relative_max_tolerance = 1.e10;
   /* Ridge detection */
   mesh->info.dhd      = MMG5_ANGEDG;
+  /* to adapt more thoroughly close to boundaries */
+  mesh->info.bdy_adaptation = MMG5_OFF;
 }
 
 
@@ -147,6 +155,9 @@ int MMG2D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam, MMG5_int va
       mesh->info.dhd    = MMG5_ANGEDG;
     }
     break;
+  case MMG2D_IPARAM_bdy_adaptation:
+    mesh->info.bdy_adaptation = val;
+    break;
   case MMG2D_IPARAM_nofem :
     mesh->info.setfem = (val==1)? 0 : 1;
     break;
@@ -162,6 +173,9 @@ int MMG2D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam, MMG5_int va
   case MMG2D_IPARAM_isosurf :
     mesh->info.isosurf = val;
     break;
+  case MMG2D_IPARAM_isotropic_smoothing :
+    mesh->info.isotropic_pt_relocation = val;
+    break;
   case MMG2D_IPARAM_lag :
 #ifdef USE_ELAS
     if ( val < 0 || val > 2 )
@@ -349,6 +363,15 @@ int MMG2D_Set_dparameter(MMG5_pMesh mesh, MMG5_pSol sol, int dparam, double val)
     else
       mesh->info.hausd    = val;
     break;
+  case MMG2D_DPARAM_hmin_factor :
+    mesh->info.relative_min_tolerance = val;
+    break;
+  case MMG2D_DPARAM_hmax_factor :
+    mesh->info.relative_max_tolerance = val;
+    break;
+  case MMG2D_DPARAM_limit_angle :
+    mesh->info.limit_angle = val;
+    break;
   case MMG2D_DPARAM_ls :
     mesh->info.ls       = val;
     break;

src/mmg2d/libmmg2d.h

@@ -138,11 +138,16 @@ extern "C" {
     MMG2D_DPARAM_hausd,             /*!< [val], Global Hausdorff distance (on all the boundary surfaces of the mesh) */
     MMG2D_DPARAM_hgrad,             /*!< [val], Global gradation */
     MMG2D_DPARAM_hgradreq,          /*!< [val], Global gradation on required entites (advanced usage) */
+    MMG2D_DPARAM_hmin_factor,       /*!< [val], Factor on minimal edge length when limit_angle is active */
+    MMG2D_DPARAM_hmax_factor,       /*!< [val], Factor on maximal edge length when limit_angle is active */
     MMG2D_DPARAM_ls,                /*!< [val], Function value where the level set is to be discretized */
     MMG2D_DPARAM_xreg,              /*!< [val], Relaxation parameter for coordinate regularization (0<val<1) */
     MMG2D_DPARAM_rmc,               /*!< [-1/val], Remove small disconnected components in level-set mode */
     MMG2D_IPARAM_nofem,             /*!< [1/0], Do not attempt to make the mesh suitable for finite-element computations */
     MMG2D_IPARAM_isoref,            /*!< [0/n], Iso-surface boundary material reference */
+    MMG2D_IPARAM_isotropic_smoothing,       /*!< [0/1], Avoid/enforce isotropic smoothing even with anisotropic metric */
+    MMG2D_DPARAM_limit_angle,       /*!< [val], Minimal angle in triangles under which remeshing is achieved */
+    MMG2D_IPARAM_bdy_adaptation,    /*!< [1/0], Enable thorough adaptation close to the boundaries (if limit_angle < Pi) */
   };
 
 /*----------------------------- function headers -----------------------------*/

src/mmg2d/libmmg2d_tools.c

@@ -148,6 +148,11 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
         if ( !MMG2D_Set_solSize(mesh,met,MMG5_Vertex,0,MMG5_Tensor) )
           return 0;
         break;
+      case 'b':
+        if ( !strcmp(argv[i],"-bdyadaptation") ) {
+          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_bdy_adaptation,1) )
+            return 0;
+        }
       case 'd':
         if ( !strcmp(argv[i],"-default") ) {
           mesh->mark=1;
@@ -194,6 +199,12 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
           else if ( !strcmp(argv[i],"-hgrad") ) {
             param = MMG2D_DPARAM_hgrad;
           }
+          else if ( !strcmp(argv[i],"-hminfactor") ) {
+            param = MMG2D_DPARAM_hmin_factor;
+          }
+          else if ( !strcmp(argv[i],"-hmaxfactor") ) {
+            param = MMG2D_DPARAM_hmax_factor;
+          }
           else {
             /* Arg unknown by Mmg: arg starts with -h but is not known */
             MMG2D_usage(argv[0]);
@@ -235,6 +246,10 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
                                      atoi(argv[i])) )
             return 0;
         }
+        else if ( !strcmp(argv[i],"-isotropicsmoothing") ) {
+          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_isotropic_smoothing,1) )
+            return 0;
+        }
         else {
           MMG2D_usage(argv[0]);
           return 0;
@@ -278,6 +293,10 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
             }
           }
         }
+        else if ( !strcmp(argv[i],"-limitangle") && ++i < argc ) {
+          if ( !MMG2D_Set_dparameter(mesh,met,MMG2D_DPARAM_limit_angle,atof(argv[i])) )
+            return 0;
+        }
         break;
       case 'm':  /* memory */
         if ( !strcmp(argv[i],"-met") ) {

src/mmg2d/mmg2d1.c

@@ -50,7 +50,7 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
     if ( typchk == 2 && it == 0 ) {
 // #warning Luca: check consistency with 3D
     }
-    if ( !mesh->info.noinsert ) {
+    if ( !mesh->info.noinsert && (mesh->info.limit_angle >= M_PI/3. || mesh->info.limit_angle <= 0.)) {
       /* Memory free */
       MMG5_DEL_MEM(mesh,mesh->adja);
       mesh->adja = 0;
@@ -107,6 +107,158 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
   return 1;
 }
 
+/**
+ * \param mesh pointer to the mesh structure.
+ * \param k triangle index.
+ * \return the minimal angle of a triangle or 0 if the triangle size
+ * is not conform to the tolerance sizes.
+ *
+ * Compute the minimal angle of triangle k. If the velocity of each
+ * node is given, the extrapolated minimal angle of the triangle
+ * is also calculated. If the edge lengths are too small of too
+ * large, 0 is returned as we want this triangle to be remeshed.
+ *
+ */
+double MMG2D_minangle(MMG5_pMesh mesh, int k) {
+
+  MMG5_pTria pt = &mesh->tria[k];
+  MMG5_pPoint pp1, pp2, pp3;
+  int* adja = &mesh->adja[3*(k-1)+1];
+  int p1, p2, p3;
+  double length[6];
+  double c1 = 1.35*mesh->info.limit_angle;
+  if (c1 > 1.) c1 = 1.;
+  double coef = 1.;
+  double factor_max = mesh->info.relative_max_tolerance;
+  double factor_min = mesh->info.relative_min_tolerance;
+
+  p1 = pt->v[0];
+  p2 = pt->v[1];
+  p3 = pt->v[2];
+
+  pp1 = &mesh->point[p1];
+  pp2 = &mesh->point[p2];
+  pp3 = &mesh->point[p3];
+
+  double square_root_area = sqrt(0.5*fabs((pp2->c[0]-pp1->c[0])*(pp3->c[1]-pp1->c[1]) -
+                                          (pp2->c[1]-pp1->c[1])*(pp3->c[0]-pp1->c[0])));
+  double mean_velocity = 0.;
+  double* velocity = &mesh->lagrangian_velocity[0];
+
+  // Velocity is only considered for already existing points
+  if (velocity && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark
+               && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1)
+    mean_velocity = ( sqrt(pow(velocity[p1-1],2) + pow(velocity[p1-1+mesh->mark],2)) + 
+                      sqrt(pow(velocity[p2-1],2) + pow(velocity[p2-1+mesh->mark],2)) + 
+                      sqrt(pow(velocity[p3-1],2) + pow(velocity[p3-1+mesh->mark],2))  )/3.;
+
+  double dt, x1, x2, x3, y1, y2, y3;
+  if (mean_velocity < 1.e12) {
+    dt = 0.;
+    x1 = pp1->c[0];
+    x2 = pp2->c[0];
+    x3 = pp3->c[0];
+    y1 = pp1->c[1];
+    y2 = pp2->c[1];
+    y3 = pp3->c[1];
+  }
+  else {
+    // Extrapolated triangle location based on a fictitious time step
+    dt = 0.25*square_root_area / mean_velocity;
+    x1 = pp1->c[0] + velocity[p1-1]*dt;
+    x2 = pp2->c[0] + velocity[p2-1]*dt;
+    x3 = pp3->c[0] + velocity[p3-1]*dt;
+    y1 = pp1->c[1] + velocity[p1-1+mesh->mark]*dt;
+    y2 = pp2->c[1] + velocity[p2-1+mesh->mark]*dt;
+    y3 = pp3->c[1] + velocity[p3-1+mesh->mark]*dt;
+  }
+
+  // Compute the length of the triangle edges
+  length[0] = pow(x1-x2,2.) + pow(y1-y2,2.);
+  length[1] = pow(x2-x3,2.) + pow(y2-y3,2.);
+  length[2] = pow(x3-x1,2.) + pow(y3-y1,2.);
+
+  // Compute the length of the fictitious triangle edges
+  length[3] = pow(pp1->c[0] - pp2->c[0],2.) + pow(pp1->c[1] - pp2->c[1],2.);
+  length[4] = pow(pp2->c[0] - pp3->c[0],2.) + pow(pp2->c[1] - pp3->c[1],2.);
+  length[5] = pow(pp3->c[0] - pp1->c[0],2.) + pow(pp3->c[1] - pp1->c[1],2.);
+
+  // Check if the edge lengths satisfy the extended tolerances
+  double min = length[0];
+  double max = length[0];
+  int min_index = 0;
+
+  for (int i = 1; i < 3; i++) {
+    if (length[i] < min) {
+      min = length[i];
+      min_index = i;
+    }
+    if (length[i] > max) max = length[i];
+  }
+
+  if (sqrt(min) < mesh->info.hmin/factor_min) return 0.;
+  if (sqrt(max) > factor_max*mesh->info.hmax) return 0.;
+
+  double cosA = ( length[(min_index+1)%3] + length[(min_index+2)%3] - length[min_index] ) / 
+                (2*pow(length[(min_index+1)%3],0.5)*pow(length[(min_index+2)%3],0.5));
+
+  min = length[3];
+  max = length[3];
+  min_index = 0;
+  for (int i = 1; i < 3; i++) {
+    if (length[i+3] < min) {
+      min = length[i+3];
+      min_index = i;
+    }
+    if (length[i+3] > max) max = length[i+3];
+  }
+
+  if (sqrt(min) < mesh->info.hmin/factor_min) return 0.;
+  if (sqrt(max) > factor_max*mesh->info.hmax) return 0.;
+
+  // If the triangle is on a boundary, the angle tolerance is divided by two
+  // to extend remeshing close to boundaries where the mesh can be often squeezed
+  if (mesh->info.bdy_adaptation) {
+    if ((pp1->tag & MG_BDY) || (pp2->tag & MG_BDY) || (pp3->tag & MG_BDY)) coef = c1;
+    else {
+      // The neighbours of boundary triangles also have a reduced tolerance
+      for(int i = 0; i < 3; i++ ) {
+        int k1 = adja[i]/3;
+
+        if (!k1) continue;
+
+        MMG5_pTria pt_neigh = &mesh->tria[k1];
+
+        for (int j = 0; j < 3; j++) {
+          MMG5_pPoint pp = &mesh->point[pt_neigh->v[j]];
+          if (pp->tag & MG_BDY) {
+            coef = c1;
+            break;
+          }
+        }
+
+        if (coef < 0.999) break;
+      }
+    }
+  }
+
+  // Minimal angle
+  double cosA2 = ( length[(min_index+1)%3+3] + length[(min_index+2)%3+3] - length[min_index+3] ) / 
+                 (2*pow(length[(min_index+1)%3+3],0.5)*pow(length[(min_index+2)%3+3],0.5));
+
+  if (cosA2 > cosA) cosA = cosA2;
+
+  if (cosA <= -1.) {
+    return 3.14;
+  }
+  else if (cosA >= 1) {
+    return 0.;
+  }
+  else return coef*acos(cosA);
+
+}
+
+
 /* Travel triangles and split long edges according to patterns */
 MMG5_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
   MMG5_pTria      pt;
@@ -130,6 +282,8 @@ MMG5_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
     if ( !MG_EOK(pt) || (pt->ref < 0) ) continue;
     if ( MG_SIN(pt->tag[0]) || MG_SIN(pt->tag[1]) || MG_SIN(pt->tag[2]) )  continue;
 
+    if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue;
+
     /* Check if pt should be cut */
     pt->flag = 0;
 
@@ -450,6 +604,8 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax) {
     pt = &mesh->tria[k];
     if ( !MG_EOK(pt) || pt->ref < 0 ) continue;
 
+    if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue;
+
     /* Travel 3 edges of the triangle and decide whether to collapse p1->p2, based on length criterion */
     pt->flag = 0; // was here before, but probably serves for nothing
 
@@ -523,12 +679,19 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
   maxit = 2;
   mesh->base++;
 
+  // If some triangles are not adapted using the limiting option on the triangle minimal angle, 
+  // using fem = 1 here can lead to very poor triangles that remain in the final mesh.
+  // So, in this case, we enforce the use of fem = 0.
+  if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3.) mesh->info.fem = 0;
+
   do {
     ns = 0;
     for (k=1; k<=mesh->nt; k++) {
       pt = &mesh->tria[k];
       if ( !MG_EOK(pt) || pt->ref < 0 ) continue;
 
+      if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue;
+
       for (i=0; i<3; i++) {
         if ( MG_SIN(pt->tag[i]) || MG_EDG(pt->tag[i]) ) continue;
         else if ( MMG2D_chkswp(mesh,met,k,i,typchk) ) {
@@ -542,6 +705,10 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
   while ( ns > 0 && ++it<maxit );
   if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nns > 0 )
     fprintf(stdout,"     %8" MMG5_PRId " edge swapped\n",nns);
+
+  // The fem value is restored after the swaps
+  if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3.) mesh->info.fem = mesh->info.setfem;
+
   return nns;
 }
 
@@ -653,6 +820,8 @@ MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) {
     pt = &mesh->tria[k];
     if ( !MG_EOK(pt) || pt->ref < 0 ) continue;
 
+    if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue;
+
     imax = -1;
     lmax = -1.0;
     for (i=0; i<3; i++) {
@@ -715,6 +884,8 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
       pt = &mesh->tria[k];
       if ( !MG_EOK(pt) || pt->ref < 0 ) continue;
 
+      if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue;
+
       for (i=0; i<3; i++) {
         p0 = &mesh->point[pt->v[i]];
         if ( p0->flag == base || MG_SIN(p0->tag) || p0->tag & MG_NOM ) continue;
@@ -726,6 +897,9 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
           ier = MMG2D_movedgpt(mesh,met,ilist,list,improve);
           if ( ier ) ns++;
         }
+        else if (mesh->info.isotropic_pt_relocation) {
+            ier = MMG2D_movintpt(mesh,met,ilist,list,improve);
+        }
         else {
           if ( met->size == 3 && met->m )
             ier = MMG2D_movintpt_ani(mesh,met,ilist,list,improve);
@@ -760,6 +934,9 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
  **/
 int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) {
 
+  mesh->mark = mesh->np;
+  for (int i = 1; i <= mesh->np; i++) mesh->point[i].tmp = -1;
+
   /* Stage 1: creation of a geometric mesh */
   if ( abs(mesh->info.imprim) > 4 || mesh->info.ddebug )
     fprintf(stdout,"  ** GEOMETRIC MESH\n");