diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index d64c66379b212df43d57302e21c3f0716f8b1982..e786431eb919d35a8d63ed1acd13c5a11d6ba4bd 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -186,7 +186,8 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
     p.adaptBlobLayerFact = (int) o->value[10]->value();
     p.adaptBlobDistFact = o->value[11]->value();
     p.optPrimSurfMesh = false;
-    HighOrderMeshOptimizer(GModel::current(), p);
+    // HighOrderMeshOptimizer(GModel::current(), p);
+    HighOrderMeshOptimizerNew(GModel::current(), p);
     break;
   }
   case 1: {                                                               // Elastic analogy
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index e384e1d1494519bdac59825dc7b5c19f5e5311b4..a12635cdf81aa724110d396c5bd977fda06209f8 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -9,6 +9,7 @@
 #include <list>
 #include <string>
 #include <vector>
+#include <set>
 #include "Range.h"
 #include "SPoint3.h"
 #include "SBoundingBox3d.h"
@@ -66,6 +67,9 @@ class GEntity {
   // vertex arrays to draw the mesh efficiently
   VertexArray *va_lines, *va_triangles;
 
+  // Set of high-order elements fixed by "fast curving"
+  std::set<MElement*> curvedBLElements;
+
  public:
   // all known native model types
   enum ModelType {
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 62336f7c023f4785fe265d3c6ae43b04a7a0d99c..445e6dfdfe290c21ecc6913f8bf423915aa1efa8 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -926,8 +926,9 @@ double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt,
 
 
 
-void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
-                 int metaElType, std::vector<MVertex*> &baseVert,
+void curveColumn(const FastCurvingParameters &p, GEntity *ent,
+                 GEntity *bndEnt, int metaElType,
+                 std::vector<MVertex*> &baseVert,
                  const std::vector<MVertex*> &topPrimVert, MElement *aboveElt,
                  std::vector<MElement*> &blob, std::set<MVertex*> &movedVert,
                  DbgOutputMeta &dbgOut)
@@ -939,7 +940,7 @@ void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
 
   // If 2D P2 and allowed, modify base vertices to minimize distance between wall edge and CAD
   if ((p.optimizeGeometry) && (metaElType == TYPE_QUA) && (order == 2))
-    optimizeCADDist2DP2(geomEnt, baseVert);
+    optimizeCADDist2DP2(bndEnt, baseVert);
 
   // Create meta-element
   MetaEl metaElt(metaElType, order, baseVert, topPrimVert);
@@ -973,15 +974,17 @@ void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
 
   // Curve elements
   dbgOut.addMetaEl(metaElt);
-  for (int iEl = 0; iEl < blob.size(); iEl++)
+  for (int iEl = 0; iEl < blob.size(); iEl++) {
     curveElement(metaElt, movedVert, blob[iEl]);
+    ent->curvedBLElements.insert(blob[iEl]);
+  }
 }
 
 
 
 void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
-                         GEntity *bndEnt, MElement *bndElt,
-                         std::set<MVertex*> movedVert,
+                         GEntity *ent, GEntity *bndEnt,
+                         MElement *bndElt, std::set<MVertex*> movedVert,
                          const FastCurvingParameters &p,
                          DbgOutputMeta &dbgOut)
 {
@@ -1021,15 +1024,16 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
   dbgOutCol.addBlob(blob);
   dbgOutCol.write("col_KO", bndElt->getNum());
   if (aboveElt == 0) std::cout << "DBGTT: aboveElt = 0 for bnd. elt. " << bndElt->getNum() << std::endl;
-  curveColumn(p, bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
-              movedVert, dbgOut);
+  curveColumn(p, ent, bndEnt, metaElType, baseVert, topPrimVert,
+              aboveElt, blob, movedVert, dbgOut);
   dbgOutCol.write("col_OK", bndElt->getNum());
 }
 
 
 
 void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
-                      GEntity *bndEnt, const FastCurvingParameters &p)
+                      GEntity *ent, GEntity *bndEnt,
+                      const FastCurvingParameters &p)
 {
   // Build list of bnd. elements to consider
   std::list<MElement*> bndEl;
@@ -1054,7 +1058,8 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
   std::set<MVertex*> movedVert;
   for(std::list<MElement*>::iterator itBE = bndEl.begin();
       itBE != bndEl.end(); itBE++) // Loop over bnd. elements
-    curveMeshFromBndElt(ed2el, face2el, bndEnt, *itBE, movedVert, p, dbgOut);
+    curveMeshFromBndElt(ed2el, face2el, ent, bndEnt,
+                        *itBE, movedVert, p, dbgOut);
   dbgOut.write("meta-elements", bndEnt->tag());
 }
 
@@ -1118,7 +1123,7 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p,
       if ((bndType == GEntity::Line) || (bndType == GEntity::Plane)) continue;  // Skip if boundary is straight
       Msg::Info("Curving elements in entity %d for boundary entity %d...",
                 gEnt->tag(), bndEnt->tag());
-      curveMeshFromBnd(ed2el, face2el, bndEnt, p);
+      curveMeshFromBnd(ed2el, face2el, allGEnt[iEnt], bndEnt, p);
     }
   }
 
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
index b709bc08b663805334b5022c630d767da5589996..67fb56dc0cfe8b3f9d4acf0d275a9ae9934c85b6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
@@ -45,7 +45,7 @@ struct FastCurvingParameters {
 
   FastCurvingParameters () :
     dim(3), onlyVisible(true), optimizeGeometry(false),
-    curveOuterBL(OUTER_NOCURVE), maxNumLayers(100), maxRho(0.5),
+    curveOuterBL(OUTER_NOCURVE), maxNumLayers(100), maxRho(0.3),
     maxAngle(3.1415927*10./180.), maxAngleInner(3.1415927*30./180.)
   {
   }
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 96cd8c848934598c4d8f916a73873d54ddd2bf4d..8a648eca2fe1e1d74d734be7fcef79cac3fdace0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -727,7 +727,7 @@ struct HOPatchDefParameters : public MeshOptPatchDef
 private:
   double jacMin, jacMax;
   double distanceFactor;
-  bool optCAD;
+  bool optCAD, lockCurvedBLElts;
   double optCADDistMax, optCADWeight;
 };
 
@@ -751,6 +751,7 @@ HOPatchDefParameters::HOPatchDefParameters(const OptHomParameters &p)
   optCAD = p.optCAD;
   optCADDistMax = p.optCADDistMax;
   optCADWeight = p.optCADWeight;
+  lockCurvedBLElts = p.lockCurvedBLElts;
 }
 
 
@@ -791,6 +792,10 @@ int HOPatchDefParameters::inPatch(const SPoint3 &badBary,
                                   double limDist, MElement *el,
                                   GEntity* gEnt) const
 {
+  if (lockCurvedBLElts && (gEnt != 0)) {
+    const std::set<MElement*> &lockedElts = gEnt->curvedBLElements;
+    if (lockedElts.find(el) != lockedElts.end()) return -1;
+  }
   return testElInDist(badBary, limDist, el) ? 1 : 0;
 }
 
@@ -803,7 +808,7 @@ void HighOrderMeshOptimizerNew(std::vector<GEntity*> &entities, OptHomParameters
   par.dim = p.dim;
   par.onlyVisible = p.onlyVisible;
   par.fixBndNodes = p.fixBndNodes;
-  par.useGeomForPatches = false;
+  par.useGeomForPatches = p.lockCurvedBLElts;
   par.useGeomForOpt = false;
   par.useBoundaries = p.optCAD;
   HOPatchDefParameters patchDef(p);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index f7d12da60e71a578fd08d3e0e01beef46f1820d1..79608afe799c699ed23754dfcab25ec27e1c09e0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -57,6 +57,7 @@ struct OptHomParameters {
   double optCADWeight;//Weight
   double optCADDistMax;//Maximum allowed distance from the CAD
   double discrTolerance;
+  bool lockCurvedBLElts; // Do not include in optimization elements already fixed by "fast curving"
 
   // OUTPUT ------>
   int SUCCESS; // 0 --> success , 1 --> Not converged
@@ -70,7 +71,7 @@ struct OptHomParameters {
       distanceFactor(12), fixBndNodes(false), strategy(0), maxAdaptBlob(3),
       adaptBlobLayerFact(2.), adaptBlobDistFact(2.), optPrimSurfMesh(false),
       optCAD(false), optCADWeight(1000.), optCADDistMax(1.e22), discrTolerance(1.e-4),
-      SUCCESS(0), numBlobs(0), minJac(0.), maxJac(0.), CPU(0.)
+      lockCurvedBLElts(true), SUCCESS(0), numBlobs(0), minJac(0.), maxJac(0.), CPU(0.)
   {
   }
 };