From b1b0c04c3d88731184a63910e585d27d4b3bea10 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 15 Feb 2009 15:14:04 +0000
Subject: [PATCH] - moved writeGEO() into the entity classes so that we can
 exploit   kernel-specific information - better documentation of lists
 returned by extrusions

---
 Geo/GEdge.cpp         |  27 ++++++
 Geo/GEdge.h           |   4 +
 Geo/GFace.cpp         |  32 +++++++
 Geo/GFace.h           |   3 +
 Geo/GModelIO_Geo.cpp  | 197 ++----------------------------------------
 Geo/GRegion.cpp       |  17 ++++
 Geo/GRegion.h         |   4 +
 Geo/GVertex.cpp       |   6 ++
 Geo/GVertex.h         |   8 +-
 Geo/OCCEdge.cpp       |  60 ++++++-------
 Geo/OCCEdge.h         |   1 +
 Geo/OCCFace.cpp       |   2 +-
 Geo/gmshEdge.cpp      |  66 ++++++++++++++
 Geo/gmshEdge.h        |   1 +
 Geo/gmshVertex.cpp    |   6 ++
 Geo/gmshVertex.h      |   1 +
 doc/texinfo/gmsh.texi |  25 +++++-
 17 files changed, 232 insertions(+), 228 deletions(-)

diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index dcec6e5897..635c06edcd 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -122,6 +122,33 @@ std::string GEdge::getAdditionalInfoString()
   return sstream.str();
 }
 
+void GEdge::writeGEO(FILE *fp)
+{
+  if(!getBeginVertex() || !getEndVertex() || geomType() == DiscreteCurve) return;
+
+  if(geomType() == Line){
+    fprintf(fp, "Line(%d) = {%d, %d};\n", 
+            tag(), getBeginVertex()->tag(), getEndVertex()->tag());
+  }
+  else{
+    // approximate other curves by splines
+    Range<double> bounds = parBounds(0);
+    double umin = bounds.low();
+    double umax = bounds.high();
+    fprintf(fp, "p%d = newp;\n", tag());
+    for(int i = 1; i < minimumDrawSegments(); i++){
+      double u = umin + (double)i / minimumDrawSegments() * (umax - umin);
+      GPoint p = point(u);
+      fprintf(fp, "Point(p%d + %d) = {%.16g, %.16g, %.16g};\n", 
+              tag(), i, p.x(), p.y(), p.z());
+    }
+    fprintf(fp, "Spline(%d) = {%d", tag(), getBeginVertex()->tag());
+    for(int i = 1; i < minimumDrawSegments(); i++)
+      fprintf(fp, ", p%d + %d", tag(), i);
+    fprintf(fp, ", %d};\n", getEndVertex()->tag());
+  }
+}
+
 bool GEdge::containsParam(double pt) const
 {
   Range<double> rg = parBounds(0);
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 3dbbefb2b5..1942eb0084 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -6,6 +6,7 @@
 #ifndef _GEDGE_H_
 #define _GEDGE_H_
 
+#include <stdio.h>
 #include "GEntity.h"
 #include "GVertex.h"
 #include "SVector3.h"
@@ -95,6 +96,9 @@ class GEdge : public GEntity {
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
+  // export in GEO format
+  virtual void writeGEO(FILE *fp);
+
   // tell if the edge is a 3D edge (in opposition with a trimmed curve on a surface)
   virtual bool is3D() const { return true; }
 
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index f0454bf757..0bb307b78a 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -196,6 +196,38 @@ std::string GFace::getAdditionalInfoString()
   return sstream.str();
 }
 
+void GFace::writeGEO(FILE *fp)
+{
+  if(geomType() == DiscreteSurface) return;
+
+  std::list<GEdge*> edg = edges();
+  std::list<int> dir = orientations();
+  if(edg.size() && dir.size() == edg.size()){
+    std::vector<int> num, ori;
+    for(std::list<GEdge*>::iterator it = edg.begin(); it != edg.end(); it++)
+      num.push_back((*it)->tag());
+    for(std::list<int>::iterator it = dir.begin(); it != dir.end(); it++)
+      ori.push_back((*it) > 0 ? 1 : -1);
+    fprintf(fp, "Line Loop(%d) = ", tag());
+    for(unsigned int i = 0; i < num.size(); i++){
+      if(i)
+        fprintf(fp, ", %d", num[i] * ori[i]);
+      else
+        fprintf(fp, "{%d", num[i] * ori[i]);
+    }
+    fprintf(fp, "};\n");
+    if(geomType() == GEntity::Plane){
+      fprintf(fp, "Plane Surface(%d) = {%d};\n", tag(), tag());
+    }
+    else if(edg.size() == 3 || edg.size() == 4){
+      fprintf(fp, "Ruled Surface(%d) = {%d};\n", tag(), tag());
+    }
+    else{
+      Msg::Error("Skipping surface %d in export", tag());
+    }
+  }
+}
+
 void GFace::computeMeanPlane()
 {
   std::vector<SPoint3> pts;
diff --git a/Geo/GFace.h b/Geo/GFace.h
index e82ffd58c1..6ed4fea375 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -148,6 +148,9 @@ class GFace : public GEntity
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
+  // export in GEO format
+  virtual void writeGEO(FILE *fp);
+
   // fill the crude representation cross
   virtual bool buildRepresentationCross();
 
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index d5b034f7a1..2345fdb38d 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -194,191 +194,6 @@ class writeFieldGEO {
   }
 };
 
-class writeGVertexGEO {
- private :
-  FILE *geo;
- public :
-  writeGVertexGEO(FILE *fp) { geo = fp ? fp : stdout; }
-  void operator() (GVertex *gv)
-  {
-    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());
-    }
-  }
-};
-
-class writeGEdgeGEO {
- private :
-  FILE *geo;
- public :
-  writeGEdgeGEO(FILE *fp) { geo = fp ? fp : stdout; }
-  void operator () (GEdge *ge)
-  {
-    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, "Spline(%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::Error("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
-          fprintf(geo, "{%d", v->Num);
-        if(i % 6 == 7)
-          fprintf(geo, "\n");
-      }
-      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, "Spline(%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());
-        }
-      }
-    }
-  }
-};
-
-class writeGFaceGEO {
- private :
-  FILE *geo;
- public :
-  writeGFaceGEO(FILE *fp) { geo = fp ? fp : stdout; }
-  void operator () (GFace *gf)
-  {
-    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);
-      fprintf(geo, "Line Loop(%d) = ", gf->tag());
-      for(unsigned int i = 0; i < num.size(); i++){
-        if(i)
-          fprintf(geo, ", %d", num[i] * ori[i]);
-        else
-          fprintf(geo, "{%d", num[i] * ori[i]);
-      }
-      fprintf(geo, "};\n");
-      if(gf->geomType() == GEntity::Plane){
-        fprintf(geo, "Plane Surface(%d) = {%d};\n", gf->tag(), gf->tag());
-      }
-      else if(edges.size() == 3 || edges.size() == 4){
-        fprintf(geo, "Ruled Surface(%d) = {%d};\n", gf->tag(), gf->tag());
-      }
-      else{
-        Msg::Error("Skipping surface %d in export", gf->tag());
-      }
-    }
-  }
-};
-
-class writeGRegionGEO {
- private :
-  FILE *geo;
- public :
-  writeGRegionGEO(FILE *fp) { geo = fp ? fp : stdout; }
-  void operator () (GRegion *gr)
-  {
-    if(gr->geomType() == GEntity::DiscreteVolume) return;
-
-    std::list<GFace*> faces = gr->faces();
-    if(faces.size()){
-      fprintf(geo, "Surface Loop(%d) = ", gr->tag());
-      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(), gr->tag());
-    }
-  }
-};
-
 class writePhysicalGroupGEO {
  private :
   FILE *geo;
@@ -431,10 +246,14 @@ int GModel::writeGEO(const std::string &name, bool printLabels)
 {
   FILE *fp = fopen(name.c_str(), "w");
 
-  std::for_each(firstVertex(), lastVertex(), writeGVertexGEO(fp));
-  std::for_each(firstEdge(), lastEdge(), writeGEdgeGEO(fp));
-  std::for_each(firstFace(), lastFace(), writeGFaceGEO(fp));
-  std::for_each(firstRegion(), lastRegion(), writeGRegionGEO(fp));
+  for(viter it = firstVertex(); it != lastVertex(); it++)
+    (*it)->writeGEO(fp);
+  for(eiter it = firstEdge(); it != lastEdge(); it++)
+    (*it)->writeGEO(fp);
+  for(fiter it = firstFace(); it != lastFace(); it++)
+    (*it)->writeGEO(fp);
+  for(riter it = firstRegion(); it != lastRegion(); it++)
+    (*it)->writeGEO(fp);
 
   std::map<int, std::string> labels;
 #if !defined(HAVE_NO_PARSER)
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 7a9f406e54..dcb1950406 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -136,6 +136,23 @@ std::string GRegion::getAdditionalInfoString()
   return sstream.str();
 }
 
+void GRegion::writeGEO(FILE *fp)
+{
+  if(geomType() == DiscreteVolume) return;
+
+  if(l_faces.size()){
+    fprintf(fp, "Surface Loop(%d) = ", tag());
+    for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); it++) {
+      if(it != l_faces.begin())
+        fprintf(fp, ", %d", (*it)->tag());
+      else
+        fprintf(fp, "{%d", (*it)->tag());
+    }
+    fprintf(fp, "};\n");
+    fprintf(fp, "Volume(%d) = {%d};\n", tag(), tag());
+  }
+}
+
 std::list<GEdge*> GRegion::edges() const
 {
   std::list<GEdge*> e;
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index e97b5057fe..9aeb4ccd38 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -6,6 +6,7 @@
 #ifndef _GREGION_H_
 #define _GREGION_H_
 
+#include <stdio.h>
 #include "GEntity.h"
 
 class MElement;
@@ -51,6 +52,9 @@ class GRegion : public GEntity {
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
+  // export in GEO format
+  virtual void writeGEO(FILE *fp);
+
   // number of types of elements
   int getNumElementTypes() const { return 4; }
 
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index f7796d8249..821a95404f 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -56,6 +56,12 @@ std::string GVertex::getAdditionalInfoString()
   return sstream.str();
 }
 
+void GVertex::writeGEO(FILE *fp)
+{
+  fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+          tag(), x(), y(), z(), prescribedMeshSizeAtVertex());
+}
+
 unsigned int GVertex::getNumMeshElements()
 {
   return points.size(); 
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 5d525c251c..27db6204b6 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -6,6 +6,7 @@
 #ifndef _GVERTEX_H_
 #define _GVERTEX_H_
 
+#include <stdio.h>
 #include "GEntity.h"
 #include "GPoint.h"
 #include "SPoint2.h"
@@ -39,6 +40,9 @@ class GVertex : public GEntity
   void addEdge(GEdge *e);
   void delEdge(GEdge *e);
 
+  // get the edges that this vertex bounds
+  virtual std::list<GEdge*> edges() const{ return l_edges; }
+
   // get the dimension of the vertex (0)
   virtual int dim() const { return 0; }
 
@@ -58,8 +62,8 @@ class GVertex : public GEntity
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
-  // get the edges that this vertex bounds
-   virtual std::list<GEdge*> edges() const{ return l_edges; }
+  // export in GEO format
+  virtual void writeGEO(FILE *fp);
 
   // get number of elements in the mesh
   unsigned int getNumMeshElements();
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 4a0988d0da..2211896ada 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -30,23 +30,6 @@ OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex
   // build the reverse curve
   c_rev = c;
   c_rev.Reverse();
-
-
-//   if (v0 == v1){
-//     const int JJ = 52;
-//     for (int i=1;i<JJ;i++){
-//       const double t = i/((double) JJ);
-//       const double xi = s0 + (s1-s0) * t; 
-//       GPoint p = point(xi);
-//       MEdgeVertex *v = new MEdgeVertex (p.x(),p.y(),p.z(),this,xi);
-//       mesh_vertices.push_back(v);
-//       meshAttributes.Method = MESH_NONE;
-//       if (i == 1)lines.push_back(new MLine(v1->mesh_vertices[0],v));
-//       else if (i == JJ-1)lines.push_back(new MLine(v,v2->mesh_vertices[0]));
-//       else lines.push_back(new MLine(mesh_vertices[i-1],v));
-//     }
-//   }
-
 }
 
 Range<double> OCCEdge::parBounds(int i) const
@@ -54,7 +37,7 @@ Range<double> OCCEdge::parBounds(int i) const
   return Range<double>(s0, s1);
 }
 
-void OCCEdge::setTrimmed (OCCFace *f)
+void OCCEdge::setTrimmed(OCCFace *f)
 {
   if (!trimmed){
     trimmed = f;
@@ -128,7 +111,7 @@ GPoint OCCEdge::point(double par) const
     return trimmed->point(u, v);
   }
   else if(!curve.IsNull()){
-    gp_Pnt pnt = curve->Value (par);
+    gp_Pnt pnt = curve->Value(par);
     return GPoint(pnt.X(), pnt.Y(), pnt.Z());
   }
   else{
@@ -141,7 +124,7 @@ SVector3 OCCEdge::firstDer(double par) const
 {  
   BRepAdaptor_Curve brepc(c);
   BRepLProp_CLProps prop(brepc, 1, 1e-5);
-  prop.SetParameter (par);
+  prop.SetParameter(par);
   gp_Vec d1 = prop.D1();
   return SVector3(d1.X(), d1.Y(), d1.Z());
 }
@@ -243,19 +226,32 @@ double OCCEdge::curvature(double par) const
       Crv = prop.Curvature();
   }
   if(Crv <= eps) Crv = eps;
-  
-  // std::list<GFace*> ff = faces();
-  // std::list<GFace *>::iterator it =  ff.begin();
-  // while (it != ff.end()){
-  //   SPoint2 par2 = reparamOnFace((*it),par,1);
-  //   const double cc = (*it)->curvature ( par2 );
-  //   if (cc > 0)
-  //     Crv = std::max( Crv, cc);  
-  //   ++it;
-  // }  
-  // printf("curvature = %12.5E\n",Crv); 
-
   return Crv;
 }
 
+void OCCEdge::writeGEO(FILE *fp)
+{
+  if(geomType() == Circle){
+    gp_Pnt center;
+    if(curve.IsNull()){
+      center = Handle(Geom_Circle)::DownCast(curve2d)->Location();
+    }
+    else{
+      center = Handle(Geom_Circle)::DownCast(curve)->Location();
+    }
+    // GEO supports only circle arcs < Pi
+    if(s1 - s0 < M_PI){
+      fprintf(fp, "p%d = newp;\n", tag());
+      fprintf(fp, "Point(p%d + 1) = {%.16g, %.16g, %.16g};\n", 
+              tag(), center.X(), center.Y(), center.Z());
+      fprintf(fp, "Circle(%d) = {%d, p%d + 1, %d};\n", 
+              tag(), getBeginVertex()->tag(), tag(), getEndVertex()->tag());
+    }
+    else
+      GEdge::writeGEO(fp);
+  }
+  else
+    GEdge::writeGEO(fp);
+}
+
 #endif
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index d04b7476a9..d757651f26 100644
--- a/Geo/OCCEdge.h
+++ b/Geo/OCCEdge.h
@@ -41,6 +41,7 @@ class OCCEdge : public GEdge {
   bool is3D() const { return !curve.IsNull(); }
   void setTrimmed(OCCFace *);
   bool isSeam(const GFace *) const;
+  virtual void writeGEO(FILE *fp);
 };
 
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 2f7b56e330..9fda93d3a9 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -187,7 +187,7 @@ GEntity::GeomType OCCFace::geomType() const
   return Unknown;
 }
 
-double OCCFace::curvature (const SPoint2 &param) const
+double OCCFace::curvature(const SPoint2 &param) const
 {
   const double eps = 1.e-12;
   BRepAdaptor_Surface sf(s, Standard_True);
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index a476b89d88..2a7f48ae17 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -258,3 +258,69 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const
   else
     return GEdge::reparamOnFace(face, epar, dir);
 }
+
+void gmshEdge::writeGEO(FILE *fp)
+{
+  if(!c || c->Num < 0 || c->Typ == MSH_SEGM_DISCRETE) return;
+  switch (c->Typ) {
+  case MSH_SEGM_LINE:
+    fprintf(fp, "Line(%d) = ", c->Num);
+    break;
+  case MSH_SEGM_CIRC:
+  case MSH_SEGM_CIRC_INV:
+    fprintf(fp, "Circle(%d) = ", c->Num);
+    break;
+  case MSH_SEGM_ELLI:
+  case MSH_SEGM_ELLI_INV:
+    fprintf(fp, "Ellipse(%d) = ", c->Num);
+    break;
+  case MSH_SEGM_NURBS:
+    fprintf(fp, "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(fp, "%d", v->Num);
+      else
+        fprintf(fp, ", %d", v->Num);
+      if(i % 8 == 7 && i != List_Nbr(c->Control_Points) - 1)
+        fprintf(fp, "\n");
+    }
+    fprintf(fp, "}\n");
+    fprintf(fp, "  Knots {");
+    for(int j = 0; j < List_Nbr(c->Control_Points) + c->degre + 1; j++) {
+      if(!j)
+        fprintf(fp, "%.16g", c->k[j]);
+      else
+        fprintf(fp, ", %.16g", c->k[j]);
+      if(j % 5 == 4 && j != List_Nbr(c->Control_Points) + c->degre)
+        fprintf(fp, "\n        ");
+    }
+    fprintf(fp, "}\n");
+    fprintf(fp, "  Order %d;\n", c->degre);
+    return;
+  case MSH_SEGM_SPLN:
+    fprintf(fp, "Spline(%d) = ", c->Num);
+    break;
+  case MSH_SEGM_BSPLN:
+    fprintf(fp, "BSpline(%d) = ", c->Num);
+    break;
+  case MSH_SEGM_BEZIER:
+    fprintf(fp, "Bezier(%d) = ", c->Num);
+    break;
+  default:
+    Msg::Error("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(fp, ", %d", v->Num);
+    else
+      fprintf(fp, "{%d", v->Num);
+    if(i % 6 == 7)
+      fprintf(fp, "\n");
+  }
+  fprintf(fp, "};\n");
+}
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index 0c3d2daf29..72155dbd33 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -26,6 +26,7 @@ class gmshEdge : public GEdge {
   virtual int minimumDrawSegments() const;
   virtual void resetMeshAttributes();
   virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
+  virtual void writeGEO(FILE *fp);
 };
 
 #endif
diff --git a/Geo/gmshVertex.cpp b/Geo/gmshVertex.cpp
index dfa94f0360..7a324e33d1 100644
--- a/Geo/gmshVertex.cpp
+++ b/Geo/gmshVertex.cpp
@@ -97,3 +97,9 @@ SPoint2 gmshVertex::reparamOnFace(const GFace *face, int dir) const
     return GVertex::reparamOnFace(face, dir);
   }
 }
+
+void gmshVertex::writeGEO(FILE *fp)
+{
+  fprintf(fp, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+          v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, v->lc);
+}
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
index afdd9b00c3..fcb3774c4f 100644
--- a/Geo/gmshVertex.h
+++ b/Geo/gmshVertex.h
@@ -33,6 +33,7 @@ class gmshVertex : public GVertex {
     v->lc = meshSize;
   }
   virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
+  virtual void writeGEO(FILE *fp);
 };
 
 #endif
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 3ebbf0f0f9..3dcdbda9df 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -1001,9 +1001,6 @@ of a given geometry point (@pxref{Points}). The last two cases permit to
 retrieve the indices of entities created through geometrical transformations
 and extrusions (see @ref{Transformations}, and @ref{Extrusions}).
 
-@c todo: explain in detail what numbers are returned once we decide what to
-@c do (chapeau et body pour extrude, etc.)...
-
 To see the practical use of such expressions, have a look at the first
 couple of examples in @ref{Tutorial}. Note that, in order to lighten the
 syntax, you can always omit the braces @code{@{@}} enclosing an
@@ -1971,6 +1968,26 @@ the X, Y and Z components of any point on this axis; the last
   Point | Line | Surface @{ @var{expression-list} @}; @dots{}
 @end example
 
+As explained in @ref{Floating point expressions}, @var{extrude} can be
+used in an expression, in which case it returns a list of identification
+numbers. By default, the list contains the ``top'' of the extruded
+entity at index 0 and the extruded entity at index 1, followed by the
+``sides'' of the extruded entity at indices 2, 3, etc.  For example:
+
+@example
+  Point(1) = @{0,0,0@};
+  Point(2) = @{1,0,0@};
+  Line(1) = @{1, 2@};
+  out[] = Extrude@{0,1,0@}@{ Line@{1@}; @};
+  Printf("top line = %g", out[0]);
+  Printf("surface = %g", out[1]);
+  Printf("side lines = %g and %g", out[2], out[3]);
+@end example
+
+This behaviour can be changed with the
+@code{Geometry.ExtrudeReturnLateralEntities} option (@pxref{Geometry
+options list}).
+
 @c .........................................................................
 @c Transformations
 @c .........................................................................
@@ -3155,7 +3172,7 @@ the nodes is given in @ref{Node ordering}.
 @item @var{number-of-string-tags}
 gives the number of string tags that follow. By default the first
 @var{string-tag} is interpreted as the name of the post-processing view.
-@c FIXME and the second as the name of the interpolation scheme.
+@c todo: and the second as the name of the interpolation scheme?
 
 @item @var{number-of-real-tags}
 gives the number of real number tags that follow. By default the first
-- 
GitLab