diff --git a/cmake/testing/mmg3d_tests.cmake b/cmake/testing/mmg3d_tests.cmake
index 0c64adeff..ccd81ec34 100644
--- a/cmake/testing/mmg3d_tests.cmake
+++ b/cmake/testing/mmg3d_tests.cmake
@@ -853,7 +853,7 @@ endif()
 
 ADD_TEST(NAME test_para_tria
   COMMAND ${EXECUT_MMG3D}
-  -ar 0.02 -nofem -nosizreq -hgradreq -1 -hgrad -1
+  -ar 0.02 -fem 0 -nosizreq -hgradreq -1 -hgrad -1
   ${MMG3D_CI_TESTS}/test_para_tria/proc0.mesh
   -sol ${MMG3D_CI_TESTS}/test_para_tria/proc0.sol
   ${CTEST_OUTPUT_DIR}/proc0.o.mesh
diff --git a/src/mmg2d/API_functions_2d.c b/src/mmg2d/API_functions_2d.c
index 3205094b3..18ac9d346 100644
--- a/src/mmg2d/API_functions_2d.c
+++ b/src/mmg2d/API_functions_2d.c
@@ -147,8 +147,8 @@ int MMG2D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam, MMG5_int va
       mesh->info.dhd    = MMG5_ANGEDG;
     }
     break;
-  case MMG2D_IPARAM_nofem :
-    mesh->info.setfem = (val==1)? 0 : 1;
+  case MMG2D_IPARAM_fem :
+    mesh->info.setfem = val;
     break;
   case MMG2D_IPARAM_opnbdy :
     mesh->info.opnbdy = val;
diff --git a/src/mmg2d/libmmg2d.h b/src/mmg2d/libmmg2d.h
index f0097d566..93473d256 100644
--- a/src/mmg2d/libmmg2d.h
+++ b/src/mmg2d/libmmg2d.h
@@ -141,7 +141,7 @@ extern "C" {
     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_fem,               /*!< [0/1], make the mesh suitable for finite-element computations */
     MMG2D_IPARAM_isoref,            /*!< [0/n], Iso-surface boundary material reference */
   };
 
diff --git a/src/mmg2d/libmmg2d_tools.c b/src/mmg2d/libmmg2d_tools.c
index 0ebf7e9ab..84da804ed 100644
--- a/src/mmg2d/libmmg2d_tools.c
+++ b/src/mmg2d/libmmg2d_tools.c
@@ -83,7 +83,8 @@ int MMG2D_usage(char *name) {
   fprintf(stdout,"-3dMedit val read and write for gmsh visu: output only if val=1, input and output if val=2, input if val=3\n");
   fprintf(stdout,"\n");
 
-  fprintf(stdout,"-nofem       do not force Mmg to create a finite element mesh \n");
+  fprintf(stdout,"-fem     n   n = 0 : do not force mmg to create a finite element mesh \n");
+  fprintf(stdout,"             n = 1 : create a finite element mesh (default)\n");
   fprintf(stdout,"-nosurf      no surface modifications\n");
 
   /* Common parameters (second section) */
@@ -168,6 +169,17 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
             return 0;
           }
         }
+        else if ( !strcmp(argv[i],"-fem") ) {
+          if ( ++i < argc && isdigit(argv[i][0]) && (atoi(argv[i]) == 0 || atoi(argv[i]) == 1) ) {
+            if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_fem,atoi(argv[i])) )
+              return 0;
+          }
+          else {
+            fprintf(stderr,"\nMissing or unexpected argument option %s\n",argv[i-1]);
+            MMG2D_usage(argv[0]);
+            return 0;
+          }
+        }
         break;
       case 'h':
         param = MMG5_UNSET;
@@ -309,10 +321,6 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
         }
         break;
       case 'n':
-        if ( !strcmp(argv[i],"-nofem") ) {
-          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_nofem,1) )
-            return 0;
-        }
         if ( !strcmp(argv[i],"-nreg") ) {
           if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_nreg,1) )
             return 0;
diff --git a/src/mmg3d/API_functions_3d.c b/src/mmg3d/API_functions_3d.c
index d86290777..c94c86a5f 100644
--- a/src/mmg3d/API_functions_3d.c
+++ b/src/mmg3d/API_functions_3d.c
@@ -2153,8 +2153,8 @@ int MMG3D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam,MMG5_int val
       mesh->info.dhd    = MMG5_ANGEDG;
     }
     break;
-  case MMG3D_IPARAM_nofem :
-    mesh->info.setfem = (val==1)? 0 : 1;
+  case MMG3D_IPARAM_fem :
+    mesh->info.setfem = val;
     break;
   case MMG3D_IPARAM_opnbdy :
     mesh->info.opnbdy = val;
diff --git a/src/mmg3d/colver_3d.c b/src/mmg3d/colver_3d.c
index 30aed7edd..8d71067c1 100644
--- a/src/mmg3d/colver_3d.c
+++ b/src/mmg3d/colver_3d.c
@@ -76,7 +76,7 @@ int MMG5_chkcol_int(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t iface,
     /* Update edges tag for pt0 (needed by lenedg). */
 
     /* prevent from recreating internal edge between boundaries */
-    if ( mesh->info.fem==typchk ) {
+    if ( mesh->info.fem ) {
       p0 = &mesh->point[nq];
       if ( (p0->tag & MG_BDY) && !(p0->tag & MG_PARBDY) ) {
         i = ip;
diff --git a/src/mmg3d/libmmg3d.c b/src/mmg3d/libmmg3d.c
index aab14438e..62bfc6679 100644
--- a/src/mmg3d/libmmg3d.c
+++ b/src/mmg3d/libmmg3d.c
@@ -1073,8 +1073,8 @@ int MMG3D_mmg3dlib(MMG5_pMesh mesh,MMG5_pSol met) {
     fprintf(stdout,"\n  -- PHASE 1 : ANALYSIS\n");
   }
 
-  /* reset fem value to user setting (needed for multiple library call) */
-  mesh->info.fem = mesh->info.setfem;
+  /* reset fem value to 0 (needed for multiple library call) */
+  mesh->info.fem = 0;
 
   /* scaling mesh */
   if ( !MMG5_scaleMesh(mesh,met,NULL) )   _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
@@ -1312,8 +1312,8 @@ int MMG3D_mmg3dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
     fprintf(stdout,"\n  -- PHASE 1 : ISOSURFACE DISCRETIZATION\n");
   }
 
-  /* reset fem value to user setting (needed for multiple library call) */
-  mesh->info.fem = mesh->info.setfem;
+  /* reset fem value to 0 (needed for multiple library call) */
+  mesh->info.fem = 0;
 
   /* scaling mesh */
   if ( !MMG5_scaleMesh(mesh,met,sol) ) {
@@ -1578,7 +1578,7 @@ int MMG3D_mmg3dmov(MMG5_pMesh mesh,MMG5_pSol met, MMG5_pSol disp) {
   }
 
   /* reset fem value to user setting (needed for multiple library call) */
-  mesh->info.fem = mesh->info.setfem;
+  mesh->info.fem = 0;
 
   /* scaling mesh */
   if ( !MMG5_scaleMesh(mesh,met,disp) )   _LIBMMG5_RETURN(mesh,met,disp,MMG5_STRONGFAILURE);
diff --git a/src/mmg3d/libmmg3d.h b/src/mmg3d/libmmg3d.h
index 26d3c265c..6641b4c5e 100644
--- a/src/mmg3d/libmmg3d.h
+++ b/src/mmg3d/libmmg3d.h
@@ -147,7 +147,7 @@ enum MMG3D_Param {
   MMG3D_IPARAM_angle,                     /*!< [1/0], Turn on/off angle detection */
   MMG3D_IPARAM_iso,                       /*!< [1/0], Enable level-set discretization (volume and surfaces) */
   MMG3D_IPARAM_isosurf,                   /*!< [1/0], Enable level-set discretization on the surfaces only */
-  MMG3D_IPARAM_nofem,                     /*!< [1/0], Do not attempt to make the mesh suitable for finite-element computations */
+  MMG3D_IPARAM_fem,                       /*!< [0/1/2], make the mesh suitable for finite-element computations with two different levels of requirements */
   MMG3D_IPARAM_opnbdy,                    /*!< [1/0], Preserve triangles at interface of 2 domains with the same reference */
   MMG3D_IPARAM_lag,                       /*!< [-1/0/1/2], Enable Lagrangian motion */
   MMG3D_IPARAM_optim,                     /*!< [1/0], Optimize mesh keeping its initial edge sizes */
diff --git a/src/mmg3d/libmmg3d_private.h b/src/mmg3d/libmmg3d_private.h
index 19d1fc75a..13ea75519 100644
--- a/src/mmg3d/libmmg3d_private.h
+++ b/src/mmg3d/libmmg3d_private.h
@@ -149,6 +149,7 @@ extern "C" {
 #define MMG3D_BADKAL    0.2
 #define MMG3D_MAXKAL     1.
 
+#define MMG3D_STRONGFEM 2
 
 #define MMG3D_NPMAX  1000000 //200000
 #define MMG3D_NAMAX   200000 //40000
diff --git a/src/mmg3d/libmmg3d_tools.c b/src/mmg3d/libmmg3d_tools.c
index 43fa2b7d0..5b1c385ad 100644
--- a/src/mmg3d/libmmg3d_tools.c
+++ b/src/mmg3d/libmmg3d_tools.c
@@ -184,7 +184,9 @@ int MMG3D_usage(char *prog) {
 #endif
   fprintf(stdout,"\n");
 
-  fprintf(stdout,"-nofem       do not force Mmg to create a finite element mesh \n");
+  fprintf(stdout,"-fem     n   n = 0 : do not force mmg to create a finite element mesh \n");
+  fprintf(stdout,"             n = 1 : create a finite element mesh (default)\n");
+  fprintf(stdout,"             n = 2 : forbid regular edges joining feature points\n");
   fprintf(stdout,"-nosurf      no surface modifications\n");
 
   fprintf(stdout,"\n");
@@ -308,6 +310,17 @@ int MMG3D_storeknownar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,
             return 0;
           }
         }
+        else if ( !strcmp(argv[i],"-fem") ) {
+          if ( ++i < argc && isdigit(argv[i][0]) ) {
+            if ( !MMG3D_Set_iparameter(mesh,met,MMG3D_IPARAM_fem,atoi(argv[i])) )
+              return 0;
+          }
+          else {
+            fprintf(stderr,"\nMissing or unexpected argument option %s\n",argv[i-1]);
+            MMG3D_usage(argv[0]);
+            return 0;
+          }
+        }
         else {
           /* Arg unknown by Mmg: arg starts with -f but is not known */
           MMG_ARGV_APPEND(argv, mmgArgv, i, *mmgArgc,return 0);
@@ -458,11 +471,7 @@ int MMG3D_storeknownar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,
         }
         break;
       case 'n':
-        if ( !strcmp(argv[i],"-nofem") ) {
-          if ( !MMG3D_Set_iparameter(mesh,met,MMG3D_IPARAM_nofem,1) )
-            return 0;
-        }
-        else if ( !strcmp(argv[i],"-nreg") ) {
+        if ( !strcmp(argv[i],"-nreg") ) {
           if ( !MMG3D_Set_iparameter(mesh,met,MMG3D_IPARAM_nreg,1) )
             return 0;
         }
diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c
index dc005e3c7..d930bee33 100644
--- a/src/mmg3d/mmg3d1.c
+++ b/src/mmg3d/mmg3d1.c
@@ -487,9 +487,17 @@ int8_t MMG5_chkedg(MMG5_pMesh mesh,MMG5_Tria *pt,int8_t ori, double hmax,
     /*   } */
     /* } */
 
-    hma2 = MMG3D_LLONG*MMG3D_LLONG*hmax*hmax;
+    if ( mesh->info.setfem == MMG3D_STRONGFEM) {
+      /* Split regular boundary edges connecting ridge or non-manifold points */
+      if ( (MG_GEO_OR_NOM(p[i1]->tag) && MG_GEO_OR_NOM(p[i2]->tag)) && !MG_GEO_OR_NOM(pt->tag[i]) ) {
+        MG_SET(pt->flag,i);
+        continue;
+      }
+    }
 
     /* check length */
+    hma2 = MMG3D_LLONG*MMG3D_LLONG*hmax*hmax;
+
     ux = p[i2]->c[0] - p[i1]->c[0];
     uy = p[i2]->c[1] - p[i1]->c[1];
     uz = p[i2]->c[2] - p[i1]->c[2];
@@ -2282,7 +2290,7 @@ static MMG5_int MMG3D_anatets_ani(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
  * \remark ridge points creation: with fem (finite element method) mode a tetra
  * cannot have 2 boundary faces. Thus, the ridge point is created from a given
  * tetra and it is seen a second time from another tetra, which allows to update
- * its second normal. With nofem mode, a ref edge or ridge can be at the
+ * its second normal. With fem 0 mode, a ref edge or ridge can be at the
  * interface of 2 boundary faces belonging to the same tetra.
  *
  */
@@ -2428,7 +2436,7 @@ MMG3D_anatets_iso(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
 
           memcpy(ppt->n,to,3*sizeof(double));
 
-          if ( mesh->info.fem<typchk ) {
+          if ( !mesh->info.fem ) {
             /* A ridge can be at the interface of 2 boundary faces of the same
              * tetra: second normal has to be computed */
             if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
@@ -3176,13 +3184,22 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
   minit = 3;
   maxit = 6;
   mesh->gap = 0.5;
-  do {
-    if ( typchk==2 && lastit==1 )  ++mesh->info.fem;
 
+  if ( mesh->info.setfem && typchk == 1 ) {
+    mesh->info.fem = 1;
+  }
+  else if ( mesh->info.setfem && typchk == 2) {
+    mesh->info.fem = 0;
+  }
+
+  do {
+    if ( mesh->info.setfem && typchk == 2 && lastit == 1 ) {
+      mesh->info.fem = 1;
+    }
 
     /* split or swap tetra with more than 2 bdry faces */
     nf = ier = 0;
-    if ( mesh->info.fem == typchk ) {
+    if ( mesh->info.fem ) {
       ier = MMG5_anatet4(mesh,met,&nf,typchk);
       if ( ier < 0 )  return 0;
     }
@@ -3272,9 +3289,6 @@ int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
       ++lastit;
     }
     else if ( lastit ) {
-      /* Avoid the incrementation of mesh->info.fem if we have detected a last
-         iteration but anatet4 leads to have nc, nf or ns != 0 so we perform a last
-         iter */
       ++lastit;
     }
   }
diff --git a/src/mmg3d/mmg3d1_delone.c b/src/mmg3d/mmg3d1_delone.c
index 62911114b..de831afde 100644
--- a/src/mmg3d/mmg3d1_delone.c
+++ b/src/mmg3d/mmg3d1_delone.c
@@ -203,9 +203,8 @@ int MMG3D_mmg3d1_delone_split(MMG5_pMesh mesh, MMG5_pSol met,
      * only non bdy edge are considered */
 
     int8_t force_splt = 0;
-    const int8_t fem_mode = 2; // value of info.fem in case of fem mode
 
-    if ( mesh->info.fem == fem_mode ) {
+    if ( mesh->info.fem ) {
       /* Force splitting of internal edges connecting bdy points */
       if ( MG_TRUE_BDY(p0->tag) && MG_TRUE_BDY(p1->tag) ) {
         force_splt = 1;
diff --git a/src/mmg3d/swap_3d.c b/src/mmg3d/swap_3d.c
index 92cf12872..81db2cb5b 100644
--- a/src/mmg3d/swap_3d.c
+++ b/src/mmg3d/swap_3d.c
@@ -137,6 +137,20 @@ int MMG5_chkswpbdy(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list,int ilist,
            (tt2.v[MMG5_inxt2[ia2]] == nq && tt2.v[MMG5_iprv2[ia2]] == np) );
   na2 = tt2.v[ia2];
 
+  if ( mesh->info.setfem == MMG3D_STRONGFEM ) {
+#ifndef NDEBUG
+    /* Check that we don't have a regular bdy edge connecting ridge or nm points */
+    if ( MG_GEO_OR_NOM(mesh->point[np].tag) && MG_GEO_OR_NOM(mesh->point[nq].tag) ) {
+      assert ( ! MG_GEO_OR_NOM(tt1.tag[ia1]) );
+      assert ( ! MG_GEO_OR_NOM(tt2.tag[ia2]) );
+    }
+#endif
+    /* No swap if it creates a regular boundary edge connecting ridge or non-manifold points */
+    if ( MG_GEO_OR_NOM(mesh->point[na1].tag) && MG_GEO_OR_NOM(mesh->point[na2].tag) ) {
+      return 0;
+    }
+  }
+
   /* Check non convexity (temporarily use b0,b1)*/
   MMG5_norpts(mesh,tt1.v[ia1],tt1.v[MMG5_inxt2[ia1]],tt2.v[ia2],b0);
   MMG5_norpts(mesh,tt2.v[ia2],tt2.v[MMG5_inxt2[ia2]],tt1.v[ia1],b1);