diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index aaf2ec1072b769e086e3de2d5785c7f268c59fbd..c23a0de91c903815cc503e6902e7efa73150a764 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -412,7 +412,7 @@ void GFace::computeMeanPlane()
   std::vector<SPoint3> pts;
 
   if(geomType() == Plane) {
-    // Special case for planar surfaces: we first try to compute the
+    // Special case for planar CAD surfaces: we first try to compute the
     // parametrization in a way that is more robust with respect to
     // perturbations of the boundary. This is crucial for sensitivity analyses
     // and GFace::relocateMeshVertices(): after perturbation of the boundary, we
@@ -421,6 +421,11 @@ void GFace::computeMeanPlane()
     std::list<GEdge*> edg = edges();
     for(std::list<GEdge*>::const_iterator ite = edg.begin(); ite != edg.end(); ite++){
       const GEdge *e = *ite;
+      if(e->geomType() == GEntity::DiscreteCurve ||
+         e->geomType() == GEntity::BoundaryLayerCurve){
+        pts.clear();
+        break;
+      }
       Range<double> b = e->parBounds(0);
       GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
       pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));