From 4c0f47654fdea3a82e7e1de958f450b6aadfd0c5 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 14 Sep 2009 10:36:51 +0000
Subject: [PATCH] added GModel::getMeshElementByTag (uses a dynamically created
 cache, like what we do for mesh vertices)

---
 Geo/GModel.cpp             | 36 ++++++++++++++++++++++++++++++++++++
 Geo/GModel.h               |  9 +++++++--
 Geo/MElement.h             |  1 +
 Post/PViewDataGModelIO.cpp | 29 ++++++++++++++++-------------
 4 files changed, 60 insertions(+), 15 deletions(-)

diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 535636c16c..c8fac9805f 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 cf98e4298b..a5969f4792 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 d768db0050..c259731866 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 9e951e587a..92e5a926e3 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");
             }
-- 
GitLab