From 6ec04bdb4dc563ec72533976f27bc569206e82d9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 11 Mar 2007 20:19:06 +0000
Subject: [PATCH] - preliminary work to interface more parser functions with
 general   GEntities

- we can now export simple GModels as.geo files (points, straight line
  segments, plane surfaces, ruled surfaces and volumes are treated
  correctly; all non-straight-line curves are exported as splines; all
  non-plane/ruled surfaces are not exported)

- tentative fix for the unv export problem reported on the mailing
  list (fprintf() buggy on some machines??)

- new IO (both read and write) for structured grids using the NASA
  Plot3D format (for Olivier)
---
 Common/CommandLine.cpp          |   7 +-
 Fltk/Callbacks.cpp              |  24 +-
 Geo/GFace.h                     |  13 +-
 Geo/GModel.h                    |  43 ++-
 Geo/GModelIO_Fourier.cpp        |  26 +-
 Geo/GModelIO_Geo.cpp            | 242 ++++++------
 Geo/GModelIO_Mesh.cpp           | 177 ++++++++-
 Geo/GRegion.h                   |  10 +-
 Geo/Geo.cpp                     |   8 +-
 Geo/Geo.h                       |  71 ++--
 Geo/GeoInterpolation.cpp        |  57 +--
 Geo/MVertex.cpp                 |   6 +-
 Graphics/Mesh.cpp               |   4 +-
 Makefile                        |   4 +-
 Mesh/meshGFace.cpp              |   3 +-
 Mesh/meshGFaceExtruded.cpp      |  11 +-
 Mesh/meshGFaceTransfinite.cpp   |  85 +++--
 Mesh/meshGRegion.cpp            |   3 +-
 Mesh/meshGRegionExtruded.cpp    |   5 +-
 Mesh/meshGRegionTransfinite.cpp | 248 ++++++------
 Parser/CreateFile.cpp           |   8 +-
 Parser/Gmsh.tab.cpp             | 653 +++++++++++++++++---------------
 Parser/Gmsh.y                   |  73 ++--
 Parser/Gmsh.yy.cpp              |   4 +-
 Parser/OpenFile.cpp             |   5 +-
 benchmarks/3d/sph.geo           |   8 +-
 benchmarks/3d/transfinite.geo   |   7 +-
 doc/TODO                        |  14 +-
 doc/VERSIONS                    |  10 +-
 doc/gmsh.1                      |   4 +-
 doc/texinfo/command_line.texi   |   2 +-
 31 files changed, 1072 insertions(+), 763 deletions(-)

diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 5054b7d693..39bc1bb08e 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.94 2007-01-30 08:56:36 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.95 2007-03-11 20:18:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -81,7 +81,8 @@ void Print_Usage(char *name){
   Msg(DIRECT, "  -1, -2, -3            Perform 1D, 2D or 3D mesh generation, then exit");
   Msg(DIRECT, "  -saveall              Save all elements (discard physical group definitions)");
   Msg(DIRECT, "  -o file               Specify mesh output file name");
-  Msg(DIRECT, "  -format string        Set output mesh format (msh, unv, vrml, stl, mesh, bdf, cgns, med)");
+  Msg(DIRECT, "  -format string        Set output mesh format (msh, msh1, msh2, unv, vrml, stl, mesh");
+  Msg(DIRECT, "                        bdf, p3d, cgns, med)");
   Msg(DIRECT, "  -algo string          Select mesh algorithm (iso, netgen, tetgen)");
   Msg(DIRECT, "  -smooth int           Set number of mesh smoothing steps");
   Msg(DIRECT, "  -optimize             Optimize quality of tetrahedral elements");
@@ -395,6 +396,8 @@ void Get_Options(int argc, char *argv[])
             CTX.mesh.format = FORMAT_MESH;
           else if(!strcmp(argv[i], "bdf"))
             CTX.mesh.format = FORMAT_BDF;
+          else if(!strcmp(argv[i], "p3d"))
+            CTX.mesh.format = FORMAT_P3D;
           else if(!strcmp(argv[i], "cgns"))
             CTX.mesh.format = FORMAT_CGNS;
           else if(!strcmp(argv[i], "med"))
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 13518c70ec..0b27d00dd6 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.516 2007-02-28 06:58:46 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.517 2007-03-11 20:18:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -607,6 +607,7 @@ static char *input_formats =
   "\tI-deas universal mesh (*.unv)"
   "\tMedit mesh (*.mesh)"
   "\tNastran bulk data file (*.bdf)"
+  "\tPlot3D structured mesh (*.p3d)"
   "\tSTL surface mesh (*.stl)"
   "\tVRML surface mesh (*.wrl)"
 #if defined(HAVE_LIBJPEG)
@@ -649,13 +650,14 @@ int _save_msh(char *name){ return msh_dialog(name); }
 int _save_pos(char *name){ return generic_mesh_dialog(name, "POS Options", FORMAT_POS); }
 int _save_options(char *name){ return options_dialog(name); }
 int _save_geo(char *name){ CreateOutputFile(name, FORMAT_GEO); return 1; }
+int _save_cgns(char *name){ CreateOutputFile(name, FORMAT_CGNS); return 1; }
 int _save_unv(char *name){ return generic_mesh_dialog(name, "UNV Options", FORMAT_UNV); }
+int _save_med(char *name){ return generic_mesh_dialog(name, "MED Options", FORMAT_MED); }
 int _save_mesh(char *name){ return generic_mesh_dialog(name, "MESH Options", FORMAT_MESH); }
 int _save_bdf(char *name){ return bdf_dialog(name); }
+int _save_p3d(char *name){ return generic_mesh_dialog(name, "P3D Options", FORMAT_P3D); }
 int _save_stl(char *name){ return stl_dialog(name); }
 int _save_vrml(char *name){ return generic_mesh_dialog(name, "VRML Options", FORMAT_VRML); }
-int _save_cgns(char *name){ CreateOutputFile(name, FORMAT_CGNS); return 1; }
-int _save_med(char *name){ return generic_mesh_dialog(name, "MED Options", FORMAT_MED); }
 int _save_eps(char *name){ return gl2ps_dialog(name, "EPS Options", FORMAT_EPS); }
 int _save_gif(char *name){ return gif_dialog(name); }
 int _save_jpeg(char *name){ return jpeg_dialog(name); }
@@ -674,13 +676,14 @@ int _save_auto(char *name)
   case FORMAT_POS  : return _save_pos(name);
   case FORMAT_OPT  : return _save_options(name);
   case FORMAT_GEO  : return _save_geo(name);
+  case FORMAT_CGNS : return _save_cgns(name);
   case FORMAT_UNV  : return _save_unv(name);
+  case FORMAT_MED  : return _save_med(name);
   case FORMAT_MESH : return _save_mesh(name);
   case FORMAT_BDF  : return _save_bdf(name);
+  case FORMAT_P3D  : return _save_p3d(name);
   case FORMAT_STL  : return _save_stl(name);
   case FORMAT_VRML : return _save_vrml(name);
-  case FORMAT_CGNS : return _save_cgns(name);
-  case FORMAT_MED  : return _save_med(name);
   case FORMAT_EPS  : return _save_eps(name);
   case FORMAT_GIF  : return _save_gif(name);
   case FORMAT_JPEG : return _save_jpeg(name);
@@ -712,17 +715,18 @@ void file_save_as_cb(CALLBACK_ARGS)
     {"Gmsh mesh statistics (*.pos)", _save_pos},
     {"Gmsh options (*.opt)", _save_options},
     {"Gmsh unrolled geometry (*.geo)", _save_geo},
-    {"I-deas universal mesh (*.unv)", _save_unv},
-    {"Medit mesh (*.mesh)", _save_mesh},
-    {"Nastran bulk data file (*.bdf)", _save_bdf},
-    {"STL surface mesh (*.stl)", _save_stl},
-    {"VRML surface mesh (*.wrl)", _save_vrml},
 #if defined(HAVE_LIBCGNS)
     {"CGNS (*.cgns)", _save_cgns},
 #endif
+    {"I-deas universal mesh (*.unv)", _save_unv},
 #if defined(HAVE_MED)
     {"MED (*.med)", _save_med},
 #endif
+    {"Medit mesh (*.mesh)", _save_mesh},
+    {"Nastran bulk data file (*.bdf)", _save_bdf},
+    {"Plot3D structured mesh (*.p3d)", _save_p3d},
+    {"STL surface mesh (*.stl)", _save_stl},
+    {"VRML surface mesh (*.wrl)", _save_vrml},
     {"Encapsulated PostScript (*.eps)", _save_eps},
     {"GIF (*.gif)", _save_gif},
 #if defined(HAVE_LIBJPEG)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 7fd70ce8d4..03c453cd5f 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -67,15 +67,15 @@ class GFace : public GEntity
   void delRegion(GRegion *r){ if(r1 == r) r1 = r2; r2=0; }
 
   // edge orientations.
-  virtual std::list<int> orientations() const{return l_dirs;}
+  virtual std::list<int> orientations() const { return l_dirs; }
   // Edges that bound this entity or that this entity bounds.
-  virtual std::list<GEdge*> edges() const{return l_edges;}
+  virtual std::list<GEdge*> edges() const { return l_edges; }
   // Edges that are embedded in this face.
-  virtual std::list<GEdge*> emb_edges() const{return embedded_edges;}
+  virtual std::list<GEdge*> emb_edges() const { return embedded_edges; }
   // Vertices that bound this entity or that this entity bounds.
   virtual std::list<GVertex*> vertices() const;
 
-  virtual int dim() const {return 2;}
+  virtual int dim() const { return 2; }
   virtual void setVisibility(char val, bool recursive=false);
 
   // compute the parameters UV from a point XYZ
@@ -162,8 +162,9 @@ class GFace : public GEntity
   // of start/end points
   std::vector<SPoint3> cross;
 
-  // a map for accessing the transfinite vertices using a pair of indices
-  std::map<std::pair<int, int>, MVertex*> transfinite_vertices;
+  // a array for accessing the transfinite vertices using a pair of
+  // indices
+  std::vector<std::vector<MVertex*> > transfinite_vertices;
 
   std::vector<MTriangle*> triangles;
   std::vector<MQuadrangle*> quadrangles;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 3cd1759f8a..8a5410c502 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -52,7 +52,7 @@ class GModel
  public:
   GModel() : modelName("Untitled"), normals(0) {}
   GModel(const std::string &name) : modelName(name), normals(0) {}
-  virtual ~GModel(){ deleteOCCInternals(); destroy(); }
+  ~GModel(){ deleteOCCInternals(); destroy(); }
   
   typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
   typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
@@ -61,10 +61,10 @@ class GModel
   typedef std::map<int, std::string>::iterator piter;
 
   // Deletes everything in a GModel 
-  virtual void destroy();
+  void destroy();
 
   // Returns the geometric tolerance for the entire model.
-  virtual double tolerance() const { return 1.e-14; }
+  double tolerance() const { return 1.e-14; }
 
   // Get the number of regions in this model.
   int numRegion() const { return regions.size(); }
@@ -83,10 +83,10 @@ class GModel
   viter lastVertex() { return vertices.end(); }
 
   // Find the region with the given tag.
-  virtual GRegion *regionByTag(int n) const;
-  virtual GFace *faceByTag(int n) const;
-  virtual GEdge *edgeByTag(int n) const;
-  virtual GVertex *vertexByTag(int n) const;
+  GRegion *regionByTag(int n) const;
+  GFace *faceByTag(int n) const;
+  GEdge *edgeByTag(int n) const;
+  GVertex *vertexByTag(int n) const;
 
   void add(GRegion *r) { regions.insert(r); }
   void add(GFace *f) { faces.insert(f); }
@@ -131,26 +131,26 @@ class GModel
   std::string getPhysicalName(int number);
 
   // The bounding box
-  virtual SBoundingBox3d bounds();
+  SBoundingBox3d bounds();
 
   // Returns the mesh status for the entire model
-  virtual int getMeshStatus(bool countDiscrete=true);
+  int getMeshStatus(bool countDiscrete=true);
 
   // Returns the total number of vertices in the mesh
-  virtual int numVertices();
+  int numVertices();
 
   // Returns the total number of elements in the mesh
-  virtual int numElements();
+  int numElements();
 
   // The list of partitions
-  virtual std::set<int> &getMeshPartitions() { return meshPartitions; }
-  virtual std::set<int> &recomputeMeshPartitions();
+  std::set<int> &getMeshPartitions() { return meshPartitions; }
+  std::set<int> &recomputeMeshPartitions();
 
   // Deletes all the partitions
-  virtual void deleteMeshPartitions();
+  void deleteMeshPartitions();
 
   // Performs various coherence tests on the mesh
-  virtual void checkMeshCoherence();
+  void checkMeshCoherence();
 
   // A container for smooth normals
   smooth_normals *normals;
@@ -159,12 +159,12 @@ class GModel
   // =========================================
 
   // Gmsh native CAD format
-  virtual int importTHEM();
-  virtual int readGEO(const std::string &name);
-  virtual int writeGEO(const std::string &name);
+  int importTHEM();
+  int readGEO(const std::string &name);
+  int writeGEO(const std::string &name);
 
   // Fourier model
-  virtual int readFourier(const std::string &name);
+  int readFourier(const std::string &name);
 
   // OCC model
   int readOCCBREP(const std::string &name);
@@ -209,6 +209,11 @@ class GModel
   int writeBDF(const std::string &name, int format=0,
 	       bool saveAll=false, double scalingFactor=1.0);
 
+  // Plot3D structured mesh format
+  int readP3D(const std::string &name);
+  int writeP3D(const std::string &name, 
+	       bool saveAll=false, double scalingFactor=1.0);
+
   // CFD General Notation System files
   int readCGNS(const std::string &name);
   int writeCGNS(const std::string &name, double scalingFactor=1.0);
diff --git a/Geo/GModelIO_Fourier.cpp b/Geo/GModelIO_Fourier.cpp
index 5286d21cb4..e70198ca8f 100644
--- a/Geo/GModelIO_Fourier.cpp
+++ b/Geo/GModelIO_Fourier.cpp
@@ -89,8 +89,8 @@ public:
 
 #if 1
     switch(gf->tag()){
-    case 0: M = 30; N = 30; break; // falcon front
-    case 1: M = 30; N = 120; break; // falcon back
+    case 0: M = 30; N = 22; break; // falcon front
+    case 1: M = 30; N = 70; break; // falcon back
     case 2: break; // wing front
     case 3: M = 25; N = 30; break; // wing top (N along the fuselage)
     case 4: M = 25; N = 30; break; // wing bottom
@@ -103,6 +103,9 @@ public:
       for(int j = 0; j < N; j++){
 	double u = i/(double)(M - 1);
 	double v = j/(double)(N - 1);
+	if(gf->tag()==0){
+	  v = sin(M_PI/2.*v - M_PI/2.) + 1.; 
+	}
 
 	//if(gf->tag() == 2){
 	  // hack for sttr report 2 (ship model, top right patch)
@@ -123,12 +126,13 @@ public:
     }
     for(int i = 0; i < M - 1; i++){
       for(int j = 0; j < N - 1; j++){
-#if 0
+#if 1
 	MQuadrangle *q = new MQuadrangle(gf->mesh_vertices[i * N + j],
 					 gf->mesh_vertices[(i + 1) * N + j],
 					 gf->mesh_vertices[(i + 1) * N + (j + 1)],
 					 gf->mesh_vertices[i * N + (j + 1)]);
-	//if(FM->GetOrientation(gf->tag()) < 0) q->revert();
+	//if(FM->GetOrientation(gf->tag()) < 0) 
+	q->revert();
 	gf->quadrangles.push_back(q);
 #else
 	MVertex *v1 = gf->mesh_vertices[i * N + j];
@@ -872,6 +876,20 @@ void cleanUpAndMergeAllFaces(GModel *m)
       if(v[0] != v[1] && v[0] != v[2] && v[1] != v[2])
 	newgf->triangles.push_back(new MTriangle(v));
     }
+    for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
+      std::vector<MVertex*> v(4);
+      for(int j = 0; j < 4; j++){
+	itp = pos.find(gf->quadrangles[i]->getVertex(j));
+	if(itp == pos.end())
+	  Msg(GERROR, "Could not find vertex");
+	else
+	  v[j] = *itp;
+      }
+      delete gf->quadrangles[i];
+      //if(v[0] != v[1] && v[0] != v[2] && v[0] != v[3] && 
+      // v[1] != v[2] && v[1] != v[3] && v[2] != v[3])
+	newgf->quadrangles.push_back(new MQuadrangle(v));
+    }
     m->remove(gf);
   }
   m->add(newgf);
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 3b4687cc2c..cd81571d38 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Geo.cpp,v 1.8 2007-01-25 15:50:57 geuzaine Exp $
+// $Id: GModelIO_Geo.cpp,v 1.9 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -146,12 +146,17 @@ class writeGVertexGEO {
   writeGVertexGEO(FILE *fp) { geo = fp ? fp : stdout; }
   void operator() (GVertex *gv)
   {
-    if(gv->getNativeType() != GEntity::GmshModel) return;
-    Vertex *v = (Vertex*)gv->getNativePtr();
-    if(!v) return;
-
-    fprintf(geo, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
-	    v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, v->lc);
+    if(gv->getNativeType() == GEntity::GmshModel){
+      Vertex *v = (Vertex*)gv->getNativePtr();
+      if(!v) return;
+      fprintf(geo, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+	      v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, v->lc);
+    }
+    else{
+      fprintf(geo, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+	      gv->tag(), gv->x(), gv->y(), gv->z(), 
+	      gv->prescribedMeshSizeAtVertex());
+    }
   }
 };
 
@@ -162,73 +167,98 @@ class writeGEdgeGEO {
   writeGEdgeGEO(FILE *fp) { geo = fp ? fp : stdout; }
   void operator () (GEdge *ge)
   {
-    if(ge->getNativeType() != GEntity::GmshModel) return;
-    Curve *c = (Curve *)ge->getNativePtr();
-    if(!c || c->Num < 0 || c->Typ == MSH_SEGM_DISCRETE) return;
-
-    switch (c->Typ) {
-    case MSH_SEGM_LINE:
-      fprintf(geo, "Line (%d) = ", c->Num);
-      break;
-    case MSH_SEGM_CIRC:
-    case MSH_SEGM_CIRC_INV:
-      fprintf(geo, "Circle (%d) = ", c->Num);
-      break;
-    case MSH_SEGM_ELLI:
-    case MSH_SEGM_ELLI_INV:
-      fprintf(geo, "Ellipse (%d) = ", c->Num);
-      break;
-    case MSH_SEGM_NURBS:
-      fprintf(geo, "Nurbs (%d) = {", c->Num);
+    if(ge->geomType() == GEntity::DiscreteCurve) return;
+    
+    if(ge->getNativeType() == GEntity::GmshModel){
+      Curve *c = (Curve *)ge->getNativePtr();
+      if(!c || c->Num < 0) return;
+      switch (c->Typ) {
+      case MSH_SEGM_LINE:
+	fprintf(geo, "Line (%d) = ", c->Num);
+	break;
+      case MSH_SEGM_CIRC:
+      case MSH_SEGM_CIRC_INV:
+	fprintf(geo, "Circle (%d) = ", c->Num);
+	break;
+      case MSH_SEGM_ELLI:
+      case MSH_SEGM_ELLI_INV:
+	fprintf(geo, "Ellipse (%d) = ", c->Num);
+	break;
+      case MSH_SEGM_NURBS:
+	fprintf(geo, "Nurbs (%d) = {", c->Num);
+	for(int i = 0; i < List_Nbr(c->Control_Points); i++) {
+	  Vertex *v;
+	  List_Read(c->Control_Points, i, &v);
+	  if(!i)
+	    fprintf(geo, "%d", v->Num);
+	  else
+	    fprintf(geo, ", %d", v->Num);
+	  if(i % 8 == 7 && i != List_Nbr(c->Control_Points) - 1)
+	    fprintf(geo, "\n");
+	}
+	fprintf(geo, "}\n");
+	fprintf(geo, "  Knots {");
+	for(int j = 0; j < List_Nbr(c->Control_Points) + c->degre + 1; j++) {
+	  if(!j)
+	    fprintf(geo, "%.16g", c->k[j]);
+	  else
+	    fprintf(geo, ", %.16g", c->k[j]);
+	  if(j % 5 == 4 && j != List_Nbr(c->Control_Points) + c->degre)
+	    fprintf(geo, "\n        ");
+	}
+	fprintf(geo, "}\n");
+	fprintf(geo, "  Order %d;\n", c->degre);
+	return;
+      case MSH_SEGM_SPLN:
+	fprintf(geo, "CatmullRom (%d) = ", c->Num);
+	break;
+      case MSH_SEGM_BSPLN:
+	fprintf(geo, "BSpline (%d) = ", c->Num);
+	break;
+      case MSH_SEGM_BEZIER:
+	fprintf(geo, "Bezier (%d) = ", c->Num);
+	break;
+      default:
+	Msg(GERROR, "Unknown curve type %d", c->Typ);
+	return;
+      }
       for(int i = 0; i < List_Nbr(c->Control_Points); i++) {
 	Vertex *v;
 	List_Read(c->Control_Points, i, &v);
-	if(!i)
-	  fprintf(geo, "%d", v->Num);
-	else
+	if(i)
 	  fprintf(geo, ", %d", v->Num);
-	if(i % 8 == 7 && i != List_Nbr(c->Control_Points) - 1)
-	  fprintf(geo, "\n");
-      }
-      fprintf(geo, "}\n");
-      fprintf(geo, "  Knots {");
-      for(int j = 0; j < List_Nbr(c->Control_Points) + c->degre + 1; j++) {
-	if(!j)
-	  fprintf(geo, "%.16g", c->k[j]);
 	else
-	  fprintf(geo, ", %.16g", c->k[j]);
-	if(j % 5 == 4 && j != List_Nbr(c->Control_Points) + c->degre)
-	  fprintf(geo, "\n        ");
+	  fprintf(geo, "{%d", v->Num);
+	if(i % 6 == 7)
+	  fprintf(geo, "\n");
       }
-      fprintf(geo, "}\n");
-      fprintf(geo, "  Order %d;\n", c->degre);
-      return;
-    case MSH_SEGM_SPLN:
-      fprintf(geo, "CatmullRom (%d) = ", c->Num);
-      break;
-    case MSH_SEGM_BSPLN:
-      fprintf(geo, "BSpline (%d) = ", c->Num);
-      break;
-    case MSH_SEGM_BEZIER:
-      fprintf(geo, "Bezier (%d) = ", c->Num);
-      break;
-    default:
-      Msg(GERROR, "Unknown curve type %d", c->Typ);
-      return;
+      fprintf(geo, "};\n");
     }
-
-    for(int i = 0; i < List_Nbr(c->Control_Points); i++) {
-      Vertex *v;
-      List_Read(c->Control_Points, i, &v);
-      if(i)
-	fprintf(geo, ", %d", v->Num);
-      else
-	fprintf(geo, "{%d", v->Num);
-      if(i % 6 == 7)
-	fprintf(geo, "\n");
+    else{
+      if(ge->getBeginVertex() && ge->getEndVertex()){
+	if(ge->geomType() == GEntity::Line){
+	  fprintf(geo, "Line (%d) = {%d, %d};\n", 
+		  ge->tag(), ge->getBeginVertex()->tag(), ge->getEndVertex()->tag());
+	}
+	else{
+	  // approximate all other curves by splines
+	  Range<double> bounds = ge->parBounds(0);
+	  double umin = bounds.low();
+	  double umax = bounds.high();
+	  fprintf(geo, "p%d = newp;\n", ge->tag());
+	  for(int i = 1; i < ge->minimumDrawSegments(); i++){
+	    double u = umin + (double)i / ge->minimumDrawSegments() * (umax - umin);
+	    GPoint p = ge->point(u);
+	    fprintf(geo, "Point(p%d + %d) = {%.16g, %.16g, %.16g, 1.e+22};\n", 
+		    ge->tag(), i, p.x(), p.y(), p.z());
+	  }
+	  fprintf(geo, "CatmullRom (%d) = {%d", ge->tag(), ge->getBeginVertex()->tag());
+	  for(int i = 1; i < ge->minimumDrawSegments(); i++)
+	    fprintf(geo, ", p%d + %d", ge->tag(), i);
+	  fprintf(geo, ", %d};\n", ge->getEndVertex()->tag());
+	}
+      }
     }
-    
-    fprintf(geo, "};\n");
   }
 };
 
@@ -239,32 +269,34 @@ class writeGFaceGEO {
   writeGFaceGEO(FILE *fp) { geo = fp ? fp : stdout; }
   void operator () (GFace *gf)
   {
-    if(gf->getNativeType() != GEntity::GmshModel) return;
-    Surface *s = (Surface *)gf->getNativePtr();
-    if(!s || s->Typ == MSH_SURF_DISCRETE) return;
-    
-    int NUMLOOP = s->Num + 1000000;
-    if(List_Nbr(s->Generatrices)){
+    if(gf->geomType() == GEntity::DiscreteSurface) return;
+
+    std::list<GEdge*> edges = gf->edges();
+    std::list<int> orientations = gf->orientations();
+    if(edges.size() && orientations.size() == edges.size()){
+      std::vector<int> num, ori;
+      for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++)
+	num.push_back((*it)->tag());
+      for(std::list<int>::iterator it = orientations.begin(); it != orientations.end(); it++)
+	ori.push_back((*it) > 0 ? 1 : -1);
+      int NUMLOOP = gf->tag() + 1000000;
       fprintf(geo, "Line Loop (%d) = ", NUMLOOP);
-      for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
-	Curve *c;
-	List_Read(s->Generatrices, i, &c);
+      for(unsigned int i = 0; i < num.size(); i++){
 	if(i)
-	  fprintf(geo, ", %d", c->Num);
+	  fprintf(geo, ", %d", num[i] * ori[i]);
 	else
-	  fprintf(geo, "{%d", c->Num);
+	  fprintf(geo, "{%d", num[i] * ori[i]);
       }
       fprintf(geo, "};\n");
-    }
-    
-    switch (s->Typ) {
-    case MSH_SURF_REGL:
-    case MSH_SURF_TRIC:
-      fprintf(geo, "Ruled Surface (%d) = {%d};\n", s->Num, NUMLOOP);
-      break;
-    case MSH_SURF_PLAN:
-      fprintf(geo, "Plane Surface (%d) = {%d};\n", s->Num, NUMLOOP);
-      break;
+      if(gf->geomType() == GEntity::Plane){
+	fprintf(geo, "Plane Surface (%d) = {%d};\n", gf->tag(), NUMLOOP);
+      }
+      else if(edges.size() == 3 || edges.size() == 4){
+	fprintf(geo, "Ruled Surface (%d) = {%d};\n", gf->tag(), NUMLOOP);
+      }
+      else{
+	Msg(GERROR, "Skipping surface %d in export", gf->tag());
+      }
     }
   }
 };
@@ -276,28 +308,20 @@ class writeGRegionGEO {
   writeGRegionGEO(FILE *fp) { geo = fp ? fp : stdout; }
   void operator () (GRegion *gr)
   {
-    if(gr->getNativeType() != GEntity::GmshModel) return;
-    Volume *vol = (Volume *)gr->getNativePtr();
-    if(!vol || vol->Typ == MSH_VOLUME_DISCRETE) return;
-    
-    int NUMLOOP = vol->Num + 1000000;
-    
-    fprintf(geo, "Surface Loop (%d) = ", NUMLOOP);
-    
-    for(int i = 0; i < List_Nbr(vol->Surfaces); i++) {
-      Surface *s;
-      List_Read(vol->Surfaces, i, &s);
-      if(i)
-	fprintf(geo, ", %d", s->Num);
-      else
-	fprintf(geo, "{%d", s->Num);
-    }
-    fprintf(geo, "};\n");
-    
-    switch (vol->Typ) {
-    case MSH_VOLUME:
-      fprintf(geo, "Volume (%d) = {%d};\n", vol->Num, NUMLOOP);
-      break;
+    if(gr->geomType() == GEntity::DiscreteVolume) return;
+
+    std::list<GFace*> faces = gr->faces();
+    if(faces.size()){
+      int NUMLOOP = gr->tag() + 1000000;
+      fprintf(geo, "Surface Loop (%d) = ", NUMLOOP);
+      for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++) {
+	if(it != faces.begin())
+	  fprintf(geo, ", %d", (*it)->tag());
+	else
+	  fprintf(geo, "{%d", (*it)->tag());
+      }
+      fprintf(geo, "};\n");
+      fprintf(geo, "Volume (%d) = {%d};\n", gr->tag(), NUMLOOP);
     }
   }
 };
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index e29ccc51bb..c1ea5a461a 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.9 2007-02-27 22:09:44 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.10 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1894,3 +1894,178 @@ int GModel::writeBDF(const std::string &name, int format, bool saveAll,
   fclose(fp);
   return 1;
 }
+
+int GModel::readP3D(const std::string &name)
+{
+  FILE *fp = fopen(name.c_str(), "r");
+  if(!fp){
+    Msg(GERROR, "Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+
+  int numBlocks = 0;
+  if(fscanf(fp, "%d", &numBlocks) != 1 || numBlocks <= 0) return 0;
+
+  std::vector<int> Ni(numBlocks), Nj(numBlocks), Nk(numBlocks);
+  for(int n = 0; n < numBlocks; n++)
+    if(fscanf(fp, "%d %d %d", &Ni[n], &Nj[n], &Nk[n]) != 3) return 0;
+
+  for(int n = 0; n < numBlocks; n++){
+    if(Nk[n] == 1){
+      GFace *gf = new gmshFace(this, numFace() + 1);
+      add(gf);
+      gf->transfinite_vertices.resize(Ni[n]);
+      for(int i = 0; i < Ni[n]; i++)
+	gf->transfinite_vertices[i].resize(Nj[n]);
+      for(int coord = 0; coord < 3; coord++){
+	for(int i = 0; i < Ni[n]; i++){
+	  for(int j = 0; j < Nj[n]; j++){
+	    double d;
+	    if(fscanf(fp, "%lf", &d) != 1) return 0;
+	    if(coord == 0){
+	      MVertex *v = new MVertex(d, 0., 0., gf);
+	      gf->transfinite_vertices[i][j] = v;
+	      gf->mesh_vertices.push_back(v);
+	    }
+	    else if(coord == 1){
+	      gf->transfinite_vertices[i][j]->y() = d;
+	    }
+	    else if(coord == 2){
+	      gf->transfinite_vertices[i][j]->z() = d;
+	    }
+	  }
+	}
+      }
+      for(unsigned int i = 0; i < gf->transfinite_vertices.size() - 1; i++)
+	for(unsigned int j = 0; j < gf->transfinite_vertices[0].size() - 1; j++)
+	  gf->quadrangles.push_back
+	    (new MQuadrangle(gf->transfinite_vertices[i    ][j    ],
+			     gf->transfinite_vertices[i + 1][j    ],
+			     gf->transfinite_vertices[i + 1][j + 1],
+			     gf->transfinite_vertices[i    ][j + 1]));
+    }
+    else{
+      GRegion *gr = new gmshRegion(this, numRegion() + 1);
+      add(gr);
+      gr->transfinite_vertices.resize(Ni[n]);
+      for(int i = 0; i < Ni[n]; i++){
+	gr->transfinite_vertices[i].resize(Nj[n]);
+	for(int j = 0; j < Nj[n]; j++){
+	  gr->transfinite_vertices[i][j].resize(Nk[n]);
+	}
+      }
+      for(int coord = 0; coord < 3; coord++){
+	for(int i = 0; i < Ni[n]; i++){
+	  for(int j = 0; j < Nj[n]; j++){
+	    for(int k = 0; k < Nk[n]; k++){
+	      double d;
+	      if(fscanf(fp, "%lf", &d) != 1){printf("aaa\n"); return 0;}
+	      if(coord == 0){
+		MVertex *v = new MVertex(d, 0., 0., gr);
+		gr->transfinite_vertices[i][j][k] = v;
+		gr->mesh_vertices.push_back(v);
+	      }
+	      else if(coord == 1){
+		gr->transfinite_vertices[i][j][k]->y() = d;
+	      }
+	      else if(coord == 2){
+		gr->transfinite_vertices[i][j][k]->z() = d;
+	      }
+	    }
+	  }
+	}
+      }
+      for(unsigned int i = 0; i < gr->transfinite_vertices.size() - 1; i++)
+	for(unsigned int j = 0; j < gr->transfinite_vertices[0].size() - 1; j++)
+	  for(unsigned int k = 0; k < gr->transfinite_vertices[0][0].size() - 1; k++)
+	    gr->hexahedra.push_back
+	      (new MHexahedron(gr->transfinite_vertices[i    ][j    ][k    ],
+			       gr->transfinite_vertices[i + 1][j    ][k    ],
+			       gr->transfinite_vertices[i + 1][j + 1][k    ],
+			       gr->transfinite_vertices[i    ][j + 1][k    ],
+			       gr->transfinite_vertices[i    ][j    ][k + 1],
+			       gr->transfinite_vertices[i + 1][j    ][k + 1],
+			       gr->transfinite_vertices[i + 1][j + 1][k + 1],
+			       gr->transfinite_vertices[i    ][j + 1][k + 1]));
+    }
+  }  
+  
+  fclose(fp);
+  return 1;
+}
+
+int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor)
+{
+  FILE *fp = fopen(name.c_str(), "w");
+  if(!fp){
+    Msg(GERROR, "Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+
+  if(noPhysicalGroups()) saveAll = true;
+
+  std::vector<GFace*> faces;
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+    if((*it)->transfinite_vertices.size() && 
+       (*it)->transfinite_vertices[0].size() &&
+       ((*it)->physicals.size() || saveAll)) faces.push_back(*it);
+
+  std::vector<GRegion*> regions;
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    if((*it)->transfinite_vertices.size() && 
+       (*it)->transfinite_vertices[0].size() && 
+       (*it)->transfinite_vertices[0][0].size() && 
+       ((*it)->physicals.size() || saveAll)) regions.push_back(*it);
+  
+  if(faces.empty() && regions.empty()){
+    Msg(WARNING, "No structured grids to save");
+    fclose(fp);
+    return 0;
+  }
+
+  fprintf(fp, "%d\n", faces.size() + regions.size());
+  
+  for(unsigned int i = 0; i < faces.size(); i++)
+    fprintf(fp, "%d %d 1\n", 
+	    faces[i]->transfinite_vertices.size(),
+	    faces[i]->transfinite_vertices[0].size());
+
+  for(unsigned int i = 0; i < regions.size(); i++)
+    fprintf(fp, "%d %d %d\n", 
+	    regions[i]->transfinite_vertices.size(),
+	    regions[i]->transfinite_vertices[0].size(),
+	    regions[i]->transfinite_vertices[0][0].size());
+  
+  for(unsigned int i = 0; i < faces.size(); i++){
+    GFace *gf = faces[i];
+    for(int coord = 0; coord < 3; coord++){
+      for(int j = 0; j < gf->transfinite_vertices.size(); j++){
+	for(int k = 0; k < gf->transfinite_vertices[j].size(); k++){
+	  MVertex *v = gf->transfinite_vertices[j][k];
+	  double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
+	  fprintf(fp, "%g ", d * scalingFactor);
+	}
+	fprintf(fp, "\n");
+      }
+    }
+  }
+  
+  for(unsigned int i = 0; i < regions.size(); i++){
+    GRegion *gr = regions[i];
+    for(int coord = 0; coord < 3; coord++){
+      for(int j = 0; j < gr->transfinite_vertices.size(); j++){
+	for(int k = 0; k < gr->transfinite_vertices[j].size(); k++){
+	  for(int l = 0; l < gr->transfinite_vertices[j][k].size(); l++){
+	    MVertex *v = gr->transfinite_vertices[j][k][l];
+	    double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
+	    fprintf(fp, "%g ", d * scalingFactor);
+	  }
+	  fprintf(fp, "\n");
+	}
+      }
+    }
+  }
+  
+  fclose(fp);
+  return 1;
+}
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index fb8e1e19a7..2d93c161dd 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -34,11 +34,11 @@ class GRegion : public GEntity {
   GRegion(GModel *model, int tag);
   virtual ~GRegion();
 
-  virtual int dim() const {return 3;}
+  virtual int dim() const { return 3; }
   virtual void setVisibility(char val, bool recursive=false);
-  virtual std::list<GFace*> faces() const{return l_faces;}
+  virtual std::list<GFace*> faces() const{ return l_faces; }
   virtual std::list<GEdge*> edges() const;
-  void set(std::list<GFace*> &f) {l_faces= f;}
+  void set(std::list<GFace*> &f) { l_faces= f; }
 
   // The bounding box
   virtual SBoundingBox3d bounds() const; 
@@ -60,6 +60,10 @@ class GRegion : public GEntity {
     std::vector<GVertex*> corners;
   } meshAttributes ;
 
+  // a array for accessing the transfinite vertices using a triplet of
+  // indices
+  std::vector<std::vector<std::vector<MVertex*> > > transfinite_vertices;
+
   std::vector<MTetrahedron*> tetrahedra;
   std::vector<MHexahedron*> hexahedra;
   std::vector<MPrism*> prisms;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 20cfadc15a..a29c9743cd 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.85 2007-03-05 15:00:28 geuzaine Exp $
+// $Id: Geo.cpp,v 1.86 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1127,6 +1127,12 @@ void DeleteShape(int Type, int Num)
   case MSH_VOLUME:
     DeleteVolume(Num);
     break;
+  case MSH_POINT_FROM_GMODEL:
+  case MSH_SEGM_FROM_GMODEL:
+  case MSH_SURF_FROM_GMODEL:
+  case MSH_VOLUME_FROM_GMODEL:
+    Msg(GERROR, "Deletion of external CAD entities not implemented yet");
+    break;
   default:
     Msg(GERROR, "Impossible to delete entity %d (of type %d)", Num, Type);
     break;
diff --git a/Geo/Geo.h b/Geo/Geo.h
index 75776b757d..87ed999f3a 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -29,41 +29,42 @@
 #include "SPoint2.h"
 #include "ExtrudeParams.h"
 
-#define MSH_POINT            1
-#define MSH_POINT_BND_LAYER  20
-#define MSH_POINT_DISCRETE   100
-
-#define MSH_SEGM_LINE        2
-#define MSH_SEGM_SPLN        3
-#define MSH_SEGM_CIRC        4
-#define MSH_SEGM_CIRC_INV    5
-#define MSH_SEGM_ELLI        6
-#define MSH_SEGM_ELLI_INV    7
-#define MSH_SEGM_LOOP        8
-#define MSH_SEGM_BSPLN       15
-#define MSH_SEGM_BND_LAYER   16
-#define MSH_SEGM_NURBS       17
-#define MSH_SEGM_BEZIER      18
-#define MSH_SEGM_PARAMETRIC  19
-#define MSH_SEGM_DISCRETE    101
-
-#define MSH_SURF_PLAN        9
-#define MSH_SURF_REGL        10
-#define MSH_SURF_TRIC        11
-#define MSH_SURF_BND_LAYER   12
-#define MSH_SURF_LOOP        13
-#define MSH_SURF_DISCRETE    102
-
-#define MSH_VOLUME           14
-#define MSH_VOLUME_DISCRETE  103
-
-#define MSH_PHYSICAL_POINT   300
-#define MSH_PHYSICAL_LINE    310
-#define MSH_PHYSICAL_SURFACE 320
-#define MSH_PHYSICAL_VOLUME  330
-
-#define MSH_POINT_ATTRACTOR  400
-#define MSH_LINE_ATTRACTOR   401
+#define MSH_POINT              100
+#define MSH_POINT_BND_LAYER    101
+#define MSH_POINT_DISCRETE     102
+#define MSH_POINT_FROM_GMODEL  103
+
+#define MSH_SEGM_LINE          200
+#define MSH_SEGM_SPLN          201
+#define MSH_SEGM_CIRC          202
+#define MSH_SEGM_CIRC_INV      203
+#define MSH_SEGM_ELLI          204
+#define MSH_SEGM_ELLI_INV      205
+#define MSH_SEGM_LOOP          206
+#define MSH_SEGM_BSPLN         207
+#define MSH_SEGM_NURBS         208
+#define MSH_SEGM_BEZIER        209
+#define MSH_SEGM_PARAMETRIC    210
+#define MSH_SEGM_BND_LAYER     211
+#define MSH_SEGM_DISCRETE      212
+#define MSH_SEGM_FROM_GMODEL   213
+
+#define MSH_SURF_PLAN          300
+#define MSH_SURF_REGL          301
+#define MSH_SURF_TRIC          302
+#define MSH_SURF_BND_LAYER     303
+#define MSH_SURF_LOOP          304
+#define MSH_SURF_DISCRETE      305
+#define MSH_SURF_FROM_GMODEL   306
+
+#define MSH_VOLUME             400
+#define MSH_VOLUME_DISCRETE    401
+#define MSH_VOLUME_FROM_GMODEL 402
+
+#define MSH_PHYSICAL_POINT     500
+#define MSH_PHYSICAL_LINE      501
+#define MSH_PHYSICAL_SURFACE   502
+#define MSH_PHYSICAL_VOLUME    503
 
 struct Coord{
   double X, Y, Z;
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index b619459c85..eff7449192 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.23 2007-03-05 11:07:14 remacle Exp $
+// $Id: GeoInterpolation.cpp,v 1.24 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -347,42 +347,43 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
     List_Read(c->Control_Points, i, &v[1]);
     List_Read(c->Control_Points, i + 1, &v[2]);
     if(!i) {
-			if(c->beg==c->end){
-				List_Read(c->Control_Points, N - 2, &v[0]);
-			}else{
-				v[0] = &temp1;
-				v[0]->Pos.X = 2. * v[1]->Pos.X - v[2]->Pos.X;
-				v[0]->Pos.Y = 2. * v[1]->Pos.Y - v[2]->Pos.Y;
-				v[0]->Pos.Z = 2. * v[1]->Pos.Z - v[2]->Pos.Z;
-				v[0]->pntOnGeometry = v[1]->pntOnGeometry * 2. - v[2]->pntOnGeometry;
-			}
+      if(c->beg == c->end){
+	List_Read(c->Control_Points, N - 2, &v[0]);
+      }
+      else{
+	v[0] = &temp1;
+	v[0]->Pos.X = 2. * v[1]->Pos.X - v[2]->Pos.X;
+	v[0]->Pos.Y = 2. * v[1]->Pos.Y - v[2]->Pos.Y;
+	v[0]->Pos.Z = 2. * v[1]->Pos.Z - v[2]->Pos.Z;
+	v[0]->pntOnGeometry = v[1]->pntOnGeometry * 2. - v[2]->pntOnGeometry;
+      }
     }
     else {
       List_Read(c->Control_Points, i - 1, &v[0]);
     }
     if(i == N - 2) {
-			if(c->beg==c->end){
-				List_Read(c->Control_Points, 1, &v[3]);
-			}else{
-				v[3] = &temp2;
-				v[3]->Pos.X = 2. * v[2]->Pos.X - v[1]->Pos.X;
-				v[3]->Pos.Y = 2. * v[2]->Pos.Y - v[1]->Pos.Y;
-				v[3]->Pos.Z = 2. * v[2]->Pos.Z - v[1]->Pos.Z;
-				v[3]->pntOnGeometry = v[2]->pntOnGeometry * 2. - v[1]->pntOnGeometry;
-			}
+      if(c->beg == c->end){
+	List_Read(c->Control_Points, 1, &v[3]);
+      }
+      else{
+	v[3] = &temp2;
+	v[3]->Pos.X = 2. * v[2]->Pos.X - v[1]->Pos.X;
+	v[3]->Pos.Y = 2. * v[2]->Pos.Y - v[1]->Pos.Y;
+	v[3]->Pos.Z = 2. * v[2]->Pos.Z - v[1]->Pos.Z;
+	v[3]->pntOnGeometry = v[2]->pntOnGeometry * 2. - v[1]->pntOnGeometry;
+      }
     }
     else {
       List_Read(c->Control_Points, i + 2, &v[3]);
     }
-    if (c->geometry)
-      {
-	SPoint2 pp =  InterpolateCubicSpline(v, t, c->mat, 0, t1, t2,c->geometry);
-	SPoint3 pt = c->geometry->point(pp);
-	V.Pos.X = pt.x();
-	V.Pos.Y = pt.y();
-	V.Pos.Z = pt.z();
-	return V;
-      }
+    if(c->geometry){
+      SPoint2 pp = InterpolateCubicSpline(v, t, c->mat, 0, t1, t2,c->geometry);
+      SPoint3 pt = c->geometry->point(pp);
+      V.Pos.X = pt.x();
+      V.Pos.Y = pt.y();
+      V.Pos.Z = pt.z();
+      return V;
+    }
     else
       return InterpolateCubicSpline(v, t, c->mat, 0, t1, t2);
 
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 44ef6dd0cd..99756033e6 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.10 2006-11-27 22:22:14 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.11 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -75,9 +75,7 @@ void MVertex::writeUNV(FILE *fp, double scalingFactor)
   char tmp[128];
   sprintf(tmp, "%25.16E%25.16E%25.16E\n", x() * scalingFactor, 
 	  y() * scalingFactor, z() * scalingFactor);
-  tmp[21] = 'D';
-  tmp[46] = 'D';
-  tmp[71] = 'D';
+  for(int i = 0; i < strlen(tmp); i++) if(tmp[i] == 'E') tmp[i] = 'D';
   fprintf(fp, tmp);
 }
 
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 7317807c0c..6eb2e0d8d4 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.197 2007-02-05 08:59:31 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.198 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -224,7 +224,7 @@ static void drawVertexLabel(GEntity *e, MVertex *v, int partition=-1)
   int physical = np ? e->physicals[np - 1] : 0;
   char str[256];
   if(CTX.mesh.label_type == 4)
-    sprintf(str, "(%g,%g,%g)", v->x(), v->y(), v->z());
+    sprintf(str, "(%.16g,%.16g,%.16g)", v->x(), v->y(), v->z());
   else if(CTX.mesh.label_type == 3){
     if(partition < 0)
       sprintf(str, "NA");
diff --git a/Makefile b/Makefile
index f0ed8386f3..e38bbb80ef 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.429 2007-02-28 06:58:46 geuzaine Exp $
+# $Id: Makefile,v 1.430 2007-03-11 20:18:57 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -23,7 +23,7 @@ include variables
 
 GMSH_MAJOR_VERSION = 2
 GMSH_MINOR_VERSION = 0
-GMSH_PATCH_VERSION = 4
+GMSH_PATCH_VERSION = 5
 GMSH_EXTRA_VERSION = "-cvs"
 
 GMSH_VERSION = ${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}.${GMSH_PATCH_VERSION}${GMSH_EXTRA_VERSION}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index c7b7b94d81..70e02e7090 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.65 2007-03-10 13:42:05 remacle Exp $
+// $Id: meshGFace.cpp,v 1.66 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1487,6 +1487,7 @@ void deMeshGFace::operator() (GFace *gf)
 
   for (unsigned int i=0;i<gf->mesh_vertices.size();i++) delete gf->mesh_vertices[i];
   gf->mesh_vertices.clear();
+  gf->transfinite_vertices.clear();
   for (unsigned int i=0;i<gf->triangles.size();i++) delete gf->triangles[i];
   gf->triangles.clear();
   for (unsigned int i=0;i<gf->quadrangles.size();i++) delete gf->quadrangles[i];
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index b7d32decd2..eca6f0556c 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceExtruded.cpp,v 1.17 2007-02-26 08:25:39 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.18 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -105,7 +105,8 @@ void extrudeMesh(GEdge *from, GFace *to,
 	  MVertex tmp(x[p], y[p], z[p], 0, -1);
 	  itp = pos.find(&tmp);
 	  if(itp == pos.end()) {
-	    Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
+	    Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
+		tmp.x(), tmp.y(), tmp.z(), to->tag());
 	    return;
 	  }
 	  verts.push_back(*itp);
@@ -143,7 +144,8 @@ void copyMesh(GFace *from, GFace *to,
 		  tmp.x(), tmp.y(), tmp.z());
       itp = pos.find(&tmp);
       if(itp == pos.end()) {
-	Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
+	Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
+	    tmp.x(), tmp.y(), tmp.z(), to->tag());
 	return;
       }
       verts.push_back(*itp);
@@ -159,7 +161,8 @@ void copyMesh(GFace *from, GFace *to,
 		  tmp.x(), tmp.y(), tmp.z());
       itp = pos.find(&tmp);
       if(itp == pos.end()) {
-	Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
+	Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d", 
+	    tmp.x(), tmp.y(), tmp.z(), to->tag());
 	return;
       }
       verts.push_back(*itp);
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index 13f7d8792e..116a3e7ea9 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.18 2007-02-21 08:17:16 geuzaine Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.19 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -195,36 +195,38 @@ int MeshTransfiniteSurface(GFace *gf)
            0            L
   */
 
-  std::map<std::pair<int,int>, MVertex*> &tab(gf->transfinite_vertices);
+  std::vector<std::vector<MVertex*> > &tab(gf->transfinite_vertices);
+  tab.resize(L + 1);
+  for(int i = 0; i <= L; i++) tab[i].resize(H + 1);
 
   if(corners.size () == 4){
-    tab[std::make_pair(0,0)] = m_vertices[0];
-    tab[std::make_pair(L,0)] = m_vertices[L];
-    tab[std::make_pair(L,H)] = m_vertices[L+H];
-    tab[std::make_pair(0,H)] = m_vertices[2*L+H];
+    tab[0][0] = m_vertices[0];
+    tab[L][0] = m_vertices[L];
+    tab[L][H] = m_vertices[L+H];
+    tab[0][H] = m_vertices[2*L+H];
     for (int i = 1; i < L; i++){
-      tab[std::make_pair(i,0)] = m_vertices[i];
-      tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
+      tab[i][0] = m_vertices[i];
+      tab[i][H] = m_vertices[2*L+H-i];
     }
     for(int i = 1; i < H; i++){
-      tab[std::make_pair(L,i)] = m_vertices[L+i];
-      tab[std::make_pair(0,i)] = m_vertices[2*L+2*H-i];
+      tab[L][i] = m_vertices[L+i];
+      tab[0][i] = m_vertices[2*L+2*H-i];
     }
   }
   else{
-    tab[std::make_pair(0,0)] = m_vertices[0];
-    tab[std::make_pair(L,0)] = m_vertices[L];
-    tab[std::make_pair(L,H)] = m_vertices[L+H];
+    tab[0][0] = m_vertices[0];
+    tab[L][0] = m_vertices[L];
+    tab[L][H] = m_vertices[L+H];
     // degenerated, only necessary for transfinite volume algo
-    tab[std::make_pair(0,H)] = m_vertices[0]; 
+    tab[0][H] = m_vertices[0]; 
     for (int i = 1; i < L; i++){
-      tab[std::make_pair(i,0)] = m_vertices[i];
-      tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
+      tab[i][0] = m_vertices[i];
+      tab[i][H] = m_vertices[2*L+H-i];
     }
     for(int i = 1; i < H;i++){
-      tab[std::make_pair(L,i)] = m_vertices[L+i];
+      tab[L][i] = m_vertices[L+i];
       // degenerated, only necessary for transfinite volume algo
-      tab[std::make_pair(0,i)] = m_vertices[0];
+      tab[0][i] = m_vertices[0];
     }
   }
 
@@ -252,7 +254,7 @@ int MeshTransfiniteSurface(GFace *gf)
 	GPoint gp = gf->point(SPoint2(Up, Vp));
 	MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp);
 	gf->mesh_vertices.push_back(newv);
-	tab[std::make_pair(i,j)] = newv;
+	tab[i][j] = newv;
       }
     }
   }
@@ -289,7 +291,7 @@ int MeshTransfiniteSurface(GFace *gf)
  	GPoint gp = gf->point(SPoint2(Up, Vp));
 	MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp);
 	gf->mesh_vertices.push_back(newv);
-	tab[std::make_pair(i,j)] = newv;
+	tab[i][j] = newv;
       }
     }
   }  
@@ -299,16 +301,15 @@ int MeshTransfiniteSurface(GFace *gf)
     for (int IT = 0; IT< CTX.mesh.nb_smoothing; IT++){
       for(int i = 1; i < L; i++){
 	for(int j = 1; j < H; j++){
-	  MVertex *v11 = tab[std::make_pair(i-1,j-1)];
-	  MVertex *v12 = tab[std::make_pair(i-1,j)];
-	  MVertex *v13 = tab[std::make_pair(i-1,j+1)];	      
-	  MVertex *v21 = tab[std::make_pair(i,j-1)];
-	  MVertex *v22 = tab[std::make_pair(i,j)];
-	  MVertex *v23 = tab[std::make_pair(i,j+1)];
-	  MVertex *v31 = tab[std::make_pair(i+1,j-1)];
-	  MVertex *v32 = tab[std::make_pair(i+1,j)];
-	  MVertex *v33 = tab[std::make_pair(i+1,j+1)];
-	  
+	  MVertex *v11 = tab[i - 1][j - 1];
+	  MVertex *v12 = tab[i - 1][j    ];
+	  MVertex *v13 = tab[i - 1][j + 1];	      
+	  MVertex *v21 = tab[i    ][j - 1];
+	  MVertex *v22 = tab[i    ][j    ];
+	  MVertex *v23 = tab[i    ][j + 1];
+	  MVertex *v31 = tab[i + 1][j - 1];
+	  MVertex *v32 = tab[i + 1][j    ];
+	  MVertex *v33 = tab[i + 1][j + 1];
 	  double alpha = 0.25 * (DSQR(v23->x() - v21->x()) +
 				 DSQR(v23->y() - v21->y()) +
 				 DSQR(v23->z() - v21->z()));
@@ -336,7 +337,7 @@ int MeshTransfiniteSurface(GFace *gf)
     // recompute corresponding u,v coordinates (necessary e.g. for 2nd order algo)
     for(int i = 1; i < L; i++){
       for(int j = 1; j < H; j++){
-	MVertex *v = tab[std::make_pair(i,j)];
+	MVertex *v = tab[i][j];
 	SPoint2 param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
 	v->setParameter(0, param[0]);
 	v->setParameter(1, param[1]);
@@ -348,10 +349,10 @@ int MeshTransfiniteSurface(GFace *gf)
     // create elements
     for(int i = 0; i < L ; i++){
       for(int j = 0; j < H; j++){
-	MVertex *v1 = tab[std::make_pair(i,j)];
-	MVertex *v2 = tab[std::make_pair(i+1,j)];
-	MVertex *v3 = tab[std::make_pair(i+1,j+1)];
-	MVertex *v4 = tab[std::make_pair(i,j+1)];
+	MVertex *v1 = tab[i    ][j    ];
+	MVertex *v2 = tab[i + 1][j    ];
+	MVertex *v3 = tab[i + 1][j + 1];
+	MVertex *v4 = tab[i    ][j + 1];
 	if(gf->meshAttributes.recombine)
 	  gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
 	else if(gf->meshAttributes.transfiniteArrangement == 1 ||
@@ -370,17 +371,17 @@ int MeshTransfiniteSurface(GFace *gf)
   }
   else{      
     for(int j = 0; j < H; j++){
-      MVertex *v1 = tab[std::make_pair(0,0)];
-      MVertex *v2 = tab[std::make_pair(1,j)];
-      MVertex *v3 = tab[std::make_pair(1,j+1)];
+      MVertex *v1 = tab[0    ][0    ];
+      MVertex *v2 = tab[1    ][j    ];
+      MVertex *v3 = tab[1    ][j + 1];
       gf->triangles.push_back(new MTriangle(v1, v2, v3));
     }
     for(int i = 1; i < L ; i++){
       for(int j = 0; j < H; j++){
-	MVertex *v1 = tab[std::make_pair(i,j)];
-	MVertex *v2 = tab[std::make_pair(i+1,j)];
-	MVertex *v3 = tab[std::make_pair(i+1,j+1)];
-	MVertex *v4 = tab[std::make_pair(i,j+1)];
+	MVertex *v1 = tab[i    ][j    ];
+	MVertex *v2 = tab[i + 1][j    ];
+	MVertex *v3 = tab[i + 1][j + 1];
+	MVertex *v4 = tab[i    ][j + 1];
 	if(gf->meshAttributes.recombine)
 	  gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
 	else if(gf->meshAttributes.transfiniteArrangement == 1 ||
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 15fd98ea14..e6e611f819 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.28 2007-02-26 08:25:39 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.29 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -345,6 +345,7 @@ void deMeshGRegion::operator() (GRegion *gr)
   for(unsigned int i = 0; i < gr->mesh_vertices.size(); i++)
     delete gr->mesh_vertices[i];
   gr->mesh_vertices.clear();
+  gr->transfinite_vertices.clear();
   for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
     delete gr->tetrahedra[i];
   gr->tetrahedra.clear();
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index a14e81ff61..0acc05dd76 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionExtruded.cpp,v 1.11 2007-03-05 09:30:53 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.12 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -112,7 +112,8 @@ int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
     MVertex tmp(x[p], y[p], z[p], 0, -1);
     itp = pos.find(&tmp);
     if(itp == pos.end())
-      Msg(GERROR, "Could not find extruded vertex");
+      Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g)",
+	  tmp.x(), tmp.y(), tmp.z());
     else
       verts.push_back(*itp);
   }
diff --git a/Mesh/meshGRegionTransfinite.cpp b/Mesh/meshGRegionTransfinite.cpp
index fb4cc6cdc9..fb8df72906 100644
--- a/Mesh/meshGRegionTransfinite.cpp
+++ b/Mesh/meshGRegionTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionTransfinite.cpp,v 1.4 2007-01-16 11:31:42 geuzaine Exp $
+// $Id: meshGRegionTransfinite.cpp,v 1.5 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -59,58 +59,58 @@
     match.
 */
 
-#define CREATE_HEX new MHexahedron(tab[(i)   + N_i*(j)   + N_i*N_j*(k)],   \
-				   tab[(i+1) + N_i*(j)   + N_i*N_j*(k)],   \
-				   tab[(i+1) + N_i*(j+1) + N_i*N_j*(k)],   \
-				   tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				   tab[(i)   + N_i*(j)   + N_i*N_j*(k+1)], \
-				   tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				   tab[(i+1) + N_i*(j+1) + N_i*N_j*(k+1)], \
-				   tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)])
-
-#define CREATE_PRISM_1 new MPrism(tab[(i)   + N_i*(j)   + N_i*N_j*(k)],   \
-				  tab[(i+1) + N_i*(j)   + N_i*N_j*(k)],   \
-				  tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				  tab[(i)   + N_i*(j)   + N_i*N_j*(k+1)], \
-				  tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				  tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)])
-
-#define CREATE_PRISM_2 new MPrism(tab[(i+1) + N_i*(j+1) + N_i*N_j*(k)],   \
-				  tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				  tab[(i+1) + N_i*(j)   + N_i*N_j*(k)],   \
-				  tab[(i+1) + N_i*(j+1) + N_i*N_j*(k+1)], \
-				  tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)], \
-				  tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)])
-
-#define CREATE_SIM_1 new MTetrahedron(tab[(i)   + N_i*(j)   + N_i*N_j*(k)], \
-				      tab[(i+1) + N_i*(j)   + N_i*N_j*(k)], \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k)], \
-				      tab[(i)   + N_i*(j)   + N_i*N_j*(k+1)])
-
-#define CREATE_SIM_2 new MTetrahedron(tab[(i+1) + N_i*(j)   + N_i*N_j*(k)],   \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				      tab[(i)   + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)])
-
-#define CREATE_SIM_3 new MTetrahedron(tab[(i)   + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)])
-
-#define CREATE_SIM_4 new MTetrahedron(tab[(i+1) + N_i*(j)   + N_i*N_j*(k)],   \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				      tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j+1) + N_i*N_j*(k)])
-
-#define CREATE_SIM_5 new MTetrahedron(tab[(i)   + N_i*(j+1) + N_i*N_j*(k)],   \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j+1) + N_i*N_j*(k)])
-
-#define CREATE_SIM_6 new MTetrahedron(tab[(i+1) + N_i*(j)   + N_i*N_j*(k+1)], \
-				      tab[(i)   + N_i*(j+1) + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j+1) + N_i*N_j*(k+1)], \
-				      tab[(i+1) + N_i*(j+1) + N_i*N_j*(k)])
+#define CREATE_HEX new MHexahedron(tab[i    ][j    ][k    ], \
+				   tab[i + 1][j    ][k    ], \
+				   tab[i + 1][j + 1][k    ], \
+				   tab[i    ][j + 1][k    ], \
+				   tab[i    ][j    ][k + 1], \
+				   tab[i + 1][j    ][k + 1], \
+				   tab[i + 1][j + 1][k + 1], \
+				   tab[i    ][j + 1][k + 1])
+
+#define CREATE_PRISM_1 new MPrism(tab[i    ][j    ][k    ], \
+				  tab[i + 1][j    ][k    ], \
+				  tab[i    ][j + 1][k    ], \
+				  tab[i    ][j    ][k + 1], \
+				  tab[i + 1][j    ][k + 1], \
+				  tab[i    ][j + 1][k + 1])
+
+#define CREATE_PRISM_2 new MPrism(tab[i + 1][j + 1][k    ], \
+				  tab[i    ][j + 1][k    ], \
+				  tab[i + 1][j    ][k    ], \
+				  tab[i + 1][j + 1][k + 1], \
+				  tab[i    ][j + 1][k + 1], \
+				  tab[i + 1][j    ][k + 1])
+
+#define CREATE_SIM_1 new MTetrahedron(tab[i    ][j    ][k    ], \
+				      tab[i + 1][j    ][k    ], \
+				      tab[i    ][j + 1][k    ], \
+				      tab[i    ][j    ][k + 1])
+
+#define CREATE_SIM_2 new MTetrahedron(tab[i + 1][j    ][k    ], \
+				      tab[i    ][j + 1][k    ], \
+				      tab[i    ][j    ][k + 1], \
+				      tab[i + 1][j    ][k + 1])
+
+#define CREATE_SIM_3 new MTetrahedron(tab[i    ][j    ][k + 1], \
+				      tab[i + 1][j    ][k + 1], \
+				      tab[i    ][j + 1][k    ], \
+				      tab[i    ][j + 1][k + 1])
+
+#define CREATE_SIM_4 new MTetrahedron(tab[i + 1][j    ][k    ], \
+				      tab[i    ][j + 1][k    ], \
+				      tab[i + 1][j    ][k + 1], \
+				      tab[i + 1][j + 1][k    ])
+
+#define CREATE_SIM_5 new MTetrahedron(tab[i    ][j + 1][k    ], \
+				      tab[i    ][j + 1][k + 1], \
+				      tab[i + 1][j    ][k + 1], \
+				      tab[i + 1][j + 1][k    ])
+
+#define CREATE_SIM_6 new MTetrahedron(tab[i + 1][j    ][k + 1], \
+				      tab[i    ][j + 1][k + 1], \
+				      tab[i + 1][j + 1][k + 1], \
+				      tab[i + 1][j + 1][k    ])
 
 double transfiniteHex(double f1, double f2, double f3, double f4, 
 		      double f5, double f6,
@@ -174,15 +174,12 @@ public:
   GOrientedTransfiniteFace(GFace *gf, std::vector<GVertex*> &corners) 
     : _gf(gf), _L(0), _H(0), _permutation(-1), _index(-1)
   { 
-    std::map<std::pair<int, int>, MVertex*>::iterator it;
-    for(it = gf->transfinite_vertices.begin(); it != gf->transfinite_vertices.end(); ++it){
-      _L = std::max(_L, it->first.first);
-      _H = std::max(_H, it->first.second);
-    }
+    _L = gf->transfinite_vertices.size() - 1;
+    if(_L <= 0) return;
+    _H = gf->transfinite_vertices[0].size() - 1;
+    if(_H <= 0) return;
     Msg(DEBUG, "Face %d: L = %d  H = %d", gf->tag(), _L, _H);
 
-    if(!_L || !_H) return;
-
     // get the corners of the transfinite volume interpolation
     std::vector<MVertex*> s(8);
     if(corners.size() == 8){
@@ -205,16 +202,16 @@ public:
     // get the corners of the transfinite surface mesh
     std::vector<MVertex*> c(4);
     if(_gf->meshAttributes.corners.size() == 4){
-      c[0] = _gf->transfinite_vertices[std::make_pair(0, 0)];
-      c[1] = _gf->transfinite_vertices[std::make_pair(_L, 0)];
-      c[2] = _gf->transfinite_vertices[std::make_pair(_L, _H)];
-      c[3] = _gf->transfinite_vertices[std::make_pair(0, _H)];
+      c[0] = _gf->transfinite_vertices[0][0];
+      c[1] = _gf->transfinite_vertices[_L][0];
+      c[2] = _gf->transfinite_vertices[_L][_H];
+      c[3] = _gf->transfinite_vertices[0][_H];
     }
     else if(_gf->meshAttributes.corners.size() == 3){
-      c[0] = _gf->transfinite_vertices[std::make_pair(0, 0)];
-      c[1] = _gf->transfinite_vertices[std::make_pair(_L, 0)];
-      c[2] = _gf->transfinite_vertices[std::make_pair(_L, _H)];
-      c[3] = _gf->transfinite_vertices[std::make_pair(0, 0)];
+      c[0] = _gf->transfinite_vertices[0][0];
+      c[1] = _gf->transfinite_vertices[_L][0];
+      c[2] = _gf->transfinite_vertices[_L][_H];
+      c[3] = _gf->transfinite_vertices[0][0];
     }
     else
       return;
@@ -239,7 +236,7 @@ public:
     Msg(DEBUG, "Found face index %d  (permutation = %d)", _index, _permutation);
     for(int i = 0; i <= _L; i++)
       for(int j = 0; j <= _H; j++)
-	_list.push_back(_gf->transfinite_vertices[std::make_pair(i, j)]);
+	_list.push_back(_gf->transfinite_vertices[i][j]);
   }
 
   // returns the index of the face in the reference hexahedron
@@ -340,7 +337,14 @@ int MeshTransfiniteVolume(GRegion *gr)
   MVertex *s6 = orientedFaces[5].getVertex(N_i - 1, N_j - 1);
   MVertex *s7 = orientedFaces[5].getVertex(0, N_j - 1);
 
-  MVertex **tab = new MVertex*[N_i * N_j * N_k];
+  std::vector<std::vector<std::vector<MVertex*> > > &tab(gr->transfinite_vertices);
+  tab.resize(N_i);
+  for(int i = 0; i < N_i; i++){
+    tab[i].resize(N_j);
+    for(int j = 0; j < N_j; j++){
+      tab[i][j].resize(N_k);
+    }
+  }
 
   for(int i = 0; i < N_i; i++) {
     double u = lengths_i[i] / L_i;
@@ -378,33 +382,31 @@ int MeshTransfiniteVolume(GRegion *gr)
 	else
           f3 = c8;
 
-	int index = i + N_i * j + N_i * N_j * k;
-
         if(i && j && k && i != N_i - 1 && j != N_j - 1 && k != N_k - 1) {
           MVertex *newv = transfiniteHex(f0, f1, f2, f3, f4, f5,
 					 c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11,
 					 s0, s1, s2, s3, s4, s5, s6, s7,
 					 u, v, w);
 	  gr->mesh_vertices.push_back(newv);
-          tab[index] = newv;
+          tab[i][j][k] = newv;
         }
         else if(!i) {
-          tab[index] = f3;
+          tab[i][j][k] = f3;
         }
         else if(!j) {
-          tab[index] = f0;
+          tab[i][j][k] = f0;
         }
         else if(!k) {
-          tab[index] = f4;
+          tab[i][j][k] = f4;
         }
         else if(i == N_i - 1) {
-          tab[index] = f1;
+          tab[i][j][k] = f1;
         }
         else if(j == N_j - 1) {
-          tab[index] = f2;
+          tab[i][j][k] = f2;
         }
         else if(k == N_k - 1) {
-          tab[index] = f5;
+          tab[i][j][k] = f5;
         }
       }
     }
@@ -424,34 +426,34 @@ int MeshTransfiniteVolume(GRegion *gr)
           else if(!orientedFaces[0].recombined() && orientedFaces[1].recombined() && 
 		  !orientedFaces[2].recombined() && orientedFaces[3].recombined() && 
 		  orientedFaces[4].recombined() && orientedFaces[5].recombined()) {
-            gr->prisms.push_back(new MPrism(tab[(i) + N_i * (j) + N_i * N_j * (k)],
-					    tab[(i + 1) + N_i * (j) + N_i * N_j * (k)],
-					    tab[(i) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k)],
-					    tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k + 1)]));
-	    gr->prisms.push_back(new MPrism(tab[(i + 1) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i + 1) + N_i * (j) + N_i * N_j * (k)],
-					    tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k + 1)],
-					    tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k)]));
+            gr->prisms.push_back(new MPrism(tab[i    ][j    ][k    ],
+					    tab[i + 1][j    ][k    ],
+					    tab[i    ][j    ][k + 1],
+					    tab[i    ][j + 1][k    ],
+					    tab[i + 1][j + 1][k    ],
+					    tab[i    ][j + 1][k + 1]));
+	    gr->prisms.push_back(new MPrism(tab[i + 1][j    ][k + 1],
+					    tab[i    ][j    ][k + 1],
+					    tab[i + 1][j    ][k    ],
+					    tab[i + 1][j + 1][k + 1],
+					    tab[i    ][j + 1][k + 1],
+					    tab[i + 1][j + 1][k    ]));
           }
           else if(orientedFaces[0].recombined() && !orientedFaces[1].recombined() && 
 		  orientedFaces[2].recombined() && !orientedFaces[3].recombined() && 
 		  orientedFaces[4].recombined() && orientedFaces[5].recombined()) {
-            gr->prisms.push_back(new MPrism(tab[(i + 1) + N_i * (j) + N_i * N_j * (k)],
-					    tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k)],
-					    tab[(i + 1) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j) + N_i * N_j * (k)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k)],
-					    tab[(i) + N_i * (j) + N_i * N_j * (k + 1)]));
-            gr->prisms.push_back(new MPrism(tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k + 1)],
-					    tab[(i + 1) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i + 1) + N_i * (j + 1) + N_i * N_j * (k)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j) + N_i * N_j * (k + 1)],
-					    tab[(i) + N_i * (j + 1) + N_i * N_j * (k)]));
+            gr->prisms.push_back(new MPrism(tab[i + 1][j    ][k    ],
+					    tab[i + 1][j + 1][k    ],
+					    tab[i + 1][j    ][k + 1],
+					    tab[i    ][j    ][k    ],
+					    tab[i    ][j + 1][k    ],
+					    tab[i    ][j    ][k + 1]));
+            gr->prisms.push_back(new MPrism(tab[i + 1][j + 1][k + 1],
+					    tab[i + 1][j    ][k + 1],
+					    tab[i + 1][j + 1][k    ],
+					    tab[i    ][j + 1][k + 1],
+					    tab[i    ][j    ][k + 1],
+					    tab[i    ][j + 1][k    ]));
           }
           else if(orientedFaces[0].recombined() && orientedFaces[1].recombined() && 
 		  orientedFaces[2].recombined() && orientedFaces[3].recombined() && 
@@ -471,7 +473,6 @@ int MeshTransfiniteVolume(GRegion *gr)
           }
           else {
             Msg(GERROR, "Wrong surface recombination in transfinite volume %d", gr->tag());
-	    delete [] tab;
             return 0;
           }
         }
@@ -487,32 +488,31 @@ int MeshTransfiniteVolume(GRegion *gr)
            (orientedFaces[0].recombined() && orientedFaces[1].recombined() && 
 	    orientedFaces[2].recombined() && !orientedFaces[4].recombined() && 
 	    !orientedFaces[5].recombined())) {
-          gr->prisms.push_back(new MPrism(tab[N_i * (j) + N_i * N_j * (k)],
-					  tab[1 + N_i * (j) + N_i * N_j * (k)],
-					  tab[1 + N_i * (j + 1) + N_i * N_j * (k)],
-					  tab[N_i * (j) + N_i * N_j * (k + 1)],
-					  tab[1 + N_i * (j) + N_i * N_j * (k + 1)],
-					  tab[1 + N_i * (j + 1) + N_i * N_j * (k + 1)]));
+          gr->prisms.push_back(new MPrism(tab[0    ][j    ][k    ],
+					  tab[1    ][j    ][k    ],
+					  tab[1    ][j + 1][k    ],
+					  tab[0    ][j    ][k + 1],
+					  tab[1    ][j    ][k + 1],
+					  tab[1    ][j + 1][k + 1]));
         }
         else if(!orientedFaces[0].recombined() && !orientedFaces[1].recombined() && 
 		!orientedFaces[2].recombined() && !orientedFaces[4].recombined() && 
 		!orientedFaces[5].recombined()) {
-          gr->tetrahedra.push_back(new MTetrahedron(tab[+N_i * (j) + N_i * N_j * (k)],
-						    tab[1 + N_i * (j) + N_i * N_j * (k)],
-						    tab[1 + N_i * (j + 1) + N_i * N_j * (k)],
-						    tab[+N_i * (j) + N_i * N_j * (k + 1)]));
-          gr->tetrahedra.push_back(new MTetrahedron(tab[1 + N_i * (j) + N_i * N_j * (k)],
-						    tab[1 + N_i * (j + 1) + N_i * N_j * (k)],
-						    tab[+N_i * (j) + N_i * N_j * (k + 1)],
-						    tab[1 + N_i * (j) + N_i * N_j * (k + 1)]));
-          gr->tetrahedra.push_back(new MTetrahedron(tab[+N_i * (j) + N_i * N_j * (k + 1)],
-						    tab[1 + N_i * (j + 1) + N_i * N_j * (k + 1)],
-						    tab[1 + N_i * (j) + N_i * N_j * (k + 1)],
-						    tab[1 + N_i * (j + 1) + N_i * N_j * (k)]));
+          gr->tetrahedra.push_back(new MTetrahedron(tab[0    ][j    ][k    ],
+						    tab[1    ][j    ][k    ],
+						    tab[1    ][j + 1][k    ],
+						    tab[0    ][j    ][k + 1]));
+          gr->tetrahedra.push_back(new MTetrahedron(tab[1    ][j    ][k    ],
+						    tab[1    ][j + 1][k    ],
+						    tab[0    ][j    ][k + 1],
+						    tab[1    ][j    ][k + 1]));
+          gr->tetrahedra.push_back(new MTetrahedron(tab[0    ][j    ][k + 1],
+						    tab[1    ][j + 1][k + 1],
+						    tab[1    ][j    ][k + 1],
+						    tab[1    ][j + 1][k    ]));
         }
         else {
           Msg(GERROR, "Wrong surface recombination in transfinite volume %d", gr->tag());
-	  delete [] tab;
           return 0;
         }
       }
@@ -543,7 +543,6 @@ int MeshTransfiniteVolume(GRegion *gr)
           }
           else {
             Msg(GERROR, "Wrong surface recombination in transfinite volume %d", gr->tag());
-	    delete [] tab;
             return 0;
           }
         }
@@ -551,6 +550,5 @@ int MeshTransfiniteVolume(GRegion *gr)
     }
   }
 
-  delete [] tab;
   return 1;
 }
diff --git a/Parser/CreateFile.cpp b/Parser/CreateFile.cpp
index 5483fc6570..520b58c549 100644
--- a/Parser/CreateFile.cpp
+++ b/Parser/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.14 2007-01-29 14:52:58 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.15 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -64,6 +64,7 @@ int GuessFileFormatFromFileName(char *name)
   else if(!strcmp(ext, ".mesh")) return FORMAT_MESH;
   else if(!strcmp(ext, ".bdf"))  return FORMAT_BDF;
   else if(!strcmp(ext, ".nas"))  return FORMAT_BDF;
+  else if(!strcmp(ext, ".p3d"))  return FORMAT_P3D;
   else if(!strcmp(ext, ".wrl"))  return FORMAT_VRML;
   else if(!strcmp(ext, ".vrml")) return FORMAT_VRML;
   else if(!strcmp(ext, ".gif"))  return FORMAT_GIF;
@@ -95,6 +96,7 @@ void GetDefaultFileName(int format, char *name)
   case FORMAT_MED:  strcpy(ext, ".med"); break;
   case FORMAT_MESH: strcpy(ext, ".mesh"); break;
   case FORMAT_BDF:  strcpy(ext, ".bdf"); break;
+  case FORMAT_P3D:  strcpy(ext, ".p3d"); break;
   case FORMAT_VRML: strcpy(ext, ".wrl"); break;
   case FORMAT_GIF:  strcpy(ext, ".gif"); break;
   case FORMAT_JPEG: strcpy(ext, ".jpg"); break;
@@ -172,6 +174,10 @@ void CreateOutputFile(char *filename, int format)
 		     CTX.mesh.save_all, CTX.mesh.scaling_factor);
     break;
 
+  case FORMAT_P3D:
+    GMODEL->writeP3D(name, CTX.mesh.save_all, CTX.mesh.scaling_factor);
+    break;
+
   case FORMAT_CGNS:
     GMODEL->writeCGNS(name, CTX.mesh.scaling_factor);
     break;
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 03384fba52..ad7aab1549 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -125,7 +125,7 @@
 
 #line 1 "Gmsh.y"
 
-// $Id: Gmsh.tab.cpp,v 1.312 2007-03-09 14:57:06 remacle Exp $
+// $Id: Gmsh.tab.cpp,v 1.313 2007-03-11 20:18:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -523,33 +523,33 @@ static const short yyrline[] = { 0,
    574,   579,   583,   592,   594,   595,   596,   597,   600,   602,
    605,   640,   679,   733,   750,   768,   779,   798,   812,   829,
    855,   882,   896,   913,   927,   944,   964,   987,   997,  1011,
-  1016,  1024,  1051,  1067,  1095,  1125,  1147,  1165,  1183,  1201,
-  1227,  1245,  1271,  1291,  1315,  1339,  1365,  1382,  1401,  1419,
-  1458,  1463,  1468,  1474,  1496,  1518,  1534,  1554,  1571,  1588,
-  1608,  1614,  1619,  1624,  1631,  1633,  1634,  1637,  1642,  1646,
-  1662,  1678,  1694,  1714,  1727,  1733,  1739,  1750,  1760,  1770,
-  1786,  1804,  1818,  1825,  1831,  1840,  1853,  1904,  1919,  1930,
-  1950,  1960,  1982,  1986,  1991,  1996,  2006,  2023,  2039,  2065,
-  2092,  2124,  2131,  2136,  2142,  2146,  2154,  2163,  2171,  2179,
-  2183,  2191,  2195,  2203,  2207,  2215,  2219,  2226,  2230,  2242,
-  2249,  2256,  2263,  2270,  2277,  2284,  2291,  2298,  2305,  2309,
-  2316,  2320,  2327,  2331,  2338,  2342,  2349,  2353,  2360,  2364,
-  2371,  2375,  2382,  2386,  2393,  2397,  2407,  2411,  2416,  2426,
-  2448,  2472,  2476,  2499,  2518,  2536,  2554,  2583,  2618,  2623,
-  2650,  2664,  2682,  2689,  2695,  2698,  2706,  2716,  2718,  2719,
-  2720,  2721,  2722,  2723,  2724,  2725,  2732,  2733,  2734,  2735,
-  2736,  2737,  2738,  2739,  2740,  2741,  2742,  2743,  2744,  2745,
-  2746,  2747,  2748,  2749,  2750,  2751,  2752,  2753,  2754,  2755,
-  2756,  2757,  2758,  2759,  2760,  2761,  2762,  2763,  2765,  2766,
-  2767,  2768,  2769,  2770,  2771,  2772,  2773,  2774,  2775,  2776,
-  2777,  2778,  2779,  2780,  2781,  2782,  2783,  2784,  2785,  2790,
-  2795,  2796,  2797,  2798,  2799,  2800,  2804,  2820,  2835,  2855,
-  2869,  2882,  2905,  2923,  2941,  2959,  2977,  2984,  2989,  2993,
-  2997,  3001,  3007,  3012,  3016,  3020,  3026,  3032,  3039,  3045,
-  3049,  3054,  3058,  3066,  3076,  3085,  3093,  3099,  3110,  3130,
-  3140,  3150,  3160,  3177,  3203,  3209,  3213,  3217,  3229,  3234,
-  3246,  3253,  3274,  3279,  3293,  3299,  3305,  3310,  3318,  3326,
-  3340,  3354,  3358,  3377,  3399
+  1016,  1024,  1050,  1066,  1095,  1126,  1148,  1166,  1184,  1202,
+  1228,  1246,  1272,  1292,  1316,  1340,  1366,  1383,  1402,  1420,
+  1459,  1464,  1469,  1475,  1497,  1519,  1535,  1555,  1572,  1589,
+  1609,  1615,  1620,  1625,  1632,  1634,  1635,  1638,  1643,  1647,
+  1670,  1693,  1716,  1743,  1756,  1762,  1768,  1779,  1789,  1799,
+  1815,  1833,  1847,  1854,  1860,  1869,  1882,  1933,  1948,  1959,
+  1979,  1989,  2011,  2015,  2020,  2025,  2035,  2052,  2068,  2094,
+  2121,  2153,  2160,  2165,  2171,  2175,  2183,  2192,  2200,  2208,
+  2212,  2220,  2224,  2232,  2236,  2244,  2248,  2255,  2259,  2271,
+  2278,  2285,  2292,  2299,  2306,  2313,  2320,  2327,  2334,  2338,
+  2345,  2349,  2356,  2360,  2367,  2371,  2378,  2382,  2389,  2393,
+  2400,  2404,  2411,  2415,  2422,  2426,  2436,  2440,  2445,  2455,
+  2477,  2501,  2505,  2528,  2547,  2565,  2583,  2612,  2647,  2652,
+  2679,  2693,  2711,  2718,  2724,  2727,  2735,  2745,  2747,  2748,
+  2749,  2750,  2751,  2752,  2753,  2754,  2761,  2762,  2763,  2764,
+  2765,  2766,  2767,  2768,  2769,  2770,  2771,  2772,  2773,  2774,
+  2775,  2776,  2777,  2778,  2779,  2780,  2781,  2782,  2783,  2784,
+  2785,  2786,  2787,  2788,  2789,  2790,  2791,  2792,  2794,  2795,
+  2796,  2797,  2798,  2799,  2800,  2801,  2802,  2803,  2804,  2805,
+  2806,  2807,  2808,  2809,  2810,  2811,  2812,  2813,  2814,  2819,
+  2824,  2825,  2826,  2827,  2828,  2829,  2833,  2849,  2864,  2884,
+  2898,  2911,  2934,  2952,  2970,  2988,  3006,  3013,  3018,  3022,
+  3026,  3030,  3036,  3041,  3045,  3049,  3055,  3061,  3068,  3074,
+  3078,  3083,  3087,  3095,  3105,  3114,  3122,  3128,  3139,  3159,
+  3169,  3179,  3189,  3206,  3232,  3238,  3242,  3246,  3258,  3263,
+  3275,  3282,  3303,  3308,  3322,  3328,  3334,  3339,  3347,  3355,
+  3369,  3383,  3387,  3406,  3428
 };
 #endif
 
@@ -3700,20 +3700,19 @@ case 82:
 	double z = CTX.geom.scaling_factor * yyvsp[-1].v[2];
 	double lc = CTX.geom.scaling_factor * yyvsp[-1].v[3];
 	Vertex *v;
-	if (!myGmshSurface)
+	if(!myGmshSurface)
 	  v = Create_Vertex(num, x, y, z, lc, 1.0);
 	else
 	  v = Create_Vertex(num, x, y, myGmshSurface, lc);
-
 	Tree_Add(THEM->Points, &v);
-	AddToTemporaryBoundingBox(v->Pos.X,v->Pos.Y,v->Pos.Z);
+	AddToTemporaryBoundingBox(v->Pos.X, v->Pos.Y, v->Pos.Z);
       }
       yyval.s.Type = MSH_POINT;
       yyval.s.Num = num;
     ;
     break;}
 case 83:
-#line 1052 "Gmsh.y"
+#line 1051 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].i;
       if(FindPhysicalGroup(num, MSH_PHYSICAL_POINT)){
@@ -3731,7 +3730,7 @@ case 83:
     ;
     break;}
 case 84:
-#line 1068 "Gmsh.y"
+#line 1067 "Gmsh.y"
 {
       double pars[] = { CTX.lc/10, CTX.lc/100., CTX.lc/20, 10, 3 };
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
@@ -3748,15 +3747,16 @@ case 84:
 	List_Read(yyvsp[-3].l, i, &d);
 	Vertex *v = FindPoint((int)d); 
 	if(v)
-	  att->addPoint(v->Pos.X,v->Pos.Y,v->Pos.Z);
+	  att->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
 	else{
 	  GVertex *gv = GMODEL->vertexByTag((int)d);
 	  if(gv) 
-	    att->addPoint(gv->x(),gv->y(),gv->z());
+	    att->addPoint(gv->x(), gv->y(), gv->z());
 	}
       }
       att->buildFastSearchStructures();
-      yyval.s.Type = MSH_POINT_ATTRACTOR;
+      // dummy values
+      yyval.s.Type = 0;
       yyval.s.Num = 0;
     ;
     break;}
@@ -3778,22 +3778,23 @@ case 85:
 	List_Read(yyvsp[-3].l, i, &d);
 	Curve *c = FindCurve((int)d); 
 	if(c){
-	  buildListOfPoints( att , c , (int) pars[3] );
+	  buildListOfPoints(att, c, (int)pars[3]);
 	}
 	else{
 	  GEdge *ge = GMODEL->edgeByTag((int)d);
 	  if(ge){
-	    buildListOfPoints( att , ge , (int) pars[3] );
+	    buildListOfPoints(att, ge, (int)pars[3]);
 	  }
 	}
       }
       att->buildFastSearchStructures();
-      yyval.s.Type = MSH_LINE_ATTRACTOR;
+      // dummy values
+      yyval.s.Type = 0;
       yyval.s.Num = 0;
     ;
     break;}
 case 86:
-#line 1126 "Gmsh.y"
+#line 1127 "Gmsh.y"
 {      
       for(int i = 0; i < List_Nbr(yyvsp[-3].l); i++){
 	double d;
@@ -3814,7 +3815,7 @@ case 86:
     ;
     break;}
 case 87:
-#line 1148 "Gmsh.y"
+#line 1149 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindCurve(num)){
@@ -3834,7 +3835,7 @@ case 87:
     ;
     break;}
 case 88:
-#line 1166 "Gmsh.y"
+#line 1167 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindCurve(num)){
@@ -3854,7 +3855,7 @@ case 88:
     ;
     break;}
 case 89:
-#line 1184 "Gmsh.y"
+#line 1185 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindCurve(num)){
@@ -3874,7 +3875,7 @@ case 89:
     ;
     break;}
 case 90:
-#line 1202 "Gmsh.y"
+#line 1203 "Gmsh.y"
 {
       int num = (int)yyvsp[-6].d;
       if(FindCurve(num)){
@@ -3902,7 +3903,7 @@ case 90:
     ;
     break;}
 case 91:
-#line 1228 "Gmsh.y"
+#line 1229 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindCurve(num)){
@@ -3922,7 +3923,7 @@ case 91:
     ;
     break;}
 case 92:
-#line 1246 "Gmsh.y"
+#line 1247 "Gmsh.y"
 {
       int num = (int)yyvsp[-6].d;
       if(FindCurve(num)){
@@ -3950,7 +3951,7 @@ case 92:
     ;
     break;}
 case 93:
-#line 1273 "Gmsh.y"
+#line 1274 "Gmsh.y"
 {
       int num = (int)yyvsp[-14].d;
       if(FindCurve(num)){
@@ -3971,7 +3972,7 @@ case 93:
     ;
     break;}
 case 94:
-#line 1292 "Gmsh.y"
+#line 1293 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(List_Nbr(yyvsp[-1].l) < 4){
@@ -3997,7 +3998,7 @@ case 94:
     ;
     break;}
 case 95:
-#line 1316 "Gmsh.y"
+#line 1317 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(List_Nbr(yyvsp[-1].l) < 4){
@@ -4023,7 +4024,7 @@ case 95:
     ;
     break;}
 case 96:
-#line 1340 "Gmsh.y"
+#line 1341 "Gmsh.y"
 {
       int num = (int)yyvsp[-8].d;
       if(List_Nbr(yyvsp[-5].l) + (int)yyvsp[-1].d + 1 != List_Nbr(yyvsp[-3].l)){
@@ -4051,7 +4052,7 @@ case 96:
     ;
     break;}
 case 97:
-#line 1366 "Gmsh.y"
+#line 1367 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindEdgeLoop(num)){
@@ -4070,7 +4071,7 @@ case 97:
     ;
     break;}
 case 98:
-#line 1383 "Gmsh.y"
+#line 1384 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].i;
       if(FindPhysicalGroup(num, MSH_PHYSICAL_LINE)){
@@ -4088,7 +4089,7 @@ case 98:
     ;
     break;}
 case 99:
-#line 1402 "Gmsh.y"
+#line 1403 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindSurface(num)){
@@ -4108,7 +4109,7 @@ case 99:
     ;
     break;}
 case 100:
-#line 1420 "Gmsh.y"
+#line 1421 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d, type = 0;
       if(FindSurface(num)){
@@ -4148,26 +4149,26 @@ case 100:
     ;
     break;}
 case 101:
-#line 1459 "Gmsh.y"
+#line 1460 "Gmsh.y"
 {
     myGmshSurface = 0;
   ;
     break;}
 case 102:
-#line 1464 "Gmsh.y"
+#line 1465 "Gmsh.y"
 {
     myGmshSurface = gmshSurface :: surfaceByTag ( (int) yyvsp[-1].d);
   ;
     break;}
 case 103:
-#line 1469 "Gmsh.y"
+#line 1470 "Gmsh.y"
 {
     int num = (int)yyvsp[-6].d, type = 0;
     myGmshSurface = gmshParametricSurface::NewParametricSurface ((int)yyvsp[-6].d,yyvsp[-3].c,yyvsp[-2].c,yyvsp[-1].c);
   ;
     break;}
 case 104:
-#line 1475 "Gmsh.y"
+#line 1476 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d, type = 0;
       if (List_Nbr(yyvsp[-1].l) != 2){
@@ -4191,7 +4192,7 @@ case 104:
     ;
     break;}
 case 105:
-#line 1497 "Gmsh.y"
+#line 1498 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d, type = 0;
       if (List_Nbr(yyvsp[-1].l) != 2){
@@ -4215,7 +4216,7 @@ case 105:
     ;
     break;}
 case 106:
-#line 1519 "Gmsh.y"
+#line 1520 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindSurfaceLoop(num)){
@@ -4233,7 +4234,7 @@ case 106:
     ;
     break;}
 case 107:
-#line 1535 "Gmsh.y"
+#line 1536 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].i;
       if(FindPhysicalGroup(num, MSH_PHYSICAL_SURFACE)){
@@ -4251,7 +4252,7 @@ case 107:
     ;
     break;}
 case 108:
-#line 1555 "Gmsh.y"
+#line 1556 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindVolume(num)){
@@ -4270,7 +4271,7 @@ case 108:
     ;
     break;}
 case 109:
-#line 1572 "Gmsh.y"
+#line 1573 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].d;
       if(FindVolume(num)){
@@ -4289,7 +4290,7 @@ case 109:
     ;
     break;}
 case 110:
-#line 1589 "Gmsh.y"
+#line 1590 "Gmsh.y"
 {
       int num = (int)yyvsp[-4].i;
       if(FindPhysicalGroup(num, MSH_PHYSICAL_VOLUME)){
@@ -4307,59 +4308,59 @@ case 110:
     ;
     break;}
 case 111:
-#line 1610 "Gmsh.y"
+#line 1611 "Gmsh.y"
 {
       TranslateShapes(yyvsp[-3].v[0], yyvsp[-3].v[1], yyvsp[-3].v[2], yyvsp[-1].l);
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 112:
-#line 1615 "Gmsh.y"
+#line 1616 "Gmsh.y"
 {
       RotateShapes(yyvsp[-8].v[0], yyvsp[-8].v[1], yyvsp[-8].v[2], yyvsp[-6].v[0], yyvsp[-6].v[1], yyvsp[-6].v[2], yyvsp[-4].d, yyvsp[-1].l);
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 113:
-#line 1620 "Gmsh.y"
+#line 1621 "Gmsh.y"
 {
       SymmetryShapes(yyvsp[-3].v[0], yyvsp[-3].v[1], yyvsp[-3].v[2], yyvsp[-3].v[3], yyvsp[-1].l);
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 114:
-#line 1625 "Gmsh.y"
+#line 1626 "Gmsh.y"
 {
       DilatShapes(yyvsp[-6].v[0], yyvsp[-6].v[1], yyvsp[-6].v[2], yyvsp[-4].d, yyvsp[-1].l);
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 115:
-#line 1632 "Gmsh.y"
+#line 1633 "Gmsh.y"
 { yyval.l = yyvsp[0].l; ;
     break;}
 case 116:
-#line 1633 "Gmsh.y"
+#line 1634 "Gmsh.y"
 { yyval.l = yyvsp[0].l; ;
     break;}
 case 117:
-#line 1634 "Gmsh.y"
+#line 1635 "Gmsh.y"
 { yyval.l = yyvsp[0].l; ;
     break;}
 case 118:
-#line 1639 "Gmsh.y"
+#line 1640 "Gmsh.y"
 {
       yyval.l = List_Create(3, 3, sizeof(Shape));
     ;
     break;}
 case 119:
-#line 1643 "Gmsh.y"
+#line 1644 "Gmsh.y"
 {
       List_Add(yyval.l, &yyvsp[0].s);
     ;
     break;}
 case 120:
-#line 1647 "Gmsh.y"
+#line 1648 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-2].l); i++){
 	double d;
@@ -4367,17 +4368,24 @@ case 120:
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Vertex *v = FindPoint(TheShape.Num);
-	if(!v)
-	  yymsg(WARNING, "Unknown point %d", TheShape.Num);
-	else{
+	if(v){
 	  TheShape.Type = MSH_POINT;
 	  List_Add(yyval.l, &TheShape);
 	}
+	else{
+	  GVertex *gv = GMODEL->vertexByTag(TheShape.Num);
+	  if(gv){
+	    TheShape.Type = MSH_POINT_FROM_GMODEL;
+	    List_Add(yyval.l, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown point %d", TheShape.Num);
+	}
       }
     ;
     break;}
 case 121:
-#line 1663 "Gmsh.y"
+#line 1671 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-2].l); i++){
 	double d;
@@ -4385,17 +4393,24 @@ case 121:
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Curve *c = FindCurve(TheShape.Num);
-	if(!c)
-	  yymsg(WARNING, "Unknown curve %d", TheShape.Num);
-	else{
+	if(c){
 	  TheShape.Type = c->Typ;
 	  List_Add(yyval.l, &TheShape);
 	}
+	else{
+	  GEdge *ge = GMODEL->edgeByTag(TheShape.Num);
+	  if(ge){
+	    TheShape.Type = MSH_SEGM_FROM_GMODEL;
+	    List_Add(yyval.l, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown curve %d", TheShape.Num);
+	}
       }
     ;
     break;}
 case 122:
-#line 1679 "Gmsh.y"
+#line 1694 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-2].l); i++){
 	double d;
@@ -4403,17 +4418,24 @@ case 122:
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Surface *s = FindSurface(TheShape.Num);
-	if(!s)
-	  yymsg(WARNING, "Unknown surface %d", TheShape.Num);
-	else{
+	if(s){
 	  TheShape.Type = s->Typ;
 	  List_Add(yyval.l, &TheShape);
 	}
+	else{
+	  GFace *gf = GMODEL->faceByTag(TheShape.Num);
+	  if(gf){
+	    TheShape.Type = MSH_SURF_FROM_GMODEL;
+	    List_Add(yyval.l, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown surface %d", TheShape.Num);
+	}
       }
     ;
     break;}
 case 123:
-#line 1695 "Gmsh.y"
+#line 1717 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-2].l); i++){
 	double d;
@@ -4421,17 +4443,24 @@ case 123:
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Volume *v = FindVolume(TheShape.Num);
-	if(!v)
-	  yymsg(WARNING, "Unknown volume %d", TheShape.Num);
-	else{
+	if(v){
 	  TheShape.Type = v->Typ;
 	  List_Add(yyval.l, &TheShape);
 	}
+	else{
+	  GRegion *gr = GMODEL->regionByTag(TheShape.Num);
+	  if(gr){
+	    TheShape.Type = MSH_VOLUME_FROM_GMODEL;
+	    List_Add(yyval.l, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown volume %d", TheShape.Num);
+	}
       }
     ;
     break;}
 case 124:
-#line 1716 "Gmsh.y"
+#line 1745 "Gmsh.y"
 {
       yyval.l = List_Create(3, 3, sizeof(Shape));
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
@@ -4444,7 +4473,7 @@ case 124:
     ;
     break;}
 case 125:
-#line 1728 "Gmsh.y"
+#line 1757 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-4].c, "View")) AliasView((int)yyvsp[-2].d, 0);
       Free(yyvsp[-4].c);
@@ -4452,7 +4481,7 @@ case 125:
     ;
     break;}
 case 126:
-#line 1734 "Gmsh.y"
+#line 1763 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-4].c, "View")) AliasView((int)yyvsp[-2].d, 0);
       Free(yyvsp[-4].c);
@@ -4460,7 +4489,7 @@ case 126:
     ;
     break;}
 case 127:
-#line 1740 "Gmsh.y"
+#line 1769 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-4].c, "View")) AliasView((int)yyvsp[-2].d, 1);
       Free(yyvsp[-4].c);
@@ -4468,7 +4497,7 @@ case 127:
     ;
     break;}
 case 128:
-#line 1752 "Gmsh.y"
+#line 1781 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
 	Shape TheShape;
@@ -4479,7 +4508,7 @@ case 128:
     ;
     break;}
 case 129:
-#line 1761 "Gmsh.y"
+#line 1790 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-4].c, "View")){
 	RemoveViewByIndex((int)yyvsp[-2].d);
@@ -4491,7 +4520,7 @@ case 129:
     ;
     break;}
 case 130:
-#line 1771 "Gmsh.y"
+#line 1800 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-1].c, "Meshes") || !strcmp(yyvsp[-1].c, "All")){
 	GMODEL->destroy();
@@ -4509,7 +4538,7 @@ case 130:
     ;
     break;}
 case 131:
-#line 1787 "Gmsh.y"
+#line 1816 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-2].c, "Empty") && !strcmp(yyvsp[-1].c, "Views")){
 	for(int i = List_Nbr(CTX.post.list) - 1; i >= 0; i--){
@@ -4525,7 +4554,7 @@ case 131:
     ;
     break;}
 case 132:
-#line 1806 "Gmsh.y"
+#line 1835 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
 	Shape TheShape;
@@ -4536,7 +4565,7 @@ case 132:
     ;
     break;}
 case 133:
-#line 1820 "Gmsh.y"
+#line 1849 "Gmsh.y"
 {
       for(int i = 0; i < 4; i++)
 	VisibilityShape(yyvsp[-1].c, i, 1);
@@ -4544,7 +4573,7 @@ case 133:
     ;
     break;}
 case 134:
-#line 1826 "Gmsh.y"
+#line 1855 "Gmsh.y"
 {
       for(int i = 0; i < 4; i++)
 	VisibilityShape(yyvsp[-1].c, i, 0);
@@ -4552,7 +4581,7 @@ case 134:
     ;
     break;}
 case 135:
-#line 1832 "Gmsh.y"
+#line 1861 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
 	Shape TheShape;
@@ -4563,7 +4592,7 @@ case 135:
     ;
     break;}
 case 136:
-#line 1841 "Gmsh.y"
+#line 1870 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
 	Shape TheShape;
@@ -4574,7 +4603,7 @@ case 136:
     ;
     break;}
 case 137:
-#line 1855 "Gmsh.y"
+#line 1884 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-2].c, "Include")){
 	char tmpstring[1024];
@@ -4626,7 +4655,7 @@ case 137:
     ;
     break;}
 case 138:
-#line 1905 "Gmsh.y"
+#line 1934 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-6].c, "Save") && !strcmp(yyvsp[-5].c, "View")){
 	Post_View **vv = (Post_View **)List_Pointer_Test(CTX.post.list, (int)yyvsp[-3].d);
@@ -4643,7 +4672,7 @@ case 138:
     ;
     break;}
 case 139:
-#line 1920 "Gmsh.y"
+#line 1949 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-6].c, "Background") && !strcmp(yyvsp[-5].c, "Mesh")  && !strcmp(yyvsp[-4].c, "View")){
 	Post_View **vv = (Post_View **)List_Pointer_Test(CTX.post.list, (int)yyvsp[-2].d);
@@ -4656,7 +4685,7 @@ case 139:
     ;
     break;}
 case 140:
-#line 1931 "Gmsh.y"
+#line 1960 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-2].c, "Sleep")){
 	SleepInSeconds(yyvsp[-1].d);
@@ -4678,7 +4707,7 @@ case 140:
     ;
     break;}
 case 141:
-#line 1951 "Gmsh.y"
+#line 1980 "Gmsh.y"
 {
        try {
 	 GMSH_PluginManager::instance()->action(yyvsp[-4].c, yyvsp[-1].c, 0);
@@ -4690,7 +4719,7 @@ case 141:
      ;
     break;}
 case 142:
-#line 1961 "Gmsh.y"
+#line 1990 "Gmsh.y"
 {
       if(!strcmp(yyvsp[-1].c, "ElementsFromAllViews"))
 	CombineViews(0, 1, CTX.post.combine_remove_orig);
@@ -4714,27 +4743,27 @@ case 142:
     ;
     break;}
 case 143:
-#line 1983 "Gmsh.y"
+#line 2012 "Gmsh.y"
 {
       exit(0);
     ;
     break;}
 case 144:
-#line 1987 "Gmsh.y"
+#line 2016 "Gmsh.y"
 {
       CTX.forced_bbox = 0;
       SetBoundingBox();
     ;
     break;}
 case 145:
-#line 1992 "Gmsh.y"
+#line 2021 "Gmsh.y"
 {
       CTX.forced_bbox = 1;
       SetBoundingBox(yyvsp[-12].d, yyvsp[-10].d, yyvsp[-8].d, yyvsp[-6].d, yyvsp[-4].d, yyvsp[-2].d);
     ;
     break;}
 case 146:
-#line 1997 "Gmsh.y"
+#line 2026 "Gmsh.y"
 {
 #if defined(HAVE_FLTK)
       Draw();
@@ -4742,7 +4771,7 @@ case 146:
     ;
     break;}
 case 147:
-#line 2009 "Gmsh.y"
+#line 2038 "Gmsh.y"
 {
       LoopControlVariablesTab[ImbricatedLoop][0] = yyvsp[-3].d;
       LoopControlVariablesTab[ImbricatedLoop][1] = yyvsp[-1].d;
@@ -4759,7 +4788,7 @@ case 147:
     ;
     break;}
 case 148:
-#line 2024 "Gmsh.y"
+#line 2053 "Gmsh.y"
 {
       LoopControlVariablesTab[ImbricatedLoop][0] = yyvsp[-5].d;
       LoopControlVariablesTab[ImbricatedLoop][1] = yyvsp[-3].d;
@@ -4777,7 +4806,7 @@ case 148:
     ;
     break;}
 case 149:
-#line 2040 "Gmsh.y"
+#line 2069 "Gmsh.y"
 {
       LoopControlVariablesTab[ImbricatedLoop][0] = yyvsp[-3].d;
       LoopControlVariablesTab[ImbricatedLoop][1] = yyvsp[-1].d;
@@ -4805,7 +4834,7 @@ case 149:
     ;
     break;}
 case 150:
-#line 2066 "Gmsh.y"
+#line 2095 "Gmsh.y"
 {
       LoopControlVariablesTab[ImbricatedLoop][0] = yyvsp[-5].d;
       LoopControlVariablesTab[ImbricatedLoop][1] = yyvsp[-3].d;
@@ -4834,7 +4863,7 @@ case 150:
     ;
     break;}
 case 151:
-#line 2093 "Gmsh.y"
+#line 2122 "Gmsh.y"
 {
       if(ImbricatedLoop <= 0){
 	yymsg(GERROR, "Invalid For/EndFor loop");
@@ -4868,7 +4897,7 @@ case 151:
     ;
     break;}
 case 152:
-#line 2125 "Gmsh.y"
+#line 2154 "Gmsh.y"
 {
       if(!FunctionManager::Instance()->createFunction(yyvsp[0].c, yyin, yyname, yylineno))
 	yymsg(GERROR, "Redefinition of function %s", yyvsp[0].c);
@@ -4877,14 +4906,14 @@ case 152:
     ;
     break;}
 case 153:
-#line 2132 "Gmsh.y"
+#line 2161 "Gmsh.y"
 {
       if(!FunctionManager::Instance()->leaveFunction(&yyin, yyname, yylineno))
 	yymsg(GERROR, "Error while exiting function");
     ;
     break;}
 case 154:
-#line 2137 "Gmsh.y"
+#line 2166 "Gmsh.y"
 {
       if(!FunctionManager::Instance()->enterFunction(yyvsp[-1].c, &yyin, yyname, yylineno))
 	yymsg(GERROR, "Unknown function %s", yyvsp[-1].c);
@@ -4892,18 +4921,18 @@ case 154:
     ;
     break;}
 case 155:
-#line 2143 "Gmsh.y"
+#line 2172 "Gmsh.y"
 {
       if(!yyvsp[-1].d) skip_until("If", "EndIf");
     ;
     break;}
 case 156:
-#line 2147 "Gmsh.y"
+#line 2176 "Gmsh.y"
 {
     ;
     break;}
 case 157:
-#line 2156 "Gmsh.y"
+#line 2185 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, yyvsp[-1].l, 
@@ -4913,7 +4942,7 @@ case 157:
     ;
     break;}
 case 158:
-#line 2164 "Gmsh.y"
+#line 2193 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, yyvsp[-1].l, 
@@ -4923,7 +4952,7 @@ case 158:
     ;
     break;}
 case 159:
-#line 2172 "Gmsh.y"
+#line 2201 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, yyvsp[-1].l, 
@@ -4933,13 +4962,13 @@ case 159:
     ;
     break;}
 case 160:
-#line 2180 "Gmsh.y"
+#line 2209 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 161:
-#line 2184 "Gmsh.y"
+#line 2213 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, yyvsp[-3].l, 
@@ -4949,13 +4978,13 @@ case 161:
     ;
     break;}
 case 162:
-#line 2192 "Gmsh.y"
+#line 2221 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 163:
-#line 2196 "Gmsh.y"
+#line 2225 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, yyvsp[-3].l, 
@@ -4965,13 +4994,13 @@ case 163:
     ;
     break;}
 case 164:
-#line 2204 "Gmsh.y"
+#line 2233 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 165:
-#line 2208 "Gmsh.y"
+#line 2237 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, yyvsp[-3].l, 
@@ -4981,13 +5010,13 @@ case 165:
     ;
     break;}
 case 166:
-#line 2216 "Gmsh.y"
+#line 2245 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 167:
-#line 2220 "Gmsh.y"
+#line 2249 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(BOUNDARY_LAYER, yyvsp[-3].l, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -4996,13 +5025,13 @@ case 167:
     ;
     break;}
 case 168:
-#line 2227 "Gmsh.y"
+#line 2256 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 169:
-#line 2231 "Gmsh.y"
+#line 2260 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       extr.mesh.ViewIndex = yyvsp[-6].d;
@@ -5014,7 +5043,7 @@ case 169:
     ;
     break;}
 case 170:
-#line 2243 "Gmsh.y"
+#line 2272 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)yyvsp[-4].d, 
@@ -5023,7 +5052,7 @@ case 170:
     ;
     break;}
 case 171:
-#line 2250 "Gmsh.y"
+#line 2279 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)yyvsp[-4].d, 
@@ -5032,7 +5061,7 @@ case 171:
     ;
     break;}
 case 172:
-#line 2257 "Gmsh.y"
+#line 2286 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)yyvsp[-4].d, 
@@ -5041,7 +5070,7 @@ case 172:
     ;
     break;}
 case 173:
-#line 2264 "Gmsh.y"
+#line 2293 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)yyvsp[-8].d, 
@@ -5050,7 +5079,7 @@ case 173:
     ;
     break;}
 case 174:
-#line 2271 "Gmsh.y"
+#line 2300 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)yyvsp[-8].d, 
@@ -5059,7 +5088,7 @@ case 174:
     ;
     break;}
 case 175:
-#line 2278 "Gmsh.y"
+#line 2307 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)yyvsp[-8].d, 
@@ -5068,7 +5097,7 @@ case 175:
     ;
     break;}
 case 176:
-#line 2285 "Gmsh.y"
+#line 2314 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)yyvsp[-10].d, 
@@ -5077,7 +5106,7 @@ case 176:
     ;
     break;}
 case 177:
-#line 2292 "Gmsh.y"
+#line 2321 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)yyvsp[-10].d, 
@@ -5086,7 +5115,7 @@ case 177:
     ;
     break;}
 case 178:
-#line 2299 "Gmsh.y"
+#line 2328 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)yyvsp[-10].d, 
@@ -5095,13 +5124,13 @@ case 178:
     ;
     break;}
 case 179:
-#line 2306 "Gmsh.y"
+#line 2335 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 180:
-#line 2310 "Gmsh.y"
+#line 2339 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)yyvsp[-8].d, 
@@ -5110,13 +5139,13 @@ case 180:
     ;
     break;}
 case 181:
-#line 2317 "Gmsh.y"
+#line 2346 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 182:
-#line 2321 "Gmsh.y"
+#line 2350 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)yyvsp[-8].d, 
@@ -5125,13 +5154,13 @@ case 182:
     ;
     break;}
 case 183:
-#line 2328 "Gmsh.y"
+#line 2357 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 184:
-#line 2332 "Gmsh.y"
+#line 2361 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)yyvsp[-8].d, 
@@ -5140,13 +5169,13 @@ case 184:
     ;
     break;}
 case 185:
-#line 2339 "Gmsh.y"
+#line 2368 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 186:
-#line 2343 "Gmsh.y"
+#line 2372 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)yyvsp[-12].d, 
@@ -5155,13 +5184,13 @@ case 186:
     ;
     break;}
 case 187:
-#line 2350 "Gmsh.y"
+#line 2379 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 188:
-#line 2354 "Gmsh.y"
+#line 2383 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)yyvsp[-12].d, 
@@ -5170,13 +5199,13 @@ case 188:
     ;
     break;}
 case 189:
-#line 2361 "Gmsh.y"
+#line 2390 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 190:
-#line 2365 "Gmsh.y"
+#line 2394 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)yyvsp[-12].d, 
@@ -5185,13 +5214,13 @@ case 190:
     ;
     break;}
 case 191:
-#line 2372 "Gmsh.y"
+#line 2401 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 192:
-#line 2376 "Gmsh.y"
+#line 2405 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)yyvsp[-14].d, 
@@ -5200,13 +5229,13 @@ case 192:
     ;
     break;}
 case 193:
-#line 2383 "Gmsh.y"
+#line 2412 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 194:
-#line 2387 "Gmsh.y"
+#line 2416 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)yyvsp[-14].d, 
@@ -5215,13 +5244,13 @@ case 194:
     ;
     break;}
 case 195:
-#line 2394 "Gmsh.y"
+#line 2423 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;
     break;}
 case 196:
-#line 2398 "Gmsh.y"
+#line 2427 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)yyvsp[-14].d, 
@@ -5230,17 +5259,17 @@ case 196:
     ;
     break;}
 case 197:
-#line 2409 "Gmsh.y"
+#line 2438 "Gmsh.y"
 {
     ;
     break;}
 case 198:
-#line 2412 "Gmsh.y"
+#line 2441 "Gmsh.y"
 {
     ;
     break;}
 case 199:
-#line 2418 "Gmsh.y"
+#line 2447 "Gmsh.y"
 {
       extr.mesh.ExtrudeMesh = true;
       extr.mesh.NbLayer = 1;
@@ -5251,7 +5280,7 @@ case 199:
     ;
     break;}
 case 200:
-#line 2427 "Gmsh.y"
+#line 2456 "Gmsh.y"
 {
       double d;
       extr.mesh.ExtrudeMesh = true;
@@ -5275,7 +5304,7 @@ case 200:
     ;
     break;}
 case 201:
-#line 2449 "Gmsh.y"
+#line 2478 "Gmsh.y"
 {
       yymsg(WARNING, "Explicit region numbers in layers are deprecated");
       double d;
@@ -5301,13 +5330,13 @@ case 201:
     ;
     break;}
 case 202:
-#line 2473 "Gmsh.y"
+#line 2502 "Gmsh.y"
 {
       extr.mesh.Recombine = true;
     ;
     break;}
 case 203:
-#line 2477 "Gmsh.y"
+#line 2506 "Gmsh.y"
 {
       int num = (int)yyvsp[-6].d;
       if(FindSurface(num)){
@@ -5328,7 +5357,7 @@ case 203:
     ;
     break;}
 case 204:
-#line 2501 "Gmsh.y"
+#line 2530 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-3].l); i++){
 	double d;
@@ -5348,7 +5377,7 @@ case 204:
     ;
     break;}
 case 205:
-#line 2519 "Gmsh.y"
+#line 2548 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-6].l); i++){
 	double d;
@@ -5368,7 +5397,7 @@ case 205:
     ;
     break;}
 case 206:
-#line 2537 "Gmsh.y"
+#line 2566 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-6].l); i++){
 	double d;
@@ -5388,7 +5417,7 @@ case 206:
     ;
     break;}
 case 207:
-#line 2555 "Gmsh.y"
+#line 2584 "Gmsh.y"
 {
       Surface *s = FindSurface((int)yyvsp[-4].d);
       if(!s)
@@ -5419,7 +5448,7 @@ case 207:
     ;
     break;}
 case 208:
-#line 2584 "Gmsh.y"
+#line 2613 "Gmsh.y"
 {
       Surface *s = FindSurface((int)yyvsp[-5].d);
       if(!s)
@@ -5456,14 +5485,14 @@ case 208:
     ;
     break;}
 case 209:
-#line 2619 "Gmsh.y"
+#line 2648 "Gmsh.y"
 {
       yymsg(WARNING, "Elliptic Surface is deprecated: use Transfinite instead (with smoothing)");
       List_Delete(yyvsp[-1].l);
     ;
     break;}
 case 210:
-#line 2624 "Gmsh.y"
+#line 2653 "Gmsh.y"
 {
       Volume *v = FindVolume((int)yyvsp[-4].d);
       if(!v)
@@ -5492,7 +5521,7 @@ case 210:
     ;
     break;}
 case 211:
-#line 2651 "Gmsh.y"
+#line 2680 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-3].l); i++){
 	double d;
@@ -5508,7 +5537,7 @@ case 211:
     ;
     break;}
 case 212:
-#line 2665 "Gmsh.y"
+#line 2694 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[-1].l); i++){
 	double d;
@@ -5523,7 +5552,7 @@ case 212:
     ;
     break;}
 case 213:
-#line 2684 "Gmsh.y"
+#line 2713 "Gmsh.y"
 { 
       Surface *s = FindSurface((int)yyvsp[-2].d);
       if(s)
@@ -5531,7 +5560,7 @@ case 213:
     ;
     break;}
 case 214:
-#line 2690 "Gmsh.y"
+#line 2719 "Gmsh.y"
 {
       Surface *s = FindSurface((int)yyvsp[-2].d);
       if(s)
@@ -5539,55 +5568,55 @@ case 214:
     ;
     break;}
 case 215:
-#line 2696 "Gmsh.y"
+#line 2725 "Gmsh.y"
 {
     ;
     break;}
 case 216:
-#line 2699 "Gmsh.y"
+#line 2728 "Gmsh.y"
 {
     ;
     break;}
 case 217:
-#line 2708 "Gmsh.y"
+#line 2737 "Gmsh.y"
 { 
       ReplaceAllDuplicates();
     ;
     break;}
 case 218:
-#line 2717 "Gmsh.y"
+#line 2746 "Gmsh.y"
 { yyval.d = yyvsp[0].d;           ;
     break;}
 case 219:
-#line 2718 "Gmsh.y"
+#line 2747 "Gmsh.y"
 { yyval.d = yyvsp[-1].d;           ;
     break;}
 case 220:
-#line 2719 "Gmsh.y"
+#line 2748 "Gmsh.y"
 { yyval.d = -yyvsp[0].d;          ;
     break;}
 case 221:
-#line 2720 "Gmsh.y"
+#line 2749 "Gmsh.y"
 { yyval.d = yyvsp[0].d;           ;
     break;}
 case 222:
-#line 2721 "Gmsh.y"
+#line 2750 "Gmsh.y"
 { yyval.d = !yyvsp[0].d;          ;
     break;}
 case 223:
-#line 2722 "Gmsh.y"
+#line 2751 "Gmsh.y"
 { yyval.d = yyvsp[-2].d - yyvsp[0].d;      ;
     break;}
 case 224:
-#line 2723 "Gmsh.y"
+#line 2752 "Gmsh.y"
 { yyval.d = yyvsp[-2].d + yyvsp[0].d;      ;
     break;}
 case 225:
-#line 2724 "Gmsh.y"
+#line 2753 "Gmsh.y"
 { yyval.d = yyvsp[-2].d * yyvsp[0].d;      ;
     break;}
 case 226:
-#line 2726 "Gmsh.y"
+#line 2755 "Gmsh.y"
 { 
       if(!yyvsp[0].d)
 	yymsg(GERROR, "Division by zero in '%g / %g'", yyvsp[-2].d, yyvsp[0].d);
@@ -5596,247 +5625,247 @@ case 226:
     ;
     break;}
 case 227:
-#line 2732 "Gmsh.y"
+#line 2761 "Gmsh.y"
 { yyval.d = (int)yyvsp[-2].d % (int)yyvsp[0].d;  ;
     break;}
 case 228:
-#line 2733 "Gmsh.y"
+#line 2762 "Gmsh.y"
 { yyval.d = pow(yyvsp[-2].d, yyvsp[0].d);  ;
     break;}
 case 229:
-#line 2734 "Gmsh.y"
+#line 2763 "Gmsh.y"
 { yyval.d = yyvsp[-2].d < yyvsp[0].d;      ;
     break;}
 case 230:
-#line 2735 "Gmsh.y"
+#line 2764 "Gmsh.y"
 { yyval.d = yyvsp[-2].d > yyvsp[0].d;      ;
     break;}
 case 231:
-#line 2736 "Gmsh.y"
+#line 2765 "Gmsh.y"
 { yyval.d = yyvsp[-2].d <= yyvsp[0].d;     ;
     break;}
 case 232:
-#line 2737 "Gmsh.y"
+#line 2766 "Gmsh.y"
 { yyval.d = yyvsp[-2].d >= yyvsp[0].d;     ;
     break;}
 case 233:
-#line 2738 "Gmsh.y"
+#line 2767 "Gmsh.y"
 { yyval.d = yyvsp[-2].d == yyvsp[0].d;     ;
     break;}
 case 234:
-#line 2739 "Gmsh.y"
+#line 2768 "Gmsh.y"
 { yyval.d = yyvsp[-2].d != yyvsp[0].d;     ;
     break;}
 case 235:
-#line 2740 "Gmsh.y"
+#line 2769 "Gmsh.y"
 { yyval.d = yyvsp[-2].d && yyvsp[0].d;     ;
     break;}
 case 236:
-#line 2741 "Gmsh.y"
+#line 2770 "Gmsh.y"
 { yyval.d = yyvsp[-2].d || yyvsp[0].d;     ;
     break;}
 case 237:
-#line 2742 "Gmsh.y"
+#line 2771 "Gmsh.y"
 { yyval.d = yyvsp[-4].d? yyvsp[-2].d : yyvsp[0].d;  ;
     break;}
 case 238:
-#line 2743 "Gmsh.y"
+#line 2772 "Gmsh.y"
 { yyval.d = exp(yyvsp[-1].d);      ;
     break;}
 case 239:
-#line 2744 "Gmsh.y"
+#line 2773 "Gmsh.y"
 { yyval.d = log(yyvsp[-1].d);      ;
     break;}
 case 240:
-#line 2745 "Gmsh.y"
+#line 2774 "Gmsh.y"
 { yyval.d = log10(yyvsp[-1].d);    ;
     break;}
 case 241:
-#line 2746 "Gmsh.y"
+#line 2775 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-1].d);     ;
     break;}
 case 242:
-#line 2747 "Gmsh.y"
+#line 2776 "Gmsh.y"
 { yyval.d = sin(yyvsp[-1].d);      ;
     break;}
 case 243:
-#line 2748 "Gmsh.y"
+#line 2777 "Gmsh.y"
 { yyval.d = asin(yyvsp[-1].d);     ;
     break;}
 case 244:
-#line 2749 "Gmsh.y"
+#line 2778 "Gmsh.y"
 { yyval.d = cos(yyvsp[-1].d);      ;
     break;}
 case 245:
-#line 2750 "Gmsh.y"
+#line 2779 "Gmsh.y"
 { yyval.d = acos(yyvsp[-1].d);     ;
     break;}
 case 246:
-#line 2751 "Gmsh.y"
+#line 2780 "Gmsh.y"
 { yyval.d = tan(yyvsp[-1].d);      ;
     break;}
 case 247:
-#line 2752 "Gmsh.y"
+#line 2781 "Gmsh.y"
 { yyval.d = atan(yyvsp[-1].d);     ;
     break;}
 case 248:
-#line 2753 "Gmsh.y"
+#line 2782 "Gmsh.y"
 { yyval.d = atan2(yyvsp[-3].d, yyvsp[-1].d);;
     break;}
 case 249:
-#line 2754 "Gmsh.y"
+#line 2783 "Gmsh.y"
 { yyval.d = sinh(yyvsp[-1].d);     ;
     break;}
 case 250:
-#line 2755 "Gmsh.y"
+#line 2784 "Gmsh.y"
 { yyval.d = cosh(yyvsp[-1].d);     ;
     break;}
 case 251:
-#line 2756 "Gmsh.y"
+#line 2785 "Gmsh.y"
 { yyval.d = tanh(yyvsp[-1].d);     ;
     break;}
 case 252:
-#line 2757 "Gmsh.y"
+#line 2786 "Gmsh.y"
 { yyval.d = fabs(yyvsp[-1].d);     ;
     break;}
 case 253:
-#line 2758 "Gmsh.y"
+#line 2787 "Gmsh.y"
 { yyval.d = floor(yyvsp[-1].d);    ;
     break;}
 case 254:
-#line 2759 "Gmsh.y"
+#line 2788 "Gmsh.y"
 { yyval.d = ceil(yyvsp[-1].d);     ;
     break;}
 case 255:
-#line 2760 "Gmsh.y"
+#line 2789 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 256:
-#line 2761 "Gmsh.y"
+#line 2790 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 257:
-#line 2762 "Gmsh.y"
+#line 2791 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-3].d*yyvsp[-3].d+yyvsp[-1].d*yyvsp[-1].d); ;
     break;}
 case 258:
-#line 2763 "Gmsh.y"
+#line 2792 "Gmsh.y"
 { yyval.d = yyvsp[-1].d*(double)rand()/(double)RAND_MAX; ;
     break;}
 case 259:
-#line 2765 "Gmsh.y"
+#line 2794 "Gmsh.y"
 { yyval.d = exp(yyvsp[-1].d);      ;
     break;}
 case 260:
-#line 2766 "Gmsh.y"
+#line 2795 "Gmsh.y"
 { yyval.d = log(yyvsp[-1].d);      ;
     break;}
 case 261:
-#line 2767 "Gmsh.y"
+#line 2796 "Gmsh.y"
 { yyval.d = log10(yyvsp[-1].d);    ;
     break;}
 case 262:
-#line 2768 "Gmsh.y"
+#line 2797 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-1].d);     ;
     break;}
 case 263:
-#line 2769 "Gmsh.y"
+#line 2798 "Gmsh.y"
 { yyval.d = sin(yyvsp[-1].d);      ;
     break;}
 case 264:
-#line 2770 "Gmsh.y"
+#line 2799 "Gmsh.y"
 { yyval.d = asin(yyvsp[-1].d);     ;
     break;}
 case 265:
-#line 2771 "Gmsh.y"
+#line 2800 "Gmsh.y"
 { yyval.d = cos(yyvsp[-1].d);      ;
     break;}
 case 266:
-#line 2772 "Gmsh.y"
+#line 2801 "Gmsh.y"
 { yyval.d = acos(yyvsp[-1].d);     ;
     break;}
 case 267:
-#line 2773 "Gmsh.y"
+#line 2802 "Gmsh.y"
 { yyval.d = tan(yyvsp[-1].d);      ;
     break;}
 case 268:
-#line 2774 "Gmsh.y"
+#line 2803 "Gmsh.y"
 { yyval.d = atan(yyvsp[-1].d);     ;
     break;}
 case 269:
-#line 2775 "Gmsh.y"
+#line 2804 "Gmsh.y"
 { yyval.d = atan2(yyvsp[-3].d, yyvsp[-1].d);;
     break;}
 case 270:
-#line 2776 "Gmsh.y"
+#line 2805 "Gmsh.y"
 { yyval.d = sinh(yyvsp[-1].d);     ;
     break;}
 case 271:
-#line 2777 "Gmsh.y"
+#line 2806 "Gmsh.y"
 { yyval.d = cosh(yyvsp[-1].d);     ;
     break;}
 case 272:
-#line 2778 "Gmsh.y"
+#line 2807 "Gmsh.y"
 { yyval.d = tanh(yyvsp[-1].d);     ;
     break;}
 case 273:
-#line 2779 "Gmsh.y"
+#line 2808 "Gmsh.y"
 { yyval.d = fabs(yyvsp[-1].d);     ;
     break;}
 case 274:
-#line 2780 "Gmsh.y"
+#line 2809 "Gmsh.y"
 { yyval.d = floor(yyvsp[-1].d);    ;
     break;}
 case 275:
-#line 2781 "Gmsh.y"
+#line 2810 "Gmsh.y"
 { yyval.d = ceil(yyvsp[-1].d);     ;
     break;}
 case 276:
-#line 2782 "Gmsh.y"
+#line 2811 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 277:
-#line 2783 "Gmsh.y"
+#line 2812 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 278:
-#line 2784 "Gmsh.y"
+#line 2813 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-3].d*yyvsp[-3].d+yyvsp[-1].d*yyvsp[-1].d); ;
     break;}
 case 279:
-#line 2785 "Gmsh.y"
+#line 2814 "Gmsh.y"
 { yyval.d = yyvsp[-1].d*(double)rand()/(double)RAND_MAX; ;
     break;}
 case 280:
-#line 2794 "Gmsh.y"
+#line 2823 "Gmsh.y"
 { yyval.d = yyvsp[0].d; ;
     break;}
 case 281:
-#line 2795 "Gmsh.y"
+#line 2824 "Gmsh.y"
 { yyval.d = 3.141592653589793; ;
     break;}
 case 282:
-#line 2796 "Gmsh.y"
+#line 2825 "Gmsh.y"
 { yyval.d = ParUtil::Instance()->rank(); ;
     break;}
 case 283:
-#line 2797 "Gmsh.y"
+#line 2826 "Gmsh.y"
 { yyval.d = ParUtil::Instance()->size(); ;
     break;}
 case 284:
-#line 2798 "Gmsh.y"
+#line 2827 "Gmsh.y"
 { yyval.d = Get_GmshMajorVersion(); ;
     break;}
 case 285:
-#line 2799 "Gmsh.y"
+#line 2828 "Gmsh.y"
 { yyval.d = Get_GmshMinorVersion(); ;
     break;}
 case 286:
-#line 2800 "Gmsh.y"
+#line 2829 "Gmsh.y"
 { yyval.d = Get_GmshPatchVersion(); ;
     break;}
 case 287:
-#line 2805 "Gmsh.y"
+#line 2834 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[0].c;
@@ -5851,7 +5880,7 @@ case 287:
     ;
     break;}
 case 288:
-#line 2821 "Gmsh.y"
+#line 2850 "Gmsh.y"
 {
       char tmpstring[1024];
       sprintf(tmpstring, "%s_%d", yyvsp[-4].c, (int)yyvsp[-1].d) ;
@@ -5868,7 +5897,7 @@ case 288:
     ;
     break;}
 case 289:
-#line 2836 "Gmsh.y"
+#line 2865 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-3].c;
@@ -5890,7 +5919,7 @@ case 289:
     ;
     break;}
 case 290:
-#line 2856 "Gmsh.y"
+#line 2885 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-2].c;
@@ -5906,7 +5935,7 @@ case 290:
     ;
     break;}
 case 291:
-#line 2870 "Gmsh.y"
+#line 2899 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-1].c;
@@ -5921,7 +5950,7 @@ case 291:
     ;
     break;}
 case 292:
-#line 2883 "Gmsh.y"
+#line 2912 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-4].c;
@@ -5943,7 +5972,7 @@ case 292:
     ;
     break;}
 case 293:
-#line 2906 "Gmsh.y"
+#line 2935 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -5963,7 +5992,7 @@ case 293:
     ;
     break;}
 case 294:
-#line 2924 "Gmsh.y"
+#line 2953 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -5983,7 +6012,7 @@ case 294:
     ;
     break;}
 case 295:
-#line 2942 "Gmsh.y"
+#line 2971 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -6003,7 +6032,7 @@ case 295:
     ;
     break;}
 case 296:
-#line 2960 "Gmsh.y"
+#line 2989 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -6023,107 +6052,107 @@ case 296:
     ;
     break;}
 case 297:
-#line 2978 "Gmsh.y"
+#line 3007 "Gmsh.y"
 { 
       yyval.d = GetValue(yyvsp[-3].c, yyvsp[-1].d);
       Free(yyvsp[-3].c);
     ;
     break;}
 case 298:
-#line 2986 "Gmsh.y"
+#line 3015 "Gmsh.y"
 {
       memcpy(yyval.v, yyvsp[0].v, 5*sizeof(double));
     ;
     break;}
 case 299:
-#line 2990 "Gmsh.y"
+#line 3019 "Gmsh.y"
 {
       for(int i = 0; i < 5; i++) yyval.v[i] = -yyvsp[0].v[i];
     ;
     break;}
 case 300:
-#line 2994 "Gmsh.y"
+#line 3023 "Gmsh.y"
 { 
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[0].v[i];
     ;
     break;}
 case 301:
-#line 2998 "Gmsh.y"
+#line 3027 "Gmsh.y"
 { 
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[-2].v[i] - yyvsp[0].v[i];
     ;
     break;}
 case 302:
-#line 3002 "Gmsh.y"
+#line 3031 "Gmsh.y"
 {
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[-2].v[i] + yyvsp[0].v[i];
     ;
     break;}
 case 303:
-#line 3009 "Gmsh.y"
+#line 3038 "Gmsh.y"
 { 
       yyval.v[0] = yyvsp[-9].d;  yyval.v[1] = yyvsp[-7].d;  yyval.v[2] = yyvsp[-5].d;  yyval.v[3] = yyvsp[-3].d; yyval.v[4] = yyvsp[-1].d;
     ;
     break;}
 case 304:
-#line 3013 "Gmsh.y"
+#line 3042 "Gmsh.y"
 { 
       yyval.v[0] = yyvsp[-7].d;  yyval.v[1] = yyvsp[-5].d;  yyval.v[2] = yyvsp[-3].d;  yyval.v[3] = yyvsp[-1].d; yyval.v[4] = 1.0;
     ;
     break;}
 case 305:
-#line 3017 "Gmsh.y"
+#line 3046 "Gmsh.y"
 {
       yyval.v[0] = yyvsp[-5].d;  yyval.v[1] = yyvsp[-3].d;  yyval.v[2] = yyvsp[-1].d;  yyval.v[3] = 0.0; yyval.v[4] = 1.0;
     ;
     break;}
 case 306:
-#line 3021 "Gmsh.y"
+#line 3050 "Gmsh.y"
 {
       yyval.v[0] = yyvsp[-5].d;  yyval.v[1] = yyvsp[-3].d;  yyval.v[2] = yyvsp[-1].d;  yyval.v[3] = 0.0; yyval.v[4] = 1.0;
     ;
     break;}
 case 307:
-#line 3028 "Gmsh.y"
+#line 3057 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(List_T*));
       List_Add(yyval.l, &(yyvsp[0].l));
     ;
     break;}
 case 308:
-#line 3033 "Gmsh.y"
+#line 3062 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].l));
     ;
     break;}
 case 309:
-#line 3041 "Gmsh.y"
+#line 3070 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 310:
-#line 3046 "Gmsh.y"
+#line 3075 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
     ;
     break;}
 case 311:
-#line 3050 "Gmsh.y"
+#line 3079 "Gmsh.y"
 {
       // creates an empty list
       yyval.l = List_Create(2, 1, sizeof(double));
     ;
     break;}
 case 312:
-#line 3055 "Gmsh.y"
+#line 3084 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 313:
-#line 3059 "Gmsh.y"
+#line 3088 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
       for(int i = 0; i < List_Nbr(yyval.l); i++){
@@ -6133,7 +6162,7 @@ case 313:
     ;
     break;}
 case 314:
-#line 3067 "Gmsh.y"
+#line 3096 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
       for(int i = 0; i < List_Nbr(yyval.l); i++){
@@ -6143,7 +6172,7 @@ case 314:
     ;
     break;}
 case 315:
-#line 3078 "Gmsh.y"
+#line 3107 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
       for(int i = 0; i < List_Nbr(yyval.l); i++){
@@ -6153,7 +6182,7 @@ case 315:
     ;
     break;}
 case 316:
-#line 3086 "Gmsh.y"
+#line 3115 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
       for(int i = 0; i < List_Nbr(yyval.l); i++){
@@ -6163,7 +6192,7 @@ case 316:
     ;
     break;}
 case 317:
-#line 3094 "Gmsh.y"
+#line 3123 "Gmsh.y"
 { 
       yyval.l = List_Create(2, 1, sizeof(double)); 
       for(double d = yyvsp[-2].d; (yyvsp[-2].d < yyvsp[0].d) ? (d <= yyvsp[0].d) : (d >= yyvsp[0].d); (yyvsp[-2].d < yyvsp[0].d) ? (d += 1.) : (d -= 1.)) 
@@ -6171,7 +6200,7 @@ case 317:
     ;
     break;}
 case 318:
-#line 3100 "Gmsh.y"
+#line 3129 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double)); 
       if(!yyvsp[0].d || (yyvsp[-4].d < yyvsp[-2].d && yyvsp[0].d < 0) || (yyvsp[-4].d > yyvsp[-2].d && yyvsp[0].d > 0)){
@@ -6184,7 +6213,7 @@ case 318:
    ;
     break;}
 case 319:
-#line 3111 "Gmsh.y"
+#line 3140 "Gmsh.y"
 {
       // Returns the coordinates of a point and fills a list with it.
       // This allows to ensure e.g. that relative point positions are
@@ -6206,7 +6235,7 @@ case 319:
     ;
     break;}
 case 320:
-#line 3131 "Gmsh.y"
+#line 3160 "Gmsh.y"
 {
       yyval.l = List_Create(List_Nbr(yyvsp[0].l), 1, sizeof(double));
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
@@ -6218,7 +6247,7 @@ case 320:
     ;
     break;}
 case 321:
-#line 3141 "Gmsh.y"
+#line 3170 "Gmsh.y"
 {
       yyval.l = List_Create(List_Nbr(yyvsp[0].l), 1, sizeof(double));
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
@@ -6230,7 +6259,7 @@ case 321:
     ;
     break;}
 case 322:
-#line 3151 "Gmsh.y"
+#line 3180 "Gmsh.y"
 {
       yyval.l = List_Create(List_Nbr(yyvsp[0].l), 1, sizeof(double));
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
@@ -6242,7 +6271,7 @@ case 322:
     ;
     break;}
 case 323:
-#line 3161 "Gmsh.y"
+#line 3190 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -6261,7 +6290,7 @@ case 323:
     ;
     break;}
 case 324:
-#line 3178 "Gmsh.y"
+#line 3207 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -6287,26 +6316,26 @@ case 324:
     ;
     break;}
 case 325:
-#line 3205 "Gmsh.y"
+#line 3234 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 326:
-#line 3210 "Gmsh.y"
+#line 3239 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
     ;
     break;}
 case 327:
-#line 3214 "Gmsh.y"
+#line 3243 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 328:
-#line 3218 "Gmsh.y"
+#line 3247 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
 	double d;
@@ -6317,19 +6346,19 @@ case 328:
     ;
     break;}
 case 329:
-#line 3231 "Gmsh.y"
+#line 3260 "Gmsh.y"
 {
       yyval.u = CTX.PACK_COLOR((int)yyvsp[-7].d, (int)yyvsp[-5].d, (int)yyvsp[-3].d, (int)yyvsp[-1].d);
     ;
     break;}
 case 330:
-#line 3235 "Gmsh.y"
+#line 3264 "Gmsh.y"
 {
       yyval.u = CTX.PACK_COLOR((int)yyvsp[-5].d, (int)yyvsp[-3].d, (int)yyvsp[-1].d, 255);
     ;
     break;}
 case 331:
-#line 3247 "Gmsh.y"
+#line 3276 "Gmsh.y"
 {
       int flag;
       yyval.u = Get_ColorForString(ColorString, -1, yyvsp[0].c, &flag);
@@ -6338,7 +6367,7 @@ case 331:
     ;
     break;}
 case 332:
-#line 3254 "Gmsh.y"
+#line 3283 "Gmsh.y"
 {
       unsigned int (*pColOpt)(int num, int action, unsigned int value);
       StringXColor *pColCat;
@@ -6359,13 +6388,13 @@ case 332:
     ;
     break;}
 case 333:
-#line 3276 "Gmsh.y"
+#line 3305 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 334:
-#line 3280 "Gmsh.y"
+#line 3309 "Gmsh.y"
 {
       yyval.l = List_Create(256, 10, sizeof(unsigned int));
       GmshColorTable *ct = Get_ColorTable((int)yyvsp[-3].d);
@@ -6379,26 +6408,26 @@ case 334:
     ;
     break;}
 case 335:
-#line 3295 "Gmsh.y"
+#line 3324 "Gmsh.y"
 {
       yyval.l = List_Create(256, 10, sizeof(unsigned int));
       List_Add(yyval.l, &(yyvsp[0].u));
     ;
     break;}
 case 336:
-#line 3300 "Gmsh.y"
+#line 3329 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].u));
     ;
     break;}
 case 337:
-#line 3307 "Gmsh.y"
+#line 3336 "Gmsh.y"
 {
       yyval.c = yyvsp[0].c;
     ;
     break;}
 case 338:
-#line 3311 "Gmsh.y"
+#line 3340 "Gmsh.y"
 {
       yyval.c = (char *)Malloc(32*sizeof(char));
       time_t now;
@@ -6408,7 +6437,7 @@ case 338:
     ;
     break;}
 case 339:
-#line 3319 "Gmsh.y"
+#line 3348 "Gmsh.y"
 {
       yyval.c = (char *)Malloc((strlen(yyvsp[-3].c)+strlen(yyvsp[-1].c)+1)*sizeof(char));
       strcpy(yyval.c, yyvsp[-3].c);
@@ -6418,7 +6447,7 @@ case 339:
     ;
     break;}
 case 340:
-#line 3327 "Gmsh.y"
+#line 3356 "Gmsh.y"
 {
       yyval.c = (char *)Malloc((strlen(yyvsp[-1].c)+1)*sizeof(char));
       int i;
@@ -6434,7 +6463,7 @@ case 340:
     ;
     break;}
 case 341:
-#line 3341 "Gmsh.y"
+#line 3370 "Gmsh.y"
 {
       yyval.c = (char *)Malloc((strlen(yyvsp[-1].c)+1)*sizeof(char));
       int i;
@@ -6450,13 +6479,13 @@ case 341:
     ;
     break;}
 case 342:
-#line 3355 "Gmsh.y"
+#line 3384 "Gmsh.y"
 {
       yyval.c = yyvsp[-1].c;
     ;
     break;}
 case 343:
-#line 3359 "Gmsh.y"
+#line 3388 "Gmsh.y"
 {
       char tmpstring[1024];
       int i = PrintListOfDouble(yyvsp[-3].c, yyvsp[-1].l, tmpstring);
@@ -6477,7 +6506,7 @@ case 343:
     ;
     break;}
 case 344:
-#line 3378 "Gmsh.y"
+#line 3407 "Gmsh.y"
 { 
       char* (*pStrOpt)(int num, int action, char *value);
       StringXString *pStrCat;
@@ -6501,7 +6530,7 @@ case 344:
     ;
     break;}
 case 345:
-#line 3400 "Gmsh.y"
+#line 3429 "Gmsh.y"
 { 
       char* (*pStrOpt)(int num, int action, char *value);
       StringXString *pStrCat;
@@ -6746,7 +6775,7 @@ yyerrhandle:
     }
   return 1;
 }
-#line 3423 "Gmsh.y"
+#line 3452 "Gmsh.y"
 
 
 void DeleteSymbol(void *a, void *b){
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 0e6f2cb435..e517b37f5d 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1,5 +1,5 @@
 %{
-// $Id: Gmsh.y,v 1.267 2007-03-05 11:01:21 geuzaine Exp $
+// $Id: Gmsh.y,v 1.268 2007-03-11 20:19:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1037,13 +1037,12 @@ Shape :
 	double z = CTX.geom.scaling_factor * $6[2];
 	double lc = CTX.geom.scaling_factor * $6[3];
 	Vertex *v;
-	if (!myGmshSurface)
+	if(!myGmshSurface)
 	  v = Create_Vertex(num, x, y, z, lc, 1.0);
 	else
 	  v = Create_Vertex(num, x, y, myGmshSurface, lc);
-
 	Tree_Add(THEM->Points, &v);
-	AddToTemporaryBoundingBox(v->Pos.X,v->Pos.Y,v->Pos.Z);
+	AddToTemporaryBoundingBox(v->Pos.X, v->Pos.Y, v->Pos.Z);
       }
       $$.Type = MSH_POINT;
       $$.Num = num;
@@ -1081,15 +1080,16 @@ Shape :
 	List_Read($3, i, &d);
 	Vertex *v = FindPoint((int)d); 
 	if(v)
-	  att->addPoint(v->Pos.X,v->Pos.Y,v->Pos.Z);
+	  att->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
 	else{
 	  GVertex *gv = GMODEL->vertexByTag((int)d);
 	  if(gv) 
-	    att->addPoint(gv->x(),gv->y(),gv->z());
+	    att->addPoint(gv->x(), gv->y(), gv->z());
 	}
       }
       att->buildFastSearchStructures();
-      $$.Type = MSH_POINT_ATTRACTOR;
+      // dummy values
+      $$.Type = 0;
       $$.Num = 0;
     }
   | tAttractor tLine ListOfDouble tAFFECT ListOfDouble tEND
@@ -1109,17 +1109,18 @@ Shape :
 	List_Read($3, i, &d);
 	Curve *c = FindCurve((int)d); 
 	if(c){
-	  buildListOfPoints( att , c , (int) pars[3] );
+	  buildListOfPoints(att, c, (int)pars[3]);
 	}
 	else{
 	  GEdge *ge = GMODEL->edgeByTag((int)d);
 	  if(ge){
-	    buildListOfPoints( att , ge , (int) pars[3] );
+	    buildListOfPoints(att, ge, (int)pars[3]);
 	  }
 	}
       }
       att->buildFastSearchStructures();
-      $$.Type = MSH_LINE_ATTRACTOR;
+      // dummy values
+      $$.Type = 0;
       $$.Num = 0;
     }
   | tCharacteristic tLength ListOfDouble tAFFECT FExpr tEND
@@ -1651,12 +1652,19 @@ ListOfShapes :
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Vertex *v = FindPoint(TheShape.Num);
-	if(!v)
-	  yymsg(WARNING, "Unknown point %d", TheShape.Num);
-	else{
+	if(v){
 	  TheShape.Type = MSH_POINT;
 	  List_Add($$, &TheShape);
 	}
+	else{
+	  GVertex *gv = GMODEL->vertexByTag(TheShape.Num);
+	  if(gv){
+	    TheShape.Type = MSH_POINT_FROM_GMODEL;
+	    List_Add($$, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown point %d", TheShape.Num);
+	}
       }
     }
   | ListOfShapes tLine '{' RecursiveListOfDouble '}' tEND
@@ -1667,12 +1675,19 @@ ListOfShapes :
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Curve *c = FindCurve(TheShape.Num);
-	if(!c)
-	  yymsg(WARNING, "Unknown curve %d", TheShape.Num);
-	else{
+	if(c){
 	  TheShape.Type = c->Typ;
 	  List_Add($$, &TheShape);
 	}
+	else{
+	  GEdge *ge = GMODEL->edgeByTag(TheShape.Num);
+	  if(ge){
+	    TheShape.Type = MSH_SEGM_FROM_GMODEL;
+	    List_Add($$, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown curve %d", TheShape.Num);
+	}
       }
     }
   | ListOfShapes tSurface '{' RecursiveListOfDouble '}' tEND
@@ -1683,12 +1698,19 @@ ListOfShapes :
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Surface *s = FindSurface(TheShape.Num);
-	if(!s)
-	  yymsg(WARNING, "Unknown surface %d", TheShape.Num);
-	else{
+	if(s){
 	  TheShape.Type = s->Typ;
 	  List_Add($$, &TheShape);
 	}
+	else{
+	  GFace *gf = GMODEL->faceByTag(TheShape.Num);
+	  if(gf){
+	    TheShape.Type = MSH_SURF_FROM_GMODEL;
+	    List_Add($$, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown surface %d", TheShape.Num);
+	}
       }
     }
   | ListOfShapes tVolume '{' RecursiveListOfDouble '}' tEND
@@ -1699,12 +1721,19 @@ ListOfShapes :
 	Shape TheShape;
 	TheShape.Num = (int)d;
 	Volume *v = FindVolume(TheShape.Num);
-	if(!v)
-	  yymsg(WARNING, "Unknown volume %d", TheShape.Num);
-	else{
+	if(v){
 	  TheShape.Type = v->Typ;
 	  List_Add($$, &TheShape);
 	}
+	else{
+	  GRegion *gr = GMODEL->regionByTag(TheShape.Num);
+	  if(gr){
+	    TheShape.Type = MSH_VOLUME_FROM_GMODEL;
+	    List_Add($$, &TheShape);
+	  }
+	  else
+	    yymsg(WARNING, "Unknown volume %d", TheShape.Num);
+	}
       }
     }
 ;
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index 89abe07bb4..a403cabe43 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -2,7 +2,7 @@
 /* A lexical scanner generated by flex */
 
 /* Scanner skeleton version:
- * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.311 2007-03-09 14:57:08 remacle Exp $
+ * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.312 2007-03-11 20:19:02 geuzaine Exp $
  */
 
 #define FLEX_SCANNER
@@ -740,7 +740,7 @@ char *yytext;
 #line 1 "Gmsh.l"
 #define INITIAL 0
 #line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.311 2007-03-09 14:57:08 remacle Exp $
+// $Id: Gmsh.yy.cpp,v 1.312 2007-03-11 20:19:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 57492fc97e..2680a7b672 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.143 2007-03-05 09:30:57 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.144 2007-03-11 20:19:05 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -350,6 +350,9 @@ int MergeFile(char *name, int warn_if_missing)
 	  !strcmp(ext, ".nas") || !strcmp(ext, ".NAS")){
     status = GMODEL->readBDF(name);
   }
+  else if(!strcmp(ext, ".p3d") || !strcmp(ext, ".P3D")){
+    status = GMODEL->readP3D(name);
+  }
 #if defined(HAVE_FLTK)
   else if(!strcmp(ext, ".pnm") || !strcmp(ext, ".PNM") ||
 	  !strcmp(ext, ".pbm") || !strcmp(ext, ".PBM") ||
diff --git a/benchmarks/3d/sph.geo b/benchmarks/3d/sph.geo
index 0f4a835045..398fa06723 100644
--- a/benchmarks/3d/sph.geo
+++ b/benchmarks/3d/sph.geo
@@ -41,7 +41,7 @@ Complex Volume(28) = {27};
 Transfinite Line {1,2,3,4,5,6} = 10 ;
 Transfinite Line {-8,-10,-12} = 5 Using Power 1.6 ;
 
-Transfinite Surface {22} = {4,7,6,3};
+Transfinite Surface {22} = {4,7,6,3} Right;
 Transfinite Surface {20} = {3,2,5,6};
 Transfinite Surface {18} = {4,2,5,7};
 
@@ -49,8 +49,6 @@ Transfinite Surface {24} = {3,4,2};
 Transfinite Surface {26} = {6,7,5};
 
 Recombine Surface {18,20,22};
-/*
-Recombine Surface {24,26};
-*/
-Transfinite Volume {28} = {3,4,2,6,7,5};
+//Recombine Surface {24,26};
 
+Transfinite Volume {28} = {3,4,2,6,7,5};
diff --git a/benchmarks/3d/transfinite.geo b/benchmarks/3d/transfinite.geo
index 737a6f7217..7d168cc0d2 100644
--- a/benchmarks/3d/transfinite.geo
+++ b/benchmarks/3d/transfinite.geo
@@ -59,9 +59,4 @@ Transfinite Volume {26} = {7,10,2,1,8,9,6,5};
 
 Recombine Surface {6,20,18,22,16,24};
 
-
-
-
-
-
-
+Physical Volume(3333) = 26;
diff --git a/doc/TODO b/doc/TODO
index 2133dd435c..fb5dd37255 100644
--- a/doc/TODO
+++ b/doc/TODO
@@ -1,4 +1,4 @@
-$Id: TODO,v 1.52 2007-03-02 13:59:00 geuzaine Exp $
+$Id: TODO,v 1.53 2007-03-11 20:19:06 geuzaine Exp $
 
 ********************************************************************
 
@@ -15,8 +15,8 @@ reinterface Triangle for plane surfaces
 
 ********************************************************************
 
-Bug: need to better understand how Netgen deals with char lenghts: a
-larger geometry (with appropriately scaled lcs) leads to a larger
+Bug: need to better understand how Netgen deals with charact. lengths:
+a larger geometry (with appropriately scaled lcs) leads to a larger
 mesh!
 
 ********************************************************************
@@ -75,8 +75,9 @@ a surface should add it with reverse orientation?
 
 ********************************************************************
 
-fill transfinite_vertices in meshGFaceExtrude so that we can use
-extruded+recombined surfaces to create Transfinite Volumes?
+fill transfinite_vertices in meshGFaceExtrude/meshGRegionExtrude when
+it makes sense (so that we can use extruded+recombined surfaces to
+create Transfinite Volumes, or use the P3D output format)
 
 ********************************************************************
 
@@ -108,8 +109,7 @@ implement better algo to determine which axes to draw
 
 ********************************************************************
 
-extrude along a given function (e.g. radius to extrude a sphere) or
-along a parametric curve
+extrude along a given function or along a parametric curve
 
 ********************************************************************
 
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 00f59b186f..39d0167d93 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,8 +1,12 @@
-$Id: VERSIONS,v 1.380 2007-02-26 08:25:46 geuzaine Exp $
+$Id: VERSIONS,v 1.381 2007-03-11 20:19:06 geuzaine Exp $
 
 new since 2.0: volumes can now be defined from external CAD surfaces;
-Delaunay/Tetgen algorithm is now used by default when available; fixed
-various bugs.
+Delaunay/Tetgen algorithm is now used by default when available;
+re-added support for Plot3D structured mesh format; added ability to
+export external CAD models as GEO files (this only works for the
+limited set of geometrical primitives available in the GEO language,
+of course--so trying to convert e.g. a trimmed NURBS from a STEP file
+into a GEO file will fail); fixed various bugs.
 
 2.0 (February 5, 2007): new geometry and mesh databases, with support
 for STEP and IGES import via OpenCascade; complete rewrite of geometry
diff --git a/doc/gmsh.1 b/doc/gmsh.1
index 68192cfe82..0e0c0f1dfd 100644
--- a/doc/gmsh.1
+++ b/doc/gmsh.1
@@ -1,4 +1,4 @@
-.\" $Id: gmsh.1,v 1.78 2007-02-12 08:36:14 geuzaine Exp $
+.\" $Id: gmsh.1,v 1.79 2007-03-11 20:19:06 geuzaine Exp $
 .TH Gmsh 1 "09 March 2006" "Gmsh 2.0" "Gmsh Manual Pages"
 .UC 4
 .\" ********************************************************************
@@ -49,7 +49,7 @@ save all elements (and discard all physical group definitions).
 specify mesh output file name.
 .TP 4
 .B \-format string
-set output mesh format (msh, unv, vrml, stl, mesh, bdf, cgns, med).
+set output mesh format (msh, unv, vrml, stl, mesh, bdf, p3d, cgns, med).
 .TP 4
 .B \-algo string
 select mesh algorithm (iso, netgen, tetgen).
diff --git a/doc/texinfo/command_line.texi b/doc/texinfo/command_line.texi
index eb3f229a69..3d8b1c8496 100644
--- a/doc/texinfo/command_line.texi
+++ b/doc/texinfo/command_line.texi
@@ -19,7 +19,7 @@ Save all elements (discard physical group definitions)
 @item -o file
 Specify mesh output file name
 @item -format string
-Set output mesh format (msh, unv, vrml, stl, mesh, bdf, cgns, med)
+Set output mesh format (msh, unv, vrml, stl, mesh, bdf, p3d, cgns, med)
 @item -algo string
 Select mesh algorithm (iso, netgen, tetgen)
 @item -smooth int
-- 
GitLab