diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 279bad89447044c4b115753094f724d1b28c542d..be1f5cc07c59946c0c9599f6cf2115619df3ce1c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -861,7 +861,7 @@ StringXNumber GeometryOptions_Number[] = {
     "Display geometry surfaces?" },
   { F|O, "SurfaceNumbers" , opt_geometry_surfaces_num , 0. , 
     "Display surface numbers?" },
-  { F|O, "SurfaceType" , opt_geometry_surface_type , 0. , 
+  { F|O, "SurfaceType" , opt_geometry_surface_type , 1. , 
     "Display surfaces as wireframe (0) or solid (1)" },
 
   { F|O, "Tangents" , opt_geometry_tangents , 0. ,
diff --git a/Fltk/GUI_Projection.cpp b/Fltk/GUI_Projection.cpp
index 86d930e054e5fe3df3206b8218ce1ba01a483809..b94db3c61411a0e7bd6850e36d8721d2f7209c9b 100644
--- a/Fltk/GUI_Projection.cpp
+++ b/Fltk/GUI_Projection.cpp
@@ -603,7 +603,7 @@ void update_cb(Fl_Widget *w, void *data)
     for (int i = 0; i < ps->GetNumParameters(); i++)
       ps->SetParameter(i, p->parameters[i + 10]->value());
 
-    p->face->computeGraphicsRep(64, 64); // FIXME: hardcoded for now
+    p->face->buildSTLTriangulation();
 
     // project selected points and elements and update u,v display
     std::vector<double> u, v, dist;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 1eb708fba284425d0dac375b4ef04b121c5f76c6..ab292ff588c5d4595cf5eb9e6e6160ab2f83116a 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.44 2008-02-05 15:59:34 geuzaine Exp $
+// $Id: GFace.cpp,v 1.45 2008-02-05 23:16:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -541,10 +541,25 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p) const
   return SPoint2(U,V);
 }
 
-void GFace::computeGraphicsRep(int nu, int nv)
+struct graphics_point{
+  double xyz[3];
+  SVector3 n;
+};
+
+bool GFace::buildSTLTriangulation()
 {
-  _graphicsRep.resize(nu);
-  for(int i = 0; i < nu; i++) _graphicsRep[i].resize(nv);
+  // Build a simple triangulation for surfaces we know are not
+  // trimmed. Do nothing by default for complex surfaces (we might
+  // want to change this is the future)
+
+  if(geomType() != ParametricSurface && geomType() != ProjectionFace)
+    return false;
+
+  const int nu = 64, nv = 64;
+  graphics_point p[nu][nv];
+
+  if(va_geom_triangles) delete va_geom_triangles;
+  va_geom_triangles = new VertexArray(3, 2 * (nu - 1) * (nv - 1));
 
   Range<double> ubounds = parBounds(0);
   Range<double> vbounds = parBounds(1);
@@ -553,18 +568,35 @@ void GFace::computeGraphicsRep(int nu, int nv)
 
   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);
-      struct graphics_point gp;
-      GPoint p = point(u, v);
-      gp.xyz[0] = p.x();
-      gp.xyz[1] = p.y();
-      gp.xyz[2] = p.z();
-      SVector3 n = normal(SPoint2(u, v));
-      gp.n[0] = n.x();
-      gp.n[1] = n.y();
-      gp.n[2] = n.z();
-      _graphicsRep[i][j] = gp;
+      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));
+    }
+  }
+
+  // i,j+1 *---* i+1,j+1
+  //       | / |
+  //   i,j *---* i+1,j
+  unsigned int c = CTX.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);
     }
   }
+  va_geom_triangles->finalize();
+  return true;
 }
diff --git a/Geo/GFace.h b/Geo/GFace.h
index b985de3e0c2ee5b615b295dea4690a59d5986de1..4a4c4786b28b5f3b343f2a89a7c883c6061802c7 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -39,12 +39,6 @@ struct mean_plane
   double x, y, z;
 };
 
-struct graphics_point
-{
-  float xyz[3];
-  float n[3];
-};
-
 struct surface_params
 {
   double radius, radius2, height, cx, cy, cz;
@@ -55,11 +49,6 @@ class GRegion;
 // A model face. 
 class GFace : public GEntity 
 {
- private:
-  // a graphical representation for topologically simple surfaces: a
-  // 2D array of points/normals
-  std::vector<std::vector<graphics_point> > _graphicsRep;
-
  protected: 
   // edge loops, will replace what follows list of al the edges of the
   // face
@@ -73,9 +62,7 @@ class GFace : public GEntity
   // i.e. closed edge loops.  the first wire is the one that is the
   // outer contour of the face.
   void resolveWires();
-  // builds a STL triangulation and fill the vertex array
-  // va_geom_triangles (by default, do nothing)
-  virtual bool buildSTLTriangulation () { return false; }
+
  public:
   GFace(GModel *model, int tag);
   virtual ~GFace();
@@ -150,6 +137,10 @@ class GFace : public GEntity
   // Returns a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
+  // Builds a STL triangulation and fills the vertex array
+  // va_geom_triangles
+  virtual bool buildSTLTriangulation();
+
   // 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);
@@ -195,11 +186,8 @@ class GFace : public GEntity
   // of start/end points
   std::vector<SPoint3> cross;
 
-  // fill the graphics representation
-  void computeGraphicsRep(int nu, int nv);
-
-  // fill the graphics representation
-  std::vector<std::vector<graphics_point> > &getGraphicsRep(){ return _graphicsRep;}
+  // a vertex array containing an STL representation of the surface
+  VertexArray *va_geom_triangles;
 
   // a array for accessing the transfinite vertices using a pair of
   // indices
@@ -207,9 +195,6 @@ class GFace : public GEntity
 
   std::vector<MTriangle*> triangles;
   std::vector<MQuadrangle*> quadrangles;
-
-  // Vertex array to draw the geometry efficiently
-  VertexArray *va_geom_triangles;
 };
 
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index ca999a9f9081b110d39f73d459fa722d826ade5c..87512adffa29a2840b6d30dd6acfb5a8bfccba69 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.31 2008-02-05 18:58:04 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.32 2008-02-05 23:16:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -321,7 +321,7 @@ bool OCCFace::buildSTLTriangulation()
   unsigned int col[4] = {c, c, c, c};
   for (int j = 1; j <= ntriangles; j++){
     Poly_Triangle triangle = (triangulation->Triangles())(j);
-    triangle.Get(p1,p2,p3);
+    triangle.Get(p1, p2, p3);
     gp_Pnt2d x1 = (triangulation->UVNodes())(p1);
     gp_Pnt2d x2 = (triangulation->UVNodes())(p2);
     gp_Pnt2d x3 = (triangulation->UVNodes())(p3);
diff --git a/Geo/fourierFace.cpp b/Geo/fourierFace.cpp
index 405baf18c2595f7806e77903ae87f0765f373c13..ead0e16fd10ea852f16f1d999cec075e11cf5644 100644
--- a/Geo/fourierFace.cpp
+++ b/Geo/fourierFace.cpp
@@ -14,18 +14,7 @@ fourierFace::fourierFace(GModel *m, FM::TopoFace *face_, int tag,
     l_edges.push_back((*it));
     l_dirs.push_back(1);   
   }
-}
-
-fourierFace::fourierFace(GModel *m, FM::TopoFace *face_, int tag, 
-			 std::list<GEdge*> l_edges_, std::list<int> l_dirs_) 
-  : GFace(m,tag), face(face_)
-{
-  for (std::list<GEdge*>::iterator it = l_edges_.begin();
-       it != l_edges_.end(); it++)
-    l_edges.push_back((*it));  
-  for (std::list<int>::iterator it = l_dirs_.begin();
-       it != l_dirs_.end(); it++)
-    l_dirs.push_back((*it));  
+  buildSTLTriangulation();
 }
 
 Range<double> fourierFace::parBounds(int i) const
diff --git a/Geo/fourierFace.h b/Geo/fourierFace.h
index f4325b21b0c14dbb9989e030e8e20298775d6be7..538f7e55205c6196b0e96688620012b90368b200 100644
--- a/Geo/fourierFace.h
+++ b/Geo/fourierFace.h
@@ -15,8 +15,6 @@ class fourierFace : public GFace {
   FM::TopoFace *face;
  public:
   fourierFace(GModel *m, FM::TopoFace *face_, int tag, std::list<GEdge*> l_edges_);
-  fourierFace(GModel *m, FM::TopoFace *face_, int tag, std::list<GEdge*> l_edges_,
-	std::list<int> l_dirs_);
   virtual ~fourierFace() {}
   Range<double> parBounds(int i) const; 
   virtual GPoint point(double par1, double par2) const; 
diff --git a/Geo/fourierProjectionFace.cpp b/Geo/fourierProjectionFace.cpp
index 70dc6f316e258e319a4babac2d3bdd31f1fe97a0..2053b25ca5c3ea4053eda1d3f8e04fded590eda1 100644
--- a/Geo/fourierProjectionFace.cpp
+++ b/Geo/fourierProjectionFace.cpp
@@ -1,4 +1,5 @@
 #include "fourierProjectionFace.h"
+#include "VertexArray.h"
 
 #if defined(HAVE_FOURIER_MODEL)
 
@@ -6,7 +7,10 @@ fourierProjectionFace::fourierProjectionFace(GModel *m, int num)
   : GFace(m,num), ps_(0) {}
 
 fourierProjectionFace::fourierProjectionFace(GModel *m, int num, FM::ProjectionSurface* ps)
-  : GFace(m,num), ps_(ps) {}
+  : GFace(m,num), ps_(ps) 
+{
+  buildSTLTriangulation();
+}
 
 fourierProjectionFace::~fourierProjectionFace() {}
 
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 434ba30a5e4ffecbb52ced67deb627ab86f0925a..430a84fce36aad965cc85555443f973e6201a730 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.148 2008-02-05 21:26:42 geuzaine Exp $
+// $Id: Geom.cpp,v 1.149 2008-02-05 23:16:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -203,7 +203,7 @@ class drawGEdge {
 
 class drawGFace {
  private:
-  void _drawNonPlaneGFace(GFace *f)
+  void _drawParametricGFace(GFace *f)
   {
     if(CTX.geom.surfaces && !f->va_geom_triangles) {
       glEnable(GL_LINE_STIPPLE);
@@ -254,64 +254,6 @@ class drawGFace {
 		  p.x(), p.y(), p.z(), n[0], n[1], n[2], CTX.geom.light);
     }
   }
-  
-  void _drawParametricGFace(GFace *f)
-  {
-    std::vector<std::vector<graphics_point> > &gr(f->getGraphicsRep());
-
-    const unsigned int N = 64;
-
-    // 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
-   
-    static bool busy = false;
-    if(gr.size() != N && !busy) {
-      busy = true; 
-      f->computeGraphicsRep(N, N);
-      busy = false;
-    }
-
-    if(gr.size() != N) return;
-
-    if(f->geomType() == GEntity::ProjectionFace)
-      glColor4ubv((GLubyte *) & CTX.color.geom.projection);
-
-    if(CTX.geom.surface_type > 0)
-      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
-    else
-      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
-
-    if(CTX.geom.light) glEnable(GL_LIGHTING);
-    glBegin(GL_QUADS);
-    for(unsigned int i = 1; i < gr.size(); i++){
-      for(unsigned int j = 1; j < gr[0].size(); j++){
-        glNormal3fv(gr[i - 1][j - 1].n);
-        glVertex3fv(gr[i - 1][j - 1].xyz);
-        glNormal3fv(gr[i    ][j - 1].n);
-        glVertex3fv(gr[i    ][j - 1].xyz);
-        glNormal3fv(gr[i    ][j    ].n);
-        glVertex3fv(gr[i    ][j    ].xyz);
-        glNormal3fv(gr[i - 1][j    ].n);
-        glVertex3fv(gr[i - 1][j    ].xyz);
-      }
-    }
-    glEnd();
-    glDisable(GL_LIGHTING);
-    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
-
-    if(CTX.geom.normals) {
-      GPoint p = f->point(0.5, 0.5);
-      SVector3 n = f->normal(SPoint2(0.5, 0.5));
-      for(int i = 0; i < 3; i++)
-	n[i] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[i];
-      glColor4ubv((GLubyte *) & CTX.color.geom.normals);
-      Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
-		  CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
-		  p.x(), p.y(), p.z(), n[0], n[1], n[2], CTX.geom.light);
-    }
-  }
-
   void _drawPlaneGFace(GFace *f)
   {
     // We create data here and the routine is not designed to be
@@ -442,9 +384,14 @@ class drawGFace {
       glEnableClientState(GL_COLOR_ARRAY);
     }
     if(CTX.polygon_offset) glEnable(GL_POLYGON_OFFSET_FILL);
+    if(CTX.geom.surface_type > 0)
+      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+    else
+      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
     glDrawArrays(GL_TRIANGLES, 0, va->getNumVertices());
     glDisable(GL_POLYGON_OFFSET_FILL);
     glDisable(GL_LIGHTING);
+    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
     glDisableClientState(GL_VERTEX_ARRAY);
     glDisableClientState(GL_NORMAL_ARRAY);
     glDisableClientState(GL_COLOR_ARRAY);
@@ -475,16 +422,16 @@ public :
     }
 
     if(CTX.geom.surfaces && f->va_geom_triangles)
-      _drawVertexArray(f->va_geom_triangles, CTX.geom.light, f->getSelection(), 
-		       CTX.color.geom.selection);
+      _drawVertexArray
+	(f->va_geom_triangles, CTX.geom.light, 
+	 (f->geomType() == GEntity::ProjectionFace) ? true : f->getSelection(), 
+	 (f->geomType() == GEntity::ProjectionFace) ? CTX.color.geom.projection : 
+	 CTX.color.geom.selection);
     
     if(f->geomType() == GEntity::Plane)
       _drawPlaneGFace(f);
-    else if(f->geomType() == GEntity::ProjectionFace ||
-	    f->geomType() == GEntity::ParametricSurface)
-      _drawParametricGFace(f);
     else
-      _drawNonPlaneGFace(f);
+      _drawParametricGFace(f);
     
     if(CTX.draw_bbox) drawBBox(f);
     
@@ -596,12 +543,13 @@ void Draw_Geom()
     gl2psLineWidth(CTX.line_width * CTX.print.eps_line_width_factor);
     if(!CTX.axes_auto_position){
       Draw_Axes(CTX.axes, CTX.axes_tics, CTX.axes_format, CTX.axes_label, 
-		CTX.axes_position,CTX.axes_mikado);
+		CTX.axes_position, CTX.axes_mikado);
     }
     else if(geometryExists){
       double bb[6] = {CTX.min[0], CTX.max[0], CTX.min[1], 
 		      CTX.max[1], CTX.min[2], CTX.max[2]};
-      Draw_Axes(CTX.axes, CTX.axes_tics, CTX.axes_format, CTX.axes_label, bb,CTX.axes_mikado);
+      Draw_Axes(CTX.axes, CTX.axes_tics, CTX.axes_format, CTX.axes_label, 
+		bb, CTX.axes_mikado);
     }
   }