From 5be5ad6418d020a05999d0a72063df445945be44 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 29 Mar 2008 21:36:30 +0000
Subject: [PATCH] ****************************** *   WARNING - PLEASE TEST    *
 ******************************

- changed I/O so that we NEVER change vertex number ids once a mesh has
been generated (or loaded)

- removed dataIndex from MVertex: stepData in PVewDataGModel now directly
uses the vertex/element number ids. This can lead to increased memory
use in corner cases (for which PViewDataList works fine, so it's not a
big deal), but it is simpler and faster
---
 Geo/GModel.cpp             | 58 ++++++++++++--------------
 Geo/GModel.h               |  9 +---
 Geo/GModelIO_MED.cpp       | 84 +++++++++++++++++++-------------------
 Geo/GModelIO_Mesh.cpp      | 68 +++++-------------------------
 Geo/MElement.cpp           | 14 +++----
 Geo/MElement.h             |  5 +++
 Geo/MVertex.cpp            | 30 +++++++-------
 Geo/MVertex.h              | 20 ++++++---
 Geo/gmshVertex.h           |  2 +-
 Mesh/meshGFaceBDS.cpp      |  4 +-
 Post/OctreePost.cpp        |  8 ++--
 Post/PViewDataGModel.cpp   | 27 ++++++------
 Post/PViewDataGModel.h     | 18 ++++++--
 Post/PViewDataGModelIO.cpp | 81 ++++++++++++------------------------
 14 files changed, 183 insertions(+), 245 deletions(-)

diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 1db423cd45..8a177d5803 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.79 2008-03-29 15:36:02 geuzaine Exp $
+// $Id: GModel.cpp,v 1.80 2008-03-29 21:36:29 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -42,7 +42,7 @@ extern Context_T CTX;
 std::vector<GModel*> GModel::list;
 
 GModel::GModel(std::string name)
-  : _maxVertexDataIndex(-1), _geo_internals(0), _occ_internals(0), _fields(0),
+  : _geo_internals(0), _occ_internals(0), _fields(0),
     modelName(name), normals(0)
 {
   list.push_back(this);
@@ -108,7 +108,6 @@ void GModel::destroy()
 
   _vertexVectorCache.clear();
   _vertexMapCache.clear();
-  _maxVertexDataIndex = -1;
 
   MVertex::resetGlobalNumber();
   MElement::resetGlobalNumber();
@@ -493,82 +492,79 @@ void GModel::removeInvisibleElements()
   }
 }
 
-int GModel::renumberMeshVertices(bool saveAll)
+int GModel::indexMeshVertices(bool all)
 {
-  _vertexVectorCache.clear();
-  _vertexMapCache.clear();
-
   // 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]->setNum(-1);
+      (*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]->setNum(-1);
+      (*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]->setNum(-1);
+      (*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]->setNum(-1);
+      (*it)->mesh_vertices[i]->setIndex(-1);
 
   // tag all mesh vertices belonging to elements that need to be saved
   // with 0
   for(viter it = firstVertex(); it != lastVertex(); ++it)
-    if(saveAll || (*it)->physicals.size()){
+    if(all || (*it)->physicals.size()){
       for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-        (*it)->mesh_vertices[i]->setNum(0);
+        (*it)->mesh_vertices[i]->setIndex(0);
     }
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    if(saveAll || (*it)->physicals.size()){
+    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)->setNum(0);
+          (*it)->lines[i]->getVertex(j)->setIndex(0);
     }
   for(fiter it = firstFace(); it != lastFace(); ++it)
-    if(saveAll || (*it)->physicals.size()){
+    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)->setNum(0);
+          (*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)->setNum(0);
+          (*it)->quadrangles[i]->getVertex(j)->setIndex(0);
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    if(saveAll || (*it)->physicals.size()){
+    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)->setNum(0);
+          (*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)->setNum(0);
+          (*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)->setNum(0);
+          (*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)->setNum(0);
+          (*it)->pyramids[i]->getVertex(j)->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]->getNum())
-        (*it)->mesh_vertices[i]->setNum(++numVertices);
+      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]->getNum())
-        (*it)->mesh_vertices[i]->setNum(++numVertices);
+      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]->getNum())
-        (*it)->mesh_vertices[i]->setNum(++numVertices);
+      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]->getNum())
-        (*it)->mesh_vertices[i]->setNum(++numVertices);
+      if(!(*it)->mesh_vertices[i]->getIndex())
+        (*it)->mesh_vertices[i]->setIndex(++numVertices);
 
   return numVertices;
 }
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 94702b258a..1add954174 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -42,7 +42,6 @@ class GModel
   // post-processing I/O)
   std::vector<MVertex*> _vertexVectorCache;
   std::map<int, MVertex*> _vertexMapCache;
-  int _maxVertexDataIndex;
 
   GEO_Internals *_geo_internals;
   void _createGEOInternals();
@@ -176,12 +175,8 @@ class GModel
   // Access a mesh vertex by tag, using the vertex cache
   MVertex *getMeshVertexByTag(int n);
 
-  // Renumber all the (used) mesh vertices in a continuous sequence
-  int renumberMeshVertices(bool saveAll);
-
-  // Get/set the maximum data index used in mesh vertices
-  int getMaxVertexDataIndex(){ return _maxVertexDataIndex; }
-  void setMaxVertexDataIndex(int n){ _maxVertexDataIndex = n; }
+  // index all the (used) mesh vertices in a continuous sequence
+  int indexMeshVertices(bool all);
 
   // Deletes all invisble mesh elements
   void removeInvisibleElements();
diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index ea531a3181..9d25a2058b 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_MED.cpp,v 1.19 2008-03-29 15:36:02 geuzaine Exp $
+// $Id: GModelIO_MED.cpp,v 1.20 2008-03-29 21:36:29 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -228,25 +228,23 @@ int GModel::readMED(const std::string &name)
 
 template<class T>
 static void fillElementsMED(med_int family, std::vector<T*> &elements, med_int &ele, 
-			    std::vector<med_int> &num, std::vector<med_int> &conn,
-			    std::vector<med_int> &fam, med_geometrie_element &type)
+			    std::vector<med_int> &conn, std::vector<med_int> &fam, 
+			    med_geometrie_element &type)
 {
   for(unsigned int i = 0; i < elements.size(); i++){
-    num.push_back(++ele);
     for(int j = 0; j < elements[i]->getNumVertices(); j++)
-      conn.push_back(elements[i]->getVertexMED(j)->getNum());
+      conn.push_back(elements[i]->getVertexMED(j)->getIndex());
     fam.push_back(family);
     if(!i) getElementTypeForMED(elements[i]->getTypeForMSH(), type);
   }
 }
 
-static void writeElementsMED(med_idt &fid, char *meshName,
-			     std::vector<med_int> &num, std::vector<med_int> &conn, 
+static void writeElementsMED(med_idt &fid, char *meshName, std::vector<med_int> &conn, 
 			     std::vector<med_int> &fam, med_geometrie_element type)
 {
-  if(num.empty()) return;
+  if(fam.empty()) return;
   if(MEDelementsEcr(fid, meshName, (med_int)3, &conn[0], MED_FULL_INTERLACE,
-		    0, MED_FAUX, 0, MED_FAUX, &fam[0], (med_int)num.size(),
+		    0, MED_FAUX, 0, MED_FAUX, &fam[0], (med_int)fam.size(),
 		    MED_MAILLE, type, MED_NOD) < 0)
     Msg(GERROR, "Could not write elements");
 }
@@ -276,11 +274,13 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
   // if there are no physicals we save all the elements
   if(noPhysicalGroups()) saveAll = true;
 
-  // renumber the vertices we will save in a continuous sequence (MED
+  // index the vertices we save in a continuous sequence (MED
   // connectivity is given in terms of vertex indices)
-  renumberMeshVertices(saveAll);
+  indexMeshVertices(saveAll);
 
-  // fill a vector containing all the geometrical entities in the model
+  // fill 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;
   entities.insert(entities.end(), vertices.begin(), vertices.end());
   entities.insert(entities.end(), edges.begin(), edges.end());
@@ -324,20 +324,19 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
   // write the nodes
   {
     std::vector<med_float> coord;
-    std::vector<med_int> num, fam;
+    std::vector<med_int> fam;
     for(unsigned int i = 0; i < entities.size(); i++){
       for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
 	MVertex *v = entities[i]->mesh_vertices[j];
-	if(v->getNum() >= 0){
+	if(v->getIndex() >= 0){
 	  coord.push_back(v->x() * scalingFactor);
 	  coord.push_back(v->y() * scalingFactor);
 	  coord.push_back(v->z() * scalingFactor);
-	  num.push_back(v->getNum());
 	  fam.push_back(0); // we never create node families
 	}
       }
     }
-    if(num.empty()){
+    if(fam.empty()){
       Msg(GERROR, "No nodes to write in MED mesh");
       return 0;
     }
@@ -346,8 +345,8 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
     char coordUnit[3 * MED_TAILLE_PNOM + 1] = 
       "unknown         unknown         unknown         ";
     if(MEDnoeudsEcr(fid, meshName, (med_int)3, &coord[0], MED_FULL_INTERLACE, 
-		    MED_CART, coordName, coordUnit, 0, MED_FAUX, &num[0], 
-		    MED_VRAI, &fam[0], (med_int)num.size()) < 0)
+		    MED_CART, coordName, coordUnit, 0, MED_FAUX, 0, MED_FAUX, 
+		    &fam[0], (med_int)fam.size()) < 0)
       Msg(GERROR, "Could not write nodes");
   }
   
@@ -357,64 +356,63 @@ int GModel::writeMED(const std::string &name, bool saveAll, double scalingFactor
     med_int ele = 0;
 
     { // points
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(viter it = firstVertex(); it != lastVertex(); it++){
 	if(saveAll || (*it)->physicals.size()){
-	  num.push_back(++ele);
-	  conn.push_back((*it)->mesh_vertices[0]->getNum());
+	  conn.push_back((*it)->mesh_vertices[0]->getIndex());
 	  fam.push_back(families[*it]);
 	}
       }
-      writeElementsMED(fid, meshName, num, conn, fam, MED_POINT1);
+      writeElementsMED(fid, meshName, conn, fam, MED_POINT1);
     }
     { // lines
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(eiter it = firstEdge(); it != lastEdge(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->lines, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->lines, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // triangles
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(fiter it = firstFace(); it != lastFace(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->triangles, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->triangles, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // quads
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(fiter it = firstFace(); it != lastFace(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->quadrangles, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->quadrangles, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // tets
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(riter it = firstRegion(); it != lastRegion(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->tetrahedra, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->tetrahedra, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // hexas
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(riter it = firstRegion(); it != lastRegion(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->hexahedra, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->hexahedra, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // prisms
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(riter it = firstRegion(); it != lastRegion(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->prisms, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->prisms, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
     { // pyramids
-      std::vector<med_int> num, conn, fam;
+      std::vector<med_int> conn, fam;
       for(riter it = firstRegion(); it != lastRegion(); it++)
 	if(saveAll || (*it)->physicals.size())
-	  fillElementsMED(families[*it], (*it)->pyramids, ele, num, conn, fam, typ);
-      writeElementsMED(fid, meshName, num, conn, fam, typ);
+	  fillElementsMED(families[*it], (*it)->pyramids, ele, conn, fam, typ);
+      writeElementsMED(fid, meshName, conn, fam, typ);
     }
   }
   
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index f14d2bb905..a16c539e20 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.47 2008-03-25 20:48:32 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.48 2008-03-29 21:36:29 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -412,25 +412,9 @@ 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.
-  bool noElements = true;
-  for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++){
-    noElements &= elements[i].empty();
+  for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
     _storeElementsInEntities(elements[i]);
-  }
 
-  // special case: if there are no elements, create one geometry
-  // vertex for each mesh vertex
-  if(noElements){
-    Msg(INFO, "No elements in mesh: creating geometry vertices");
-    for(unsigned int i = 0; i < vertexVector.size(); i++){
-      MVertex *v = vertexVector[i];
-      if(v) points[v->getNum()].push_back(v);
-    }
-    for(std::map<int, MVertex*>::iterator it = vertexMap.begin(); 
-        it != vertexMap.end(); ++it) 
-      points[it->second->getNum()].push_back(it->second);
-  }
-  
   // treat points separately
   for(std::map<int, std::vector<MVertex*> >::iterator it = points.begin(); 
       it != points.end(); ++it){
@@ -533,9 +517,9 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   // if there are no physicals we save all the elements
   if(noPhysicalGroups()) saveAll = true;
 
-  // get the number of vertices and renumber the vertices in a
-  // continuous sequence
-  int numVertices = renumberMeshVertices(saveAll);
+  // get the number of vertices and index the vertices in a continuous
+  // sequence
+  int numVertices = indexMeshVertices(saveAll);
   
   // binary format exists only in version 2
   if(binary) version = 2.0;
@@ -672,38 +656,6 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
     fprintf(fp, "$ENDELM\n");
   }
 
-#if 0 // test NodeData
-  std::vector<MVertex*> allVertices;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      allVertices.push_back((*it)->mesh_vertices[i]);
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      allVertices.push_back((*it)->mesh_vertices[i]);
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      allVertices.push_back((*it)->mesh_vertices[i]);
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      allVertices.push_back((*it)->mesh_vertices[i]);
-  fprintf(fp, "$NodeData\n");
-  fprintf(fp, "\"test\"\n");
-  fprintf(fp, "0 0 0 0 1 %d\n", allVertices.size());
-  for(unsigned int i = 0; i < allVertices.size(); i++){
-    int tag = allVertices[i]->getNum();
-    double val = allVertices[i]->x() * allVertices[i]->y();
-    if(binary){
-      fwrite(&tag, sizeof(int), 1, fp);
-      fwrite(&val, sizeof(double), 1, fp);
-    }
-    else{
-      fprintf(fp, "%d %.16g\n", tag, val);
-    }
-  }
-  if(binary) fprintf(fp, "\n");
-  fprintf(fp, "$EndNodeData\n");
-#endif
-
   fclose(fp);
   return 1;
 }
@@ -1130,7 +1082,7 @@ int GModel::writeVRML(const std::string &name, bool saveAll, double scalingFacto
 
   if(noPhysicalGroups()) saveAll = true;
 
-  renumberMeshVertices(saveAll);
+  indexMeshVertices(saveAll);
 
   fprintf(fp, "#VRML V1.0 ascii\n");
   fprintf(fp, "#created by Gmsh\n");
@@ -1348,7 +1300,7 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
 
   if(noPhysicalGroups()) saveAll = true;
 
-  renumberMeshVertices(saveAll);
+  indexMeshVertices(saveAll);
 
   // nodes
   fprintf(fp, "%6d\n", -1);
@@ -1437,7 +1389,7 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
             fprintf(fp, "\n");
             row = 0;
           }
-          fprintf(fp, "%10d%10d%10d%10d", 7, (*it2)->getNum(), 0, 0);
+          fprintf(fp, "%10d%10d%10d%10d", 7, (*it2)->getIndex(), 0, 0);
           row++;
         }
         fprintf(fp, "\n");
@@ -1561,7 +1513,7 @@ int GModel::writeMESH(const std::string &name, bool saveAll, double scalingFacto
 
   if(noPhysicalGroups()) saveAll = true;
 
-  int numVertices = renumberMeshVertices(saveAll);
+  int numVertices = indexMeshVertices(saveAll);
 
   fprintf(fp, " MeshVersionFormatted 1\n");
   fprintf(fp, " Dimension\n");
@@ -1906,7 +1858,7 @@ int GModel::writeBDF(const std::string &name, int format, bool saveAll,
 
   if(noPhysicalGroups()) saveAll = true;
 
-  renumberMeshVertices(saveAll);
+  indexMeshVertices(saveAll);
 
   fprintf(fp, "$ Created by Gmsh\n");
 
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 98d6d7b279..38ac7a2fd5 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.62 2008-03-23 21:42:57 geuzaine Exp $
+// $Id: MElement.cpp,v 1.63 2008-03-29 21:36:29 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -411,7 +411,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
 
   int verts[30];
   for(int i = 0; i < n; i++)
-    verts[i] = getVertex(i)->getNum();
+    verts[i] = getVertex(i)->getIndex();
 
   if(!binary){
     for(int i = 0; i < n; i++)
@@ -531,7 +531,7 @@ void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
 void MElement::writeVRML(FILE *fp)
 {
   for(int i = 0; i < getNumVertices(); i++)
-    fprintf(fp, "%d,", getVertex(i)->getNum() - 1);
+    fprintf(fp, "%d,", getVertex(i)->getIndex() - 1);
   fprintf(fp, "-1,\n");
 }
 
@@ -553,7 +553,7 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
   if(physical < 0) revert();
 
   for(int k = 0; k < n; k++) {
-    fprintf(fp, "%10d", getVertexUNV(k)->getNum());
+    fprintf(fp, "%10d", getVertexUNV(k)->getIndex());
     if(k % 8 == 7)
       fprintf(fp, "\n");
   }
@@ -566,7 +566,7 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
 void MElement::writeMESH(FILE *fp, int elementary)
 {
   for(int i = 0; i < getNumVertices(); i++)
-    fprintf(fp, " %d", getVertex(i)->getNum());
+    fprintf(fp, " %d", getVertex(i)->getIndex());
   fprintf(fp, " %d\n", elementary);
 }
 
@@ -583,7 +583,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   if(format == 0){ // free field format
     fprintf(fp, "%s,%d,%d", str, _num, elementary);
     for(int i = 0; i < n; i++){
-      fprintf(fp, ",%d", getVertex(i)->getNum());
+      fprintf(fp, ",%d", getVertex(i)->getIndex());
       if(i != n - 1 && !((i + 3) % 8)){
         fprintf(fp, ",+%s%d\n+%s%d", cont[ncont], _num, cont[ncont], _num);
         ncont++;
@@ -596,7 +596,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   else{ // small or large field format
     fprintf(fp, "%-8s%-8d%-8d", str, _num, elementary);
     for(int i = 0; i < n; i++){
-      fprintf(fp, "%-8d", getVertex(i)->getNum());
+      fprintf(fp, "%-8d", getVertex(i)->getIndex());
       if(i != n - 1 && !((i + 3) % 8)){
         fprintf(fp, "+%s%-6d\n+%s%-6d", cont[ncont], _num, cont[ncont], _num);
         ncont++;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index f0896b2418..90b654f9c3 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -34,9 +34,14 @@ class GFace;
 class MElement 
 {
  private:
+  // the maximum element id number in the mesh
   static int _globalNum;
+  // the id number of the element (this number is unique and is
+  // guaranteed never to change once a mesh has been generated)
   int _num;
+  // the number of the mesh partition the element belongs to
   short _partition;
+  // a visibility flag
   char _visible;
  protected:
   void _getEdgeRep(MVertex *v0, MVertex *v1, 
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 503707c27a..47e0964efb 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.22 2008-03-20 11:44:06 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.23 2008-03-29 21:36:29 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -46,16 +46,16 @@ bool MVertexLessThanLexicographic::operator()(const MVertex *v1, const MVertex *
 
 void MVertex::writeMSH(FILE *fp, bool binary, double scalingFactor)
 {
-  if(_num < 0) return; // negative vertices are never saved
+  if(_index < 0) return; // negative index vertices are never saved
 
   if(!binary){
-    fprintf(fp, "%d %.16g %.16g %.16g\n", _num, 
+    fprintf(fp, "%d %.16g %.16g %.16g\n", _index, 
             x() * scalingFactor, 
             y() * scalingFactor,
             z() * scalingFactor);
   }
   else{
-    fwrite(&_num, sizeof(int), 1, fp);
+    fwrite(&_index, sizeof(int), 1, fp);
     double data[3] = {x() * scalingFactor, y() * scalingFactor, z() * scalingFactor};
     fwrite(data, sizeof(double), 3, fp);
   }
@@ -70,19 +70,19 @@ void MVertex::writeMSH(FILE *fp, double version, bool binary, int num,
       fprintf(fp, " %d %d 1", physical, elementary);
     else
       fprintf(fp, " 2 %d %d", physical, elementary);
-    fprintf(fp, " %d\n", _num);
+    fprintf(fp, " %d\n", _index);
   }
   else{
     int tags[4] = {num, physical, elementary, 0};
     fwrite(tags, sizeof(int), 4, fp);
-    int verts[1] = {_num};
+    int verts[1] = {_index};
     fwrite(verts, sizeof(int), 1, fp);
   }
 }
 
 void MVertex::writeVRML(FILE *fp, double scalingFactor)
 {
-  if(_num < 0) return; // negative vertices are never saved
+  if(_index < 0) return; // negative index vertices are never saved
 
   fprintf(fp, "%.16g %.16g %.16g,\n",
           x() * scalingFactor, y() * scalingFactor, z() * scalingFactor);
@@ -90,12 +90,12 @@ void MVertex::writeVRML(FILE *fp, double scalingFactor)
 
 void MVertex::writeUNV(FILE *fp, double scalingFactor)
 {
-  if(_num < 0) return; // negative vertices are never saved
+  if(_index < 0) return; // negative index vertices are never saved
 
   int coord_sys = 1;
   int displacement_coord_sys = 1;
   int color = 11;
-  fprintf(fp, "%10d%10d%10d%10d\n", _num, coord_sys, displacement_coord_sys, color);
+  fprintf(fp, "%10d%10d%10d%10d\n", _index, coord_sys, displacement_coord_sys, color);
   // hack to print the numbers with "D+XX" exponents
   char tmp[128];
   sprintf(tmp, "%25.16E%25.16E%25.16E\n", x() * scalingFactor, 
@@ -106,7 +106,7 @@ void MVertex::writeUNV(FILE *fp, double scalingFactor)
 
 void MVertex::writeMESH(FILE *fp, double scalingFactor)
 {
-  if(_num < 0) return; // negative vertices are never saved
+  if(_index < 0) return; // negative index vertices are never saved
 
   fprintf(fp, " %20.14G      %20.14G      %20.14G      %d\n", 
           x() * scalingFactor, y() * scalingFactor, z() * scalingFactor, 0);
@@ -141,7 +141,7 @@ static void double_to_char8(double val, char *str)
 
 void MVertex::writeBDF(FILE *fp, int format, double scalingFactor)
 {
-  if(_num < 0) return; // negative vertices are never saved
+  if(_index < 0) return; // negative index vertices are never saved
 
   char xs[17], ys[17], zs[17];
   double x1 = x() * scalingFactor;
@@ -150,17 +150,17 @@ void MVertex::writeBDF(FILE *fp, int format, double scalingFactor)
   if(format == 0){
     // free field format (max 8 char per field, comma separated, 10 per line)
     double_to_char8(x1, xs); double_to_char8(y1, ys); double_to_char8(z1, zs);
-    fprintf(fp, "GRID,%d,%d,%s,%s,%s\n", _num, 0, xs, ys, zs);
+    fprintf(fp, "GRID,%d,%d,%s,%s,%s\n", _index, 0, xs, ys, zs);
   }
   else if(format == 1){ 
     // small field format (8 char par field, 10 per line)
     double_to_char8(x1, xs); double_to_char8(y1, ys); double_to_char8(z1, zs);
-    fprintf(fp, "GRID    %-8d%-8d%-8s%-8s%-8s\n", _num, 0, xs, ys, zs);
+    fprintf(fp, "GRID    %-8d%-8d%-8s%-8s%-8s\n", _index, 0, xs, ys, zs);
   }
   else{ 
     // large field format (8 char first/last field, 16 char middle, 6 per line)
-    fprintf(fp, "GRID*   %-16d%-16d%-16.9G%-16.9G*N%-6d\n", _num, 0, x1, y1, _num);
-    fprintf(fp, "*N%-6d%-16.9G\n", _num, z1);
+    fprintf(fp, "GRID*   %-16d%-16d%-16.9G%-16.9G*N%-6d\n", _index, 0, x1, y1, _index);
+    fprintf(fp, "*N%-6d%-16.9G\n", _index, z1);
   }
 }
 
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 142f03e887..1fcc9375f4 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -37,16 +37,25 @@ class MVertexLessThanLexicographic{
 // A mesh vertex.
 class MVertex{
  private:
+  // the maximum vertex id number in the mesh
   static int _globalNum;
+  // the id number of the vertex (this number is unique and is
+  // guaranteed never to change once a mesh has been generated)
   int _num;
+  // a vertex index, used for example when saving a mesh (this index
+  // is not necessarily unique, can change after mesh renumbering,
+  // etc.)
+  int _index;
+  // a visibility and polynomial order flags
   char _visible, _order;
+  // the cartesian coordinates of the vertex
   double _x, _y, _z;
+  // the geometrical entity the vertex is associated with
   GEntity *_ge;
-  int _dataIndex;
 
  public :
   MVertex(double x, double y, double z, GEntity *ge=0, int num=0) 
-    : _visible(true), _order(1), _x(x), _y(y), _z(z), _ge(ge), _dataIndex(-1)
+    : _visible(true), _order(1), _x(x), _y(y), _z(z), _ge(ge)
   {
     if(num){
       _num = num;
@@ -55,6 +64,7 @@ class MVertex{
     else{
       _num = ++_globalNum;
     }
+    _index = num;
   }
   virtual ~MVertex(){}
 
@@ -87,9 +97,9 @@ class MVertex{
   inline int getNum() const { return _num; }
   inline void setNum(int num) { _num = num; }
 
-  // get/set the data index
-  inline int getDataIndex() { return _dataIndex; }
-  inline void setDataIndex(int index) { _dataIndex = index; }
+  // get/set the index
+  inline int getIndex() const { return _index; }
+  inline void setIndex(int index) { _index = index; }
 
   // get/set ith parameter
   virtual bool getParameter(int i, double &par) const{ return false; }
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
index 43ed184f50..2c4e0240b5 100644
--- a/Geo/gmshVertex.h
+++ b/Geo/gmshVertex.h
@@ -29,7 +29,7 @@ class gmshVertex : public GVertex {
   Vertex *v;
 
  public:
-  gmshVertex(GModel *m, Vertex *_v) : GVertex(m, _v->Num,_v->lc), v(_v)
+  gmshVertex(GModel *m, Vertex *_v) : GVertex(m, _v->Num, _v->lc), v(_v)
   {
     mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
   }
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 18e216bb89..603bd9afa7 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.10 2008-03-20 11:44:08 geuzaine Exp $
+// $Id: meshGFaceBDS.cpp,v 1.11 2008-03-29 21:36:30 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -802,7 +802,7 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
 void gmshCollapseSmallEdges(GModel &gm)
 {
   return;
-  gm.renumberMeshVertices(true);
+  // gm.renumberMeshVertices(true);
   std::list<GFace*> faces;
   for (GModel::fiter fit = gm.firstFace(); fit != gm.lastFace(); fit++){
     faces.push_back(*fit);
diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp
index 743b63907c..855bd24514 100644
--- a/Post/OctreePost.cpp
+++ b/Post/OctreePost.cpp
@@ -1,4 +1,4 @@
-// $Id: OctreePost.cpp,v 1.8 2008-03-25 20:48:32 geuzaine Exp $
+// $Id: OctreePost.cpp,v 1.9 2008-03-29 21:36:30 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -420,10 +420,8 @@ bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep, doub
   MElement *e = (MElement*)in;
 
   int dataIndex[8];
-  for(int i = 0; i < e->getNumVertices(); i++){
-    dataIndex[i] = e->getVertex(i)->getDataIndex();
-    if(dataIndex[i] < 0) return false; // no data
-  }
+  for(int i = 0; i < e->getNumVertices(); i++)
+    dataIndex[i] = e->getVertex(i)->getNum();
   
   double U[3];
   e->xyz2uvw(P, U);
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 288ad42fe9..cc43db5930 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModel.cpp,v 1.36 2008-03-29 15:36:02 geuzaine Exp $
+// $Id: PViewDataGModel.cpp,v 1.37 2008-03-29 21:36:30 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -130,17 +130,20 @@ int PViewDataGModel::getNumElements(int step, int ent)
 
 int PViewDataGModel::getDimension(int step, int ent, int ele)
 {
+  // no sanity checks (assumed to be guarded by skipElement)
   return _steps[step]->getEntity(ent)->getMeshElement(ele)->getDim();
 }
 
 int PViewDataGModel::getNumNodes(int step, int ent, int ele)
 {
+  // no sanity checks (assumed to be guarded by skipElement)
   return _steps[step]->getEntity(ent)->getMeshElement(ele)->getNumVertices();
 }
 
 void PViewDataGModel::getNode(int step, int ent, int ele, int nod, 
                               double &x, double &y, double &z)
 {
+  // no sanity checks (assumed to be guarded by skipElement)
   MVertex *v = _steps[step]->getEntity(ent)->getMeshElement(ele)->getVertex(nod);
   x = v->x();
   y = v->y();
@@ -149,18 +152,20 @@ void PViewDataGModel::getNode(int step, int ent, int ele, int nod,
 
 int PViewDataGModel::getNumComponents(int step, int ent, int ele)
 {
+  // no sanity checks (assumed to be guarded by skipElement)
   return _steps[step]->getNumComponents();
 }
 
 void PViewDataGModel::getValue(int step, int ent, int ele, int nod, int comp, double &val)
 {
+  // no sanity checks (assumed to be guarded by skipElement)
   MVertex *v = _steps[step]->getEntity(ent)->getMeshElement(ele)->getVertex(nod);
-  int index = v->getDataIndex();
-  val = _steps[step]->getData(index)[comp];
+  val = _steps[step]->getData(v->getNum())[comp];
 }
 
 int PViewDataGModel::getNumEdges(int step, int ent, int ele)
 { 
+  // no sanity checks (assumed to be guarded by skipElement)
   return _steps[step]->getEntity(ent)->getMeshElement(ele)->getNumEdges();
 }
 
@@ -171,20 +176,19 @@ bool PViewDataGModel::skipEntity(int step, int ent)
 
 bool PViewDataGModel::skipElement(int step, int ent, int ele)
 {
-  if(step >= (int)_steps.size() || !_steps[step]->getNumData()) return true;
-  MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
+  if(step >= getNumTimeSteps()) return true;
+  stepData<double> *data = _steps[step];
+  if(!_steps[step]->getNumData()) return true;
+  MElement *e = data->getEntity(ent)->getMeshElement(ele);
   if(!e->getVisibility()) return true;
-  for(int i = 0; i < e->getNumVertices(); i++){
-    int index = e->getVertex(i)->getDataIndex();
-    if(index < 0 || index >= (int)_steps[step]->getNumData()) return true;
-    if(!_steps[step]->getData(index)) return true;
-  }
+  for(int i = 0; i < e->getNumVertices(); i++)
+    if(!data->getData(e->getVertex(i)->getNum())) return true;
   return false;
 }
 
 bool PViewDataGModel::hasTimeStep(int step)
 {
-  if(step < (int)_steps.size() && _steps[step]->getNumData()) return true;
+  if(step < getNumTimeSteps() && _steps[step]->getNumData()) return true;
   return false;
 }
 
@@ -209,7 +213,6 @@ GEntity *PViewDataGModel::getEntity(int step, int ent)
 
 bool PViewDataGModel::getValue(int step, int dataIndex, int comp, double &val)
 {
-  if(dataIndex < 0 || dataIndex >= (int)_steps[step]->getNumData()) return false;
   double *d = _steps[step]->getData(dataIndex);
   if(!d) return false;
   val = d[comp];
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 5a4c811c68..7c2b7d2a12 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -51,7 +51,12 @@ class stepData{
   // the number of components in the data (one stepData contains only
   // a single field type)
   int _numComp;
-  // the values, indexed by dataIndex in MVertex or MElement
+  // the values, indexed by MVertex or MElement id numbers (If the
+  // numbering is sparse, or if we only have data for high-id
+  // entities, the vector has zero entries and is thus not
+  // optimal. This is the price to pay if we want 1) rapid access to
+  // the data and 2) not to store any additional info in MVertex or
+  // MElement)
   std::vector<real*> *_data;
  public:
   stepData(GModel *model, DataType type, int numComp, 
@@ -99,8 +104,13 @@ class stepData{
   }
   real *getData(int index, bool allocIfNeeded=false)
   {
-    if(!_data || index >= (int)_data->size()) resizeData(index + 100); // optimize this
-    if(allocIfNeeded && !(*_data)[index]) (*_data)[index] = new real[_numComp];
+    if(allocIfNeeded){
+      if(index >= getNumData()) resizeData(index + 100); // optimize this
+      if(!(*_data)[index]) (*_data)[index] = new real[_numComp];
+    }
+    else{
+      if(index >= getNumData()) return 0;
+    }
     return (*_data)[index];
   }
   void destroyData()
@@ -154,7 +164,7 @@ class PViewDataGModel : public PViewData {
 
   // direct access to GModel entities
   GEntity *getEntity(int step, int ent);
-  // direct access to value by dataIndex
+  // direct access to value by index
   bool getValue(int step, int dataIndex, int comp, double &val);
 
   // I/O routines
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 7f30032a1d..01035b21c4 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModelIO.cpp,v 1.16 2008-03-29 15:36:02 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.17 2008-03-29 21:36:30 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -61,18 +61,7 @@ bool PViewDataGModel::readMSH(std::string fileName, int fileIndex, FILE *fp,
     else{
       if(fscanf(fp, "%d", &num) != 1) return false;
     }
-    MVertex *v = _steps[step]->getModel()->getMeshVertexByTag(num);
-    if(!v){
-      Msg(GERROR, "Unknown vertex %d in data", num);
-      return false;
-    }
-    if(v->getDataIndex() < 0){
-      int max = _steps[step]->getModel()->getMaxVertexDataIndex();
-      _steps[step]->getModel()->setMaxVertexDataIndex(max + 1);
-      v->setDataIndex(max + 1);
-    }
-    int index = v->getDataIndex();
-    double *d = _steps[step]->getData(index, true);
+    double *d = _steps[step]->getData(num, true);
     if(binary){
       if((int)fread(d, sizeof(double), numComp, fp) != numComp) return false;
       if(swap) swapBytes((char*)d, sizeof(double), numComp);
@@ -114,15 +103,6 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
     return false;
   }
 
-  // map data index to vertex tags
-  std::vector<int> tags(model->getMaxVertexDataIndex() + 1, 0);
-  for(int i = 0; i < _steps[0]->getNumEntities(); i++){
-    for(unsigned int j = 0; j < _steps[0]->getEntity(i)->mesh_vertices.size(); j++){
-      MVertex *v = _steps[0]->getEntity(i)->mesh_vertices[j];
-      if(v->getDataIndex() >= 0) tags[v->getDataIndex()] = v->getNum();
-    }
-  }
-
   for(unsigned int step = 0; step < _steps.size(); step++){
     int numNodes = 0, numComp = _steps[step]->getNumComponents();
     for(int i = 0; i < _steps[step]->getNumData(); i++)
@@ -134,12 +114,18 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
               numComp, numNodes);
       for(int i = 0; i < _steps[step]->getNumData(); i++){
         if(_steps[step]->getData(i)){
+	  MVertex *v = _steps[step]->getModel()->getMeshVertexByTag(i);
+	  if(!v){
+	    Msg(GERROR, "Unknown vertex %d in data", i);
+	    return false;
+	  }
+	  int num = v->getIndex();
           if(binary){
-            fwrite(&tags[i], sizeof(int), 1, fp);
+            fwrite(&num, sizeof(int), 1, fp);
             fwrite(_steps[step]->getData(i), sizeof(double), numComp, fp);
           }
           else{
-            fprintf(fp, "%d", tags[i]);
+            fprintf(fp, "%d", num);
             for(int k = 0; k < numComp; k++)
               fprintf(fp, " %.16g", _steps[step]->getData(i)[k]);
             fprintf(fp, "\n");
@@ -294,17 +280,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	  nodeTags.clear();
 	for(unsigned int i = 0; i < profile.size(); i++){
 	  int num = nodeTags.empty() ? profile[i] : nodeTags[profile[i] - 1];
-	  MVertex *v = _steps[step]->getModel()->getMeshVertexByTag(num);
-	  if(!v){
-	    Msg(GERROR, "Unknown vertex %d in data", num);
-	    return false;
-	  }
-	  if(v->getDataIndex() < 0){
-	    int max = _steps[step]->getModel()->getMaxVertexDataIndex();
-	    _steps[step]->getModel()->setMaxVertexDataIndex(max + 1);
-	    v->setDataIndex(max + 1);
-	  }
-	  fillData(_steps[step], v->getDataIndex(), val, i, numComp);
+	  fillData(_steps[step], num, val, profile[i] - 1, numComp);
 	}
       }
       else{
@@ -346,23 +322,21 @@ bool PViewDataGModel::writeMED(std::string fileName)
     return false;
   }
 
-  // map data index to vertex tags
-  std::vector<int> tags(model->getMaxVertexDataIndex() + 1, 0);
-  for(int i = 0; i < _steps[0]->getNumEntities(); i++){
-    for(unsigned int j = 0; j < _steps[0]->getEntity(i)->mesh_vertices.size(); j++){
-      MVertex *v = _steps[0]->getEntity(i)->mesh_vertices[j];
-      if(v->getDataIndex() >= 0) tags[v->getDataIndex()] = v->getNum();
-    }
-  }
-
   // compute profile
-  std::vector<med_int> prof;
+  std::vector<med_int> profile, nums;
   for(int i = 0; i < _steps[0]->getNumData(); i++){
-    if(_steps[0]->getData(i))
-      prof.push_back(tags[i]);
+    if(_steps[0]->getData(i)){
+      MVertex *v = _steps[0]->getModel()->getMeshVertexByTag(i);
+      if(!v){
+	Msg(GERROR, "Unknown vertex %d in data", i);
+	return false;
+      }
+      profile.push_back(v->getIndex());
+      nums.push_back(i);
+    }
   }
   char *profileName = (char*)"nodeProfile";
-  if(MEDprofilEcr(fid, &prof[0], (med_int)prof.size(), profileName) < 0){
+  if(MEDprofilEcr(fid, &profile[0], (med_int)profile.size(), profileName) < 0){
     Msg(GERROR, "Could not create MED profile");
     return false;
   }
@@ -385,18 +359,15 @@ bool PViewDataGModel::writeMED(std::string fileName)
     unsigned int n = 0;
     for(int i = 0; i < _steps[step]->getNumData(); i++)
       if(_steps[step]->getData(i)) n++;
-    if(n != prof.size() || numComp != _steps[step]->getNumComponents()){
+    if(n != profile.size() || numComp != _steps[step]->getNumComponents()){
       Msg(GERROR, "Skipping incompatible step");
       continue;
     }
     double time = _steps[step]->getTime();
     std::vector<double> val(numNodes * numComp);
-    for(int i = 0; i < _steps[step]->getNumData(); i++){
-      if(_steps[step]->getData(i)){
-	for(int k = 0; k < numComp; k++)
-	  val[i * numComp + k] = _steps[step]->getData(i)[k];
-      }
-    }
+    for(unsigned int i = 0; i < profile.size(); i++)
+      for(int k = 0; k < numComp; k++)
+	val[(profile[i] - 1) * numComp + k] = _steps[step]->getData(nums[i])[k];
     if(MEDchampEcr(fid, meshName, fieldName, (unsigned char*)&val[0], 
 		   MED_FULL_INTERLACE, numNodes, MED_NOGAUSS, MED_ALL,
 		   profileName, MED_COMPACT, MED_NOEUD, MED_NONE, (med_int)step,
-- 
GitLab