diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 6db55cc2bcac00bc2dc9b7f37b43e99120b6cf33..5b6ba9dec75b96f5111f736011b2fd6b222176f5 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -235,7 +235,7 @@ int MergeFile(std::string fileName, bool warnIfMissing)
}
char header[256];
- fgets(header, sizeof(header), fp);
+ if(!fgets(header, sizeof(header), fp)) return 0;
fclose(fp);
Msg::StatusBar(2, true, "Reading '%s'", fileName.c_str());
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 8b9b89c58d0549322bc11cba49a1ac7209bb18d3..b19caf2898cd4d417288dc2a1610e58f18c9278c 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -80,17 +80,17 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
return true;
}
-static void createElementMSH(GModel *m, int num, int typeMSH, int physical,
- int reg, int part, std::vector<MVertex*> &v,
- std::map<int, std::vector<MElement*> > elements[10],
- std::map<int, std::map<int, std::string> > physicals[4])
+static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
+ int reg, int part, std::vector<MVertex*> &v,
+ std::map<int, std::vector<MElement*> > elements[10],
+ std::map<int, std::map<int, std::string> > physicals[4])
{
MElementFactory factory;
MElement *e = factory.create(typeMSH, v, num, part);
if(!e){
Msg::Error("Unknown type of element %d", typeMSH);
- return;
+ return NULL;
}
switch(e->getType()){
@@ -114,7 +114,7 @@ static void createElementMSH(GModel *m, int num, int typeMSH, int physical,
elements[8][reg].push_back(e); break;
case TYPE_POLYH :
elements[9][reg].push_back(e); break;
- default : Msg::Error("Wrong type of element"); return;
+ default : Msg::Error("Wrong type of element"); return NULL;
}
int dim = e->getDim();
@@ -122,8 +122,48 @@ static void createElementMSH(GModel *m, int num, int typeMSH, int physical,
physicals[dim][reg][physical] = "unnamed";
if(part) m->getMeshPartitions().insert(part);
+ return e;
}
+static MElement *getParent(int parentNum, std::map<int, std::vector<MElement*> > &elements)
+{
+ std::map<int, std::vector<MElement*> >::iterator it = elements.begin();
+ for(; it != elements.end(); ++it) {
+ std::vector<MElement*>::iterator itE = it->second.begin();
+ for(; itE != it->second.end(); itE++)
+ if((*itE)->getNum() == parentNum) {
+ MElement *p = (*itE);
+ it->second.erase(itE);
+ return p;
+ }
+ }
+ return NULL;
+}
+
+static MElement *getParent(int parentNum, int dim,
+ std::map<int, std::vector<MElement*> > elements[10])
+{
+ MElement *parent = NULL;
+ switch(dim){
+ case 0 :
+ return getParent(parentNum, elements[0]);
+ case 1 :
+ return getParent(parentNum, elements[1]);
+ case 2 :
+ if(parent = getParent(parentNum, elements[2])) return parent;
+ if(parent = getParent(parentNum, elements[3])) return parent;
+ if(parent = getParent(parentNum, elements[8])) return parent;
+ return parent;
+ case 3 :
+ if(parent = getParent(parentNum, elements[4])) return parent;
+ if(parent = getParent(parentNum, elements[5])) return parent;
+ if(parent = getParent(parentNum, elements[6])) return parent;
+ if(parent = getParent(parentNum, elements[7])) return parent;
+ if(parent = getParent(parentNum, elements[9])) return parent;
+ return parent;
+ default : return parent;
+ }
+}
int GModel::readMSH(const std::string &name)
{
@@ -293,34 +333,38 @@ int GModel::readMSH(const std::string &name)
if(!fgets(str, sizeof(str), fp)) return 0;
int numElements;
+ std::map<int, MElement*> parents;
sscanf(str, "%d", &numElements);
Msg::Info("%d elements", numElements);
Msg::ResetProgressMeter();
if(!binary){
for(int i = 0; i < numElements; i++) {
- int num, type, physical = 0, elementary = 0, partition = 0, numVertices;
+ int num, type, physical = 0, elementary = 0, partition = 0, parentN = 0, numVertices;
if(version <= 1.0){
- fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices);
+ if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices) != 5)
+ return 0;
if(numVertices != MElement::getInfoMSH(type)) return 0;
}
else{
int numTags;
- fscanf(fp, "%d %d %d", &num, &type, &numTags);
+ if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
for(int j = 0; j < numTags; j++){
int tag;
- fscanf(fp, "%d", &tag);
+ if(fscanf(fp, "%d", &tag) != 1) return 0;
if(j == 0) physical = tag;
else if(j == 1) elementary = tag;
else if(j == 2) partition = tag;
+ else if(j == 3) parentN = tag;
// ignore any other tags for now
}
if(!(numVertices = MElement::getInfoMSH(type))) {
if(type != MSH_POLYG_ && type != MSH_POLYH_) return 0;
- fscanf(fp, "%d", &numVertices);
+ if(fscanf(fp, "%d", &numVertices) != 1) return 0;
}
}
int *indices = new int[numVertices];
- for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]);
+ for(int j = 0; j < numVertices; j++)
+ if(fscanf(fp, "%d", &indices[j]) != 1) return 0;
std::vector<MVertex*> vertices;
if(vertexVector.size()){
if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0;
@@ -328,8 +372,18 @@ int GModel::readMSH(const std::string &name)
else{
if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
}
- createElementMSH(this, num, type, physical, elementary, partition,
- vertices, elements, physicals);
+ MElement *e = createElementMSH(this, num, type, physical, elementary,
+ partition, vertices, elements, physicals);
+ if(parentN) {
+ MElement *parent = NULL;
+ if(parents.find(parentN) == parents.end()){
+ parent = getParent(parentN, e->getDim(), elements);
+ if(parent) parents[parentN] = parent;
+ else Msg::Error("Parent element %d not found!",parentN);
+ }
+ else parent = parents.find(parentN)->second;
+ e->setParent(parent);
+ }
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, "Reading elements");
delete [] indices;
@@ -372,7 +426,6 @@ int GModel::readMSH(const std::string &name)
numElementsPartial += numElms;
}
}
-
}
else if(!strncmp(&str[1], "NodeData", 8)) {
@@ -384,7 +437,6 @@ int GModel::readMSH(const std::string &name)
_vertexMapCache = vertexMap;
postpro = true;
break;
-
}
else if(!strncmp(&str[1], "ElementData", 11) ||
!strncmp(&str[1], "ElementNodeData", 15)) {
@@ -392,7 +444,6 @@ int GModel::readMSH(const std::string &name)
// there's some element post-processing data to read later on
postpro = true;
break;
-
}
do {
@@ -467,19 +518,42 @@ static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elem
}
template<class T>
-static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, bool saveAll,
+static void writeElementsMSH(FILE *fp, T *ele, bool saveAll,
double version, bool binary, int &num, int elementary,
+ std::vector<int> &physicals, int parentNum)
+{
+ if(saveAll)
+ ele->writeMSH(fp, version, binary, ++num, elementary, 0, parentNum);
+ else
+ for(unsigned int j = 0; j < physicals.size(); j++)
+ ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j], parentNum);
+}
+
+template<class T>
+static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, bool saveAll,
+ double version, bool binary, int &num, int elementary,
std::vector<int> &physicals)
{
for(unsigned int i = 0; i < ele.size(); i++)
- if(saveAll)
- ele[i]->writeMSH(fp, version, binary, ++num, elementary, 0);
- else
- for(unsigned int j = 0; j < physicals.size(); j++)
- ele[i]->writeMSH(fp, version, binary, ++num, elementary, physicals[j]);
+ writeElementsMSH(fp, ele[i], saveAll, version, binary, num,
+ elementary, physicals, 0);
+}
+
+template<class T>
+static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, bool saveAll,
+ double version, bool binary, int &num, int elementary,
+ std::vector<int> &physicals, std::map<MElement *, int> &parentsNum)
+{
+ for(unsigned int i = 0; i < ele.size(); i++) {
+ int parentNum = (ele[i]->getParent() != NULL &&
+ parentsNum.find(ele[i]->getParent()) != parentsNum.end()) ?
+ parentsNum.find(ele[i]->getParent())->second : 0;
+ writeElementsMSH(fp, ele[i], saveAll, version, binary, num,
+ elementary, physicals, parentNum);
+ }
}
-int GModel::writeMSH(const std::string &name, double version, bool binary,
+int GModel::writeMSH(const std::string &name, double version, bool binary,
bool saveAll, bool saveParametric, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
@@ -491,8 +565,6 @@ 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 index the vertices in a continuous
// sequence
int numVertices = indexMeshVertices(saveAll);
@@ -504,6 +576,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
// list have the same type, i.e., they are all of the same
// polynomial order)
std::map<int, int> elements;
+ std::map<MElement *, GEntity *> parents[2];
+ int numParents = 0;
for(viter it = firstVertex(); it != lastVertex(); ++it){
int p = (saveAll ? 1 : (*it)->physicals.size());
int n = p * (*it)->points.size();
@@ -522,6 +596,11 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
if(n) elements[(*it)->quadrangles[0]->getTypeForMSH()] += n;
n = p * (*it)->polygons.size();
if(n) elements[(*it)->polygons[0]->getTypeForMSH()] += n;
+ for(int i = 0; i < (*it)->polygons.size(); i++)
+ if((*it)->polygons[i]->ownsParent()) {
+ parents[0][(*it)->polygons[i]->getParent()] = (*it);
+ numParents += p;
+ }
}
for(riter it = firstRegion(); it != lastRegion(); ++it){
int p = (saveAll ? 1 : (*it)->physicals.size());
@@ -535,14 +614,18 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
if(n) elements[(*it)->pyramids[0]->getTypeForMSH()] += n;
n = p * (*it)->polyhedra.size();
if(n) elements[(*it)->polyhedra[0]->getTypeForMSH()] += n;
+ for(int i = 0; i < (*it)->polyhedra.size(); i++)
+ if((*it)->polyhedra[i]->ownsParent()) {
+ parents[1][(*it)->polyhedra[i]->getParent()] = (*it);
+ numParents += p;
+ }
}
-
-
int numElements = 0;
std::map<int, int>::const_iterator it = elements.begin();
for(; it != elements.end(); ++it)
numElements += it->second;
+ numElements += numParents;
if(version >= 2.0){
fprintf(fp, "$MeshFormat\n");
@@ -567,7 +650,6 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
fprintf(fp, "$ParametricNodes\n");
else
fprintf(fp, "$Nodes\n");
-
}
else
fprintf(fp, "$NOD\n");
@@ -580,8 +662,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric,
scalingFactor);
-
-
+
if(binary) fprintf(fp, "\n");
if(version >= 2.0){
@@ -599,50 +680,102 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
fprintf(fp, "%d\n", numElements);
int num = 0;
+ std::map<MElement *, int> parentsNum;
+ std::map<MElement *, GEntity *>::const_iterator itP;
+
writeElementHeaderMSH(binary, fp, elements, MSH_PNT);
for(viter it = firstVertex(); it != lastVertex(); ++it)
writeElementsMSH(fp, (*it)->points, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
- writeElementHeaderMSH
- (binary, fp, elements, MSH_LIN_2, MSH_LIN_3, MSH_LIN_4, MSH_LIN_5);
+
+ writeElementHeaderMSH(binary, fp, elements, MSH_LIN_2, MSH_LIN_3, MSH_LIN_4, MSH_LIN_5);
for(eiter it = firstEdge(); it != lastEdge(); ++it)
writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
- writeElementHeaderMSH(binary, fp, elements, MSH_TRI_3, MSH_TRI_6, MSH_TRI_9,
+
+ writeElementHeaderMSH(binary, fp, elements, MSH_TRI_3, MSH_TRI_6, MSH_TRI_9,
MSH_TRI_10, MSH_TRI_12, MSH_TRI_15, MSH_TRI_15I, MSH_TRI_21);
+ for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
+ if(itP->first->getType() == TYPE_TRI) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
-
writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8);
+ for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
+ if(itP->first->getType() == TYPE_QUA) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
writeElementHeaderMSH(binary, fp, elements, MSH_POLYG_);
+ for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
+ if(itP->first->getType() == TYPE_POLYG) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(fiter it = firstFace(); it != lastFace(); it++)
writeElementsMSH(fp, (*it)->polygons, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, parentsNum);
+
writeElementHeaderMSH(binary, fp, elements, MSH_TET_4, MSH_TET_10, MSH_TET_20,
MSH_TET_35, MSH_TET_56, MSH_TET_52);
+ for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
+ if(itP->first->getType() == TYPE_TET) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
writeElementHeaderMSH(binary, fp, elements, MSH_HEX_8, MSH_HEX_27, MSH_HEX_20);
+ for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
+ if(itP->first->getType() == TYPE_HEX) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
writeElementHeaderMSH(binary, fp, elements, MSH_PRI_6, MSH_PRI_18, MSH_PRI_15);
+ for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
+ if(itP->first->getType() == TYPE_PRI) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
writeElementHeaderMSH(binary, fp, elements, MSH_PYR_5, MSH_PYR_14, MSH_PYR_13);
+ for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
+ if(itP->first->getType() == TYPE_PYR) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
writeElementHeaderMSH(binary, fp, elements, MSH_POLYH_);
+ for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
+ if(itP->first->getType() == TYPE_POLYH) {
+ writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
+ itP->second->tag(), itP->second->physicals, 0);
+ parentsNum[itP->first] = num;
+ }
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->polyhedra, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, parentsNum);
if(binary) fprintf(fp, "\n");
@@ -732,7 +865,7 @@ int GModel::readSTL(const std::string &name, double tolerance)
// "solid", or binary data header
char buffer[256];
- fgets(buffer, sizeof(buffer), fp);
+ if(!fgets(buffer, sizeof(buffer), fp)) return 0;
// workaround for stupid tools which use "solid" to start their
// binary files
@@ -946,7 +1079,7 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
const char *format;
fpos_t position;
fgetpos(fp, &position);
- fgets(tmp, sizeof(tmp), fp);
+ if(!fgets(tmp, sizeof(tmp), fp)) return 0;;
fsetpos(fp, &position);
if(strstr(tmp, ","))
format = " , %d";
@@ -1366,7 +1499,7 @@ int GModel::readMESH(const std::string &name)
}
char buffer[256];
- fgets(buffer, sizeof(buffer), fp);
+ if(!fgets(buffer, sizeof(buffer), fp)) return 0;
char str[256];
int format;
@@ -2153,10 +2286,11 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
char buffer[256], buffer2[256];
- fgets(buffer, sizeof(buffer), fp); // version line
- fgets(buffer, sizeof(buffer), fp); // title
+ if(!fgets(buffer, sizeof(buffer), fp)) return 0; // version line
+ if(!fgets(buffer, sizeof(buffer), fp)) return 0; // title
- fscanf(fp, "%s", buffer); // ASCII or BINARY
+ if(fscanf(fp, "%s", buffer) != 1) // ASCII or BINARY
+ Msg::Error("Failed reading buffer");
bool binary = false;
if(!strcmp(buffer, "BINARY")) binary = true;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 80e8823a6c15dff2e829fc3b07698c699bea13eb..69ee7515c946b8191fd4a146ee5b38a3f4654dcd 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -382,7 +382,7 @@ double MElement::interpolateDiv(double val[], double u, double v, double w, int
}
void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical)
+ int elementary, int physical, int parentNum)
{
int type = getTypeForMSH();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index c259731866bae9ac59da2ac560504b0b8b43d07e..d7daec74470df49a46c9c14aba64ac7df0cb02d5 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -132,9 +132,11 @@ class MElement
}
// get parent and children for hierarchial grids
- virtual MElement *getParent() const { return 0; }
+ virtual MElement *getParent() const { return NULL; }
+ virtual void setParent(MElement *p) {}
virtual int getNumChildren() const { return 0; }
- virtual MElement *getChild(int i) const { return 0; }
+ virtual MElement *getChild(int i) const { return NULL; }
+ virtual bool ownsParent() const { return false; }
//get the type of the element
virtual int getType() const = 0;
@@ -216,7 +218,7 @@ class MElement
// IO routines
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
- int num=0, int elementary=1, int physical=1);
+ int num=0, int elementary=1, int physical=1, int parentNum=0);
virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
bool printGamma, bool printEta, bool printRho,
bool printDisto,double scalingFactor=1.0, int elementary=1);
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 068374cc223609d14ab64cb03a474db913268658..2bd688d6fac21c34fcb676e4c74b6a2dadabcb5c 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -134,7 +134,7 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
}
void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical)
+ int elementary, int physical, int parentNum)
{
int type = getTypeForMSH();
@@ -151,8 +151,12 @@ void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
fprintf(fp, "%d %d", num ? num : numE, type);
if(version < 2.0)
fprintf(fp, " %d %d", abs(physical), elementary);
- else
- fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+ else {
+ if(parentNum)
+ fprintf(fp, " 4 %d %d %d %d", abs(physical), elementary, partE, parentNum);
+ else
+ fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+ }
}
else{
int tags[4] = {num ? num : numE, abs(physical), elementary, partE};
@@ -306,7 +310,7 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
}
void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical)
+ int elementary, int physical, int parentNum)
{
int type = getTypeForMSH();
@@ -323,8 +327,12 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
fprintf(fp, "%d %d", num ? num : numE, type);
if(version < 2.0)
fprintf(fp, " %d %d", abs(physical), elementary);
- else
- fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+ else {
+ if(parentNum)
+ fprintf(fp, " 4 %d %d %d %d", abs(physical), elementary, partE, parentNum);
+ else
+ fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+ }
}
else{
int tags[4] = {num ? num : numE, abs(physical), elementary, partE};
@@ -495,7 +503,6 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
unsigned int nbQ = quads.size();
unsigned int nbTe = tetras.size();
unsigned int nbH = hexas.size();
- std::map<int, double> nodeLs[8];
// copy element
std::vector<MVertex*> vmv;
@@ -528,7 +535,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
MElementFactory factory;
- MElement *copy = factory.create(eType, vmv, e->getNum(), ePart);
+ MElement *copy = factory.create(eType, vmv, ++numEle, ePart);
switch (eType) {
case MSH_TET_4 :
@@ -542,6 +549,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
+ std::map<int, double> nodeLs[4];
for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
@@ -555,6 +563,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z(),
e->getVertex(6)->x(), e->getVertex(6)->y(), e->getVertex(6)->z(),
e->getVertex(7)->x(), e->getVertex(7)->y(), e->getVertex(7)->z());
+ std::map<int, double> nodeLs[8];
for(int i = 0; i < 8; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
isCut = H.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
hexas, tetras, quads, triangles, lines, 0, nodeLs);
@@ -564,6 +573,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
+ std::map<int, double> nodeLs[4];
for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
nodeLs[3] = verticesLs[e->getVertex(5)->getNum()];
bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
@@ -593,6 +603,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
+ std::map<int, double> nodeLs[4];
for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
nodeLs[3] = verticesLs[e->getVertex(4)->getNum()];
bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
@@ -608,13 +619,14 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
isCut = iC1 || iC2;
}
else if(eType == MSH_POLYH_){
+ std::map<int, double> nodeLs[4];
for(int i = 0; i < e->getNumChildren(); i++) {
MTetrahedron *t = (MTetrahedron*) e->getChild(i);
DI_Tetra Tet(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
t->getVertex(3)->x(), t->getVertex(3)->y(), t->getVertex(3)->z());
- for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
+ for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[t->getVertex(i)->getNum()];
bool iC = Tet.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
isCut = (isCut || iC);
@@ -657,13 +669,13 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
poly[1].push_back(mt);
}
if(poly[0].size()) {
- p1 = new MPolyhedron(poly[0], copy, true, ++numEle, ePart);
+ p1 = new MPolyhedron(poly[0], ++numEle, ePart, true, copy);
int reg = getElementaryTag(-1, elementary, newtags[3]);
elements[9][reg].push_back(p1);
assignPhysicals(GM, gePhysicals, reg, 3, physicals);
}
if(poly[1].size()) {
- p2 = new MPolyhedron(poly[1], copy, false, ++numEle, ePart);
+ p2 = new MPolyhedron(poly[1], ++numEle, ePart, false, copy);
elements[9][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 3, physicals);
}
@@ -728,6 +740,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
DI_Triangle T(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
+ std::map<int, double> nodeLs[3];
for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, 0, nodeLs);
@@ -737,17 +750,19 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
+ std::map<int, double> nodeLs[4];
for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
isCut = Q.cut(*ls, ipV, ipS, cp, integOrder,integOrder,integOrder,
quads, triangles, lines, 0, nodeLs);
}
else if(eType == MSH_POLYG_){
+ std::map<int, double> nodeLs[3];
for(int i = 0; i < e->getNumChildren(); i++) {
MTriangle *t = (MTriangle*) e->getChild(i);
DI_Triangle Tri(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());
- for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
+ for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[t->getVertex(i)->getNum()];
bool iC = Tri.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, 0, nodeLs);
isCut = (isCut || iC);
@@ -790,13 +805,13 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
poly[1].push_back(mt);
}
if(poly[0].size()) {
- p1 = new MPolygon(poly[0], copy, true, ++numEle, ePart);
+ p1 = new MPolygon(poly[0], ++numEle, ePart, true, copy);
int reg = getElementaryTag(-1, elementary, newtags[2]);
elements[8][reg].push_back(p1);
assignPhysicals(GM, gePhysicals, reg, 2, physicals);
}
if(poly[1].size()) {
- p2 = new MPolygon(poly[1], copy, false, ++numEle, ePart);
+ p2 = new MPolygon(poly[1], ++numEle, ePart, false, copy);
elements[8][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
}
@@ -852,6 +867,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
{
DI_Line L(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z());
+ std::map<int, double> nodeLs[2];
for(int i = 0; i < 2; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
isCut = L.cut(*ls, ipV, cp, integOrder, lines, 0, nodeLs);
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index c52f486ab741108c2c2fbd03065b9dd31cdb1066..7277f0fb4c1672e6800277d1fa3d1adb9ceeda63 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -26,8 +26,9 @@ class MPolyhedron : public MElement {
std::vector<MFace> _faces;
void _init();
public:
- MPolyhedron(std::vector<MVertex*> v, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0)
+ MPolyhedron(std::vector<MVertex*> v, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
+ : MElement(num, part), _owner(owner), _orig(orig)
{
if(v.size() % 4){
Msg::Error("Got %d vertices for polyhedron", (int)v.size());
@@ -37,14 +38,8 @@ class MPolyhedron : public MElement {
_parts.push_back(new MTetrahedron(v[i], v[i + 1], v[i + 2], v[i + 3]));
_init();
}
- MPolyhedron(std::vector<MTetrahedron*> vT, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0)
- {
- for(unsigned int i = 0; i < vT.size(); i++)
- _parts.push_back(vT[i]);
- _init();
- }
- MPolyhedron(std::vector<MTetrahedron*> vT, MElement* orig, bool owner, int num=0, int part=0)
+ MPolyhedron(std::vector<MTetrahedron*> vT, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
: MElement(num, part), _owner(owner), _orig(orig)
{
for(unsigned int i = 0; i < vT.size(); i++)
@@ -131,10 +126,12 @@ class MPolyhedron : public MElement {
return _parts.size() * getNGQTetPts(pOrder);
}
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
- int num=0, int elementary=1, int physical=1);
+ int num=0, int elementary=1, int physical=1, int parentNum=0);
virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p) { _orig = p; }
virtual int getNumChildren() const { return _parts.size(); }
virtual MElement *getChild(int i) const { return _parts[i]; }
+ virtual bool ownsParent() const { return _owner; }
};
class MPolygon : public MElement {
@@ -147,15 +144,17 @@ class MPolygon : public MElement {
std::vector<MEdge> _edges;
void _initVertices();
public:
- MPolygon(std::vector<MVertex*> v, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0)
+ MPolygon(std::vector<MVertex*> v, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
+ : MElement(num, part), _owner(owner), _orig(orig)
{
for(unsigned int i = 0; i < v.size() / 3; i++)
_parts.push_back(new MTriangle(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]));
_initVertices();
}
- MPolygon(std::vector<MElement*> vT, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0)
+ MPolygon(std::vector<MTriangle*> vT, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
+ : MElement(num, part), _owner(owner), _orig(orig)
{
for(unsigned int i = 0; i < vT.size(); i++){
MTriangle *t = (MTriangle*) vT[i];
@@ -163,20 +162,6 @@ class MPolygon : public MElement {
}
_initVertices();
}
- MPolygon(std::vector<MTriangle*> vT, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0)
- {
- for(unsigned int i = 0; i < vT.size(); i++)
- _parts.push_back(vT[i]);
- _initVertices();
- }
- MPolygon(std::vector<MTriangle*> vT, MElement* orig, bool owner, int num=0, int part=0)
- : MElement(num, part), _owner(owner), _orig(orig)
- {
- for(unsigned int i = 0; i < vT.size(); i++)
- _parts.push_back(vT[i]);
- _initVertices();
- }
~MPolygon()
{
if(_owner)
@@ -245,10 +230,12 @@ class MPolygon : public MElement {
return _parts.size() * getNGQTPts(pOrder);
}
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
- int num=0, int elementary=1, int physical=1);
+ int num=0, int elementary=1, int physical=1, int parentNum=0);
virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p) { _orig = p; }
virtual int getNumChildren() const { return _parts.size(); }
virtual MElement *getChild(int i) const { return _parts[i]; }
+ virtual bool ownsParent() const { return _owner; }
};
class MTriangleBorder : public MTriangle {
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 5cc6391926fc28bd621b2b9e6911a7abd2da0d56..e7288d6fd7b419ab93a297df83527a349fe8803a 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -662,7 +662,7 @@ void DI_Element::addLs (const double *ls) {
for(int i = 0; i < nbMid(); i++)
mid_[i]->addLs(ls[nbVert()+i]);
}
-void DI_Element::addLs (int primTag, std::map<int, double> nodeLs[8]) { printf("p");
+void DI_Element::addLs (int primTag, std::map<int, double> nodeLs[8]) {
for(int i = 0; i < nbVert(); i++)
pts_[i]->addLs(nodeLs[i][primTag]);
}
@@ -673,7 +673,7 @@ void DI_Element::addLs (const DI_Element *e) {
for(int i = 0; i < nbMid(); i++)
mid_[i]->addLs(e);
}
-void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) { printf("e");
+void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
if(type() != e->type()) {
printf("Error : addLs with element of different type\n");
}