more work on MSH3

Common/GmshDefines.h

@@ -197,7 +197,6 @@
 #define MSH_PYR_245 131
 // Additional types
 #define MSH_PYR_1   132
-
 #define MSH_PNT_SUB 133
 #define MSH_LIN_SUB 134
 #define MSH_TRI_SUB 135

Geo/GModelIO_MSH.cpp

@@ -129,8 +129,8 @@ int GModel::readMSH(const std::string &name)
       if(sscanf(str, "%d", &numVertices) != 1) return 0;
       Msg::Info("%d vertices", numVertices);
       Msg::ResetProgressMeter();
-      std::map<int, MVertex*> vertexMap;
-      std::vector<MVertex*> vertexVector;
+      _vertexMapCache.clear();
+      _vertexVectorCache.clear();
       int maxVertex = -1;
       minVertex = numVertices + 1;
       for(int i = 0; i < numVertices; i++) {
@@ -210,32 +210,28 @@ int GModel::readMSH(const std::string &name)
         }
         minVertex = std::min(minVertex, num);
         maxVertex = std::max(maxVertex, num);
-        if(vertexMap.count(num))
+        if(_vertexMapCache.count(num))
           Msg::Warning("Skipping duplicate vertex %d", num);
-        vertexMap[num] = vertex;
+        _vertexMapCache[num] = vertex;
         if(numVertices > 100000)
           Msg::ProgressMeter(i + 1, numVertices, true, "Reading nodes");
       }
       // if the vertex numbering is dense, transfer the map into a vector to
       // speed up element creation
-      if((int)vertexMap.size() == numVertices &&
+      if((int)_vertexMapCache.size() == numVertices &&
          ((minVertex == 1 && maxVertex == numVertices) ||
           (minVertex == 0 && maxVertex == numVertices - 1))){
         Msg::Info("Vertex numbering is dense");
-        vertexVector.resize(vertexMap.size() + 1);
+        _vertexVectorCache.resize(_vertexMapCache.size() + 1);
         if(minVertex == 1)
-          vertexVector[0] = 0;
+          _vertexVectorCache[0] = 0;
         else
-          vertexVector[numVertices] = 0;
-        std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
-        for(; it != vertexMap.end(); ++it)
-          vertexVector[it->first] = it->second;
-        vertexMap.clear();
+          _vertexVectorCache[numVertices] = 0;
+        for(std::map<int, MVertex*>::const_iterator it = _vertexMapCache.begin();
+            it != _vertexMapCache.end(); ++it)
+          _vertexVectorCache[it->first] = it->second;
+        _vertexMapCache.clear();
       }
-
-      // cache the vertex indexing data
-      _vertexVectorCache = vertexVector;
-      _vertexMapCache = vertexMap;
     }
 
     // $Elements section
@@ -245,79 +241,53 @@ int GModel::readMSH(const std::string &name)
       if(sscanf(str, "%d", &numElements) != 1) return 0;
       Msg::Info("%d elements", numElements);
       Msg::ResetProgressMeter();
-      std::map<int, MElement*> elementMap;
+      _elementMapCache.clear();
       for(int i = 0; i < numElements; i++) {
-        int num, type, numTags, numVertices;
+        int num, type, entity, numData;
         if(!binary){
-          if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
+          if(fscanf(fp, "%d %d %d %d", &num, &type, &entity, &numData) != 4) return 0;
         }
         else{
           if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
           if(swap) SwapBytes((char*)&num, sizeof(int), 1);
           if(fread(&type, sizeof(int), 1, fp) != 1) return 0;
           if(swap) SwapBytes((char*)&type, sizeof(int), 1);
-          if(fread(&numTags, sizeof(int), 1, fp) != 1) return 0;
-          if(swap) SwapBytes((char*)&numTags, sizeof(int), 1);
+          if(fread(&entity, sizeof(int), 1, fp) != 1) return 0;
+          if(swap) SwapBytes((char*)&entity, sizeof(int), 1);
+          if(fread(&numData, sizeof(int), 1, fp) != 1) return 0;
+          if(swap) SwapBytes((char*)&numData, sizeof(int), 1);
         }
-        std::vector<int> tags;
-        if(numTags > 0){
-          tags.resize(numTags);
+        std::vector<int> data;
+        if(numData > 0){
+          data.resize(numData);
           if(!binary){
-            for(int j = 0; j < numTags; j++){
-              if(fscanf(fp, "%d", &tags[j]) != 1) return 0;
+            for(int j = 0; j < numData; j++){
+              if(fscanf(fp, "%d", &data[j]) != 1) return 0;
             }
           }
           else{
-            if(fread(&tags[0], sizeof(int), numTags, fp) != numTags) return 0;
-            if(swap) SwapBytes((char*)&tags[0], sizeof(int), numTags);
-          }
+            if(fread(&data[0], sizeof(int), numData, fp) != numData) return 0;
+            if(swap) SwapBytes((char*)&data[0], sizeof(int), numData);
           }
-        if(!(numVertices = MElement::getInfoMSH(type))) {
-          return 0;
-        }
-        std::vector<int> indices(numVertices);
-        if(!binary){
-          for(int j = 0; j < numVertices; j++)
-            if(fscanf(fp, "%d", &indices[j]) != 1) return 0;
         }
-        else{
-          if(fread(&indices[0], sizeof(int), numVertices, fp) != numVertices) return 0;
-          if(swap) SwapBytes((char*)&indices[0], sizeof(int), numVertices);
-        }
-        std::vector<MVertex*> vertices;
-        if(_vertexVectorCache.size()){
-          if(!getVertices(numVertices, indices, _vertexVectorCache, vertices, minVertex))
-            return 0;
-        }
-        else{
-          if(!getVertices(numVertices, indices, _vertexMapCache, vertices))
-            return 0;
-        }
-        std::vector<short> ghosts;
         MElementFactory factory;
-        MElement *element = factory.create(num, type, tags, vertices, elementMap,
-                                           ghosts);
-        elementMap[num] = element;
-        int part = element->getPartition();
-        if(part) getMeshPartitions().insert(part);
-        int elementary = tags.size() ? tags[0] : 0;
+        MElement *element = factory.create(num, type, data, this);
         switch(element->getType()){
-        case TYPE_PNT: elements[0][elementary].push_back(element); break;
-        case TYPE_LIN: elements[1][elementary].push_back(element); break;
-        case TYPE_TRI: elements[2][elementary].push_back(element); break;
-        case TYPE_QUA: elements[3][elementary].push_back(element); break;
-        case TYPE_TET: elements[4][elementary].push_back(element); break;
-        case TYPE_HEX: elements[5][elementary].push_back(element); break;
-        case TYPE_PRI: elements[6][elementary].push_back(element); break;
-        case TYPE_PYR: elements[7][elementary].push_back(element); break;
-        }
-        for(unsigned int j = 0; j < ghosts.size(); j++)
-          _ghostCells.insert(std::pair<MElement*, short>(element, ghosts[j]));
+        case TYPE_PNT: elements[0][entity].push_back(element); break;
+        case TYPE_LIN: elements[1][entity].push_back(element); break;
+        case TYPE_TRI: elements[2][entity].push_back(element); break;
+        case TYPE_QUA: elements[3][entity].push_back(element); break;
+        case TYPE_TET: elements[4][entity].push_back(element); break;
+        case TYPE_HEX: elements[5][entity].push_back(element); break;
+        case TYPE_PRI: elements[6][entity].push_back(element); break;
+        case TYPE_PYR: elements[7][entity].push_back(element); break;
+        case TYPE_POLYG: elements[8][entity].push_back(element); break;
+        case TYPE_POLYH: elements[9][entity].push_back(element); break;
+        }
+        _elementMapCache[num] = element;
         if(numElements > 100000)
           Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
       }
-      // cache the element map
-      _elementMapCache = elementMap;
     }
 
     // Post-processing sections

Geo/MElement.cpp

@@ -7,6 +7,7 @@
 #include <math.h>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
+#include "GModel.h"
 #include "MElement.h"
 #include "MPoint.h"
 #include "MLine.h"
@@ -710,7 +711,7 @@ double MElement::integrateFlux(double val[], int face, int pOrder, int order)
   return result;
 }
 
-void MElement::writeMSH(FILE *fp, bool binary, int elementary,
+void MElement::writeMSH(FILE *fp, bool binary, int entity,
                         std::vector<short> *ghosts)
 {
   int num = getNum();
@@ -720,37 +721,40 @@ void MElement::writeMSH(FILE *fp, bool binary, int elementary,
   // if necessary, change the ordering of the vertices to get positive volume
   setVolumePositive();
 
-  std::vector<int> info;
-  info.push_back(0);
-  info.push_back(elementary);
+  std::vector<int> verts;
+  getVerticesIdForMSH(verts);
+
+  // FIXME: once we create elements using their own interpretion of data, we
+  // should move this also into each element base class
+  std::vector<int> data;
+  data.insert(data.end(), verts.begin(), verts.end());
   if(getParent())
-    info.push_back(getParent()->getNum());
+    data.push_back(getParent()->getNum());
   if(getPartition()){
     if(ghosts){
-      info.push_back(1 + ghosts->size());
-      info.push_back(getPartition());
-      info.insert(info.end(), ghosts->begin(), ghosts->end());
+      data.push_back(1 + ghosts->size());
+      data.push_back(getPartition());
+      data.insert(data.end(), ghosts->begin(), ghosts->end());
     }
     else{
-      info.push_back(1);
-      info.push_back(getPartition());
+      data.push_back(1);
+      data.push_back(getPartition());
     }
   }
-  info[0] = info.size() - 1;
-  std::vector<int> verts;
-  getVerticesIdForMSH(verts);
-  info.insert(info.end(), verts.begin(), verts.end());
+  int numData = data.size();
 
   if(!binary){
-    fprintf(fp, "%d %d", num, type);
-    for(unsigned int i = 0; i < info.size(); i++)
-      fprintf(fp, " %d", info[i]);
+    fprintf(fp, "%d %d %d %d", num, type, entity, numData);
+    for(unsigned int i = 0; i < numData; i++)
+      fprintf(fp, " %d", data[i]);
     fprintf(fp, "\n");
   }
   else{
     fwrite(&num, sizeof(int), 1, fp);
     fwrite(&type, sizeof(int), 1, fp);
-    fwrite(&info[0], sizeof(int), info.size(), fp);
+    fwrite(&entity, sizeof(int), 1, fp);
+    fwrite(&numData, sizeof(int), 1, fp);
+    fwrite(&data[0], sizeof(int), numData, fp);
   }
 }
 
@@ -1389,30 +1393,55 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
   }
 }
 
-MElement *MElementFactory::create(int num, int type, const std::vector<int> &tags,
-                                  std::vector<MVertex*> &v,
-                                  std::map<int, MElement*> &elementCache,
-                                  std::vector<short> &ghosts)
+MElement *MElementFactory::create(int num, int type, const std::vector<int> &data,
+                                  GModel *model)
 {
-  MElement *parent = 0;
+  // This should be rewritten: each element should register itself in a static
+  // factory owned e.g. directly by MElement, and interpret its data by
+  // itself. This would remove the ugly switch in the routine above.
+
+  int numVertices = MElement::getInfoMSH(type), startVertices = 0;
+  if(data.size() && !numVertices){
+    startVertices = 1;
+    numVertices = data[0];
+  }
+
+  std::vector<MVertex*> vertices(numVertices);
+  if(data.size() > startVertices + numVertices - 1){
+    for(int i = 0; i < numVertices; i++)
+      vertices[i] = model->getMeshVertexByTag(data[startVertices + i]);
+  }
+  else{
+    return 0;
+  }
+
   int part = 0;
-  if(tags.size() > 2 && (type == MSH_PNT_SUB || type == MSH_LIN_SUB ||
+  int startPartitions = startVertices + numVertices;
+
+  MElement *parent = 0;
+  if((type == MSH_PNT_SUB || type == MSH_LIN_SUB ||
       type == MSH_TRI_SUB || type == MSH_TET_SUB)){
-    parent = elementCache[tags[1]];
-    if(tags.size() > 3 && tags[2]){ // num partitions
-      part = tags[3];
-      for(int i = 0; i < tags[2] - 1; i++)
-        ghosts.push_back(tags[4 + i]);
+    parent = model->getMeshElementByTag(data[startPartitions]);
+    startPartitions += 1;
   }
-  }
-  else if(tags.size() > 1){
-    if(tags[1]){ // num partitions
-      part = tags[2];
-      for(int i = 0; i < tags[1] - 1; i++)
-        ghosts.push_back(tags[3 + i]);
+
+  std::vector<short> ghosts;
+  if(data.size() > startPartitions){
+    int numPartitions = data[startPartitions];
+    if(numPartitions > 0 && data.size() > startPartitions + numPartitions - 1){
+      part = data[startPartitions + 1];
+      for(int i = 1; i < numPartitions; i++)
+        ghosts.push_back(data[startPartitions + 1 + i]);
     }
   }
-  return create(type, v, num, part, false, parent);
+
+  MElement *element = create(type, vertices, num, part, false, parent);
+
+  for(unsigned int j = 0; j < ghosts.size(); j++)
+    model->getGhostCells().insert(std::pair<MElement*, short>(element, ghosts[j]));
+  if(part) model->getMeshPartitions().insert(part);
+
+  return element;
 }
 
 void MElement::xyzTouvw(fullMatrix<double> *xu)

Geo/MElement.h

@@ -18,6 +18,8 @@
 #include "JacobianBasis.h"
 #include "GaussIntegration.h"
 
+class GModel;
+
 // A mesh element.
 class MElement
 {
@@ -357,10 +359,7 @@ class MElementFactory{
  public:
   MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0,
                    bool owner=false, MElement *parent=0, MElement *d1=0, MElement *d2=0);
-  MElement *create(int num, int type, const std::vector<int> &tags,
-                   std::vector<MVertex*> &v,
-                   std::map<int, MElement*> &elementCache,
-                   std::vector<short> &ghosts);
+  MElement *create(int num, int type, const std::vector<int> &data, GModel *model);
 };
 
 // Traits of various elements based on the dimension.  These generally define

Geo/MSubElement.h

@@ -63,7 +63,8 @@ class MSubTriangle : public MTriangle
   MSubTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0,
                bool owner=false, MElement* orig=NULL)
     : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubTriangle(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+  MSubTriangle(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
+               MElement* orig=NULL)
     : MTriangle(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
   ~MSubTriangle();
   virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
@@ -94,9 +95,11 @@ class MSubLine : public MLine
   std::vector<MElement*> _parents;
   IntPt *_intpt;
  public:
-  MSubLine(MVertex *v0, MVertex *v1, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+  MSubLine(MVertex *v0, MVertex *v1, int num=0, int part=0, bool owner=false,
+           MElement* orig=NULL)
     : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubLine(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+  MSubLine(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
+           MElement* orig=NULL)
     : MLine(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
   ~MSubLine();
   virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
@@ -129,7 +132,8 @@ class MSubPoint : public MPoint
  public:
   MSubPoint(MVertex *v0, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
     : MPoint(v0, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubPoint(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+  MSubPoint(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
+            MElement* orig=NULL)
     : MPoint(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
   ~MSubPoint();
   virtual int getTypeForMSH() const { return MSH_PNT_SUB; }