diff --git a/Box/Box.cpp b/Box/Box.cpp
index 61fde2e4a39a7384e511204f91515c874bcb3bf3..53a63669fd002e670ce018bde583e705637dde39 100644
--- a/Box/Box.cpp
+++ b/Box/Box.cpp
@@ -1,4 +1,4 @@
-// $Id: Box.cpp,v 1.49 2008-03-20 11:44:01 geuzaine Exp $
+// $Id: Box.cpp,v 1.50 2008-03-23 21:42:56 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -87,13 +87,10 @@ int GMSHBOX(int argc, char *argv[])
       MergeFile(CTX.files[i].c_str());
     if(CTX.bgm_filename) {
       MergeFile(CTX.bgm_filename);
-      if(PView::list.size()){
-        GModel::current()->getFields()->set_background_mesh(PView::list.back()->getNum() - 1);
-      }
-      else{
+      if(PView::list.size())
+        GModel::current()->getFields()->set_background_mesh(PView::list.size() - 1);
+      else
         fprintf(stderr, ERROR_STR "Invalid background mesh (no view)\n");
-        exit(1);
-      }
     }
     if(CTX.batch > 0) {
       GModel::current()->mesh(CTX.batch);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index e7ce1930f3829b8678af054ddfc95cffad989d5f..9259a39df75b681c91f39c665ed800fc95f254bb 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.572 2008-03-21 18:27:38 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.573 2008-03-23 21:42:56 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -664,7 +664,7 @@ int _save_options(const char *name){ return options_dialog(name); }
 int _save_geo(const char *name){ return geo_dialog(name); }
 int _save_cgns(const char *name){ CreateOutputFile(name, FORMAT_CGNS); return 1; }
 int _save_unv(const char *name){ return unv_dialog(name); }
-int _save_med(const char *name){ return generic_mesh_dialog(name, "MED Options", FORMAT_MED); }
+int _save_med(const char *name){ CreateOutputFile(name, FORMAT_MED); return 1; }
 int _save_mesh(const char *name){ return generic_mesh_dialog(name, "MESH Options", FORMAT_MESH); }
 int _save_bdf(const char *name){ return bdf_dialog(name); }
 int _save_p3d(const char *name){ return generic_mesh_dialog(name, "P3D Options", FORMAT_P3D); }
diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp
index 9af7058daa640e6d290ee02acbdacaa9461b0e11..a939bf2a689f44a67289ca1ba11a9b197fc74093 100644
--- a/Fltk/Main.cpp
+++ b/Fltk/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.125 2008-03-20 11:44:03 geuzaine Exp $
+// $Id: Main.cpp,v 1.126 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -87,9 +87,8 @@ int main(int argc, char *argv[])
         PView::combine(true, 2, CTX.post.combine_remove_orig);
       if(CTX.bgm_filename) {
         MergeFile(CTX.bgm_filename);
-        if(PView::list.size()){
-          GModel::current()->getFields()->set_background_mesh(PView::list.back()->getNum() - 1);
-        }
+        if(PView::list.size())
+          GModel::current()->getFields()->set_background_mesh(PView::list.size() - 1);
         else
           Msg(GERROR, "Invalid background mesh (no view)");
       }
@@ -168,7 +167,7 @@ int main(int argc, char *argv[])
   if(CTX.bgm_filename) {
     MergeFile(CTX.bgm_filename);
     if(PView::list.size())
-      GModel::current()->getFields()->set_background_mesh(PView::list.back()->getNum() - 1);
+      GModel::current()->getFields()->set_background_mesh(PView::list.size() - 1);
     else
       Msg(GERROR, "Invalid background mesh (no view)");
   }
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 8c86186d7ba77777d42f9b2d6b7a0c6b7e199c9d..0d420918fa1011fe18a6e478472f753d31794055 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.76 2008-03-20 11:44:04 geuzaine Exp $
+// $Id: GModel.cpp,v 1.77 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -22,6 +22,10 @@
 #include <string.h>
 #include "GModel.h"
 #include "MElement.h"
+#include "discreteRegion.h"
+#include "discreteFace.h"
+#include "discreteEdge.h"
+#include "discreteVertex.h"
 
 #if defined(HAVE_GMSH_EMBEDDED)
 #  include "GmshEmbedded.h"
@@ -43,7 +47,7 @@ GModel::GModel(std::string name)
 {
   list.push_back(this);
   // at the moment we always create (at least an empty) GEO model
-  createGEOInternals();
+  _createGEOInternals();
 
 #if !defined(HAVE_GMSH_EMBEDDED)
   _fields = new FieldManager();
@@ -54,8 +58,8 @@ GModel::~GModel()
 {
   std::vector<GModel*>::iterator it = std::find(list.begin(), list.end(), this);
   if(it != list.end()) list.erase(it);
-  deleteGEOInternals();
-  deleteOCCInternals();
+  _deleteGEOInternals();
+  _deleteOCCInternals();
   destroy();
 #if !defined(HAVE_GMSH_EMBEDDED)
   delete _fields;
@@ -485,38 +489,6 @@ void GModel::removeInvisibleElements()
   }
 }
 
-template<class T>
-static void associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &elements)
-{
-  for(unsigned int i = 0; i < elements.size(); i++)
-    for(int j = 0; j < elements[i]->getNumVertices(); j++)
-      elements[i]->getVertex(j)->setEntity(ge);
-}
-
-void GModel::associateEntityWithMeshVertices()
-{
-  // loop on regions, then on faces, edges and vertices and store the
-  // entity pointer in the the elements' vertices (this way we
-  // associate the entity of lowest geometrical degree with each
-  // vertex)
-  for(riter it = firstRegion(); it != lastRegion(); ++it){
-    associateEntityWithElementVertices(*it, (*it)->tetrahedra);
-    associateEntityWithElementVertices(*it, (*it)->hexahedra);
-    associateEntityWithElementVertices(*it, (*it)->prisms);
-    associateEntityWithElementVertices(*it, (*it)->pyramids);
-  }
-  for(fiter it = firstFace(); it != lastFace(); ++it){
-    associateEntityWithElementVertices(*it, (*it)->triangles);
-    associateEntityWithElementVertices(*it, (*it)->quadrangles);
-  }
-  for(eiter it = firstEdge(); it != lastEdge(); ++it){
-    associateEntityWithElementVertices(*it, (*it)->lines);
-  }
-  for(viter it = firstVertex(); it != lastVertex(); ++it){
-    (*it)->mesh_vertices[0]->setEntity(*it);
-  }
-}
-
 int GModel::renumberMeshVertices(bool saveAll)
 {
   _vertexVectorCache.clear();
@@ -724,3 +696,86 @@ void GModel::checkMeshCoherence()
     MElementLessThanLexicographic::tolerance = old_tol;
   }
 }
+
+template<class T>
+static void _addElements(std::vector<T*> &dst, const std::vector<MElement*> &src)
+{
+  for(unsigned int i = 0; i < src.size(); i++) dst.push_back((T*)src[i]);
+}
+
+void GModel::_storeElementsInEntities(std::map<int, std::vector<MElement*> > &map)
+{
+  std::map<int, std::vector<MElement*> >::const_iterator it = map.begin();
+  for(; it != map.end(); ++it){
+    if(!it->second.size()) continue;
+    int numEdges = it->second[0]->getNumEdges();
+    switch(numEdges){
+    case 1: 
+      {
+        GEdge *e = getEdgeByTag(it->first);
+        if(!e){
+          e = new discreteEdge(this, it->first);
+          add(e);
+        }
+        _addElements(e->lines, it->second);
+      }
+      break;
+    case 3: case 4: 
+      {
+        GFace *f = getFaceByTag(it->first);
+        if(!f){
+          f = new discreteFace(this, it->first);
+          add(f);
+        }
+        if(numEdges == 3) _addElements(f->triangles, it->second);
+        else _addElements(f->quadrangles, it->second);
+      }
+      break;
+    case 6: case 12: case 9: case 8:
+      {
+        GRegion *r = getRegionByTag(it->first);
+        if(!r){
+          r = new discreteRegion(this, it->first);
+          add(r);
+        }
+        if(numEdges == 6) _addElements(r->tetrahedra, it->second);
+        else if(numEdges == 12) _addElements(r->hexahedra, it->second);
+        else if(numEdges == 9) _addElements(r->prisms, it->second);
+        else _addElements(r->pyramids, it->second);
+      }
+      break;
+    }
+  }
+}
+
+template<class T>
+static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &elements)
+{
+  for(unsigned int i = 0; i < elements.size(); i++)
+    for(int j = 0; j < elements[i]->getNumVertices(); j++)
+      elements[i]->getVertex(j)->setEntity(ge);
+}
+
+void GModel::_associateEntityWithMeshVertices()
+{
+  // loop on regions, then on faces, edges and vertices and store the
+  // entity pointer in the the elements' vertices (this way we
+  // associate the entity of lowest geometrical degree with each
+  // vertex)
+  for(riter it = firstRegion(); it != lastRegion(); ++it){
+    _associateEntityWithElementVertices(*it, (*it)->tetrahedra);
+    _associateEntityWithElementVertices(*it, (*it)->hexahedra);
+    _associateEntityWithElementVertices(*it, (*it)->prisms);
+    _associateEntityWithElementVertices(*it, (*it)->pyramids);
+  }
+  for(fiter it = firstFace(); it != lastFace(); ++it){
+    _associateEntityWithElementVertices(*it, (*it)->triangles);
+    _associateEntityWithElementVertices(*it, (*it)->quadrangles);
+  }
+  for(eiter it = firstEdge(); it != lastEdge(); ++it){
+    _associateEntityWithElementVertices(*it, (*it)->lines);
+  }
+  for(viter it = firstVertex(); it != lastVertex(); ++it){
+    (*it)->mesh_vertices[0]->setEntity(*it);
+  }
+}
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 11c96342d9e21d2a91ca8f663a589e91db80eb3e..f075f0e7f9821dec19a1203528a82efc3604b566 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -45,15 +45,23 @@ class GModel
   int _maxVertexDataIndex;
 
   GEO_Internals *_geo_internals;
-  void createGEOInternals();
-  void deleteGEOInternals();
+  void _createGEOInternals();
+  void _deleteGEOInternals();
 
   OCC_Internals *_occ_internals;
-  void deleteOCCInternals();
+  void _deleteOCCInternals();
 
   // Characteristic Lengths fields
   FieldManager *_fields;
 
+  // Store the elements in the map (indexed by elementary region
+  // number) into the model, creating discrete geometrical entities on
+  // the fly if needed
+  void _storeElementsInEntities(std::map<int, std::vector<MElement*> > &map);
+
+  // loop over all vertices connected to elements and associate geo entity
+  void _associateEntityWithMeshVertices();
+
  protected:
   std::string modelName;
   std::set<GRegion*, GEntityLessThan> regions;
@@ -83,6 +91,10 @@ class GModel
   // Access characteristic length fields
   FieldManager *getFields(){ return _fields; }
 
+  // Get/set the model name
+  void setName(std::string name){ modelName = name; }
+  std::string getName(){ return modelName; }
+
   // Get the number of regions in this model.
   int getNumRegions() const { return regions.size(); }
   int getNumFaces() const { return faces.size(); }
@@ -164,9 +176,6 @@ class GModel
   // Access a mesh vertex by tag, using the vertex cache
   MVertex *getMeshVertexByTag(int n);
 
-  // loop over all vertices connected to elements and associate geo entity
-  void associateEntityWithMeshVertices();
-
   // Renumber all the (used) mesh vertices in a continuous sequence
   int renumberMeshVertices(bool saveAll);
 
@@ -251,7 +260,8 @@ class GModel
   int writeCGNS(const std::string &name, double scalingFactor=1.0);
 
   // Med interface ("Modele d'Echange de Donnees")
-  int writeMED(const std::string &name);
+  int readMED(const std::string &name);
+  int writeMED(const std::string &name, double scalingFactor=1.0);
 };
 
 #endif
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 937b9e15b83a34f010075ec36d203ea1faa650b3..9ac30c2e24dc1ee0ebd1e46d8bd286dccca7c819 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Geo.cpp,v 1.19 2008-03-20 11:44:05 geuzaine Exp $
+// $Id: GModelIO_Geo.cpp,v 1.20 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -32,23 +32,23 @@
 #include "Parser.h" // for Symbol_T
 #include "Field.h"
 
-int GModel::readGEO(const std::string &name)
-{
-  ParseFile((char*)name.c_str(), 1);
-  return importGEOInternals();
-}
-
-void GModel::createGEOInternals()
+void GModel::_createGEOInternals()
 {
   _geo_internals = new GEO_Internals;
 }
 
-void GModel::deleteGEOInternals()
+void GModel::_deleteGEOInternals()
 {
   delete _geo_internals;
   _geo_internals = 0;
 }
 
+int GModel::readGEO(const std::string &name)
+{
+  ParseFile((char*)name.c_str(), 1);
+  return importGEOInternals();
+}
+
 int GModel::importGEOInternals()
 {
   if(Tree_Nbr(_geo_internals->Points)) {
diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index 6dd5463f6d719fa527ccd3560cb378c428883143..00086493b833ff855607a5cbc840f944271fb84c 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_MED.cpp,v 1.10 2008-03-21 17:09:06 geuzaine Exp $
+// $Id: GModelIO_MED.cpp,v 1.11 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -19,373 +19,247 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+#include <string>
 #include "GModel.h"
-#include "MElement.h"
-#include "MVertex.h"
-#include "MEdge.h"
 #include "Message.h"
 
 #if defined(HAVE_MED)
 
 #include <map>
-#include <string>
-#include <vector>
-#include <iostream>
 #include <sstream>
-#include "GModelIO_MED.h"
+#include <vector>
+#include "MElement.h"
+#include "MVertex.h"
 
 extern "C" {
 #include <med.h>
 }
 
-ConversionData Data::MyConversionData;
-
-typedef std::list<int>::const_iterator listIter;
-
-// _____________________________________________ //
-//                                               //
-// Implementation  de la classe  ConversionData  //
-//  ____________________________________________ //
-
-ConversionData::ConversionData()
-{ 
-        // Correspondance des types GMSH et MEd
-        static const med_geometrie_element ValuesTypesOfElts[] = {  MED_SEG2,    MED_TRIA3,  MED_QUAD4,   MED_TETRA4,
-                                                                    MED_HEXA8,   MED_PENTA6, MED_PYRA5,   MED_SEG3, 
-                                                                    MED_TRIA6,   MED_QUAD8,  MED_TETRA10, MED_HEXA20, 
-                                                                    MED_PENTA15, MED_PYRA13, MED_POINT1,  MED_QUAD8,
-                                                                    MED_HEXA20,  MED_PENTA15, MED_PYRA13,  MED_NONE};
-                
-        int i=0;
-        while (ValuesTypesOfElts[i] != MED_NONE)
-        {
-           typesOfElts[i+1]=ValuesTypesOfElts[i];
-           ++i;
-        }
-
-        // **************************
-        // Numérotation des Noeuds
-        // **************************
-        //
-        // line ( 1er et 2nd Ordre)
-        i=1;
-        while ( i!=3 ) { 
-           medVertexOrder[1].push_back(i);
-           medVertexOrder[8].push_back(i);
-           ++i;
-        }
-        medVertexOrder[8].push_back(3);
-
-        // tria
-        i=1;
-        while ( i!=4 ) { 
-           medVertexOrder[2].push_back(i);
-           medVertexOrder[9].push_back(i);
-           ++i;
-        }
-        medVertexOrder[9].push_back(4);
-        medVertexOrder[9].push_back(5);
-        medVertexOrder[9].push_back(6);
-
-        // quad
-        i=1;
-        while ( i!=5 ) { 
-           medVertexOrder[3].push_back(i);
-           medVertexOrder[10].push_back(i);
-           medVertexOrder[16].push_back(i);
-           ++i;
-        }
-        medVertexOrder[10].push_back(5);
-        medVertexOrder[16].push_back(5);
-        medVertexOrder[10].push_back(6);
-        medVertexOrder[16].push_back(6);
-        medVertexOrder[10].push_back(7);
-        medVertexOrder[16].push_back(7);
-        medVertexOrder[10].push_back(8);
-        medVertexOrder[16].push_back(8);
-
-        // tetra
-        static const int OrderforTetra[] = { 1, 3, 2, 4, 0 };
-        i=0;
-        while ( OrderforTetra[i] != 0 ) { 
-           medVertexOrder[4].push_back(OrderforTetra[i]);
-           medVertexOrder[11].push_back(OrderforTetra[i]);
-           ++i;
-        }
-        static const int OrderforTetra2[] = { 7, 6 , 5, 8 , 9 , 10, 0 }; 
-        i=0;
-        while ( OrderforTetra2[i] != 0 ) { 
-           medVertexOrder[11].push_back(OrderforTetra2[i]);
-           ++i;
-        }
-
-        // hexa
-        // le 1 MED est le 1er Gmsh,   le 2nd Med est le 4ieme Gmsh , le 3 le 6 ...
-        static const int OrderforHexa[] = { 1, 5, 8,  4, 2, 6, 7, 3, 0};
-        i=0;
-        while (OrderforHexa[i] != 0)
-        {
-           medVertexOrder[5].push_back(OrderforHexa[i]);
-           medVertexOrder[12].push_back(OrderforHexa[i]);
-           medVertexOrder[17].push_back(OrderforHexa[i]);
-           ++i;
-        };
-        static const int OrderforHexa2[] = { 11,18,16,10,13,19,15,16,9,17,20,14, 0 }; 
-        i=0;
-        while ( OrderforHexa2[i] != 0 ) { 
-           medVertexOrder[12].push_back(OrderforHexa2[i]);
-           medVertexOrder[17].push_back(OrderforHexa2[i]);
-           ++i;
-        }
-
-        // Prism
-        static const int OrderforPrism[] = { 1, 3, 2, 4, 6, 5, 0};
-        i=0;
-        while (OrderforPrism[i] != 0)
-        {
-           medVertexOrder[6].push_back(OrderforPrism[i]);
-           medVertexOrder[13].push_back(OrderforPrism[i]);
-           medVertexOrder[18].push_back(OrderforPrism[i]);
-           ++i;
-        };
-        static const int OrderforPrism2[] = { 8, 10, 7, 14, 13, 15, 9, 12, 11, 0};
-        i=0;
-        while ( OrderforPrism2[i] != 0 ) { 
-           medVertexOrder[13].push_back(OrderforPrism2[i]);
-           medVertexOrder[18].push_back(OrderforPrism2[i]);
-           ++i;
-        }
-
-
-        //Pyra
-        static const int OrderforPyra[] = { 1, 4, 3, 2, 5, 0};
-        i=0;
-        while (OrderforPyra[i] != 0)
-        {
-           medVertexOrder[7].push_back(OrderforPyra[i]);
-           medVertexOrder[14].push_back(OrderforPyra[i]);
-           medVertexOrder[19].push_back(OrderforPyra[i]);
-           ++i;
-        };
-        static const int OrderforPyra2[] = { 7, 11 ,9 , 6, 8, 13, 12 , 10, 0};
-        i=0;
-        while ( OrderforPyra2[i] != 0 ) { 
-           medVertexOrder[14].push_back(OrderforPyra2[i]);
-           medVertexOrder[19].push_back(OrderforPyra2[i]);
-           ++i;
-        }
-};
-
-
-// ____________________________________ //
-//                                      //
-// Implementation  de la classe  MedIO  //
-//  ___________________________________ //
-
-MedIO::MedIO() : _numOfNode(1), _boolOpen(0), _meshName("monMaillage")
-{ 
-};
-
-// --------------------------------------------
-int MedIO::SetFile(const std::string &FileName)
-// --------------------------------------------
-// Open Med File
-// Headers Creation in Med File 
-// Mesh Creation in Med File
+static int getTypeForMED(int msh, med_geometrie_element &med)
 {
-   _FileName = FileName;
-   // All Meshes are 3D and unstructured Meshes
-   med_err ret = 0;
-   char des[MED_TAILLE_DESC+1]="Med file generated by  Gmsh";
-
-   _fid = MEDouvrir((char *) _FileName.c_str(),MED_CREATION);
-   if (_fid < 0) {
-       Msg(GERROR, "Unable to open file '%s'", _FileName.c_str());
-       return 0;
-   }
-
-   if (MEDfichDesEcr(_fid,des) < 0) {
-       Msg(GERROR, "Unable to write descriptor in file '%s'", _FileName.c_str());
-       return 0;
-   }
-
-   std::string description = "gmsh";
-   if (MEDmaaCr(_fid,(char *) _meshName.c_str(),3,MED_NON_STRUCTURE,(char *) description.c_str()) < 0) 
-   {
-       Msg(GERROR, "Error in Mesh Creation '%s'", _FileName.c_str());
-       return 0;
-   }
-   _boolOpen=1;
-    Msg(INFO, "File Open");
-   return 1;
+  switch(msh) {
+  case MSH_LIN_2: med = MED_SEG2; return 2; 
+  case MSH_TRI_3: med = MED_TRIA3; return 3; 
+  case MSH_QUA_4: med = MED_QUAD4; return 4; 
+  case MSH_TET_4: med = MED_TETRA4; return 4; 
+  case MSH_HEX_8: med = MED_HEXA8; return 8; 
+  case MSH_PRI_6: med = MED_PENTA6; return 6; 
+  case MSH_PYR_5: med = MED_PYRA5; return 5; 
+  case MSH_LIN_3: med = MED_SEG3; return 3; 
+  case MSH_TRI_6: med = MED_TRIA6; return 6; 
+  case MSH_TET_10: med = MED_TETRA10; return 10;
+  case MSH_PNT: med = MED_POINT1; return 1; 
+  case MSH_QUA_8: med = MED_QUAD8; return 8; 
+  case MSH_HEX_20: med = MED_HEXA20; return 20;
+  case MSH_PRI_15: med = MED_PENTA15; return 15;
+  case MSH_PYR_13: med = MED_PYRA13; return 13;
+  default: med = MED_NONE; return 0; 
+  }
 }
- 
-// -------------------------------------------------
-int MedIO::AddNode(MVertex* v, const int famille)
-// -------------------------------------------------
-{ 
-    //Msg(INFO, "Creation %d", v->getNum());
-   
-   if ( elements.find(v->getNum()) != elements.end() ) return 1;
-   coordonnees.push_back((med_float)v->x());
-   coordonnees.push_back((med_float)v->y());
-   coordonnees.push_back((med_float)v->z());
 
-   numOpt.push_back(v->getNum());
-   elements[v->getNum()] = _numOfNode ;
-   _numOfNode = _numOfNode + 1;
-
-   families.push_back(famille);
-   numFamilles.insert(famille);
-   return 1;
+template<class T>
+static void fillElementsMED(med_int family, std::vector<T*> &elements, med_int &ele, 
+			    std::vector<med_int> &num, std::vector<med_int> &conn,
+			    std::vector<med_int> &fam, med_geometrie_element &type)
+{
+  for(unsigned int i = 0; i < elements.size(); i++){
+    num.push_back(++ele);
+    for(int j = 0; j < elements[i]->getNumVertices(); j++)
+      conn.push_back(elements[i]->getVertexMED(j)->getNum());
+    fam.push_back(family);
+    if(!i) getTypeForMED(elements[i]->getTypeForMSH(), type);
+  }
 }
 
-// -------------------------------
-int MedIO::CreateFamilles( )
-// -------------------------------
+static void writeElementsMED(med_idt &fid, char *meshName,
+			     std::vector<med_int> &num, std::vector<med_int> &conn, 
+			     std::vector<med_int> &fam, med_geometrie_element type)
 {
-   numFamilles.insert(0);
-   std::set<int>::const_iterator itFam;
-   for (itFam = numFamilles.begin(); itFam != numFamilles.end(); ++itFam) {
-        med_err CR;
-        if (*itFam != 0 )
-        { std::ostringstream oss;
-          oss << *itFam;
-          std::string fam = "F_" + oss.str();
-          std::string group = "G_" + oss.str();
-          while (group.size() < 80) group = group + " ";
-          CR = MEDfamCr (_fid, (char *) _meshName.c_str(),(char *)fam.c_str(),*itFam, 0,0,0,0,(char *)group.c_str(),1);
-          CR=0;
-        }
-        else
-        {
-          std::string fam = "Famille0";
-          CR = MEDfamCr (_fid, (char *) _meshName.c_str(),(char *)fam.c_str(),*itFam, 0,0,0,0,0,0);
-        }
-        if ( CR < 0 )
-        {
-            Msg(GERROR, "Error in Family Creation '%d'", *itFam );
-            return 0;
-         }
-   }
-   return 1;
+  if(num.empty()) return;
+  if(MEDelementsEcr(fid, meshName, (med_int)3, &conn[0], MED_FULL_INTERLACE,
+		    0, MED_FAUX, 0, MED_FAUX, &fam[0], (med_int)num.size(),
+		    MED_MAILLE, type, MED_NOD) < 0)
+    Msg(GERROR, "Could not write elements");
 }
 
-
-// -------------------------------
-int MedIO::CreateElemt()
-// -------------------------------
+int GModel::writeMED(const std::string &name, double scalingFactor)
 {
-    std::map<int,med_geometrie_element>::const_iterator eltIter = Data::MyConversionData.typesOfElts.begin();
-    for ( ; eltIter != Data::MyConversionData.typesOfElts.end(); ++eltIter ) {
-       med_geometrie_element typemed =(*eltIter).second;
-       if (typemed == MED_POINT1) continue;
-       int nbNoeudElt = typemed % 100 ;
-       int nbElements = LesConn[typemed].size() / nbNoeudElt;
-       if (nbElements != 0 )
-           med_err CR = MEDelementsEcr (_fid, (char*) _meshName.c_str(),(med_int) 3, 
-                         &LesConn[typemed][0], MED_FULL_INTERLACE,
-                         NULL, MED_FAUX, NULL, MED_FAUX,
-                         &famElts[typemed][0],nbElements,
-                           MED_MAILLE,typemed,MED_NOD);
-
+  med_idt fid = MEDouvrir((char*)name.c_str(), MED_CREATION);
+  if(fid < 0) {
+    Msg(GERROR, "Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+
+  // write header
+  char des[MED_TAILLE_DESC + 1] = "MED file generated by Gmsh";
+  if(MEDfichDesEcr(fid, des) < 0) {
+    Msg(GERROR, "Unable to write MED descriptor");
+    return 0;
+  }
+
+  char *meshName = (char*)getName().c_str();
+  
+  // create 3D unstructured mesh
+  if(MEDmaaCr(fid, meshName, 3, MED_NON_STRUCTURE, (char*)"gmsh") < 0){
+    Msg(GERROR, "Could not create MED mesh");
+    return 0;
+  }
+  
+  // renumber all the vertices in a continuous sequence (MED
+  // connectivity is given in terms of vertex indices)
+  renumberMeshVertices(true);
+
+  // fill a vector containing all the geometrical entities in the model
+  std::vector<GEntity*> entities;
+  entities.insert(entities.end(), vertices.begin(), vertices.end());
+  entities.insert(entities.end(), edges.begin(), edges.end());
+  entities.insert(entities.end(), faces.begin(), faces.end());
+  entities.insert(entities.end(), regions.begin(), regions.end());
+
+  std::map<GEntity*, int> families;
+
+  // write the families
+  {
+    // always create a "0" family, with no groups or attributes
+    if(MEDfamCr(fid, meshName, (char*)"F_0", 0, 0, 0, 0, 0, 0, 0) < 0)
+      Msg(GERROR, "Could not create MED family 0");
+
+    // create one family per elementary entity, with one group per
+    // physical entity and no attributes
+    for(unsigned int i = 0; i < entities.size(); i++){
+      int num = - (families.size() + 1);
+      families[entities[i]] = num;
+      std::ostringstream fs;
+      fs << entities[i]->dim() << "D_" << entities[i]->tag();
+      std::string familyName = "F_" + fs.str();
+      std::string groupName;
+      for(unsigned j = 0; j < entities[i]->physicals.size(); j++){
+	std::string tmp = getPhysicalName(entities[i]->physicals[j]);
+	std::ostringstream gs;
+	gs << entities[i]->dim() << "D_" << tmp;
+	if(tmp.empty()) gs << entities[i]->physicals[j];
+	groupName += "G_" + gs.str();
+	groupName.resize((j + 1) * 80, ' ');
+      }
+      if(MEDfamCr(fid, meshName, (char*)familyName.c_str(), 
+		  (med_int)num, 0, 0, 0, 0, (char*)groupName.c_str(),
+		  (med_int)entities[i]->physicals.size()) < 0)
+	Msg(GERROR, "Could not create MED family %d", num);
     }
-};
-// ------------------------------------------
-int MedIO::Ecrit()
-// ------------------------------------------
-{
-    if (_boolOpen != 1)
-    {
-        Msg(GERROR, "File not Open");
-        return 0;
+  }
+  
+  // write the nodes
+  {
+    std::vector<med_float> coord;
+    std::vector<med_int> num, fam;
+    for(unsigned int i = 0; i < entities.size(); i++){
+      for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+	MVertex *v = entities[i]->mesh_vertices[j];
+	if(v->getNum() >= 0){
+	  coord.push_back(v->x() * scalingFactor);
+	  coord.push_back(v->y() * scalingFactor);
+	  coord.push_back(v->z() * scalingFactor);
+	  num.push_back(v->getNum());
+	  fam.push_back(0); // we never create node families
+	}
+      }
     }
-
-    int nbNoeuds=coordonnees.size() / 3;
-    if (nbNoeuds != families.size())
-    {
-        Msg(GERROR, "bad Vectors");
-        return 0;
+    char coordName[3 * MED_TAILLE_PNOM + 1] = 
+      "x               y               z               ";
+    char coordUnit[3 * MED_TAILLE_PNOM + 1] = 
+      "inconnu         inconnu         inconnu         ";
+    if(MEDnoeudsEcr(fid, meshName, (med_int)3, &coord[0], MED_FULL_INTERLACE, 
+		    MED_CART, coordName, coordUnit, 0, MED_FAUX, &num[0], 
+		    MED_VRAI, &fam[0], (med_int)num.size()) < 0)
+      Msg(GERROR, "Could not write nodes");
+  }
+  
+  // write the elements
+  {
+    med_geometrie_element typ;
+    med_int ele = 0;
+
+    { // points
+      std::vector<med_int> num, conn, fam;
+      for(viter it = firstVertex(); it != lastVertex(); it++){
+	num.push_back(++ele);
+	conn.push_back((*it)->mesh_vertices[0]->getNum());
+	fam.push_back(families[*it]);
+      }
+      writeElementsMED(fid, meshName, num, conn, fam, MED_POINT1);
     }
-
-    // *********************
-    // Creation des Familles
-    // *********************
-    int CRFam = CreateFamilles();
-    if ( CRFam < 0 )
-    {
-       Msg(GERROR, "Error in Nodes Families Creation ");
-       return 0;
+    { // lines
+      std::vector<med_int> num, conn, fam;
+      for(eiter it = firstEdge(); it != lastEdge(); it++)
+	fillElementsMED(families[*it], (*it)->lines, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
     }
-
-    // *************************
-    // Creation des Coordonnees
-    // *************************
-    char nomcoo[3*MED_TAILLE_PNOM+1] = "x               y               z               ";
-    char unicoo[3*MED_TAILLE_PNOM+1] = "inconnu         inconnu         inconnu         ";
-
-    med_err CR = MEDnoeudsEcr(_fid, (char*) _meshName.c_str(),(med_int) 3, 
-                              &coordonnees[0], MED_FULL_INTERLACE, MED_CART,
-                              nomcoo,unicoo, NULL, MED_FAUX, 
-                              &numOpt[0], MED_VRAI, 
-                              &families[0], nbNoeuds);
-    Msg(INFO, "%d ", CR);
-    if ( CR < 0 )
-    {
-       Msg(GERROR, "Error in Nodes Creation ");
-       return 0;
+    { // triangles
+      std::vector<med_int> num, conn, fam;
+      for(fiter it = firstFace(); it != lastFace(); it++)
+	fillElementsMED(families[*it], (*it)->triangles, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
     }
-    // ***************************
-    // Creation des connectivités
-    // ***************************
-    int CRElt = CreateElemt();
-    if ( CRElt < 0 )
-    {
-       Msg(GERROR, "Error in Elements Creation ");
-       return 0;
+    { // quads
+      std::vector<med_int> num, conn, fam;
+      for(fiter it = firstFace(); it != lastFace(); it++)
+	fillElementsMED(families[*it], (*it)->quadrangles, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
     }
-    return 1;
-}
-
-//--------------------
-int MedIO::CloseFile()
-//--------------------
-{
-   if (MEDfermer(_fid) < 0) 
-   {
-       Msg(GERROR, "Unable to close file '%s'",(char * ) _FileName.c_str());
-       return 0;
-   }
+    { // tets
+      std::vector<med_int> num, conn, fam;
+      for(riter it = firstRegion(); it != lastRegion(); it++)
+	fillElementsMED(families[*it], (*it)->tetrahedra, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
+    }
+    { // hexas
+      std::vector<med_int> num, conn, fam;
+      for(riter it = firstRegion(); it != lastRegion(); it++)
+	fillElementsMED(families[*it], (*it)->hexahedra, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
+    }
+    { // prisms
+      std::vector<med_int> num, conn, fam;
+      for(riter it = firstRegion(); it != lastRegion(); it++)
+	fillElementsMED(families[*it], (*it)->prisms, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
+    }
+    { // pyramids
+      std::vector<med_int> num, conn, fam;
+      for(riter it = firstRegion(); it != lastRegion(); it++)
+	fillElementsMED(families[*it], (*it)->pyramids, ele, num, conn, fam, typ);
+      writeElementsMED(fid, meshName, num, conn, fam, typ);
+    }
+  }
+  
+  if(MEDfermer(fid) < 0){
+    Msg(GERROR, "Unable to close file '%s'", (char*)name.c_str());
+    return 0;
+  }
+  
   return 1;
 }
 
-
-int GModel::writeMED(const std::string &name)
+int GModel::readMED(const std::string &name)
 {
 
-   MedIO MedDriver = MedIO();
-   int CR1 = MedDriver.SetFile(name);
-
-   renumberMeshVertices(noPhysicalGroups());
-   MedDriver.TraiteMed(vertices);
-   MedDriver.TraiteMed(edges);
-   MedDriver.TraiteMed(faces);
-   MedDriver.TraiteMed(regions);
-
-   int CR2 = MedDriver.Ecrit();
-   int CR3 = MedDriver.CloseFile();
-
-   return CR1 * CR2 * CR3 ;
+  return 1;
 }
 
 #else
 
-int GModel::writeMED(const std::string &name)
+int GModel::writeMED(const std::string &name, double scalingFactor)
 {
   Msg(GERROR, "Gmsh has to be compiled with MED support to write '%s'",
       name.c_str());
   return 0;
 }
 
-#endif                // du HAVE_LIBMED
+int GModel::readMED(const std::string &name)
+{
+  Msg(GERROR, "Gmsh has to be compiled with MED support to read '%s'",
+      name.c_str());
+  return 0;
+}
 
+#endif
diff --git a/Geo/GModelIO_MED.h b/Geo/GModelIO_MED.h
deleted file mode 100644
index d840c50bf8ffcdfc0e8e44a5a2fef0903fbffc56..0000000000000000000000000000000000000000
--- a/Geo/GModelIO_MED.h
+++ /dev/null
@@ -1,245 +0,0 @@
-#ifndef _GMODELIO_MED_H_
-#define _GMODELIO_MED_H_
-
-#include <string>
-#include <map>
-#include <vector>
-#include <set>
-#include <typeinfo>
-
-#include "GEntity.h"
-#include "GRegion.h"
-#include "GEdge.h"
-#include "GFace.h"
-#include "Message.h"
-
-#if defined(HAVE_MED)
-
-extern "C"
-{
-           #include "med.h"
-}
-
-typedef std::map<med_geometrie_element,std::vector<int> >  connectivities;
-typedef std::set<int>                  setFamille;
-
-class MedIO ;
-
-// _________________________________________ //
-//                                           //
-// Declaration  de la classe  ConversionData // 
-// (implementation ds le cpp)                //
-// ________________________________________  //
-//
-class ConversionData 
-{
-        public :
-
-           ConversionData();
-
-           std::map<int,med_geometrie_element> typesOfElts;
-           std::map<int,std::list<int> > medVertexOrder;
-
-           std::map<int,int> familleParDimension;
-           std::map<int,int> famillefamille;
-
-};
-
-// _________________________________________________ //
-//                                                   //
-// Declaration  de la classe  TraiteMailledeBase     // 
-// (implementations en fin de fichier)               //
-// ________________________________________________  //
-
-class Data 
-{
-    public :
-    static ConversionData MyConversionData ;
-};
-
-
-template<class T> 
-class TraiteMailledeBase : Data
-{
-  protected:
-
-  static inline int RecupereFamille(const T &ele, const int dimension, MedIO & monDriver);
-  template<class G> static void RecupereElement(const std::vector<G*> &ele, const int famille, MedIO& monDriver);
-
-};
-
-
-// _________________________________________________ //
-//                                                   //
-// Specialisations de la classe  TraiteMailledeBase  // 
-// selon le type de maille                           //
-// ________________________________________________  //
-
-template <class T>
-struct TraiteMaille : public TraiteMailledeBase<T> {
-  static inline  void AddElement(const T& ele, MedIO & monDriver)
-  { Msg(GERROR, "Wrong Call : basic Method AddElt"); };
-};
-
-
-template <>
-struct TraiteMaille<GVertex> : public TraiteMailledeBase<GVertex> {
-  static inline  void AddElement(const GVertex & ele, MedIO& monDriver) ;
-  // l'implementation est  en fin de fichier pour pouvoir utiliser la methode AddNode de MedIO
-};
-
-template <>
-struct TraiteMaille<GEdge> : public TraiteMailledeBase<GEdge> {
-  static inline void AddElement(const GEdge & ele, MedIO& monDriver) 
-  {
-           RecupereElement (ele.lines,  RecupereFamille(ele, 1, monDriver), monDriver);
-  }
-};
-
-template <>
-struct TraiteMaille<GFace> : public TraiteMailledeBase<GFace> {
-  static inline void AddElement(const GFace & ele,  MedIO&  monDriver) 
-  {
-           RecupereElement (ele.triangles,   RecupereFamille(ele, 2, monDriver), monDriver);
-           RecupereElement (ele.quadrangles, RecupereFamille(ele, 2, monDriver), monDriver);
-  }
-};
-
-template <>
-struct TraiteMaille<GRegion> : public TraiteMailledeBase<GRegion>{
-  static inline void AddElement(const GRegion & ele, MedIO &  monDriver) 
-  {
-           RecupereElement (ele.tetrahedra, RecupereFamille(ele, 3, monDriver), monDriver );
-           RecupereElement (ele.hexahedra,  RecupereFamille(ele, 3, monDriver), monDriver);
-           RecupereElement (ele.prisms,     RecupereFamille(ele, 3, monDriver), monDriver);
-           RecupereElement (ele.pyramids,   RecupereFamille(ele, 3, monDriver), monDriver);
-     }
-};
-
-// ______________________________________ //
-//                                        //
-// Declaration  de la classe  MedIO       // 
-// (implementation ds le cpp)             //
-// _____________________________________  //
-
-class MVertex;
-class MEdge;
-class MedIO 
-{
-   public:
-      MedIO();
-      int SetFile   (const std::string& theFileName);
-      int AddNode   (MVertex* const v, const int famille );
-      int Ecrit     ();
-      int CloseFile ();
-
-   private :
-      int CreateFamilles();
-      int CreateElemt();
-
-   public :
-      connectivities LesConn;
-      setFamille     numFamilles;
-      std::vector<int>       families;
-      std::map<int,int>      elements;
-      std::map<med_geometrie_element,std::vector<int> > famElts;
-
-   private :
-
-      std::vector<med_float> coordonnees;
-      std::vector<med_int>   numOpt;
-  
-      int _numOfNode;
-      int _boolOpen;
-
-      med_idt _fid;
-      std::string  _FileName;
-      std::string  _meshName;
-
-   public : 
-     static ConversionData MyConversionData ;
-    
-      
-   template<class T>
-   void TraiteMed(const std::set<T*,GEntityLessThan> & seq) 
-   {
-        typename std::set<T*,GEntityLessThan>::const_iterator iter=seq.begin();
-
-        for ( ; iter != seq.end() ; ++iter)
-        {
-             const T& element= *(*iter);
-             TraiteMaille<T>::AddElement(element, *this);
-        }
-   }
-};
-
-// _______________________________________________________________ //
-//                                                                 //
-// Implementation des Methodes des classes de type TraiteElement   // 
-// qui utilisent MedIO                                             //
-// _______________________________________________________________ //
-
-void TraiteMaille<GVertex>::AddElement(const GVertex & ele, MedIO& monDriver)
-{
-     int famille = 0;
-     if (ele.physicals.size() != 0 ) famille = ele.physicals[0];
-     monDriver.AddNode(ele.mesh_vertices[0],famille);
-}
-
-//                                      *-*-*-*-*-*-*-*-*-*
-
-template<class T> int TraiteMailledeBase<T>::RecupereFamille(const T &ele, const int dimension, MedIO & monDriver)
-// 
-// Dans le format MED, Les familles de mailles sont negatives (contrairement au physical group de GMSH)
-// Les familles ne contiennent des mailles de la meme dimension 
-{
-     int famille = 0;
-     if (ele.physicals.size() != 0) famille=ele.physicals[0];
-     if (famille == 0) return 0; 
-
-     if (famille > 0)  famille = -1 * famille;
-
-     int familleInitiale=famille;
-     while (MyConversionData.familleParDimension[famille] != dimension)
-     {
-        if (monDriver.numFamilles.find(famille) == monDriver.numFamilles.end())
-        {   monDriver.numFamilles.insert(famille);
-            MyConversionData.familleParDimension[famille]=dimension;
-        }
-         else famille=famille-1000;
-     }
-
-     MyConversionData.famillefamille[familleInitiale]=famille;
-     return famille;
-}
-
-//                                      *-*-*-*-*-*-*-*-*-*
-
-template<class T>
-template<class G> void TraiteMailledeBase<T>::RecupereElement(const std::vector<G*> &ele, const int famille, MedIO& monDriver) 
-// l'implementation aurait pu etre laissee avec la declaration (ok pour GCC 3.3.5)
-// mais pour eviter des eventuels pb de portage il semble preferable de la separer
-{
-
-     typedef std::list<int>::const_iterator listIter;
-
-     if (ele.size() == 0) return;
-     const int type=ele[0]->getTypeForMSH();
-     const med_geometrie_element medType= Data::MyConversionData.typesOfElts[type];
-
-     for (unsigned int elt = 0; elt < ele.size(); ++elt) {
-        const int monNum=ele[elt]->getNum();
-        for (listIter monIter  = Data::MyConversionData.medVertexOrder[type].begin(); 
-                      monIter != Data::MyConversionData.medVertexOrder[type].end(); 
-                      ++monIter) {
-              const int NoeudATraiter = *monIter - 1;
-              const int Noeud = ele[elt]->getVertex(NoeudATraiter)->getNum();
-              monDriver.AddNode(ele[elt]->getVertex(NoeudATraiter),0);
-              monDriver.LesConn[medType].push_back(monDriver.elements[Noeud]);
-           }
-        monDriver.famElts[medType].push_back(famille);
-        }
-}
-
-#endif               // du HAVE_LIBMED
-#endif               // du define de _GMODELIO_MED_H_
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 9165e18780ecfc962e140d8c78ea9a53a02893bf..296373508b725d59b03c26e0d7b7384a2ab61ea9 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.45 2008-03-21 07:21:05 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.46 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -38,58 +38,6 @@
 #  include "Message.h"
 #endif
 
-template<class T>
-static void addElements(std::vector<T*> &dst, const std::vector<MElement*> &src)
-{
-  for(unsigned int i = 0; i < src.size(); i++) dst.push_back((T*)src[i]);
-}
-
-static void storeElementsInEntities(GModel *m, 
-                                    std::map<int, std::vector<MElement*> > &map)
-{
-  std::map<int, std::vector<MElement*> >::const_iterator it = map.begin();
-  for(; it != map.end(); ++it){
-    if(!it->second.size()) continue;
-    int numEdges = it->second[0]->getNumEdges();
-    switch(numEdges){
-    case 1: 
-      {
-        GEdge *e = m->getEdgeByTag(it->first);
-        if(!e){
-          e = new discreteEdge(m, it->first);
-          m->add(e);
-        }
-        addElements(e->lines, it->second);
-      }
-      break;
-    case 3: case 4: 
-      {
-        GFace *f = m->getFaceByTag(it->first);
-        if(!f){
-          f = new discreteFace(m, it->first);
-          m->add(f);
-        }
-        if(numEdges == 3) addElements(f->triangles, it->second);
-        else addElements(f->quadrangles, it->second);
-      }
-      break;
-    case 6: case 12: case 9: case 8:
-      {
-        GRegion *r = m->getRegionByTag(it->first);
-        if(!r){
-          r = new discreteRegion(m, it->first);
-          m->add(r);
-        }
-        if(numEdges == 6) addElements(r->tetrahedra, it->second);
-        else if(numEdges == 12) addElements(r->hexahedra, it->second);
-        else if(numEdges == 9) addElements(r->prisms, it->second);
-        else addElements(r->pyramids, it->second);
-      }
-      break;
-    }
-  }
-}
-
 static void storeVerticesInEntities(std::map<int, MVertex*> &vertices)
 {
   std::map<int, MVertex*>::const_iterator it = vertices.begin();
@@ -211,38 +159,31 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
                              std::map<int, std::vector<MElement*> > elem[7],
                              std::map<int, std::map<int, std::string> > physicals[4])
 {
-  int dim = 0;
-
-  switch (type) {
-  case MSH_PNT:    points[reg].push_back(v[0]); dim = 0; break;
-  case MSH_LIN_2:  elem[0][reg].push_back(new MLine(v, num, part)); dim = 1; break;
-  case MSH_LIN_3:  elem[0][reg].push_back(new MLine3(v, num, part)); dim = 1; break;
-  case MSH_LIN_4:  elem[0][reg].push_back(new MLineN(v, num, part)); dim = 1; break;
-  case MSH_LIN_5:  elem[0][reg].push_back(new MLineN(v, num, part)); dim = 1; break;
-  case MSH_LIN_6:  elem[0][reg].push_back(new MLineN(v, num, part)); dim = 1; break;
-  case MSH_TRI_3:  elem[1][reg].push_back(new MTriangle(v, num, part)); dim = 2; break;
-  case MSH_TRI_6:  elem[1][reg].push_back(new MTriangle6(v, num, part)); dim = 2; break;
-  case MSH_TRI_9:  elem[1][reg].push_back(new MTriangleN(v, 3, num, part)); dim = 2; break;
-  case MSH_TRI_10: elem[1][reg].push_back(new MTriangleN(v, 3, num, part)); dim = 2; break;
-  case MSH_TRI_12: elem[1][reg].push_back(new MTriangleN(v, 4, num, part)); dim = 2; break;
-  case MSH_TRI_15: elem[1][reg].push_back(new MTriangleN(v, 4, num, part)); dim = 2; break;
-  case MSH_TRI_15I:elem[1][reg].push_back(new MTriangleN(v, 5, num, part)); dim = 2; break;
-  case MSH_TRI_21: elem[1][reg].push_back(new MTriangleN(v, 5, num, part)); dim = 2; break;
-  case MSH_QUA_4:  elem[2][reg].push_back(new MQuadrangle(v, num, part)); dim = 2; break;
-  case MSH_QUA_8:  elem[2][reg].push_back(new MQuadrangle8(v, num, part)); dim = 2; break;
-  case MSH_QUA_9:  elem[2][reg].push_back(new MQuadrangle9(v, num, part)); dim = 2; break;
-  case MSH_TET_4:  elem[3][reg].push_back(new MTetrahedron(v, num, part)); dim = 3; break;
-  case MSH_TET_10: elem[3][reg].push_back(new MTetrahedron10(v, num, part)); dim = 3; break;
-  case MSH_HEX_8:  elem[4][reg].push_back(new MHexahedron(v, num, part)); dim = 3; break;
-  case MSH_HEX_20: elem[4][reg].push_back(new MHexahedron20(v, num, part)); dim = 3; break;
-  case MSH_HEX_27: elem[4][reg].push_back(new MHexahedron27(v, num, part)); dim = 3; break;
-  case MSH_PRI_6:  elem[5][reg].push_back(new MPrism(v, num, part)); dim = 3; break;
-  case MSH_PRI_15: elem[5][reg].push_back(new MPrism15(v, num, part)); dim = 3; break;
-  case MSH_PRI_18: elem[5][reg].push_back(new MPrism18(v, num, part)); dim = 3; break;
-  case MSH_PYR_5:  elem[6][reg].push_back(new MPyramid(v, num, part)); dim = 3; break;
-  case MSH_PYR_13: elem[6][reg].push_back(new MPyramid13(v, num, part)); dim = 3; break;
-  case MSH_PYR_14: elem[6][reg].push_back(new MPyramid14(v, num, part)); dim = 3; break;
-  default: Msg(GERROR, "Unknown type (%d) for element %d", type, num); break;
+  int dim;
+  if(type == MSH_POINT){
+    dim = 0;
+    points[reg].push_back(v[0]);
+  }
+  else{
+    MElementFactory factory;
+    MElement *e = factory.create(type, v, num, part);
+    if(!e){
+      Msg(GERROR, "Unknown type of element %d", type);
+      return;
+    }
+    dim = e->getDim();
+    int idx;
+    switch(e->getNumEdges()){
+    case 1 : idx = 0; break;
+    case 3 : idx = 1; break;
+    case 4 : idx = 2; break;
+    case 6 : idx = 3; break;
+    case 12 : idx = 4; break;
+    case 9 : idx = 5; break;
+    case 8 : idx = 6; break;
+    default : Msg(GERROR, "Wrong number of edges in element"); return;
+    }
+    elem[idx][reg].push_back(e);
   }
   
   if(physical && (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
@@ -474,7 +415,7 @@ int GModel::readMSH(const std::string &name)
   bool noElements = true;
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++){
     noElements &= elements[i].empty();
-    storeElementsInEntities(this, elements[i]);
+    _storeElementsInEntities(elements[i]);
   }
 
   // special case: if there are no elements, create one geometry
@@ -503,7 +444,7 @@ int GModel::readMSH(const std::string &name)
   }
 
   // associate the correct geometrical entity with each mesh vertex
-  associateEntityWithMeshVertices();
+  _associateEntityWithMeshVertices();
 
   // special case for geometry vertices: now that the correct
   // geometrical entity has been associated with the vertices, we
@@ -1171,8 +1112,8 @@ int GModel::readVRML(const std::string &name)
   }
 
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
-    storeElementsInEntities(this, elements[i]);
-  associateEntityWithMeshVertices();
+    _storeElementsInEntities(elements[i]);
+  _associateEntityWithMeshVertices();
   storeVerticesInEntities(allVertexVector);
 
   fclose(fp);
@@ -1374,8 +1315,8 @@ int GModel::readUNV(const std::string &name)
   }
   
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
-    storeElementsInEntities(this, elements[i]);
-  associateEntityWithMeshVertices();
+    _storeElementsInEntities(elements[i]);
+  _associateEntityWithMeshVertices();
   storeVerticesInEntities(vertexMap);
   for(int i = 0; i < 4; i++)  
     storePhysicalTagsInEntities(this, i, physicals[i]);
@@ -1602,8 +1543,8 @@ int GModel::readMESH(const std::string &name)
   }
 
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
-    storeElementsInEntities(this, elements[i]);
-  associateEntityWithMeshVertices();
+    _storeElementsInEntities(elements[i]);
+  _associateEntityWithMeshVertices();
   storeVerticesInEntities(vertexVector);
 
   fclose(fp);
@@ -1946,8 +1887,8 @@ int GModel::readBDF(const std::string &name)
   }
   
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
-    storeElementsInEntities(this, elements[i]);
-  associateEntityWithMeshVertices();
+    _storeElementsInEntities(elements[i]);
+  _associateEntityWithMeshVertices();
   storeVerticesInEntities(vertexMap);
 
   fclose(fp);
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index a503a4aafd94e2366d5c96f0e4b023e0baebaf0c..0fcd376fd811ef2e4901d0edb96cf7dd8fd0ff51 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.30 2008-03-20 11:44:05 geuzaine Exp $
+// $Id: GModelIO_OCC.cpp,v 1.31 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -431,6 +431,12 @@ void OCC_Internals::buildGModel(GModel *model)
   }
 }
 
+void GModel::_deleteOCCInternals()
+{
+  if(_occ_internals) delete _occ_internals;
+  _occ_internals = 0;
+}
+
 int GModel::readOCCSTEP(const std::string &fn)
 {
   _occ_internals = new OCC_Internals;
@@ -459,12 +465,6 @@ int GModel::readOCCBREP(const std::string &fn)
   return 1;
 }
 
-void GModel::deleteOCCInternals()
-{
-  if(_occ_internals) delete _occ_internals;
-  _occ_internals = 0;
-}
-
 /*
   OCC Creation routines
 */
@@ -584,6 +584,10 @@ void OCC_Internals::Sphere(const SPoint3 & center, const double & radius,
 
 #else
 
+void GModel::_deleteOCCInternals()
+{
+}
+
 int GModel::readOCCSTEP(const std::string &fn)
 {
   Msg(GERROR, "Gmsh has to be compiled with OpenCascade support to load '%s'",
@@ -605,8 +609,4 @@ int GModel::readOCCBREP(const std::string &fn)
   return 0;
 }
 
-void GModel::deleteOCCInternals()
-{
-}
-
 #endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index ebc5446bce0c49f4ec2bd4b9a08f9a5fd16da9c2..98d6d7b2793bc91172069329025adab1acae4c14 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.61 2008-03-20 11:44:06 geuzaine Exp $
+// $Id: MElement.cpp,v 1.62 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -764,3 +764,39 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
     z[0] = pnt1.z(); z[1] = pnt2.z();
   }
 }
+
+MElement *MElementFactory::create(int type, std::vector<MVertex*> &v, 
+				  int num, int part)
+{
+  switch (type) {
+  case MSH_PNT:    return 0;
+  case MSH_LIN_2:  return new MLine(v, num, part);
+  case MSH_LIN_3:  return new MLine3(v, num, part);
+  case MSH_LIN_4:  return new MLineN(v, num, part);
+  case MSH_LIN_5:  return new MLineN(v, num, part);
+  case MSH_LIN_6:  return new MLineN(v, num, part);
+  case MSH_TRI_3:  return new MTriangle(v, num, part);
+  case MSH_TRI_6:  return new MTriangle6(v, num, part);
+  case MSH_TRI_9:  return new MTriangleN(v, 3, num, part);
+  case MSH_TRI_10: return new MTriangleN(v, 3, num, part);
+  case MSH_TRI_12: return new MTriangleN(v, 4, num, part);
+  case MSH_TRI_15: return new MTriangleN(v, 4, num, part);
+  case MSH_TRI_15I:return new MTriangleN(v, 5, num, part);
+  case MSH_TRI_21: return new MTriangleN(v, 5, num, part);
+  case MSH_QUA_4:  return new MQuadrangle(v, num, part);
+  case MSH_QUA_8:  return new MQuadrangle8(v, num, part);
+  case MSH_QUA_9:  return new MQuadrangle9(v, num, part);
+  case MSH_TET_4:  return new MTetrahedron(v, num, part);
+  case MSH_TET_10: return new MTetrahedron10(v, num, part);
+  case MSH_HEX_8:  return new MHexahedron(v, num, part);
+  case MSH_HEX_20: return new MHexahedron20(v, num, part);
+  case MSH_HEX_27: return new MHexahedron27(v, num, part);
+  case MSH_PRI_6:  return new MPrism(v, num, part);
+  case MSH_PRI_15: return new MPrism15(v, num, part);
+  case MSH_PRI_18: return new MPrism18(v, num, part);
+  case MSH_PYR_5:  return new MPyramid(v, num, part);
+  case MSH_PYR_13: return new MPyramid13(v, num, part);
+  case MSH_PYR_14: return new MPyramid14(v, num, part);
+  default:         return 0;
+  }
+}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 69001261d0bbf17fbb4a0f21e4b6f3c9b3180d1e..f0896b2418c003d8eee730df1fcf765dc171d0eb 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -89,6 +89,9 @@ class MElement
   // get the vertex using the Nastran BDF ordering
   virtual MVertex *getVertexBDF(int num){ return getVertex(num); }
 
+  // get the vertex using MED ordering
+  virtual MVertex *getVertexMED(int num){ return getVertex(num); }
+
   // get the number of vertices associated with edges, faces and
   // volumes (nonzero only for higher order elements)
   virtual int getNumEdgeVertices(){ return 0; }
@@ -208,6 +211,11 @@ class MElementLessThanLexicographic{
   }
 };
 
+class MElementFactory{
+ public:
+  MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0);
+};
+
 class MLine : public MElement {
  protected:
   MVertex *_v[2];
@@ -372,6 +380,11 @@ class MTriangle : public MElement {
   virtual double gammaShapeMeasure();
   virtual int getNumVertices(){ return 3; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[3] = {0, 2, 1};
+    return getVertex(map[num]); 
+  }
   virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
   { 
     if(_v[0] != v1 && _v[0] != v2) return _v[0];
@@ -475,6 +488,11 @@ class MTriangle6 : public MTriangle {
     static const int map[6] = {0, 3, 1, 4, 2, 5};
     return getVertex(map[num]); 
   }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[6] = {0, 2, 1, 5, 4, 3};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 3; }
   virtual int getNumEdgesRep();
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
@@ -612,6 +630,11 @@ class MQuadrangle : public MElement {
   virtual int getDim(){ return 2; }
   virtual int getNumVertices(){ return 4; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[4] = {0, 3, 2, 1};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdges(){ return 4; }
   virtual MEdge getEdge(int num)
   {
@@ -702,6 +725,11 @@ class MQuadrangle8 : public MQuadrangle {
     static const int map[8] = {0, 4, 1, 5, 2, 6, 3, 7};
     return getVertex(map[num]); 
   }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[8] = {0, 3, 2, 1, 7, 6, 5, 4};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 4; }
   virtual int getNumEdgesRep(){ return 8; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -811,6 +839,11 @@ class MTetrahedron : public MElement {
   virtual int getDim(){ return 3; }
   virtual int getNumVertices(){ return 4; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[4] = {0, 2, 1, 3};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdges(){ return 6; }
   virtual MEdge getEdge(int num)
   {
@@ -935,6 +968,11 @@ class MTetrahedron10 : public MTetrahedron {
     static const int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8};
     return getVertex(map[num]); 
   }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[10] = {0, 2, 1, 3, 6, 5, 4, 7, 8, 9};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 6; }
   virtual int getNumEdgesRep(){ return 12; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/Makefile b/Geo/Makefile
index c2108d34c1bc8c9d3cbea85c1f83387c62f1515d..47687d6eb30f59eab7c1b4b2cfb7a9a41b70a2b4 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.196 2008-03-19 20:13:45 geuzaine Exp $
+# $Id: Makefile,v 1.197 2008-03-23 21:42:57 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -202,8 +202,7 @@ GModelIO_CGNS.o: GModelIO_CGNS.cpp GModel.h GVertex.h GEntity.h Range.h \
   MEdgeHash.h ../Common/Hash.h MFaceHash.h
 GModelIO_MED.o: GModelIO_MED.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
-  GFace.h GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
-  MVertex.h MEdge.h MFace.h ../Common/Message.h
+  GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h
 ExtrudeParams.o: ExtrudeParams.cpp ../Common/Message.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 7c8f41abb36e2613235fda640ac23bd936c5abd1..20e4085b7bd48352cee77703fbedce7a852bf92c 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1,4 +1,4 @@
-// $Id: Field.cpp,v 1.25 2008-03-21 09:55:43 geuzaine Exp $
+// $Id: Field.cpp,v 1.26 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -74,7 +74,7 @@ class FieldOptionDouble : public FieldOption {
 };
 
 class FieldOptionInt : public FieldOption {
-public:
+ public:
   int &val;
   FieldOptionType get_type()
   {
@@ -643,7 +643,8 @@ class ParametricField : public Field {
 
 class PostViewField : public Field {
   OctreePost *octree;
- public:int view_index;
+ public:
+  int view_index;
   double operator() (double x, double y, double z)
   {
     if(view_index < 0 || view_index >= (int)PView::list.size())
diff --git a/Parser/CreateFile.cpp b/Parser/CreateFile.cpp
index d53ff7cbfabcceed232fe49d07a012c05abd50b8..6891fb0d389553b58cfbf4468eda391ee8746105 100644
--- a/Parser/CreateFile.cpp
+++ b/Parser/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.26 2008-03-20 11:44:09 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.27 2008-03-23 21:42:57 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -186,7 +186,7 @@ void CreateOutputFile(const char *filename, int format)
     break;
 
   case FORMAT_MED:
-    GModel::current()->writeMED(name);
+    GModel::current()->writeMED(name, CTX.mesh.scaling_factor);
     break;
 
   case FORMAT_POS:
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index b852c7839bab4609e319922da538e34c72adee2a..6530b9d3f69cd117617c9a5899cc04b4daba5d78 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -324,7 +324,7 @@
 /* Copy the first part of user declarations.  */
 #line 1 "Gmsh.y"
 
-// $Id: Gmsh.tab.cpp,v 1.355 2008-03-21 09:39:39 geuzaine Exp $
+// $Id: Gmsh.tab.cpp,v 1.356 2008-03-23 21:42:58 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -993,28 +993,28 @@ static const yytype_uint16 yyrline[] =
     1498,  1514,  1534,  1551,  1568,  1589,  1594,  1599,  1604,  1609,
     1620,  1626,  1635,  1636,  1641,  1644,  1648,  1671,  1694,  1717,
     1745,  1754,  1758,  1771,  1787,  1802,  1816,  1822,  1828,  1837,
-    1851,  1899,  1915,  1929,  1948,  1958,  1980,  1984,  1989,  1994,
-    2006,  2023,  2040,  2067,  2094,  2125,  2133,  2139,  2146,  2150,
-    2159,  2167,  2175,  2184,  2183,  2196,  2195,  2208,  2207,  2220,
-    2219,  2231,  2230,  2246,  2253,  2260,  2267,  2274,  2281,  2288,
-    2295,  2302,  2310,  2309,  2321,  2320,  2332,  2331,  2343,  2342,
-    2354,  2353,  2365,  2364,  2376,  2375,  2387,  2386,  2398,  2397,
-    2412,  2415,  2421,  2430,  2450,  2473,  2477,  2501,  2519,  2537,
-    2555,  2584,  2619,  2624,  2651,  2665,  2678,  2695,  2701,  2707,
-    2710,  2719,  2729,  2730,  2731,  2732,  2733,  2734,  2735,  2736,
-    2737,  2744,  2745,  2746,  2747,  2748,  2749,  2750,  2751,  2752,
-    2753,  2754,  2755,  2756,  2757,  2758,  2759,  2760,  2761,  2762,
-    2763,  2764,  2765,  2766,  2767,  2768,  2769,  2770,  2771,  2772,
-    2773,  2774,  2775,  2777,  2778,  2779,  2780,  2781,  2782,  2783,
-    2784,  2785,  2786,  2787,  2788,  2789,  2790,  2791,  2792,  2793,
-    2794,  2795,  2796,  2797,  2806,  2807,  2808,  2809,  2810,  2811,
-    2812,  2816,  2832,  2847,  2867,  2880,  2893,  2916,  2934,  2952,
-    2970,  2988,  2996,  3000,  3004,  3008,  3012,  3019,  3023,  3027,
-    3031,  3038,  3043,  3051,  3056,  3060,  3065,  3069,  3077,  3088,
-    3096,  3104,  3110,  3121,  3141,  3151,  3161,  3178,  3205,  3210,
-    3214,  3218,  3231,  3235,  3247,  3254,  3275,  3279,  3294,  3299,
-    3306,  3310,  3317,  3321,  3329,  3337,  3351,  3365,  3369,  3388,
-    3411
+    1851,  1899,  1915,  1928,  1947,  1957,  1979,  1983,  1988,  1993,
+    2005,  2022,  2039,  2066,  2093,  2124,  2132,  2138,  2145,  2149,
+    2158,  2166,  2174,  2183,  2182,  2195,  2194,  2207,  2206,  2219,
+    2218,  2230,  2229,  2245,  2252,  2259,  2266,  2273,  2280,  2287,
+    2294,  2301,  2309,  2308,  2320,  2319,  2331,  2330,  2342,  2341,
+    2353,  2352,  2364,  2363,  2375,  2374,  2386,  2385,  2397,  2396,
+    2411,  2414,  2420,  2429,  2449,  2472,  2476,  2500,  2518,  2536,
+    2554,  2583,  2618,  2623,  2650,  2664,  2677,  2694,  2700,  2706,
+    2709,  2718,  2728,  2729,  2730,  2731,  2732,  2733,  2734,  2735,
+    2736,  2743,  2744,  2745,  2746,  2747,  2748,  2749,  2750,  2751,
+    2752,  2753,  2754,  2755,  2756,  2757,  2758,  2759,  2760,  2761,
+    2762,  2763,  2764,  2765,  2766,  2767,  2768,  2769,  2770,  2771,
+    2772,  2773,  2774,  2776,  2777,  2778,  2779,  2780,  2781,  2782,
+    2783,  2784,  2785,  2786,  2787,  2788,  2789,  2790,  2791,  2792,
+    2793,  2794,  2795,  2796,  2805,  2806,  2807,  2808,  2809,  2810,
+    2811,  2815,  2831,  2846,  2866,  2879,  2892,  2915,  2933,  2951,
+    2969,  2987,  2995,  2999,  3003,  3007,  3011,  3018,  3022,  3026,
+    3030,  3037,  3042,  3050,  3055,  3059,  3064,  3068,  3076,  3087,
+    3095,  3103,  3109,  3120,  3140,  3150,  3160,  3177,  3204,  3209,
+    3213,  3217,  3230,  3234,  3246,  3253,  3274,  3278,  3293,  3298,
+    3305,  3309,  3316,  3320,  3328,  3336,  3350,  3364,  3368,  3387,
+    3410
 };
 #endif
 
@@ -5707,9 +5707,8 @@ yyreduce:
     {
       if(!strcmp((yyvsp[(1) - (7)].c), "Background") && !strcmp((yyvsp[(2) - (7)].c), "Mesh")  && !strcmp((yyvsp[(3) - (7)].c), "View")){
 	int index = (int)(yyvsp[(5) - (7)].d);
-	if(index >= 0 && index < (int)PView::list.size()){
+	if(index >= 0 && index < (int)PView::list.size())
 	  GModel::current()->getFields()->set_background_mesh(index);
-	}
 	else
 	  yymsg(GERROR, "Unknown view %d", index);
       }
@@ -5720,7 +5719,7 @@ yyreduce:
     break;
 
   case 143:
-#line 1930 "Gmsh.y"
+#line 1929 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(1) - (3)].c), "Sleep")){
 	SleepInSeconds((yyvsp[(2) - (3)].d));
@@ -5742,7 +5741,7 @@ yyreduce:
     break;
 
   case 144:
-#line 1949 "Gmsh.y"
+#line 1948 "Gmsh.y"
     {
        try {
 	 GMSH_PluginManager::instance()->action((yyvsp[(3) - (7)].c), (yyvsp[(6) - (7)].c), 0);
@@ -5755,7 +5754,7 @@ yyreduce:
     break;
 
   case 145:
-#line 1959 "Gmsh.y"
+#line 1958 "Gmsh.y"
     {
       if(!strcmp((yyvsp[(2) - (3)].c), "ElementsFromAllViews"))
 	PView::combine(false, 1, CTX.post.combine_remove_orig);
@@ -5780,14 +5779,14 @@ yyreduce:
     break;
 
   case 146:
-#line 1981 "Gmsh.y"
+#line 1980 "Gmsh.y"
     {
       exit(0);
     ;}
     break;
 
   case 147:
-#line 1985 "Gmsh.y"
+#line 1984 "Gmsh.y"
     {
       CTX.forced_bbox = 0;
       SetBoundingBox();
@@ -5795,7 +5794,7 @@ yyreduce:
     break;
 
   case 148:
-#line 1990 "Gmsh.y"
+#line 1989 "Gmsh.y"
     {
       CTX.forced_bbox = 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));
@@ -5803,7 +5802,7 @@ yyreduce:
     break;
 
   case 149:
-#line 1995 "Gmsh.y"
+#line 1994 "Gmsh.y"
     {
 #if defined(HAVE_FLTK)
       Draw();
@@ -5812,7 +5811,7 @@ yyreduce:
     break;
 
   case 150:
-#line 2007 "Gmsh.y"
+#line 2006 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(3) - (6)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(5) - (6)].d);
@@ -5832,7 +5831,7 @@ yyreduce:
     break;
 
   case 151:
-#line 2024 "Gmsh.y"
+#line 2023 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(3) - (8)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(5) - (8)].d);
@@ -5852,7 +5851,7 @@ yyreduce:
     break;
 
   case 152:
-#line 2041 "Gmsh.y"
+#line 2040 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(5) - (8)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(7) - (8)].d);
@@ -5882,7 +5881,7 @@ yyreduce:
     break;
 
   case 153:
-#line 2068 "Gmsh.y"
+#line 2067 "Gmsh.y"
     {
       LoopControlVariablesTab[ImbricatedLoop][0] = (yyvsp[(5) - (10)].d);
       LoopControlVariablesTab[ImbricatedLoop][1] = (yyvsp[(7) - (10)].d);
@@ -5912,7 +5911,7 @@ yyreduce:
     break;
 
   case 154:
-#line 2095 "Gmsh.y"
+#line 2094 "Gmsh.y"
     {
       if(ImbricatedLoop <= 0){
 	yymsg(GERROR, "Invalid For/EndFor loop");
@@ -5946,7 +5945,7 @@ yyreduce:
     break;
 
   case 155:
-#line 2126 "Gmsh.y"
+#line 2125 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->createFunction((yyvsp[(2) - (2)].c), gmsh_yyin, gmsh_yyname,
 						      gmsh_yylineno))
@@ -5957,7 +5956,7 @@ yyreduce:
     break;
 
   case 156:
-#line 2134 "Gmsh.y"
+#line 2133 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->leaveFunction(&gmsh_yyin, gmsh_yyname,
 						     gmsh_yylineno))
@@ -5966,7 +5965,7 @@ yyreduce:
     break;
 
   case 157:
-#line 2140 "Gmsh.y"
+#line 2139 "Gmsh.y"
     {
       if(!FunctionManager::Instance()->enterFunction((yyvsp[(2) - (3)].c), &gmsh_yyin, gmsh_yyname,
 						     gmsh_yylineno))
@@ -5976,20 +5975,20 @@ yyreduce:
     break;
 
   case 158:
-#line 2147 "Gmsh.y"
+#line 2146 "Gmsh.y"
     {
       if(!(yyvsp[(3) - (4)].d)) skip_until("If", "EndIf");
     ;}
     break;
 
   case 159:
-#line 2151 "Gmsh.y"
+#line 2150 "Gmsh.y"
     {
     ;}
     break;
 
   case 160:
-#line 2160 "Gmsh.y"
+#line 2159 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, (yyvsp[(4) - (5)].l), 
@@ -6000,7 +5999,7 @@ yyreduce:
     break;
 
   case 161:
-#line 2168 "Gmsh.y"
+#line 2167 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, (yyvsp[(10) - (11)].l), 
@@ -6011,7 +6010,7 @@ yyreduce:
     break;
 
   case 162:
-#line 2176 "Gmsh.y"
+#line 2175 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[(12) - (13)].l), 
@@ -6022,14 +6021,14 @@ yyreduce:
     break;
 
   case 163:
-#line 2184 "Gmsh.y"
+#line 2183 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 164:
-#line 2188 "Gmsh.y"
+#line 2187 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE, (yyvsp[(4) - (7)].l), 
@@ -6040,14 +6039,14 @@ yyreduce:
     break;
 
   case 165:
-#line 2196 "Gmsh.y"
+#line 2195 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 166:
-#line 2200 "Gmsh.y"
+#line 2199 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(ROTATE, (yyvsp[(10) - (13)].l), 
@@ -6058,14 +6057,14 @@ yyreduce:
     break;
 
   case 167:
-#line 2208 "Gmsh.y"
+#line 2207 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 168:
-#line 2212 "Gmsh.y"
+#line 2211 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShapes(TRANSLATE_ROTATE, (yyvsp[(12) - (15)].l), 
@@ -6076,14 +6075,14 @@ yyreduce:
     break;
 
   case 169:
-#line 2220 "Gmsh.y"
+#line 2219 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 170:
-#line 2224 "Gmsh.y"
+#line 2223 "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.,
@@ -6093,14 +6092,14 @@ yyreduce:
     break;
 
   case 171:
-#line 2231 "Gmsh.y"
+#line 2230 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 172:
-#line 2235 "Gmsh.y"
+#line 2234 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       extr.mesh.ViewIndex = (int)(yyvsp[(4) - (10)].d);
@@ -6113,7 +6112,7 @@ yyreduce:
     break;
 
   case 173:
-#line 2247 "Gmsh.y"
+#line 2246 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[(4) - (8)].d), 
@@ -6123,7 +6122,7 @@ yyreduce:
     break;
 
   case 174:
-#line 2254 "Gmsh.y"
+#line 2253 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (8)].d), 
@@ -6133,7 +6132,7 @@ yyreduce:
     break;
 
   case 175:
-#line 2261 "Gmsh.y"
+#line 2260 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (8)].d), 
@@ -6143,7 +6142,7 @@ yyreduce:
     break;
 
   case 176:
-#line 2268 "Gmsh.y"
+#line 2267 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[(4) - (12)].d), 
@@ -6153,7 +6152,7 @@ yyreduce:
     break;
 
   case 177:
-#line 2275 "Gmsh.y"
+#line 2274 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (12)].d), 
@@ -6163,7 +6162,7 @@ yyreduce:
     break;
 
   case 178:
-#line 2282 "Gmsh.y"
+#line 2281 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (12)].d), 
@@ -6173,7 +6172,7 @@ yyreduce:
     break;
 
   case 179:
-#line 2289 "Gmsh.y"
+#line 2288 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[(4) - (14)].d), 
@@ -6183,7 +6182,7 @@ yyreduce:
     break;
 
   case 180:
-#line 2296 "Gmsh.y"
+#line 2295 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (14)].d), 
@@ -6193,7 +6192,7 @@ yyreduce:
     break;
 
   case 181:
-#line 2303 "Gmsh.y"
+#line 2302 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (14)].d), 
@@ -6203,14 +6202,14 @@ yyreduce:
     break;
 
   case 182:
-#line 2310 "Gmsh.y"
+#line 2309 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 183:
-#line 2314 "Gmsh.y"
+#line 2313 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_POINT, (int)(yyvsp[(4) - (12)].d), 
@@ -6220,14 +6219,14 @@ yyreduce:
     break;
 
   case 184:
-#line 2321 "Gmsh.y"
+#line 2320 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 185:
-#line 2325 "Gmsh.y"
+#line 2324 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (12)].d), 
@@ -6237,14 +6236,14 @@ yyreduce:
     break;
 
   case 186:
-#line 2332 "Gmsh.y"
+#line 2331 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 187:
-#line 2336 "Gmsh.y"
+#line 2335 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (12)].d), 
@@ -6254,14 +6253,14 @@ yyreduce:
     break;
 
   case 188:
-#line 2343 "Gmsh.y"
+#line 2342 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 189:
-#line 2347 "Gmsh.y"
+#line 2346 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_POINT, (int)(yyvsp[(4) - (16)].d), 
@@ -6271,14 +6270,14 @@ yyreduce:
     break;
 
   case 190:
-#line 2354 "Gmsh.y"
+#line 2353 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 191:
-#line 2358 "Gmsh.y"
+#line 2357 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (16)].d), 
@@ -6288,14 +6287,14 @@ yyreduce:
     break;
 
   case 192:
-#line 2365 "Gmsh.y"
+#line 2364 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 193:
-#line 2369 "Gmsh.y"
+#line 2368 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (16)].d), 
@@ -6305,14 +6304,14 @@ yyreduce:
     break;
 
   case 194:
-#line 2376 "Gmsh.y"
+#line 2375 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 195:
-#line 2380 "Gmsh.y"
+#line 2379 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_POINT, (int)(yyvsp[(4) - (18)].d), 
@@ -6322,14 +6321,14 @@ yyreduce:
     break;
 
   case 196:
-#line 2387 "Gmsh.y"
+#line 2386 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 197:
-#line 2391 "Gmsh.y"
+#line 2390 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SEGM_LINE, (int)(yyvsp[(4) - (18)].d), 
@@ -6339,14 +6338,14 @@ yyreduce:
     break;
 
   case 198:
-#line 2398 "Gmsh.y"
+#line 2397 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = extr.mesh.Recombine = false;
     ;}
     break;
 
   case 199:
-#line 2402 "Gmsh.y"
+#line 2401 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(Shape));
       ExtrudeShape(TRANSLATE_ROTATE, MSH_SURF_PLAN, (int)(yyvsp[(4) - (18)].d), 
@@ -6356,19 +6355,19 @@ yyreduce:
     break;
 
   case 200:
-#line 2413 "Gmsh.y"
+#line 2412 "Gmsh.y"
     {
     ;}
     break;
 
   case 201:
-#line 2416 "Gmsh.y"
+#line 2415 "Gmsh.y"
     {
     ;}
     break;
 
   case 202:
-#line 2422 "Gmsh.y"
+#line 2421 "Gmsh.y"
     {
       extr.mesh.ExtrudeMesh = true;
       extr.mesh.NbLayer = 1;
@@ -6380,7 +6379,7 @@ yyreduce:
     break;
 
   case 203:
-#line 2431 "Gmsh.y"
+#line 2430 "Gmsh.y"
     {
       double d;
       extr.mesh.ExtrudeMesh = true;
@@ -6403,7 +6402,7 @@ yyreduce:
     break;
 
   case 204:
-#line 2451 "Gmsh.y"
+#line 2450 "Gmsh.y"
     {
       yymsg(GERROR, "Explicit region numbers in layers are deprecated");
       double d;
@@ -6429,14 +6428,14 @@ yyreduce:
     break;
 
   case 205:
-#line 2474 "Gmsh.y"
+#line 2473 "Gmsh.y"
     {
       extr.mesh.Recombine = true;
     ;}
     break;
 
   case 206:
-#line 2478 "Gmsh.y"
+#line 2477 "Gmsh.y"
     {
       int num = (int)(yyvsp[(3) - (9)].d);
       if(FindSurface(num)){
@@ -6458,7 +6457,7 @@ yyreduce:
     break;
 
   case 207:
-#line 2502 "Gmsh.y"
+#line 2501 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (6)].l)); i++){
 	double d;
@@ -6479,7 +6478,7 @@ yyreduce:
     break;
 
   case 208:
-#line 2520 "Gmsh.y"
+#line 2519 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (9)].l)); i++){
 	double d;
@@ -6500,7 +6499,7 @@ yyreduce:
     break;
 
   case 209:
-#line 2538 "Gmsh.y"
+#line 2537 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (9)].l)); i++){
 	double d;
@@ -6521,7 +6520,7 @@ yyreduce:
     break;
 
   case 210:
-#line 2556 "Gmsh.y"
+#line 2555 "Gmsh.y"
     {
       Surface *s = FindSurface((int)(yyvsp[(4) - (8)].d));
       if(!s)
@@ -6553,7 +6552,7 @@ yyreduce:
     break;
 
   case 211:
-#line 2585 "Gmsh.y"
+#line 2584 "Gmsh.y"
     {
       Surface *s = FindSurface((int)(yyvsp[(4) - (9)].d));
       if(!s)
@@ -6591,7 +6590,7 @@ yyreduce:
     break;
 
   case 212:
-#line 2620 "Gmsh.y"
+#line 2619 "Gmsh.y"
     {
       yymsg(WARNING, "Elliptic Surface is deprecated: use Transfinite instead (with smoothing)");
       List_Delete((yyvsp[(7) - (8)].l));
@@ -6599,7 +6598,7 @@ yyreduce:
     break;
 
   case 213:
-#line 2625 "Gmsh.y"
+#line 2624 "Gmsh.y"
     {
       Volume *v = FindVolume((int)(yyvsp[(4) - (8)].d));
       if(!v)
@@ -6629,7 +6628,7 @@ yyreduce:
     break;
 
   case 214:
-#line 2652 "Gmsh.y"
+#line 2651 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (6)].l)); i++){
 	double d;
@@ -6646,7 +6645,7 @@ yyreduce:
     break;
 
   case 215:
-#line 2666 "Gmsh.y"
+#line 2665 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (4)].l)); i++){
 	double d;
@@ -6662,7 +6661,7 @@ yyreduce:
     break;
 
   case 216:
-#line 2679 "Gmsh.y"
+#line 2678 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (6)].l)); i++){
 	double d;
@@ -6676,7 +6675,7 @@ yyreduce:
     break;
 
   case 217:
-#line 2696 "Gmsh.y"
+#line 2695 "Gmsh.y"
     { 
       Surface *s = FindSurface((int)(yyvsp[(8) - (10)].d));
       if(s)
@@ -6685,7 +6684,7 @@ yyreduce:
     break;
 
   case 218:
-#line 2702 "Gmsh.y"
+#line 2701 "Gmsh.y"
     {
       Surface *s = FindSurface((int)(yyvsp[(8) - (10)].d));
       if(s)
@@ -6694,66 +6693,66 @@ yyreduce:
     break;
 
   case 219:
-#line 2708 "Gmsh.y"
+#line 2707 "Gmsh.y"
     {
     ;}
     break;
 
   case 220:
-#line 2711 "Gmsh.y"
+#line 2710 "Gmsh.y"
     {
     ;}
     break;
 
   case 221:
-#line 2720 "Gmsh.y"
+#line 2719 "Gmsh.y"
     { 
       ReplaceAllDuplicates();
     ;}
     break;
 
   case 222:
-#line 2729 "Gmsh.y"
+#line 2728 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (1)].d);           ;}
     break;
 
   case 223:
-#line 2730 "Gmsh.y"
+#line 2729 "Gmsh.y"
     { (yyval.d) = (yyvsp[(2) - (3)].d);           ;}
     break;
 
   case 224:
-#line 2731 "Gmsh.y"
+#line 2730 "Gmsh.y"
     { (yyval.d) = -(yyvsp[(2) - (2)].d);          ;}
     break;
 
   case 225:
-#line 2732 "Gmsh.y"
+#line 2731 "Gmsh.y"
     { (yyval.d) = (yyvsp[(2) - (2)].d);           ;}
     break;
 
   case 226:
-#line 2733 "Gmsh.y"
+#line 2732 "Gmsh.y"
     { (yyval.d) = !(yyvsp[(2) - (2)].d);          ;}
     break;
 
   case 227:
-#line 2734 "Gmsh.y"
+#line 2733 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) - (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 228:
-#line 2735 "Gmsh.y"
+#line 2734 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) + (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 229:
-#line 2736 "Gmsh.y"
+#line 2735 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) * (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 230:
-#line 2738 "Gmsh.y"
+#line 2737 "Gmsh.y"
     { 
       if(!(yyvsp[(3) - (3)].d))
 	yymsg(GERROR, "Division by zero in '%g / %g'", (yyvsp[(1) - (3)].d), (yyvsp[(3) - (3)].d));
@@ -6763,307 +6762,307 @@ yyreduce:
     break;
 
   case 231:
-#line 2744 "Gmsh.y"
+#line 2743 "Gmsh.y"
     { (yyval.d) = (int)(yyvsp[(1) - (3)].d) % (int)(yyvsp[(3) - (3)].d);  ;}
     break;
 
   case 232:
-#line 2745 "Gmsh.y"
+#line 2744 "Gmsh.y"
     { (yyval.d) = pow((yyvsp[(1) - (3)].d), (yyvsp[(3) - (3)].d));  ;}
     break;
 
   case 233:
-#line 2746 "Gmsh.y"
+#line 2745 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) < (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 234:
-#line 2747 "Gmsh.y"
+#line 2746 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) > (yyvsp[(3) - (3)].d);      ;}
     break;
 
   case 235:
-#line 2748 "Gmsh.y"
+#line 2747 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) <= (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 236:
-#line 2749 "Gmsh.y"
+#line 2748 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) >= (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 237:
-#line 2750 "Gmsh.y"
+#line 2749 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) == (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 238:
-#line 2751 "Gmsh.y"
+#line 2750 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) != (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 239:
-#line 2752 "Gmsh.y"
+#line 2751 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) && (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 240:
-#line 2753 "Gmsh.y"
+#line 2752 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (3)].d) || (yyvsp[(3) - (3)].d);     ;}
     break;
 
   case 241:
-#line 2754 "Gmsh.y"
+#line 2753 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (5)].d)? (yyvsp[(3) - (5)].d) : (yyvsp[(5) - (5)].d);  ;}
     break;
 
   case 242:
-#line 2755 "Gmsh.y"
+#line 2754 "Gmsh.y"
     { (yyval.d) = exp((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 243:
-#line 2756 "Gmsh.y"
+#line 2755 "Gmsh.y"
     { (yyval.d) = log((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 244:
-#line 2757 "Gmsh.y"
+#line 2756 "Gmsh.y"
     { (yyval.d) = log10((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 245:
-#line 2758 "Gmsh.y"
+#line 2757 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 246:
-#line 2759 "Gmsh.y"
+#line 2758 "Gmsh.y"
     { (yyval.d) = sin((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 247:
-#line 2760 "Gmsh.y"
+#line 2759 "Gmsh.y"
     { (yyval.d) = asin((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 248:
-#line 2761 "Gmsh.y"
+#line 2760 "Gmsh.y"
     { (yyval.d) = cos((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 249:
-#line 2762 "Gmsh.y"
+#line 2761 "Gmsh.y"
     { (yyval.d) = acos((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 250:
-#line 2763 "Gmsh.y"
+#line 2762 "Gmsh.y"
     { (yyval.d) = tan((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 251:
-#line 2764 "Gmsh.y"
+#line 2763 "Gmsh.y"
     { (yyval.d) = atan((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 252:
-#line 2765 "Gmsh.y"
+#line 2764 "Gmsh.y"
     { (yyval.d) = atan2((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d));;}
     break;
 
   case 253:
-#line 2766 "Gmsh.y"
+#line 2765 "Gmsh.y"
     { (yyval.d) = sinh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 254:
-#line 2767 "Gmsh.y"
+#line 2766 "Gmsh.y"
     { (yyval.d) = cosh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 255:
-#line 2768 "Gmsh.y"
+#line 2767 "Gmsh.y"
     { (yyval.d) = tanh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 256:
-#line 2769 "Gmsh.y"
+#line 2768 "Gmsh.y"
     { (yyval.d) = fabs((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 257:
-#line 2770 "Gmsh.y"
+#line 2769 "Gmsh.y"
     { (yyval.d) = floor((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 258:
-#line 2771 "Gmsh.y"
+#line 2770 "Gmsh.y"
     { (yyval.d) = ceil((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 259:
-#line 2772 "Gmsh.y"
+#line 2771 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 260:
-#line 2773 "Gmsh.y"
+#line 2772 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 261:
-#line 2774 "Gmsh.y"
+#line 2773 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (6)].d)*(yyvsp[(3) - (6)].d)+(yyvsp[(5) - (6)].d)*(yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 262:
-#line 2775 "Gmsh.y"
+#line 2774 "Gmsh.y"
     { (yyval.d) = (yyvsp[(3) - (4)].d)*(double)rand()/(double)RAND_MAX; ;}
     break;
 
   case 263:
-#line 2777 "Gmsh.y"
+#line 2776 "Gmsh.y"
     { (yyval.d) = exp((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 264:
-#line 2778 "Gmsh.y"
+#line 2777 "Gmsh.y"
     { (yyval.d) = log((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 265:
-#line 2779 "Gmsh.y"
+#line 2778 "Gmsh.y"
     { (yyval.d) = log10((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 266:
-#line 2780 "Gmsh.y"
+#line 2779 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 267:
-#line 2781 "Gmsh.y"
+#line 2780 "Gmsh.y"
     { (yyval.d) = sin((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 268:
-#line 2782 "Gmsh.y"
+#line 2781 "Gmsh.y"
     { (yyval.d) = asin((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 269:
-#line 2783 "Gmsh.y"
+#line 2782 "Gmsh.y"
     { (yyval.d) = cos((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 270:
-#line 2784 "Gmsh.y"
+#line 2783 "Gmsh.y"
     { (yyval.d) = acos((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 271:
-#line 2785 "Gmsh.y"
+#line 2784 "Gmsh.y"
     { (yyval.d) = tan((yyvsp[(3) - (4)].d));      ;}
     break;
 
   case 272:
-#line 2786 "Gmsh.y"
+#line 2785 "Gmsh.y"
     { (yyval.d) = atan((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 273:
-#line 2787 "Gmsh.y"
+#line 2786 "Gmsh.y"
     { (yyval.d) = atan2((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d));;}
     break;
 
   case 274:
-#line 2788 "Gmsh.y"
+#line 2787 "Gmsh.y"
     { (yyval.d) = sinh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 275:
-#line 2789 "Gmsh.y"
+#line 2788 "Gmsh.y"
     { (yyval.d) = cosh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 276:
-#line 2790 "Gmsh.y"
+#line 2789 "Gmsh.y"
     { (yyval.d) = tanh((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 277:
-#line 2791 "Gmsh.y"
+#line 2790 "Gmsh.y"
     { (yyval.d) = fabs((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 278:
-#line 2792 "Gmsh.y"
+#line 2791 "Gmsh.y"
     { (yyval.d) = floor((yyvsp[(3) - (4)].d));    ;}
     break;
 
   case 279:
-#line 2793 "Gmsh.y"
+#line 2792 "Gmsh.y"
     { (yyval.d) = ceil((yyvsp[(3) - (4)].d));     ;}
     break;
 
   case 280:
-#line 2794 "Gmsh.y"
+#line 2793 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 281:
-#line 2795 "Gmsh.y"
+#line 2794 "Gmsh.y"
     { (yyval.d) = fmod((yyvsp[(3) - (6)].d), (yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 282:
-#line 2796 "Gmsh.y"
+#line 2795 "Gmsh.y"
     { (yyval.d) = sqrt((yyvsp[(3) - (6)].d)*(yyvsp[(3) - (6)].d)+(yyvsp[(5) - (6)].d)*(yyvsp[(5) - (6)].d)); ;}
     break;
 
   case 283:
-#line 2797 "Gmsh.y"
+#line 2796 "Gmsh.y"
     { (yyval.d) = (yyvsp[(3) - (4)].d)*(double)rand()/(double)RAND_MAX; ;}
     break;
 
   case 284:
-#line 2806 "Gmsh.y"
+#line 2805 "Gmsh.y"
     { (yyval.d) = (yyvsp[(1) - (1)].d); ;}
     break;
 
   case 285:
-#line 2807 "Gmsh.y"
+#line 2806 "Gmsh.y"
     { (yyval.d) = 3.141592653589793; ;}
     break;
 
   case 286:
-#line 2808 "Gmsh.y"
+#line 2807 "Gmsh.y"
     { (yyval.d) = ParUtil::Instance()->rank(); ;}
     break;
 
   case 287:
-#line 2809 "Gmsh.y"
+#line 2808 "Gmsh.y"
     { (yyval.d) = ParUtil::Instance()->size(); ;}
     break;
 
   case 288:
-#line 2810 "Gmsh.y"
+#line 2809 "Gmsh.y"
     { (yyval.d) = Get_GmshMajorVersion(); ;}
     break;
 
   case 289:
-#line 2811 "Gmsh.y"
+#line 2810 "Gmsh.y"
     { (yyval.d) = Get_GmshMinorVersion(); ;}
     break;
 
   case 290:
-#line 2812 "Gmsh.y"
+#line 2811 "Gmsh.y"
     { (yyval.d) = Get_GmshPatchVersion(); ;}
     break;
 
   case 291:
-#line 2817 "Gmsh.y"
+#line 2816 "Gmsh.y"
     {
       Symbol TheSymbol;
       TheSymbol.Name = (yyvsp[(1) - (1)].c);
@@ -7079,7 +7078,7 @@ yyreduce:
     break;
 
   case 292:
-#line 2833 "Gmsh.y"
+#line 2832 "Gmsh.y"
     {
       char tmpstring[1024];
       sprintf(tmpstring, "%s_%d", (yyvsp[(1) - (5)].c), (int)(yyvsp[(4) - (5)].d)) ;
@@ -7097,7 +7096,7 @@ yyreduce:
     break;
 
   case 293:
-#line 2848 "Gmsh.y"
+#line 2847 "Gmsh.y"
     {
       Symbol TheSymbol;
       TheSymbol.Name = (yyvsp[(1) - (4)].c);
@@ -7120,7 +7119,7 @@ yyreduce:
     break;
 
   case 294:
-#line 2868 "Gmsh.y"
+#line 2867 "Gmsh.y"
     {
       Symbol TheSymbol;
       TheSymbol.Name = (yyvsp[(2) - (4)].c);
@@ -7136,7 +7135,7 @@ yyreduce:
     break;
 
   case 295:
-#line 2881 "Gmsh.y"
+#line 2880 "Gmsh.y"
     {
       Symbol TheSymbol;
       TheSymbol.Name = (yyvsp[(1) - (2)].c);
@@ -7152,7 +7151,7 @@ yyreduce:
     break;
 
   case 296:
-#line 2894 "Gmsh.y"
+#line 2893 "Gmsh.y"
     {
       Symbol TheSymbol;
       TheSymbol.Name = (yyvsp[(1) - (5)].c);
@@ -7175,7 +7174,7 @@ yyreduce:
     break;
 
   case 297:
-#line 2917 "Gmsh.y"
+#line 2916 "Gmsh.y"
     {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -7196,7 +7195,7 @@ yyreduce:
     break;
 
   case 298:
-#line 2935 "Gmsh.y"
+#line 2934 "Gmsh.y"
     {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -7217,7 +7216,7 @@ yyreduce:
     break;
 
   case 299:
-#line 2953 "Gmsh.y"
+#line 2952 "Gmsh.y"
     {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -7238,7 +7237,7 @@ yyreduce:
     break;
 
   case 300:
-#line 2971 "Gmsh.y"
+#line 2970 "Gmsh.y"
     {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -7259,7 +7258,7 @@ yyreduce:
     break;
 
   case 301:
-#line 2989 "Gmsh.y"
+#line 2988 "Gmsh.y"
     { 
       (yyval.d) = GetValue((yyvsp[(3) - (6)].c), (yyvsp[(5) - (6)].d));
       Free((yyvsp[(3) - (6)].c));
@@ -7267,70 +7266,70 @@ yyreduce:
     break;
 
   case 302:
-#line 2997 "Gmsh.y"
+#line 2996 "Gmsh.y"
     {
       memcpy((yyval.v), (yyvsp[(1) - (1)].v), 5*sizeof(double));
     ;}
     break;
 
   case 303:
-#line 3001 "Gmsh.y"
+#line 3000 "Gmsh.y"
     {
       for(int i = 0; i < 5; i++) (yyval.v)[i] = -(yyvsp[(2) - (2)].v)[i];
     ;}
     break;
 
   case 304:
-#line 3005 "Gmsh.y"
+#line 3004 "Gmsh.y"
     { 
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(2) - (2)].v)[i];
     ;}
     break;
 
   case 305:
-#line 3009 "Gmsh.y"
+#line 3008 "Gmsh.y"
     { 
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(1) - (3)].v)[i] - (yyvsp[(3) - (3)].v)[i];
     ;}
     break;
 
   case 306:
-#line 3013 "Gmsh.y"
+#line 3012 "Gmsh.y"
     {
       for(int i = 0; i < 5; i++) (yyval.v)[i] = (yyvsp[(1) - (3)].v)[i] + (yyvsp[(3) - (3)].v)[i];
     ;}
     break;
 
   case 307:
-#line 3020 "Gmsh.y"
+#line 3019 "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 308:
-#line 3024 "Gmsh.y"
+#line 3023 "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 309:
-#line 3028 "Gmsh.y"
+#line 3027 "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 310:
-#line 3032 "Gmsh.y"
+#line 3031 "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 311:
-#line 3039 "Gmsh.y"
+#line 3038 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(List_T*));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].l)));
@@ -7338,14 +7337,14 @@ yyreduce:
     break;
 
   case 312:
-#line 3044 "Gmsh.y"
+#line 3043 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].l)));
     ;}
     break;
 
   case 313:
-#line 3052 "Gmsh.y"
+#line 3051 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].d)));
@@ -7353,14 +7352,14 @@ yyreduce:
     break;
 
   case 314:
-#line 3057 "Gmsh.y"
+#line 3056 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(1) - (1)].l);
     ;}
     break;
 
   case 315:
-#line 3061 "Gmsh.y"
+#line 3060 "Gmsh.y"
     {
       // creates an empty list
       (yyval.l) = List_Create(2, 1, sizeof(double));
@@ -7368,14 +7367,14 @@ yyreduce:
     break;
 
   case 316:
-#line 3066 "Gmsh.y"
+#line 3065 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (3)].l);
     ;}
     break;
 
   case 317:
-#line 3070 "Gmsh.y"
+#line 3069 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(3) - (4)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7386,7 +7385,7 @@ yyreduce:
     break;
 
   case 318:
-#line 3078 "Gmsh.y"
+#line 3077 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(4) - (5)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7397,7 +7396,7 @@ yyreduce:
     break;
 
   case 319:
-#line 3089 "Gmsh.y"
+#line 3088 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (2)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7408,7 +7407,7 @@ yyreduce:
     break;
 
   case 320:
-#line 3097 "Gmsh.y"
+#line 3096 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(3) - (3)].l);
       for(int i = 0; i < List_Nbr((yyval.l)); i++){
@@ -7419,7 +7418,7 @@ yyreduce:
     break;
 
   case 321:
-#line 3105 "Gmsh.y"
+#line 3104 "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)); ((yyvsp[(1) - (3)].d) < (yyvsp[(3) - (3)].d)) ? (d += 1.) : (d -= 1.)) 
@@ -7428,7 +7427,7 @@ yyreduce:
     break;
 
   case 322:
-#line 3111 "Gmsh.y"
+#line 3110 "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)){
@@ -7442,7 +7441,7 @@ yyreduce:
     break;
 
   case 323:
-#line 3122 "Gmsh.y"
+#line 3121 "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
@@ -7465,7 +7464,7 @@ yyreduce:
     break;
 
   case 324:
-#line 3142 "Gmsh.y"
+#line 3141 "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++){
@@ -7478,7 +7477,7 @@ yyreduce:
     break;
 
   case 325:
-#line 3152 "Gmsh.y"
+#line 3151 "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++){
@@ -7491,7 +7490,7 @@ yyreduce:
     break;
 
   case 326:
-#line 3162 "Gmsh.y"
+#line 3161 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -7511,7 +7510,7 @@ yyreduce:
     break;
 
   case 327:
-#line 3179 "Gmsh.y"
+#line 3178 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -7538,7 +7537,7 @@ yyreduce:
     break;
 
   case 328:
-#line 3206 "Gmsh.y"
+#line 3205 "Gmsh.y"
     {
       (yyval.l) = List_Create(2, 1, sizeof(double));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].d)));
@@ -7546,21 +7545,21 @@ yyreduce:
     break;
 
   case 329:
-#line 3211 "Gmsh.y"
+#line 3210 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(1) - (1)].l);
     ;}
     break;
 
   case 330:
-#line 3215 "Gmsh.y"
+#line 3214 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].d)));
     ;}
     break;
 
   case 331:
-#line 3219 "Gmsh.y"
+#line 3218 "Gmsh.y"
     {
       for(int i = 0; i < List_Nbr((yyvsp[(3) - (3)].l)); i++){
 	double d;
@@ -7572,21 +7571,21 @@ yyreduce:
     break;
 
   case 332:
-#line 3232 "Gmsh.y"
+#line 3231 "Gmsh.y"
     {
       (yyval.u) = CTX.PACK_COLOR((int)(yyvsp[(2) - (9)].d), (int)(yyvsp[(4) - (9)].d), (int)(yyvsp[(6) - (9)].d), (int)(yyvsp[(8) - (9)].d));
     ;}
     break;
 
   case 333:
-#line 3236 "Gmsh.y"
+#line 3235 "Gmsh.y"
     {
       (yyval.u) = CTX.PACK_COLOR((int)(yyvsp[(2) - (7)].d), (int)(yyvsp[(4) - (7)].d), (int)(yyvsp[(6) - (7)].d), 255);
     ;}
     break;
 
   case 334:
-#line 3248 "Gmsh.y"
+#line 3247 "Gmsh.y"
     {
       int flag;
       (yyval.u) = Get_ColorForString(ColorString, -1, (yyvsp[(1) - (1)].c), &flag);
@@ -7596,7 +7595,7 @@ yyreduce:
     break;
 
   case 335:
-#line 3255 "Gmsh.y"
+#line 3254 "Gmsh.y"
     {
       unsigned int (*pColOpt)(int num, int action, unsigned int value);
       StringXColor *pColCat;
@@ -7617,14 +7616,14 @@ yyreduce:
     break;
 
   case 336:
-#line 3276 "Gmsh.y"
+#line 3275 "Gmsh.y"
     {
       (yyval.l) = (yyvsp[(2) - (3)].l);
     ;}
     break;
 
   case 337:
-#line 3280 "Gmsh.y"
+#line 3279 "Gmsh.y"
     {
       (yyval.l) = List_Create(256, 10, sizeof(unsigned int));
       GmshColorTable *ct = Get_ColorTable((int)(yyvsp[(3) - (6)].d));
@@ -7639,7 +7638,7 @@ yyreduce:
     break;
 
   case 338:
-#line 3295 "Gmsh.y"
+#line 3294 "Gmsh.y"
     {
       (yyval.l) = List_Create(256, 10, sizeof(unsigned int));
       List_Add((yyval.l), &((yyvsp[(1) - (1)].u)));
@@ -7647,35 +7646,35 @@ yyreduce:
     break;
 
   case 339:
-#line 3300 "Gmsh.y"
+#line 3299 "Gmsh.y"
     {
       List_Add((yyval.l), &((yyvsp[(3) - (3)].u)));
     ;}
     break;
 
   case 340:
-#line 3307 "Gmsh.y"
+#line 3306 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(1) - (1)].c);
     ;}
     break;
 
   case 341:
-#line 3311 "Gmsh.y"
+#line 3310 "Gmsh.y"
     {
       Msg(WARNING, "Named string expressions not implemented yet");
     ;}
     break;
 
   case 342:
-#line 3318 "Gmsh.y"
+#line 3317 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(1) - (1)].c);
     ;}
     break;
 
   case 343:
-#line 3322 "Gmsh.y"
+#line 3321 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc(32*sizeof(char));
       time_t now;
@@ -7686,7 +7685,7 @@ yyreduce:
     break;
 
   case 344:
-#line 3330 "Gmsh.y"
+#line 3329 "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));
@@ -7697,7 +7696,7 @@ yyreduce:
     break;
 
   case 345:
-#line 3338 "Gmsh.y"
+#line 3337 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc((strlen((yyvsp[(3) - (4)].c))+1)*sizeof(char));
       int i;
@@ -7714,7 +7713,7 @@ yyreduce:
     break;
 
   case 346:
-#line 3352 "Gmsh.y"
+#line 3351 "Gmsh.y"
     {
       (yyval.c) = (char *)Malloc((strlen((yyvsp[(3) - (4)].c))+1)*sizeof(char));
       int i;
@@ -7731,14 +7730,14 @@ yyreduce:
     break;
 
   case 347:
-#line 3366 "Gmsh.y"
+#line 3365 "Gmsh.y"
     {
       (yyval.c) = (yyvsp[(3) - (4)].c);
     ;}
     break;
 
   case 348:
-#line 3370 "Gmsh.y"
+#line 3369 "Gmsh.y"
     {
       char tmpstring[1024];
       int i = PrintListOfDouble((yyvsp[(3) - (6)].c), (yyvsp[(5) - (6)].l), tmpstring);
@@ -7760,7 +7759,7 @@ yyreduce:
     break;
 
   case 349:
-#line 3389 "Gmsh.y"
+#line 3388 "Gmsh.y"
     { 
       const char* (*pStrOpt)(int num, int action, const char *value);
       StringXString *pStrCat;
@@ -7786,7 +7785,7 @@ yyreduce:
     break;
 
   case 350:
-#line 3412 "Gmsh.y"
+#line 3411 "Gmsh.y"
     { 
       const char* (*pStrOpt)(int num, int action, const char *value);
       StringXString *pStrCat;
@@ -7813,7 +7812,7 @@ yyreduce:
 
 
 /* Line 1267 of yacc.c.  */
-#line 7817 "Gmsh.tab.cpp"
+#line 7816 "Gmsh.tab.cpp"
       default: break;
     }
   YY_SYMBOL_PRINT ("-> $$ =", yyr1[yyn], &yyval, &yyloc);
@@ -8027,7 +8026,7 @@ yyreturn:
 }
 
 
-#line 3436 "Gmsh.y"
+#line 3435 "Gmsh.y"
 
 
 void DeleteSymbol(void *a, void *b)
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index e88819f6e58774a804ae6e2a38cf808cda0dca4b..dff43642403c66f3fbaed5314affc79f854a5244 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1,5 +1,5 @@
 %{
-// $Id: Gmsh.y,v 1.306 2008-03-21 07:21:07 geuzaine Exp $
+// $Id: Gmsh.y,v 1.307 2008-03-23 21:43:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -1916,9 +1916,8 @@ Command :
     {
       if(!strcmp($1, "Background") && !strcmp($2, "Mesh")  && !strcmp($3, "View")){
 	int index = (int)$5;
-	if(index >= 0 && index < (int)PView::list.size()){
+	if(index >= 0 && index < (int)PView::list.size())
 	  GModel::current()->getFields()->set_background_mesh(index);
-	}
 	else
 	  yymsg(GERROR, "Unknown view %d", index);
       }
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index 97d69ab9f3ab7a8943dcd7f85fcfb87437d0b29a..fdefba39819f98a8aacf3e32bfaa01917ac2406d 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -835,7 +835,7 @@ int gmsh_yy_flex_debug = 0;
 char *gmsh_yytext;
 #line 1 "Gmsh.l"
 #line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.355 2008-03-21 09:39:41 geuzaine Exp $
+// $Id: Gmsh.yy.cpp,v 1.356 2008-03-23 21:43:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 6daebad69e9a37c64a104b5fcc6afd1cfd06a0ba..0e5dc7527f9b0856112ca10cb0ae762d0d160b97 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.179 2008-03-20 11:44:12 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.180 2008-03-23 21:43:03 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -274,6 +274,8 @@ void SetProjectName(const char *name)
   strncpy(CTX.no_ext_filename, no_ext, 255);
   strncpy(CTX.base_filename, base, 255);
 
+  GModel::current()->setName(base);
+    
 #if defined(HAVE_FLTK)
   if(!CTX.batch) WID->set_title(CTX.filename);
 #endif
@@ -350,6 +352,9 @@ int MergeFile(const char *name, int warn_if_missing)
   else if(!strcmp(ext, ".mesh") || !strcmp(ext, ".MESH")){
     status = m->readMESH(name);
   }
+  else if(!strcmp(ext, ".med") || !strcmp(ext, ".MED")){
+    status = m->readMED(name);
+  }
   else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF") ||
           !strcmp(ext, ".nas") || !strcmp(ext, ".NAS")){
     status = m->readBDF(name);
diff --git a/benchmarks/2d/function_field.geo b/benchmarks/2d/function_field.geo
index e6f6b54de32580495e6a4e0a132825890d24b78e..6f9744a8838a98bfa2027f005b5e10d15cbe1ced 100644
--- a/benchmarks/2d/function_field.geo
+++ b/benchmarks/2d/function_field.geo
@@ -11,5 +11,6 @@ Line(4) = {4,1} ;
 Line Loop(5) = {4,1,-2,3} ;
 Plane Surface(6) = {5} ;
 
-Function Field(1) = "Cos(2*3.14*x)/5 + 0.21";
-Characteristic Length Field{1};
+Field[1] = MathEval;
+Field[1].F = "Cos(2*3.14*x)/5 + 0.21";
+Background Field = 1;
diff --git a/benchmarks/3d/function_field.geo b/benchmarks/3d/function_field.geo
index 1ed4680ffe137b09951cfe43b3ef3e9de5ee76ab..b3b95d7a9f1d0860c3b9938b803b65a4ef0aeeaa 100644
--- a/benchmarks/3d/function_field.geo
+++ b/benchmarks/3d/function_field.geo
@@ -12,5 +12,6 @@ Line Loop(5) = {4,1,-2,3} ;
 Plane Surface(6) = {5} ;
 Extrude {0,0,1} { Surface{6}; }
 
-Function Field(1) = "Cos(2*3.14*x)/5 + 0.21";
-Characteristic Length Field{1};
+Field[1] = MathEval;
+Field[1].F = "Cos(2*3.14*x)/5 + 0.21";
+Background Field = 1;