From 429dd4024f92b6ee1007baa756badd1cfaee6040 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-Fran=C3=A7ois=20Remacle=20=28students=29?=
 <jean-francois.remacle@uclouvain.be>
Date: Thu, 30 Apr 2009 15:25:33 +0000
Subject: [PATCH] added createTopology when Merging a MSH File and implemented
 discreteEdge

---
 Common/OpenFile.cpp           |   1 +
 Common/OpenFile.h             |   1 +
 Fltk/classificationEditor.cpp |  10 +-
 Geo/GEdge.h                   |   2 +
 Geo/GEdgeCompound.cpp         |  41 +--
 Geo/GModel.cpp                |   4 +-
 Geo/GModel.h                  |   3 +
 Geo/GModelIO_Geo.cpp          |  14 +-
 Geo/GModelIO_Mesh.cpp         | 115 +++++++-
 Geo/MVertex.h                 |   3 +-
 Geo/discreteEdge.cpp          | 187 ++++++++++++-
 Geo/discreteEdge.h            |  11 +-
 Geo/gmshEdge.cpp              |   2 +-
 Mesh/meshGEdge.cpp            |  28 +-
 Parser/Gmsh.tab.cpp           | 479 +++++++++++++++++-----------------
 Parser/Gmsh.y                 |   1 +
 16 files changed, 618 insertions(+), 284 deletions(-)

diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index bc5a3d4ebf..3f2adb4ccb 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -426,6 +426,7 @@ void OpenProject(std::string fileName)
   // merge the file
   MergeFile(fileName);
 
+
   // merge the associated option file if there is one
   if(!StatFile(fileName + ".opt"))
     MergeFile(fileName + ".opt");
diff --git a/Common/OpenFile.h b/Common/OpenFile.h
index ec2806c49e..c247f9f2d6 100644
--- a/Common/OpenFile.h
+++ b/Common/OpenFile.h
@@ -13,6 +13,7 @@ void ParseString(std::string str);
 void OpenProject(std::string filename);
 void OpenProjectMacFinder(const char *fileName);
 int MergeFile(std::string fileName, bool warnIfMissing=false);
+
 void ClearProject();
 void SetBoundingBox(double xmin, double xmax,
                     double ymin, double ymax, 
diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index 7dd6d1ad6a..c11971d98f 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -273,7 +273,7 @@ static GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
   std::map<std::pair<int, int>, GEdge*>::iterator it = 
     newEdges.find(std::make_pair<int, int>(i1, i2));
   if(it == newEdges.end()){
-    discreteEdge *temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1);
+    discreteEdge *temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
     GModel::current()->add(temporary);
     newEdges[std::make_pair<int, int>(i1, i2)] = temporary;
     return temporary;
@@ -383,7 +383,7 @@ static void updateedges_cb(Fl_Widget* w, void* data)
 {
   classificationEditor *e = (classificationEditor*)data;
  
-  // printf("%d edges detected\n", e->edges_detected.size());
+  printf("%d inside edges detected\n", e->edges_detected.size());
 
   for(unsigned int i = 0; i < e->temporary->lines.size(); i++){
     delete e->temporary->lines[i];
@@ -398,10 +398,12 @@ static void updateedges_cb(Fl_Widget* w, void* data)
     e->temporary->lines.push_back(new MLine(ea.v1, ea.v2));
   } 
 
+  printf("%d boundary edges detected\n",  e->edges_lonly.size());
   if(e->_togbuttons[CLASSTOGBUTTON_CLOS]->value()){
     for(unsigned int i = 0 ; i < e->edges_lonly.size(); i++){
       edge_angle ea = e->edges_lonly[i];
       e->temporary->lines.push_back(new MLine(ea.v1, ea.v2));
+      //check if closed loop
     } 
   }
   
@@ -620,9 +622,9 @@ classificationEditor::classificationEditor()
   // saved for the ones that have been saved by the user
   // and that will be used for next step
 
-  temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1);
+  temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
   GModel::current()->add(temporary);
-  saved = new discreteEdge(GModel::current(), maxEdgeNum() + 1);
+  saved = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
   GModel::current()->add(saved);
   
   _window->end();
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 02c3779ff2..cf1e2072ba 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -12,6 +12,7 @@
 #include "SVector3.h"
 #include "SPoint3.h"
 #include "SPoint2.h"
+#include "MLine.h"
 
 class MElement;
 class MLine;
@@ -36,6 +37,7 @@ class GEdge : public GEntity {
   GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1);
   virtual ~GEdge();
 
+
   // delete mesh data
   virtual void deleteMesh();
 
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 4c27a6fc55..a7130f6dda 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -16,7 +16,7 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
   v1 = _orientation[N-1] ? _compound[N-1]->getEndVertex() :   _compound[N-1]->getBeginVertex();
   v0->addEdge(this);
   v1->addEdge(this);
-  //  printf("%d -> %d\n",v0->tag(),v1->tag());
+  //printf("%d -> %d\n",v0->tag(),v1->tag());
   parametrize();
 }
 
@@ -30,19 +30,21 @@ void GEdgeCompound::orderEdges()
     edges.push_back(_compound[i]);
   }
 
-  // find a lonly edge
+ // find a lonly edge
   std::map<GVertex*,GEdge*> tempv;
   for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
     GVertex *v1 = (*it)->getBeginVertex();
-    GVertex *v2 = (*it)->getBeginVertex();
+    GVertex *v2 = (*it)->getEndVertex();
+    //printf("BEFORE ORDERING segment vert=%d, %d y =%g %g\n", v1->tag(), v2->tag() , v1->y(), v2->y());
     std::map<GVertex*,GEdge*>::iterator it1 = tempv.find(v1);
-    if (it1==tempv.end())tempv.insert(std::make_pair(v1,*it));
+    if (it1==tempv.end()) tempv.insert(std::make_pair(v1,*it));
     else tempv.erase(it1);
     std::map<GVertex*,GEdge*>::iterator it2 = tempv.find(v2);
-    if (it2==tempv.end())tempv.insert(std::make_pair(v2,*it));
+    if (it2==tempv.end()) tempv.insert(std::make_pair(v2,*it));
     else tempv.erase(it2);
   }
 
+ //find the first GEdge and erase it from the list edges
   GEdge *firstEdge;
   if (tempv.size()==2){   // non periodic
     firstEdge = (tempv.begin())->second;
@@ -60,8 +62,10 @@ void GEdgeCompound::orderEdges()
     }
   else{
     Msg::Error("EdgeCompound %d is wrong (it has %d end points)",tag(),tempv.size());
+    exit(1);
   }
 
+  //loop over all segments to order segments and store it in the list _c
   _c.push_back(firstEdge); 
   _orientation.push_back(1);
   GVertex *first = _c[0]->getBeginVertex();
@@ -72,9 +76,7 @@ void GEdgeCompound::orderEdges()
     bool found = false;
     for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
       GEdge *e = *it;
-      //      printf("last %d edge %d %d\n",last->tag(),e->getBeginVertex()->tag(),
-      //      	     e->getEndVertex()->tag());
-      if (e->getBeginVertex() == last){
+     if (e->getBeginVertex() == last){
 	_c.push_back(e); 
 	edges.erase(it);
 	_orientation.push_back(1);
@@ -97,16 +99,18 @@ void GEdgeCompound::orderEdges()
 	first = last;
 	last = temp;
 	_orientation[0] = 0;
-	printf("coucou\n");
-      }
+     }
       else {
 	Msg::Error("Compound Edge %d is wrong",tag());
 	return;
       }
     }
   }  
+
+  //edges is now a list of ordered GEdges
   _compound = _c;
 
+  //special case reverse orientation 
   if (_compound.size() < 2)return;
   if (_orientation[0] && _compound[0]->getEndVertex() != _compound[1]->getEndVertex() 
       && _compound[0]->getEndVertex() != _compound[1]->getBeginVertex()){  
@@ -115,14 +119,17 @@ void GEdgeCompound::orderEdges()
       _orientation[i] = !_orientation[i] ;
     }
   }
+
+//   for (int i=0;i<_compound.size();i++){
+//     printf("i = %d , o %d e %d (%d,%d)\n",
+// 	   i, (int)_orientation[i],
+// 	   _compound[i]->tag(),
+// 	   _compound[i]->getBeginVertex()->tag(),
+// 	   _compound[i]->getEndVertex()->tag());
+//  }
+//   exit(1);
+  
   return;
-   for (int i=0;i<_compound.size();i++){
-     printf("o %d e %d (%d,%d)\n",
-	    (int)_orientation[i],
-	    _compound[i]->tag(),
-	    _compound[i]->getBeginVertex()->tag(),
-	    _compound[i]->getEndVertex()->tag());
-   }
 
 }
 
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index c7039168d2..67051916df 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -659,9 +659,9 @@ void GModel::_storeElementsInEntities(std::map<int, std::vector<MElement*> > &ma
       {
         GEdge *e = getEdgeByTag(it->first);
         if(!e){
-          e = new discreteEdge(this, it->first);
+          e = new discreteEdge(this, it->first, 0, 0);
           add(e);
-        }
+       }
         _addElements(e->lines, it->second);
       }
       break;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 400f2ffcd5..a5c4fb2cff 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -297,6 +297,9 @@ class GModel
   int writeMSH(const std::string &name, double version=1.0, bool binary=false,
                bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0);
 
+  //Create topology from mesh
+  void createTopologyFromMSH();
+
   // Mesh statistics (as Gmsh post-processing views)
   int writePOS(const std::string &name, bool printElementary,
                bool printElementNumber, bool printGamma, bool printEta, bool printRho,
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 072ea30450..fb4aa77443 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -42,6 +42,8 @@ int GModel::readGEO(const std::string &name)
 
 int GModel::importGEOInternals()
 {
+
+  //printf("Dans import GEO internals \n");
   if(Tree_Nbr(_geo_internals->Points)) {
     List_T *points = Tree2List(_geo_internals->Points);
     for(int i = 0; i < List_Nbr(points); i++){
@@ -67,7 +69,7 @@ int GModel::importGEOInternals()
           e = new gmshEdge(this, c,
                            getVertexByTag(c->beg->Num),
                            getVertexByTag(c->end->Num));
-          add(e);
+         add(e);
         }
         else
           e->resetMeshAttributes();
@@ -181,6 +183,16 @@ int GModel::importGEOInternals()
   Msg::Debug("%d Faces", faces.size());
   Msg::Debug("%d Regions", regions.size());
 
+  for ( std::set<GVertex*, MVertexLessThanLexicographic>::iterator it  = vertices.begin(); it != vertices.end(); it++){
+    printf("WARNING:import GEO vert of Type: %s \n", (*it)->getTypeString().c_str());
+   }
+  for ( std::set<GEdge*, MVertexLessThanLexicographic>::iterator it  = edges.begin(); it != edges.end(); it++){
+    printf("WARNING:import GEO edge of Type: %s \n", (*it)->getTypeString().c_str());
+  }
+  for ( std::set<GFace*, MVertexLessThanLexicographic>::iterator it  = faces.begin(); it != faces.end(); it++){
+    printf("WARNING:import GEO face of Type:  %s \n", (*it)->getTypeString().c_str());
+  }
+
   return 1;
 }
 
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 83e2a59435..6728d8fb54 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -22,6 +22,14 @@
 #include "discreteFace.h"
 #include "StringUtils.h"
 #include "GmshMessage.h"
+#include "discreteVertex.h"
+#include "discreteEdge.h"
+#include "discreteFace.h"
+#include "discreteRegion.h"
+#include "MElement.h"
+#include "GEdgeCompound.h"
+
+#include <iostream> // DBG
 
 static void storePhysicalTagsInEntities(GModel *m, int dim,
                                         std::map<int, std::map<int, std::string> > &map)
@@ -106,6 +114,107 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
   if(part) m->getMeshPartitions().insert(part);
 }
 
+void GModel::createTopologyFromMSH(){
+
+  //printf("Dans createTopologyFromMSH \n");
+
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
+  std::vector<discreteVertex*> vertices;
+  std::vector<discreteEdge*> edges;
+  std::vector<discreteFace*> faces;
+  std::vector<discreteRegion*> regions;
+
+  for (std::vector<GEntity*>::iterator entity = entities.begin(); entity != entities.end(); entity++) {
+    switch ((*entity)->dim()) {
+    case 0:
+      vertices.push_back((discreteVertex*) *entity);
+      break;
+    case 1:
+      edges.push_back((discreteEdge*) *entity);
+      break;
+    case 2:
+      faces.push_back((discreteFace*) *entity);
+      break;
+    case 3:
+      regions.push_back((discreteRegion*) *entity);
+      break;
+    }
+  }
+  //printf("vertices size =%d \n", vertices.size());
+  //printf("edges size =%d \n", edges.size());
+  //printf("faces size =%d \n", faces.size());
+  //printf("regions size =%d \n", regions.size());
+
+  int tag = 100;
+  for (std::vector<discreteEdge*>::iterator edge = edges.begin(); edge != edges.end(); edge++){
+    if (tag < (*edge)->tag() ) tag = (*edge)->tag() + 1;
+  }
+
+ //For each discreteEdge, build a new GEdgeCompound
+  for (std::vector<discreteEdge*>::iterator edge = edges.begin(); edge != edges.end(); edge++){
+
+    //printf("createTopology: %d  EDGES, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
+
+    //create a map with the tags of the mesh vertices
+    std::map<int, GVertex*> myMap;
+    for (std::vector<MLine*>::const_iterator it = (*edge)->lines.begin() ; it != (*edge)->lines.end() ; ++it){  
+      int tagB = (*it)->getVertex(0)->getNum();
+      int tagE = (*it)->getVertex(1)->getNum();
+
+      std::map<int, GVertex*>::iterator it1 = myMap.find(tagB);
+      std::map<int, GVertex*>::iterator it2 = myMap.find(tagE);
+      if (it1 == myMap.end()){
+	GVertex *gvB = new discreteVertex(this,tagB);
+	gvB->mesh_vertices.push_back((*it)->getVertex(0)); 
+	gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
+	myMap.insert(std::make_pair(tagB, gvB));
+      }
+      if (it2 == myMap.end()){
+	GVertex *gvE = new discreteVertex(this,tagE);
+	gvE->mesh_vertices.push_back((*it)->getVertex(1)); 
+	gvE->points.push_back(new MPoint(gvE->mesh_vertices.back()));
+	myMap.insert(std::make_pair(tagE, gvE));
+      }
+    }
+
+ //    for(std::map<int, GVertex*>::const_iterator it = myMap.begin(); it != myMap.end(); ++it){
+//       printf(" tag=%d tagsize=%d\n", it->first, myMap.size());
+//     }
+
+    //create a vector composed of plenty of discreteEdges from the Mlines of the original discreteVertex
+    std::vector<GEdge*> e_compound;
+
+   for (std::vector<MLine*>::const_iterator it = (*edge)->lines.begin() ; it != (*edge)->lines.end() ; ++it){  
+     //printf("MLine =%d %d \n", (*it)->getVertex(0)->getNum(), (*it)->getVertex(1)->getNum());
+
+      int tagB = (*it)->getVertex(0)->getNum();
+      int tagE = (*it)->getVertex(1)->getNum();
+      std::map<int,GVertex*>::iterator it1 = myMap.find(tagB);
+      std::map<int,GVertex*>::iterator it2 = myMap.find(tagE);
+      GVertex *gvB = it1->second;
+      GVertex *gvE = it2->second;
+
+      GEdge *temp = new discreteEdge(this, tag, gvB, gvE); //new GEdge corresponding to the MLine
+      gvB->addEdge(temp);
+      gvE->addEdge(temp);
+
+      e_compound.push_back(temp); //add the compound to the GEdge
+      tag ++;
+
+    }
+
+   //now, we can create the GEdgeCompound
+    GEdge *gec = new GEdgeCompound(this, tag, e_compound);
+    add(gec);
+
+  }
+
+  return;
+
+}
+
 int GModel::readMSH(const std::string &name)
 {
   FILE *fp = fopen(name.c_str(), "rb");
@@ -569,7 +678,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   writeElementHeaderMSH
     (binary, fp, elements, MSH_LIN_2, MSH_LIN_3, MSH_LIN_4, MSH_LIN_5);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
+   writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
   writeElementHeaderMSH(binary, fp, elements, MSH_TRI_3, MSH_TRI_6, MSH_TRI_9, 
                         MSH_TRI_10, MSH_TRI_12, MSH_TRI_15, MSH_TRI_15I, MSH_TRI_21);
@@ -950,7 +1059,7 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
     return 0;
   }
   Msg::Info("%d elements", elements[0][region].size() + 
-      elements[1][region].size() + elements[2][region].size());
+	    elements[1][region].size() + elements[2][region].size());
   return 1;
 }
 
@@ -2279,7 +2388,7 @@ int GModel::readDIFF(const std::string &name)
     while(strstr(str, "Max number of nodes in an element:")==NULL){
       if(!fgets(str, sizeof(str), fp) || feof(fp))
         break;
-      }
+    }
     if(sscanf(str, "%*s %*s %*s %*s %*s %*s %*s %d", &numVerticesPerElement) != 1)
       return 0;
     Msg::Info("numVerticesPerElement %d",numVerticesPerElement); 
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 3873f65405..e73821e730 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -70,8 +70,7 @@ class MVertex{
   inline void setPolynomialOrder(char order){ _order = order; }
 
   // get/set the coordinates
-  inline double x() const { return _x; }
-  inline double y() const { return _y; }
+  inline double x() const { return _x; }  inline double y() const { return _y; }
   inline double z() const { return _z; }
   inline double & x() { return _x; }
   inline double & y() { return _y; }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index cb966f103d..1cb0ee5c64 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -6,33 +6,204 @@
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "discreteEdge.h"
+#include "MLine.h"
+#include "Numeric.h"
+
+#include <vector>
+#include <list>
 
 #if !defined(HAVE_GMSH_EMBEDDED)
 #include "Geo.h"
 #endif
 
-discreteEdge::discreteEdge(GModel *model, int num) : GEdge(model, num, 0, 0) 
+discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
+  : GEdge(model, num, _v0, _v1)
 {
 #if !defined(HAVE_GMSH_EMBEDDED)
   Curve *c = Create_Curve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1.);
   Tree_Add(model->getGEOInternals()->Curves, &c);
   CreateReversedCurve(c);
 #endif
+  
 }
 
-GPoint discreteEdge::point(double p) const 
+void discreteEdge::orderEdge() 
 {
-  Msg::Error("Cannot evaluate point on discrete edge");
-  return GPoint();
+
+  printf(" *** ORDERING DISCRETE EDGE %d of size %d \n", this->tag(), lines.size());
+
+  std::vector<MLine*> _m ;  
+  std::list<MLine*> segments;
+
+  //store all lines in a list : segments
+  for (int i=0;i<lines.size();i++){
+    //printf("BEFORE ORDERING segment i=%d, vert=%d, %d\n", i, lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
+    segments.push_back(lines[i]);
+  }
+
+  // find a lonly MLine
+  std::map<MVertex*,MLine*> boundv;
+  for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+    MVertex *vL = (*it)->getVertex(0);
+    MVertex *vR = (*it)->getVertex(1);
+    //printf("MLIne %d %d = (%g, %g, %g) (%g, %g, %g)\n",vL->getNum(), vR->getNum(), vL->x(),vL->y(),vL->z(), vR->x(),vR->y(),vR->z() );
+    std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
+    if (it1==boundv.end()) boundv.insert(std::make_pair(vL,*it));
+    else boundv.erase(it1);
+    std::map<MVertex*,MLine*>::iterator it2 = boundv.find(vR);
+    if (it2==boundv.end())       boundv.insert(std::make_pair(vR,*it));
+    else boundv.erase(it2);
+  }   
+
+  //find the first MLine and erase it from the list segments
+  MLine *firstLine;
+  if (boundv.size()==2){   // non periodic
+    firstLine = (boundv.begin())->second;
+    for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+      if (*it == firstLine){
+	segments.erase(it);
+	break;
+      }
+    }    
+  }
+  else if (boundv.size()==0) // periodic
+    {
+      firstLine = *(segments.begin());
+      segments.erase(segments.begin());
+    }
+  else{
+    Msg::Error("EdgeCompound %d is wrong (it has %d end points)",tag(),boundv.size());
+  }
+
+  //loop over all segments to order segments and store it in the list _m
+  _m.push_back(firstLine); 
+  _orientation.push_back(1);
+  MVertex *first = _m[0]->getVertex(0);
+  MVertex *last = _m[0]->getVertex(1);   
+  //printf("first =%d last =%d \n", first->getNum(), last->getNum());
+  while (first != last){
+    if (segments.empty())break;
+    bool found = false;
+    for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+      MLine *e = *it;
+      if (e->getVertex(0) == last){
+	printf("orientation 1: beginV=%d \n", e->getVertex(0)->getNum());
+	_m.push_back(e); 
+	segments.erase(it);
+	_orientation.push_back(1);
+	last = e->getVertex(1);
+	found = true;
+	break;
+      }
+      else if (e->getVertex(1) == last){
+	printf("orientation 0: endV=%d \n", e->getVertex(1)->getNum());
+	_m.push_back(e); 
+	segments.erase(it);
+	_orientation.push_back(0);
+	last = e->getVertex(0);
+	found = true;
+	break;
+      }
+    }
+    if (!found  && _orientation[0]){ //reverse orientation of first Line
+      if (_m.size() == 1 ){
+	MVertex *temp = first;
+	first = last;
+	last = temp;
+	_orientation[0] = 0;
+     }
+      else {
+	Msg::Error("Discrete Edge %d is wrong",tag());
+	return;
+      }
+    }
+  }  
+
+  //lines is now a list of ordered MLines
+  lines = _m;
+
+  //special case reverse orientation 
+  if (lines.size() < 2) return;
+  if (_orientation[0] && lines[0]->getVertex(1) != lines[1]->getVertex(1) 
+      && lines[0]->getVertex(1) != lines[1]->getVertex(0)){  
+    printf("coucou here \n");
+    for (int i=0;i<lines.size();i++) _orientation[i] = !_orientation[i] ;
+  }
+ 
+  for (int i=0;i<lines.size();i++){
+    printf("AFTER ORDERING segment or=%d, vert=%d, %d\n", (int)_orientation[i], lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
+  }
+
+ return;
+
+}
+void discreteEdge::parametrize() 
+{
+
+  for (int i=0;i<lines.size()+1;i++){
+    _pars.push_back(i);
+  }   
+
+/*
+  +---------------+--------------+----------- ... ----------+
+  _pars[0]=0   _pars[1]=1    _pars[2]=2             _pars[N+1]=N+1
+*/
+  
+//    for (int i=0;i<lines.size()+1;i++){
+//      printf("_pars[%d]=%g\n",i, _pars[i] );
+//    }
+
+}
+
+void discreteEdge::getLocalParameter ( const double &t,
+					int &iLine,
+					double & tLoc) const
+{
+
+  for (iLine=0 ; iLine<lines.size() ;iLine++){
+   double tmin = _pars[iLine];
+    double tmax = _pars[iLine+1];
+    if (t >= tmin && t <= tmax){      
+      tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin)  :  1 - (t-tmin)/(tmax-tmin)  ;
+      //printf("getlocal param t=%g, iLine=%d, tLoc=%g \n", t, iLine, tLoc);
+     return;
+    }
+  }
+}
+
+GPoint discreteEdge::point(double par) const 
+{
+  
+  //printf("dans point v0 =%d (%g,%g,%g), v1=%d (%g,%g,%g)\n", v0->tag(), v0->x(),v0->y(), v0->z(), v1->tag(), v1->x(),v1->y(), v1->z());
+
+  double x, y, z;
+  
+  x = v0->x() + par * (v1->x()- v0->x());
+  y = v0->y() + par * (v1->y()- v0->y());
+  z = v0->z() + par * (v1->z()- v0->z());
+
+  //printf("Discrete Edge POint par=%g, x= %g, y = %g, z = %g \n", par, x,y,z);
+  return GPoint(x,y,z);
 }
 
 SVector3 discreteEdge::firstDer(double par) const 
 {
-  Msg::Error("Cannot evaluate derivative on discrete edge");
-  return SVector3();
+
+  double dx,dy,dz;
+  double dt = sqrt((v1->x()-v0->x())*(v1->x()-v0->x()) + (v1->y()-v0->y())*(v1->y()-v0->y()) + (v1->z()-v0->z())*(v1->z()-v0->z()));
+  dx = (v1->x() - v0->x() ) /dt;
+  dy = (v1->y() - v0->y() ) /dt;
+  dz = (v1->z() - v0->z() ) /dt;
+
+  //printf("firstDer discreteEdge  par=%g, dx=%g, dy=%g dz=%g dt=%g \n", par,dx,dy,dz, dt);
+  return SVector3(dx,dy,dz);
+
 }
 
 Range<double> discreteEdge::parBounds(int i) const {
-  Msg::Error("Cannot specify bounds for parametric coordinate on discrete edge");
-  return Range<double>(0,0);
+
+  return Range<double>(0, 1);
+  //return Range<double>(0, lines.size()-1); 
+  //return Range<double>(0, _pars[lines.size()-1]); 
+  
 }
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index fc24c51275..c69e9b0564 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -10,13 +10,22 @@
 #include "GEdge.h"
 
 class discreteEdge : public GEdge {
+ protected:
+  std::vector<double> _pars;
+  std::vector<int> _orientation;
  public:
-  discreteEdge(GModel *model, int num);
+  discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
   virtual ~discreteEdge() {}
+  void getLocalParameter ( const double &t,
+			   int &iEdge,
+			   double & tLoc) const;
   virtual GeomType geomType() const { return DiscreteCurve; }
   virtual GPoint point(double p) const;
   virtual SVector3 firstDer(double par) const;
   virtual Range<double> parBounds(int) const;
+  void setVertices(GVertex *_v0, GVertex *_v1){ v0 = _v0; v1 = _v1; }
+  void orderEdge();
+  void parametrize();
 };
 
 #endif
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 6ce4647ce9..e6bd953fa4 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -30,7 +30,7 @@ void gmshEdge::resetMeshAttributes()
 
 Range<double> gmshEdge::parBounds(int i) const
 { 
-  return Range<double>(c->ubeg, c->uend);
+ return Range<double>(c->ubeg, c->uend);
 }
 
 GPoint gmshEdge::point(double par) const
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index ce1ec9bd7f..bba8021fd6 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -237,7 +237,9 @@ void meshGEdge::operator() (GEdge *ge)
 
   ge->model()->setCurrentMeshEntity(ge);
 
-  if(ge->geomType() == GEntity::DiscreteCurve) return;
+  if(ge->geomType() == GEntity::DiscreteCurve) {
+    return;
+  }
   if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
   if(ge->meshAttributes.Method == MESH_NONE) return;
   if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;
@@ -276,6 +278,7 @@ void meshGEdge::operator() (GEdge *ge)
   else if(ge->meshAttributes.Method == MESH_TRANSFINITE){
     a = Integration(ge, t_begin, t_end, F_Transfinite, Points, 1.e-8);
     N = ge->meshAttributes.nbPointsTransfinite;
+    printf("Mesh transfinite N=%d, a=%g \n", N, a);
   }
   else{
     if(CTX::instance()->mesh.lcIntegrationPrecision > 1.e-8){
@@ -290,8 +293,9 @@ void meshGEdge::operator() (GEdge *ge)
 		      CTX::instance()->mesh.lcIntegrationPrecision);
     }
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
+    //printf("Mesh NOT transfinite N=%d, a=%g \n", N, a);
   }
-  
+
   // if the curve is periodic and if the begin vertex is identical to
   // the end vertex and if this vertex has only one model curve
   // adjacent to it, then the vertex is not connecting any other
@@ -301,6 +305,7 @@ void meshGEdge::operator() (GEdge *ge)
   if(ge->getBeginVertex() == ge->getEndVertex() && 
      ge->getBeginVertex()->edges().size() == 1){
     end_p = beg_p = ge->point(t_begin);
+    printf("Meshing periodic closed curve \n");
   }
   else{
     MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
@@ -341,11 +346,10 @@ void meshGEdge::operator() (GEdge *ge)
   }
 
   for(unsigned int i = 0; i < ge->mesh_vertices.size() + 1; i++){
-    MVertex *v0 = (i == 0) ? 
-      ge->getBeginVertex()->mesh_vertices[0] : ge->mesh_vertices[i - 1];
-    MVertex *v1 = (i == ge->mesh_vertices.size()) ? 
-      ge->getEndVertex()->mesh_vertices[0] : ge->mesh_vertices[i];
+    MVertex *v0 = (i == 0) ?       ge->getBeginVertex()->mesh_vertices[0] : ge->mesh_vertices[i - 1];
+    MVertex *v1 = (i == ge->mesh_vertices.size()) ?       ge->getEndVertex()->mesh_vertices[0] : ge->mesh_vertices[i];
     ge->lines.push_back(new MLine(v0, v1));
+    //printf("New Line v0=%g v1=%g \n", v0->y(), v1->y());
   }
 
   if(ge->getBeginVertex() == ge->getEndVertex() && 
@@ -355,4 +359,16 @@ void meshGEdge::operator() (GEdge *ge)
     v0->y() = beg_p.y();
     v0->z() = beg_p.z();
   }
+
+  //for (std::vector<MLine*>::iterator it = ge->lines.begin() ; it != ge->lines.end() ; ++it){
+  //  printf("line with : first = %d last=%d \n",  (*it)->getVertex(0)->getIndex(), (*it)->getVertex(1)->getIndex() );
+  //}
+
+//   //if DiscreteCurve, erase all old MLines
+//   if(ge->geomType() == GEntity::DiscreteCurve) {
+//     ge->lines.erase(ge->lines.begin(), ge->lines.end()-(N-1));
+//   }
+
+ 
+
 }
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 8296588f78..4063e836eb 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -981,28 +981,28 @@ static const yytype_uint16 yyrline[] =
     1380,  1405,  1430,  1446,  1462,  1490,  1510,  1528,  1545,  1566,
     1571,  1576,  1581,  1586,  1606,  1612,  1623,  1624,  1629,  1632,
     1636,  1659,  1682,  1705,  1733,  1742,  1746,  1761,  1788,  1805,
-    1819,  1825,  1831,  1840,  1854,  1902,  1920,  1935,  1954,  1966,
-    1990,  1994,  1999,  2004,  2015,  2032,  2049,  2068,  2087,  2115,
-    2123,  2129,  2136,  2140,  2149,  2157,  2165,  2174,  2173,  2186,
-    2185,  2198,  2197,  2210,  2209,  2222,  2229,  2236,  2243,  2250,
-    2257,  2264,  2271,  2278,  2286,  2285,  2297,  2296,  2308,  2307,
-    2319,  2318,  2330,  2329,  2341,  2340,  2352,  2351,  2363,  2362,
-    2374,  2373,  2388,  2391,  2397,  2406,  2426,  2449,  2453,  2477,
-    2480,  2496,  2499,  2512,  2515,  2521,  2524,  2531,  2587,  2657,
-    2662,  2729,  2772,  2798,  2821,  2844,  2847,  2856,  2860,  2876,
-    2877,  2878,  2879,  2880,  2881,  2882,  2883,  2884,  2891,  2892,
-    2893,  2894,  2895,  2896,  2897,  2898,  2899,  2900,  2901,  2902,
-    2903,  2904,  2905,  2906,  2907,  2908,  2909,  2910,  2911,  2912,
-    2913,  2914,  2915,  2916,  2917,  2918,  2919,  2920,  2921,  2922,
-    2924,  2925,  2926,  2927,  2928,  2929,  2930,  2931,  2932,  2933,
-    2934,  2935,  2936,  2937,  2938,  2939,  2940,  2941,  2942,  2943,
-    2944,  2953,  2954,  2955,  2956,  2957,  2958,  2959,  2963,  2976,
-    2988,  3003,  3013,  3023,  3041,  3046,  3051,  3061,  3071,  3079,
-    3083,  3087,  3091,  3095,  3102,  3106,  3110,  3114,  3121,  3126,
-    3133,  3138,  3142,  3147,  3151,  3159,  3170,  3174,  3186,  3194,
-    3202,  3209,  3220,  3240,  3250,  3260,  3270,  3290,  3295,  3299,
-    3303,  3315,  3319,  3331,  3338,  3348,  3352,  3367,  3372,  3379,
-    3383,  3396,  3404,  3415,  3419,  3427,  3435,  3449,  3463,  3467
+    1819,  1825,  1831,  1840,  1854,  1903,  1921,  1936,  1955,  1967,
+    1991,  1995,  2000,  2005,  2016,  2033,  2050,  2069,  2088,  2116,
+    2124,  2130,  2137,  2141,  2150,  2158,  2166,  2175,  2174,  2187,
+    2186,  2199,  2198,  2211,  2210,  2223,  2230,  2237,  2244,  2251,
+    2258,  2265,  2272,  2279,  2287,  2286,  2298,  2297,  2309,  2308,
+    2320,  2319,  2331,  2330,  2342,  2341,  2353,  2352,  2364,  2363,
+    2375,  2374,  2389,  2392,  2398,  2407,  2427,  2450,  2454,  2478,
+    2481,  2497,  2500,  2513,  2516,  2522,  2525,  2532,  2588,  2658,
+    2663,  2730,  2773,  2799,  2822,  2845,  2848,  2857,  2861,  2877,
+    2878,  2879,  2880,  2881,  2882,  2883,  2884,  2885,  2892,  2893,
+    2894,  2895,  2896,  2897,  2898,  2899,  2900,  2901,  2902,  2903,
+    2904,  2905,  2906,  2907,  2908,  2909,  2910,  2911,  2912,  2913,
+    2914,  2915,  2916,  2917,  2918,  2919,  2920,  2921,  2922,  2923,
+    2925,  2926,  2927,  2928,  2929,  2930,  2931,  2932,  2933,  2934,
+    2935,  2936,  2937,  2938,  2939,  2940,  2941,  2942,  2943,  2944,
+    2945,  2954,  2955,  2956,  2957,  2958,  2959,  2960,  2964,  2977,
+    2989,  3004,  3014,  3024,  3042,  3047,  3052,  3062,  3072,  3080,
+    3084,  3088,  3092,  3096,  3103,  3107,  3111,  3115,  3122,  3127,
+    3134,  3139,  3143,  3148,  3152,  3160,  3171,  3175,  3187,  3195,
+    3203,  3210,  3221,  3241,  3251,  3261,  3271,  3291,  3296,  3300,
+    3304,  3316,  3320,  3332,  3339,  3349,  3353,  3368,  3373,  3380,
+    3384,  3397,  3405,  3416,  3420,  3428,  3436,  3450,  3464,  3468
 };
 #endif
 
@@ -5723,6 +5723,7 @@ yyreduce:
 	char tmpstring[1024];
 	FixRelativePath((yyvsp[(2) - (3)].c), tmpstring);
 	MergeFile(tmpstring, true);
+	GModel::current()->createTopologyFromMSH();
       }
       else if(!strcmp((yyvsp[(1) - (3)].c), "System"))
 	SystemCall((yyvsp[(2) - (3)].c));
@@ -5733,7 +5734,7 @@ yyreduce:
     break;
 
   case 145:
-#line 1903 "Gmsh.y"
+#line 1904 "Gmsh.y"
     {
 #if !defined(HAVE_NO_POST)
       if(!strcmp((yyvsp[(1) - (7)].c), "Save") && !strcmp((yyvsp[(2) - (7)].c), "View")){
@@ -5754,7 +5755,7 @@ yyreduce:
     break;
 
   case 146:
-#line 1921 "Gmsh.y"
+#line 1922 "Gmsh.y"
     {
 #if !defined(HAVE_NO_POST)
       if(!strcmp((yyvsp[(1) - (7)].c), "Background") && !strcmp((yyvsp[(2) - (7)].c), "Mesh")  && !strcmp((yyvsp[(3) - (7)].c), "View")){
@@ -5772,7 +5773,7 @@ yyreduce:
     break;
 
   case 147:
-#line 1936 "Gmsh.y"
+#line 1937 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(1) - (3)].c), "Sleep")){
 	SleepInSeconds((yyvsp[(2) - (3)].d));
@@ -5794,7 +5795,7 @@ yyreduce:
     break;
 
   case 148:
-#line 1955 "Gmsh.y"
+#line 1956 "Gmsh.y"
     {
 #if !defined(HAVE_NO_POST)
        try {
@@ -5809,7 +5810,7 @@ yyreduce:
     break;
 
   case 149:
-#line 1967 "Gmsh.y"
+#line 1968 "Gmsh.y"
     {
 #if !defined(HAVE_NO_POST)
       if(!strcmp((yyvsp[(2) - (3)].c), "ElementsFromAllViews"))
@@ -5836,14 +5837,14 @@ yyreduce:
     break;
 
   case 150:
-#line 1991 "Gmsh.y"
+#line 1992 "Gmsh.y"
     {
       exit(0);
     ;}
     break;
 
   case 151:
-#line 1995 "Gmsh.y"
+#line 1996 "Gmsh.y"
     {
       CTX::instance()->forcedBBox = 0;
       SetBoundingBox();
@@ -5851,7 +5852,7 @@ yyreduce:
     break;
 
   case 152:
-#line 2000 "Gmsh.y"
+#line 2001 "Gmsh.y"
     {
       CTX::instance()->forcedBBox = 1;
       SetBoundingBox((yyvsp[(3) - (15)].d), (yyvsp[(5) - (15)].d), (yyvsp[(7) - (15)].d), (yyvsp[(9) - (15)].d), (yyvsp[(11) - (15)].d), (yyvsp[(13) - (15)].d));
@@ -5859,7 +5860,7 @@ yyreduce:
     break;
 
   case 153:
-#line 2005 "Gmsh.y"
+#line 2006 "Gmsh.y"
     {
 #if defined(HAVE_FLTK)
       Draw();
@@ -5868,7 +5869,7 @@ yyreduce:
     break;
 
   case 154:
-#line 2016 "Gmsh.y"
+#line 2017 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(3) - (6)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(5) - (6)].d);
@@ -5888,7 +5889,7 @@ yyreduce:
     break;
 
   case 155:
-#line 2033 "Gmsh.y"
+#line 2034 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(3) - (8)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(5) - (8)].d);
@@ -5908,7 +5909,7 @@ yyreduce:
     break;
 
   case 156:
-#line 2050 "Gmsh.y"
+#line 2051 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(5) - (8)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(7) - (8)].d);
@@ -5930,7 +5931,7 @@ yyreduce:
     break;
 
   case 157:
-#line 2069 "Gmsh.y"
+#line 2070 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(5) - (10)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(7) - (10)].d);
@@ -5952,7 +5953,7 @@ yyreduce:
     break;
 
   case 158:
-#line 2088 "Gmsh.y"
+#line 2089 "Gmsh.y"
     {
       if(ImbricatedLoop <= 0){
 	yymsg(0, "Invalid For/EndFor loop");
@@ -5983,7 +5984,7 @@ yyreduce:
     break;
 
   case 159:
-#line 2116 "Gmsh.y"
+#line 2117 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->createFunction
          ((yyvsp[(2) - (2)].c), gmsh_yyin, gmsh_yyname, gmsh_yylineno))
@@ -5994,7 +5995,7 @@ yyreduce:
     break;
 
   case 160:
-#line 2124 "Gmsh.y"
+#line 2125 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->leaveFunction
          (&gmsh_yyin, gmsh_yyname, gmsh_yylineno))
@@ -6003,7 +6004,7 @@ yyreduce:
     break;
 
   case 161:
-#line 2130 "Gmsh.y"
+#line 2131 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->enterFunction
          ((yyvsp[(2) - (3)].c), &gmsh_yyin, gmsh_yyname, gmsh_yylineno))
@@ -6013,20 +6014,20 @@ yyreduce:
     break;
 
   case 162:
-#line 2137 "Gmsh.y"
+#line 2138 "Gmsh.y"
     {
       if(!(yyvsp[(3) - (4)].d)) skip_until("If", "EndIf");
     ;}
     break;
 
   case 163:
-#line 2141 "Gmsh.y"
+#line 2142 "Gmsh.y"
     {
     ;}
     break;
 
   case 164:
-#line 2150 "Gmsh.y"
+#line 2151 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, (yyvsp[(4) - (5)].l), 
@@ -6037,7 +6038,7 @@ yyreduce:
     break;
 
   case 165:
-#line 2158 "Gmsh.y"
+#line 2159 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, (yyvsp[(10) - (11)].l), 
@@ -6048,7 +6049,7 @@ yyreduce:
     break;
 
   case 166:
-#line 2166 "Gmsh.y"
+#line 2167 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[(12) - (13)].l), 
@@ -6059,14 +6060,14 @@ yyreduce:
     break;
 
   case 167:
-#line 2174 "Gmsh.y"
+#line 2175 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 168:
-#line 2178 "Gmsh.y"
+#line 2179 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, (yyvsp[(4) - (7)].l), 
@@ -6077,14 +6078,14 @@ yyreduce:
     break;
 
   case 169:
-#line 2186 "Gmsh.y"
+#line 2187 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 170:
-#line 2190 "Gmsh.y"
+#line 2191 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, (yyvsp[(10) - (13)].l), 
@@ -6095,14 +6096,14 @@ yyreduce:
     break;
 
   case 171:
-#line 2198 "Gmsh.y"
+#line 2199 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 172:
-#line 2202 "Gmsh.y"
+#line 2203 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[(12) - (15)].l), 
@@ -6113,14 +6114,14 @@ yyreduce:
     break;
 
   case 173:
-#line 2210 "Gmsh.y"
+#line 2211 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 174:
-#line 2214 "Gmsh.y"
+#line 2215 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(BOUNDARY_LAYER, (yyvsp[(3) - (6)].l), 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -6130,7 +6131,7 @@ yyreduce:
     break;
 
   case 175:
-#line 2223 "Gmsh.y"
+#line 2224 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[(4) - (8)].d), 
@@ -6140,7 +6141,7 @@ yyreduce:
     break;
 
   case 176:
-#line 2230 "Gmsh.y"
+#line 2231 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (8)].d), 
@@ -6150,7 +6151,7 @@ yyreduce:
     break;
 
   case 177:
-#line 2237 "Gmsh.y"
+#line 2238 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (8)].d), 
@@ -6160,7 +6161,7 @@ yyreduce:
     break;
 
   case 178:
-#line 2244 "Gmsh.y"
+#line 2245 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[(4) - (12)].d), 
@@ -6170,7 +6171,7 @@ yyreduce:
     break;
 
   case 179:
-#line 2251 "Gmsh.y"
+#line 2252 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (12)].d), 
@@ -6180,7 +6181,7 @@ yyreduce:
     break;
 
   case 180:
-#line 2258 "Gmsh.y"
+#line 2259 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (12)].d), 
@@ -6190,7 +6191,7 @@ yyreduce:
     break;
 
   case 181:
-#line 2265 "Gmsh.y"
+#line 2266 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[(4) - (14)].d), 
@@ -6200,7 +6201,7 @@ yyreduce:
     break;
 
   case 182:
-#line 2272 "Gmsh.y"
+#line 2273 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (14)].d), 
@@ -6210,7 +6211,7 @@ yyreduce:
     break;
 
   case 183:
-#line 2279 "Gmsh.y"
+#line 2280 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (14)].d), 
@@ -6220,14 +6221,14 @@ yyreduce:
     break;
 
   case 184:
-#line 2286 "Gmsh.y"
+#line 2287 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 185:
-#line 2290 "Gmsh.y"
+#line 2291 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[(4) - (12)].d), 
@@ -6237,14 +6238,14 @@ yyreduce:
     break;
 
   case 186:
-#line 2297 "Gmsh.y"
+#line 2298 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 187:
-#line 2301 "Gmsh.y"
+#line 2302 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (12)].d), 
@@ -6254,14 +6255,14 @@ yyreduce:
     break;
 
   case 188:
-#line 2308 "Gmsh.y"
+#line 2309 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 189:
-#line 2312 "Gmsh.y"
+#line 2313 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (12)].d), 
@@ -6271,14 +6272,14 @@ yyreduce:
     break;
 
   case 190:
-#line 2319 "Gmsh.y"
+#line 2320 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 191:
-#line 2323 "Gmsh.y"
+#line 2324 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[(4) - (16)].d), 
@@ -6288,14 +6289,14 @@ yyreduce:
     break;
 
   case 192:
-#line 2330 "Gmsh.y"
+#line 2331 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 193:
-#line 2334 "Gmsh.y"
+#line 2335 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (16)].d), 
@@ -6305,14 +6306,14 @@ yyreduce:
     break;
 
   case 194:
-#line 2341 "Gmsh.y"
+#line 2342 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 195:
-#line 2345 "Gmsh.y"
+#line 2346 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (16)].d), 
@@ -6322,14 +6323,14 @@ yyreduce:
     break;
 
   case 196:
-#line 2352 "Gmsh.y"
+#line 2353 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 197:
-#line 2356 "Gmsh.y"
+#line 2357 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[(4) - (18)].d), 
@@ -6339,14 +6340,14 @@ yyreduce:
     break;
 
   case 198:
-#line 2363 "Gmsh.y"
+#line 2364 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 199:
-#line 2367 "Gmsh.y"
+#line 2368 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (18)].d), 
@@ -6356,14 +6357,14 @@ yyreduce:
     break;
 
   case 200:
-#line 2374 "Gmsh.y"
+#line 2375 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 201:
-#line 2378 "Gmsh.y"
+#line 2379 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (18)].d), 
@@ -6373,19 +6374,19 @@ yyreduce:
     break;
 
   case 202:
-#line 2389 "Gmsh.y"
+#line 2390 "Gmsh.y"
     {
     ;}
     break;
 
   case 203:
-#line 2392 "Gmsh.y"
+#line 2393 "Gmsh.y"
     {
     ;}
     break;
 
   case 204:
-#line 2398 "Gmsh.y"
+#line 2399 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = true;
       extr.mesh.NbLayer = 1;
@@ -6397,7 +6398,7 @@ yyreduce:
     break;
 
   case 205:
-#line 2407 "Gmsh.y"
+#line 2408 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = true;
       extr.mesh.NbLayer = List_Nbr((yyvsp[(3) - (7)].l));
@@ -6420,7 +6421,7 @@ yyreduce:
     break;
 
   case 206:
-#line 2427 "Gmsh.y"
+#line 2428 "Gmsh.y"
     {
       yymsg(0, "Explicit region numbers in layers are deprecated");
       extr.mesh.ExtrudeMesh = true;
@@ -6446,14 +6447,14 @@ yyreduce:
     break;
 
   case 207:
-#line 2450 "Gmsh.y"
+#line 2451 "Gmsh.y"
     {
       extr.mesh.Recombine = true;
     ;}
     break;
 
   case 208:
-#line 2454 "Gmsh.y"
+#line 2455 "Gmsh.y"
     {
       int num = (int)(yyvsp[(3) - (9)].d);
       if(FindSurface(num)){
@@ -6475,14 +6476,14 @@ yyreduce:
     break;
 
   case 209:
-#line 2477 "Gmsh.y"
+#line 2478 "Gmsh.y"
     {
       (yyval.v)[0] = (yyval.v)[1] = 1.;
     ;}
     break;
 
   case 210:
-#line 2481 "Gmsh.y"
+#line 2482 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(2) - (3)].c), "Progression") || !strcmp((yyvsp[(2) - (3)].c), "Power"))
         (yyval.v)[0] = 1.;
@@ -6498,14 +6499,14 @@ yyreduce:
     break;
 
   case 211:
-#line 2496 "Gmsh.y"
+#line 2497 "Gmsh.y"
     {
       (yyval.i) = -1; // left
     ;}
     break;
 
   case 212:
-#line 2500 "Gmsh.y"
+#line 2501 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(1) - (1)].c), "Right"))
         (yyval.i) = 1;
@@ -6518,35 +6519,35 @@ yyreduce:
     break;
 
   case 213:
-#line 2512 "Gmsh.y"
+#line 2513 "Gmsh.y"
     {
      (yyval.l) = List_Create(1, 1, sizeof(double));
    ;}
     break;
 
   case 214:
-#line 2516 "Gmsh.y"
+#line 2517 "Gmsh.y"
     {
      (yyval.l) = (yyvsp[(2) - (2)].l);
    ;}
     break;
 
   case 215:
-#line 2521 "Gmsh.y"
+#line 2522 "Gmsh.y"
     {
       (yyval.i) = 45;
     ;}
     break;
 
   case 216:
-#line 2525 "Gmsh.y"
+#line 2526 "Gmsh.y"
     {
       (yyval.i) = (int)(yyvsp[(2) - (2)].d);
     ;}
     break;
 
   case 217:
-#line 2532 "Gmsh.y"
+#line 2533 "Gmsh.y"
     {
       int type = (int)(yyvsp[(6) - (7)].v)[0];
       double coef = fabs((yyvsp[(6) - (7)].v)[1]);
@@ -6605,7 +6606,7 @@ yyreduce:
     break;
 
   case 218:
-#line 2588 "Gmsh.y"
+#line 2589 "Gmsh.y"
     {
       int k = List_Nbr((yyvsp[(4) - (6)].l));
       if(k != 0 && k != 3 && k != 4){
@@ -6678,7 +6679,7 @@ yyreduce:
     break;
 
   case 219:
-#line 2658 "Gmsh.y"
+#line 2659 "Gmsh.y"
     {
       yymsg(1, "Elliptic Surface is deprecated: use Transfinite instead (with smoothing)");
       List_Delete((yyvsp[(7) - (8)].l));
@@ -6686,7 +6687,7 @@ yyreduce:
     break;
 
   case 220:
-#line 2663 "Gmsh.y"
+#line 2664 "Gmsh.y"
     {
       int k = List_Nbr((yyvsp[(4) - (5)].l));
       if(k != 0 && k != 6 && k != 8){
@@ -6756,7 +6757,7 @@ yyreduce:
     break;
 
   case 221:
-#line 2730 "Gmsh.y"
+#line 2731 "Gmsh.y"
     {
       if(!(yyvsp[(3) - (5)].l)){
 	List_T *tmp = Tree2List(GModel::current()->getGEOInternals()->Surfaces);
@@ -6802,7 +6803,7 @@ yyreduce:
     break;
 
   case 222:
-#line 2773 "Gmsh.y"
+#line 2774 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (6)].l)); i++){
 	double d;
@@ -6825,7 +6826,7 @@ yyreduce:
     break;
 
   case 223:
-#line 2799 "Gmsh.y"
+#line 2800 "Gmsh.y"
     { 
       Surface *s = FindSurface((int)(yyvsp[(8) - (10)].d));
       if(s){
@@ -6851,7 +6852,7 @@ yyreduce:
     break;
 
   case 224:
-#line 2822 "Gmsh.y"
+#line 2823 "Gmsh.y"
     {
       Surface *s = FindSurface((int)(yyvsp[(8) - (10)].d));
       if(s){
@@ -6877,26 +6878,26 @@ yyreduce:
     break;
 
   case 225:
-#line 2845 "Gmsh.y"
+#line 2846 "Gmsh.y"
     {
     ;}
     break;
 
   case 226:
-#line 2848 "Gmsh.y"
+#line 2849 "Gmsh.y"
     {
     ;}
     break;
 
   case 227:
-#line 2857 "Gmsh.y"
+#line 2858 "Gmsh.y"
     { 
       ReplaceAllDuplicates();
     ;}
     break;
 
   case 228:
-#line 2861 "Gmsh.y"
+#line 2862 "Gmsh.y"
     { 
       if(!strcmp((yyvsp[(2) - (3)].c), "Geometry"))
         ReplaceAllDuplicates();
@@ -6909,47 +6910,47 @@ yyreduce:
     break;
 
   case 229:
-#line 2876 "Gmsh.y"
+#line 2877 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (1)].d);           ;}
     break;
 
   case 230:
-#line 2877 "Gmsh.y"
+#line 2878 "Gmsh.y"
     { (yyval.d) = (yyvsp[(2) - (3)].d);           ;}
     break;
 
   case 231:
-#line 2878 "Gmsh.y"
+#line 2879 "Gmsh.y"
     { (yyval.d) = -(yyvsp[(2) - (2)].d);          ;}
     break;
 
   case 232:
-#line 2879 "Gmsh.y"
+#line 2880 "Gmsh.y"
     { (yyval.d) = (yyvsp[(2) - (2)].d);           ;}
     break;
 
   case 233:
-#line 2880 "Gmsh.y"
+#line 2881 "Gmsh.y"
     { (yyval.d) = !(yyvsp[(2) - (2)].d);          ;}
     break;
 
   case 234:
-#line 2881 "Gmsh.y"
+#line 2882 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) - (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 235:
-#line 2882 "Gmsh.y"
+#line 2883 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) + (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 236:
-#line 2883 "Gmsh.y"
+#line 2884 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) * (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 237:
-#line 2885 "Gmsh.y"
+#line 2886 "Gmsh.y"
     { 
       if(!(yyvsp[(3) - (3)].d))
 	yymsg(0, "Division by zero in '%g / %g'", (yyvsp[(1) - (3)].d), (yyvsp[(3) - (3)].d));
@@ -6959,307 +6960,307 @@ yyreduce:
     break;
 
   case 238:
-#line 2891 "Gmsh.y"
+#line 2892 "Gmsh.y"
     { (yyval.d) = (int)(yyvsp[(1) - (3)].d) % (int)(yyvsp[(3) - (3)].d);  ;}
     break;
 
   case 239:
-#line 2892 "Gmsh.y"
+#line 2893 "Gmsh.y"
     { (yyval.d) = pow((yyvsp[(1) - (3)].d), (yyvsp[(3) - (3)].d));  ;}
     break;
 
   case 240:
-#line 2893 "Gmsh.y"
+#line 2894 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) < (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 241:
-#line 2894 "Gmsh.y"
+#line 2895 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) > (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 242:
-#line 2895 "Gmsh.y"
+#line 2896 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) <= (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 243:
-#line 2896 "Gmsh.y"
+#line 2897 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) >= (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 244:
-#line 2897 "Gmsh.y"
+#line 2898 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) == (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 245:
-#line 2898 "Gmsh.y"
+#line 2899 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) != (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 246:
-#line 2899 "Gmsh.y"
+#line 2900 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) && (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 247:
-#line 2900 "Gmsh.y"
+#line 2901 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) || (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 248:
-#line 2901 "Gmsh.y"
+#line 2902 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (5)].d) ? (yyvsp[(3) - (5)].d) : (yyvsp[(5) - (5)].d); ;}
     break;
 
   case 249:
-#line 2902 "Gmsh.y"
+#line 2903 "Gmsh.y"
     { (yyval.d) = exp((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 250:
-#line 2903 "Gmsh.y"
+#line 2904 "Gmsh.y"
     { (yyval.d) = log((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 251:
-#line 2904 "Gmsh.y"
+#line 2905 "Gmsh.y"
     { (yyval.d) = log10((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 252:
-#line 2905 "Gmsh.y"
+#line 2906 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 253:
-#line 2906 "Gmsh.y"
+#line 2907 "Gmsh.y"
     { (yyval.d) = sin((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 254:
-#line 2907 "Gmsh.y"
+#line 2908 "Gmsh.y"
     { (yyval.d) = asin((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 255:
-#line 2908 "Gmsh.y"
+#line 2909 "Gmsh.y"
     { (yyval.d) = cos((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 256:
-#line 2909 "Gmsh.y"
+#line 2910 "Gmsh.y"
     { (yyval.d) = acos((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 257:
-#line 2910 "Gmsh.y"
+#line 2911 "Gmsh.y"
     { (yyval.d) = tan((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 258:
-#line 2911 "Gmsh.y"
+#line 2912 "Gmsh.y"
     { (yyval.d) = atan((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 259:
-#line 2912 "Gmsh.y"
+#line 2913 "Gmsh.y"
     { (yyval.d) = atan2((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d));;}
     break;
 
   case 260:
-#line 2913 "Gmsh.y"
+#line 2914 "Gmsh.y"
     { (yyval.d) = sinh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 261:
-#line 2914 "Gmsh.y"
+#line 2915 "Gmsh.y"
     { (yyval.d) = cosh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 262:
-#line 2915 "Gmsh.y"
+#line 2916 "Gmsh.y"
     { (yyval.d) = tanh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 263:
-#line 2916 "Gmsh.y"
+#line 2917 "Gmsh.y"
     { (yyval.d) = fabs((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 264:
-#line 2917 "Gmsh.y"
+#line 2918 "Gmsh.y"
     { (yyval.d) = floor((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 265:
-#line 2918 "Gmsh.y"
+#line 2919 "Gmsh.y"
     { (yyval.d) = ceil((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 266:
-#line 2919 "Gmsh.y"
+#line 2920 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 267:
-#line 2920 "Gmsh.y"
+#line 2921 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 268:
-#line 2921 "Gmsh.y"
+#line 2922 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (6)].d) * (yyvsp[(3) - (6)].d) + (yyvsp[(5) - (6)].d) * (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 269:
-#line 2922 "Gmsh.y"
+#line 2923 "Gmsh.y"
     { (yyval.d) = (yyvsp[(3) - (4)].d) * (double)rand() / (double)RAND_MAX; ;}
     break;
 
   case 270:
-#line 2924 "Gmsh.y"
+#line 2925 "Gmsh.y"
     { (yyval.d) = exp((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 271:
-#line 2925 "Gmsh.y"
+#line 2926 "Gmsh.y"
     { (yyval.d) = log((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 272:
-#line 2926 "Gmsh.y"
+#line 2927 "Gmsh.y"
     { (yyval.d) = log10((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 273:
-#line 2927 "Gmsh.y"
+#line 2928 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 274:
-#line 2928 "Gmsh.y"
+#line 2929 "Gmsh.y"
     { (yyval.d) = sin((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 275:
-#line 2929 "Gmsh.y"
+#line 2930 "Gmsh.y"
     { (yyval.d) = asin((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 276:
-#line 2930 "Gmsh.y"
+#line 2931 "Gmsh.y"
     { (yyval.d) = cos((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 277:
-#line 2931 "Gmsh.y"
+#line 2932 "Gmsh.y"
     { (yyval.d) = acos((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 278:
-#line 2932 "Gmsh.y"
+#line 2933 "Gmsh.y"
     { (yyval.d) = tan((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 279:
-#line 2933 "Gmsh.y"
+#line 2934 "Gmsh.y"
     { (yyval.d) = atan((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 280:
-#line 2934 "Gmsh.y"
+#line 2935 "Gmsh.y"
     { (yyval.d) = atan2((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d));;}
     break;
 
   case 281:
-#line 2935 "Gmsh.y"
+#line 2936 "Gmsh.y"
     { (yyval.d) = sinh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 282:
-#line 2936 "Gmsh.y"
+#line 2937 "Gmsh.y"
     { (yyval.d) = cosh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 283:
-#line 2937 "Gmsh.y"
+#line 2938 "Gmsh.y"
     { (yyval.d) = tanh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 284:
-#line 2938 "Gmsh.y"
+#line 2939 "Gmsh.y"
     { (yyval.d) = fabs((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 285:
-#line 2939 "Gmsh.y"
+#line 2940 "Gmsh.y"
     { (yyval.d) = floor((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 286:
-#line 2940 "Gmsh.y"
+#line 2941 "Gmsh.y"
     { (yyval.d) = ceil((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 287:
-#line 2941 "Gmsh.y"
+#line 2942 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 288:
-#line 2942 "Gmsh.y"
+#line 2943 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 289:
-#line 2943 "Gmsh.y"
+#line 2944 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (6)].d) * (yyvsp[(3) - (6)].d) + (yyvsp[(5) - (6)].d) * (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 290:
-#line 2944 "Gmsh.y"
+#line 2945 "Gmsh.y"
     { (yyval.d) = (yyvsp[(3) - (4)].d) * (double)rand() / (double)RAND_MAX; ;}
     break;
 
   case 291:
-#line 2953 "Gmsh.y"
+#line 2954 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (1)].d); ;}
     break;
 
   case 292:
-#line 2954 "Gmsh.y"
+#line 2955 "Gmsh.y"
     { (yyval.d) = 3.141592653589793; ;}
     break;
 
   case 293:
-#line 2955 "Gmsh.y"
+#line 2956 "Gmsh.y"
     { (yyval.d) = Msg::GetCommRank(); ;}
     break;
 
   case 294:
-#line 2956 "Gmsh.y"
+#line 2957 "Gmsh.y"
     { (yyval.d) = Msg::GetCommSize(); ;}
     break;
 
   case 295:
-#line 2957 "Gmsh.y"
+#line 2958 "Gmsh.y"
     { (yyval.d) = GetGmshMajorVersion(); ;}
     break;
 
   case 296:
-#line 2958 "Gmsh.y"
+#line 2959 "Gmsh.y"
     { (yyval.d) = GetGmshMinorVersion(); ;}
     break;
 
   case 297:
-#line 2959 "Gmsh.y"
+#line 2960 "Gmsh.y"
     { (yyval.d) = GetGmshPatchVersion(); ;}
     break;
 
   case 298:
-#line 2964 "Gmsh.y"
+#line 2965 "Gmsh.y"
     {
       if(!gmsh_yysymbols.count((yyvsp[(1) - (1)].c))){
 	yymsg(0, "Unknown variable '%s'", (yyvsp[(1) - (1)].c));
@@ -7272,7 +7273,7 @@ yyreduce:
     break;
 
   case 299:
-#line 2977 "Gmsh.y"
+#line 2978 "Gmsh.y"
     {
       char tmpstring[1024];
       sprintf(tmpstring, "%s_%d", (yyvsp[(1) - (5)].c), (int)(yyvsp[(4) - (5)].d)) ;
@@ -7287,7 +7288,7 @@ yyreduce:
     break;
 
   case 300:
-#line 2989 "Gmsh.y"
+#line 2990 "Gmsh.y"
     {
       int index = (int)(yyvsp[(3) - (4)].d);
       if(!gmsh_yysymbols.count((yyvsp[(1) - (4)].c))){
@@ -7305,7 +7306,7 @@ yyreduce:
     break;
 
   case 301:
-#line 3004 "Gmsh.y"
+#line 3005 "Gmsh.y"
     {
       if(!gmsh_yysymbols.count((yyvsp[(2) - (4)].c))){
 	yymsg(0, "Unknown variable '%s'", (yyvsp[(2) - (4)].c));
@@ -7318,7 +7319,7 @@ yyreduce:
     break;
 
   case 302:
-#line 3014 "Gmsh.y"
+#line 3015 "Gmsh.y"
     {
       if(!gmsh_yysymbols.count((yyvsp[(1) - (2)].c))){
 	yymsg(0, "Unknown variable '%s'", (yyvsp[(1) - (2)].c));
@@ -7331,7 +7332,7 @@ yyreduce:
     break;
 
   case 303:
-#line 3024 "Gmsh.y"
+#line 3025 "Gmsh.y"
     {
       int index = (int)(yyvsp[(3) - (5)].d);
       if(!gmsh_yysymbols.count((yyvsp[(1) - (5)].c))){
@@ -7349,7 +7350,7 @@ yyreduce:
     break;
 
   case 304:
-#line 3042 "Gmsh.y"
+#line 3043 "Gmsh.y"
     {
       NumberOption(GMSH_GET, (yyvsp[(1) - (3)].c), 0, (yyvsp[(3) - (3)].c), (yyval.d));
       Free((yyvsp[(1) - (3)].c)); Free((yyvsp[(3) - (3)].c));
@@ -7357,7 +7358,7 @@ yyreduce:
     break;
 
   case 305:
-#line 3047 "Gmsh.y"
+#line 3048 "Gmsh.y"
     {
       NumberOption(GMSH_GET, (yyvsp[(1) - (6)].c), (int)(yyvsp[(3) - (6)].d), (yyvsp[(6) - (6)].c), (yyval.d));
       Free((yyvsp[(1) - (6)].c)); Free((yyvsp[(6) - (6)].c));
@@ -7365,7 +7366,7 @@ yyreduce:
     break;
 
   case 306:
-#line 3052 "Gmsh.y"
+#line 3053 "Gmsh.y"
     {
       double d = 0.;
       if(NumberOption(GMSH_GET, (yyvsp[(1) - (4)].c), 0, (yyvsp[(3) - (4)].c), d)){
@@ -7378,7 +7379,7 @@ yyreduce:
     break;
 
   case 307:
-#line 3062 "Gmsh.y"
+#line 3063 "Gmsh.y"
     {
       double d = 0.;
       if(NumberOption(GMSH_GET, (yyvsp[(1) - (7)].c), (int)(yyvsp[(3) - (7)].d), (yyvsp[(6) - (7)].c), d)){
@@ -7391,7 +7392,7 @@ yyreduce:
     break;
 
   case 308:
-#line 3072 "Gmsh.y"
+#line 3073 "Gmsh.y"
     { 
       (yyval.d) = Msg::GetValue((yyvsp[(3) - (6)].c), (yyvsp[(5) - (6)].d));
       Free((yyvsp[(3) - (6)].c));
@@ -7399,70 +7400,70 @@ yyreduce:
     break;
 
   case 309:
-#line 3080 "Gmsh.y"
+#line 3081 "Gmsh.y"
     {
       memcpy((yyval.v), (yyvsp[(1) - (1)].v), 5*sizeof(double));
     ;}
     break;
 
   case 310:
-#line 3084 "Gmsh.y"
+#line 3085 "Gmsh.y"
     {
       for(int i = 0; i < 5; i++) (yyval.v)[i] = -(yyvsp[(2) - (2)].v)[i];
     ;}
     break;
 
   case 311:
-#line 3088 "Gmsh.y"
+#line 3089 "Gmsh.y"
     { 
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(2) - (2)].v)[i];
     ;}
     break;
 
   case 312:
-#line 3092 "Gmsh.y"
+#line 3093 "Gmsh.y"
     { 
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(1) - (3)].v)[i] - (yyvsp[(3) - (3)].v)[i];
     ;}
     break;
 
   case 313:
-#line 3096 "Gmsh.y"
+#line 3097 "Gmsh.y"
     {
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(1) - (3)].v)[i] + (yyvsp[(3) - (3)].v)[i];
     ;}
     break;
 
   case 314:
-#line 3103 "Gmsh.y"
+#line 3104 "Gmsh.y"
     { 
       (yyval.v)[0] = (yyvsp[(2) - (11)].d);  (yyval.v)[1] = (yyvsp[(4) - (11)].d);  (yyval.v)[2] = (yyvsp[(6) - (11)].d);  (yyval.v)[3] = (yyvsp[(8) - (11)].d); (yyval.v)[4] = (yyvsp[(10) - (11)].d);
     ;}
     break;
 
   case 315:
-#line 3107 "Gmsh.y"
+#line 3108 "Gmsh.y"
     { 
       (yyval.v)[0] = (yyvsp[(2) - (9)].d);  (yyval.v)[1] = (yyvsp[(4) - (9)].d);  (yyval.v)[2] = (yyvsp[(6) - (9)].d);  (yyval.v)[3] = (yyvsp[(8) - (9)].d); (yyval.v)[4] = 1.0;
     ;}
     break;
 
   case 316:
-#line 3111 "Gmsh.y"
+#line 3112 "Gmsh.y"
     {
       (yyval.v)[0] = (yyvsp[(2) - (7)].d);  (yyval.v)[1] = (yyvsp[(4) - (7)].d);  (yyval.v)[2] = (yyvsp[(6) - (7)].d);  (yyval.v)[3] = 0.0; (yyval.v)[4] = 1.0;
     ;}
     break;
 
   case 317:
-#line 3115 "Gmsh.y"
+#line 3116 "Gmsh.y"
     {
       (yyval.v)[0] = (yyvsp[(2) - (7)].d);  (yyval.v)[1] = (yyvsp[(4) - (7)].d);  (yyval.v)[2] = (yyvsp[(6) - (7)].d);  (yyval.v)[3] = 0.0; (yyval.v)[4] = 1.0;
     ;}
     break;
 
   case 318:
-#line 3122 "Gmsh.y"
+#line 3123 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(List_T*));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].l)));
@@ -7470,14 +7471,14 @@ yyreduce:
     break;
 
   case 319:
-#line 3127 "Gmsh.y"
+#line 3128 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].l)));
     ;}
     break;
 
   case 320:
-#line 3134 "Gmsh.y"
+#line 3135 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].d)));
@@ -7485,14 +7486,14 @@ yyreduce:
     break;
 
   case 321:
-#line 3139 "Gmsh.y"
+#line 3140 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(1) - (1)].l);
     ;}
     break;
 
   case 322:
-#line 3143 "Gmsh.y"
+#line 3144 "Gmsh.y"
     {
       // creates an empty list
       (yyval.l) = List_Create(2, 1, sizeof(double));
@@ -7500,14 +7501,14 @@ yyreduce:
     break;
 
   case 323:
-#line 3148 "Gmsh.y"
+#line 3149 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (3)].l);
     ;}
     break;
 
   case 324:
-#line 3152 "Gmsh.y"
+#line 3153 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(3) - (4)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7518,7 +7519,7 @@ yyreduce:
     break;
 
   case 325:
-#line 3160 "Gmsh.y"
+#line 3161 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(4) - (5)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7529,14 +7530,14 @@ yyreduce:
     break;
 
   case 326:
-#line 3171 "Gmsh.y"
+#line 3172 "Gmsh.y"
     { 
       (yyval.l) = (yyvsp[(1) - (1)].l); 
     ;}
     break;
 
   case 327:
-#line 3175 "Gmsh.y"
+#line 3176 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(1) - (1)].c), "*") || !strcmp((yyvsp[(1) - (1)].c), "all"))
         (yyval.l) = 0;
@@ -7548,7 +7549,7 @@ yyreduce:
     break;
 
   case 328:
-#line 3187 "Gmsh.y"
+#line 3188 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (2)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7559,7 +7560,7 @@ yyreduce:
     break;
 
   case 329:
-#line 3195 "Gmsh.y"
+#line 3196 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(3) - (3)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7570,7 +7571,7 @@ yyreduce:
     break;
 
   case 330:
-#line 3203 "Gmsh.y"
+#line 3204 "Gmsh.y"
     { 
       (yyval.l) = List_Create(2, 1, sizeof(double)); 
       for(double d = (yyvsp[(1) - (3)].d); ((yyvsp[(1) - (3)].d) < (yyvsp[(3) - (3)].d)) ? (d <= (yyvsp[(3) - (3)].d)) : (d >= (yyvsp[(3) - (3)].d)); 
@@ -7580,7 +7581,7 @@ yyreduce:
     break;
 
   case 331:
-#line 3210 "Gmsh.y"
+#line 3211 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double)); 
       if(!(yyvsp[(5) - (5)].d) || ((yyvsp[(1) - (5)].d) < (yyvsp[(3) - (5)].d) && (yyvsp[(5) - (5)].d) < 0) || ((yyvsp[(1) - (5)].d) > (yyvsp[(3) - (5)].d) && (yyvsp[(5) - (5)].d) > 0)){
@@ -7594,7 +7595,7 @@ yyreduce:
     break;
 
   case 332:
-#line 3221 "Gmsh.y"
+#line 3222 "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
@@ -7617,7 +7618,7 @@ yyreduce:
     break;
 
   case 333:
-#line 3241 "Gmsh.y"
+#line 3242 "Gmsh.y"
     {
       (yyval.l) = List_Create(List_Nbr((yyvsp[(1) - (1)].l)), 1, sizeof(double));
       for(int i = 0; i < List_Nbr((yyvsp[(1) - (1)].l)); i++){
@@ -7630,7 +7631,7 @@ yyreduce:
     break;
 
   case 334:
-#line 3251 "Gmsh.y"
+#line 3252 "Gmsh.y"
     {
       (yyval.l) = List_Create(List_Nbr((yyvsp[(1) - (1)].l)), 1, sizeof(double));
       for(int i = 0; i < List_Nbr((yyvsp[(1) - (1)].l)); i++){
@@ -7643,7 +7644,7 @@ yyreduce:
     break;
 
   case 335:
-#line 3261 "Gmsh.y"
+#line 3262 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       if(!gmsh_yysymbols.count((yyvsp[(1) - (3)].c)))
@@ -7656,7 +7657,7 @@ yyreduce:
     break;
 
   case 336:
-#line 3271 "Gmsh.y"
+#line 3272 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       if(!gmsh_yysymbols.count((yyvsp[(1) - (6)].c)))
@@ -7676,7 +7677,7 @@ yyreduce:
     break;
 
   case 337:
-#line 3291 "Gmsh.y"
+#line 3292 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].d)));
@@ -7684,21 +7685,21 @@ yyreduce:
     break;
 
   case 338:
-#line 3296 "Gmsh.y"
+#line 3297 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(1) - (1)].l);
     ;}
     break;
 
   case 339:
-#line 3300 "Gmsh.y"
+#line 3301 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].d)));
     ;}
     break;
 
   case 340:
-#line 3304 "Gmsh.y"
+#line 3305 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (3)].l)); i++){
 	double d;
@@ -7710,21 +7711,21 @@ yyreduce:
     break;
 
   case 341:
-#line 3316 "Gmsh.y"
+#line 3317 "Gmsh.y"
     {
       (yyval.u) = CTX::instance()->packColor((int)(yyvsp[(2) - (9)].d), (int)(yyvsp[(4) - (9)].d), (int)(yyvsp[(6) - (9)].d), (int)(yyvsp[(8) - (9)].d));
     ;}
     break;
 
   case 342:
-#line 3320 "Gmsh.y"
+#line 3321 "Gmsh.y"
     {
       (yyval.u) = CTX::instance()->packColor((int)(yyvsp[(2) - (7)].d), (int)(yyvsp[(4) - (7)].d), (int)(yyvsp[(6) - (7)].d), 255);
     ;}
     break;
 
   case 343:
-#line 3332 "Gmsh.y"
+#line 3333 "Gmsh.y"
     {
       int flag;
       (yyval.u) = GetColorForString(ColorString, -1, (yyvsp[(1) - (1)].c), &flag);
@@ -7734,7 +7735,7 @@ yyreduce:
     break;
 
   case 344:
-#line 3339 "Gmsh.y"
+#line 3340 "Gmsh.y"
     {
       unsigned int val = 0;
       ColorOption(GMSH_GET, (yyvsp[(1) - (5)].c), 0, (yyvsp[(5) - (5)].c), val);
@@ -7744,14 +7745,14 @@ yyreduce:
     break;
 
   case 345:
-#line 3349 "Gmsh.y"
+#line 3350 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (3)].l);
     ;}
     break;
 
   case 346:
-#line 3353 "Gmsh.y"
+#line 3354 "Gmsh.y"
     {
       (yyval.l) = List_Create(256, 10, sizeof(unsigned int));
       GmshColorTable *ct = GetColorTable((int)(yyvsp[(3) - (6)].d));
@@ -7766,7 +7767,7 @@ yyreduce:
     break;
 
   case 347:
-#line 3368 "Gmsh.y"
+#line 3369 "Gmsh.y"
     {
       (yyval.l) = List_Create(256, 10, sizeof(unsigned int));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].u)));
@@ -7774,21 +7775,21 @@ yyreduce:
     break;
 
   case 348:
-#line 3373 "Gmsh.y"
+#line 3374 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].u)));
     ;}
     break;
 
   case 349:
-#line 3380 "Gmsh.y"
+#line 3381 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(1) - (1)].c);
     ;}
     break;
 
   case 350:
-#line 3384 "Gmsh.y"
+#line 3385 "Gmsh.y"
     {
       if(!gmsh_yystringsymbols.count((yyvsp[(1) - (1)].c))){
 	yymsg(0, "Unknown string variable '%s'", (yyvsp[(1) - (1)].c));
@@ -7804,7 +7805,7 @@ yyreduce:
     break;
 
   case 351:
-#line 3397 "Gmsh.y"
+#line 3398 "Gmsh.y"
     { 
       std::string out;
       StringOption(GMSH_GET, (yyvsp[(1) - (3)].c), 0, (yyvsp[(3) - (3)].c), out);
@@ -7815,7 +7816,7 @@ yyreduce:
     break;
 
   case 352:
-#line 3405 "Gmsh.y"
+#line 3406 "Gmsh.y"
     { 
       std::string out;
       StringOption(GMSH_GET, (yyvsp[(1) - (6)].c), (int)(yyvsp[(3) - (6)].d), (yyvsp[(6) - (6)].c), out);
@@ -7826,14 +7827,14 @@ yyreduce:
     break;
 
   case 353:
-#line 3416 "Gmsh.y"
+#line 3417 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(1) - (1)].c);
     ;}
     break;
 
   case 354:
-#line 3420 "Gmsh.y"
+#line 3421 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc(32 * sizeof(char));
       time_t now;
@@ -7844,7 +7845,7 @@ yyreduce:
     break;
 
   case 355:
-#line 3428 "Gmsh.y"
+#line 3429 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc((strlen((yyvsp[(3) - (6)].c)) + strlen((yyvsp[(5) - (6)].c)) + 1) * sizeof(char));
       strcpy((yyval.c), (yyvsp[(3) - (6)].c));
@@ -7855,7 +7856,7 @@ yyreduce:
     break;
 
   case 356:
-#line 3436 "Gmsh.y"
+#line 3437 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc((strlen((yyvsp[(3) - (4)].c)) + 1) * sizeof(char));
       int i;
@@ -7872,7 +7873,7 @@ yyreduce:
     break;
 
   case 357:
-#line 3450 "Gmsh.y"
+#line 3451 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc((strlen((yyvsp[(3) - (4)].c)) + 1) * sizeof(char));
       int i;
@@ -7889,14 +7890,14 @@ yyreduce:
     break;
 
   case 358:
-#line 3464 "Gmsh.y"
+#line 3465 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(3) - (4)].c);
     ;}
     break;
 
   case 359:
-#line 3468 "Gmsh.y"
+#line 3469 "Gmsh.y"
     {
       char tmpstring[1024];
       int i = PrintListOfDouble((yyvsp[(3) - (6)].c), (yyvsp[(5) - (6)].l), tmpstring);
@@ -7919,7 +7920,7 @@ yyreduce:
 
 
 /* Line 1267 of yacc.c.  */
-#line 7923 "Gmsh.tab.cpp"
+#line 7924 "Gmsh.tab.cpp"
       default: break;
     }
   YY_SYMBOL_PRINT ("-> $$ =", yyr1[yyn], &yyval, &yyloc);
@@ -8133,7 +8134,7 @@ yyreturn:
 }
 
 
-#line 3488 "Gmsh.y"
+#line 3489 "Gmsh.y"
 
 
 int PrintListOfDouble(char *format, List_T *list, char *buffer)
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index f0e19e4546..43768c092b 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1892,6 +1892,7 @@ Command :
 	char tmpstring[1024];
 	FixRelativePath($2, tmpstring);
 	MergeFile(tmpstring, true);
+	GModel::current()->createTopologyFromMSH();
       }
       else if(!strcmp($1, "System"))
 	SystemCall($2);
-- 
GitLab