diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 535636c16c963af5fff89fde5cf7f318bc10f3f0..c8fac9805ff88e5eceaf499ecea9242727de7885 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -116,6 +116,8 @@ void GModel::destroyMeshCaches()
 {
   _vertexVectorCache.clear();
   _vertexMapCache.clear();
+  _elementVectorCache.clear();
+  _elementMapCache.clear();
   if(_octree) Octree_Delete(_octree);
   _octree = 0;
 }
@@ -545,6 +547,40 @@ void GModel::getMeshVerticesForPhysicalGroup(int dim, int num, std::vector<MVert
   v.insert(v.begin(), sv.begin(), sv.end());
 }
 
+MElement *GModel::getMeshElementByTag(int n)
+{
+  if(_elementVectorCache.empty() && _elementMapCache.empty()){
+    Msg::Debug("Rebuilding mesh element cache");
+    _elementVectorCache.clear();
+    _elementMapCache.clear();
+    bool dense = (getNumMeshElements() == MElement::getGlobalNumber());
+    std::vector<GEntity*> entities;
+    getEntities(entities);
+    if(dense){
+      Msg::Debug("Good: we have a dense element numbering in the cache");
+      // numbering starts at 1
+      _elementVectorCache.resize(MElement::getGlobalNumber() + 1);
+      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);
+          _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);
+          _elementMapCache[e->getNum()] = e;
+        }
+    }
+  }
+
+  if(n < (int)_elementVectorCache.size())
+    return _elementVectorCache[n];
+  else
+    return _elementMapCache[n];
+}
+
 template <class T>
 static void removeInvisible(std::vector<T*> &elements, bool all)
 {
diff --git a/Geo/GModel.h b/Geo/GModel.h
index cf98e4298bff2a97f1d5d8a6d063202766ea1777..a5969f4792cfe7629d7ddd271afb062f552e5989 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -40,10 +40,12 @@ class GModel
   // the visibility flag
   char _visible;
 
-  // vertex cache to speed-up direct access by vertex number (used for
-  // post-processing I/O)
+  // vertex and element caches to speed-up direct access by tag (only
+  // used for post-processing I/O)
   std::vector<MVertex*> _vertexVectorCache;
   std::map<int, MVertex*> _vertexMapCache;
+  std::vector<MElement*> _elementVectorCache;
+  std::map<int, MElement*> _elementMapCache;
 
   // an octree for fast mesh element lookup
   Octree *_octree;
@@ -256,6 +258,9 @@ class GModel
   // access a mesh element by coordinates (using an octree search)
   MElement *getMeshElementByCoord(SPoint3 &p);
 
+  // access a mesh element by tag, using the element cache
+  MElement *getMeshElementByTag(int n);
+
   // return the total number of vertices in the mesh
   int getNumMeshVertices();
 
diff --git a/Geo/MElement.h b/Geo/MElement.h
index d768db005070d8cc6917e97d46075076d1e53d54..c259731866bae9ac59da2ac560504b0b8b43d07e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -46,6 +46,7 @@ class MElement
   virtual ~MElement(){}
 
   // reset the global node number
+  static int getGlobalNumber(){ return _globalNum; }
   static void resetGlobalNumber(){ _globalNum = 0; }
 
   // set/get the tolerance for isInside() test
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 9e951e587a064b714a6d1b3425dbc0218ca9e532..92e5a926e3f3258f4a3f7944d10de7e64ae505ca 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -157,11 +157,8 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
         fprintf(fp, "$EndNodeData\n");
       }
       else{
-        int mult = 1;
-        if(_type == ElementNodeData){
+        if(_type == ElementNodeData)
           fprintf(fp, "$ElementNodeData\n");
-          mult = 1; // FIXME
-        }
         else
           fprintf(fp, "$ElementData\n");
         fprintf(fp, "1\n\"%s\"\n", getName().c_str());
@@ -169,20 +166,26 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
         fprintf(fp, "3\n%d\n%d\n%d\n", step, numComp, numEnt);
         for(int i = 0; i < _steps[step]->getNumData(); i++){
           if(_steps[step]->getData(i)){
-            // FIXME 
-            // MElement *e = model->getMeshElementByTag(i);
-            // if(!e){
-            //   
-            // }
-            // int num = e->getNum();
-            int num = i;
+            MElement *e = model->getMeshElementByTag(i);
+            if(!e){
+              Msg::Error("Unknown element %d in data", i);
+              return false;
+            }
+            int mult = 1;
+            if(_type == ElementNodeData)
+              mult = e->getNumVertices();
+            int num = e->getNum();
             if(binary){
               fwrite(&num, sizeof(int), 1, fp);
-              fwrite(_steps[step]->getData(i), sizeof(double), numComp, fp);
+              if(_type == ElementNodeData)
+                fwrite(&mult, sizeof(int), 1, fp);
+              fwrite(_steps[step]->getData(i), sizeof(double), numComp * mult, fp);
             }
             else{
               fprintf(fp, "%d", num);
-              for(int k = 0; k < numComp; k++)
+              if(_type == ElementNodeData)
+                fprintf(fp, " %d", mult);
+              for(int k = 0; k < numComp * mult; k++)
                 fprintf(fp, " %.16g", _steps[step]->getData(i)[k]);
               fprintf(fp, "\n");
             }