diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 8770d1b851ed929e06e3cddd3ca0907946cfc867..89f1a4e77622aa1e2ec66651d81dd4596ee0f607 100644
--- a/Common/GmshDefines.h
+++ b/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
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 73232135000dca622c650db97de970e1068ca60a..aa1210406db94c14df499eca95540312f02b93aa 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/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++) {
@@ -157,7 +157,7 @@ int GModel::readMSH(const std::string &name)
case 0:
{
GVertex *gv = getVertexByTag(entity);
- if (gv) gv->deleteMesh();
+ if(gv) gv->deleteMesh();
vertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
}
break;
@@ -166,7 +166,7 @@ int GModel::readMSH(const std::string &name)
GEdge *ge = getEdgeByTag(entity);
double u;
if(!binary){
- if (fscanf(fp, "%lf", &u) != 1) return 0;
+ if(fscanf(fp, "%lf", &u) != 1) return 0;
}
else{
if(fread(&u, sizeof(double), 1, fp) != 1) return 0;
@@ -180,7 +180,7 @@ int GModel::readMSH(const std::string &name)
GFace *gf = getFaceByTag(entity);
double uv[2];
if(!binary){
- if (fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) return 0;
+ if(fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) return 0;
}
else{
if(fread(uv, sizeof(double), 2, fp) != 2) return 0;
@@ -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;
+ 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;
}
- for(unsigned int j = 0; j < ghosts.size(); j++)
- _ghostCells.insert(std::pair<MElement*, short>(element, ghosts[j]));
+ _elementMapCache[num] = element;
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
}
- // cache the element map
- _elementMapCache = elementMap;
}
// Post-processing sections
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index c1e8c2282d1f5a6618926d64a1a5721d5a82d537..18e6dc3b0d7f92a5a07a7a12594c943883a01adf 100644
--- a/Geo/MElement.cpp
+++ b/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;
- int part = 0;
- if(tags.size() > 2 && (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]);
- }
+ // 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 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]);
+ else{
+ return 0;
+ }
+
+ int part = 0;
+ 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 = model->getMeshElementByTag(data[startPartitions]);
+ startPartitions += 1;
+ }
+
+ 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)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7630b6664d4e0a5a298284b15e630c8fd94abcbf..67d48a6ca1d891c3923d1950a89f5a242d0c20b4 100644
--- a/Geo/MElement.h
+++ b/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
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index 4f241d1facb678e4976873714dc7c45b58566f5a..53ba9bb60b63272a8c84e09bda06b3ce220cdb28 100644
--- a/Geo/MSubElement.h
+++ b/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; }