diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index ad92be02e4bc74df115d69f5fc35cb0e15be226b..1d18a5b081d7ef663f180f294878169fd7e8b259 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -108,20 +108,6 @@ void GEdge::setVisibility(char val, bool recursive)
   }
 }
 
-void GEdge::recomputeMeshPartitions()
-{
-  for(unsigned int i = 0; i < lines.size(); i++) {
-    int part = lines[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-}
-
-void GEdge::deleteMeshPartitions()
-{
-  for(unsigned int i = 0; i < lines.size(); i++)
-    lines[i]->setPartition(0);
-}
-
 std::string GEdge::getAdditionalInfoString()
 {
   if(!v0 || !v1) return std::string("");
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 4c9a8978054e314ff26efa3ef54bab196bd20098..d9a01a4577a3d667db86a4ac833902f6fe10edab 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -49,7 +49,7 @@ class GEdge : public GEntity {
   // get the bounding box
   virtual SBoundingBox3d bounds() const;
 
-  // faces that bound this entity bounds
+  // faces that this entity bounds
   virtual std::list<GFace*> faces() const { return l_faces; }
 
   // get the parameter location for a point in space on the edge
@@ -76,13 +76,7 @@ class GEdge : public GEntity {
   virtual double curvature(double par) const;
 
   // reparmaterize the point onto the given face
-  virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const;
-
-  // recompute the mesh partitions defined on this edge
-  void recomputeMeshPartitions();
-
-  // delete the mesh partitions defined on this edge
-  void deleteMeshPartitions();
+  virtual SPoint2 reparamOnFace(GFace *face, double epar, int dir) const;
 
   // return the minimum number of segments used for meshing the edge
   virtual int minimumMeshSegments() const { return 1; }
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index b4adc573d9a6a0af03c0df7913d40615ba78d73b..582bb18b9a3218a0f4fe54c35d20a06c9659d0cb 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -4,7 +4,9 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include <stdio.h>
+#include "GModel.h"
 #include "GEntity.h"
+#include "MElement.h"
 
 #if defined(HAVE_GMSH_EMBEDDED)
 #  include "GmshEmbedded.h"
@@ -75,3 +77,18 @@ std::string GEntity::getInfoString()
 
   return out;
 }
+
+void GEntity::recomputeMeshPartitions()
+{
+  for(unsigned int i = 0; i < getNumMeshElements(); i++) {
+    int part = getMeshElement(i)->getPartition();
+    if(part) model()->getMeshPartitions().insert(part);
+  }
+}
+
+void GEntity::deleteMeshPartitions()
+{
+  for(unsigned int i = 0; i < getNumMeshElements(); i++)
+    getMeshElement(i)->setPartition(0);
+}
+
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 5476df2cc76743168e27147a15c5af73b4c5be5f..78634c68b3a41959228bec9a1c90d0454b248c47 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -245,6 +245,12 @@ class GEntity {
 
   // get the mesh vertex at the given index
   MVertex *getMeshVertex(unsigned int index) { return mesh_vertices[index]; }
+
+  // recompute the mesh partitions defined on this entity
+  void recomputeMeshPartitions();
+
+  // delete the mesh partitions defined on this entity
+  void deleteMeshPartitions();
 };
 
 class GEntityLessThan {
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 4b155a4ac4d9e8a5d377f9ca049df7cac2009a68..0ca129f65bab80b86dd5548e8361424c263aeaf8 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -150,26 +150,6 @@ void GFace::setVisibility(char val, bool recursive)
   }
 }
 
-void GFace::recomputeMeshPartitions()
-{
-  for(unsigned int i = 0; i < triangles.size(); i++) {
-    int part = triangles[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-  for(unsigned int i = 0; i < quadrangles.size(); i++) {
-    int part = quadrangles[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-}
-
-void GFace::deleteMeshPartitions()
-{
-  for(unsigned int i = 0; i < triangles.size(); i++)
-    triangles[i]->setPartition(0);
-  for(unsigned int i = 0; i < quadrangles.size(); i++)
-    quadrangles[i]->setPartition(0);
-}
-
 std::string GFace::getAdditionalInfoString()
 {
   if(l_edges.empty()) return std::string("");
diff --git a/Geo/GFace.h b/Geo/GFace.h
index f20f3d9f09223018fc26128ca8e8d3c03d3ad1c8..3f7948bea1725917882bcedfdfb67f587c365cfe 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -126,12 +126,6 @@ class GFace : public GEntity
   // return the curvature i.e. the divergence of the normal
   virtual double curvature(const SPoint2 &param) const;
 
-  // recompute the mesh partitions defined on this face
-  void recomputeMeshPartitions();
-
-  // delete the mesh partitions defined on this face
-  void deleteMeshPartitions();
-
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 48af9dad50e66c0678a0a27038cee1a56ce09015..c9965077cd1aa4203a740bf866672fcfce3f7a86 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -221,123 +221,74 @@ void GModel::snapVertices()
   }
 }
 
-std::vector<GEntity*> GModel::getEntities()
+void GModel::getEntities(std::vector<GEntity*> &entities)
 {
-  std::vector<GEntity*> entities;
+  entities.clear();
   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());
-  return entities;
 }
 
 bool GModel::noPhysicalGroups()
 {
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    if((*it)->physicals.size()) return false;
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    if((*it)->physicals.size()) return false;
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    if((*it)->physicals.size()) return false;
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    if((*it)->physicals.size()) return false;
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    if(entities[i]->physicals.size()) return false;
   return true;
 }
 
-static void addInGroup(GEntity *ge, std::map<int, std::vector<GEntity*> > &group)
-{
-  for(unsigned int i = 0; i < ge->physicals.size(); i++){
-    // physicals can be stored with negative signs when the entity
-    // should be "reversed"
-    int p = std::abs(ge->physicals[i]);
-    if(std::find(group[p].begin(), group[p].end(), ge) == group[p].end())
-      group[p].push_back(ge);
-  }
-}
-
 void GModel::getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4])
 {
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    addInGroup(*it, groups[0]);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    addInGroup(*it, groups[1]);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    addInGroup(*it, groups[2]);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    addInGroup(*it, groups[3]);
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++){
+    std::map<int, std::vector<GEntity*> > &group(groups[entities[i]->dim()]);
+    for(unsigned int j = 0; j < entities[i]->physicals.size(); j++){
+      // physicals can be stored with negative signs when the entity
+      // should be "reversed"
+      int p = std::abs(entities[i]->physicals[j]);
+      if(std::find(group[p].begin(), group[p].end(), entities[i]) == group[p].end())
+	group[p].push_back(entities[i]);
+    }
+  }
 }
 
 void GModel::deletePhysicalGroups()
 {
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    (*it)->physicals.clear();
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    (*it)->physicals.clear();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    (*it)->physicals.clear();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    (*it)->physicals.clear();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    entities[i]->physicals.clear();
 }
 
 void GModel::deletePhysicalGroup(int dim, int num)
 {
-  switch(dim){
-  case 0:
-    for(viter it = firstVertex(); it != lastVertex(); ++it){
-      std::vector<int> p;
-      for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-        if((*it)->physicals[i] != num) p.push_back((*it)->physicals[i]);
-      (*it)->physicals = p;
-    }
-    break;
-  case 1:
-    for(eiter it = firstEdge(); it != lastEdge(); ++it){
-      std::vector<int> p;
-      for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-        if((*it)->physicals[i] != num) p.push_back((*it)->physicals[i]);
-      (*it)->physicals = p;
-    }
-    break;
-  case 2:
-    for(fiter it = firstFace(); it != lastFace(); ++it){
-      std::vector<int> p;
-      for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-        if((*it)->physicals[i] != num) p.push_back((*it)->physicals[i]);
-      (*it)->physicals = p;
-    }
-    break;
-  case 3:
-    for(riter it = firstRegion(); it != lastRegion(); ++it){
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++){
+    if(dim == entities[i]->dim()){
       std::vector<int> p;
-      for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-        if((*it)->physicals[i] != num) p.push_back((*it)->physicals[i]);
-      (*it)->physicals = p;
+      for(unsigned int j = 0; j < entities[i]->physicals.size(); j++)
+        if(entities[i]->physicals[j] != num) 
+	  p.push_back(entities[i]->physicals[j]);
+      entities[i]->physicals = p;
     }
-    break;
   }
 }
 
 int GModel::getMaxPhysicalNumber()
 {
+  std::vector<GEntity*> entities;
+  getEntities(entities);
   int num = 0;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-      num = std::max(num, std::abs((*it)->physicals[i]));
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-      num = std::max(num, std::abs((*it)->physicals[i]));
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-      num = std::max(num, std::abs((*it)->physicals[i]));
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->physicals.size(); i++)
-      num = std::max(num, std::abs((*it)->physicals[i]));
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->physicals.size(); j++)
+      num = std::max(num, std::abs(entities[i]->physicals[j]));
   return num;
 }
 
-
-
-
 int GModel::setPhysicalName(std::string name, int number)
 {
   // check if the name is already used
@@ -359,21 +310,13 @@ std::string GModel::getPhysicalName(int number)
 
 SBoundingBox3d GModel::bounds()
 {
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  // using the mesh vertices for now; should use entities[i]->bounds() instead
   SBoundingBox3d bb;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    bb += (*it)->bounds();
-
-  // use the mesh instead of the geometry for now
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      bb += (*it)->mesh_vertices[i]->point();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      bb += (*it)->mesh_vertices[i]->point();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      bb += (*it)->mesh_vertices[i]->point();
-
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+      bb += entities[i]->mesh_vertices[j]->point();
   return bb;
 }
 
@@ -410,27 +353,21 @@ int GModel::getMeshStatus(bool countDiscrete)
 
 int GModel::getNumMeshVertices()
 {
-  int n = 0;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    n += (*it)->mesh_vertices.size();
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    n += (*it)->mesh_vertices.size();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    n += (*it)->mesh_vertices.size();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    n += (*it)->mesh_vertices.size();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  unsigned int n = 0;
+  for(unsigned int i = 0; i < entities.size(); i++)
+    n += entities[i]->mesh_vertices.size();
   return n;
 }
 
 int GModel::getNumMeshElements()
 {
+  std::vector<GEntity*> entities;
+  getEntities(entities);
   unsigned int n = 0;
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    n += (*it)->getNumMeshElements();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    n += (*it)->getNumMeshElements();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    n += (*it)->getNumMeshElements();
+  for(unsigned int i = 0; i < entities.size(); i++)
+    n += entities[i]->getNumMeshElements();
   return n;
 }
 
@@ -501,10 +438,11 @@ MElement *GModel::getMeshElementByCoord(SPoint3 &p)
     const int maxElePerBucket = 100; // memory vs. speed trade-off
     _octree = Octree_Create(maxElePerBucket, min, size, 
 			    MElementBB, MElementCentroid, MElementInEle);
-    std::vector<GEntity*> ent = getEntities();
-    for(unsigned int i = 0; i < ent.size(); i++)
-      for(unsigned int j = 0; j < ent[i]->getNumMeshElements(); j++)
-	Octree_Insert(ent[i]->getMeshElement(j), _octree);
+    std::vector<GEntity*> entities;
+    getEntities(entities);
+    for(unsigned int i = 0; i < entities.size(); i++)
+      for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+	Octree_Insert(entities[i]->getMeshElement(j), _octree);
     Octree_Arrange(_octree);
   }
   double P[3] = {p.x(), p.y(), p.z()};
@@ -515,13 +453,6 @@ MElement *GModel::getMeshElementByCoord(SPoint3 &p)
 #endif
 }
 
-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];
-}
-
 MVertex *GModel::getMeshVertexByTag(int n)
 {
   if(_vertexVectorCache.empty() && _vertexMapCache.empty()){
@@ -529,30 +460,24 @@ MVertex *GModel::getMeshVertexByTag(int n)
     _vertexVectorCache.clear();
     _vertexMapCache.clear();
     bool dense = (getNumMeshVertices() == MVertex::getGlobalNumber());
+    std::vector<GEntity*> entities;
+    getEntities(entities);
     if(dense){
       // numbering starts at 1
       _vertexVectorCache.resize(MVertex::getGlobalNumber() + 1);
-      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);
+      for(unsigned int i = 0; i < entities.size(); i++)
+	for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+	  _vertexVectorCache[entities[i]->mesh_vertices[j]->getNum()] = 
+	    entities[i]->mesh_vertices[j];
     }
     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);
+      for(unsigned int i = 0; i < entities.size(); i++)
+	for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+	  _vertexMapCache[entities[i]->mesh_vertices[j]->getNum()] = 
+	    entities[i]->mesh_vertices[j];
     }
   }
-
+  
   if(n < (int)_vertexVectorCache.size())
     return _vertexVectorCache[n];
   else
@@ -623,93 +548,37 @@ void GModel::removeInvisibleElements()
 
 int GModel::indexMeshVertices(bool all)
 {
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
   // tag all mesh vertices with -1 (negative vertices will not be
   // saved)
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->setIndex(-1);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->setIndex(-1);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->setIndex(-1);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->setIndex(-1);
-
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+      entities[i]->mesh_vertices[j]->setIndex(-1);
+  
   // tag all mesh vertices belonging to elements that need to be saved
   // with 0
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    if(all || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-        (*it)->mesh_vertices[i]->setIndex(0);
-    }
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    if(all || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->lines.size(); i++)
-        for(int j = 0; j < (*it)->lines[i]->getNumVertices(); j++)
-          (*it)->lines[i]->getVertex(j)->setIndex(0);
-    }
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    if(all || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-        for(int j = 0; j < (*it)->triangles[i]->getNumVertices(); j++)
-          (*it)->triangles[i]->getVertex(j)->setIndex(0);
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-        for(int j = 0; j < (*it)->quadrangles[i]->getNumVertices(); j++)
-          (*it)->quadrangles[i]->getVertex(j)->setIndex(0);
-    }
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    if(all || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
-        for(int j = 0; j < (*it)->tetrahedra[i]->getNumVertices(); j++)
-          (*it)->tetrahedra[i]->getVertex(j)->setIndex(0);
-      for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
-        for(int j = 0; j < (*it)->hexahedra[i]->getNumVertices(); j++)
-          (*it)->hexahedra[i]->getVertex(j)->setIndex(0);
-      for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
-        for(int j = 0; j < (*it)->prisms[i]->getNumVertices(); j++)
-          (*it)->prisms[i]->getVertex(j)->setIndex(0);
-      for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
-        for(int j = 0; j < (*it)->pyramids[i]->getNumVertices(); j++)
-          (*it)->pyramids[i]->getVertex(j)->setIndex(0);
-    }
+  for(unsigned int i = 0; i < entities.size(); i++)
+    if(all || entities[i]->physicals.size())
+      for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+	for(int k = 0; k < entities[i]->getMeshElement(j)->getNumVertices(); k++)
+	  entities[i]->getMeshElement(j)->getVertex(k)->setIndex(0);
 
   // renumber all the mesh vertices tagged with 0
   int numVertices = 0;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      if(!(*it)->mesh_vertices[i]->getIndex())
-        (*it)->mesh_vertices[i]->setIndex(++numVertices);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      if(!(*it)->mesh_vertices[i]->getIndex())
-        (*it)->mesh_vertices[i]->setIndex(++numVertices);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      if(!(*it)->mesh_vertices[i]->getIndex())
-        (*it)->mesh_vertices[i]->setIndex(++numVertices);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      if(!(*it)->mesh_vertices[i]->getIndex())
-        (*it)->mesh_vertices[i]->setIndex(++numVertices);
-
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+      if(!entities[i]->mesh_vertices[j]->getIndex())
+        entities[i]->mesh_vertices[j]->setIndex(++numVertices);
+  
   return numVertices;
 }
 
-/*
-int GModel::unIndexMeshVerticesNotInPartition(int partition)
-{
-  // tag all the vertices not belonging to the partition with index -1
-  // (so they won't be saved)
-  Msg::Error("unIndexMeshVerticesNotInPartition not implemented");
-}
-*/
-
 void GModel::scaleMesh(double factor)
 {
-  std::vector<GEntity*> entities = getEntities();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
   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];
@@ -722,79 +591,21 @@ void GModel::scaleMesh(double factor)
 void GModel::recomputeMeshPartitions()
 {
   meshPartitions.clear();
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    (*it)->recomputeMeshPartitions();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    (*it)->recomputeMeshPartitions();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    (*it)->recomputeMeshPartitions();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    entities[i]->recomputeMeshPartitions();
 }
 
 void GModel::deleteMeshPartitions()
 {
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    (*it)->deleteMeshPartitions();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    (*it)->deleteMeshPartitions();
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    (*it)->deleteMeshPartitions();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    entities[i]->deleteMeshPartitions();
   meshPartitions.clear();
 }
 
-static int checkVertices(std::vector<MVertex*> &vertices,
-                         std::set<MVertex*, MVertexLessThanLexicographic> &pos)
-{
-  int num = 0;
-  for(unsigned int i = 0; i < vertices.size(); i++){
-    MVertex *v = vertices[i];
-    std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(v);
-    if(it == pos.end()){
-      pos.insert(v);
-    }
-    else{
-      Msg::Info("Vertices %d and %d have identical position (%g, %g, %g)",
-          (*it)->getNum(), v->getNum(), v->x(), v->y(), v->z());
-      num++;
-    }
-  }
-  return num;
-}
-
-template <class T>
-static int checkElements(int tag,
-                         std::vector<T*> &elements,
-                         std::set<MElement*, MElementLessThanLexicographic> &pos)
-{
-  int num = 0;
-  for(unsigned int i = 0; i < elements.size(); i++){
-    MElement *e = elements[i];
-    std::set<MElement*, MElementLessThanLexicographic>::iterator it = pos.find(e);
-    if(it == pos.end()){
-      pos.insert(e);
-    }
-    else{
-      char temp[256], temp2[256];
-      sprintf(temp, "Elements %d ( ", (*it)->getNum());
-      for (int i = 0; i < (*it)->getNumVertices();i++){
-        sprintf(temp2, "%d ", (*it)->getVertex(i)->getNum());
-        strcat(temp, temp2);
-      }
-      sprintf(temp2, ") on entity %d has same barycenter as element %d ( ",
-              tag, e->getNum());
-      strcat(temp, temp2);
-      for (int i = 0; i < e->getNumVertices(); i++){
-        sprintf(temp2, "%d ", e->getVertex(i)->getNum());
-        strcat(temp, temp2);
-      }
-      sprintf(temp2, ")");
-      strcat(temp, temp2);
-      Msg::Info("%s", temp);
-      num++;
-    }
-  }
-  return num;
-}
-
 void GModel::checkMeshCoherence()
 {
   int numEle = getNumMeshElements();
@@ -806,20 +617,29 @@ void GModel::checkMeshCoherence()
   double lc = bb.empty() ? 1. : norm(SVector3(bb.max(), bb.min()));
   double tol = CTX.geom.tolerance * lc;
 
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
   // check for duplicate mesh vertices
   {
     double old_tol = MVertexLessThanLexicographic::tolerance;
     MVertexLessThanLexicographic::tolerance = tol;
     std::set<MVertex*, MVertexLessThanLexicographic> pos;
     int num = 0;
-    for(viter it = firstVertex(); it != lastVertex(); ++it)
-      num += checkVertices((*it)->mesh_vertices, pos);
-    for(eiter it = firstEdge(); it != lastEdge(); ++it)
-      num += checkVertices((*it)->mesh_vertices, pos);
-    for(fiter it = firstFace(); it != lastFace(); ++it)
-      num += checkVertices((*it)->mesh_vertices, pos);
-    for(riter it = firstRegion(); it != lastRegion(); ++it)
-      num += checkVertices((*it)->mesh_vertices, pos);
+    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];
+	std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(v);
+	if(it == pos.end()){
+	  pos.insert(v);
+	}
+	else{
+	  Msg::Info("Vertices %d and %d have identical position (%g, %g, %g)",
+		    (*it)->getNum(), v->getNum(), v->x(), v->y(), v->z());
+	  num++;
+	}
+      }
+    }
     if(num) Msg::Warning("%d duplicate vertices", num);
     MVertexLessThanLexicographic::tolerance = old_tol;
   }
@@ -830,17 +650,33 @@ void GModel::checkMeshCoherence()
     MElementLessThanLexicographic::tolerance = tol;
     std::set<MElement*, MElementLessThanLexicographic> pos;
     int num = 0;
-    for(eiter it = firstEdge(); it != lastEdge(); ++it)
-      num += checkElements((*it)->tag(), (*it)->lines, pos);
-    for(fiter it = firstFace(); it != lastFace(); ++it){
-      num += checkElements((*it)->tag(), (*it)->triangles, pos);
-      num += checkElements((*it)->tag(), (*it)->quadrangles, pos);
-    }
-    for(riter it = firstRegion(); it != lastRegion(); ++it){
-      num += checkElements((*it)->tag(), (*it)->tetrahedra, pos);
-      num += checkElements((*it)->tag(), (*it)->hexahedra, pos);
-      num += checkElements((*it)->tag(), (*it)->prisms, pos);
-      num += checkElements((*it)->tag(), (*it)->pyramids, pos);
+    for(unsigned int i = 0; i < entities.size(); i++){
+      for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+	MElement *e = entities[i]->getMeshElement(j);
+	std::set<MElement*, MElementLessThanLexicographic>::iterator it = pos.find(e);
+	if(it == pos.end()){
+	  pos.insert(e);
+	}
+	else{
+	  char temp[256], temp2[256];
+	  sprintf(temp, "Element %d ( ", (*it)->getNum());
+	  for (int k = 0; k < (*it)->getNumVertices(); k++){
+	    sprintf(temp2, "%d ", (*it)->getVertex(k)->getNum());
+	    strcat(temp, temp2);
+	  }
+	  sprintf(temp2, ") on entity %d has same barycenter as element %d ( ",
+		  entities[i]->tag(), e->getNum());
+	  strcat(temp, temp2);
+	  for (int k = 0; k < e->getNumVertices(); k++){
+	    sprintf(temp2, "%d ", e->getVertex(k)->getNum());
+	    strcat(temp, temp2);
+	  }
+	  sprintf(temp2, ")");
+	  strcat(temp, temp2);
+	  Msg::Info("%s", temp);
+	  num++;
+	}
+      }
     }
     if(num) Msg::Warning("%d duplicate elements", num);
     MElementLessThanLexicographic::tolerance = old_tol;
@@ -860,6 +696,17 @@ void GModel::_storeElementsInEntities(std::map<int, std::vector<MElement*> > &ma
     if(!it->second.size()) continue;
     int numEdges = it->second[0]->getNumEdges();
     switch(numEdges){
+    case 0: 
+      {
+        GVertex *v = getVertexByTag(it->first);
+        if(!v){
+          v = new discreteVertex(this, it->first);
+          add(v);
+        }
+	if(v->points.empty()) // CAD points already have one by default
+	  _addElements(v->points, it->second);
+      }
+      break;
     case 1: 
       {
         GEdge *e = getEdgeByTag(it->first);
@@ -926,7 +773,37 @@ void GModel::_associateEntityWithMeshVertices()
     _associateEntityWithElementVertices(*it, (*it)->lines);
   }
   for(viter it = firstVertex(); it != lastVertex(); ++it){
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->setEntity(*it);
+    _associateEntityWithElementVertices(*it, (*it)->points);
+  }
+}
+
+void GModel::_storeVerticesInEntities(std::map<int, MVertex*> &vertices)
+{
+  std::map<int, MVertex*>::const_iterator it = vertices.begin();
+  for(; it != vertices.end(); ++it){
+    MVertex *v = it->second;
+    GEntity *ge = v->onWhat();
+    if(ge){
+      if(ge->dim() || ge->mesh_vertices.empty()) // special case for points
+	ge->mesh_vertices.push_back(v);
+    }
+    else
+      delete v; // we delete all unused vertices
+  }
+}
+
+void GModel::_storeVerticesInEntities(std::vector<MVertex*> &vertices)
+{
+  for(unsigned int i = 0; i < vertices.size(); i++){
+    MVertex *v = vertices[i];
+    if(v){ // the vector is allowed to have null entries
+      GEntity *ge = v->onWhat();
+      if(ge) {
+	if(ge->dim() || ge->mesh_vertices.empty()) // special case for points
+	  ge->mesh_vertices.push_back(v);
+      }
+      else
+        delete v; // we delete all unused vertices
+    }
   }
 }
diff --git a/Geo/GModel.h b/Geo/GModel.h
index a98a560b8a030e57e698a18c4844da6cea1713a1..e5ce013e2fa76d6ef980f3b0178ccb60523aed5b 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -62,6 +62,11 @@ class GModel
   // entity
   void _associateEntityWithMeshVertices();
 
+  // store the vertices in the entity they are associated with, and
+  // delete those that are not associated with any geo entity
+  void _storeVerticesInEntities(std::map<int, MVertex*> &vertices);
+  void _storeVerticesInEntities(std::vector<MVertex*> &vertices);
+
   // entity that is currently being meshed (used for error reporting)
   GEntity *_currentMeshEntity;
 
@@ -153,8 +158,8 @@ class GModel
   // Snap vertices on model edges by using geometry tolerance
   void snapVertices();
 
-  // Get a vector containing all the entities in the model
-  std::vector<GEntity*> getEntities();
+  // Fill a vector containing all the entities in the model
+  void getEntities(std::vector<GEntity*> &entities);
 
   // Checks if there are no physical entities in the model
   bool noPhysicalGroups();
diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index e93e562a9d9fb5160951aa8326f6e024b7562d64..213b04295cddca25ca3450acd869722861accdab 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -234,6 +234,7 @@ int GModel::readMED(const std::string &name, int meshIndex)
 			   (meshDim > 1) ? coord[meshDim * i + 1] : 0., 
 			   (meshDim > 2) ? coord[meshDim * i + 2] : 0.,
 			   0, nodeTags.empty() ? 0 : nodeTags[i]);
+
   // read elements
   for(int mshType = 0; mshType < 50; mshType++){ // loop over all possible MSH types
     med_geometrie_element type = msh2medElementType(mshType);
@@ -255,37 +256,19 @@ int GModel::readMED(const std::string &name, int meshIndex)
     std::vector<med_int> eleTags(numEle);
     if(MEDnumLire(fid, meshName, &eleTags[0], numEle, MED_MAILLE, type) < 0)
       eleTags.clear();
-    if(numNodPerEle == 1){ // special case for points
-      for(int j = 0; j < numEle; j++){    
-	GVertex *v = getVertexByTag(-fam[j]);
-	if(!v){
-	  v = new discreteVertex(this, -fam[j]);
-	  add(v);
-	}
-	v->mesh_vertices.push_back(verts[conn[j] - 1]);
-      }
-    }
-    else{
-      std::map<int, std::vector<MElement*> > elements;
-      MElementFactory factory;
-      for(int j = 0; j < numEle; j++){    
-	std::vector<MVertex*> v(numNodPerEle);
-	for(int k = 0; k < numNodPerEle; k++)
-	  v[k] = verts[conn[numNodPerEle * j + med2mshNodeIndex(type, k)] - 1];
-	MElement *e = factory.create(mshType, v, eleTags.empty() ? 0 : eleTags[j]);
-	if(e) elements[-fam[j]].push_back(e);
-      }
-      _storeElementsInEntities(elements);
+    std::map<int, std::vector<MElement*> > elements;
+    MElementFactory factory;
+    for(int j = 0; j < numEle; j++){    
+      std::vector<MVertex*> v(numNodPerEle);
+      for(int k = 0; k < numNodPerEle; k++)
+	v[k] = verts[conn[numNodPerEle * j + med2mshNodeIndex(type, k)] - 1];
+      MElement *e = factory.create(mshType, v, eleTags.empty() ? 0 : eleTags[j]);
+      if(e) elements[-fam[j]].push_back(e);
     }
+    _storeElementsInEntities(elements);
   }
   _associateEntityWithMeshVertices();
-  for(unsigned int i = 0; i < verts.size(); i++){
-    GEntity *ge = verts[i]->onWhat();
-    // store vertices (except for points, which are already ok)
-    if(ge && ge->dim() > 0) ge->mesh_vertices.push_back(verts[i]);
-    // delete unused vertices
-    if(!ge) delete verts[i];
-  }
+  _storeVerticesInEntities(verts);
 
   // read family info
   med_int numFamilies = MEDnFam(fid, meshName);
@@ -397,7 +380,8 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
   // get a vector containing all the geometrical entities in the
   // model (the ordering of the entities must be the same as the one
   // used during the indexing of the vertices)
-  std::vector<GEntity*> entities = getEntities();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
 
   std::map<GEntity*, int> families;
 
@@ -470,13 +454,10 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
     med_geometrie_element typ;
     { // points
       std::vector<med_int> conn, fam;
-      for(viter it = firstVertex(); it != lastVertex(); it++){
-	if(saveAll || (*it)->physicals.size()){
-	  conn.push_back((*it)->mesh_vertices[0]->getIndex());
-	  fam.push_back(families[*it]);
-	}
-      }
-      writeElementsMED(fid, meshName, conn, fam, MED_POINT1);
+      for(viter it = firstVertex(); it != lastVertex(); it++)
+	if(saveAll || (*it)->physicals.size())
+	  fillElementsMED(families[*it], (*it)->points, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // lines
       std::vector<med_int> conn, fam;
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 7499a99a2729a850e16e0c59375ed5e76f254e60..8b88b8622d14c3661ece5584a05672f1ef79e915 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -23,33 +23,6 @@
 #  include "Message.h"
 #endif
 
-static void storeVerticesInEntities(std::map<int, MVertex*> &vertices)
-{
-  std::map<int, MVertex*>::const_iterator it = vertices.begin();
-  for(; it != vertices.end(); ++it){
-    MVertex *v = it->second;
-    GEntity *ge = v->onWhat();
-    if(ge) 
-      ge->mesh_vertices.push_back(v);
-    else
-      delete v; // we delete all unused vertices
-  }
-}
-
-static void storeVerticesInEntities(std::vector<MVertex*> &vertices)
-{
-  for(unsigned int i = 0; i < vertices.size(); i++){
-    MVertex *v = vertices[i];
-    if(v){ // the vector can have null entries (first or last element)
-      GEntity *ge = v->onWhat();
-      if(ge) 
-        ge->mesh_vertices.push_back(v);
-      else
-        delete v; // we delete all unused vertices
-    }
-  }
-}
-
 static void storePhysicalTagsInEntities(GModel *m, int dim,
                                         std::map<int, std::map<int, std::string> > &map)
 {
@@ -103,37 +76,30 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
 
 static void createElementMSH(GModel *m, int num, int type, int physical, 
                              int reg, int part, std::vector<MVertex*> &v, 
-                             std::map<int, std::vector<MVertex*> > &points,
-                             std::map<int, std::vector<MElement*> > elem[7],
+                             std::map<int, std::vector<MElement*> > elements[7],
                              std::map<int, std::map<int, std::string> > physicals[4])
 {
-  int dim;
-  if(type == MSH_PNT){
-    dim = 0;
-    points[reg].push_back(v[0]);
-  }
-  else{
-    MElementFactory factory;
-    MElement *e = factory.create(type, v, num, part);
-
-    if(!e){
-      Msg::Error("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::Error("Wrong number of edges in element"); return;
-    }
-    elem[idx][reg].push_back(e);
-  }
+  MElementFactory factory;
+  MElement *e = factory.create(type, v, num, part);
+  if(!e){
+    Msg::Error("Unknown type of element %d", type);
+    return;
+  }
+
+  int dim = e->getDim();
+  int idx;
+  switch(e->getNumEdges()){
+  case 0 : idx = 0; break;
+  case 1 : idx = 1; break;
+  case 3 : idx = 2; break;
+  case 4 : idx = 3; break;
+  case 6 : idx = 4; break;
+  case 12 : idx = 5; break;
+  case 9 : idx = 6; break;
+  case 8 : idx = 7; break;
+  default : Msg::Error("Wrong number of edges in element"); return;
+  }
+  elements[idx][reg].push_back(e);
   
   if(physical && (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
     physicals[dim][reg][physical] = "unnamed";
@@ -152,8 +118,7 @@ int GModel::readMSH(const std::string &name)
   double version = 1.0;
   bool binary = false, swap = false;
   char str[256] = "XXX";
-  std::map<int, std::vector<MVertex*> > points;
-  std::map<int, std::vector<MElement*> > elements[7];
+  std::map<int, std::vector<MElement*> > elements[8];
   std::map<int, std::map<int, std::string> > physicals[4];
   bool postpro = false;
 
@@ -289,7 +254,7 @@ int GModel::readMSH(const std::string &name)
             if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
           }
           createElementMSH(this, num, type, physical, elementary, partition, 
-                           vertices, points, elements, physicals);
+                           vertices, elements, physicals);
 	  if(numElements > 100000) 
 	    Msg::ProgressMeter(i + 1, numElements, "Reading elements");
         }
@@ -322,7 +287,7 @@ int GModel::readMSH(const std::string &name)
               if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
             }
             createElementMSH(this, num, type, physical, elementary, partition, 
-                             vertices, points, elements, physicals);
+                             vertices, elements, physicals);
 	    if(numElements > 100000) 
 	      Msg::ProgressMeter(numElementsPartial + i + 1, numElements, 
 				 "Reading elements");
@@ -361,36 +326,18 @@ int GModel::readMSH(const std::string &name)
   }
 
   // store the elements in their associated elementary entity. If the
-  // entity does not exist, create a new one.
+  // entity does not exist, create a new (discrete) one.
   for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
     _storeElementsInEntities(elements[i]);
 
-  // treat points separately
-  for(std::map<int, std::vector<MVertex*> >::iterator it = points.begin(); 
-      it != points.end(); ++it){
-    GVertex *v = getVertexByTag(it->first);
-    if(!v){
-      v = new discreteVertex(this, it->first);
-      add(v);
-    }
-    for(unsigned int i = 0; i < it->second.size(); i++) 
-      v->mesh_vertices.push_back(it->second[i]);
-  }
-
   // associate the correct geometrical entity with each mesh vertex
   _associateEntityWithMeshVertices();
 
-  // special case for geometry vertices: now that the correct
-  // geometrical entity has been associated with the vertices, we
-  // reset mesh_vertices so that it can be filled again below
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    (*it)->mesh_vertices.clear();
-
   // store the vertices in their associated geometrical entity
   if(vertexVector.size())
-    storeVerticesInEntities(vertexVector);
+    _storeVerticesInEntities(vertexVector);
   else
-    storeVerticesInEntities(vertexMap);
+    _storeVerticesInEntities(vertexMap);
 
   // store the physical tags
   for(int i = 0; i < 4; i++)
@@ -480,8 +427,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   std::map<int,int> elements;
   for(viter it = firstVertex(); it != lastVertex(); ++it){
     int p = (saveAll ? 1 : (*it)->physicals.size());
-    int n = p * (*it)->mesh_vertices.size();
-    if(n) elements[MSH_PNT] += n;
+    int n = p * (*it)->points.size();
+    if(n) elements[(*it)->points[0]->getTypeForMSH()] += n;
   }
   for(eiter it = firstEdge(); it != lastEdge(); ++it){
     int p = (saveAll ? 1 : (*it)->physicals.size());
@@ -536,18 +483,12 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
     fprintf(fp, "$NOD\n");
 
   fprintf(fp, "%d\n", numVertices);
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
+
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
+      entities[i]->mesh_vertices[j]->writeMSH(fp, binary, scalingFactor);
 
   if(binary) fprintf(fp, "\n");
 
@@ -565,7 +506,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
 
   writeElementHeaderMSH(binary, fp, elements, MSH_PNT);
   for(viter it = firstVertex(); it != lastVertex(); ++it)
-    writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, binary, num,
+    writeElementsMSH(fp, (*it)->points, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
   writeElementHeaderMSH(binary, fp, elements, MSH_LIN_2, MSH_LIN_3, MSH_LIN_4, MSH_LIN_5);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
@@ -657,42 +598,14 @@ int GModel::writePOS(const std::string &name, bool printElementary,
   fprintf(fp, "View \"Statistics\" {\n");
   fprintf(fp, "T2(1.e5,30,%d){%s};\n", (1<<16)|(4<<8), names.c_str());
 
-  for(eiter it = firstEdge(); it != lastEdge(); ++it) {
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->lines.size(); i++)
-        (*it)->lines[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                  scalingFactor, (*it)->tag());
-    }
-  }
-
-  for(fiter it = firstFace(); it != lastFace(); ++it) {
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-        (*it)->triangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                      scalingFactor, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-        (*it)->quadrangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                        scalingFactor, (*it)->tag());
-    }
-  }
-
-  for(riter it = firstRegion(); it != lastRegion(); ++it) {
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
-        (*it)->tetrahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
-                                       scalingFactor, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
-        (*it)->hexahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                      scalingFactor, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
-        (*it)->prisms[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                   scalingFactor, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
-        (*it)->pyramids[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5], 
-                                     scalingFactor, (*it)->tag());
-    }
-  }
-
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    if(saveAll || entities[i]->physicals.size())
+      for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+	entities[i]->getMeshElement(j)->writePOS(fp, f[0], f[1], f[2], f[3],
+						 f[4], f[5], scalingFactor, 
+						 entities[i]->tag());
   fprintf(fp, "};\n");
 
   fclose(fp);
@@ -1021,7 +934,7 @@ int GModel::readVRML(const std::string &name)
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     _storeElementsInEntities(elements[i]);
   _associateEntityWithMeshVertices();
-  storeVerticesInEntities(allVertexVector);
+  _storeVerticesInEntities(allVertexVector);
 
   fclose(fp);
   return 1;
@@ -1224,7 +1137,8 @@ int GModel::readUNV(const std::string &name)
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     _storeElementsInEntities(elements[i]);
   _associateEntityWithMeshVertices();
-  storeVerticesInEntities(vertexMap);
+  _storeVerticesInEntities(vertexMap);
+
   for(int i = 0; i < 4; i++)  
     storePhysicalTagsInEntities(this, i, physicals[i]);
 
@@ -1232,18 +1146,6 @@ int GModel::readUNV(const std::string &name)
   return 1;
 }
 
-template<class T>
-static void writeElementsUNV(FILE *fp, const std::vector<T*> &ele, int saveAll, 
-                             int &num, int elementary, std::vector<int> &physicals)
-{
-  for(unsigned int i = 0; i < ele.size(); i++)
-    if(saveAll)
-      ele[i]->writeUNV(fp, ++num, elementary, 0);
-    else
-      for(unsigned int j = 0; j < physicals.size(); j++)
-        ele[i]->writeUNV(fp, ++num, elementary, physicals[j]);
-}
-
 int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNodes, 
                      double scalingFactor)
 {
@@ -1257,38 +1159,30 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
 
   indexMeshVertices(saveAll);
 
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
   // nodes
   fprintf(fp, "%6d\n", -1);
   fprintf(fp, "%6d\n", 2411);
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
+      entities[i]->mesh_vertices[j]->writeUNV(fp, scalingFactor);
   fprintf(fp, "%6d\n", -1);  
 
   // elements
   fprintf(fp, "%6d\n", -1);
   fprintf(fp, "%6d\n", 2412);
   int num = 0;
-  for(eiter it = firstEdge(); it != lastEdge(); ++it){
-    writeElementsUNV(fp, (*it)->lines, saveAll, num, (*it)->tag(), (*it)->physicals);
-  }
-  for(fiter it = firstFace(); it != lastFace(); ++it){
-    writeElementsUNV(fp, (*it)->triangles, saveAll, num, (*it)->tag(), (*it)->physicals);
-    writeElementsUNV(fp, (*it)->quadrangles, saveAll, num, (*it)->tag(), (*it)->physicals);
-  }
-  for(riter it = firstRegion(); it != lastRegion(); ++it){
-    writeElementsUNV(fp, (*it)->tetrahedra, saveAll, num, (*it)->tag(), (*it)->physicals);
-    writeElementsUNV(fp, (*it)->hexahedra, saveAll, num, (*it)->tag(), (*it)->physicals);
-    writeElementsUNV(fp, (*it)->prisms, saveAll, num, (*it)->tag(), (*it)->physicals);
+  for(unsigned int i = 0; i < entities.size(); i++){
+    for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+      MElement *e = entities[i]->getMeshElement(j);
+      if(saveAll)
+	e->writeUNV(fp, ++num, entities[i]->tag(), 0);
+      else
+	for(unsigned int k = 0; k < entities[i]->physicals.size(); k++)
+	  e->writeUNV(fp, ++num, entities[i]->tag(), entities[i]->physicals[k]);
+    }
   }
   fprintf(fp, "%6d\n", -1);
 
@@ -1474,7 +1368,7 @@ int GModel::readMESH(const std::string &name)
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     _storeElementsInEntities(elements[i]);
   _associateEntityWithMeshVertices();
-  storeVerticesInEntities(vertexVector);
+  _storeVerticesInEntities(vertexVector);
 
   fclose(fp);
   return 1;
@@ -1498,18 +1392,11 @@ int GModel::writeMESH(const std::string &name, bool saveAll, double scalingFacto
 
   fprintf(fp, " Vertices\n");
   fprintf(fp, " %d\n", numVertices);
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
+      entities[i]->mesh_vertices[j]->writeMESH(fp, scalingFactor);
   
   int numTriangles = 0, numQuadrangles = 0, numTetrahedra = 0;
   for(fiter it = firstFace(); it != lastFace(); ++it){
@@ -1831,7 +1718,7 @@ int GModel::readBDF(const std::string &name)
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     _storeElementsInEntities(elements[i]);
   _associateEntityWithMeshVertices();
-  storeVerticesInEntities(vertexMap);
+  _storeVerticesInEntities(vertexMap);
 
   fclose(fp);
   return 1;
@@ -1852,46 +1739,20 @@ int GModel::writeBDF(const std::string &name, int format, bool saveAll,
 
   fprintf(fp, "$ Created by Gmsh\n");
 
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
-  
-  for(eiter it = firstEdge(); it != lastEdge(); ++it){
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->lines.size(); i++)
-        (*it)->lines[i]->writeBDF(fp, format, (*it)->tag());
-    }
-  }
-  for(fiter it = firstFace(); it != lastFace(); ++it){
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-        (*it)->triangles[i]->writeBDF(fp, format, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-        (*it)->quadrangles[i]->writeBDF(fp, format, (*it)->tag());
-    }
-  }
-  for(riter it = firstRegion(); it != lastRegion(); ++it){
-    if(saveAll || (*it)->physicals.size()){
-      for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
-        (*it)->tetrahedra[i]->writeBDF(fp, format, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
-        (*it)->hexahedra[i]->writeBDF(fp, format, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
-        (*it)->prisms[i]->writeBDF(fp, format, (*it)->tag());
-      for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
-        (*it)->pyramids[i]->writeBDF(fp, format, (*it)->tag());
-    }
-  }
-  
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
+  // nodes
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
+      entities[i]->mesh_vertices[j]->writeBDF(fp, scalingFactor);
+
+  // elements
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+      if(saveAll || entities[i]->physicals.size())
+	entities[i]->getMeshElement(j)->writeBDF(fp, format, entities[i]->tag());
+
   fprintf(fp, "ENDDATA\n");
    
   fclose(fp);
@@ -2103,7 +1964,8 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
   fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
 
   // get all the entities in the model
-  std::vector<GEntity*> entities = getEntities();
+  std::vector<GEntity*> entities;
+  getEntities(entities);
 
   // write mesh vertices
   fprintf(fp, "POINTS %d double\n", numVertices);
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index b4f55b540dbadd88d0e7ac4cba6dbf4be7b4ce9e..cfe39c7161b70ea2d96de56c9f0124a05ca6443b 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -116,38 +116,6 @@ void GRegion::setVisibility(char val, bool recursive)
   }
 }
 
-void GRegion::recomputeMeshPartitions()
-{
-  for(unsigned int i = 0; i < tetrahedra.size(); i++) {
-    int part = tetrahedra[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-  for(unsigned int i = 0; i < hexahedra.size(); i++) {
-    int part = hexahedra[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-  for(unsigned int i = 0; i < prisms.size(); i++) {
-    int part = prisms[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-  for(unsigned int i = 0; i < pyramids.size(); i++) {
-    int part = pyramids[i]->getPartition();
-    if(part) model()->getMeshPartitions().insert(part);
-  }
-}
-
-void GRegion::deleteMeshPartitions()
-{
-  for(unsigned int i = 0; i < tetrahedra.size(); i++)
-    tetrahedra[i]->setPartition(0);
-  for(unsigned int i = 0; i < hexahedra.size(); i++)
-    hexahedra[i]->setPartition(0);
-  for(unsigned int i = 0; i < prisms.size(); i++)
-    prisms[i]->setPartition(0);
-  for(unsigned int i = 0; i < pyramids.size(); i++)
-    pyramids[i]->setPartition(0);
-}
-
 std::string GRegion::getAdditionalInfoString()
 {
   if(l_faces.empty()) return std::string("");
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index ccb61be903e37100eb5a658e815f024cec268dbb..3fb32d99953c936f56d823ada8eaf6a5837d0421 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -44,12 +44,6 @@ class GRegion : public GEntity {
   // check if the region is connected to another region by an edge
   bool edgeConnected(GRegion *r) const;
 
-  // recompute the mesh partitions defined on this region
-  void recomputeMeshPartitions();
-
-  // delete the mesh partitions defined on this region
-  void deleteMeshPartitions();
-
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index 64a797cf8a5c5c421b0955f9abd6916d64a64a36..f1b52a26effcc404e24454df0db55da29bdeab76 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -7,7 +7,7 @@
 #include <algorithm>
 #include "GVertex.h"
 #include "GFace.h"
-#include "MVertex.h"
+#include "MElement.h"
 
 #if defined(HAVE_GMSH_EMBEDDED)
 #  include "GmshEmbedded.h"
@@ -23,6 +23,9 @@ GVertex::~GVertex()
 {
   for(unsigned int i = 0; i < mesh_vertices.size(); i++)
     delete mesh_vertices[i];
+
+  for(unsigned int i = 0; i < points.size(); i++)
+    delete points[i];
 }
 
 void GVertex::setPosition(GPoint &p)
@@ -57,3 +60,15 @@ std::string GVertex::getAdditionalInfoString()
   }
   return std::string(str);
 }
+
+unsigned int GVertex::getNumMeshElements()
+{
+  return points.size(); 
+}
+
+MElement *GVertex::getMeshElement(unsigned int index)
+{ 
+  if(index < points.size())
+    return points[index]; 
+  return 0;
+}
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 9d477a00e7494c3b5eeb693471013058d5a52625..61b6ee2adcff2148a75ffbc56784b7e985a412e3 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -10,6 +10,9 @@
 #include "GPoint.h"
 #include "SPoint2.h"
 
+class MElement;
+class MPoint;
+
 // A model vertex.
 class GVertex : public GEntity 
 {
@@ -19,22 +22,47 @@ class GVertex : public GEntity
  public:
   GVertex(GModel *m, int tag, double ms=1.e22);
   virtual ~GVertex();
+
+  // get/set the coordinates of the vertex
   virtual GPoint point() const = 0;
   virtual double x() const = 0;
   virtual double y() const = 0;
   virtual double z() const = 0;
   virtual void setPosition(GPoint &p);
+
+  // add/delete an edge bounded by this vertex
   void addEdge(GEdge *e);
   void delEdge(GEdge *e);
+
+  // get the dimension of the vertex (0)
   virtual int dim() const { return 0; }
+
+  // get the geometric type of the vertex
   virtual GeomType geomType() const { return Point; }
+
+  // get/set the prescribed mesh size at the vertex
   inline double prescribedMeshSizeAtVertex() const { return meshSize; }
   virtual void setPrescribedMeshSizeAtVertex(double l) { meshSize = l; }
+
+  // get the bounding box
   virtual SBoundingBox3d bounds() const { return SBoundingBox3d(SPoint3(x(), y(), z())); }
+
+  // reparmaterize the point onto the given face
   virtual SPoint2 reparamOnFace(GFace *gf, int) const;
+
+  // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
-  virtual std::list<GEdge*> edges() const{ return l_edges; }
+
+  // get the edges that this vertex bounds
+   virtual std::list<GEdge*> edges() const{ return l_edges; }
+
+  // get number of elements in the mesh
+  unsigned int getNumMeshElements();
+
+  // get the element at the given index
+  MElement *getMeshElement(unsigned int index);
+
+  std::vector<MPoint*> points;
 };
 
 #endif
-
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index d93d5297a871f548963de5fdc9fe5f0caab929b0..2e75b164e234c10368d4bab88d68a55103efa3fe 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -17,6 +17,11 @@ class MEdge {
   char _si[2]; // sorted indices
 
  public:
+  MEdge()
+  {
+    _v[0] = _v[1] = 0;
+    _si[0] = _si[1] = 0;
+  }
   MEdge(MVertex *v0, MVertex *v1)
   {
     _v[0] = v0; _v[1] = v1;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 20a7f3f10e039e754dd044808495a942c6e89cb2..6b95f2cb19af5c8f093cb275cd285595c1a863be 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -996,7 +996,7 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
 				  int num, int part)
 {
   switch (type) {
-  case MSH_PNT:    return 0;
+  case MSH_PNT:    return new MPoint(v, num, part);
   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);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index c976e651db2d9c28b2e6923c0343c65483171b7f..63cf9e5eda31a63b3ec686ff0b84c69259f1d782 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -114,7 +114,9 @@ class MElement
 
   // get all the vertices on an edge
   virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
-    = 0;
+  {
+    v.resize(0);
+  }
 
   // get the faces
   virtual int getNumFaces() = 0;
@@ -144,7 +146,7 @@ class MElement
   virtual SPoint3 barycenter();
 
   // revert the orientation of the element
-  virtual void revert() = 0;
+  virtual void revert(){}
 
   // compute and change the orientation of 3D elements to get
   // positive volume
@@ -236,6 +238,52 @@ class MElementFactory{
   MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0);
 };
 
+/*
+ * MPoint
+ *
+ */
+class MPoint : public MElement {
+ protected:
+  MVertex *_v[1];
+ public :
+  MPoint(MVertex *v0, int num=0, int part=0) 
+    : MElement(num, part)
+  {
+    _v[0] = v0;
+  }
+  MPoint(std::vector<MVertex*> &v, int num=0, int part=0) 
+    : MElement(num, part)
+  {
+    _v[0] = v[0];
+  }
+  ~MPoint(){}
+  virtual int getDim(){ return 0; }
+  virtual int getNumVertices(){ return 1; }
+  virtual MVertex *getVertex(int num){ return _v[0]; }
+  virtual int getNumEdges(){ return 0; }
+  virtual MEdge getEdge(int num){ return MEdge(); }
+  virtual int getNumEdgesRep(){ return 0; }
+  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n){}
+  virtual int getNumFaces(){ return 0; }
+  virtual MFace getFace(int num){ return MFace(); }
+  virtual int getNumFacesRep(){ return 0; }
+  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){}
+  virtual int getTypeForMSH(){ return MSH_PNT; }
+  virtual const char *getStringForPOS(){ return "SP"; }
+  virtual void getShapeFunction(int num, double u, double v, double w, double &s) 
+  {
+    s = 1.;
+  }
+  virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]) 
+  {
+    s[0] = s[1] = s[2] = 0.;
+  }
+  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  {
+    return true;
+  }
+};
+
 /*
  * MLine
  *
@@ -2876,6 +2924,7 @@ template <> struct DimTr<2>
     element->getEdgeVertices(iFace, v);
   }
 };
+
 template <> struct DimTr<3>
 {
   typedef MFace FaceT;
@@ -2893,4 +2942,5 @@ template <> struct DimTr<3>
     element->getFaceVertices(iFace, v);
   }
 };
+
 #endif
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index fbbe266c80141efabe309c84bb4c6cc2540d86ea..e03508e9e3177823aac70a56cddb5aef4e01c0ab 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -45,25 +45,6 @@ void MVertex::writeMSH(FILE *fp, bool binary, double scalingFactor)
   }
 }
 
-void MVertex::writeMSH(FILE *fp, double version, bool binary, int num, 
-                       int elementary, int physical)
-{
-  if(!binary){
-    fprintf(fp, "%d 15", num);
-    if(version < 2.0)
-      fprintf(fp, " %d %d 1", physical, elementary);
-    else
-      fprintf(fp, " 2 %d %d", physical, elementary);
-    fprintf(fp, " %d\n", _index);
-  }
-  else{
-    int tags[4] = {num, physical, elementary, 0};
-    fwrite(tags, sizeof(int), 4, fp);
-    int verts[1] = {_index};
-    fwrite(verts, sizeof(int), 1, fp);
-  }
-}
-
 void MVertex::writeVRML(FILE *fp, double scalingFactor)
 {
   if(_index < 0) return; // negative index vertices are never saved
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 3d663fc0274ff8b2c2389bb916b3ffda56ff173a..9f894222a709a0a3207a40c50c7ac31c3ff544d8 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -107,8 +107,6 @@ class MVertex{
 
   // IO routines
   void writeMSH(FILE *fp, bool binary=false, double scalingFactor=1.0);
-  void writeMSH(FILE *fp, double version, bool binary, int num, 
-                int elementary, int physical);
   void writeVRML(FILE *fp, double scalingFactor=1.0);
   void writeUNV(FILE *fp, double scalingFactor=1.0);
   void writeVTK(FILE *fp, bool binary=false, double scalingFactor=1.0);
diff --git a/Geo/Makefile b/Geo/Makefile
index 68c9c1103f2f0ebfba460f20102eb911f86d6cd0..be8292cd1b6070365be591e95ac5852d19518ef5 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -67,12 +67,15 @@ depend:
 	rm -f Makefile.new
 
 # DO NOT DELETE THIS LINE
-GEntity.o: GEntity.cpp GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
-  ../Common/VertexArray.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+GEntity.o: GEntity.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/VertexArray.h ../Geo/SVector3.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
 GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GFace.h GEdgeLoop.h GEdge.h \
-  SVector3.h Pair.h MVertex.h ../Common/Message.h
+  SVector3.h Pair.h MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h \
+  MFace.h ../Common/Message.h
 GEdge.o: GEdge.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 \
@@ -97,8 +100,9 @@ gmshVertex.o: gmshVertex.cpp GFace.h GEntity.h Range.h SPoint3.h \
   SVector3.h Pair.h gmshVertex.h Geo.h ../Common/GmshDefines.h \
   gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h MVertex.h \
-  GeoInterpolation.h ../Common/Message.h
+  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
+  GeoInterpolation.h ../Common/Message.h MVertex.h MElement.h MEdge.h \
+  MFace.h
 gmshEdge.o: gmshEdge.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 gmshEdge.h Geo.h ../Common/GmshDefines.h \
@@ -127,23 +131,23 @@ gmshSurface.o: gmshSurface.cpp gmshSurface.h Pair.h Range.h SPoint2.h \
   ../Numeric/NumericEmbedded.h ../Common/Message.h
 OCCVertex.o: OCCVertex.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 OCCVertex.h OCCIncludes.h MVertex.h \
-  OCCEdge.h OCCFace.h
+  GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h OCCEdge.h \
+  OCCFace.h MVertex.h MElement.h ../Common/GmshDefines.h MEdge.h MFace.h
 OCCEdge.o: OCCEdge.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 ../Common/Message.h OCCEdge.h OCCVertex.h \
-  OCCIncludes.h MVertex.h OCCFace.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
+  OCCIncludes.h OCCFace.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/PartitionOptions.h
 OCCFace.o: OCCFace.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 OCCVertex.h OCCIncludes.h MVertex.h \
-  OCCEdge.h OCCFace.h ../Common/Message.h ../Numeric/Numeric.h \
+  GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h OCCEdge.h \
+  OCCFace.h ../Common/Message.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h ../Common/VertexArray.h ../Geo/SVector3.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
 OCCRegion.o: OCCRegion.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 OCCVertex.h OCCIncludes.h MVertex.h \
-  OCCEdge.h OCCFace.h OCCRegion.h ../Common/Message.h
+  GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h OCCEdge.h \
+  OCCFace.h OCCRegion.h ../Common/Message.h
 discreteEdge.o: discreteEdge.cpp discreteEdge.h 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 Geo.h \
@@ -195,9 +199,9 @@ GModelIO_Geo.o: GModelIO_Geo.cpp GModel.h GVertex.h GEntity.h Range.h \
   gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
   ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  ../Common/OpenFile.h ../Common/Message.h gmshVertex.h MVertex.h \
-  gmshFace.h gmshEdge.h gmshRegion.h ../Mesh/Field.h ../Geo/Geo.h \
-  ../Post/PView.h ../Geo/SPoint3.h ../Parser/Parser.h
+  ../Common/OpenFile.h ../Common/Message.h gmshVertex.h gmshFace.h \
+  gmshEdge.h gmshRegion.h ../Mesh/Field.h ../Geo/Geo.h ../Post/PView.h \
+  ../Geo/SPoint3.h ../Parser/Parser.h
 GModelIO_Mesh.o: GModelIO_Mesh.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 ../Common/GmshDefines.h MElement.h \
@@ -208,8 +212,8 @@ GModelIO_OCC.o: GModelIO_OCC.cpp GModelIO_OCC.h 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 OCCIncludes.h \
   ../Common/Message.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/PartitionOptions.h OCCVertex.h MVertex.h OCCEdge.h OCCFace.h \
-  OCCRegion.h MElement.h ../Common/GmshDefines.h MEdge.h MFace.h \
+  ../Mesh/PartitionOptions.h OCCVertex.h OCCEdge.h OCCFace.h OCCRegion.h \
+  MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
   ../Common/OpenFile.h
 GModelIO_Fourier.o: GModelIO_Fourier.cpp GModel.h GVertex.h GEntity.h \
   Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
@@ -218,10 +222,7 @@ GModelIO_Fourier.o: GModelIO_Fourier.cpp GModel.h GVertex.h GEntity.h \
   GModelIO_Fourier.h
 GModelIO_CGNS.o: GModelIO_CGNS.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 ../Common/Message.h MZone.h \
-  MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
-  MEdgeHash.h ../Common/Hash.h MFaceHash.h CustomContainer.h \
-  MZoneBoundary.h CGNSOptions.h
+  GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h CGNSOptions.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 ../Common/Message.h
@@ -289,14 +290,8 @@ MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
   ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
   ../Numeric/GaussLegendre1D.h ../Common/Message.h ../Common/Context.h \
   ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h \
-  ../Mesh/qualityMeasures.h
-MZone.o: MZone.cpp MZone.h MElement.h ../Common/GmshDefines.h MVertex.h \
-  SPoint3.h MEdge.h SVector3.h MFace.h GFace.h GEntity.h Range.h \
-  SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h \
-  Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h MFaceHash.h \
-  CustomContainer.h
-MZoneBoundary.o: MZoneBoundary.cpp MZoneBoundary.h MZone.h MElement.h \
-  ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
-  GFace.h GEntity.h Range.h SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h \
-  GVertex.h SPoint2.h Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h \
-  MFaceHash.h CustomContainer.h ../Common/Message.h
+  ../Mesh/qualityMeasures.h ../Mesh/meshGFaceDelaunayInsertion.h \
+  ../Geo/MElement.h ../Mesh/meshGRegionDelaunayInsertion.h \
+  ../Mesh/BackgroundMesh.h ../Mesh/qualityMeasures.h
+MZone.o: MZone.cpp
+MZoneBoundary.o: MZoneBoundary.cpp
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 3be07f4d83981d3072bfe3e5f9761afc19272fa9..f5c142b5be0d6ae40660056251abb037790b15e7 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -7,9 +7,35 @@
 #include "OCCVertex.h"
 #include "OCCEdge.h"
 #include "OCCFace.h"
+#include "MVertex.h"
+#include "MElement.h"
 
 #if defined(HAVE_OCC)
 
+OCCVertex::OCCVertex(GModel *m, int num, TopoDS_Vertex _v) 
+  : GVertex(m, num), v(_v)
+{
+  max_curvature = -1;
+  gp_Pnt pnt = BRep_Tool::Pnt(v);
+  _x = pnt.X();
+  _y = pnt.Y();
+  _z = pnt.Z();
+  mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
+  points.push_back(new MPoint(mesh_vertices.back()));
+}
+
+void OCCVertex::setPosition(GPoint &p)
+{
+  _x = p.x();
+  _y = p.y();
+  _z = p.z();
+  if(mesh_vertices.size()){
+    mesh_vertices[0]->x() = p.x();
+    mesh_vertices[0]->y() = p.y();
+    mesh_vertices[0]->z() = p.z();
+  }
+}
+
 double max_surf_curvature(const GVertex *gv, double x, double y, double z, 
                           const GEdge *_myGEdge)
 {
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index fe32155b3ad9d3d48856e470817c73686d8902be..4c359b4eba49bc968d357814146035dd41b70cbe 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -9,7 +9,6 @@
 #include "GModel.h"
 #include "OCCIncludes.h"
 #include "GVertex.h"
-#include "MVertex.h"
 
 #if defined(HAVE_OCC)
 
@@ -20,32 +19,16 @@ class OCCVertex : public GVertex {
   mutable double max_curvature;
   double max_curvature_of_surfaces() const;
  public:
-  OCCVertex(GModel *m, int num, TopoDS_Vertex _v) : GVertex(m, num), v(_v)
-  {
-    max_curvature = -1;
-    gp_Pnt pnt = BRep_Tool::Pnt(v);
-    _x = pnt.X();
-    _y = pnt.Y();
-    _z = pnt.Z();
-    mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
-  }
+  OCCVertex(GModel *m, int num, TopoDS_Vertex _v);
   virtual ~OCCVertex() {}
   virtual GPoint point() const { return GPoint(x(), y(), z()); }
   virtual double x() const { return _x; }
   virtual double y() const { return _y; }
   virtual double z() const { return _z; }
-  virtual void setPosition(GPoint &p)
-  {
-    _x = p.x();
-    _y = p.y();
-    _z = p.z();
-    delete mesh_vertices[0];
-    mesh_vertices.clear();
-    mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
-  }
+  virtual void setPosition(GPoint &p);
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*)&v; }
-  virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
+  virtual SPoint2 reparamOnFace(GFace *gf, int) const;
 };
 
 #endif
diff --git a/Geo/discreteVertex.h b/Geo/discreteVertex.h
index 88861961e5317bc7262bcd6d26a5bdfacc278b5d..885f2db7694bbf7d2ae484ecd7de9711f33d9864 100644
--- a/Geo/discreteVertex.h
+++ b/Geo/discreteVertex.h
@@ -12,8 +12,8 @@
 
 class discreteVertex : public GVertex {
  public:
-  discreteVertex(GModel *m, int num) : GVertex(m, num) {}
-  virtual ~discreteVertex() {}
+  discreteVertex(GModel *m, int num) : GVertex(m, num){}
+  virtual ~discreteVertex(){}
   virtual GPoint point() const { return GPoint(x(), y(), z(), this); }
   virtual double x() const 
   { 
diff --git a/Geo/gmshVertex.cpp b/Geo/gmshVertex.cpp
index 2392b94013b0f802775ef283679cee8b381ee353..99f5d889ba756ba9801e2195652485df22fdaffa 100644
--- a/Geo/gmshVertex.cpp
+++ b/Geo/gmshVertex.cpp
@@ -9,6 +9,35 @@
 #include "Geo.h"
 #include "GeoInterpolation.h"
 #include "Message.h"
+#include "MVertex.h"
+#include "MElement.h"
+
+gmshVertex::gmshVertex(GModel *m, Vertex *_v)
+  : GVertex(m, _v->Num, _v->lc), v(_v)
+{
+  mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
+  points.push_back(new MPoint(mesh_vertices.back()));
+}
+
+void gmshVertex::setPosition(GPoint &p)
+{
+  v->Pos.X = p.x();
+  v->Pos.Y = p.y();
+  v->Pos.Z = p.z();
+  if(mesh_vertices.size()){
+    mesh_vertices[0]->x() = p.x();
+    mesh_vertices[0]->y() = p.y();
+    mesh_vertices[0]->z() = p.z();
+  }
+}
+
+GEntity::GeomType gmshVertex::geomType() const
+{
+  if(v->Typ == MSH_POINT_BND_LAYER)
+    return BoundaryLayerPoint;
+  else
+    return Point;
+}
 
 SPoint2 gmshVertex::reparamOnFace(GFace *face, int dir) const
 {
@@ -68,11 +97,3 @@ SPoint2 gmshVertex::reparamOnFace(GFace *face, int dir) const
     return GVertex::reparamOnFace(face, dir);
   }
 }
-
-GEntity::GeomType gmshVertex::geomType() const
-{
-  if(v->Typ == MSH_POINT_BND_LAYER)
-    return BoundaryLayerPoint;
-  else
-    return Point;
-}
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
index 319fa57aab7d8265610318c2a0b09fe0264696b6..f7082a1ef0010c24f46767ac70d2046c22121cbf 100644
--- a/Geo/gmshVertex.h
+++ b/Geo/gmshVertex.h
@@ -8,17 +8,13 @@
 
 #include "Geo.h"
 #include "GVertex.h"
-#include "MVertex.h"
 
 class gmshVertex : public GVertex {
  protected:
   Vertex *v;
 
  public:
-  gmshVertex(GModel *m, Vertex *_v) : GVertex(m, _v->Num, _v->lc), v(_v)
-  {
-    mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
-  }
+  gmshVertex(GModel *m, Vertex *_v);
   virtual ~gmshVertex() {}
   virtual GPoint point() const 
   {
@@ -27,17 +23,7 @@ class gmshVertex : public GVertex {
   virtual double x() const { return v->Pos.X; }
   virtual double y() const { return v->Pos.Y; }
   virtual double z() const { return v->Pos.Z; }
-  virtual void setPosition(GPoint &p)
-  {
-    v->Pos.X = p.x();
-    v->Pos.Y = p.y();
-    v->Pos.Z = p.z();
-    if(mesh_vertices.size()){
-      mesh_vertices[0]->x() = p.x();
-      mesh_vertices[0]->y() = p.y();
-      mesh_vertices[0]->z() = p.z();
-    }
-  }
+  virtual void setPosition(GPoint &p);
   virtual GeomType geomType() const;
   ModelType getNativeType() const { return GmshModel; }
   void *getNativePtr() const { return v; }
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index 5053bf84ddc7793771109159aff3860825ba9007..bccb8bd837a57722ce939020040ab21ef8944d17 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -331,7 +331,7 @@ static int BuildDelaunay(DocRecord *doc)
 
 // This routine insert the point 'newPoint' in the list dlist,
 // respecting the clock-wise orientation
-static int DListInsert(DListRecord **dlist, MPoint center, PointNumero newPoint)
+static int DListInsert(DListRecord **dlist, DPoint center, PointNumero newPoint)
 {
   DListRecord *p, *newp;
   double alpha1, alpha, beta, xx, yy;
diff --git a/Mesh/DivideAndConquer.h b/Mesh/DivideAndConquer.h
index 6f7b298427dc3be7ff5e1d0edd91b07de1710647..e75d61e3aa052325c3a25d00f8fab02ca8bfbff0 100644
--- a/Mesh/DivideAndConquer.h
+++ b/Mesh/DivideAndConquer.h
@@ -12,10 +12,10 @@ typedef int PointNumero;
 typedef struct{
   double v;
   double h;
-}MPoint;
+}DPoint;
 
 typedef struct{
-  MPoint where;
+  DPoint where;
   DListPeek adjacent;
   void *data;
 }PointRecord;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 035da470048fc683c22bf5a2839b193ddd54f806..0566d605466c6ca2ea7c7891f76e4cb32b73a961 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -391,16 +391,4 @@ HighOrder.o: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
   ../Numeric/NumericEmbedded.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/PartitionOptions.h ../Common/GmshMatrix.h \
   ../Numeric/FunctionSpace.h
-Partition.o: Partition.cpp ../Geo/MElement.h \
-  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
-  ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
-  ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/GEdge.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
-  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
-  ../Geo/GRegion.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h Partition.h \
-  ../Common/Message.h PartitionOptions.h
+Partition.o: Partition.cpp
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index ab16b5589fd529a58c670759cdf04b5d00a3fd70..c6f57cc143c8754fd24621474384b208144adc68 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -66,8 +66,8 @@ static void buildInterpLc(const std::vector<IntPoint> &lcPoints)
   IntPoint p;
   interpLc.clear();
   for(int i = 0; i < lcPoints.size(); i++){
-    p=lcPoints[i];
-    interpLc.push_back(xi2lc( p.t, p.lc));
+    p = lcPoints[i];
+    interpLc.push_back(xi2lc(p.t, p.lc));
   }
 }
 
@@ -196,8 +196,9 @@ static double trapezoidal(IntPoint * P1, IntPoint * P2)
   return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
 }
 
-static void RecursiveIntegration(GEdge *ge, IntPoint * from, IntPoint * to,
-				 double (*f) (GEdge *e, double X), std::vector<IntPoint> &Points,
+static void RecursiveIntegration(GEdge *ge, IntPoint *from, IntPoint *to,
+				 double (*f) (GEdge *e, double X), 
+				 std::vector<IntPoint> &Points,
 				 double Prec, int *depth)
 {
   IntPoint P, p1;
@@ -231,7 +232,7 @@ static void RecursiveIntegration(GEdge *ge, IntPoint * from, IntPoint * to,
 
 static double Integration(GEdge *ge, double t1, double t2, 
 			  double (*f) (GEdge *e, double X),
-			  std::vector<IntPoint>&Points, double Prec)
+			  std::vector<IntPoint> &Points, double Prec)
 {
   IntPoint from, to;
 
@@ -279,21 +280,19 @@ void meshGEdge::operator() (GEdge *ge)
 
   Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
 
-  // Create a list of integration points
-  std::vector<IntPoint> Points,lcPoints;
-
   // compute bounds
   Range<double> bounds = ge->parBounds(0);
   double t_begin = bounds.low();
   double t_end = bounds.high();
   
   // first compute the length of the curve by integrating one
+  std::vector<IntPoint> Points;
   double length = Integration(ge, t_begin, t_end, F_One, Points, 1.e-8 * CTX.lc);
   ge->setLength(length);
+  Points.clear();
 
   if(length == 0.0)
     Msg::Debug("Curve %d has a zero length", ge->tag());
-  Points.resize(0); 
 
   // Integrate detJ/lc du 
   double a;
@@ -308,6 +307,7 @@ void meshGEdge::operator() (GEdge *ge)
   }
   else{
     if(CTX.mesh.lc_integration_precision > 1.e-8){
+      std::vector<IntPoint> lcPoints;
       Integration(ge, t_begin, t_end, F_Lc_usingInterpLcBis, lcPoints, 
                   CTX.mesh.lc_integration_precision);
       buildInterpLc(lcPoints);
@@ -317,7 +317,8 @@ void meshGEdge::operator() (GEdge *ge)
       a = Integration(ge, t_begin, t_end, F_Lc_usingInterpLc, Points, 1.e-8);
     }
     else{
-      a = Integration(ge, t_begin, t_end, F_Lc, Points, CTX.mesh.lc_integration_precision);
+      a = Integration(ge, t_begin, t_end, F_Lc, Points,
+		      CTX.mesh.lc_integration_precision);
     }
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
   }
@@ -347,8 +348,8 @@ void meshGEdge::operator() (GEdge *ge)
     IntPoint P1, P2;
     ge->mesh_vertices.resize(N - 2);
     while(NUMP < N - 1) {
-      P1=Points[count-1];
-      P2=Points[count];
+      P1 = Points[count-1];
+      P2 = Points[count];
       const double d = (double)NUMP * b;
       if((fabs(P2.p) >= fabs(d)) && (fabs(P1.p) < fabs(d))) {
         double dt = P2.t - P1.t;
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 91cfc3a52a49e06a8c1e1fd56864e123d544a46e..00b68cd0be43c5c78b3d8b6c4c8afefea1129cf4 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -128,6 +128,16 @@ int PViewDataGModel::getNumTensors(int step)
   return 0;
 }
 
+int PViewDataGModel::getNumPoints(int step)
+{
+  if(_steps.empty()) return 0;
+  GModel *m = _steps[0]->getModel(); // to generalize
+  int n = 0;
+  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it)
+    n += (*it)->points.size();
+  return n;
+}
+
 int PViewDataGModel::getNumLines(int step)
 {
   if(_steps.empty()) return 0;
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 84e6b9f7ffa99a0412ac71cc0e85acabcf403304..b9df6148bdbfdfbb0c292ceaa5a09457351d0e3e 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -44,7 +44,7 @@ class stepData{
     : _model(model), _fileName(fileName), _fileIndex(fileIndex), _time(time), 
       _min(min), _max(max), _numComp(numComp), _data(0)
   {
-    _entities = _model->getEntities();
+    _model->getEntities(_entities);
     _bbox = _model->bounds();
   }
   ~stepData(){ destroyData(); }
@@ -135,7 +135,7 @@ class PViewDataGModel : public PViewData {
   int getNumScalars(int step=-1);
   int getNumVectors(int step=-1);
   int getNumTensors(int step=-1);
-  int getNumPoints(int step=-1){ return 0; }
+  int getNumPoints(int step=-1);
   int getNumLines(int step=-1);
   int getNumTriangles(int step=-1);
   int getNumQuadrangles(int step=-1);
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 0ad64cffd894cd733f7c48529849fe29833fd5b9..2465a05ef14fbb2cc9e2f51842bfa0c4d02bf54a 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -209,20 +209,18 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 
   if(numCompMsh > 9) Msg::Warning("More than 9 components in field");
 
-  // the ordering of the elements in the following lists is important:
-  // it should match the ordering of the MSH element types (when
-  // elements are saved without tags, this governs the order with
-  // which we implicitly index them in GModel::readMED)
+  // Warning! The ordering of the elements in the last two lists is
+  // important: it should match the ordering of the MSH element types
+  // (when elements are saved without tags, this governs the order
+  // with which we implicitly index them in GModel::readMED)
   const med_entite_maillage entType[] = 
     {MED_NOEUD, MED_MAILLE, MED_NOEUD_MAILLE};
-  // don't import points for now (points are not numbered in the same
-  // sequence as MElements)
   const med_geometrie_element eleType[] = 
     {MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8, 
      MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_TETRA10, 
-     /* MED_POINT1, */ MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
+     MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
   const int nodesPerEle[] = 
-    {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 10, /* 1, */ 8, 20, 15, 13};
+    {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 10, 1, 8, 20, 15, 13};
 
   med_int numSteps = 0;
   std::vector<std::pair<int, int> > pairs;
@@ -319,8 +317,8 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	    Msg::Error("Could not read Gauss points");
 	    return false;
 	  }
-	  // we should check that refcoo corresponds to our internal
-	  // reference element
+	  // FIXME: we should check that refcoo corresponds to our
+	  // internal reference element
 	  for(int i = 0; i < (int)gscoo.size(); i++){
 	    p.push_back(gscoo[i]);
 	    if(i % dim == dim - 1) for(int j = 0; j < 3 - dim; j++) p.push_back(0.);