diff --git a/Common/Context.h b/Common/Context.h
index 601aecd38377bf241c9c2d3a7b55a3e4cdd885e3..2c21f2c4ee0c3084ec16ccfd0ed7ff22986ec517 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -42,7 +42,7 @@ struct contextMeshOptions {
   int order, secondOrderLinear, secondOrderIncomplete;
   int secondOrderExperimental, meshOnlyVisible;
   int minCircPoints, minCurvPoints;
-  int hoOptimize, hoNLayers;
+  int hoOptimize, hoNLayers, hoOptPrimSurfMesh;
   double hoThresholdMin, hoThresholdMax, hoPoissonRatio;
   int saveAll, saveTri, saveGroupsOfNodes, binary, bdfFieldFormat, saveParametric;
   int smoothNormals, reverseAllNormals, zoneDefinition, clip;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index fb60147fbbbced8045649708ca0e9d413954ce59..35058f056ad4ffb32b521e71951e7e0592b37483 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -980,6 +980,8 @@ StringXNumber MeshOptions_Number[] = {
     "Minimum threshold for high order element optimization"},
   { F|O, "HighOrderThresholdMax", opt_mesh_ho_threshold_max, 2.0,
     "Maximum threshold for high order element optimization"},
+  { F|O, "HighOrderOptPrimSurfMesh", opt_mesh_ho_opt_prim_surf_mesh, 0,
+    "Try to fix flipped surface mesh elements in high-order optimizer"},
 
   { F|O, "LabelSampling" , opt_mesh_label_sampling , 1. ,
     "Label sampling rate (display one label every `LabelSampling' elements)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 796e342d5de77f54effa730343c76870b7cee35b..f980418be90295365c610069bb8f6b6ff2f7402e 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5638,6 +5638,13 @@ double opt_mesh_ho_threshold_max(OPT_ARGS_NUM)
   return CTX::instance()->mesh.hoThresholdMax;
 }
 
+double opt_mesh_ho_opt_prim_surf_mesh(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX::instance()->mesh.hoOptPrimSurfMesh = (int)val;
+  return CTX::instance()->mesh.hoOptPrimSurfMesh;
+}
+
 double opt_mesh_ho_poisson(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET) {
diff --git a/Common/Options.h b/Common/Options.h
index e2b3dd6cbd3b2769ca9ac9cc540b3db8dc69982a..1007a85cad4e21742c3068c92b12b292fb49e1f4 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -453,6 +453,7 @@ double opt_mesh_ho_nlayers(OPT_ARGS_NUM);
 double opt_mesh_ho_threshold_min(OPT_ARGS_NUM);
 double opt_mesh_ho_threshold_max(OPT_ARGS_NUM);
 double opt_mesh_ho_poisson(OPT_ARGS_NUM);
+double opt_mesh_ho_opt_prim_surf_mesh(OPT_ARGS_NUM);
 double opt_mesh_second_order_experimental(OPT_ARGS_NUM);
 double opt_mesh_second_order_linear(OPT_ARGS_NUM);
 double opt_mesh_second_order_incomplete(OPT_ARGS_NUM);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 5eb70430af6039ac8d2319304bca8d89eba3885d..1505d499969a6d7d0f5e584506b5bdf59de198a7 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -808,7 +808,7 @@ void GenerateMesh(GModel *m, int ask)
       p.BARRIER_MIN = CTX::instance()->mesh.hoThresholdMin;
       p.BARRIER_MAX = CTX::instance()->mesh.hoThresholdMax;
       p.dim = GModel::current()->getDim();
-      //p.optPrimSurfMesh = true;
+      p.optPrimSurfMesh = CTX::instance()->mesh.hoOptPrimSurfMesh;
       HighOrderMeshOptimizer(GModel::current(), p);
     }
 #else