diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 47770f1c32b2746e6892d18152a6aa90d7e550dd..deeb7fc2c6b4e62b6c6428c0f49ee80117f43ecf 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -867,18 +867,23 @@ struct graphics_point{
   SVector3 n;
 };
 
-bool GFace::buildSTLTriangulation()
+bool GFace::buildSTLTriangulation(bool force)
 {
   // Build a simple triangulation for surfaces which we know are not
   // trimmed
   if(geomType() != ParametricSurface && geomType() != ProjectionFace)
     return false;
 
-  const int nu = 64, nv = 64;
-  graphics_point p[nu][nv];
+  if(va_geom_triangles){
+    if(force)
+      delete va_geom_triangles;
+    else
+      return true;
+  }
 
-  if(va_geom_triangles) delete va_geom_triangles;
+  const int nu = 64, nv = 64;
   va_geom_triangles = new VertexArray(3, 2 * (nu - 1) * (nv - 1));
+  graphics_point p[nu][nv];
 
   Range<double> ubounds = parBounds(0);
   Range<double> vbounds = parBounds(1);
@@ -986,9 +991,31 @@ int GFace::genusGeom()
     if ((*it)->isSeam(this)){
       nSeams++;
       std::set<GEdge*>::iterator it2 = single_seams.find(*it);
-      if (it2 != single_seams.end())single_seams.erase(it2);
+      if (it2 != single_seams.end()) single_seams.erase(it2);
       else single_seams.insert(*it);
     }
   }
   return nSeams - single_seams.size();
 }
+
+void GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+                           std::vector<SVector3> *normals)
+{
+  if(!buildSTLTriangulation()){
+    Msg::Error("No STL triangulation available to fill point cloud");
+    return;
+  }
+
+  bool computePoints = points ? true : false;
+  bool computeNormals = normals ? true : false;
+  if(!computePoints && !computeNormals) return;
+  /*
+  int N = va_geom_triangles->getNumVertices();
+  for (int i = 0; i < N; i += 3) {
+    SPoint3 p((va_geom_triangles->getVertexArray(3 * i))[0],
+              (va_geom_triangles->getVertexArray(3 * i))[1],
+              (va_geom_triangles->getVertexArray(3 * i))[2]);
+  }
+  */
+  
+}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index ac78c14bdb69a681b0f4b2bf2d2110b9cffbc4b1..f6b033a93b7426e45c7dbc4e2ba14ded3ce6f16d 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -194,8 +194,9 @@ class GFace : public GEntity
   virtual bool buildRepresentationCross();
 
   // build a STL triangulation and fills the vertex array
-  // va_geom_triangles
-  virtual bool buildSTLTriangulation();
+  // va_geom_triangles (or do nothing if it already exists, unless
+  // force=true)
+  virtual bool buildSTLTriangulation(bool force=false);
 
   // recompute the mean plane of the surface from a list of points
   void computeMeanPlane(const std::vector<MVertex*> &points);
@@ -229,6 +230,11 @@ class GFace : public GEntity
   void setCompound(GFaceCompound *gfc) { compound = gfc; }
   GFaceCompound *getCompound() const { return compound; }
 
+  // add points (and optionalluy normals) in vectors so that two
+  // points are at most maxDist apart
+  void fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+                      std::vector<SVector3> *normals);
+
   struct {
     // do we recombine the triangles of the mesh?
     int recombine;
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 1e697f5160438f582659ee1caa7e536d10a1a927..96cd08ffd6523c2508ced27e932f258affcf08ca 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -303,13 +303,15 @@ surface_params OCCFace::getSurfaceParams() const
   return p;
 }
 
-bool OCCFace::buildSTLTriangulation()
+bool OCCFace::buildSTLTriangulation(bool force)
 {
   if(va_geom_triangles){
-    delete va_geom_triangles;
-    va_geom_triangles = 0;
+    if(force)
+      delete va_geom_triangles;
+    else
+      return true;
   }
-  
+
   TopLoc_Location loc;
   int p1, p2, p3;
   Bnd_Box aBox;
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index e84bed22b3f341d96717437985e155e5b84e5815..863d184ea2fe27f6f196e8bc17df096ca1172713 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -21,7 +21,7 @@ class OCCFace : public GFace {
   Handle(Geom_Surface) occface;
   double umin, umax, vmin, vmax;
   bool _periodic[2];
-  bool buildSTLTriangulation();
+  bool buildSTLTriangulation(bool force=false);
  public:
   OCCFace(GModel *m, TopoDS_Face s, int num, TopTools_IndexedMapOfShape &emap);
   virtual ~OCCFace(){}