diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index b5d96a90ea175d9e73dffb4955af1d6fcb08ec9c..ba883064818ba0f0d7f7465d938dd3d2e4e2b272 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -178,12 +178,9 @@ SOrientedBoundingBox GFace::getOBB()
       }
     } 
     else if(buildSTLTriangulation()) {
-      int N = va_geom_triangles->getNumVertices();
-      for (int i = 0; i < N; i++) {
-        SPoint3 p((va_geom_triangles->getVertexArray(3 * i))[0],
-                  (va_geom_triangles->getVertexArray(3 * i))[1],
-                  (va_geom_triangles->getVertexArray(3 * i))[2]);
-        vertices.push_back(p);
+      for (unsigned int i = 0; i < stl_vertices.size(); i++){
+        GPoint p = point(stl_vertices[i]);
+        vertices.push_back(SPoint3(p.x(), p.y(), p.z()));
       }
     }
     else {
@@ -789,10 +786,16 @@ SVector3 GFace::normal(const SPoint2 &param) const
   return n;
 }
 
-bool GFace::buildRepresentationCross()
+bool GFace::buildRepresentationCross(bool force)
 {
+  if(cross.size()){
+    if(force)
+      cross.clear();
+    else
+      return true;
+  }
+
   if(geomType() != Plane){
-    // don't try again
     cross.clear();
     cross.push_back(SPoint3(0., 0., 0.));
     return false;
@@ -804,7 +807,6 @@ bool GFace::buildRepresentationCross()
     GEdge *ge = *it;
     if(ge->geomType() == GEntity::DiscreteCurve ||
        ge->geomType() == GEntity::BoundaryLayerCurve){
-      // don't try again
       cross.clear();
       cross.push_back(SPoint3(0., 0., 0.));
       return false;
@@ -853,8 +855,8 @@ bool GFace::buildRepresentationCross()
     }
     if(end_line) cross.push_back(pt_last_inside);
   }
-  // if we couldn't determine a cross, add a dummy point so that
-  // we won't try again
+  // if we couldn't determine a cross, add a dummy point so that we
+  // won't try again unless we force the recomputation
   if(!cross.size()){
     cross.push_back(SPoint3(0., 0., 0.));
     return false;
@@ -862,18 +864,52 @@ bool GFace::buildRepresentationCross()
   return true;
 }
 
-struct graphics_point{
-  double xyz[3];
-  SVector3 n;
-};
-
 bool GFace::buildSTLTriangulation(bool force)
 {
+  if(stl_triangles.size()){
+    if(force){
+      stl_vertices.clear();
+      stl_triangles.clear();
+    }
+    else
+      return true;
+  }
+
   // Build a simple triangulation for surfaces which we know are not
   // trimmed
-  if(geomType() != ParametricSurface && geomType() != ProjectionFace)
-    return false;
+  if(geomType() == ParametricSurface || geomType() == ProjectionFace){
+    const int nu = 64, nv = 64;
+    Range<double> ubounds = parBounds(0);
+    Range<double> vbounds = parBounds(1);
+    double umin = ubounds.low(), umax = ubounds.high();
+    double vmin = vbounds.low(), vmax = vbounds.high();
+    for(int i = 0; i < nu; i++){
+      for(int j = 0; j < nv; j++){
+        double u = umin + (double)i / (double)(nu - 1) * (umax - umin);
+        double v = vmin + (double)j / (double)(nv - 1) * (vmax - vmin);
+        stl_vertices.push_back(SPoint2(u, v));
+      }
+    }
+    for(int i = 0; i < nu - 1; i++){
+      for(int j = 0; j < nv - 1; j++){
+        stl_triangles.push_back(i * nv + j);
+        stl_triangles.push_back((i + 1) * nv + j);
+        stl_triangles.push_back((i + 1) * nv + j + 1);
+        stl_triangles.push_back(i * nv + j);
+        stl_triangles.push_back((i + 1) * nv + j + 1);
+        stl_triangles.push_back(i * nv + j + 1);
+      }
+    }
+    return true;
+  }
+  
+  // build STL for general surfaces here
+
+  return false;
+}
 
+bool GFace::fillVertexArray(bool force)
+{
   if(va_geom_triangles){
     if(force)
       delete va_geom_triangles;
@@ -881,48 +917,26 @@ bool GFace::buildSTLTriangulation(bool force)
       return true;
   }
 
-  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);
-  double umin = ubounds.low(), umax = ubounds.high();
-  double vmin = vbounds.low(), vmax = vbounds.high();
-
-  for(int i = 0; i < nu; i++){
-    for(int j = 0; j < nv; j++){
-      double u = umin + (double)i / (double)(nu - 1) * (umax - umin);
-      double v = vmin + (double)j / (double)(nv - 1) * (vmax - vmin);
-      GPoint gp = point(u, v);
-      p[i][j].xyz[0] = gp.x();
-      p[i][j].xyz[1] = gp.y();
-      p[i][j].xyz[2] = gp.z();
-      p[i][j].n = normal(SPoint2(u, v));
-    }
-  }
+  if(!buildSTLTriangulation(force)) return false;
+  if(stl_triangles.empty()) return false;
 
-  // i,j+1 *---* i+1,j+1
-  //       | / |
-  //   i,j *---* i+1,j
+  va_geom_triangles = new VertexArray(3, stl_triangles.size() / 3);
   unsigned int c = CTX::instance()->color.geom.surface;
   unsigned int col[4] = {c, c, c, c};
-  for(int i = 0; i < nu - 1; i++){
-    for(int j = 0; j < nv - 1; j++){
-      double x1[3] = {p[i][j].xyz[0], p[i + 1][j].xyz[0], p[i + 1][j + 1].xyz[0]};
-      double y1[3] = {p[i][j].xyz[1], p[i + 1][j].xyz[1], p[i + 1][j + 1].xyz[1]};
-      double z1[3] = {p[i][j].xyz[2], p[i + 1][j].xyz[2], p[i + 1][j + 1].xyz[2]};
-      SVector3 n1[3] = {p[i][j].n, p[i + 1][j].n, p[i + 1][j + 1].n};
-      va_geom_triangles->add(x1, y1, z1, n1, col);
-      double x2[3] = {p[i][j].xyz[0], p[i + 1][j + 1].xyz[0], p[i][j + 1].xyz[0]};
-      double y2[3] = {p[i][j].xyz[1], p[i + 1][j + 1].xyz[1], p[i][j + 1].xyz[1]};
-      double z2[3] = {p[i][j].xyz[2], p[i + 1][j + 1].xyz[2], p[i][j + 1].xyz[2]};
-      SVector3 n2[3] = {p[i][j].n, p[i + 1][j + 1].n, p[i][j + 1].n};
-      va_geom_triangles->add(x2, y2, z2, n2, col);
-    }
+  for (unsigned int i = 0; i < stl_triangles.size(); i += 3){
+    SPoint2 &p1(stl_vertices[stl_triangles[i]]);
+    SPoint2 &p2(stl_vertices[stl_triangles[i + 1]]);
+    SPoint2 &p3(stl_vertices[stl_triangles[i + 2]]);
+    GPoint gp1 = point(p1);
+    GPoint gp2 = point(p2);
+    GPoint gp3 = point(p3);
+    double x[3] = {gp1.x(), gp2.x(), gp3.x()};
+    double y[3] = {gp1.y(), gp2.y(), gp3.y()};
+    double z[3] = {gp1.z(), gp2.z(), gp3.z()};
+    SVector3 n[3] = {normal(p1), normal(p2), normal(p3)};
+    va_geom_triangles->add(x, y, z, n, col);
   }
   va_geom_triangles->finalize();
-  return true;
 }
 
 // by default we assume that straight lines are geodesics
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 5bef52a9bf6bf2c20825e6c9a0eaed047122e915..35404222159af861122b4a1e310068ae2927717e 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -191,13 +191,15 @@ class GFace : public GEntity
   virtual void writeGEO(FILE *fp);
 
   // fill the crude representation cross
-  virtual bool buildRepresentationCross();
+  virtual bool buildRepresentationCross(bool force=false);
 
-  // build a STL triangulation and fills the vertex array
-  // va_geom_triangles (or do nothing if it already exists, unless
-  // force=true)
+  // build an STL triangulation (or do nothing if it already exists,
+  // unless force=true)
   virtual bool buildSTLTriangulation(bool force=false);
 
+  // fill the vertex array using an STL triangulation
+  bool fillVertexArray(bool force=false);
+
   // recompute the mean plane of the surface from a list of points
   void computeMeanPlane(const std::vector<MVertex*> &points);
   void computeMeanPlane(const std::vector<SPoint3> &points);
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 6034b726137f30362a0d6cdf42bb29ed3c5fc7b7..cde7ea41dd66dfa48986be9a213335a0ab7add2c 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -94,7 +94,6 @@ class GFaceCompound : public GFace {
   ModelType getNativeType() const { return GmshModel; }
   void * getNativePtr() const { return 0; }
   virtual SPoint2 getCoordinates(MVertex *v) const;
-  virtual bool buildRepresentationCross(){ return false; }
   virtual double curvatureMax(const SPoint2 &param) const;
   virtual int genusGeom ();
  private:
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 1c9bb6796e91e956e04aade59519084f3ef69aba..72844a4012c4627c1a92ed442b62bc0b13fcfb99 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -152,14 +152,11 @@ SOrientedBoundingBox GRegion::getOBB()
         }
       } 
       else if ((*b_face)->buildSTLTriangulation()) {
-        int N = (*b_face)->va_geom_triangles->getNumVertices();
-        for(int i = 0; i < N; i++) {
-          SPoint3 p(((*b_face)->va_geom_triangles->getVertexArray(3*i))[0],
-                    ((*b_face)->va_geom_triangles->getVertexArray(3*i))[1],
-                    ((*b_face)->va_geom_triangles->getVertexArray(3*i))[2]);
-          vertices.push_back(p);          
-        }
-      } 
+        for (unsigned int i = 0; i < (*b_face)->stl_vertices.size(); i++){
+          GPoint p = (*b_face)->point((*b_face)->stl_vertices[i]);
+          vertices.push_back(SPoint3(p.x(), p.y(), p.z()));
+        } 
+      }
       else {
         int N = 10;
         std::list<GEdge*> b_edges = (*b_face)->edges();
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 34dda7daf68eab635e69cd06d9366e3d491e93d9..6bf884ad87f95815f9c09ca4636b838ad9b50efe 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -11,7 +11,6 @@
 #include "OCCEdge.h"
 #include "OCCFace.h"
 #include "Numeric.h"
-#include "VertexArray.h"
 #include "Context.h"
 
 #if defined(HAVE_OCC)
@@ -84,7 +83,6 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
   umax += fabs(du) / 100.0;
   vmax += fabs(dv) / 100.0;
   occface = BRep_Tool::Surface(s);
-  if(!CTX::instance()->batch) buildSTLTriangulation();
 }
 
 Range<double> OCCFace::parBounds(int i) const
@@ -305,16 +303,15 @@ surface_params OCCFace::getSurfaceParams() const
 
 bool OCCFace::buildSTLTriangulation(bool force)
 {
-  if(va_geom_triangles){
-    if(force)
-      delete va_geom_triangles;
+  if(stl_triangles.size()){
+    if(force){
+      stl_vertices.clear();
+      stl_triangles.clear();
+    }
     else
       return true;
   }
 
-  stl_vertices.clear();
-  stl_triangles.clear();
-  
   Bnd_Box aBox;
   BRepBndLib::Add(s, aBox);
   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
@@ -355,23 +352,6 @@ bool OCCFace::buildSTLTriangulation(bool force)
     stl_triangles.push_back(p3 - 1);
   }
 
-  va_geom_triangles = new VertexArray(3, stl_triangles.size() / 3);
-  unsigned int c = CTX::instance()->color.geom.surface;
-  unsigned int col[4] = {c, c, c, c};
-  for (unsigned int i = 0; i < stl_triangles.size(); i += 3){
-    SPoint2 &p1(stl_vertices[stl_triangles[i]]);
-    SPoint2 &p2(stl_vertices[stl_triangles[i + 1]]);
-    SPoint2 &p3(stl_vertices[stl_triangles[i + 2]]);
-    GPoint gp1 = GFace::point(p1);
-    GPoint gp2 = GFace::point(p2);
-    GPoint gp3 = GFace::point(p3);
-    double x[3] = {gp1.x(), gp2.x(), gp3.x()};
-    double y[3] = {gp1.y(), gp2.y(), gp3.y()};
-    double z[3] = {gp1.z(), gp2.z(), gp3.z()};
-    SVector3 n[3] = {normal(p1), normal(p2), normal(p3)};
-    va_geom_triangles->add(x, y, z, n, col);
-  }
-  va_geom_triangles->finalize();
   return true;
 }
 
diff --git a/Graphics/drawGeom.cpp b/Graphics/drawGeom.cpp
index 07957bd1c398c9b0b153fe650f5003695fc132c0..277604f6a6f9346a99325a707be8f868091e06a6 100644
--- a/Graphics/drawGeom.cpp
+++ b/Graphics/drawGeom.cpp
@@ -236,7 +236,18 @@ class drawGFace {
   }
   void _drawParametricGFace(GFace *f)
   {
-    if (f->geomType() == GEntity::CompoundSurface)return;
+    if (f->geomType() == GEntity::CompoundSurface) return;
+
+    if(CTX::instance()->geom.surfaceType > 0){
+      // avoid reentrant calls
+      static bool busy = false;
+      if(!busy) {
+        busy = true;
+        f->fillVertexArray(f->geomType() == GEntity::ProjectionFace);
+        busy = false;
+      }
+    }
+
     Range<double> ubounds = f->parBounds(0);
     Range<double> vbounds = f->parBounds(1);
     const double uav = 0.5 * (ubounds.high() + ubounds.low());
@@ -307,20 +318,27 @@ class drawGFace {
   }
   void _drawPlaneGFace(GFace *f)
   {
+    if(CTX::instance()->geom.surfaceType > 0){
+      // avoid reentrant calls
+      static bool busy = false;
+      if(!busy) {
+        busy = true; 
+        f->fillVertexArray();
+        busy = false;
+      }
+    }
+
     if(!CTX::instance()->geom.surfaceType || !f->va_geom_triangles ||
        CTX::instance()->geom.surfacesNum || CTX::instance()->geom.normals){
-      // We create data here and the routine is not designed to be
-      // reentrant, so we must lock it to avoid race conditions when
-      // redraw events are fired in rapid succession
+      // avoid reentrant calls
       static bool busy = false;
-      if(f->cross.empty() && !busy) {
+      if(!busy) {
         busy = true; 
         f->buildRepresentationCross();
         busy = false;
       }
     }
 
-    //FIXME: cleanup buildGraphicsRep
     if(CTX::instance()->geom.surfaces) {
       if(CTX::instance()->geom.surfaceType > 0 && f->va_geom_triangles){
         _drawVertexArray(f->va_geom_triangles, CTX::instance()->geom.light, 
diff --git a/Solver/SElement.h b/Solver/SElement.h
index 518800fefac80ef95aa0130cefc28d50b477ec8c..cb59a49e38537268fd9b6cfa0ad93829d273418f 100644
--- a/Solver/SElement.h
+++ b/Solver/SElement.h
@@ -11,6 +11,38 @@
 
 // A solver element.
 
+/*
+class SFunctionSpace{
+ public:
+  
+}
+
+class SFunctionSpaceXFEM : public SFunctionSpace{
+  static simpleFunction<double> *_enrichement_s, *_enrichement_t;
+  // gradient of functions (possibly enriched)
+  void nodalFunctions ( double u, double v, double w, double s[],
+			simpleFunction<double> *_enrichement);   
+  void gradNodalFunctions ( double u, double v, double w, double invjac[3][3],
+			    double grad[][3], simpleFunction<double> *_enrichment);   
+ public:
+  static void setShapeEnrichement(simpleFunction<double>*f){_enrichement_s=f;}
+  static void setTestEnrichement(simpleFunction<double>*f){_enrichement_t=f;}
+  static const simpleFunction<double>* getShapeEnrichement(){return _enrichement_s;}
+  static const simpleFunction<double>* getTestEnrichement(){return _enrichement_t;}
+  int getNumNodalShapeFunctions() const;
+  inline MVertex *getVertex(int i) const {return _e->getParent() ? _e->getParent()->getVertex(i) : _e->getVertex(i) ; }
+  int getNumNodalTestFunctions() const;
+  void nodalShapeFunctions ( double u, double v, double w, double s[]); 
+  void gradNodalShapeFunctions ( double u, double v, double w, double invjac[3][3],
+				 double grad[][3]); 
+  void nodalTestFunctions ( double u, double v, double w, double s[]); 
+  void gradNodalTestFunctions ( double u, double v, double w, double invjac[3][3],
+			        double grad[][3]);   
+  
+}
+
+*/
+
 class SElement
 {
  private: