diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h
index fe75d6193..5dfcdb121 100644
--- a/src/common/libmmgtypes.h
+++ b/src/common/libmmgtypes.h
@@ -524,9 +524,12 @@ 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 */
   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; /*!< force the use of some isotropic functions */
+  int           bdy_adaptation;
   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 */
diff --git a/src/common/scalem.c b/src/common/scalem.c
index 33155d91b..17da1b2e8 100644
--- a/src/common/scalem.c
+++ b/src/common/scalem.c
@@ -745,7 +745,7 @@ int MMG5_unscaleMesh(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol sol) {
   mesh->info.delta = 1.;
   mesh->info.min[0]= 0.;
   mesh->info.min[1]= 0.;
-  mesh->info.min[2]= 0.;
+  if (mesh->dim == 3) mesh->info.min[2]= 0.;
 
   /* de-normalize metric */
   if ( !(met && met->np && met->m) )  return 1;
diff --git a/src/mmg2d/API_functions_2d.c b/src/mmg2d/API_functions_2d.c
index 3205094b3..44574688f 100644
--- a/src/mmg2d/API_functions_2d.c
+++ b/src/mmg2d/API_functions_2d.c
@@ -101,11 +101,16 @@ 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 istropic remeshing even with anisotropic metric */
+  mesh->info.isotropic = MMG5_OFF;
+  /* limit angle to avoid remeshing some good triangles */
+  mesh->info.limit_angle = 5.*atan(1.);
   /* Ridge detection */
   mesh->info.dhd      = MMG5_ANGEDG;
+
+  mesh->info.bdy_adaptation = MMG5_OFF;
 }
 
 
@@ -147,6 +152,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 +170,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 :
+    mesh->info.isotropic = val;
+    break;
   case MMG2D_IPARAM_lag :
 #ifdef USE_ELAS
     if ( val < 0 || val > 2 )
@@ -349,6 +360,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.min[2] = val;
+    break;
+  case MMG2D_DPARAM_hmax_factor :
+    mesh->info.max[2] = val;
+    break;
+  case MMG2D_DPARAM_limit_angle :
+    mesh->info.limit_angle = val;
+    break;
   case MMG2D_DPARAM_ls :
     mesh->info.ls       = val;
     break;
diff --git a/src/mmg2d/libmmg2d.c b/src/mmg2d/libmmg2d.c
index 892b95297..2df7224f8 100644
--- a/src/mmg2d/libmmg2d.c
+++ b/src/mmg2d/libmmg2d.c
@@ -60,7 +60,7 @@ void MMG2D_Set_commonFunc(void) {
     return;
 }
 
-int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met) {
+int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) {
   MMG5_pSol sol=NULL; // unused
   mytime    ctim[TIMEMAX];
   char      stim[32];
@@ -219,7 +219,7 @@ int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met) {
     fprintf(stdout,"\n  -- PHASE 2 : %s MESHING\n",met->size < 3 ? "ISOTROPIC" : "ANISOTROPIC");
 
   /* Mesh improvement */
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
     MMG2D_RETURN_AND_PACK(mesh,met,sol,MMG5_LOWFAILURE);
   }
@@ -306,6 +306,7 @@ int MMG2D_mmg2dmesh(MMG5_pMesh mesh,MMG5_pSol met) {
   MMG5_pSol sol=NULL; // unused
   mytime    ctim[TIMEMAX];
   char      stim[32];
+  double* velocity;
 
   assert ( mesh );
   assert ( met );
@@ -479,7 +480,7 @@ int MMG2D_mmg2dmesh(MMG5_pMesh mesh,MMG5_pSol met) {
     fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT (%s)\n",
             met->size < 3 ? "ISOTROPIC" : "ANISOTROPIC");
 
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
     MMG2D_RETURN_AND_PACK(mesh,met,sol,MMG5_LOWFAILURE);
   }
@@ -526,6 +527,7 @@ int MMG2D_mmg2dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
   mytime    ctim[TIMEMAX];
   char      stim[32];
   int8_t    mettofree = 0;
+  double*   velocity;
 
   assert ( mesh );
   assert ( sol );
@@ -758,7 +760,7 @@ int MMG2D_mmg2dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
     fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT\n");
   }
 
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( mettofree ) {
       MMG5_DEL_MEM(mesh,met->m);
       MMG5_SAFE_FREE (met);
@@ -822,6 +824,7 @@ int MMG2D_mmg2dmov(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol disp) {
   char      stim[32];
   int       ier;
   MMG5_int  k,*invalidTris;
+  double*   velocity;
 
   assert ( mesh );
   assert ( disp );
@@ -984,7 +987,7 @@ int MMG2D_mmg2dmov(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol disp) {
       fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT\n");
     }
 
-    if ( !MMG2D_mmg2d1n(mesh,met) ) {
+    if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
       if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,disp,MMG5_STRONGFAILURE);
       MMG2D_RETURN_AND_PACK(mesh,met,disp,MMG5_LOWFAILURE);
     }
diff --git a/src/mmg2d/libmmg2d.h b/src/mmg2d/libmmg2d.h
index f0097d566..45775cdb9 100644
--- a/src/mmg2d/libmmg2d.h
+++ b/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,         /*!< [0/1], Avoid/enforce isotropic remeshing 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 -----------------------------*/
@@ -2293,7 +2298,7 @@ LIBMMG2D_EXPORT int MMG2D_Free_all(const int starter,...);
  * >   END SUBROUTINE\n
  *
  */
-  LIBMMG2D_EXPORT int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol sol);
+  LIBMMG2D_EXPORT int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol sol,double* velocity);
 
 /**
  * \brief Main "program" for the mesh generation library.
diff --git a/src/mmg2d/libmmg2d_private.h b/src/mmg2d/libmmg2d_private.h
index 631690264..cab514541 100644
--- a/src/mmg2d/libmmg2d_private.h
+++ b/src/mmg2d/libmmg2d_private.h
@@ -254,9 +254,9 @@ int MMG2D_norver(MMG5_pMesh,MMG5_int );
 int MMG2D_regnor(MMG5_pMesh );
 int MMG2D_regver(MMG5_pMesh );
 int MMG2D_boulen(MMG5_pMesh , MMG5_int ,int8_t ,MMG5_int *,MMG5_int *,double *);
-int MMG2D_mmg2d1n(MMG5_pMesh ,MMG5_pSol );
-int MMG2D_anatri(MMG5_pMesh ,MMG5_pSol ,int8_t );
-int MMG2D_adptri(MMG5_pMesh ,MMG5_pSol );
+int MMG2D_mmg2d1n(MMG5_pMesh ,MMG5_pSol, double*);
+int MMG2D_anatri(MMG5_pMesh ,MMG5_pSol ,int8_t , double*);
+int MMG2D_adptri(MMG5_pMesh ,MMG5_pSol, double*);
 int MMG2D_defsiz_iso(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_defsiz_ani(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_defmetbdy_2d(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t );
@@ -265,9 +265,9 @@ MMG5_int MMG2D_grad2met_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_pTria,MMG5_int,MMG5_int)
 int MMG2D_grad2metreq_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_pTria,MMG5_int,MMG5_int);
 int MMG2D_gradsiz_ani(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_gradsizreq_ani(MMG5_pMesh ,MMG5_pSol );
-MMG5_int MMG2D_anaelt(MMG5_pMesh ,MMG5_pSol ,int );
-MMG5_int MMG2D_colelt(MMG5_pMesh ,MMG5_pSol ,int, double );
-MMG5_int MMG2D_swpmsh(MMG5_pMesh ,MMG5_pSol ,int );
+MMG5_int MMG2D_anaelt(MMG5_pMesh ,MMG5_pSol ,int ,double*);
+MMG5_int MMG2D_colelt(MMG5_pMesh ,MMG5_pSol ,int ,double, double*);
+MMG5_int MMG2D_swpmsh(MMG5_pMesh ,MMG5_pSol ,int ,double*);
 double MMG2D_lencurv_iso(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,MMG5_int );
 double MMG2D_lencurv_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,MMG5_int );
 int MMG2D_chkedg(MMG5_pMesh ,MMG5_int );
@@ -287,8 +287,8 @@ int MMG2D_swapar(MMG5_pMesh ,MMG5_int ,int8_t );
 int MMG5_interpmet22(MMG5_pMesh ,double *,double *,double ,double *);
 int MMG2D_intmet_iso(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t ,MMG5_int ,double );
 int MMG2D_intmet_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t ,MMG5_int ,double );
-MMG5_int MMG2D_adpspl(MMG5_pMesh ,MMG5_pSol );
-MMG5_int MMG2D_movtri(MMG5_pMesh ,MMG5_pSol ,int ,int8_t );
+MMG5_int MMG2D_adpspl(MMG5_pMesh ,MMG5_pSol ,double*);
+MMG5_int MMG2D_movtri(MMG5_pMesh ,MMG5_pSol ,int ,int8_t ,double*);
 MMG5_int MMG2D_chkspl(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t );
 int MMG2D_split1b(MMG5_pMesh ,MMG5_int ,int8_t ,MMG5_int );
 int MMG2D_movedgpt(MMG5_pMesh ,MMG5_pSol ,int ,MMG5_int *,int8_t );
diff --git a/src/mmg2d/libmmg2d_tools.c b/src/mmg2d/libmmg2d_tools.c
index 0ebf7e9ab..c48b53af4 100644
--- a/src/mmg2d/libmmg2d_tools.c
+++ b/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],"-bdy-adaptation") ) {
+          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],"-hmin_factor") ) {
+            param = MMG2D_DPARAM_hmin_factor;
+          }
+          else if ( !strcmp(argv[i],"-hmax_factor") ) {
+            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],"-isotropic") ) {
+          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_isotropic,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],"-limit_angle") && ++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") ) {
diff --git a/src/mmg2d/libmmg2df.c b/src/mmg2d/libmmg2df.c
index 955ccf4c7..953eb7325 100644
--- a/src/mmg2d/libmmg2df.c
+++ b/src/mmg2d/libmmg2df.c
@@ -51,7 +51,7 @@ FORTRAN_NAME(MMG2D_MMG2DLIB,mmg2d_mmg2dlib,(MMG5_pMesh *mesh,MMG5_pSol *met
                                             ,int* retval),(mesh,met
                                                            ,retval)){
 
-  *retval = MMG2D_mmg2dlib(*mesh,*met);
+  *retval = MMG2D_mmg2dlib(*mesh,*met,NULL);
 
   return;
 }
diff --git a/src/mmg2d/mmg2d.c b/src/mmg2d/mmg2d.c
index 65f8c7d89..b4b5e1d32 100644
--- a/src/mmg2d/mmg2d.c
+++ b/src/mmg2d/mmg2d.c
@@ -275,6 +275,7 @@ int MMG2D_defaultOption(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol sol) {
 int main(int argc,char *argv[]) {
   MMG5_pMesh    mesh;
   MMG5_pSol     sol,met,disp,ls;
+  double*       velocity;
   int           ier,ierSave,fmtin,fmtout;
   char          stim[32],*ptr;
 
@@ -469,7 +470,7 @@ int main(int argc,char *argv[]) {
       MMG2D_RETURN_AND_FREE(mesh,met,ls,disp,MMG5_STRONGFAILURE);
     }
 
-    ier = MMG2D_mmg2dlib(mesh,met);
+    ier = MMG2D_mmg2dlib(mesh,met,velocity);
   }
 
   if ( ier != MMG5_STRONGFAILURE ) {
diff --git a/src/mmg2d/mmg2d1.c b/src/mmg2d/mmg2d1.c
index 4b4f7704e..d0e957ac7 100644
--- a/src/mmg2d/mmg2d1.c
+++ b/src/mmg2d/mmg2d1.c
@@ -37,7 +37,7 @@
  based on patterns, collapses and swaps.
    typchk = 1 -> adaptation based on edge lengths
    typchk = 2 -> adaptation based on lengths calculated in metric met */
-int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
+int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk,double* velocity) {
   int       it,maxit;
   MMG5_int  ns,nc,nsw,nns,nnc,nnsw;
 
@@ -50,13 +50,13 @@ 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 >= 4.*atan(1.)) {
       /* Memory free */
       MMG5_DEL_MEM(mesh,mesh->adja);
       mesh->adja = 0;
 
       /* Split long edges according to patterns */
-      ns = MMG2D_anaelt(mesh,met,typchk);
+      ns = MMG2D_anaelt(mesh,met,typchk,velocity);
       if ( ns < 0 ) {
         fprintf(stderr,"  ## Unable to complete surface mesh. Exit program.\n");
         return 0;
@@ -69,7 +69,7 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
       }
 
       /* Collapse short edges */
-      nc = MMG2D_colelt(mesh,met,typchk,MMG2D_LSHRT);
+      nc = MMG2D_colelt(mesh,met,typchk,MMG2D_LSHRT,velocity);
       if ( nc < 0 ) {
         fprintf(stderr,"  ## Unable to collapse mesh. Exiting.\n");
         return 0;
@@ -82,7 +82,7 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
     }
     /* Swap edges */
     if ( !mesh->info.noswap ) {
-      nsw = MMG2D_swpmsh(mesh,met,typchk);
+      nsw = MMG2D_swpmsh(mesh,met,typchk,velocity);
       if ( nsw < 0 ) {
         fprintf(stderr,"  ## Unable to improve mesh. Exiting.\n");
         return 0;
@@ -107,8 +107,134 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
   return 1;
 }
 
+// Minimum angle function
+double MMG2D_minangle(MMG5_pMesh mesh, int k, double* velocity) {
+
+  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.max[2]; // 1.75;
+  double factor_min = mesh->info.min[2]; // 1.33;
+
+  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.;
+
+  // Velocity is only considered for already existing points
+  if (velocity != NULL && 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 {
+    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;
+  }
+
+  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.);
+
+  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.);
+
+  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 tolerance is multiplied by two
+  if (mesh->info.bdy_adaptation) {
+    if ((pp1->tag & MG_BDY) || (pp2->tag & MG_BDY) || (pp3->tag & MG_BDY)) {
+      coef = c1;
+    }
+    else {
+      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;
+      }
+    }
+  }
+
+  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_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) {
   MMG5_pTria      pt;
   MMG5_pPoint     ppt,p1,p2;
   MMG5_Hash       hash;
@@ -130,6 +256,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 <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue;
+
     /* Check if pt should be cut */
     pt->flag = 0;
 
@@ -434,7 +562,7 @@ int MMG2D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) {
 }
 
 /* Travel triangles and collapse short edges */
-MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax) {
+MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax,double* velocity) {
   MMG5_pTria   pt;
   MMG5_pPoint  p1,p2;
   double       ll;
@@ -450,6 +578,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 <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > 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
 
@@ -512,7 +642,7 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax) {
 }
 
 /* Travel triangles and swap edges to improve quality */
-MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
+MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) {
   MMG5_pTria pt;
   int        it,maxit;
   MMG5_int   ns,nns;
@@ -523,12 +653,16 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) {
   maxit = 2;
   mesh->base++;
 
+  if (mesh->info.limit_angle <= 4*atan(1.)) 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 <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > 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,13 +676,16 @@ 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);
+
+  if (mesh->info.limit_angle <= 4*atan(1.)) mesh->info.fem = mesh->info.setfem;
+
   return nns;
 }
 
 
 /* Mesh adaptation routine for the final stage of the algorithm: intertwine splitting
  based on patterns, collapses, swaps and vertex relocations.*/
-int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
+int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) {
   int                  maxit,it;
   MMG5_int             nns,ns,nnc,nc,nnsw,nsw,nnm,nm;
 
@@ -558,14 +695,14 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
   do {
 
     if ( !mesh->info.noinsert ) {
-      ns = MMG2D_adpspl(mesh,met);
+      ns = MMG2D_adpspl(mesh,met,velocity);
       if ( ns < 0 ) {
         fprintf(stderr,"  ## Problem in function adpspl."
                 " Unable to complete mesh. Exit program.\n");
         return 0;
       }
 
-      nc = MMG2D_colelt(mesh,met,2,MMG2D_LOPTS);
+      nc = MMG2D_colelt(mesh,met,2,MMG2D_LOPTS,velocity);
       if ( nc < 0 ) {
         fprintf(stderr,"  ## Problem in function colelt."
                 " Unable to complete mesh. Exit program.\n");
@@ -578,7 +715,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
     }
     
     if ( !mesh->info.noswap ) {
-      nsw = MMG2D_swpmsh(mesh,met,2);
+      nsw = MMG2D_swpmsh(mesh,met,2,velocity);
       if ( nsw < 0 ) {
         fprintf(stderr,"  ## Problem in function swpmsh."
                 " Unable to complete mesh. Exit program.\n");
@@ -589,7 +726,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
       nsw = 0;
 
     if ( !mesh->info.nomove ) {
-      nm = MMG2D_movtri(mesh,met,1,0);
+      nm = MMG2D_movtri(mesh,met,1,0,velocity);
       if ( nm < 0 ) {
         fprintf(stderr,"  ## Problem in function movtri. "
                 "Unable to complete mesh. Exit program.\n");
@@ -613,7 +750,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
 
   /* Last iterations of vertex relocation only */
   if ( !mesh->info.nomove ) {
-    nm = MMG2D_movtri(mesh,met,5,1);
+    nm = MMG2D_movtri(mesh,met,5,1,velocity);
     if ( nm < 0 ) {
       fprintf(stderr,"  ## Problem in function movtri. Unable to complete mesh."
               " Exit program.\n");
@@ -638,7 +775,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) {
  * edges are only splitted on a one-by-one basis
  *
  */
-MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) {
+MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) {
   MMG5_pTria         pt;
   double             lmax,len;
   MMG5_int           ip,ns,k,nt;
@@ -653,6 +790,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 <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue;
+
     imax = -1;
     lmax = -1.0;
     for (i=0; i<3; i++) {
@@ -694,7 +833,7 @@ MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) {
 }
 
 /* Analyze points to relocate them according to a quality criterion */
-MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
+MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve,double* velocity) {
   MMG5_pTria           pt;
   MMG5_pPoint          p0;
   MMG5_int             nnm,nm,ns,k;
@@ -715,6 +854,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 <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > 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;
@@ -727,7 +868,7 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
           if ( ier ) ns++;
         }
         else {
-          if ( met->size == 3 && met->m )
+          if ( met->size == 3 && met->m && !mesh->info.isotropic)
             ier = MMG2D_movintpt_ani(mesh,met,ilist,list,improve);
           else
             ier = MMG2D_movintpt(mesh,met,ilist,list,improve);
@@ -758,13 +899,16 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) {
  * Mesh adaptation -- new version of mmg2d1.c
  *
  **/
-int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) {
+int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) {
+
+  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");
 
-  if ( !MMG2D_anatri(mesh,met,1) ) {
+  if ( !MMG2D_anatri(mesh,met,1,velocity) ) {
     fprintf(stderr,"  ## Unable to split mesh-> Exiting.\n");
     return 0;
   }
@@ -813,7 +957,7 @@ int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) {
     return 1;
   }
 
-  if ( !MMG2D_anatri(mesh,met,2) ) {
+  if ( !MMG2D_anatri(mesh,met,2,velocity) ) {
     fprintf(stderr,"  ## Unable to proceed adaptation. Exit program.\n");
     return 0;
   }
@@ -827,7 +971,7 @@ int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) {
 
 
   /* Stage 3: fine mesh improvements */
-  if ( !MMG2D_adptri(mesh,met) ) {
+  if ( !MMG2D_adptri(mesh,met,velocity) ) {
     fprintf(stderr,"  ## Unable to make fine improvements. Exit program.\n");
     return 0;
   }