diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 6b2b35e88705cd62a227bcaa688c489f66b96562..5eb4730cb7762b59d528fc0db8184cb60d8e90d1 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -131,7 +131,6 @@ GModel::GModel(std::string name)
   list.push_back(this);
   // at the moment we always create (at least an empty) GEO model
   _createGEOInternals();
-  _newElemNumbers = NULL;
 #if defined(HAVE_OCC)
   _factory = new OCCFactory();
 #endif
@@ -147,7 +146,6 @@ GModel::~GModel()
   destroy();
   _deleteGEOInternals();
   _deleteOCCInternals();
-  if(_newElemNumbers) delete _newElemNumbers;
 #if defined(HAVE_MESH)
   delete _fields;
 #endif
@@ -247,6 +245,7 @@ void GModel::destroyMeshCaches()
   _vertexMapCache.clear();
   _elementVectorCache.clear();
   _elementMapCache.clear();
+  _elementIndexCache.clear();
   Octree_Delete(_octree);
   _octree = 0;
 }
@@ -655,7 +654,6 @@ int GModel::getNumMeshElements(unsigned c[5])
   return 0;
 }
 
-
 MElement *GModel::getMeshElementByCoord(SPoint3 &p)
 {
   if(!_octree){
@@ -723,9 +721,8 @@ void GModel::getMeshVerticesForPhysicalGroup(int dim, int num, std::vector<MVert
   v.insert(v.begin(), sv.begin(), sv.end());
 }
 
-MElement *GModel::getMeshElementByTag(int num)
+MElement *GModel::getMeshElementByTag(int n)
 {
-  int n = (_newElemNumbers) ? (*_newElemNumbers)(num) : num;
   if(_elementVectorCache.empty() && _elementMapCache.empty()){
     Msg::Debug("Rebuilding mesh element cache");
     _elementVectorCache.clear();
@@ -740,16 +737,14 @@ MElement *GModel::getMeshElementByTag(int num)
       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);
-          int ni = (_newElemNumbers) ? (*_newElemNumbers)(e->getNum()) : e->getNum();
-          _elementVectorCache[ni] = e;
+          _elementVectorCache[e->getNum()] = e;
         }
     }
     else{
       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);
-          int ni = (_newElemNumbers) ? (*_newElemNumbers)(e->getNum()) : e->getNum();
-          _elementMapCache[ni] = e;
+          _elementMapCache[e->getNum()] = e;
         }
     }
   }
@@ -760,6 +755,20 @@ MElement *GModel::getMeshElementByTag(int num)
     return _elementMapCache[n];
 }
 
+int GModel::getMeshElementIndex(MElement *e)
+{
+  int num = e->getNum();
+  if(num >= _elementIndexCache.size()) return num;
+  return _elementIndexCache[num];
+}
+
+void GModel::setMeshElementIndex(MElement *e, int index)
+{
+  int num = e->getNum();
+  if(num >= _elementIndexCache.size()) _elementIndexCache.resize(num + 1, 0);
+  _elementIndexCache[num] = index; 
+}
+
 template <class T>
 static void removeInvisible(std::vector<T*> &elements, bool all)
 {
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 6ed3ac83eb6cd030b2512f745666fc3c73d029ba..8cb2b55f11cf669dca747483ccb3e84f19fcccf4 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -17,7 +17,6 @@
 #include "GRegion.h"
 #include "SPoint3.h"
 #include "SBoundingBox3d.h"
-#include "fullMatrix.h"
 
 class Octree;
 class FM_Internals;
@@ -52,7 +51,7 @@ class GModel
   std::map<int, MVertex*> _vertexMapCache;
   std::vector<MElement*> _elementVectorCache;
   std::map<int, MElement*> _elementMapCache;
-  fullVector<double> *_newElemNumbers;
+  std::vector<int> _elementIndexCache;
 
   // ghost cell information (stores partitions for each element acting
   // as a ghost cell)
@@ -302,7 +301,10 @@ class GModel
 
   // access a mesh element by tag, using the element cache
   MElement *getMeshElementByTag(int n);
-  int getMeshElementNumByTag(int n) { return (*_newElemNumbers)(n); }
+
+  // access temporary mesh element index
+  int getMeshElementIndex(MElement *e);
+  void setMeshElementIndex(MElement *e, int index);
 
   // return the total number of vertices in the mesh
   int getNumMeshVertices();
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 8982c70c09ba5f0de3dcff1f34d5c162e7261a1e..46163460db0db4802bb710bd8f6ffd4f7b0af116 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -517,7 +517,7 @@ int GModel::readMSH(const std::string &name)
 template<class T>
 static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll, 
                             double version, bool binary, int &num, int elementary, 
-                            std::vector<int> &physicals, fullVector<double> *newElemNumbers, int parentNum=0)
+                            std::vector<int> &physicals, int parentNum=0)
 {
   std::vector<short> ghosts;
   if(model->getGhostCells().size()){
@@ -536,8 +536,7 @@ static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll,
       ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j],
                     parentNum, &ghosts);
 
-  //(*newElemNumbers)(ele->getNum()) = num;
-
+  model->setMeshElementIndex(ele, num); // should really be a multimap...
 }
 
 template<class T>
@@ -546,7 +545,6 @@ static void writeElementsMSH(FILE *fp, GModel *model, std::vector<T*> &ele,
                              bool saveAll, int saveSinglePartition, double version,
                              bool binary, int &num, int elementary,
                              std::vector<int> &physicals,
-                             fullVector<double> *newElemNumbers,
                              std::map<MElement*, int> *parentsNum=0)
 {
   for(unsigned int i = 0; i < ele.size(); i++){
@@ -561,7 +559,7 @@ static void writeElementsMSH(FILE *fp, GModel *model, std::vector<T*> &ele,
       }
     }
     writeElementMSH(fp, model, ele[i], saveAll, version, binary, num,
-                    elementary, physicals, newElemNumbers, parentNum);
+                    elementary, physicals, parentNum);
   }
 }
 
@@ -697,25 +695,24 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   fprintf(fp, "%d\n", numElements);
   int num = elementStartNum;
   std::map<MElement*, int> parentsNum;
-  if(_newElemNumbers) delete _newElemNumbers;
-  _newElemNumbers = new fullVector<double>(numElements + 1);
+  _elementIndexCache.resize(numElements + 1);
 
   // points
   for(viter it = firstVertex(); it != lastVertex(); ++it)
     writeElementsMSH(fp, this, (*it)->points, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // lines
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin(); 
       it != parents[0].end(); it++)
     if(it->first->getType() == TYPE_LIN) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
     writeElementsMSH(fp, this, (*it)->lines, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers,
+                     version, binary, num, (*it)->tag(), (*it)->physicals, 
                      &parentsNum);
 
   // triangles
@@ -723,36 +720,36 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
       it != parents[1].end(); it++)
     if(it->first->getType() == TYPE_TRI) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); ++it)
     writeElementsMSH(fp, this, (*it)->triangles, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // quads
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin(); 
       it != parents[1].end(); it++)
     if(it->first->getType() == TYPE_QUA) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); ++it)
     writeElementsMSH(fp, this, (*it)->quadrangles, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // polygons
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
       it != parents[1].end(); it++)
     if(it->first->getType() == TYPE_POLYG) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); it++)
     writeElementsMSH(fp, this, (*it)->polygons, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers, 
+                     version, binary, num, (*it)->tag(), (*it)->physicals, 
                      &parentsNum);
 
   // tets
@@ -760,60 +757,60 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
       it != parents[2].end(); it++)
     if(it->first->getType() == TYPE_TET) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, this, (*it)->tetrahedra, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // hexas
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
       it != parents[2].end(); it++)
     if(it->first->getType() == TYPE_HEX) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, this, (*it)->hexahedra, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // prisms
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
       it != parents[2].end(); it++)
     if(it->first->getType() == TYPE_PRI) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, this, (*it)->prisms, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
   
   // pyramids
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
       it != parents[2].end(); it++)
     if(it->first->getType() == TYPE_PYR) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, this, (*it)->pyramids, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers);
+                     version, binary, num, (*it)->tag(), (*it)->physicals);
 
   // polyhedra
   for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
       it != parents[2].end(); it++)
     if(it->first->getType() == TYPE_POLYH) {
       writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
-                      it->second->tag(), it->second->physicals, _newElemNumbers);
+                      it->second->tag(), it->second->physicals);
       parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, this, (*it)->polyhedra, saveAll, saveSinglePartition,
-                     version, binary, num, (*it)->tag(), (*it)->physicals, _newElemNumbers,
+                     version, binary, num, (*it)->tag(), (*it)->physicals,
                      &parentsNum);
 
   if(binary) fprintf(fp, "\n");
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 46dbc30e6e7e447dc3d5b2b1918a92b3460ef168..7aed6e2249bf7b452d0429ccbc403a4bef40e3a5 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -179,13 +179,8 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
               Msg::Error("Unknown element %d in data", i);
               return false;
             }
-            int mult = 1;
-            if(_type == ElementNodeData)
-              mult = e->getNumVertices();
-            // Warning: this will not work if the mesh we just saved
-            // above changed the renumbering of the elements (which is
-            // usually the case)...
-            int num = model->getMeshElementNumByTag(i);
+            int mult = _steps[step]->getMult(i);
+            int num = model->getMeshElementIndex(e);
             if(binary){
               fwrite(&num, sizeof(int), 1, fp);
               if(_type == ElementNodeData)