diff --git a/Fltk/GUI_Projection.cpp b/Fltk/GUI_Projection.cpp
index b94db3c61411a0e7bd6850e36d8721d2f7209c9b..7a3335edc8af8ad70fad5cf1f1946fad5e7cae52 100644
--- a/Fltk/GUI_Projection.cpp
+++ b/Fltk/GUI_Projection.cpp
@@ -26,7 +26,7 @@ extern Context_T CTX;
 static fourierProjectionFace *createProjectionFaceFromName(char *name)
 {
   GModel *m = GModel::current();
-  int tag = m->numFace() + 1;
+  int tag = m->getNumFaces() + 1;
   fourierProjectionFace *f = 0;
   if(!strcmp(name, "plane"))
     f = new fourierProjectionFace(m, tag, new FM::PlaneProjectionSurface(tag));
@@ -238,7 +238,7 @@ projectionEditor::projectionEditor()
   Fl_Group *o = new Fl_Group(WB, WB, 2 * BB, 3 * BH);
   _select[0] = new Fl_Round_Button(2 * WB + BB / 2, WB, BB, BH, "Points");
   _select[1] = new Fl_Round_Button(2 * WB + BB / 2, WB + BH, BB, BH, "Elements");
-  if(GModel::current()->numElements())
+  if(GModel::current()->getNumMeshElements())
     _select[1]->value(1);
   else
     _select[0]->value(1);
@@ -1033,7 +1033,7 @@ void action_cb(Fl_Widget *w, void *data)
     for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++)
       if((*it)->getNativeType() == GEntity::FourierModel) 
 	id = std::max(id, (*it)->tag());
-    if(id > 0) faces.push_back(m->faceByTag(id));
+    if(id > 0) faces.push_back(m->getFace(id));
   }
   else if(what == "delete_all" || what == "save_all"){
     for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index ddbc19cd5df2a356ff8e6ac1a526f0ad946ce0f6..764624c158ec0fc03b199863421195250b4501d9 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.58 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GModel.cpp,v 1.59 2008-02-20 09:20:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -37,7 +37,7 @@ extern Context_T CTX;
 std::vector<GModel*> GModel::list;
 
 GModel::GModel(std::string name)
-  : geo_internals(0), occ_internals(0), modelName(name), normals(0)
+  : _geo_internals(0), _occ_internals(0), modelName(name), normals(0)
 {
   list.push_back(this);
   // at the moment we always create (at least an empty) GEO model
@@ -92,6 +92,8 @@ void GModel::destroy()
   if(normals) delete normals;
   normals = 0;
 
+  invalidateMeshVertexCache();
+
   MVertex::resetGlobalNumber();
   MElement::resetGlobalNumber();
 #if !defined(HAVE_GMSH_EMBEDDED)
@@ -101,40 +103,40 @@ void GModel::destroy()
 #endif
 }
 
-GRegion *GModel::regionByTag(int n) const
+GRegion *GModel::getRegion(int n) const
 {
   GEntity tmp((GModel*)this, n);
-  std::set<GRegion*,GEntityLessThan>::const_iterator it =regions.find((GRegion*)&tmp);
+  std::set<GRegion*, GEntityLessThan>::const_iterator it = regions.find((GRegion*)&tmp);
   if(it != regions.end())
     return *it;
   else
     return 0;
 }
 
-GFace *GModel::faceByTag(int n) const
+GFace *GModel::getFace(int n) const
 {
   GEntity tmp((GModel*)this, n);
-  std::set<GFace*,GEntityLessThan>::const_iterator it = faces.find((GFace*)&tmp);
+  std::set<GFace*, GEntityLessThan>::const_iterator it = faces.find((GFace*)&tmp);
   if(it != faces.end())
     return *it;
   else
     return 0;
 }
 
-GEdge *GModel::edgeByTag(int n) const
+GEdge *GModel::getEdge(int n) const
 {
   GEntity tmp((GModel*)this, n);
-  std::set<GEdge*,GEntityLessThan>::const_iterator  it = edges.find((GEdge*)&tmp);
+  std::set<GEdge*, GEntityLessThan>::const_iterator it = edges.find((GEdge*)&tmp);
   if(it != edges.end())
     return *it;
   else
     return 0;
 }
 
-GVertex *GModel::vertexByTag(int n) const
+GVertex *GModel::getVertex(int n) const
 {
   GEntity tmp((GModel*)this, n);
- std::set<GVertex*,GEntityLessThan>::const_iterator  it = vertices.find((GVertex*)&tmp);
+  std::set<GVertex*, GEntityLessThan>::const_iterator it = vertices.find((GVertex*)&tmp);
   if(it != vertices.end())
     return *it;
   else
@@ -144,8 +146,7 @@ GVertex *GModel::vertexByTag(int n) const
 void GModel::remove(GRegion *r)
 {
   riter it = std::find(firstRegion(), lastRegion(), r);
-  if(it != (riter)regions.end())
-  	regions.erase(it);
+  if(it != (riter)regions.end()) regions.erase(it);
 }
 
 void GModel::remove(GFace *f)
@@ -237,6 +238,8 @@ void GModel::associateEntityWithVertices()
 
 int GModel::renumberMeshVertices(bool saveAll)
 {
+  invalidateMeshVertexCache();
+
   // tag all mesh vertices with -1 (negative vertices will not be
   // saved)
   for(viter it = firstVertex(); it != lastVertex(); ++it)
@@ -472,7 +475,7 @@ int GModel::getMeshStatus(bool countDiscrete)
   return -1;
 }
 
-int GModel::numVertices()
+int GModel::getNumMeshVertices()
 {
   int n = 0;
   for(viter it = firstVertex(); it != lastVertex(); ++it)
@@ -486,7 +489,7 @@ int GModel::numVertices()
   return n;
 }
 
-int GModel::numElements()
+int GModel::getNumMeshElements()
 {
   int n = 0;
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
@@ -504,6 +507,59 @@ int GModel::numElements()
   return n;
 }
 
+template <class T>
+static void insertMeshVertices(std::vector<MVertex*> &vertices, T &container)
+{
+  for(unsigned int i = 0; i < vertices.size(); i++)
+    container[vertices[i]->getNum()] = vertices[i];
+}
+
+void GModel::buildMeshVertexCache()
+{
+  _vertexVectorCache.clear();
+  _vertexMapCache.clear();
+  bool dense = (getNumMeshVertices() == MVertex::getGlobalNumber());
+ 
+  if(dense){
+    _vertexVectorCache.resize(MVertex::getGlobalNumber());
+    for(viter it = firstVertex(); it != lastVertex(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexVectorCache);
+    for(eiter it = firstEdge(); it != lastEdge(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexVectorCache);
+    for(fiter it = firstFace(); it != lastFace(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexVectorCache);
+    for(riter it = firstRegion(); it != lastRegion(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexVectorCache);
+  }
+  else{
+    for(viter it = firstVertex(); it != lastVertex(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexMapCache);
+    for(eiter it = firstEdge(); it != lastEdge(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexMapCache);
+    for(fiter it = firstFace(); it != lastFace(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexMapCache);
+    for(riter it = firstRegion(); it != lastRegion(); ++it)
+      insertMeshVertices((*it)->mesh_vertices, _vertexMapCache);
+  }
+}
+
+void GModel::invalidateMeshVertexCache()
+{
+  _vertexVectorCache.clear();
+  _vertexMapCache.clear();
+}
+
+MVertex *GModel::getMeshVertex(int num)
+{
+  if(num < 0) return 0;
+  if(_vertexVectorCache.empty() && _vertexMapCache.empty())
+    buildMeshVertexCache();
+  if(num < _vertexVectorCache.size())
+    return _vertexVectorCache[num];
+  else
+    return _vertexMapCache[num];
+}
+
 std::set<int> &GModel::recomputeMeshPartitions()
 {
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
@@ -582,7 +638,7 @@ static int checkElements(int tag,
 
 void GModel::checkMeshCoherence()
 {
-  int numEle = numElements();
+  int numEle = getNumMeshElements();
   if(!numEle) return;
 
   Msg(INFO, "Checking mesh coherence (%d elements)", numEle);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index b18423dfd71bb71c78168f2815142da345d35ec6..0481ddc2e510dee9dd570155fc8d796428018190 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -37,16 +37,17 @@ class smooth_normals;
 class GModel
 {
  private:
+  // A vertex cache to speed-up direct access by vertex number (used
+  // for post-processing I/O)
+  std::vector<MVertex*> _vertexVectorCache;
+  std::map<int, MVertex*> _vertexMapCache;
+
+  GEO_Internals *_geo_internals;
   void createGEOInternals();
   void deleteGEOInternals();
-  GEO_Internals *geo_internals;
 
+  OCC_Internals *_occ_internals;
   void deleteOCCInternals();
-  OCC_Internals *occ_internals;
-
- public:
-  GEO_Internals *getGEOInternals(){ return geo_internals; }
-  OCC_Internals *getOCCInternals(){ return occ_internals; }
 
  protected:
   std::string modelName;
@@ -63,26 +64,28 @@ class GModel
 
   // the static list of all loaded models
   static std::vector<GModel*> list;
+
   // returns the current model
   static GModel *current();
 
-  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
-  typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
-  typedef std::set<GEdge*, GEntityLessThan>::iterator eiter;
-  typedef std::set<GVertex*, GEntityLessThan>::iterator viter;
-  typedef std::map<int, std::string>::iterator piter;
-
   // Deletes everything in a GModel
   void destroy();
 
-  // Returns the geometric tolerance for the entire model.
-  double tolerance() const { return 1.e-14; }
+  // Access internal CAD representations
+  GEO_Internals *getGEOInternals(){ return _geo_internals; }
+  OCC_Internals *getOCCInternals(){ return _occ_internals; }
 
   // Get the number of regions in this model.
-  int numRegion() const { return regions.size(); }
-  int numFace() const { return faces.size(); }
-  int numEdge() const { return edges.size(); }
-  int numVertex() const  { return vertices.size(); }
+  int getNumRegions() const { return regions.size(); }
+  int getNumFaces() const { return faces.size(); }
+  int getNumEdges() const { return edges.size(); }
+  int getNumVertices() const  { return vertices.size(); }
+
+  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
+  typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
+  typedef std::set<GEdge*, GEntityLessThan>::iterator eiter;
+  typedef std::set<GVertex*, GEntityLessThan>::iterator viter;
+  typedef std::map<int, std::string>::iterator piter;
 
   // Get an iterator initialized to the first entity in this model.
   riter firstRegion() { return regions.begin(); }
@@ -95,10 +98,10 @@ class GModel
   viter lastVertex() { return vertices.end(); }
 
   // Find the region with the given tag.
-  GRegion *regionByTag(int n) const;
-  GFace *faceByTag(int n) const;
-  GEdge *edgeByTag(int n) const;
-  GVertex *vertexByTag(int n) const;
+  GRegion *getRegion(int n) const;
+  GFace *getFace(int n) const;
+  GEdge *getEdge(int n) const;
+  GVertex *getVertex(int n) const;
 
   void add(GRegion *r) { regions.insert(r); }
   void add(GFace *f) { faces.insert(f); }
@@ -152,10 +155,17 @@ class GModel
   int getMeshStatus(bool countDiscrete=true);
 
   // Returns the total number of vertices in the mesh
-  int numVertices();
+  int getNumMeshVertices();
 
   // Returns the total number of elements in the mesh
-  int numElements();
+  int getNumMeshElements();
+
+  // Invalidate/Rebuild the vertex cache
+  void invalidateMeshVertexCache();
+  void buildMeshVertexCache();
+
+  // Access a mesh vertex by number, using the vertex cache
+  MVertex *getMeshVertex(int num);
 
   // The list of partitions
   std::set<int> &getMeshPartitions() { return meshPartitions; }
diff --git a/Geo/GModelIO_Fourier.cpp b/Geo/GModelIO_Fourier.cpp
index 61606d9a380b544c4010b5c95a00e4290d4c5923..f81003357e96757ce9d17192cd71eec1ada80560 100644
--- a/Geo/GModelIO_Fourier.cpp
+++ b/Geo/GModelIO_Fourier.cpp
@@ -28,7 +28,7 @@ void makeGFace(FM::Patch* patch)
   
   GModel *m = GModel::current();
 
-  int tagVertex = m->numVertex();
+  int tagVertex = m->getNumVertices();
   patch->F(LL[0], LL[1], xx, yy, zz);
   FM::TopoVertex* vLL = new FM::TopoVertex(++tagVertex, xx, yy, zz);
   m->add(new fourierVertex(m, vLL->GetTag(), vLL));
@@ -47,35 +47,35 @@ void makeGFace(FM::Patch* patch)
   FM::Curve* curveT = new FM::FCurve(0, patch, UR, UL);
   FM::Curve* curveL = new FM::FCurve(0, patch, UL, LL);
   
-  int tagEdge = m->numEdge();
+  int tagEdge = m->getNumEdges();
   FM::TopoEdge* eB = new FM::TopoEdge(++tagEdge, curveB, vLL, vLR);
   i1 = eB->GetStartPoint()->GetTag();
   i2 = eB->GetEndPoint()->GetTag();
-  m->add(new fourierEdge(m, eB, eB->GetTag(), m->vertexByTag(i1),
-			 m->vertexByTag(i2)));
+  m->add(new fourierEdge(m, eB, eB->GetTag(), m->getVertex(i1),
+			 m->getVertex(i2)));
   FM::TopoEdge* eR = new FM::TopoEdge(++tagEdge, curveR, vLR, vUR); 
   i1 = eR->GetStartPoint()->GetTag();
   i2 = eR->GetEndPoint()->GetTag();
-  m->add(new fourierEdge(m, eR, eR->GetTag(), m->vertexByTag(i1),
-			 m->vertexByTag(i2))); 
+  m->add(new fourierEdge(m, eR, eR->GetTag(), m->getVertex(i1),
+			 m->getVertex(i2))); 
   FM::TopoEdge* eT = new FM::TopoEdge(++tagEdge, curveT, vUR, vUL);
   i1 = eT->GetStartPoint()->GetTag();
   i2 = eT->GetEndPoint()->GetTag();
-  m->add(new fourierEdge(m, eT, eT->GetTag(), m->vertexByTag(i1),
-			 m->vertexByTag(i2)));
+  m->add(new fourierEdge(m, eT, eT->GetTag(), m->getVertex(i1),
+			 m->getVertex(i2)));
   FM::TopoEdge* eL = new FM::TopoEdge(++tagEdge, curveL, vUL, vLL); 
   i1 = eL->GetStartPoint()->GetTag();
   i2 = eL->GetEndPoint()->GetTag();
-  m->add(new fourierEdge(m, eL, eL->GetTag(), m->vertexByTag(i1),
-			 m->vertexByTag(i2)));
+  m->add(new fourierEdge(m, eL, eL->GetTag(), m->getVertex(i1),
+			 m->getVertex(i2)));
   
-  FM::TopoFace* face = new FM::TopoFace(m->numFace() + 1, patch);
+  FM::TopoFace* face = new FM::TopoFace(m->getNumFaces() + 1, patch);
   face->AddEdge(eB); face->AddEdge(eR); 
   face->AddEdge(eT); face->AddEdge(eL);
   std::list<GEdge*> l_edges;
   for (int j = 0; j < face->GetNumEdges(); j++) {
     int tag = face->GetEdge(j)->GetTag(); 
-    l_edges.push_back(m->edgeByTag(tag));
+    l_edges.push_back(m->getEdge(tag));
   }
   m->add(new fourierFace(m, face, face->GetTag(), l_edges));
 }
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 9c96046e56a26e2a5537d01c91514a7d45aa4521..705061f9e7282fd7a6e2eaa54232e81aece7fcb5 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Geo.cpp,v 1.14 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GModelIO_Geo.cpp,v 1.15 2008-02-20 09:20:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -39,22 +39,23 @@ int GModel::readGEO(const std::string &name)
 
 void GModel::createGEOInternals()
 {
-  geo_internals = new GEO_Internals;
+  _geo_internals = new GEO_Internals;
 }
 
 void GModel::deleteGEOInternals()
 {
-  delete geo_internals;
+  delete _geo_internals;
+  _geo_internals = 0;
 }
 
 int GModel::importGEOInternals()
 {
-  if(Tree_Nbr(geo_internals->Points)) {
-    List_T *points = Tree2List(geo_internals->Points);
+  if(Tree_Nbr(_geo_internals->Points)) {
+    List_T *points = Tree2List(_geo_internals->Points);
     for(int i = 0; i < List_Nbr(points); i++){
       Vertex *p;
       List_Read(points, i, &p);
-      GVertex *v = vertexByTag(p->Num);
+      GVertex *v = getVertex(p->Num);
       if(!v){
 	v = new gmshVertex(this, p);
 	add(v);
@@ -63,17 +64,17 @@ int GModel::importGEOInternals()
     }
     List_Delete(points);
   }
-  if(Tree_Nbr(geo_internals->Curves)) {
-    List_T *curves = Tree2List(geo_internals->Curves);
+  if(Tree_Nbr(_geo_internals->Curves)) {
+    List_T *curves = Tree2List(_geo_internals->Curves);
     for(int i = 0; i < List_Nbr(curves); i++){
       Curve *c;
       List_Read(curves, i, &c);
       if(c->Num >= 0 && c->beg && c->end){
-	GEdge *e = edgeByTag(c->Num);
+	GEdge *e = getEdge(c->Num);
 	if(!e){
 	  e = new gmshEdge(this, c,
-			   vertexByTag(c->beg->Num),
-			   vertexByTag(c->end->Num));
+			   getVertex(c->beg->Num),
+			   getVertex(c->end->Num));
 	  add(e);
 	}
 	else
@@ -84,12 +85,12 @@ int GModel::importGEOInternals()
     }
     List_Delete(curves);
   }
-  if(Tree_Nbr(geo_internals->Surfaces)) {
-    List_T *surfaces = Tree2List(geo_internals->Surfaces);
+  if(Tree_Nbr(_geo_internals->Surfaces)) {
+    List_T *surfaces = Tree2List(_geo_internals->Surfaces);
     for(int i = 0; i < List_Nbr(surfaces); i++){
       Surface *s;
       List_Read(surfaces, i, &s);
-      GFace *f = faceByTag(s->Num);
+      GFace *f = getFace(s->Num);
       if(!f){
 	f = new gmshFace(this, s);
 	add(f);
@@ -101,12 +102,12 @@ int GModel::importGEOInternals()
     }
     List_Delete(surfaces);
   } 
-  if(Tree_Nbr(geo_internals->Volumes)) {
-    List_T *volumes = Tree2List(geo_internals->Volumes);
+  if(Tree_Nbr(_geo_internals->Volumes)) {
+    List_T *volumes = Tree2List(_geo_internals->Volumes);
     for(int i = 0; i < List_Nbr(volumes); i++){
       Volume *v;
       List_Read(volumes, i, &v);
-      GRegion *r = regionByTag(v->Num);
+      GRegion *r = getRegion(v->Num);
       if(!r){
 	r = new gmshRegion(this, v);
 	add(r);
@@ -118,18 +119,18 @@ int GModel::importGEOInternals()
     }
     List_Delete(volumes);
   }
-  for(int i = 0; i < List_Nbr(geo_internals->PhysicalGroups); i++){
+  for(int i = 0; i < List_Nbr(_geo_internals->PhysicalGroups); i++){
     PhysicalGroup *p;
-    List_Read(geo_internals->PhysicalGroups, i, &p);
+    List_Read(_geo_internals->PhysicalGroups, i, &p);
     for(int j = 0; j < List_Nbr(p->Entities); j++){
       int num;
       List_Read(p->Entities, j, &num);
       GEntity *ge = 0;
       switch(p->Typ){
-      case MSH_PHYSICAL_POINT:   ge = vertexByTag(abs(num)); break;
-      case MSH_PHYSICAL_LINE:    ge = edgeByTag(abs(num)); break;
-      case MSH_PHYSICAL_SURFACE: ge = faceByTag(abs(num)); break;
-      case MSH_PHYSICAL_VOLUME:  ge = regionByTag(abs(num)); break;
+      case MSH_PHYSICAL_POINT:   ge = getVertex(abs(num)); break;
+      case MSH_PHYSICAL_LINE:    ge = getEdge(abs(num)); break;
+      case MSH_PHYSICAL_SURFACE: ge = getFace(abs(num)); break;
+      case MSH_PHYSICAL_VOLUME:  ge = getRegion(abs(num)); break;
       }
       int pnum = sign(num) * p->Num;
       if(ge && std::find(ge->physicals.begin(), ge->physicals.end(), pnum) == 
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index d4a4ec0b9bbf42a21496c27450d01c6eaf8a9df4..4966ab7fb301a3d1ca4fcba515d6a3085e68e3b4 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.32 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.33 2008-02-20 09:20:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -66,7 +66,7 @@ static void storeElementsInEntities(GModel *m,
     switch(numEdges){
     case 1: 
       {
-	GEdge *e = m->edgeByTag(it->first);
+	GEdge *e = m->getEdge(it->first);
 	if(!e){
 	  e = new discreteEdge(m, it->first);
 	  m->add(e);
@@ -76,7 +76,7 @@ static void storeElementsInEntities(GModel *m,
       break;
     case 3: case 4: 
       {
-	GFace *f = m->faceByTag(it->first);
+	GFace *f = m->getFace(it->first);
 	if(!f){
 	  f = new discreteFace(m, it->first);
 	  m->add(f);
@@ -87,7 +87,7 @@ static void storeElementsInEntities(GModel *m,
       break;
     case 6: case 12: case 9: case 8:
       {
-	GRegion *r = m->regionByTag(it->first);
+	GRegion *r = m->getRegion(it->first);
 	if(!r){
 	  r = new discreteRegion(m, it->first);
 	  m->add(r);
@@ -136,10 +136,10 @@ static void storePhysicalTagsInEntities(GModel *m, int dim,
   for(; it != map.end(); ++it){
     GEntity *ge = 0;
     switch(dim){
-    case 0: ge = m->vertexByTag(it->first); break;
-    case 1: ge = m->edgeByTag(it->first); break;
-    case 2: ge = m->faceByTag(it->first); break;
-    case 3: ge = m->regionByTag(it->first); break;
+    case 0: ge = m->getVertex(it->first); break;
+    case 1: ge = m->getEdge(it->first); break;
+    case 2: ge = m->getFace(it->first); break;
+    case 3: ge = m->getRegion(it->first); break;
     }
     if(ge){
       std::map<int, std::string>::const_iterator it2 = it->second.begin();
@@ -471,7 +471,12 @@ int GModel::readMSH(const std::string &name)
 
     }
     else if(!strncmp(&str[1], "NodeData", 8)) {
-      // there's some post-processing data to read later on
+      // there's some post-processing data to read later on, so cache
+      // the vertex indexing data
+      if(vertexVector.size())
+	_vertexVectorCache = vertexVector;
+      else
+	_vertexMapCache = vertexMap;
       postpro = true;
     }
 
@@ -506,7 +511,7 @@ int GModel::readMSH(const std::string &name)
   // treat points separately
   for(std::map<int, std::vector<MVertex*> >::iterator it = points.begin(); 
       it != points.end(); ++it){
-    GVertex *v = vertexByTag(it->first);
+    GVertex *v = getVertex(it->first);
     if(!v){
       v = new discreteVertex(this, it->first);
       add(v);
@@ -919,7 +924,7 @@ int GModel::readSTL(const std::string &name, double tolerance)
   Msg(INFO, "%d facets", points.size() / 3);
 
   // create face
-  GFace *face = new discreteFace(this, numFace() + 1);
+  GFace *face = new discreteFace(this, getNumFaces() + 1);
   add(face);
 
   // create (unique) vertices and triangles
@@ -2008,7 +2013,7 @@ int GModel::readP3D(const std::string &name)
 
   for(int n = 0; n < numBlocks; n++){
     if(Nk[n] == 1){
-      GFace *gf = new discreteFace(this, numFace() + 1);
+      GFace *gf = new discreteFace(this, getNumFaces() + 1);
       add(gf);
       gf->transfinite_vertices.resize(Ni[n]);
       for(int i = 0; i < Ni[n]; i++)
@@ -2041,7 +2046,7 @@ int GModel::readP3D(const std::string &name)
 			     gf->transfinite_vertices[i    ][j + 1]));
     }
     else{
-      GRegion *gr = new discreteRegion(this, numRegion() + 1);
+      GRegion *gr = new discreteRegion(this, getNumRegions() + 1);
       add(gr);
       gr->transfinite_vertices.resize(Ni[n]);
       for(int i = 0; i < Ni[n]; i++){
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 68bfc6853e372f415b98a656b65128535654c574..991554f92cc41f4019d314d3765bbe8d77b0918a 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.26 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GModelIO_OCC.cpp,v 1.27 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -410,8 +410,8 @@ void OCC_Internals::buildGModel(GModel *model)
     TopoDS_Edge edge = TopoDS::Edge(emap(i));
     int i1 = vmap.FindIndex(TopExp::FirstVertex(edge)); 
     int i2 = vmap.FindIndex(TopExp::LastVertex(edge));
-    GVertex *v1 = model->vertexByTag(i1);
-    GVertex *v2 = model->vertexByTag(i2);
+    GVertex *v1 = model->getVertex(i1);
+    GVertex *v2 = model->getVertex(i2);
     OCCEdge *e = new OCCEdge(model, edge, i, v1, v2);
     model->add(e);
   }
@@ -433,28 +433,28 @@ void OCC_Internals::buildGModel(GModel *model)
 
 int GModel::readOCCSTEP(const std::string &fn)
 {
-  occ_internals = new OCC_Internals;
-  occ_internals->loadSTEP(fn.c_str());
-  occ_internals->buildLists();
-  occ_internals->buildGModel(this);
+  _occ_internals = new OCC_Internals;
+  _occ_internals->loadSTEP(fn.c_str());
+  _occ_internals->buildLists();
+  _occ_internals->buildGModel(this);
   return 1;
 }
 
 int GModel::readOCCIGES(const std::string &fn)
 {
-  occ_internals = new OCC_Internals;
-  occ_internals->loadIGES(fn.c_str());
-  occ_internals->buildLists();
-  occ_internals->buildGModel(this);
+  _occ_internals = new OCC_Internals;
+  _occ_internals->loadIGES(fn.c_str());
+  _occ_internals->buildLists();
+  _occ_internals->buildGModel(this);
   return 1;
 }
 
 int GModel::readOCCBREP(const std::string &fn)
 {
-  occ_internals = new OCC_Internals;
-  occ_internals->loadBREP(fn.c_str());
-  occ_internals->buildLists();
-  occ_internals->buildGModel(this);
+  _occ_internals = new OCC_Internals;
+  _occ_internals->loadBREP(fn.c_str());
+  _occ_internals->buildLists();
+  _occ_internals->buildGModel(this);
   snapGVertices();
   return 1;
 }
@@ -493,10 +493,10 @@ void GModel::snapGVertices (void)
 
 void GModel::deleteOCCInternals()
 {
-  if(occ_internals) delete occ_internals;
+  if(_occ_internals) delete _occ_internals;
+  _occ_internals = 0;
 }
 
-
 /*
   OCC Creation routines
 */
@@ -533,9 +533,10 @@ void AddSimpleShapes(TopoDS_Shape theShape, TopTools_ListOfShape& theList)
   }
 }
 
-void OCC_Internals::applyBooleanOperator ( TopoDS_Shape tool ,  const BooleanOperator & op ){
-  if (tool.IsNull())return;
-  if (shape.IsNull())shape = tool;
+void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool,  const BooleanOperator & op)
+{
+  if (tool.IsNull()) return;
+  if (shape.IsNull()) shape = tool;
   else{
     switch(op){
     case OCC_Internals::Add :
@@ -602,16 +603,17 @@ void OCC_Internals::applyBooleanOperator ( TopoDS_Shape tool ,  const BooleanOpe
   }
 }
   
-void OCC_Internals::Sphere  ( const SPoint3 & center, const double & radius, const BooleanOperator & op ){
+void OCC_Internals::Sphere(const SPoint3 & center, const double & radius,
+			   const BooleanOperator & op)
+{
   // build a sphere
-  gp_Pnt aP (center.x(), center.y(), center.z());  
+  gp_Pnt aP(center.x(), center.y(), center.z());  
   TopoDS_Shape aShape = BRepPrimAPI_MakeSphere(aP, radius).Shape(); 
   // either add it to the current shape, or use it as a tool and remove the
   // sphere from the current shape
-  applyBooleanOperator ( aShape , op );
+  applyBooleanOperator(aShape,  op);
 }
 
-
 #else
 
 int GModel::readOCCSTEP(const std::string &fn)
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 05ea81c7b6ab3163189f81e1696dcb15f2b5aae8..6cf98c00e97be5f26920dc39197ec160cb7c6a87 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.102 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: Geo.cpp,v 1.103 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -3209,7 +3209,7 @@ void setVolumeSurfaces(Volume *v, List_T *loops)
 	  List_Add(v->SurfacesOrientations, &tmp);
 	}
 	else{
-	  GFace *gf = GModel::current()->faceByTag(abs(is));
+	  GFace *gf = GModel::current()->getFace(abs(is));
 	  if(gf) {
 	    List_Add(v->SurfacesByTag, &is);
 	  }
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index c7a20cfc3bfc60f7f9b56a3806eeb6e01d2a52ce..39ccda4aec5b0f4175bbae1462ea2d793652de4b 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -42,10 +42,11 @@ class MVertex{
   char _visible, _order;
   double _x, _y, _z;
   GEntity *_ge;
+  void *_data;
 
  public :
   MVertex(double x, double y, double z, GEntity *ge=0, int num=0) 
-    : _visible(true), _order(1), _x(x), _y(y), _z(z), _ge(ge)
+    : _visible(true), _order(1), _x(x), _y(y), _z(z), _ge(ge), _data(0)
   {
     if(num){
       _num = num;
@@ -57,7 +58,8 @@ class MVertex{
   }
   virtual ~MVertex(){}
 
-  // reset the global node number
+  // get/reset the global node number
+  static int getGlobalNumber(){ return _globalNum; }
   static void resetGlobalNumber(){ _globalNum = 0; }
 
   // get/set the visibility flag
@@ -85,6 +87,10 @@ class MVertex{
   inline int getNum() const { return _num; }
   inline void setNum(int num) { _num = num; }
 
+  // get/set the data
+  inline void *getData() { return _data; }
+  inline void setData(void *data) { _data = data; }
+
   // get/set ith parameter
   virtual bool getParameter(int i, double &par) const{ return false; }
   virtual bool setParameter(int i, double par){ return false; }
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index c67f2fe13a349fb1257376de25d8c2a59cbc8237..d40b8375cb67114c7d33a71c4a438d3e5c21bbc5 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.34 2008-02-17 09:30:28 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.35 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -52,7 +52,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
     for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()){	  
       TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
       int index = emap.FindIndex(edge);
-      GEdge *e = m->edgeByTag(index);
+      GEdge *e = m->getEdge(index);
       if(!e) throw;
       l_wire.push_back(e);
       l_oris.push_back(edge.Orientation());
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 69d58e11143795e5168112dbeee64dc191693d80..cefa15c4895a1f4948630a72a378cd3ed637ef2a 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCRegion.cpp,v 1.8 2008-02-17 08:47:59 geuzaine Exp $
+// $Id: OCCRegion.cpp,v 1.9 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -38,7 +38,7 @@ OCCRegion::OCCRegion(GModel *m, TopoDS_Solid _s, int num, TopTools_IndexedMapOfS
     for(exp3.Init(shell, TopAbs_FACE); exp3.More(); exp3.Next()){	  
       TopoDS_Face face = TopoDS::Face(exp3.Current());
       int index = fmap.FindIndex(face);
-      GFace *f = m->faceByTag(index);
+      GFace *f = m->getFace(index);
       if(!f) throw;
       l_faces.push_back(f);
       f->addRegion(this);
diff --git a/Geo/findLinks.cpp b/Geo/findLinks.cpp
index 71df7faa51d0dd07897c1fa7d9ada76cf2f59404..b9272d98fc209ce2e7c91c772137fd535175b8b1 100644
--- a/Geo/findLinks.cpp
+++ b/Geo/findLinks.cpp
@@ -1,4 +1,4 @@
-// $Id: findLinks.cpp,v 1.4 2008-02-17 08:47:59 geuzaine Exp $
+// $Id: findLinks.cpp,v 1.5 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -48,7 +48,7 @@ static int complink(const void *a, const void *b)
 static void recurFindLinkedEdges(int ed, List_T *edges, Tree_T *points,
 				 Tree_T *links)
 {
-  GEdge *ge = GModel::current()->edgeByTag(ed);
+  GEdge *ge = GModel::current()->getEdge(ed);
   if(!ge){
     Msg(GERROR, "Unknown curve %d", ed);
     return;
@@ -123,7 +123,7 @@ static void orientAndSortEdges(List_T *edges, Tree_T *links)
   List_Read(temp, 0, &num);
   List_Add(edges, &num);
 
-  GEdge *ge0 = GModel::current()->edgeByTag(abs(num));
+  GEdge *ge0 = GModel::current()->getEdge(abs(num));
   if(!ge0){
     Msg(GERROR, "Unknown curve %d", abs(num));
     return;
@@ -141,7 +141,7 @@ static void orientAndSortEdges(List_T *edges, Tree_T *links)
       nxa na;
       List_Read(lk.l, j, &na);
       if(ge0->tag() != na.a && List_Search(temp, &na.a, fcmp_absint)){
-	GEdge *ge1 = GModel::current()->edgeByTag(abs(na.a));
+	GEdge *ge1 = GModel::current()->getEdge(abs(na.a));
 	if(!ge1){
 	  Msg(GERROR, "Unknown curve %d", abs(na.a));
 	  return;
@@ -176,7 +176,7 @@ int allEdgesLinked(int ed, List_T *edges)
   for(int i = 0; i < List_Nbr(edges); i++){
     int num;
     List_Read(edges, i, &num);
-    GEdge *ge = GModel::current()->edgeByTag(abs(num));
+    GEdge *ge = GModel::current()->getEdge(abs(num));
     if(!ge){
       Msg(GERROR, "Unknown curve %d", abs(num));
       return 0;
@@ -220,7 +220,7 @@ int allEdgesLinked(int ed, List_T *edges)
 static void recurFindLinkedFaces(int fac, List_T *faces, Tree_T *edges, 
 				 Tree_T *links)
 {
-  GFace *gf = GModel::current()->faceByTag(abs(fac));
+  GFace *gf = GModel::current()->getFace(abs(fac));
   if(!gf){
     Msg(GERROR, "Unknown surface %d", abs(fac));
     return;
@@ -288,7 +288,7 @@ int allFacesLinked(int fac, List_T *faces)
   for(int i = 0; i < List_Nbr(faces); i++){
     int num;
     List_Read(faces, i, &num);
-    GFace *gf = GModel::current()->faceByTag(abs(num));
+    GFace *gf = GModel::current()->getFace(abs(num));
     if(!gf){
       Msg(GERROR, "Unknown surface %d", abs(num));
       return 0;
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 4d1f92b6db36e8a6d53f88518c5b6ab80f785740..6ccc91439404f8737e18254fdfef4d8c40cea72b 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.48 2008-02-17 09:30:28 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.49 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -32,7 +32,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
   for(int i = 0 ; i < List_Nbr(s->Generatrices); i++){
     Curve *c;
     List_Read(s->Generatrices, i, &c);
-    GEdge *e = m->edgeByTag(abs(c->Num));
+    GEdge *e = m->getEdge(abs(c->Num));
     if(e){
       l_edges.push_back(e);
       e->addFace(this);
@@ -59,7 +59,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
     for(int i = 0 ; i < List_Nbr(s->EmbeddedCurves); i++){
       Curve *c;
       List_Read(s->EmbeddedCurves, i, &c);
-      GEdge *e = m->edgeByTag(abs(c->Num));
+      GEdge *e = m->getEdge(abs(c->Num));
       if(e)
 	embedded_edges.push_back(e);
       else
@@ -70,7 +70,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
     for(int i = 0 ; i < List_Nbr(s->EmbeddedPoints); i++){
       Vertex *v;
       List_Read(s->EmbeddedPoints, i, &v);
-      GVertex *gv = m->vertexByTag(v->Num);
+      GVertex *gv = m->getVertex(v->Num);
       if(gv)
 	embedded_vertices.push_back(gv);
       else
@@ -110,7 +110,7 @@ void gmshFace::resetMeshAttributes()
     for(int i = 0; i < List_Nbr(s->TrsfPoints); i++){
       Vertex *corn;
       List_Read(s->TrsfPoints, i, &corn);
-      GVertex *gv = model()->vertexByTag(corn->Num);
+      GVertex *gv = model()->getVertex(corn->Num);
       if(gv)
 	meshAttributes.corners.push_back(gv);
       else
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index b9ab5e53b62850d9fa5202068baa401de7330169..55ef08ac572a051771fa3738427381ea65da5d66 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshRegion.cpp,v 1.18 2008-02-17 08:47:59 geuzaine Exp $
+// $Id: gmshRegion.cpp,v 1.19 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -32,7 +32,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
     List_Read(v->Surfaces, i, &s);
     int ori;
     List_Read(v->SurfacesOrientations, i, &ori);
-    GFace *f = m->faceByTag(abs(s->Num));
+    GFace *f = m->getFace(abs(s->Num));
     if(f){
       l_faces.push_back(f);
       l_dirs.push_back(ori);
@@ -43,7 +43,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
   for(int i = 0; i < List_Nbr(v->SurfacesByTag); i++){
     int is;
     List_Read(v->SurfacesByTag, i, &is);
-    GFace *f = m->faceByTag(abs(is));
+    GFace *f = m->getFace(abs(is));
     if(f){
       l_faces.push_back(f);
       l_dirs.push_back(sign(is));
@@ -64,7 +64,7 @@ void gmshRegion::resetMeshAttributes()
     for(int i = 0; i < List_Nbr(v->TrsfPoints); i++){
       Vertex *corn;
       List_Read(v->TrsfPoints, i, &corn);
-      GVertex *gv = model()->vertexByTag(corn->Num);
+      GVertex *gv = model()->getVertex(corn->Num);
       if(gv)
 	meshAttributes.corners.push_back(gv);
       else
diff --git a/Geo/gmshSurface.cpp b/Geo/gmshSurface.cpp
index 5ef8640f7b0e96a1473461849e8e0eb8aa640f08..b42da18673d025904c7d2c32f282e5d979bbaf27 100644
--- a/Geo/gmshSurface.cpp
+++ b/Geo/gmshSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshSurface.cpp,v 1.10 2008-02-17 08:47:59 geuzaine Exp $
+// $Id: gmshSurface.cpp,v 1.11 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -27,7 +27,7 @@
 
 std::map<int,gmshSurface*> gmshSurface::allGmshSurfaces;
 
-gmshSurface * gmshSphere::NewSphere(int iSphere, double x, double y, double z, double r)
+gmshSurface *gmshSphere::NewSphere(int iSphere, double x, double y, double z, double r)
 {
   gmshSphere *sph = new gmshSphere(x, y, z, r);
   
@@ -39,7 +39,7 @@ gmshSurface * gmshSphere::NewSphere(int iSphere, double x, double y, double z, d
   return sph;
 }
 
-gmshSurface * gmshSurface::surfaceByTag(int iSurface)
+gmshSurface *gmshSurface::getSurface(int iSurface)
 {
   std::map<int, gmshSurface*>::iterator it = allGmshSurfaces.find(iSurface);
   if(it == allGmshSurfaces.end()){
@@ -49,7 +49,7 @@ gmshSurface * gmshSurface::surfaceByTag(int iSurface)
   return it->second;
 }
 
-SPoint3  gmshSphere::point (double par1, double par2) const
+SPoint3 gmshSphere::point(double par1, double par2) const
 {
   par2 += M_PI*.5;
   const double x = xc + r * sin(par2) * cos(par1);
@@ -59,31 +59,33 @@ SPoint3  gmshSphere::point (double par1, double par2) const
   return SPoint3(x, y, z);
 }
 
-gmshSurface * gmshPolarSphere::NewPolarSphere(int iSphere, double x, double y, double z, double r)
+gmshSurface *gmshPolarSphere::NewPolarSphere(int iSphere, double x, double y, double z,
+					     double r)
 {
   gmshPolarSphere *sph = new gmshPolarSphere(x, y, z, r);
   
   if(allGmshSurfaces.find(iSphere) != allGmshSurfaces.end()){
-    Msg(GERROR,"gmshSurface %d already exists",iSphere);
+    Msg(GERROR, "gmshSurface %d already exists", iSphere);
   }
-  
+
   allGmshSurfaces[iSphere] = sph;
   return sph;
 }
-SPoint3  gmshPolarSphere::point (double parA, double parB) const
+
+SPoint3 gmshPolarSphere::point(double parA, double parB) const
 {
-	//stereographic projection from the south pole, origin of the axis at the center of the sphere
-	//parA=2rx/(r+z)
-	//parB=2ry/(r+z)
-	double rp2=parA*parA+parB*parB;
-	double z=r*(4*r*r-rp2)/(4*r*r+rp2);
-  return SPoint3((r+z)*parA/(2*r),(r+z)*parB/(2*r),z);
+  //stereographic projection from the south pole, origin of the axis
+  //at the center of the sphere 
+  //parA=2rx/(r+z) parB=2ry/(r+z)
+  double rp2 = parA * parA + parB * parB;
+  double z = r * (4 * r * r - rp2) / (4 * r * r + rp2);
+  return SPoint3((r + z) * parA / (2 * r), (r + z) * parB / (2 * r), z);
 }
 
-
-gmshSurface * gmshParametricSurface::NewParametricSurface(int iSurf,  char *valX,  char *valY,  char *valZ)
+gmshSurface *gmshParametricSurface::NewParametricSurface(int iSurf, char *valX,
+							 char *valY, char *valZ)
 {
-  gmshParametricSurface *sph = new gmshParametricSurface(valX,valY,valZ);
+  gmshParametricSurface *sph = new gmshParametricSurface(valX, valY, valZ);
   
   if(allGmshSurfaces.find(iSurf) != allGmshSurfaces.end()){
     Msg(GERROR,"gmshSurface %d already exists",iSurf);
@@ -92,7 +94,7 @@ gmshSurface * gmshParametricSurface::NewParametricSurface(int iSurf,  char *valX
   return sph;
 }
 
-gmshParametricSurface ::gmshParametricSurface(char *valX, char *valY, char *valZ)
+gmshParametricSurface::gmshParametricSurface(char *valX, char *valY, char *valZ)
 {
 #if !defined(HAVE_MATH_EVAL)
   Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
@@ -103,7 +105,7 @@ gmshParametricSurface ::gmshParametricSurface(char *valX, char *valY, char *valZ
 #endif
 }
 
-gmshParametricSurface ::~gmshParametricSurface()
+gmshParametricSurface::~gmshParametricSurface()
 {
 #if !defined(HAVE_MATH_EVAL)
   Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
@@ -113,17 +115,18 @@ gmshParametricSurface ::~gmshParametricSurface()
   evaluator_destroy(evalZ);
 #endif
 }
-SPoint3 gmshParametricSurface ::point(double par1, double par2) const
+
+SPoint3 gmshParametricSurface::point(double par1, double par2) const
 {
 #if !defined(HAVE_MATH_EVAL)
   Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
   return SPoint3(0.,0.,0.);
 #else
-  char *names[2] = {"u","v"};
+  char *names[2] = {"u", "v"};
   double values [2] = {par1,par2};
   const double x = evaluator_evaluate(evalX, 2, names, values);
   const double y = evaluator_evaluate(evalY, 2, names, values);
   const double z = evaluator_evaluate(evalZ, 2, names, values);
-  return SPoint3(x,y,z);
+  return SPoint3(x, y, z);
 #endif
 }
diff --git a/Geo/gmshSurface.h b/Geo/gmshSurface.h
index 1c0c421cd567a4d87db0453c08e0036d3d2fe735..20ee94528d5b4c35fd82cb27a722a3e5f4edb8b8 100644
--- a/Geo/gmshSurface.h
+++ b/Geo/gmshSurface.h
@@ -36,14 +36,14 @@ protected:
   static std::map<int, gmshSurface*> allGmshSurfaces;
 public:
   virtual ~gmshSurface(){}
-  static void reset () 
+  static void reset() 
   {
-    std::map<int,gmshSurface*>::iterator it = allGmshSurfaces.begin();
+    std::map<int, gmshSurface*>::iterator it = allGmshSurfaces.begin();
     for (; it != allGmshSurfaces.end(); ++it)
       delete it->second;
     allGmshSurfaces.clear();
   };
-  static gmshSurface* surfaceByTag(int tag);
+  static gmshSurface* getSurface(int tag);
   virtual Range<double> parBounds(int i) const = 0;
   /// Underlying geometric representation of this entity.
   enum gmshSurfaceType {
@@ -64,7 +64,7 @@ public:
   virtual SVector3 normal(const SPoint2 &param) const = 0;
   // Return the first derivate of the face at the parameter location.
   virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const = 0;
-  virtual double getMetricEigenvalue ( const SPoint2 &) {throw;}
+  virtual double getMetricEigenvalue(const SPoint2 &) { throw; }
 };
 
 class gmshSphere : public gmshSurface
@@ -146,7 +146,6 @@ public:
 
 };
 
-
 class gmshParametricSurface : public gmshSurface
 {
   void *evalX, *evalY, *evalZ;
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 248d46724b8011de1963af86b53b0d1354cf6f40..12213e09d1d8992a546b80a9f295d971d58a4f04 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.152 2008-02-17 08:47:59 geuzaine Exp $
+// $Id: Geom.cpp,v 1.153 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -467,8 +467,8 @@ void Draw_Geom()
   for(int i = 0; i < 6; i++)
     glDisable((GLenum)(GL_CLIP_PLANE0 + i));
 
-  bool geometryExists = 
-    m->numVertex() || m->numEdge() || m->numFace() || m->numRegion();
+  bool geometryExists = m->getNumVertices() || m->getNumEdges() || 
+    m->getNumFaces() || m->getNumRegions();
 
   if(geometryExists && (CTX.draw_bbox || !CTX.mesh.draw)) {
     glColor4ubv((GLubyte *) & CTX.color.fg);
diff --git a/Graphics/SelectBuffer.cpp b/Graphics/SelectBuffer.cpp
index b8020e3a17645643d1b6f3976e1e85cfd2ae9013..ab1a2642a18ac067d2b556d9e9af21bf04e51edc 100644
--- a/Graphics/SelectBuffer.cpp
+++ b/Graphics/SelectBuffer.cpp
@@ -1,4 +1,4 @@
-// $Id: SelectBuffer.cpp,v 1.18 2008-02-17 08:48:00 geuzaine Exp $
+// $Id: SelectBuffer.cpp,v 1.19 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -81,9 +81,9 @@ bool ProcessSelectionBuffer(int entityType,
   // In our case the selection buffer size is equal to between 5 and 7
   // times the maximum number of possible hits
   GModel *m = GModel::current();
-  int eles = (meshSelection && CTX.pick_elements) ? 4 * m->numElements() : 0;
-  int size = 7 * (m->numVertex() + m->numEdge() + m->numFace() + 
-		  m->numRegion() + eles) + 1000 ;
+  int eles = (meshSelection && CTX.pick_elements) ? 4 * m->getNumMeshElements() : 0;
+  int size = 7 * (m->getNumVertices() + m->getNumEdges() + m->getNumFaces() + 
+		  m->getNumRegions() + eles) + 1000 ;
 
   GLuint *selectionBuffer = new GLuint[size];
   glSelectBuffer(size, selectionBuffer);
@@ -171,7 +171,7 @@ bool ProcessSelectionBuffer(int entityType,
       switch (hits[i].type) {
       case 0:
 	{
-	  GVertex *v = m->vertexByTag(hits[i].ient);
+	  GVertex *v = m->getVertex(hits[i].ient);
 	  if(!v){
 	    Msg(GERROR, "Problem in point selection processing");
 	    return false;
@@ -182,7 +182,7 @@ bool ProcessSelectionBuffer(int entityType,
 	break;
       case 1:
 	{
-	  GEdge *e = m->edgeByTag(hits[i].ient);
+	  GEdge *e = m->getEdge(hits[i].ient);
 	  if(!e){
 	    Msg(GERROR, "Problem in line selection processing");
 	    return false;
@@ -197,7 +197,7 @@ bool ProcessSelectionBuffer(int entityType,
 	break;
       case 2:
 	{
-	  GFace *f = m->faceByTag(hits[i].ient);
+	  GFace *f = m->getFace(hits[i].ient);
 	  if(!f){
 	    Msg(GERROR, "Problem in surface selection processing");
 	    return false;
@@ -212,7 +212,7 @@ bool ProcessSelectionBuffer(int entityType,
 	break;
       case 3:
 	{
-	  GRegion *r = m->regionByTag(hits[i].ient);
+	  GRegion *r = m->getRegion(hits[i].ient);
 	  if(!r){
 	    Msg(GERROR, "Problem in volume selection processing");
 	    return false;
@@ -244,19 +244,19 @@ void HighlightEntityNum(int v, int c, int s, int r)
 {
   GModel *m = GModel::current();
   if(v) {
-    GVertex *pv = m->vertexByTag(v);
+    GVertex *pv = m->getVertex(v);
     if(pv) HighlightEntity(pv);
   }
   if(c) {
-    GEdge *pc = m->edgeByTag(c);
+    GEdge *pc = m->getEdge(c);
     if(pc) HighlightEntity(pc);
   }
   if(s) {
-    GFace *ps = m->faceByTag(s);
+    GFace *ps = m->getFace(s);
     if(ps) HighlightEntity(ps);
   }
   if(r) {
-    GRegion *pr = m->regionByTag(r);
+    GRegion *pr = m->getRegion(r);
     if(pr) HighlightEntity(pr);
   }
 }
@@ -278,19 +278,19 @@ void ZeroHighlightEntityNum(int v, int c, int s, int r)
 {
   GModel *m = GModel::current();
   if(v) {
-    GVertex *pv = m->vertexByTag(v);
+    GVertex *pv = m->getVertex(v);
     if(pv) ZeroHighlightEntity(pv);
   }
   if(c) {
-    GEdge *pc = m->edgeByTag(c);
+    GEdge *pc = m->getEdge(c);
     if(pc) ZeroHighlightEntity(pc);
   }
   if(s) {
-    GFace *ps = m->faceByTag(s);
+    GFace *ps = m->getFace(s);
     if(ps) ZeroHighlightEntity(ps);
   }
   if(r) {
-    GRegion *pr = m->regionByTag(r);
+    GRegion *pr = m->getRegion(r);
     if(pr) ZeroHighlightEntity(pr);
   }
 }
diff --git a/Mesh/BoundaryLayer.cpp b/Mesh/BoundaryLayer.cpp
index 5e97156a04d5b9408a8fd5e8d19e05648fe52862..401d58ffe62c9f1c4a3b808d507fc5cc6d1ea0ad 100644
--- a/Mesh/BoundaryLayer.cpp
+++ b/Mesh/BoundaryLayer.cpp
@@ -1,4 +1,4 @@
-// $Id: BoundaryLayer.cpp,v 1.8 2008-02-17 08:48:00 geuzaine Exp $
+// $Id: BoundaryLayer.cpp,v 1.9 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -56,7 +56,7 @@ int Mesh2DWithBoundaryLayers(GModel *m)
     if(gf->geomType() == GEntity::BoundaryLayerSurface){
       ExtrudeParams *ep = gf->meshAttributes.extrude;
       if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
-	GFace *from = m->faceByTag(std::abs(ep->geo.Source));
+	GFace *from = m->getFace(std::abs(ep->geo.Source));
 	if(!from){
 	  Msg(GERROR, "Unknown source face %d for boundary layer", ep->geo.Source);
 	  return 0;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index d0131f9a8d09fd13746aacc8512cf02e30f1722a..6bf77a53598ffa046eac7bc747afbbdc314af2ca 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.135 2008-02-17 08:48:00 geuzaine Exp $
+// $Id: Generator.cpp,v 1.136 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -70,10 +70,10 @@ void GetStatistics(double stat[50], double quality[3][100])
 
   GModel *m = GModel::current();
 
-  stat[0] = m->numVertex();
-  stat[1] = m->numEdge();
-  stat[2] = m->numFace();
-  stat[3] = m->numRegion();
+  stat[0] = m->getNumVertices();
+  stat[1] = m->getNumEdges();
+  stat[2] = m->getNumFaces();
+  stat[3] = m->getNumRegions();
   
   std::map<int, std::vector<GEntity*> > physicals[4];
   m->getPhysicalGroups(physicals);
@@ -147,7 +147,7 @@ void GetStatistics(double stat[50], double quality[3][100])
 
 bool TooManyElements(GModel *m, int dim)
 {
-  if(CTX.expert_mode || !m->numVertex()) return false;
+  if(CTX.expert_mode || !m->getNumVertices()) return false;
 
   // try to detect obvious mistakes in characteristic lenghts (one of
   // the most common cause for erroneous bug reports on the mailing
@@ -155,7 +155,7 @@ bool TooManyElements(GModel *m, int dim)
   double sumAllLc = 0.;
   for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it)
     sumAllLc += (*it)->prescribedMeshSizeAtVertex();
-  sumAllLc /= (double)m->numVertex();
+  sumAllLc /= (double)m->getNumVertices();
   if(!sumAllLc || pow(CTX.lc / sumAllLc, dim) > 1.e10) 
     return !GetBinaryAnswer("Your choice of characteristic lengths will likely produce\n"
 			    "a very large mesh. Do you really want to continue?\n\n"
@@ -418,7 +418,8 @@ void GenerateMesh(int ask)
     SetOrderN(m, CTX.mesh.order, CTX.mesh.second_order_linear, 
 	      CTX.mesh.second_order_incomplete);
 
-  Msg(INFO, "%d vertices %d elements", m->numVertices(), m->numElements());
+  Msg(INFO, "%d vertices %d elements",
+      m->getNumMeshVertices(), m->getNumMeshElements());
 
   CTX.threads_lock = 0;
   CTX.mesh.changed = ENT_ALL;
diff --git a/Mesh/meshGEdgeExtruded.cpp b/Mesh/meshGEdgeExtruded.cpp
index 57cde5047c13fcb0e3bb00ef0cf593c5df37a45a..fc073838a2794e4dbd1fdd0417dc1d774e249c78 100644
--- a/Mesh/meshGEdgeExtruded.cpp
+++ b/Mesh/meshGEdgeExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdgeExtruded.cpp,v 1.9 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGEdgeExtruded.cpp,v 1.10 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -76,7 +76,7 @@ int MeshExtrudedCurve(GEdge *ge)
   }
   else {
     // curve is a copy of another curve (the "top" of the extrusion)
-    GEdge *from = ge->model()->edgeByTag(std::abs(ep->geo.Source));
+    GEdge *from = ge->model()->getEdge(std::abs(ep->geo.Source));
     if(!from){
       Msg(GERROR, "Unknown source curve %d for extrusion", ep->geo.Source);
       return 0;
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index e244aaab4db2cc129cfa8ab0d6547c76e8342c0e..5a9d59a3c6b12a351f4b08f4d5204e9931490aac 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceExtruded.cpp,v 1.25 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.26 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -217,7 +217,7 @@ int MeshExtrudedSurface(GFace *gf,
   
   if(ep->geo.Mode == EXTRUDED_ENTITY) {
     // surface is extruded from a curve
-    GEdge *from = gf->model()->edgeByTag(std::abs(ep->geo.Source));
+    GEdge *from = gf->model()->getEdge(std::abs(ep->geo.Source));
     if(!from){
       Msg(GERROR, "Unknown source curve %d for extrusion", ep->geo.Source);
       return 0;
@@ -226,7 +226,7 @@ int MeshExtrudedSurface(GFace *gf,
   }
   else {
     // surface is a copy of another surface (the "top" of the extrusion)
-    GFace *from = gf->model()->faceByTag(std::abs(ep->geo.Source));
+    GFace *from = gf->model()->getFace(std::abs(ep->geo.Source));
     if(!from){ 
       Msg(GERROR, "Unknown source surface %d for extrusion", ep->geo.Source);
       return 0;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 06543f617de81c320d5b75e58e65fbfbcda9ef2a..e5164320be84c7fb59486e72ededc6fe0602f5c7 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.41 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.42 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -158,7 +158,7 @@ void TransferTetgenMesh(GRegion *gr,
     v[0] = numberedV[out.trifacelist[i * 3 + 0] - 1];
     v[1] = numberedV[out.trifacelist[i * 3 + 1] - 1];
     v[2] = numberedV[out.trifacelist[i * 3 + 2] - 1];
-    GFace *gf = gr->model()->faceByTag(out.trifacemarkerlist[i]);
+    GFace *gf = gr->model()->getFace(out.trifacemarkerlist[i]);
     for(int j = 0; j < 3; j++){	  
       if(v[j]->onWhat()->dim() == 3){
 	v[j]->onWhat()->mesh_vertices.erase
diff --git a/Mesh/meshGRegionCarveHole.cpp b/Mesh/meshGRegionCarveHole.cpp
index b1dfdeacb75f1dc819d59944840e8107978b3479..081e785453d8c809a399eb05e97c308ed8006e9d 100644
--- a/Mesh/meshGRegionCarveHole.cpp
+++ b/Mesh/meshGRegionCarveHole.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionCarveHole.cpp,v 1.4 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGRegionCarveHole.cpp,v 1.5 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -86,7 +86,7 @@ void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces
   // add all points from carving surfaces into kdtree
   int numnodes = 0;
   for(unsigned int i = 0; i < surfaces.size(); i++){
-    GFace *gf = m->faceByTag(surfaces[i]);
+    GFace *gf = m->getFace(surfaces[i]);
     if(!gf){
       Msg(GERROR, "Unknown carving surface %d", surfaces[i]);
       return;
@@ -97,7 +97,7 @@ void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces
   ANNpointArray kdnodes = annAllocPts(numnodes, 4);
   int k = 0;
   for(unsigned int i = 0; i < surfaces.size(); i++){
-    GFace *gf = m->faceByTag(surfaces[i]);
+    GFace *gf = m->getFace(surfaces[i]);
     for(unsigned int j = 0; j < gf->mesh_vertices.size(); j++){
       kdnodes[k][0] = gf->mesh_vertices[j]->x();
       kdnodes[k][1] = gf->mesh_vertices[j]->y();
@@ -122,7 +122,7 @@ void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces
   // intersections o see who's inside)
 
   // generate discrete boundary mesh of the carved hole
-  GFace *gf = m->faceByTag(num);
+  GFace *gf = m->getFace(num);
   if(!gf) return;
   std::set<MFace, Less_Face> faces;
   std::list<GFace*> f = gr->faces();
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 3828695df9d14feea2de3453f90cc8a417069694..9465ebb985417097ba5d384a087583d9095c5940 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionExtruded.cpp,v 1.21 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.22 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -212,7 +212,7 @@ void meshGRegionExtruded::operator() (GRegion *gr)
   insertAllVertices(gr, pos);
 
   // volume is extruded from a surface
-  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  GFace *from = gr->model()->getFace(std::abs(ep->geo.Source));
   if(!from){
     Msg(GERROR, "Unknown source surface %d for extrusion", ep->geo.Source);
     return;
@@ -258,7 +258,7 @@ void phase1(GRegion *gr,
 	    std::set<std::pair<MVertex*, MVertex*> > &edges)
 {
   ExtrudeParams *ep = gr->meshAttributes.extrude;
-  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  GFace *from = gr->model()->getFace(std::abs(ep->geo.Source));
   if(!from) return;
 
   for(unsigned int i = 0; i < from->triangles.size(); i++){
@@ -295,7 +295,7 @@ void phase2(GRegion *gr,
 	    int &swap)
 {
   ExtrudeParams *ep = gr->meshAttributes.extrude;
-  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  GFace *from = gr->model()->getFace(std::abs(ep->geo.Source));
   if(!from) return;
 
   for(unsigned int i = 0; i < from->triangles.size(); i++){
@@ -361,7 +361,7 @@ void phase3(GRegion *gr,
 	    std::set<std::pair<MVertex*, MVertex*> > &edges)
 {
   ExtrudeParams *ep = gr->meshAttributes.extrude;
-  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  GFace *from = gr->model()->getFace(std::abs(ep->geo.Source));
   if(!from) return;
 
   for(unsigned int i = 0; i < from->triangles.size(); i++){
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 5ec881ceab82c3e2d985aca99ed34fcb531ef10c..5d5d44d8ed04d07deee607b974e7f88ebac76e73 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -336,7 +336,7 @@
 /* Copy the first part of user declarations.  */
 #line 1 "Gmsh.y"
 
-// $Id: Gmsh.tab.cpp,v 1.346 2008-02-17 08:48:02 geuzaine Exp $
+// $Id: Gmsh.tab.cpp,v 1.347 2008-02-20 09:20:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -4819,7 +4819,7 @@ yyreduce:
         if(v)
           att->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
         else{
-          GVertex *gv = GModel::current()->vertexByTag((int)d);
+          GVertex *gv = GModel::current()->getVertex((int)d);
           if(gv) 
             att->addPoint(gv->x(), gv->y(), gv->z());
         }
@@ -4952,7 +4952,7 @@ yyreduce:
 	if(v)
 	  attractor->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag((int)d);
+	  GVertex *gv = GModel::current()->getVertex((int)d);
 	  if(gv) 
 	    attractor->addPoint(gv->x(), gv->y(), gv->z());
 	}
@@ -4990,7 +4990,7 @@ yyreduce:
 	  att->addCurve(c, (int)pars[3]);
 	}
 	else{
-	  GEdge *ge = GModel::current()->edgeByTag((int)d);
+	  GEdge *ge = GModel::current()->getEdge((int)d);
 	  if(ge){
 	    att->addGEdge(ge, (int)pars[3]);
 	  }
@@ -5013,7 +5013,7 @@ yyreduce:
 	if(v)
 	  v->lc = (yyvsp[(5) - (6)].d);
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag((int)d);
+	  GVertex *gv = GModel::current()->getVertex((int)d);
 	  if(gv) 
 	    gv->setPrescribedMeshSizeAtVertex((yyvsp[(5) - (6)].d));
 	}
@@ -5373,7 +5373,7 @@ yyreduce:
   case 111:
 #line 1526 "Gmsh.y"
     {
-      myGmshSurface = gmshSurface::surfaceByTag((int)(yyvsp[(3) - (4)].d));
+      myGmshSurface = gmshSurface::getSurface((int)(yyvsp[(3) - (4)].d));
       (yyval.s).Type = 0;
       (yyval.s).Num = 0;
     ;}
@@ -5642,7 +5642,7 @@ yyreduce:
 	  List_Add((yyval.l), &TheShape);
 	}
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag(TheShape.Num);
+	  GVertex *gv = GModel::current()->getVertex(TheShape.Num);
 	  if(gv){
 	    TheShape.Type = MSH_POINT_FROM_GMODEL;
 	    List_Add((yyval.l), &TheShape);
@@ -5668,7 +5668,7 @@ yyreduce:
 	  List_Add((yyval.l), &TheShape);
 	}
 	else{
-	  GEdge *ge = GModel::current()->edgeByTag(TheShape.Num);
+	  GEdge *ge = GModel::current()->getEdge(TheShape.Num);
 	  if(ge){
 	    TheShape.Type = MSH_SEGM_FROM_GMODEL;
 	    List_Add((yyval.l), &TheShape);
@@ -5694,7 +5694,7 @@ yyreduce:
 	  List_Add((yyval.l), &TheShape);
 	}
 	else{
-	  GFace *gf = GModel::current()->faceByTag(TheShape.Num);
+	  GFace *gf = GModel::current()->getFace(TheShape.Num);
 	  if(gf){
 	    TheShape.Type = MSH_SURF_FROM_GMODEL;
 	    List_Add((yyval.l), &TheShape);
@@ -5720,7 +5720,7 @@ yyreduce:
 	  List_Add((yyval.l), &TheShape);
 	}
 	else{
-	  GRegion *gr = GModel::current()->regionByTag(TheShape.Num);
+	  GRegion *gr = GModel::current()->getRegion(TheShape.Num);
 	  if(gr){
 	    TheShape.Type = MSH_VOLUME_FROM_GMODEL;
 	    List_Add((yyval.l), &TheShape);
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index fe26794cf60791f75f533fb6cd08c981766da38b..ee4c1c5644ec3411b3d40953c570641d07b0f93d 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1,5 +1,5 @@
 %{
-// $Id: Gmsh.y,v 1.298 2008-02-17 08:48:04 geuzaine Exp $
+// $Id: Gmsh.y,v 1.299 2008-02-20 09:20:46 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -1040,7 +1040,7 @@ Shape :
         if(v)
           att->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
         else{
-          GVertex *gv = GModel::current()->vertexByTag((int)d);
+          GVertex *gv = GModel::current()->getVertex((int)d);
           if(gv) 
             att->addPoint(gv->x(), gv->y(), gv->z());
         }
@@ -1150,7 +1150,7 @@ Shape :
 	if(v)
 	  attractor->addPoint(v->Pos.X, v->Pos.Y, v->Pos.Z);
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag((int)d);
+	  GVertex *gv = GModel::current()->getVertex((int)d);
 	  if(gv) 
 	    attractor->addPoint(gv->x(), gv->y(), gv->z());
 	}
@@ -1185,7 +1185,7 @@ Shape :
 	  att->addCurve(c, (int)pars[3]);
 	}
 	else{
-	  GEdge *ge = GModel::current()->edgeByTag((int)d);
+	  GEdge *ge = GModel::current()->getEdge((int)d);
 	  if(ge){
 	    att->addGEdge(ge, (int)pars[3]);
 	  }
@@ -1205,7 +1205,7 @@ Shape :
 	if(v)
 	  v->lc = $5;
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag((int)d);
+	  GVertex *gv = GModel::current()->getVertex((int)d);
 	  if(gv) 
 	    gv->setPrescribedMeshSizeAtVertex($5);
 	}
@@ -1524,7 +1524,7 @@ Shape :
     }  
   | tCoordinates tSurface FExpr tEND
     {
-      myGmshSurface = gmshSurface::surfaceByTag((int)$3);
+      myGmshSurface = gmshSurface::getSurface((int)$3);
       $$.Type = 0;
       $$.Num = 0;
     }  
@@ -1746,7 +1746,7 @@ ListOfShapes :
 	  List_Add($$, &TheShape);
 	}
 	else{
-	  GVertex *gv = GModel::current()->vertexByTag(TheShape.Num);
+	  GVertex *gv = GModel::current()->getVertex(TheShape.Num);
 	  if(gv){
 	    TheShape.Type = MSH_POINT_FROM_GMODEL;
 	    List_Add($$, &TheShape);
@@ -1769,7 +1769,7 @@ ListOfShapes :
 	  List_Add($$, &TheShape);
 	}
 	else{
-	  GEdge *ge = GModel::current()->edgeByTag(TheShape.Num);
+	  GEdge *ge = GModel::current()->getEdge(TheShape.Num);
 	  if(ge){
 	    TheShape.Type = MSH_SEGM_FROM_GMODEL;
 	    List_Add($$, &TheShape);
@@ -1792,7 +1792,7 @@ ListOfShapes :
 	  List_Add($$, &TheShape);
 	}
 	else{
-	  GFace *gf = GModel::current()->faceByTag(TheShape.Num);
+	  GFace *gf = GModel::current()->getFace(TheShape.Num);
 	  if(gf){
 	    TheShape.Type = MSH_SURF_FROM_GMODEL;
 	    List_Add($$, &TheShape);
@@ -1815,7 +1815,7 @@ ListOfShapes :
 	  List_Add($$, &TheShape);
 	}
 	else{
-	  GRegion *gr = GModel::current()->regionByTag(TheShape.Num);
+	  GRegion *gr = GModel::current()->getRegion(TheShape.Num);
 	  if(gr){
 	    TheShape.Type = MSH_VOLUME_FROM_GMODEL;
 	    List_Add($$, &TheShape);
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index cd649a769f6a04c565ea299f62b6fdd9ce514f32..85468a536f362873abe4e31c11dec5654733c059 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -852,7 +852,7 @@ int gmsh_yy_flex_debug = 0;
 char *gmsh_yytext;
 #line 1 "Gmsh.l"
 #line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.346 2008-02-17 08:48:05 geuzaine Exp $
+// $Id: Gmsh.yy.cpp,v 1.347 2008-02-20 09:20:47 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
diff --git a/Plugin/Triangulate.cpp b/Plugin/Triangulate.cpp
index ec34194da57dec8ac7b3bb7b929a8c098b48e909..7083b9c2d4cfe7822957092304a4affe0d940c29 100644
--- a/Plugin/Triangulate.cpp
+++ b/Plugin/Triangulate.cpp
@@ -1,4 +1,4 @@
-// $Id: Triangulate.cpp,v 1.44 2008-02-17 08:48:08 geuzaine Exp $
+// $Id: Triangulate.cpp,v 1.45 2008-02-20 09:20:47 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -104,7 +104,7 @@ static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
     double *p = (double *)List_Pointer_Fast(inList, i);
     points.push_back(new MVertex(p[0], p[1], p[2]));
   }
-  discreteFace *s = new discreteFace(GModel::current(), GModel::current()->numFace() + 1);
+  discreteFace *s = new discreteFace(GModel::current(), GModel::current()->getNumFaces() + 1);
   s->computeMeanPlane(points);
   double plan[3][3];
   s->getMeanPlaneData(plan);
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 90f3ca3ba46bb1c9ad3e822083dcfefbbbe1c47f..866982962356af3a8e37f4ae03ba15855e79ed33 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -37,7 +37,7 @@ class PViewDataGModel : public PViewData {
   double getMin(int step=-1){ return 0.; }
   double getMax(int step=-1){ return 1.; }
   SBoundingBox3d getBoundingBox(){ return SBoundingBox3d(); }
-  int getNumElements(int type=0){ return _model->numElements(); }
+  int getNumElements(int type=0){ return _model->getNumMeshElements(); }
   int getDimension(int ele){ return 0; }
   int getNumNodes(int ele){ return 0; }
   void getNode(int ele, int nod, double &x, double &y, double &z){}