diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2811157548c3e93532c8d78237af09659ca1938b..535bb6112f5507b72a099ee90f0b1b944a2709f3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -81,7 +81,6 @@ set(GMSH_API
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
Solver/linearSystemFull.h Solver/elasticitySolver.h
-Solver/c0DgPlateSolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index b0c0db81d442b4de2caad31494b621a978241cc2..eff51558df748c26d7e1052f5ee6b7189ba3bd8e 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -350,6 +350,7 @@ int GModel::readMSH(const std::string &name)
if(!fgets(str, sizeof(str), fp)) return 0;
int numElements;
std::map<int, MElement*> parents;
+ std::set<MElement*> parentsOwned;
sscanf(str, "%d", &numElements);
Msg::Info("%d elements", numElements);
Msg::ResetProgressMeter();
@@ -367,7 +368,6 @@ int GModel::readMSH(const std::string &name)
else{
int numTags;
if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
- std::vector<int> tags(numTags);
int numPartitions = 0;
for(int j = 0; j < numTags; j++){
int tag;
@@ -375,9 +375,9 @@ int GModel::readMSH(const std::string &name)
if(j == 0) physical = tag;
else if(j == 1) elementary = tag;
else if(version < 2.2 && j == 2) partition = tag;
- else if(version >= 2.2 && j == 2) numPartitions = tag;
+ else if(version >= 2.2 && j == 2 && numTags > 3) numPartitions = tag;
else if(version >= 2.2 && j == 3) partition = tag;
- else if(j >= 4 && j < 4 + numPartitions) ghosts.push_back(-tag);
+ else if(j >= 4 && j < 4 + numPartitions - 1) ghosts.push_back(-tag);
else if(j == numTags - 1) parent = tag;
}
if(!(numVertices = MElement::getInfoMSH(type))) {
@@ -407,7 +407,12 @@ int GModel::readMSH(const std::string &name)
else Msg::Error("Parent element %d not found", parent);
}
else p = parents.find(parent)->second;
- e->setParent(p);
+ if(parentsOwned.find(p) == parentsOwned.end()) {
+ e->setParent(p, true);
+ parentsOwned.insert(p);
+ }
+ else
+ e->setParent(p, false);
}
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, "Reading elements");
@@ -432,9 +437,13 @@ int GModel::readMSH(const std::string &name)
int num = data[0];
int physical = (numTags > 0) ? data[1] : 0;
int elementary = (numTags > 1) ? data[2] : 0;
- int numPartitions = (version >= 2.2 && numTags > 2) ? data[3] : 0;
+ int numPartitions = (version >= 2.2 && numTags > 3) ? data[3] : 0;
int partition = (version < 2.2 && numTags > 2) ? data[3] :
(version >= 2.2 && numTags > 3) ? data[4] : 0;
+ int parent = (version < 2.2 && numTags > 3) ||
+ (version >= 2.2 && numPartitions && numTags > 3 + numPartitions) ||
+ (version >= 2.2 && !numPartitions && numTags > 2) ?
+ data[numTags] : 0;
int *indices = &data[numTags + 1];
std::vector<MVertex*> vertices;
if(vertexVector.size()){
@@ -448,6 +457,21 @@ int GModel::readMSH(const std::string &name)
if(numPartitions > 1)
for(int j = 0; j < numPartitions - 1; j++)
_ghostCells.insert(std::pair<MElement*, short>(e, -data[5 + j]));
+ if(parent) {
+ MElement *p = 0;
+ if(parents.find(parent) == parents.end()){
+ p = getParent(parent, e->getDim(), elements);
+ if(p) parents[parent] = p;
+ else Msg::Error("Parent element %d not found", parent);
+ }
+ else p = parents.find(parent)->second;
+ if(parentsOwned.find(p) == parentsOwned.end()) {
+ e->setParent(p, true);
+ parentsOwned.insert(p);
+ }
+ else
+ e->setParent(p, false);
+ }
if(numElements > 100000)
Msg::ProgressMeter(numElementsPartial + i + 1, numElements,
"Reading elements");
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 082a3841081be467c6fcb43c6677b341dd36045f..026b4b419ce089f07b4a087b6a00e3e7feb62302 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -403,23 +403,29 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
// if necessary, change the ordering of the vertices to get positive
// volume
setVolumePositive();
- int n = getNumVertices();
+ int n = getNumVerticesForMSH();
+ int par = (parentNum) ? 1 : 0;
+ bool poly = (type == MSH_POLYG_ || type == MSH_POLYH_);
if(!binary){
fprintf(fp, "%d %d", num ? num : _num, type);
if(version < 2.0)
fprintf(fp, " %d %d %d", abs(physical), elementary, n);
else if(!_partition)
- fprintf(fp, " 2 %d %d", abs(physical), elementary);
+ fprintf(fp, " %d %d %d", 2 + par, abs(physical), elementary);
else if(!ghosts)
- fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, _partition);
+ fprintf(fp, " %d %d %d 1 %d", 4 + par, abs(physical), elementary, _partition);
else{
int numGhosts = ghosts->size();
- fprintf(fp, " %d %d %d %d %d", 4 + numGhosts, abs(physical), elementary,
+ fprintf(fp, " %d %d %d %d %d", 4 + numGhosts + par, abs(physical), elementary,
1 + numGhosts, _partition);
for(unsigned int i = 0; i < ghosts->size(); i++)
fprintf(fp, " %d", -(*ghosts)[i]);
}
+ if(version >= 2.0 && parentNum)
+ fprintf(fp, " %d", parentNum);
+ if(version >= 2.0 && poly)
+ fprintf(fp, " %d", n);
}
else{
int numTags, numGhosts = 0;
@@ -429,6 +435,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
numGhosts = ghosts->size();
numTags = 4 + numGhosts;
}
+ numTags += par;
// we write elements in blobs of single elements; this will lead
// to suboptimal reads, but it's much simpler when the number of
// tags change from element to element (third-party codes can
@@ -438,14 +445,14 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
1 + numGhosts, _partition};
if(ghosts)
for(int i = 0; i < numGhosts; i++) blob[8 + i] = -(*ghosts)[i];
+ if(par) blob[8 + numGhosts] = par;
+ if(poly) Msg::Error("Unable to write polygons/polyhedra in binary files.");
fwrite(blob, sizeof(int), 4 + numTags, fp);
}
if(physical < 0) revert();
- int verts[60];
- for(int i = 0; i < n; i++)
- verts[i] = getVertex(i)->getIndex();
+ int *verts = getVerticesIdForMSH();
if(!binary){
for(int i = 0; i < n; i++)
@@ -457,6 +464,8 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
}
if(physical < 0) revert();
+
+ delete [] verts;
}
void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
@@ -764,6 +773,15 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
}
}
+int *MElement::getVerticesIdForMSH()
+{
+ int n = getNumVerticesForMSH();
+ int *verts = new int[n];
+ for(int i = 0; i < n; i++)
+ verts[i] = getVertex(i)->getIndex();
+ return verts;
+}
+
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 9292dee6071a44253ee83bb445d426ecf163a2fb..e026ace7c606ef147bfebcd7ce8dfc9305e244a3 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -152,7 +152,7 @@ class MElement
// get parent and children for hierarchial grids
virtual MElement *getParent() const { return NULL; }
- virtual void setParent(MElement *p) {}
+ virtual void setParent(MElement *p, bool owner = false) {}
virtual int getNumChildren() const { return 0; }
virtual MElement *getChild(int i) const { return NULL; }
virtual bool ownsParent() const { return false; }
@@ -274,6 +274,8 @@ class MElement
// return the number of vertices, as well as the element name if
// 'name' != 0
static int getInfoMSH(const int typeMSH, const char **const name=0);
+ virtual int getNumVerticesForMSH() { return getNumVertices(); }
+ virtual int *getVerticesIdForMSH();
static void registerBindings(binding *b);
};
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 6fa792afbb74ab2002d2604f53d2c4bdb82bc8da..e26e93cf5d3cf1cdfe5e6d902ba0af61533708db 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -136,57 +136,6 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = _intpt;
}
-void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical, int parentNum,
- std::vector<short> *ghosts)
-{
- int type = getTypeForMSH();
-
- if(!type) return;
-
- // if necessary, change the ordering of the vertices to get positive
- // volume
- setVolumePositive();
- int n = _parts.size() * 4;
- int numE = getNum();
- int partE = getPartition();
-
- if(!binary){
- fprintf(fp, "%d %d", num ? num : numE, type);
- if(version < 2.0)
- fprintf(fp, " %d %d", abs(physical), elementary);
- else {
- if(parentNum)
- fprintf(fp, " 5 %d %d 1 %d %d", abs(physical), elementary, partE, parentNum);
- else
- fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, partE);
- }
- }
- else{
- Msg::Error("Binary output not coded for polyhedra");
- }
-
- if(physical < 0) revert();
-
- fprintf(fp, " %d", n);
- int *verts = new int[n];
- for(unsigned int i = 0; i < _parts.size(); i++)
- for(int j = 0; j < 4; j++)
- verts[i * 4 + j] = _parts[i]->getVertex(j)->getIndex();
-
- if(!binary){
- for(int i = 0; i < n; i++)
- fprintf(fp, " %d", verts[i]);
- fprintf(fp, "\n");
- }
- else{
- fwrite(verts, sizeof(int), n, fp);
- }
-
- if(physical < 0) revert();
- delete [] verts;
-}
-
//------------------------------------------- MPolygon ------------------------------
void MPolygon::_initVertices()
@@ -274,57 +223,6 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = _intpt;
}
-void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical, int parentNum,
- std::vector<short> *ghosts)
-{
- int type = getTypeForMSH();
-
- if(!type) return;
-
- // if necessary, change the ordering of the vertices to get positive
- // volume
- setVolumePositive();
- int n = _parts.size() * 3;
- int numE = getNum();
- int partE = getPartition();
-
- if(!binary){
- fprintf(fp, "%d %d", num ? num : numE, type);
- if(version < 2.0)
- fprintf(fp, " %d %d", abs(physical), elementary);
- else {
- if(parentNum)
- fprintf(fp, " 5 %d %d 1 %d %d", abs(physical), elementary, partE, parentNum);
- else
- fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, partE);
- }
- }
- else{
- Msg::Error("Binary output not coded for polygons");
- }
-
- if(physical < 0) revert();
-
- fprintf(fp, " %d", n);
- int *verts = new int[n];
- for(unsigned int i = 0; i < _parts.size(); i++)
- for(int j = 0; j < 3; j++)
- verts[i * 3 + j] = _parts[i]->getVertex(j)->getIndex();
-
- if(!binary){
- for(int i = 0; i < n; i++)
- fprintf(fp, " %d", verts[i]);
- fprintf(fp, "\n");
- }
- else{
- fwrite(verts, sizeof(int), n, fp);
- }
-
- if(physical < 0) revert();
- delete [] verts;
-}
-
//----------------------------------- MTriangleBorder ------------------------------
void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index dc4da84e41c0ef9edcb867764aa5298f6dff67ed..438ab40c99c41b549b643189319162ba689c8efd 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -97,6 +97,16 @@ class MPolyhedron : public MElement {
}
virtual int getType() const { return TYPE_POLYH; }
virtual int getTypeForMSH() const { return MSH_POLYH_; }
+ virtual void revert()
+ {
+ for(int i = 0; i < _parts.size(); i++)
+ _parts[i]->revert();
+ _vertices.clear();
+ _innerVertices.clear();
+ _edges.clear();
+ _faces.clear();
+ _init();
+ }
virtual double getVolume()
{
double vol = 0;
@@ -123,14 +133,21 @@ class MPolyhedron : public MElement {
}
virtual bool isInside(double u, double v, double w);
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
- virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
- int num=0, int elementary=1, int physical=1,
- int parentNum=0, std::vector<short> *ghosts=0);
virtual MElement *getParent() const { return _orig; }
- virtual void setParent(MElement *p) { _orig = p; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
virtual int getNumChildren() const { return _parts.size(); }
virtual MElement *getChild(int i) const { return _parts[i]; }
virtual bool ownsParent() const { return _owner; }
+ virtual int getNumVerticesForMSH() {return _parts.size() * 4;}
+ virtual int *getVerticesIdForMSH()
+ {
+ int n = getNumVerticesForMSH();
+ int *verts = new int [n];
+ for(unsigned int i = 0; i < _parts.size(); i++)
+ for(int j = 0; j < 4; j++)
+ verts[i * 4 + j] = _parts[i]->getVertex(j)->getIndex();
+ return verts;
+ }
};
class MPolygon : public MElement {
@@ -208,6 +225,15 @@ class MPolygon : public MElement {
}
virtual int getType() const { return TYPE_POLYG; }
virtual int getTypeForMSH() const { return MSH_POLYG_; }
+ virtual void revert()
+ {
+ for(int i = 0; i < _parts.size(); i++)
+ _parts[i]->revert();
+ _vertices.clear();
+ _innerVertices.clear();
+ _edges.clear();
+ _initVertices();
+ }
virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
return _orig->getFunctionSpace(order);
@@ -226,14 +252,21 @@ class MPolygon : public MElement {
}
virtual bool isInside(double u, double v, double w);
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
- virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
- int num=0, int elementary=1, int physical=1,
- int parentNum=0, std::vector<short> *ghosts=0);
virtual MElement *getParent() const { return _orig; }
- virtual void setParent(MElement *p) { _orig = p; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
virtual int getNumChildren() const { return _parts.size(); }
virtual MElement *getChild(int i) const { return _parts[i]; }
virtual bool ownsParent() const { return _owner; }
+ virtual int getNumVerticesForMSH() {return _parts.size() * 3;}
+ virtual int *getVerticesIdForMSH()
+ {
+ int n = getNumVerticesForMSH();
+ int *verts = new int[n];
+ for(unsigned int i = 0; i < _parts.size(); i++)
+ for(int j = 0; j < 3; j++)
+ verts[i * 3 + j] = _parts[i]->getVertex(j)->getIndex();
+ return verts;
+ }
};
class MTriangleBorder : public MTriangle {