diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index b2c14347c893a9f35031e1daba804a7e6b399ce9..918a4fbabd325ba22d726f1d08f73d17556ffc35 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -275,7 +275,6 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
             normals->set(0, 0, 0);
             normals->set(0, 1, 0);
             normals->set(0, 2, 1);
-            Msg::Error("Discrete face is 2D!");
           }
         }
         break;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 1121b97601cdb31d3151828320044b3d3d3f89d3..d91c48c9c16983658134a734ba93cf66b263f63b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -128,7 +128,7 @@ void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2enti
   SVector3 geoNorm(0.,0.,0.);
   std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]);
   GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second;
-  const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
+  bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
   for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
     const int &iV = _el2V[iEl][i];
     primNodesXYZ(i,0) = _xyz[iV].x();
@@ -145,6 +145,17 @@ void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2enti
     SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
     geoNorm = ((GFace*)ge)->normal(param);
   }
+  if (!hasGeoNorm && ge && ge->geomType() == GEntity::DiscreteSurface) {
+    SBoundingBox3d bb = ge->bounds();
+    // If we don't have the CAD, check if the mesh is 2D:
+    if (!bb.empty() && bb.max().z() - bb.min().z() == .0) {
+      hasGeoNorm = true;
+      geoNorm = SVector3(0, 0, 1);
+    }
+  }
+  else if (!ge) {
+    Msg::Warning("don't have the entity");
+  }
 
   fullMatrix<double> &elNorm = _scaledNormEl[iEl];
   elNorm.resize(1,3);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 83ff11f5b0dda7560f15bc1ad1103ca423a5c2f2..3d51a4f3ae3113fa0420af01161e4cfb04b6a9a0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -644,7 +644,12 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
     Msg::Info("Computing connectivity and bad elements for entity %d...",
               entity->tag());
     calcVertex2Elements(p.dim,entity,vertex2elements);
+
     if (p.optPrimSurfMesh) calcElement2Entity(entity,element2entity);
+    // Otherwise, still compute element2entity: Hack for using the geometry
+    // normal in the computation of the Jacobian when optimizing surface meshes.
+    // Warning: Accurate for planar surface but not really for curved surface...
+    else calcElement2Entity(entity,element2entity);
   }
 
   OptHomPeriodicity periodicity = OptHomPeriodicity(entities);