diff --git a/Common/Context.h b/Common/Context.h
index 2eaa2fffb6ecd56655546673ec31bf9791b67b2f..de5eeb1826b4ce15c86b947500294242c0f04b58 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -49,6 +49,7 @@ struct contextMeshOptions {
   int nLayersPerGap;
   double gradation;
   int quadqsSizemapMethod, quadqsTopoOptimMethods;
+  double quadqsRemeshingBoldness;
   // mesh IO
   int fileFormat, firstElementTag, firstNodeTag;
   double mshFileVersion, medFileMinorVersion, scalingFactor;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 0e2c2dfbd516a2bdd75ff9d8c42f2d85ffe00c7f..a32e3af8ce875c2a52cb446015df6a893fee0b0a 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1386,7 +1386,8 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "QuadqsSizemapMethod" , opt_mesh_quadqs_sizemap_method, 0. ,
     "Size map method in QuadQuasiStructured. 0: default, 1: cross-field,"
       "2: cross-field + CAD small features adaptation,"
-      "3: from background mesh (e.g. sizes in current triangulation)"
+      "3: from background mesh (e.g. sizes in current triangulation),"
+      "4: cross-field + CAD small features adaptation (clamped by background mesh)"
   },
   { F|O, "QuadqsTopologyOptimizationMethods" , opt_mesh_quadqs_topo_optim_methods, 0. ,
     "Topology optimization methods in QuadQuasiStructured. 0: default (all),"
@@ -1395,6 +1396,11 @@ StringXNumber MeshOptions_Number[] = {
       "001: cavity remeshing,"
       "xxx: combination of multiple methods (e.g. 111 for all)"
   },
+  { F|O, "QuadqsRemeshingBoldness" , opt_mesh_quadqs_remeshing_boldness, 0.501 ,
+    "Controls how much cavity remeshing is allowed to distort"
+      " the quad mesh. From 0 (no quality decrease during remeshing) to 1"
+      " (quality can tend to 0 during remeshing)."
+  },
 
   { F|O, "Quadrangles" , opt_mesh_quadrangles , 1. ,
     "Display mesh quadrangles?" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index b0ba81367af14209baff6e2034ae340bd148bce6..f32a4fd71bea73ea6a62ad5e30f9a04571c4b5b4 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -6580,6 +6580,12 @@ double opt_mesh_quadqs_sizemap_method(OPT_ARGS_NUM)
   return CTX::instance()->mesh.quadqsSizemapMethod;
 }
 
+double opt_mesh_quadqs_remeshing_boldness(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET) CTX::instance()->mesh.quadqsRemeshingBoldness = (double)val;
+  return CTX::instance()->mesh.quadqsRemeshingBoldness;
+}
+
 double opt_mesh_quadqs_topo_optim_methods(OPT_ARGS_NUM) {
   if(action & GMSH_SET) CTX::instance()->mesh.quadqsTopoOptimMethods = (int)val;
   return CTX::instance()->mesh.quadqsTopoOptimMethods;
diff --git a/Common/Options.h b/Common/Options.h
index 59c0feb7df70309f0092e7f61c86e5d668dd840a..aac3ed6b278eb9fcc936ac74e751e74694043137 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -606,6 +606,7 @@ double opt_mesh_reparam_max_triangles(OPT_ARGS_NUM);
 double opt_mesh_ignore_parametrization(OPT_ARGS_NUM);
 double opt_mesh_quadqs_sizemap_method(OPT_ARGS_NUM);
 double opt_mesh_quadqs_topo_optim_methods(OPT_ARGS_NUM);
+double opt_mesh_quadqs_remeshing_boldness(OPT_ARGS_NUM);
 double opt_solver_listen(OPT_ARGS_NUM);
 double opt_solver_timeout(OPT_ARGS_NUM);
 double opt_solver_plugins(OPT_ARGS_NUM);
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index a1cee7104c863d396e557b450a7895a14e0891a7..1bade4c3b043b03d2430b9586cfdd9ef285af9b0 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -2361,6 +2361,14 @@ static void mesh_optimize_cb(Fl_Widget *w, void *data)
   drawContext::global()->draw();
 }
 
+static void mesh_optimize_quad_topo_cb(Fl_Widget *w, void *data)
+{
+  CTX::instance()->lock = 1;
+  GModel::current()->optimizeMesh("QuadQuasiStructured");
+  CTX::instance()->lock = 0;
+  drawContext::global()->draw();
+}
+
 static void mesh_cross_compute_cb(Fl_Widget *w, void *data)
 {
   std::vector<int> tags;
@@ -4649,6 +4657,8 @@ static menuItem static_modules[] = {
   {"0Modules/Mesh/Experimental/Convert old partitioning",
    (Fl_Callback *)mesh_convert_old_partitioning_cb},
 #endif
+  {"0Modules/Mesh/Experimental/Optimize quad topology",
+   (Fl_Callback *)mesh_optimize_quad_topo_cb},
   {"0Modules/Mesh/Reverse/Elements", (Fl_Callback *)mesh_reverse_parts_cb,
    (void *)"elements"},
   {"0Modules/Mesh/Reverse/Curves", (Fl_Callback *)mesh_reverse_parts_cb,
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 58e724951a149696ad6b70e9fd9986be794eb293..349358845229e1882d9eee5af74e1ed85d061446 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -200,7 +200,10 @@ bool discreteFace_rtree_callback(std::pair<MTriangle *, MTriangle *> *t,
 GPoint discreteFace::closestPoint(const SPoint3 &queryPoint, double maxDistance,
                                   SVector3 *normal) const
 {
-  if(_param.empty()) return GPoint();
+  if(_param.empty()) {
+    Msg::Debug("discreteFace %i: no param, no closestPoint", tag());
+    return GPoint();
+  }
 
   dfWrapper wrapper(queryPoint);
   do {
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 38ea85701d3d95db8b96f0f6287fa69e3e9fcb9d..2513719974d6cc4ce5c1bfbd8f507198ead8d6c0 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -553,6 +553,8 @@ static void Mesh2D(GModel *m)
   Msg::SetNumThreads(prevNumThreads);
 
   if(CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT) {
+    replaceBadQuadDominantMeshes(m);
+
     /* In the quasi-structured pipeline, the quad-dominant mesh
      * is subdivided into a full quad mesh */
     /* TODO: - a faster CAD projection approach (from uv)
@@ -561,6 +563,7 @@ static void Mesh2D(GModel *m)
     // RefineMesh(m, linear, true, false);
     RefineMeshWithBackgroundMeshProjection(m);
 
+
     OptimizeMesh(m, "QuadQuasiStructured");
   }
 
diff --git a/Mesh/meshQuadQuasiStructured.cpp b/Mesh/meshQuadQuasiStructured.cpp
index 761290bb8f5f1db841da968b6b5afe7290457555..2c5ee51f104912629b950f2773be2425f4ba2b0d 100644
--- a/Mesh/meshQuadQuasiStructured.cpp
+++ b/Mesh/meshQuadQuasiStructured.cpp
@@ -18,6 +18,7 @@
 #include "GModel.h"
 #include "meshGEdge.h"
 #include "meshGFace.h"
+#include "meshGFaceOptimize.h"
 #include "meshGRegion.h"
 #include "BackgroundMesh.h"
 #include "Generator.h"
@@ -68,10 +69,18 @@ template <typename Key, typename Hash = robin_hood::hash<Key>,
 using unordered_set =
   robin_hood::detail::Table<true, MaxLoadFactor100, Key, void, Hash, KeyEqual>;
 
+
+constexpr int SizeMapDefault    = 0;
+constexpr int SizeMapCrossField = 1;
+constexpr int SizeMapCrossFieldAndSmallCad = 2;
+constexpr int SizeMapBackgroundMesh = 3;
+constexpr int SizeMapCrossFieldAndBMeshOnCurves = 4;
+
 const std::string BMESH_NAME = "bmesh_quadqs";
 
 constexpr bool PARANO_QUALITY = false;
 constexpr bool PARANO_VALIDITY = false;
+constexpr bool DBG_EXPORT = false;
 
 
 /* scaling applied on integer values stored in view (background field),
@@ -85,7 +94,7 @@ int buildBackgroundField(
   const std::vector<std::array<double, 5> >& global_singularity_list,
   const std::string &viewName = "guiding_field")
 {
-  Msg::Debug("build background guiding field ...");
+  Msg::Info("Build background field (view 'guiding_field') ...");
   if(global_triangles.size() != global_triangle_directions.size()) {
     Msg::Error("build background field: incoherent sizes in input");
     return -1;
@@ -303,13 +312,15 @@ int fillSizemapFromScalarBackgroundField(
 
 std::string nameOfSizeMapMethod(int method) {
   if (method == 0) {
-    return "default (cross-field conformal scaling and CAD adaptation)";
+    return "default (" + nameOfSizeMapMethod(4) + ")";
   } else if (method == 1) {
     return "cross-field conformal scaling";
   } else if (method == 2) {
     return "cross-field conformal scaling and CAD adaptation";
   } else if (method == 3) {
     return "background mesh sizes";
+  } else if (method == 4) {
+    return "cross-field conformal scaling and CAD adaptation (clamped by background mesh)";
   }
   return "unknown";
 }
@@ -361,6 +372,42 @@ bool generateMeshWithSpecialParameters(GModel *gm,
   return true;
 }
 
+bool getGFaceBackgroundMeshLinesAndTriangles(
+    GlobalBackgroundMesh& bmesh, GFace* gf,
+    std::vector<MLine*>& lines, std::vector<MTriangle*>& triangles) {
+  lines.clear();
+  triangles.clear();
+
+  /* Collect pointers to background mesh elements */
+  std::vector<GEdge *> edges = face_edges(gf);
+  for(GEdge *ge : edges) {
+    auto it = bmesh.edgeBackgroundMeshes.find(ge);
+    if(it == bmesh.edgeBackgroundMeshes.end()) {
+      Msg::Warning("getGFaceBackgroundMeshLinesAndTriangles: GEdge %i not found "
+          "in background mesh",
+          ge->tag());
+      continue;
+    }
+    lines.reserve(lines.size() + it->second.lines.size());
+    for(size_t i = 0; i < it->second.lines.size(); ++i) {
+      lines.push_back(&(it->second.lines[i]));
+    }
+  }
+  auto it = bmesh.faceBackgroundMeshes.find(gf);
+  if(it == bmesh.faceBackgroundMeshes.end()) {
+    Msg::Warning("getGFaceBackgroundMeshLinesAndTriangles: GFace %i not found "
+        "in background mesh",
+        gf->tag());
+    return false;
+  }
+  triangles.reserve(triangles.size() + it->second.triangles.size());
+  for(size_t i = 0; i < it->second.triangles.size(); ++i) {
+    triangles.push_back(&(it->second.triangles[i]));
+  }
+
+  return true;
+}
+
 int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
                                        bool deleteGModelMesh, int N)
 {
@@ -389,7 +436,7 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
         return 0;
       }
       else if(field && field->numComponents() == 1) {
-        if(qqsSizemapMethod == 0) {
+        if(qqsSizemapMethod == SizeMapDefault) {
           Msg::Info("scalar background field exists, using it as size map");
           externalSizemap = true;
         }
@@ -495,38 +542,13 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
         continue;
 
       /* Compute a cross field on each face */
+
+      /* Get mesh elements for solver */
       std::vector<MLine *> lines;
       std::vector<MTriangle *> triangles;
-
-      /* Collect pointers to background mesh elements */
-      std::vector<GEdge *> edges = face_edges(gf);
-      for(GEdge *ge : edges) {
-        auto it = bmesh.edgeBackgroundMeshes.find(ge);
-        if(it == bmesh.edgeBackgroundMeshes.end()) {
-          Msg::Warning("BuildBackgroundMeshAndGuidingField: GEdge %i not found "
-                       "in background mesh",
-                       ge->tag());
-          continue;
-        }
-        lines.reserve(lines.size() + it->second.lines.size());
-        for(size_t i = 0; i < it->second.lines.size(); ++i) {
-          lines.push_back(&(it->second.lines[i]));
-        }
-      }
-      auto it = bmesh.faceBackgroundMeshes.find(gf);
-      if(it == bmesh.faceBackgroundMeshes.end()) {
-        Msg::Warning("BuildBackgroundMeshAndGuidingField: GFace %i not found "
-                     "in background mesh",
-                     gf->tag());
-        continue;
-      }
-      triangles.reserve(triangles.size() + it->second.triangles.size());
-      for(size_t i = 0; i < it->second.triangles.size(); ++i) {
-        triangles.push_back(&(it->second.triangles[i]));
-      }
-
-      if(triangles.size() == 0) {
-        Msg::Error("- Face %i: no triangles in background mesh", gf->tag());
+      bool oklt = getGFaceBackgroundMeshLinesAndTriangles(bmesh, gf, lines, triangles);
+      if (!oklt && triangles.size() == 0) {
+        Msg::Error("- Face %i: failed to get triangles from background mesh", gf->tag());
         continue;
       }
 
@@ -591,7 +613,7 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
         }
       }
       else if(qqsSizemapMethod ==
-              3) { /* Size map from background triangulation */
+              SizeMapBackgroundMesh) { /* Size map from background triangulation */
         int sts = fillSizemapFromTriangleSizes(triangles, localSizemap);
         if(sts != 0) {
           Msg::Warning("- Face %i: failed to fill size map from triangle sizes",
@@ -653,9 +675,9 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
         }
       }
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
+      #if defined(_OPENMP)
+      #pragma omp critical
+      #endif
       {
         append(global_triangles, triangles);
         append(global_triangle_directions, triangleDirections);
@@ -664,8 +686,10 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
           global_size_map.push_back({kv.first, kv.second});
         }
       }
+
     }
   }
+
   sort_unique(global_size_map);
 
   /* Warning: from now on, code is not optimized in terms of data structures
@@ -690,10 +714,21 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
 
   /* Minimal size on curves */
   std::unordered_map<MVertex *, double> cadMinimalSizeOnCurves;
-  if((qqsSizemapMethod == 0 && !externalSizemap) || qqsSizemapMethod == 2) {
-    int scad = computeMinimalSizeOnCurves(bmesh, cadMinimalSizeOnCurves);
-    if(scad != 0) {
-      Msg::Warning("failed to compute minimal size on CAD curves");
+  if (!externalSizemap) {
+    if (   qqsSizemapMethod == SizeMapDefault 
+        || qqsSizemapMethod == SizeMapCrossFieldAndBMeshOnCurves
+        || qqsSizemapMethod == SizeMapCrossFieldAndSmallCad)  {
+
+      bool clampMinWithTriEdges = false;
+      if (qqsSizemapMethod == SizeMapDefault 
+          || qqsSizemapMethod == SizeMapCrossFieldAndBMeshOnCurves)  {
+        clampMinWithTriEdges = true;
+      }
+
+      int scad = computeMinimalSizeOnCurves(bmesh, clampMinWithTriEdges, cadMinimalSizeOnCurves);
+      if(scad != 0) {
+        Msg::Warning("failed to compute minimal size on CAD curves");
+      }
     }
   }
 
@@ -714,6 +749,7 @@ int BuildBackgroundMeshAndGuidingField(GModel *gm, bool overwriteGModelMesh,
   /* One-way propagation of values */
   const double gradientMax =
     1.2; /* this param should be a global gmsh option */
+  Msg::Info("Sizemap smoothing (progression ratio: %.2f)", gradientMax);
   int sop = sizeMapOneWaySmoothing(global_triangles, sizeMap, gradientMax);
   if(sop != 0) { Msg::Warning("failed to compute one-way size map smoothing"); }
 
@@ -774,7 +810,7 @@ bool getSingularitiesFromBackgroundField(GFace* gf,
     if(guiding_field && guiding_field->numComponents() == 3) { field = guiding_field; }
   }
   if (field == nullptr) {
-    Msg::Error("get singularities: face %i, failed to get background field", gf->tag());
+    Msg::Debug("get singularities: face %i, failed to get background field", gf->tag());
     return false;
   }
 
@@ -804,8 +840,53 @@ bool getSingularitiesFromBackgroundField(GFace* gf,
   return true;
 }
 
+bool getSingularitiesFromNewCrossFieldComputation(
+    GlobalBackgroundMesh& bmesh,
+    GFace* gf,
+    std::vector<std::pair<SPoint3, int> >& singularities) {
+  const int N = 4;
+
+  std::vector<MLine *> lines;
+  std::vector<MTriangle *> triangles;
+  bool oklt = getGFaceBackgroundMeshLinesAndTriangles(bmesh, gf, lines, triangles);
+  if (!oklt && triangles.size() == 0) {
+    Msg::Error("- Face %i: failed to get triangles from background mesh", gf->tag());
+    return false;
+  }
+
+  /* Cross field */
+  std::vector<std::array<double, 3> > triEdgeTheta;
+  int nbDiffusionLevels = 4;
+  double thresholdNormConvergence = 1.e-2;
+  int nbBoundaryExtensionLayer = 1;
+  int verbosity = 0;
+  Msg::Info("- Face %i: compute cross field (%li triangles, %li B.C. "
+      "edges, %i diffusion levels) ...",
+      gf->tag(), triangles.size(), lines.size(),
+      nbDiffusionLevels);
+  int scf = computeCrossFieldWithHeatEquation(
+      N, triangles, lines, triEdgeTheta, nbDiffusionLevels,
+      thresholdNormConvergence, nbBoundaryExtensionLayer, verbosity);
+  if(scf != 0) {
+    Msg::Warning("- Face %i: failed to compute cross field", gf->tag());
+  }
+
+  /* Cross field singularities */
+  bool addSingularitiesAtAcuteCorners = true;
+  double thresholdInDeg = 30.;
+  int scsi = detectCrossFieldSingularities(N, triangles, triEdgeTheta,
+      addSingularitiesAtAcuteCorners,
+      thresholdInDeg, singularities);
+  if(scsi != 0) {
+    Msg::Warning("- Face %i: failed to compute cross field singularities",
+        gf->tag());
+  }
+
+  return true;
+}
+
 void computeBdrVertexIdealValence(const std::vector<MQuadrangle *> &quads,
-                                  unordered_map<MVertex *, double> &qValIdeal)
+    unordered_map<MVertex *, double> &qValIdeal)
 {
   qValIdeal.clear();
   for(MQuadrangle *f : quads) {
@@ -1554,21 +1635,20 @@ int improveInteriorValences(
 int RefineMeshWithBackgroundMeshProjection(GModel *gm)
 {
   const bool SHOW_INTERMEDIATE_VIEWS = (Msg::GetVerbosity() >= 99);
-
-  const bool SHOW_QUADTRI = SHOW_INTERMEDIATE_VIEWS;
-  if(SHOW_QUADTRI) {
+  if(SHOW_INTERMEDIATE_VIEWS) {
     std::vector<MElement *> elements;
     for(GFace *gf : model_faces(gm)) {
       append(elements,
              dynamic_cast_vector<MTriangle *, MElement *>(gf->triangles));
       append(elements,
              dynamic_cast_vector<MQuadrangle *, MElement *>(gf->quadrangles));
-      // showUVParametrization(gf,dynamic_cast_vector<MTriangle*,MElement*>(gf->triangles),"quadtri_uvs");
-      // showUVParametrization(gf,dynamic_cast_vector<MQuadrangle*,MElement*>(gf->quadrangles),"quadtri_uvs");
     }
     GeoLog::add(elements, "qqs_quadtri");
     GeoLog::flush();
   }
+  if (DBG_EXPORT) {
+    gm->writeMSH("qqs_init.msh", 4.1);
+  }
 
   Msg::Info(
     "Refine mesh (midpoint subdivision, with background projection) ...");
@@ -1576,84 +1656,59 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
   bool linear = true;
   RefineMesh(gm, linear, true, false);
 
-  if(Msg::GetVerbosity() >= 99) {
+  if(true || Msg::GetVerbosity() >= 99) {
     std::unordered_map<std::string, double> stats;
     appendQuadMeshStatistics(gm, stats, "MPS_");
     printStatistics(stats, "Quad mesh after subdivision, before projection:");
+    gm->writeMSH("qqs_subdiv_noproj.msh", 4.1);
   }
 
-  GlobalBackgroundMesh *bmesh = nullptr;
-  if(backgroudMeshExists(BMESH_NAME)) {
-    bmesh = &(getBackgroundMesh(BMESH_NAME));
-  }
-  else {
-    Msg::Warning("refine mesh with background projection: no background mesh, "
-                 "using CAD projection (slow)");
-  }
 
+  /* Convert vertex types:
+   * - MVertex on curve -> MEdgeVertex
+   * - MVertex on face  -> MFaceVertex */
   std::vector<GEdge *> edges = model_edges(gm);
-
+  std::vector<GFace *> faces = model_faces(gm);
   /* Build custom edgeToFaces because ge->faces() does not work
    * for embedded edges */
+  std::unordered_map<GEdge*,std::vector<MVertex*> > toProjectOnCurve;
+  std::unordered_map<GFace*,std::vector<MVertex*> > toProjectOnFace;
   std::unordered_map<GEdge*,std::unordered_set<GFace*> > edgeToFaces;
   for (GFace* gf: model_faces(gm)) {
     std::vector<GEdge*> fedges = face_edges(gf);
     for (GEdge* ge: fedges) edgeToFaces[ge].insert(gf);
   }
 
-#if defined(_OPENMP)
-#pragma omp parallel for schedule(dynamic)
-#endif
+  #if defined(_OPENMP)
+  #pragma omp parallel for schedule(dynamic)
+  #endif
   for(size_t e = 0; e < edges.size(); ++e) {
     GEdge *ge = edges[e];
     if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) continue;
     if(ge->lines.size() == 0 || ge->mesh_vertices.size() == 0) continue;
-
-    Msg::Debug("- Curve %i: project midpoints on curve", ge->tag());
     unordered_map<MVertex *, MVertex *> old2new_ge;
-    unordered_map<MVertex *, SPoint3> backupPositionLinear;
-    double tPrev = 0;
+    std::vector<MVertex*>& projList = toProjectOnCurve[ge];
+    projList.reserve(ge->mesh_vertices.size());
     for(size_t i = 0; i < ge->mesh_vertices.size(); ++i) {
       MVertex *v = ge->mesh_vertices[i];
       MEdgeVertex *mev = dynamic_cast<MEdgeVertex *>(v);
       if(mev == nullptr) {
-        double t = tPrev;
-        GPoint proj = ge->closestPoint(v->point(), t);
-        if(proj.succeeded()) {
-          /* Need to change the type of the MVertex */
-          MVertex *v2 =
-            new MEdgeVertex(proj.x(), proj.y(), proj.z(), ge, proj.u());
-          tPrev = proj.u();
-          ge->mesh_vertices[i] = v2;
-          old2new_ge[v] = v2;
-          backupPositionLinear[v2] = v->point();
-          delete v;
-        }
-        else {
-          Msg::Warning("- Edge %i, vertex %li: curve projection failed",
-                       ge->tag(), v->getNum());
-          MVertex *v2 = new MEdgeVertex(v->point().x(), v->point().y(),
-                                        v->point().z(), ge, tPrev);
-          ge->mesh_vertices[i] = v2;
-          old2new_ge[v] = v2;
-          backupPositionLinear[v2] = v->point();
-          delete v;
-        }
-        tPrev = t;
+        MVertex* v2 = new MEdgeVertex(v->point().x(), v->point().y(), v->point().z(), ge, 0.);
+        ge->mesh_vertices[i] = v2;
+        old2new_ge[v] = v2;
+        projList.push_back(v2);
+        delete v;
       }
     }
-
-    bool qualityOk = true;
-
-    /* Update Lines */
-    for(MLine *l : ge->lines)
+    /* Update lines */
+    for(MLine *l : ge->lines) {
       for(size_t lv = 0; lv < 2; ++lv) {
         MVertex *v = l->getVertex(lv);
         auto it = old2new_ge.find(v);
         if(it != old2new_ge.end()) { l->setVertex(lv, it->second); }
       }
-
-    /* Update adjacent faces */
+    }
+    /* Update elements in adjacent faces */
     for(GFace *gf : edgeToFaces[ge]) {
       for(size_t i = 0; i < gf->getNumMeshElements(); ++i) {
         MElement *e = gf->getMeshElement(i);
@@ -1663,36 +1718,12 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
           if(it != old2new_ge.end()) { e->setVertex(lv, it->second); }
         }
       }
-
-      /* Check quality of quads */
-      if(gf->quadrangles.size() > 0) {
-        double sicnMin, sicnAvg;
-        computeSICN(
-          dynamic_cast_vector<MQuadrangle *, MElement *>(gf->quadrangles),
-          sicnMin, sicnAvg);
-        if(sicnMin < 0.) {
-          Msg::Warning("- refine with projection: quality negative (%.3f) on "
-                       "face %i after curve %i projection, rollback",
-                       sicnMin, gf->tag(), ge->tag());
-          qualityOk = false;
-        }
-      }
-    }
-
-    if(!qualityOk) { /* restore positions on the curve */
-      for(auto &kv : backupPositionLinear) { kv.first->setXYZ(kv.second); }
-    }
-    if(PARANO_VALIDITY) {
-      errorAndAbortIfInvalidVertexInElements(
-        dynamic_cast_vector<MLine *, MElement *>(ge->lines),
-        "after curve proj");
     }
   }
 
-  std::vector<GFace *> faces = model_faces(gm);
-#if defined(_OPENMP)
-#pragma omp parallel for schedule(dynamic)
-#endif
+  #if defined(_OPENMP)
+  #pragma omp parallel for schedule(dynamic)
+  #endif
   for(size_t f = 0; f < faces.size(); ++f) {
     GFace *gf = faces[f];
     if(CTX::instance()->mesh.meshOnlyVisible && !gf->getVisibility()) continue;
@@ -1700,108 +1731,27 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
        gf->tag() != CTX::instance()->debugSurface)
       continue;
     if(gf->triangles.size() == 0 && gf->quadrangles.size() == 0) continue;
-
-    SurfaceProjector *sp = nullptr;
-    if(bmesh) {
-      auto it = bmesh->faceBackgroundMeshes.find(gf);
-      if(it == bmesh->faceBackgroundMeshes.end()) {
-        Msg::Error("background mesh not found for face %i", gf->tag());
-        continue;
-      }
-
-      /* Get pointers to triangles in the background mesh */
-      std::vector<MTriangle *> triangles(it->second.triangles.size());
-      for(size_t i = 0; i < it->second.triangles.size(); ++i) {
-        triangles[i] = &(it->second.triangles[i]);
-      }
-
-      sp = new SurfaceProjector();
-      bool oki = sp->initialize(gf, triangles);
-      if(!oki) {
-        Msg::Warning("failed to initialize surface projector");
-        delete sp;
-        sp = nullptr;
-      }
-    }
-
-    /* Project the vertices which have been introduced by the RefineMesh */
-    Msg::Debug("- Face %i: project midpoints on surface", gf->tag());
-    bool evalOnCAD = gf->haveParametrization();
-    bool projOnCad = false;
-    if(evalOnCAD && !haveNiceParametrization(gf)) {
-      /* Strong disortion in parametrization, use projection */
-      evalOnCAD = false;
-      projOnCad = true;
-    }
-
-    if(Msg::GetVerbosity() >= 999) {
-      /* In some models (S9 from MAMBO) and with some resolution (clscale 0.1),
-       * there is this issue of vertex on wrong face ...
-       * I guess this is because the CAD is garbage (duplicated corners,
-       * degenerate edges) */
-      Msg::Debug("- Face %i, verify vertices before ...", gf->tag());
-      for(size_t i = 0; i < gf->getNumMeshElements(); ++i) {
-        MElement *e = gf->getMeshElement(i);
-        for(size_t lv = 0; lv < e->getNumVertices(); ++lv) {
-          MVertex *v = e->getVertex(lv);
-          if(v->onWhat()->cast2Face() && v->onWhat() != gf) {
-            Msg::Error(
-              "- Face %i: element vertex %li is associated to entity (%i,%i)",
-              gf->tag(), v->getNum(), v->onWhat()->dim(), v->onWhat()->tag());
-            abort();
-          }
-        }
-      }
-    }
-
     unordered_map<MVertex *, MVertex *> old2new_gf;
-    unordered_map<MVertex *, SPoint3> backupPositionLinear;
-
+    std::vector<MVertex*>& projList = toProjectOnFace[gf];
+    projList.reserve(gf->mesh_vertices.size());
     for(size_t i = 0; i < gf->mesh_vertices.size(); ++i) {
       MVertex *v = gf->mesh_vertices[i];
       if(v->onWhat() != gf) {
         Msg::Error("- Face %i: vertex %li is associated to entity (%i,%i)",
-                   gf->tag(), v->getNum(), v->onWhat()->dim(),
-                   v->onWhat()->tag());
-        abort();
+            gf->tag(), v->getNum(), v->onWhat()->dim(),
+            v->onWhat()->tag());
+        continue;
       }
       MFaceVertex *mfv = dynamic_cast<MFaceVertex *>(v);
       if(mfv == nullptr) {
-        GPoint proj;
-        if(sp != nullptr) {
-          proj = sp->closestPoint(v->point().data(), evalOnCAD, projOnCad);
-
-          if(!proj.succeeded() && gf->haveParametrization()) {
-            double uvg[2] = {0., 0.};
-            proj = gf->closestPoint(v->point(), uvg);
-          }
-        }
-        else {
-          double uvg[2] = {0., 0.};
-          proj = gf->closestPoint(v->point(), uvg);
-        }
-        if(proj.succeeded()) {
-          /* Need to change the type of the MVertex */
-          MVertex *v2 = new MFaceVertex(proj.x(), proj.y(), proj.z(), gf,
-                                        proj.u(), proj.v());
-          gf->mesh_vertices[i] = v2;
-          backupPositionLinear[v2] = v->point();
-          old2new_gf[v] = v2;
-          delete v;
-        }
-        else {
-          MVertex *v2 = new MFaceVertex(v->point().x(), v->point().y(),
-                                        v->point().z(), gf, 0., 0.);
-          gf->mesh_vertices[i] = v2;
-          backupPositionLinear[v2] = v->point();
-          old2new_gf[v] = v2;
-          delete v;
-          Msg::Warning("- Face %i, vertex %li: surface projection failed",
-                       gf->tag(), v->getNum());
-        }
+        MVertex* v2 = new MFaceVertex(v->point().x(), v->point().y(),
+            v->point().z(), gf, 0., 0.);
+        gf->mesh_vertices[i] = v2;
+        old2new_gf[v] = v2;
+        projList.push_back(v2);
+        delete v;
       }
     }
-
     /* Update elements */
     for(size_t i = 0; i < gf->getNumMeshElements(); ++i) {
       MElement *e = gf->getMeshElement(i);
@@ -1811,73 +1761,123 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
         if(it2 != old2new_gf.end()) { e->setVertex(lv, it2->second); }
       }
     }
+  }
 
-    if(Msg::GetVerbosity() >= 99) {
-      Msg::Debug("- Face %i, verify vertices ...", gf->tag());
-      for(size_t i = 0; i < gf->getNumMeshElements(); ++i) {
-        MElement *e = gf->getMeshElement(i);
-        for(size_t lv = 0; lv < e->getNumVertices(); ++lv) {
-          MVertex *v = e->getVertex(lv);
-          if(v->onWhat()->cast2Edge()) {
-            MEdgeVertex *mev = dynamic_cast<MEdgeVertex *>(v);
-            if(mev == nullptr) {
-              Msg::Error("vertex attached to curve %i but not MEdgeVertex",
-                         v->onWhat()->cast2Edge()->tag());
-              abort();
-            }
-          }
-          else if(v->onWhat()->cast2Face()) {
-            MFaceVertex *mfv = dynamic_cast<MFaceVertex *>(v);
-            if(mfv == nullptr) {
-              Msg::Error("vertex attached to face %i but not MFaceVertex",
-                         v->onWhat()->cast2Face()->tag());
-              abort();
-            }
-          }
-          else if(v->onWhat()->cast2Vertex()) {
-          }
-          else {
-            Msg::Error("vertex not attach to a CAD entity");
-          }
+  /* Geometric projection on model */
+
+  /* - projections on curves */
+  #if defined(_OPENMP)
+  #pragma omp parallel for schedule(dynamic)
+  #endif
+  for (size_t e = 0; e < edges.size(); ++e) {
+    GEdge* ge = edges[e];
+    auto it = toProjectOnCurve.find(ge);
+    if (it == toProjectOnCurve.end()) continue;
+
+    double tMin = ge->parBounds(0).low();
+    double tMax = ge->parBounds(0).high();
+
+    for (size_t j = 0; j < it->second.size(); ++j) {
+      MVertex* v = it->second[j];
+      double tGuess = tMin + (tMax-tMin) * double(j+1)/double(it->second.size()+1);
+      GPoint proj = ge->closestPoint(v->point(), tGuess);
+      if(proj.succeeded()) {
+        v->setXYZ(proj.x(), proj.y(), proj.z());
+        v->setParameter(0, proj.u());
+      } else {
+        Msg::Warning("- Edge %i, vertex %li: curve projection failed",
+            ge->tag(), v->getNum());
+      }
+    }
+  }
+
+  GlobalBackgroundMesh *bmesh = nullptr;
+  if(backgroudMeshExists(BMESH_NAME)) {
+    bmesh = &(getBackgroundMesh(BMESH_NAME));
+  }
+  else {
+    Msg::Warning("refine mesh with background projection: no background mesh, "
+                 "using CAD projection (slow)");
+  }
+
+  /* - projections on faces */
+  #if defined(_OPENMP)
+  #pragma omp parallel for schedule(dynamic)
+  #endif
+  for(size_t f = 0; f < faces.size(); ++f) {
+    GFace *gf = faces[f];
+
+    auto it = toProjectOnFace.find(gf);
+    if (it == toProjectOnFace.end()) continue;
+
+    Msg::Info("- Face %i: project %li vertices and smooth the quad mesh ...", 
+        gf->tag(), it->second.size());
+
+    SurfaceProjector* sp = nullptr;
+    if(bmesh) {
+      auto it2 = bmesh->faceBackgroundMeshes.find(gf);
+      if(it2 == bmesh->faceBackgroundMeshes.end()) {
+        Msg::Error("background mesh not found for face %i", gf->tag());
+      } else {
+        /* Get pointers to triangles in the background mesh */
+        std::vector<MTriangle *> triangles(it2->second.triangles.size());
+        for(size_t i = 0; i < it2->second.triangles.size(); ++i) {
+          triangles[i] = &(it2->second.triangles[i]);
+        }
+
+        sp = new SurfaceProjector();
+        bool oki = sp->initialize(gf, triangles);
+        if(!oki) {
+          Msg::Warning("failed to initialize surface projector");
+          delete sp;
+          sp = nullptr;
         }
       }
     }
 
-    /* Check quality of quads */
-    if(gf->quadrangles.size() > 0) {
-      bool qualityOk = true;
-      double sicnMin, sicnAvg;
-      computeSICN(
-        dynamic_cast_vector<MQuadrangle *, MElement *>(gf->quadrangles),
-        sicnMin, sicnAvg);
-      if(sicnMin < 0.) {
-        Msg::Warning("- refine with projection: quality negative (%.3f) on "
-                     "face %i after projection, rollback",
-                     sicnMin, gf->tag());
-        qualityOk = false;
+    /* Project the vertices */
+    bool evalOnCAD = false;
+    bool projOnCad = false;
+    if (gf->haveParametrization() && haveNiceParametrization(gf)) evalOnCAD = true;
+    for (size_t j = 0; j < it->second.size(); ++j) {
+      MVertex* v = it->second[j];
+      GPoint proj;
+      if(sp != nullptr) {
+        proj = sp->closestPoint(v->point().data(), evalOnCAD, projOnCad);
+        if(!proj.succeeded() && gf->haveParametrization()) {
+          double uvg[2] = {0., 0.};
+          proj = gf->closestPoint(v->point(), uvg);
+        }
+      } else {
+        double uvg[2] = {0., 0.};
+        proj = gf->closestPoint(v->point(), uvg);
       }
-      if(!qualityOk) {
-        for(auto &kv : backupPositionLinear) { kv.first->setXYZ(kv.second); }
+      if(proj.succeeded()) {
+        v->setXYZ(proj.x(), proj.y(), proj.z());
+        v->setParameter(0, proj.u());
+        v->setParameter(1, proj.v());
+      } else {
+        Msg::Warning("- Face %i, vertex %li: surface projection failed",
+            gf->tag(), v->getNum());
       }
     }
 
-    /* Smooth geometry (quick) */
-    if(sp != nullptr) {
-      double timeMax = 0.3;
+    /* Quad mesh smoothing */
+    double timeMax = 0.3;
+    optimizeGeometryQuadMesh(gf, sp, timeMax);
+    double sicnMin, sicnAvg;
+    computeSICN(
+      dynamic_cast_vector<MQuadrangle *, MElement *>(gf->quadrangles),
+      sicnMin, sicnAvg);
+    if(sicnMin < 0.) {
+      double timeMax = 2.;
       optimizeGeometryQuadMesh(gf, sp, timeMax);
     }
 
-    if(sp != nullptr) delete sp;
-
-    if(PARANO_VALIDITY) {
-      errorAndAbortIfInvalidVertexInElements(
-        dynamic_cast_vector<MQuadrangle *, MElement *>(gf->quadrangles),
-        "after surf proj");
-    }
+    if (sp) delete sp;
   }
 
-  const bool SHOW_QUADINIT = SHOW_INTERMEDIATE_VIEWS;
-  if(SHOW_QUADINIT) {
+  if(SHOW_INTERMEDIATE_VIEWS) {
     std::vector<MElement *> elements;
     for(GFace *gf : model_faces(gm)) {
       append(elements,
@@ -1899,9 +1899,77 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
     errorAndAbortIfInvalidVertexInModel(gm, "after refine + proj");
   }
 
+  if (DBG_EXPORT) {
+    gm->writeMSH("qqs_subdiv.msh", 4.1);
+  }
+
+  return 0;
+}
+
+
+int checkAndReplaceQuadDominantMesh(GFace* gf, int valenceMaxBdr = 6, int valenceMaxIn = 8) {
+  /* Check valence */
+  unordered_map<MVertex *, int > valence;
+  for (MQuadrangle* q: gf->quadrangles) for (size_t lv = 0; lv < 4; ++lv) {
+    valence[q->getVertex(lv)] += 1;
+  }
+  for (MTriangle* t: gf->triangles) for (size_t lv = 0; lv < 3; ++lv) {
+    valence[t->getVertex(lv)] += 1;
+  }
+  int vMaxInt = 0;
+  int vMaxBdr = 0;
+  for (auto& kv: valence) {
+    if (kv.first->onWhat() == gf) {
+      vMaxInt = std::max(vMaxInt,kv.second);
+    } else if (kv.first->onWhat()->cast2Edge() != nullptr){
+      vMaxBdr = std::max(vMaxBdr,kv.second);
+    }
+  }
+  if (vMaxInt <= valenceMaxIn && vMaxBdr <= valenceMaxBdr) return 0;
+
+  Msg::Warning("- Face %i: quad-tri mesh have valence %i (interior) and %i (bdr), replace it with meshadapt ...",
+      gf->tag(), vMaxInt, vMaxBdr);
+
+  int algo2d = gf->getMeshingAlgo();
+  gf->setMeshingAlgo(ALGO_2D_MESHADAPT);
+  gf->mesh(false);
+  gf->setMeshingAlgo(algo2d);
+
+  if (gf->meshStatistics.status != GFace::DONE) {
+    Msg::Warning("- Face %i: failed to build triangulation", gf->tag());
+    return -1;
+  }
+
+  /* Recombine triangles in quads */
+  recombineIntoQuads(gf, false, 0, false, .1);
+
   return 0;
 }
 
+int replaceBadQuadDominantMeshes(GModel *gm) {
+  std::vector<GFace*> faces = model_faces(gm);
+
+  #if defined(_OPENMP)
+  #pragma omp parallel for schedule(dynamic)
+  #endif
+  for(size_t f = 0; f < faces.size(); ++f) {
+    GFace *gf = faces[f];
+    if(CTX::instance()->mesh.meshOnlyVisible && !gf->getVisibility()) continue;
+    if(CTX::instance()->debugSurface > 0 &&
+       gf->tag() != CTX::instance()->debugSurface)
+      continue;
+    if(gf->getNumMeshElements() == 0) continue;
+
+    int st = checkAndReplaceQuadDominantMesh(gf);
+    if (st != 0) {
+      Msg::Warning("- Face %i: failed to replace bad quad-tri mesh", gf->tag());
+    }
+  }
+
+  return 0;
+}
+
+
 int optimizeQuadMeshWithDiskQuadrangulationRemeshing(GFace *gf)
 {
   if(gf->triangles.size() > 0 || gf->quadrangles.size() == 0) return -1;
@@ -2049,6 +2117,10 @@ int quadMeshingOfSimpleFacesWithPatterns(GModel *gm,
       gm, "global check after face pattern meshing");
   }
 
+  if (DBG_EXPORT) {
+    gm->writeMSH("qqs_simplefaces.msh", 4.1);
+  }
+
   return 0;
 }
 
@@ -2106,6 +2178,10 @@ int optimizeTopologyWithDiskQuadrangulationRemeshing(GModel *gm)
       gm, "global check after disk quadrangulation remeshing");
   }
 
+  if (DBG_EXPORT) {
+    gm->writeMSH("qqs_diskrmsh.msh", 4.1);
+  }
+
   return 0;
 }
 
@@ -2139,6 +2215,8 @@ int optimizeTopologyWithCavityRemeshing(GModel *gm)
     }
   }
 
+  GlobalBackgroundMesh &bmesh = getBackgroundMesh(BMESH_NAME);
+
 #if defined(_OPENMP)
 #pragma omp parallel for schedule(dynamic)
 #endif
@@ -2156,7 +2234,10 @@ int optimizeTopologyWithCavityRemeshing(GModel *gm)
     std::vector<std::pair<SPoint3, int> > singularities;
     bool okg = getSingularitiesFromBackgroundField(gf, singularities);
     if (!okg) {
-      Msg::Warning("- Face %i: failed to get singularities from background field", gf->tag());
+      okg = getSingularitiesFromNewCrossFieldComputation(bmesh, gf, singularities);
+      if (!okg) {
+        Msg::Warning("- Face %i: failed to get singularities", gf->tag());
+      }
     }
 
     bool invertNormals = meshOrientationIsOppositeOfCadOrientation(gf);
@@ -2183,6 +2264,10 @@ int optimizeTopologyWithCavityRemeshing(GModel *gm)
 
   GeoLog::flush();
 
+  if (DBG_EXPORT) {
+    gm->writeMSH("qqs_cavrmsh.msh", 4.1);
+  }
+
   return 0;
 }
 
@@ -2227,6 +2312,11 @@ int RefineMeshWithBackgroundMeshProjection(GModel *gm)
              "RefineMeshWithBackgroundMeshProjection");
   return -10;
 }
+int replaceBadQuadDominantMeshes(GModel *gm) {
+  Msg::Error("Module QUADMESHINGTOOLS required for function "
+             "replaceBadQuadDominantMeshes");
+  return -10;
+}
 
 #endif
 /* endif: QUADMESHINGTOOLS*/
@@ -2317,6 +2407,7 @@ int transferSeamGEdgesVerticesToGFace(GModel *gm)
   return 0;
 }
 
+
 QuadqsContextUpdater::QuadqsContextUpdater()
 {
   algo2d = CTX::instance()->mesh.algo2d;
diff --git a/Mesh/meshQuadQuasiStructured.h b/Mesh/meshQuadQuasiStructured.h
index dd2c38d9f2aa9af5505c1cabca9ff7eb01a8dcdd..70b17bf7d8ec9efc4a57d6a6c5e043ca76122bd7 100644
--- a/Mesh/meshQuadQuasiStructured.h
+++ b/Mesh/meshQuadQuasiStructured.h
@@ -123,10 +123,22 @@ int quadMeshingOfSimpleFacesWithPatterns(GModel *gm,
  *        on the CAD surfaces, using the background mesh for
  *        faster projections.
  *
- * @param gm The model containg the surface meshes
+ * @param gm The model containing the surface meshes
  *
  * @return 0 if success
  */
 int RefineMeshWithBackgroundMeshProjection(GModel *gm);
 
+/**
+ * @brief The initial unstructured quad-tri mesh may contain
+ * very bad configurations (e.g. valence 50+) due to failures
+ * in algo pack. This method replaces them by meshes produced
+ * with algo meshadapt.
+ *
+ * @param gm The model containing the surface meshes
+ *
+ * @return 0 if success
+ */
+int replaceBadQuadDominantMeshes(GModel *gm);
+
 #endif
diff --git a/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.cpp b/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.cpp
index 2414658d913f36cc232f0a49351c87fecbb5ffae..a66dee088f6e0ae6b703c20827888f6c8060bc65 100644
--- a/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.cpp
+++ b/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.cpp
@@ -64,6 +64,7 @@ template <typename Key, typename Hash = robin_hood::hash<Key>, typename KeyEqual
 
 namespace QMT {
   constexpr size_t NO_SIZE_T = (size_t) -1;
+  constexpr uint32_t NO_U32 = (uint32_t) -1;
 
   bool SHOW_QUALITY = false; /* Debug only */
 
@@ -401,14 +402,15 @@ namespace QMT {
     }
 
     v2q.resize(vcount);
-    quads.resize(elements.size());
+    quads.clear();
+    quads.resize(elements.size(),{NO_U32,NO_U32,NO_U32,NO_U32});
     for (size_t f = 0; f < elements.size(); ++f) {
       MElement* q = elements[f];
-      if (q->getNumVertices() != 4) {
-        Msg::Error("buildCondensedStructure: element is not a quad");
+      if (q->getNumVertices() > 4) {
+        Msg::Error("buildCondensedStructure: element with %li vertices not supported", q->getNumVertices());
         return false;
       }
-      for (size_t lv = 0; lv < 4; ++lv) {
+      for (size_t lv = 0; lv < q->getNumVertices(); ++lv) {
         MVertex* v = q->getVertex(lv);
         auto it = old2new.find(v);
         size_t nv;
@@ -436,14 +438,14 @@ namespace QMT {
 
     /* Build the one rings for the free vertices */
     oneRings.resize(freeVertices.size());
-    vector<MElement*> adjQuads;
+    vector<MElement*> adjElts;
     for (size_t v = 0; v < freeVertices.size(); ++v) {
-      adjQuads.resize(v2q[v].size());
+      adjElts.resize(v2q[v].size());
       for (size_t lq = 0; lq < v2q[v].size(); ++lq) {
-        adjQuads[lq] = elements[v2q[v][lq]];
+        adjElts[lq] = elements[v2q[v][lq]];
       }
       std::vector<MVertex*> bnd;
-      bool okb = buildBoundary(adjQuads, bnd);
+      bool okb = buildBoundary(adjElts, bnd);
       if (!okb) {
         Msg::Error("buildCondensedStructure: failed to build boundary for stencil");
         return false;
@@ -454,7 +456,7 @@ namespace QMT {
       /* Start of the stencil */
       MVertex* vp = freeVertices[v];
       MVertex* v0 = NULL;
-      for (MElement* e: adjQuads) {
+      for (MElement* e: adjElts) {
         size_t N = e->getNumVertices();
         for (size_t j = 0; j < N; ++j) {
           MVertex* a = e->getVertex(j);
@@ -476,13 +478,15 @@ namespace QMT {
       for (size_t j = 0; j < bnd.size(); ++j) {
         if (bnd[j] == v0 && j > 0) {
           std::rotate(bnd.begin(),bnd.begin()+j,bnd.end());
+          break;
         }
       }
       if (bnd.front() != v0) {
-        Msg::Warning("buildCondensedStructure: wrong start");
+        Msg::Warning("buildCondensedStructure: wrong start (bnd size: %li)", bnd.size());
         return false;
       }
-      if (bnd.size() < 6 || bnd.size() % 2 != 0) {
+
+      if (bnd.size() < 3) {
         Msg::Warning("buildCondensedStructure: wrong boundary, size %li", bnd.size());
         return false;
       }
@@ -529,6 +533,7 @@ namespace QMT {
       std::vector<SPoint2> quad_uvs = paramOnElement(gf, elements[i]);
       for (size_t lv = 0; lv < 4; ++lv) {
         size_t v = quads[i][lv];
+        if (v == NO_U32) continue;
         if (uvs[v][0] == DBL_MAX) {
           uvs[v][0] = quad_uvs[lv][0];
           uvs[v][1] = quad_uvs[lv][1];
@@ -1147,6 +1152,8 @@ bool movePointWithKernelAndProjection(
     const std::vector<std::array<double,3> >& points,
     const std::vector<size_t>& one_ring_first,
     const std::vector<uint32_t>& one_ring_values,
+    const SmoothingKernel& kernelRegular,
+    const SmoothingKernel& kernelIrregular,
     bool project,
     SurfaceProjector* sp,
     vec3& newPos,
@@ -1160,7 +1167,7 @@ bool movePointWithKernelAndProjection(
   size_t n = one_ring_first[v+1] - one_ring_first[v];
   GPoint proj;
   double stDx = 0.;
-  if (n == 8) { /* regular vertex */
+  if (n == 8 && kernelRegular == SmoothingKernel::WinslowFDM) { /* Winslow for regular vertex */
     /* Extract geometric stencil */
     std::array<vec3,8> stencil = fillStencilRegular(v, points, one_ring_first, one_ring_values);
 
@@ -1180,12 +1187,26 @@ bool movePointWithKernelAndProjection(
   } else { /* irregular vertex */
     std::vector<vec3> stencilIrreg(8);
 
+    bool angleBased = false;
+    if (n > 4 && n % 2 == 0) {
+      if (n == 8 && kernelRegular == AngleBased) {
+        angleBased = true;
+      } else if (kernelIrregular == AngleBased) {
+        angleBased = true;
+      }
+    }
+
     /* Extract geometric stencil */
-    constexpr bool oneOverTwo = true; /* angle-based does not use diagonals */
+    bool oneOverTwo = angleBased; /* angle-based does not use diagonals */
     fillStencilIrregular(v, points, one_ring_first, one_ring_values, stencilIrreg, oneOverTwo);
 
     /* Smoothing (in 3D, not on surface) */
-    bool ok = kernelAngleBased(points[v], stencilIrreg, newPos);
+    bool ok = false;
+    if (angleBased) {
+      ok = kernelAngleBased(points[v], stencilIrreg, newPos);
+    } else {
+      ok = kernelLaplacian(stencilIrreg, newPos);
+    }
     if (!ok) return false;
 
     if (project) { /* Projection on surface */
@@ -1224,7 +1245,7 @@ bool kernelLoopWithProjection(
     bool okc = buildCondensedStructure(patch.elements,patch.intVertices,
         old2new,new2old, quads,v2q,oneRings,points);
     if (!okc) {
-      Msg::Error("buildCondensedStructure: failed to build condensed representation");
+      Msg::Warning("kernelLoopWithProjection: failed to build condensed representation");
       return false;
     }
     compress(oneRings, one_ring_first, one_ring_values);
@@ -1275,13 +1296,13 @@ bool kernelLoopWithProjection(
       if (point_uvs.size()) { newUv = point_uvs[v]; }
 
       bool ok = movePointWithKernelAndProjection(v, points, one_ring_first, one_ring_values,
-          project, opt.sp, newPos, newUv);
+          opt.kernelRegular, opt.kernelIrregular, project, opt.sp, newPos, newUv);
       if (ok) {
         double dx = length(newPos - points[v]);
         sum_dx += dx;
         /* Modify the coordinates */
         points[v] = newPos;
-        point_uvs[v] = newUv;
+        if (point_uvs.size()) point_uvs[v] = newUv;
         if (opt.localLocking && dx < opt.dxLocalMax*localAvgSize[v]) {
           locked[v] = true;
         } else {
@@ -1295,11 +1316,23 @@ bool kernelLoopWithProjection(
     if (iter == 0) {
       sum_dx0 = sum_dx;
     } else {
-      if (sum_dx < opt.dxGlobalMax * sum_dx0) break;
-      if (nMoved == 0) break;
+      if (sum_dx < opt.dxGlobalMax * sum_dx0) {
+        Msg::Debug("kernelLoopWithProjection: stop at iter %li/%li because sum_dx = %f < %f", 
+            iter, opt.outerLoopIterMax, sum_dx, opt.dxGlobalMax*sum_dx0);
+        break;
+      }
+      if (nMoved == 0) {
+        Msg::Debug("kernelLoopWithProjection: stop at iter %li/%li because no point moved", 
+            iter, opt.outerLoopIterMax);
+        break;
+      }
     }
 
-    if (Cpu() - t0 > opt.timeMax) break;
+    if (Cpu() - t0 > opt.timeMax) {
+      Msg::Debug("kernelLoopWithProjection: stop at iter %li/%li because timeout (%f sec)", 
+          iter, opt.outerLoopIterMax, Cpu()-t0);
+      break;
+    }
   }
 
   /* Update the positions */
@@ -1476,6 +1509,8 @@ bool patchProjectOnSurface(GFaceMeshPatch& patch, SurfaceProjector* sp) {
 
 bool optimizeGeometryQuadMesh(GFace* gf, SurfaceProjector* sp, double timeMax)
 {
+  // TODO FIXME: replace by optimizeGeometryQuadTriMesh when it works well
+
   if (!gf->haveParametrization() && sp == nullptr) {
     Msg::Error("optimize geometry: face %i, no CAD and no surface projector", gf->tag());
     return false;
@@ -1579,3 +1614,113 @@ bool optimizeGeometryQuadMesh(GFace* gf, SurfaceProjector* sp, double timeMax)
 
   return true;
 }
+
+bool optimizeGeometryQuadTriMesh(GFace* gf, SurfaceProjector* sp, double timeMax) {
+  // TODO FIXME: replace by optimizeGeometryQuadMesh when it works well
+
+  if (!gf->haveParametrization() && sp == nullptr) {
+    Msg::Error("optimize geometry: face %i, no CAD and no surface projector", gf->tag());
+    return false;
+  }
+
+  if (gf->quadrangles.size() == 0) {
+    Msg::Error("optimize geometry: face %i, no quads", gf->tag());
+    return false;
+  }
+
+  bool forceEvenIfBadBoundary = true;
+  GFaceMeshPatch patch;
+  std::vector<MElement*> elements = dynamic_cast_vector<MQuadrangle*,MElement*>(gf->quadrangles);
+  append(elements, dynamic_cast_vector<MTriangle*,MElement*>(gf->triangles));
+  bool okp = patchFromElements(gf, elements, patch, forceEvenIfBadBoundary);
+  if (!okp) {
+    Msg::Warning("optimize geometry: face %i, failed to build patch", gf->tag());
+    return false;
+  }
+
+  bool debugCreateViews = DBG_VIZU_G;
+  if (Msg::GetVerbosity() >= 999
+      && CTX::instance()->debugSurface == patch.gf->tag()) {
+    debugCreateViews = true;
+  }
+  const int rdi = (int)(((double)rand()/RAND_MAX)*1e4); /* only to get a random name for debugging */
+
+  int countMax = 3;
+  bool running = true;
+  size_t niter = 50;
+  int count = 0;
+  if (gf->geomType() == GEntity::Plane) {
+    /* Much faster projection, we can smooth */
+    niter = 100;
+  }
+
+  double t0 = Cpu();
+
+  while (running && count < countMax) {
+    running = false;
+    count += 1;
+    if (Cpu() - t0 > timeMax) {
+      Msg::Debug("optimize geometry: face %i, reached time limit", gf->tag());
+      break;
+    }
+
+    PatchGeometryBackup backup(patch);
+
+    double minSICNb = DBL_MAX;
+    double avgSICNb = 0.;
+    computeSICN(patch.elements, minSICNb, avgSICNb);
+
+    if (debugCreateViews) {
+      std::string method = "GFace_K";
+      GeoLog::add(patch.elements, "optim_"+method+"_IN_"+std::to_string(rdi) + "_" + std::to_string(minSICNb));
+    }
+
+    GeomOptimStats stats;
+    GeomOptimOptions opt;
+    opt.sp = sp;
+    opt.outerLoopIterMax = niter;
+    opt.timeMax = timeMax;
+    opt.localLocking = true;
+    opt.dxGlobalMax = 1.e-3;
+    opt.dxLocalMax = 1.e-5;
+    opt.force3DwithProjection = true;
+    opt.withBackup = false;
+    double t1 = Cpu();
+    bool okk = kernelLoopWithProjection(patch, opt, stats);
+    if (!okk) {
+      return false;
+    }
+    double etime = Cpu() - t1;
+
+    double minSICNa = DBL_MAX;
+    double avgSICNa = 0.;
+    computeSICN(patch.elements, minSICNa, avgSICNa);
+
+    bool keep = true;
+    if (minSICNa < minSICNb && avgSICNa < avgSICNb) keep = false;
+    if (minSICNa < 0.1 && minSICNa < minSICNb) keep = false;
+    if (minSICNa < 0.75 * minSICNb) keep = false;
+
+    if (debugCreateViews) {
+      std::string method = "GFace_K";
+      GeoLog::add(patch.elements, "optim_"+method+"_OUT_"+std::to_string(rdi) + "_" + std::to_string(minSICNa) + "_keep=" + std::to_string(keep));
+    }
+
+    if (keep) {
+      Msg::Debug("- Face %i: kernel smoothing (Mix Winslow/Angle, explicit with projections, %li vertices, %li iter max, %.3f sec), SICN min: %f -> %f, avg: %f -> %f",
+          gf->tag(),gf->mesh_vertices.size(), niter,etime,minSICNb,minSICNa,avgSICNb,avgSICNa);
+      if (minSICNa >= minSICNb && 0.99*avgSICNa > avgSICNb && etime < 1) {
+        niter *= 1.5;
+        running = true;
+      }
+    } else {
+      Msg::Debug("- Face %i: worst quality after global smoothing (Mix Winslow/Angle, explicit, %li iter max), roll back (SICN min: %f -> %f, avg: %f -> %f)",
+          gf->tag(),niter,minSICNb,minSICNa,avgSICNb,avgSICNa);
+      backup.restore();
+      break;
+    }
+  }
+
+  return true;
+
+}
diff --git a/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.h b/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.h
index cdebf5df510fc1ab732dd18bcda4dca0ef29b07f..98539518456abcb368817080753a0e89b5e51a93 100644
--- a/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.h
+++ b/contrib/QuadMeshingTools/qmtMeshGeometryOptimization.h
@@ -53,6 +53,12 @@ int patchOptimizeGeometryGlobal(
  * - Winslow if regular vertex
  * - Angle-based smoothing is irregular vertex
  */
+enum SmoothingKernel {
+  Laplacian,
+  WinslowFDM,
+  AngleBased
+};
+
 struct GeomOptimOptions {
   double useDmoIfSICNbelow = 0.1; /* use the DMO kernel if SICN quality below threshold */
   size_t outerLoopIterMax = 100; /* maximum number of loops over all vertices */
@@ -68,6 +74,8 @@ struct GeomOptimOptions {
   double qualityRangeMax = 0.8;
   bool withBackup = true; /* save the geometry before, restore if quality decreased */
   bool force3DwithProjection = false;
+  SmoothingKernel kernelRegular = SmoothingKernel::WinslowFDM;
+  SmoothingKernel kernelIrregular = SmoothingKernel::Laplacian;
 };
 
 /**
@@ -120,3 +128,15 @@ bool patchProjectOnSurface(GFaceMeshPatch& patch, SurfaceProjector* sp = nullptr
  * @return true if success
  */
 bool optimizeGeometryQuadMesh(GFace* gf, SurfaceProjector* sp = nullptr, double timeMax = DBL_MAX);
+
+/**
+ * @brief High-level function which try to make good parameter choices
+ *        automatically.
+ *
+ * @param gf The face containing the quad-tri mesh to smooth
+ * @param sp Surface projector (faster than CAD projection)
+ * @param timeMax Time budget for the smoothing
+ *
+ * @return true if success
+ */
+bool optimizeGeometryQuadTriMesh(GFace* gf, SurfaceProjector* sp = nullptr, double timeMax = DBL_MAX);
diff --git a/contrib/QuadMeshingTools/qmtMeshUtils.cpp b/contrib/QuadMeshingTools/qmtMeshUtils.cpp
index 079b85e53e14a1bef3e8fc570f1af696de1c8154..7fb25e00ad8bbcb097a8be9b1b15f9477f8310e2 100644
--- a/contrib/QuadMeshingTools/qmtMeshUtils.cpp
+++ b/contrib/QuadMeshingTools/qmtMeshUtils.cpp
@@ -810,24 +810,59 @@ std::vector<MTriangle*> trianglesFromQuads(const std::vector<MQuadrangle*>& quad
   return tris;
 }
 
-bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, bool& requireDelete) {
+bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, std::vector<MTriangle*>& requireDelete,
+    bool copyExisting) {
   triangles.clear();
-  requireDelete = false;
+  requireDelete.clear();
 
   /* Existing pure tri or pure quad mesh on the GFace */
   if (gf->triangles.size() > 0 && gf->quadrangles.size() == 0) {
-    triangles = gf->triangles;
-    requireDelete = false;
+    if (copyExisting) {
+      triangles.resize(gf->triangles.size());
+      for (size_t i = 0; i < gf->triangles.size(); ++i) {
+        triangles[i] = new MTriangle(
+            gf->triangles[i]->getVertex(0),
+            gf->triangles[i]->getVertex(1),
+            gf->triangles[i]->getVertex(2));
+      }
+      requireDelete = triangles;
+    } else {
+      triangles = gf->triangles;
+    }
     return true;
   } else if (gf->triangles.size() == 0 && gf->quadrangles.size() > 0) {
     triangles = trianglesFromQuads(gf->quadrangles);
-    requireDelete = true;
+    requireDelete = triangles;
+    return true;
+  } else if (gf->triangles.size() > 0 && gf->quadrangles.size() > 0) {
+    /* quad-dominant mesh */
+    /* - triangles from triangles */
+    if (copyExisting) {
+      triangles.resize(gf->triangles.size());
+      for (size_t i = 0; i < gf->triangles.size(); ++i) {
+        triangles[i] = new MTriangle(
+            gf->triangles[i]->getVertex(0),
+            gf->triangles[i]->getVertex(1),
+            gf->triangles[i]->getVertex(2));
+      }
+      requireDelete = triangles;
+    } else {
+      triangles = gf->triangles;
+    }
+    /* - triangles from quads */
+    std::vector<MTriangle*> quadSplit = trianglesFromQuads(gf->quadrangles);
+    append(triangles, quadSplit);
+    append(requireDelete,quadSplit);
     return true;
   }
 
   /* Check if there is the quadqs background mesh */
   const std::string BMESH_NAME = "bmesh_quadqs";
   if (backgroudMeshExists(BMESH_NAME)) {
+    if (copyExisting) {
+      Msg::Error("copyExisting not supported here (but simple to fix)");
+      return false;
+    }
     GlobalBackgroundMesh& bmesh = getBackgroundMesh(BMESH_NAME);
     auto it = bmesh.faceBackgroundMeshes.find(gf);
     if (it != bmesh.faceBackgroundMeshes.end() && it->second.triangles.size() > 0) {
@@ -836,7 +871,6 @@ bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, bool& requ
       for (size_t i = 0; i < it->second.triangles.size(); ++i) {
         triangles[i] = &(it->second.triangles[i]);
       }
-      requireDelete = false;
       return true;
     }
   }
@@ -845,12 +879,15 @@ bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, bool& requ
   for (auto& bmesh: global_bmeshes) {
     auto it = bmesh->faceBackgroundMeshes.find(gf);
     if (it != bmesh->faceBackgroundMeshes.end() && it->second.triangles.size() > 0) {
+      if (copyExisting) {
+        Msg::Error("copyExisting not supported here (but simple to fix)");
+        return false;
+      }
       /* Get pointers to triangles in the background mesh */
       triangles.resize(it->second.triangles.size());
       for (size_t i = 0; i < it->second.triangles.size(); ++i) {
         triangles[i] = &(it->second.triangles[i]);
       }
-      requireDelete = false;
       return true;
     }
   }
@@ -864,14 +901,9 @@ bool fillSurfaceProjector(GFace* gf, SurfaceProjector* sp) {
     abort();
   }
 
-  if (gf->geomType() == GFace::GeomType::Sphere) {
-    bool oka = sp->setAnalyticalProjection(gf);
-    return oka;
-  }
-
-  bool deleteTheTris = false;
   std::vector<MTriangle*> triangles;
-  bool okgt = getGFaceTriangles(gf, triangles, deleteTheTris);
+  std::vector<MTriangle*> trianglesToDel;
+  bool okgt = getGFaceTriangles(gf, triangles, trianglesToDel);
   if (!okgt) {
     Msg::Error("fillSurfaceProjector: case not supported, no triangles");
     return false;
@@ -882,9 +914,7 @@ bool fillSurfaceProjector(GFace* gf, SurfaceProjector* sp) {
     Msg::Error("failed to initialize the surface projector");
   }
 
-  if (deleteTheTris) {
-    for (MTriangle* t: triangles) delete t;
-  }
+  for (MTriangle* t: trianglesToDel) delete t;
 
   return true;
 }
@@ -1013,9 +1043,9 @@ bool fillGFaceInfo(GFace* gf, GFaceInfo& info) {
   info.intSumVal3mVal5 = 0;
 
 
-  bool deleteTheTris = false;
   std::vector<MTriangle*> triangles;
-  bool okgt = getGFaceTriangles(gf, triangles, deleteTheTris);
+  std::vector<MTriangle*> trianglesToDel;
+  bool okgt = getGFaceTriangles(gf, triangles, trianglesToDel);
   if (!okgt) {
     Msg::Error("fillSurfaceProjector: case not supported, no triangles");
     return false;
@@ -1023,9 +1053,7 @@ bool fillGFaceInfo(GFace* gf, GFaceInfo& info) {
 
   bool okr = fillGFaceInfo(gf, info, triangles);
 
-  if (deleteTheTris) {
-    for (MTriangle* t: triangles) delete t;
-  }
+  for (MTriangle* t: trianglesToDel) delete t;
 
   return okr;
 }
diff --git a/contrib/QuadMeshingTools/qmtMeshUtils.h b/contrib/QuadMeshingTools/qmtMeshUtils.h
index ebdde427dd480b8a01580ed70ac19164a4d11119..cb97c56e176b740d480ab319085837680e87928e 100644
--- a/contrib/QuadMeshingTools/qmtMeshUtils.h
+++ b/contrib/QuadMeshingTools/qmtMeshUtils.h
@@ -75,8 +75,12 @@ std::vector<SPoint2> paramOnElement(GFace* gf, MElement* e);
 /* warning: triangles are allocated, should be delete by the caller */
 std::vector<MTriangle*> trianglesFromQuads(const std::vector<MQuadrangle*>& quads);
 
-bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, bool& requireDelete);
-
+/* Find a way to get triangles associated to GFace (will look into current mesh, can split quads,
+ * will check in background mesh)
+ * If new triangles are allocated (e.g. by splitting quads), they are added to requireDelete
+ * If copyExisting: existing triangles are copied into new ones (added to requireDelete) */
+bool getGFaceTriangles(GFace* gf, std::vector<MTriangle*>& triangles, std::vector<MTriangle*>& requireDelete,
+    bool copyExisting = false);
 
 bool fillSurfaceProjector(GFace* gf, SurfaceProjector* sp);
 
diff --git a/contrib/QuadMeshingTools/qmtQuadCavityRemeshing.cpp b/contrib/QuadMeshingTools/qmtQuadCavityRemeshing.cpp
index 055ecb3d3aa0ea07577bbca1567e9836c492f52b..09f4352d1b0b3a0f8cfe029e3cc43e1c53318527 100644
--- a/contrib/QuadMeshingTools/qmtQuadCavityRemeshing.cpp
+++ b/contrib/QuadMeshingTools/qmtQuadCavityRemeshing.cpp
@@ -15,10 +15,12 @@
 // #include <cmath>
 #include <queue>
 #include <unordered_set>
+#include <cassert>
 // #include <algorithm>
 
 /* Gmsh includes */
 #include "GmshMessage.h"
+#include "Context.h"
 #include "OS.h"
 #include "GVertex.h"
 #include "GEdge.h"
@@ -678,8 +680,6 @@ namespace QMT {
       std::vector<MElement*>& newQuads,
       std::vector<MVertex*>& newVertices){
 
-    bool haveParam = gf->haveParametrization();
-
     std::vector<MVertex*> s3r = reverseVector(s3);
     std::vector<MVertex*> s4r = reverseVector(s4);
     std::vector< std::vector<MVertex*> > grid(s1.size());
@@ -711,6 +711,7 @@ namespace QMT {
           SVector3 p = (1.-v) * s1u + v * s3u + (1.-u) * s4v + u * s2v
             - ((1.-u)*(1.-v)*c00 + u*v*c11 + u * (1.-v) * c10 + (1.-u)*v*c01);
           double uv[2] = {0.,0.};
+          // TODO: interpolate uv if param available
           MVertex *vNew = new MFaceVertex(p.x(),p.y(),p.z(),gf,uv[0],uv[1]);
           newVertices.push_back(vNew);
           grid[i][j] = vNew;
@@ -748,7 +749,7 @@ namespace QMT {
 
     /* Associate exising vertices to pattern sides */
     SVector3 center(0.,0.,0.);
-    if (oldCenter) {
+    if (oldCenter != nullptr) {
       center = oldCenter->point();
     } else {
       double csum = 0.;
@@ -800,15 +801,13 @@ namespace QMT {
 
     /* Create vertices on internal points */
     for (size_t v = 0; v < P.n; ++v) if (!P.vOnBdr[v]) {
-      GPoint pp;
-      if (oldCenter) {
-        double uv[2];
-        oldCenter->getParameter(0,uv[0]);
-        oldCenter->getParameter(1,uv[1]);
-        pp = GPoint(center.x(),center.y(),center.z(),gf,uv[0],uv[1]);
-      } else {
-        pp = GPoint(center.x(),center.y(),center.z(),gf,0,0);
+      double uvc[2] = {0.,0.};
+      if (oldCenter != nullptr && gf->haveParametrization()
+          && dynamic_cast<MFaceVertex*>(oldCenter)) {
+        oldCenter->getParameter(0,uvc[0]);
+        oldCenter->getParameter(1,uvc[1]);
       }
+      GPoint pp = GPoint(center.x(),center.y(),center.z(),gf,uvc[0],uvc[1]);
 
       bool moveTowardBdr = true;
       if (moveTowardBdr) {
@@ -831,13 +830,15 @@ namespace QMT {
 
       double uv[2] = {0.,0.};
       MVertex *sing = new MFaceVertex(pp.x(),pp.y(),pp.z(),gf,pp.u(),pp.v());
-      GPoint proj = gf->closestPoint(sing->point(),uv);
-      if (proj.succeeded()) {
-        sing->x() = proj.x();
-        sing->y() = proj.y();
-        sing->z() = proj.z();
-        sing->setParameter(0,proj.u());
-        sing->setParameter(1,proj.v());
+      if (gf->haveParametrization()) {
+        GPoint proj = gf->closestPoint(sing->point(),uv);
+        if (proj.succeeded()) {
+          sing->x() = proj.x();
+          sing->y() = proj.y();
+          sing->z() = proj.z();
+          sing->setParameter(0,proj.u());
+          sing->setParameter(1,proj.v());
+        }
       }
 
       newVertices.push_back(sing);
@@ -965,6 +966,7 @@ namespace QMT {
   struct GrowingCavity {
     std::vector<bool> quadIsInside;      /* fixed size: farmer.quads.size() */
     std::vector<uint32_t> quadEdgeIsBdr; /* fixed size: 4*farmer.quads.size() */
+    std::vector<uint32_t> vertexValenceInside; /* fixed size: farmer.vertices.size() */
     unordered_set<uint32_t> verticesInt; /* growing with cavity */
     unordered_set<uint32_t> verticesBdr; /* growing with cavity */
     unordered_set<uint32_t> elements;    /* growing with cavity */
@@ -974,6 +976,7 @@ namespace QMT {
     void clear() {
       quadIsInside.clear();
       quadEdgeIsBdr.clear();
+      vertexValenceInside.clear();
       verticesInt.clear();
       verticesBdr.clear();
       elements.clear();
@@ -1004,6 +1007,9 @@ namespace QMT {
       std::vector<bool> vertexForbidden;
       std::vector<bool> quadEdgeForbidden;
 
+      /* outside allocation for functions */
+      vector<uint32_t> _vertex2pos;
+
     public:
       CavityFarmer(GFace* gf_):gf(gf_) {};
 
@@ -1234,6 +1240,7 @@ namespace QMT {
             return false;
           }
           for (size_t le = 0; le < 4; ++le) {
+            cav.vertexValenceInside[quads[f][le]] += 1;
             touched.push_back(quads[f][le]);
             uint32_t opp = adjacent[4*f+le];
             if (cav.quadEdgeIsBdr[4*f+le]) {
@@ -1258,8 +1265,9 @@ namespace QMT {
         sort_unique(touched);
 
         for (uint32_t v: touched) {
-          uint32_t valIn, valOut;
-          vertexValence(v, valIn, valOut);
+          uint32_t n = v2q_first[v+1]-v2q_first[v];
+          uint32_t valIn = cav.vertexValenceInside[v];
+          uint32_t valOut = n - valIn;
           if (valOut > 0) {
             cav.verticesBdr.insert(v);
           } else if (valOut == 0 && valIn > 0) {
@@ -1301,9 +1309,10 @@ namespace QMT {
         nConvexCorners = 0;
         double ir = 0.;
         for (uint32_t v: cav.verticesInt) {
-          uint32_t valIn, valOut;
-          vertexValence(v, valIn, valOut);
-          if (valOut) {
+          uint32_t n = v2q_first[v+1]-v2q_first[v];
+          uint32_t valIn = cav.vertexValenceInside[v];
+          uint32_t valOut = n - valIn;
+          if (valOut > 0) {
             Msg::Error("weird: cavity interior vertex has valence-out=%i",valOut);
           }
           if (valIn != 4) {
@@ -1312,8 +1321,8 @@ namespace QMT {
           }
         }
         for (uint32_t v: cav.verticesBdr) {
-          uint32_t valIn, valOut;
-          vertexValence(v, valIn, valOut);
+          uint32_t n = v2q_first[v+1]-v2q_first[v];
+          uint32_t valIn = cav.vertexValenceInside[v];
           if (valIn == 1) {
             nConvexCorners += 1;
           } else if (valIn > 2) {
@@ -1328,8 +1337,9 @@ namespace QMT {
 
       bool isValidConvexCavity() {
         for (uint32_t v: cav.verticesBdr) {
-          uint32_t valIn, valOut;
-          vertexValence(v, valIn, valOut);
+          uint32_t n = v2q_first[v+1]-v2q_first[v];
+          uint32_t valIn = cav.vertexValenceInside[v];
+          uint32_t valOut = n - valIn;
           if (valOut == 1) {
             MVertex* vp = vertices[v];
             GFace* gf = vp->onWhat()->cast2Face();
@@ -1359,12 +1369,16 @@ namespace QMT {
           running = false;
           std::vector<uint32_t> quadsAtConcavities;
           for (uint32_t v: cav.verticesBdr) {
-            uint32_t valIn, valOut;
-            uint32_t qOut = NO_U32;
-            vertexValence(v, valIn, valOut, &qOut);
+            uint32_t n = v2q_first[v+1]-v2q_first[v];
+            uint32_t valIn = cav.vertexValenceInside[v];
+            uint32_t valOut = n - valIn;
+
             /* Add exterior quad if cavity bdr vertex out-valence is one
              * and cavity bdr vertex is inside the GFace */
             if (valOut == 1) {
+              uint32_t qOut = NO_U32;
+              vertexValence(v, valIn, valOut, &qOut);
+              assert(qOut != NO_U32);
               GFace* pgf = vertices[v]->onWhat()->cast2Face();
               if (pgf == nullptr) continue;
               quadsAtConcavities.push_back(qOut);
@@ -1386,7 +1400,10 @@ namespace QMT {
         /* Prepare the cavity boundary edges to make
          * a continuous loop on them */
         vector<uint32_t> pos2vertex;
-        vector<uint32_t> vertex2pos(vertices.size(),NO_U32); /* allocation outside ? */
+        if (_vertex2pos.size() != vertices.size()) {
+          _vertex2pos.resize(vertices.size());
+        }
+        std::fill(_vertex2pos.begin(),_vertex2pos.end(), NO_U32); /* _vertex2pos allocated outside */
         vector<uint32_t> pos2nextVertex;
         vector<bool> isCorner;
         pos2vertex.reserve(cav.verticesBdr.size());
@@ -1405,15 +1422,15 @@ namespace QMT {
                 edgeFound += 1;
                 if (edgeFound == 1) { /* edgeFound >= 2 means non manifold */
                   uint32_t v2 = quads[q][(le+1)%4];
-                  if (vertex2pos[v] != NO_U32) {
+                  if (_vertex2pos[v] != NO_U32) {
                     Msg::Debug("cavity farmer: get sides, vertex already assigned, not manifold bdr loop");
                     return false;
                   }
-                  vertex2pos[v] = pos2vertex.size();
+                  _vertex2pos[v] = pos2vertex.size();
                   pos2vertex.push_back(v);
                   pos2nextVertex.push_back(v2);
-                  uint32_t valIn, valOut;
-                  vertexValence(v, valIn, valOut);
+                  uint32_t n = v2q_first[v+1]-v2q_first[v];
+                  uint32_t valIn = cav.vertexValenceInside[v];
                   if (valIn == 1) {
                     isCorner.push_back(true);
                     Ncorners += 1;
@@ -1471,7 +1488,7 @@ namespace QMT {
             // abort();
           }
 
-          uint32_t pos = vertex2pos[v];
+          uint32_t pos = _vertex2pos[v];
           if (isCorner[pos]) {
             nCornerVisited += 1;
             sides.resize(sides.size()+1);
@@ -1636,6 +1653,7 @@ namespace QMT {
         cavBackup.clear();
         cav.quadIsInside.resize(quads.size(),false);
         cav.quadEdgeIsBdr.resize(4*quads.size(),false);
+        cav.vertexValenceInside.resize(vertices.size(),0);
 
         /* constrain the growth */
         vertexForbidden.clear();
@@ -1849,6 +1867,28 @@ namespace QMT {
     return false;
   }
 
+  bool keepCavityRemeshing(
+      double sicnMinBefore, 
+      double sicnAvgBefore,
+      double sicnMinAfter, 
+      double sicnAvgAfter) {
+    double boldness = CTX::instance()->mesh.quadqsRemeshingBoldness;
+    if (boldness < 0.) boldness = 0.;
+    if (boldness > 1.) boldness = 1.;
+    if (boldness == 1. && sicnMinAfter > 0.) {
+      return true;
+    } else if (boldness == 0. && sicnMinAfter < sicnMinBefore) {
+      return false;
+    }
+    if (boldness == 0.501) { /* backward compatibility */
+      return ((sicnMinAfter > 0.3 && sicnAvgAfter > 0.5) || sicnMinAfter > sicnMinBefore);
+    }
+    double acceptMin = (1.-boldness) * sicnMinBefore;
+    double acceptAvg = (1.-boldness) * sicnAvgBefore * 0.5;
+    if (sicnMinAfter > acceptMin && sicnAvgAfter > acceptAvg) return true;
+
+    return false;
+  }
 
   bool remeshCavity(GFace* gf, RemeshableCavity& cav, QuadqsGFaceContext& ctx,
       const std::vector<MElement*>& cavityNeighborsForSmoothing) {
@@ -1862,7 +1902,7 @@ namespace QMT {
     }
 
     /* Call the remeshing code */
-    const double minSICNafer = 0.1;
+    const double minSICNafer = 0.;
     GFaceMeshDiff diff;
     int st = remeshPatchWithQuadPattern(gf, cav.patternNoAndRot, cav.sides, cav.elements,
         cav.intVertices, minSICNafer, ctx.invertNormalsForQuality, ctx.sp, diff);
@@ -1877,7 +1917,7 @@ namespace QMT {
     computeSICN(diff.before.elements, sicnMin, sicnAvg);
     double sicnMinAfter, sicnAvgAfter;
     computeSICN(diff.after.elements, sicnMinAfter, sicnAvgAfter);
-    if ((sicnMin > 0.3 && sicnAvg > 0.5) || sicnMinAfter > sicnMin) {
+    if (keepCavityRemeshing(sicnMin, sicnAvg, sicnMinAfter, sicnAvgAfter)) {
       constexpr bool verifyTopology = true;
       size_t neb = diff.before.elements.size();
       size_t nea = diff.after.elements.size();
@@ -2591,8 +2631,10 @@ int remeshPatchWithQuadPattern(
     opt.invertCADNormals = invertNormalsForQuality;
     opt.outerLoopIterMax = 100;
     opt.timeMax = 0.5;
+    opt.kernelRegular = SmoothingKernel::WinslowFDM;
+    opt.kernelIrregular = SmoothingKernel::Laplacian;
     opt.localLocking = true;
-    opt.dxGlobalMax = 1.e-3;
+    opt.dxGlobalMax = 0.;
     opt.dxLocalMax = 1.e-5;
     opt.withBackup = false;
     patchOptimizeGeometryWithKernel(diff.after, opt, stats);
@@ -2806,7 +2848,7 @@ int meshFaceWithGlobalPattern(GFace* gf, bool invertNormalsForQuality, double mi
 
   if (sides.size() == 0) {
     Msg::Debug("face %i (%li edges), failed to build sides",gf->tag(),edges.size());
-    DBG(disk);
+    // DBG(disk);
     return -1;
   }
 
diff --git a/contrib/QuadMeshingTools/qmtSizeMap.cpp b/contrib/QuadMeshingTools/qmtSizeMap.cpp
index f7f427de478e4b9929fe5774278186bf16b08893..8c7aef496765bbe20ce2d16c2f08e3d1d6886bd2 100644
--- a/contrib/QuadMeshingTools/qmtSizeMap.cpp
+++ b/contrib/QuadMeshingTools/qmtSizeMap.cpp
@@ -93,6 +93,7 @@ using namespace QMT;
 
 int computeMinimalSizeOnCurves(
     GlobalBackgroundMesh& gbm,
+    bool clampMinWithTriEdges,
     std::unordered_map<MVertex*,double>& minSize) {
   Msg::Debug("compute minimal size on curves (using background mesh) ...");
   /* Important information: all mesh elements are queried in the GlobalBackgroundMesh,
@@ -177,6 +178,27 @@ int computeMinimalSizeOnCurves(
     }
   }
 
+  if (clampMinWithTriEdges) {
+    for (GFace* gf: model_faces(gm)) {
+      auto it = gbm.faceBackgroundMeshes.find(gf);
+      if (it == gbm.faceBackgroundMeshes.end()) {
+        Msg::Error("face %i not found in background mesh", gf->tag());
+        continue;
+      }
+      for (MTriangle& t: it->second.triangles) for (size_t le = 0; le < 3; ++le) {
+        MVertex* v1 = t.getVertex(le);
+        MVertex* v2 = t.getVertex((le+1)%3);
+        double len = v1->distance(v2);
+        if (v1->onWhat()->cast2Face() == nullptr) {
+          auto itv = minSize.find(v1);
+          if (itv != minSize.end()) {
+            itv->second = std::max(itv->second,len);
+          }
+        }
+      }
+    }
+  }
+
   return 0;
 }
 
diff --git a/contrib/QuadMeshingTools/qmtSizeMap.h b/contrib/QuadMeshingTools/qmtSizeMap.h
index 115d377bb5f8d363aa25035b679a31a7484176b7..8396fcc32db6a902a718891ff42a34d43888c1a1 100644
--- a/contrib/QuadMeshingTools/qmtSizeMap.h
+++ b/contrib/QuadMeshingTools/qmtSizeMap.h
@@ -24,16 +24,18 @@ class GlobalBackgroundMesh;
  *        - the prescribed mesh size
  *
  * @param[in] gbm Background mesh on which to compute the minimal sizes
+ * @param[in] clampMinWithTriEdges If true, the minimum length is the maximum
+ * of the previously computed size (from CAD) and local background mesh triangle size.
+ * This option is useful to avoid over-refinement.
  * @param[out] minSize the minimal distance, for each MVertex of GVertex / GEdge 
  *
  * @return 0 if success
  */
 int computeMinimalSizeOnCurves(
     GlobalBackgroundMesh& gbm,
+    bool clampMinWithTriEdges,
     std::unordered_map<MVertex*,double>& minSize);
 
-
-
 /**
  * @brief One way smoothing to get a smooth scalar field where
  *        the gradient is inferior to gradientMax, the output