diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 0494ade2e52c493d60f2abdbbf24c2ecacc89b26..ccab6170755c50beefcaf1c0306abfb013dc92ea 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -935,6 +935,9 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "%d\n", count);
+ int type = -1;
+ int physical = 0;
+ int elementary = 0;
int partition = 0;
for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
@@ -942,7 +945,9 @@ int CellComplex::writeComplexMSH(const std::string &name){
if(vertex->inSubdomain()) partition = 3;
else if(vertex->onDomainBoundary()) partition = 2;
else partition = 1;
- fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getNum(), 15, 3, 0, 0, partition, vertex->getVertex(0)->getNum());
+ type = vertex->getTypeForMSH();
+ //vertex->writeMSH(fp, 2.1, false, 0, elementary, physical);
+ fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getNum(), type, 3, physical, elementary, partition, vertex->getVertex(0)->getNum());
}
std::list< std::pair<int, Cell*> > cells;
@@ -954,7 +959,9 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = edge->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+ type = cell->getTypeForMSH();
+ //cell->writeMSH(fp, 2.1, false, 0, elementary, physical);
+ fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
}
}
@@ -966,8 +973,9 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = face->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 3, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ type = cell->getTypeForMSH();
+ if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
}
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
@@ -978,50 +986,53 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = volume->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 12, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 13, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 14, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
+ type = cell->getTypeForMSH();
+ if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
+ if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
+ if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
}
}
for(int i = 0; i < _store.size(); i++){
for(citer cit = _store.at(i).begin(); cit != _store.at(i).end(); cit++){
Cell* cell = *cit;
+ type = cell->getTypeForMSH();
if(cell->inSubdomain()) partition = 3;
else if(cell->onDomainBoundary()) partition = 2;
else partition = 1;
- if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), 15, 3, 0, 0, partition, cell->getVertex(0)->getNum());
- if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+ if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum());
+ if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
if(cell->getDim() == 2){
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 3, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
}
if(cell->getDim() == 3){
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 12, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 13, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 14, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
+ if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
+ if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
+ if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
}
}
}
for(citer cit = _ecells.begin(); cit != _ecells.end(); cit++){
Cell* cell = *cit;
+ type = cell->getTypeForMSH();
if(cell->inSubdomain()) partition = 3;
else if(cell->onDomainBoundary()) partition = 2;
else partition = 1;
- if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), 15, 3, 0, 0, partition, cell->getVertex(0)->getNum());
- if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+ if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum());
+ if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
if(cell->getDim() == 2){
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 3, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
}
if(cell->getDim() == 3){
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 12, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 13, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 14, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
+ if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
+ if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
+ if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
}
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 79667eee33a96b5ccefa5ba02c78481b161f1a56..41a9bb0e85008b78759370b2b82a17a5793b1c7e 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -96,6 +96,7 @@ class Cell
virtual int getType() { return -1; }
virtual int getTypeForMSH() { return -1; }
virtual int getPartition() { return -1; }
+ virtual void setPartition(int num){}
virtual void setImmune(bool immune) { _immune = immune; };
virtual bool getImmune() const { return _immune; };
@@ -327,6 +328,7 @@ class Cell
// checks whether lower dimensional simplex tau is on the boundary of this cell
virtual int kappa(Cell* tau) const;
+ virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, int num=0, int elementary=1, int physical=1, int parentNum=0) {}
};
@@ -418,6 +420,7 @@ class OneSimplex : public Simplex, public MLine
int kappa(Cell* tau) const;
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
+
};
// Two simplex cell type.
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 5f78674b7773fcc301320f2c9ffb698bf8520d7d..1d86167fb6fb57a5e36ee1e3265f33dc4284726c 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -834,7 +834,7 @@ int Chain::writeChainMSH(const std::string &name){
//_cellComplex->writeComplexMSH(name);
- if(getSize() == 0) return 0;
+ if(getSize() == 0) return 1;
FILE *fp = fopen(name.c_str(), "a");
if(!fp){
@@ -855,26 +855,16 @@ int Chain::writeChainMSH(const std::string &name){
fprintf(fp, "%d \n", getSize());
fprintf(fp, "0 \n");
- //std::list< std::pair<int, Cell*> > cells;
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
int coeff = (*cit).second;
fprintf(fp, "%d %d \n", cell->getNum(), coeff );
- /*
- cells = cell->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell2 = (*it).second;
- fprintf(fp, "%d %d \n", cell2->getNum(), getCoeff(i)*(*it).first );
- }
- */
}
fprintf(fp, "$EndElementData\n");
-
fclose(fp);
return 1;
-
}
void Chain::createPView(){
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 06e85c8569077d5e08fa9cab280adbf8c176979c..59724dd1b2b09b85264e208103ea7212d21b0cae 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -266,12 +266,10 @@ class Chain{
// (for primary complex chains only, not for dual chains represented by primary cells (yet).)
void smoothenChain();
- // append this chain to a 2.1 MSH ASCII file as $ElementData
+ // append this chain to a MSH ASCII file as $ElementData
int writeChainMSH(const std::string &name);
// create a PView of this chain.
- // Warning: saving the PView in the GUI changes the numbering of the mesh,
- // but NOT the numbering of the post-processing ElementData accordingly! (See Post/PViewDataGModelIO.cpp writeMSH)
void createPView();
};
diff --git a/Geo/GModel.h b/Geo/GModel.h
index fc5f44e668045db87a548bfed469fcbe4da250d2..1ab199100315779e3de863b366430c76546f3c4c 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -371,7 +371,7 @@ class GModel
// Gmsh mesh file format
int readMSH(const std::string &name);
int writeMSH(const std::string &name, double version=1.0, bool binary=false,
- bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0);
+ bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0, bool eleRenumbering=true);
int writeDistanceMSH(const std::string &name, double version=1.0, bool binary=false,
bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 49baa39fa540526f7173563bed811e8ca7815f2b..edff2d5d6d62458d58fe9ca17d04bab180bfea3d 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -529,36 +529,41 @@ static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elem
template<class T>
static void writeElementsMSH(FILE *fp, T *ele, bool saveAll,
double version, bool binary, int &num, int elementary,
- std::vector<int> &physicals, int parentNum)
+ std::vector<int> &physicals, int parentNum, bool eleRenumbering=true)
{
- 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);
+ if(saveAll){
+ if(eleRenumbering) ele->writeMSH(fp, version, binary, ++num, elementary, 0, parentNum);
+ else ele->writeMSH(fp, version, binary, 0, elementary, 0, parentNum);
+ }
+ else{
+ for(unsigned int j = 0; j < physicals.size(); j++){
+ if(eleRenumbering) ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j], parentNum);
+ else ele->writeMSH(fp, version, binary, 0, elementary, physicals[j], parentNum);
+ }
+ }
}
template<class T>
static void writeElementsMSH(FILE *fp, std::vector<T*> &ele, bool saveAll,
double version, bool binary, int &num, int elementary,
- std::vector<int> &physicals)
+ std::vector<int> &physicals, bool eleRenumbering=true)
{
for(unsigned int i = 0; i < ele.size(); i++)
writeElementsMSH(fp, ele[i], saveAll, version, binary, num,
- elementary, physicals, 0);
+ elementary, physicals, 0, eleRenumbering);
}
template<class T>
static void writeElementsMSH(FILE *fp, std::vector<T*> &ele, bool saveAll,
double version, bool binary, int &num, int elementary,
- std::vector<int> &physicals, std::map<MElement *, int> &parentsNum)
+ std::vector<int> &physicals, std::map<MElement *, int> &parentsNum, bool eleRenumbering=true)
{
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);
+ elementary, physicals, parentNum, eleRenumbering);
}
}
int GModel::writeDistanceMSH(const std::string &name, double version, bool binary,
@@ -653,7 +658,7 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
}
int GModel::writeMSH(const std::string &name, double version, bool binary,
- bool saveAll, bool saveParametric, double scalingFactor)
+ bool saveAll, bool saveParametric, double scalingFactor, bool eleRenumbering)
{
FILE *fp = fopen(name.c_str(), binary ? "wb" : "w");
if(!fp){
@@ -789,99 +794,99 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
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);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8,
MSH_QUA_16, MSH_QUA_25, MSH_QUA_36, MSH_QUA_12, MSH_QUA_16I,
MSH_QUA_20);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(fiter it = firstFace(); it != lastFace(); it++)
writeElementsMSH(fp, (*it)->polygons, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals, parentsNum);
+ (*it)->tag(), (*it)->physicals, parentsNum, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals);
+ (*it)->tag(), (*it)->physicals, eleRenumbering);
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);
+ itP->second->tag(), itP->second->physicals, 0, eleRenumbering);
parentsNum[itP->first] = num;
}
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->polyhedra, saveAll, version, binary, num,
- (*it)->tag(), (*it)->physicals, parentsNum);
+ (*it)->tag(), (*it)->physicals, parentsNum, eleRenumbering);
if(binary) fprintf(fp, "\n");
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 36c986821104f3f81b6df2df419b7615615c7286..5b81afce1e155806ecd471715a27849b858ec2b7 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -74,8 +74,6 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
void Homology::findGenerators(std::string fileName){
-
-
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
//_cellComplex->printEuler();
@@ -103,8 +101,6 @@ void Homology::findGenerators(std::string fileName){
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
- //_cellComplex->writeComplexMSH(fileName);
-
Msg::Info("Computing homology groups...");
t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex);
@@ -112,8 +108,6 @@ void Homology::findGenerators(std::string fileName){
t2 = Cpu();
Msg::Info("Homology Computation complete (%g s).", t2 - t1);
- std::vector<Chain*> chainVector;
-
int HRank[4];
for(int j = 0; j < 4; j++){
HRank[j] = 0;
@@ -126,8 +120,6 @@ void Homology::findGenerators(std::string fileName){
std::string name = "H" + dimension + getDomainString() + generator;
Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, _model, name, chains->getTorsion(j,i));
- //Chain* chain2 = new Chain(chain);
- //printf("chain %d \n", i);
t1 = Cpu();
int start = chain->getSize();
chain->smoothenChain();
@@ -137,9 +129,7 @@ void Homology::findGenerators(std::string fileName){
HRank[j] = HRank[j] + 1;
if(chain->getTorsion() != 1) Msg::Warning("H%d %d has torsion coefficient %d!", j, i, chain->getTorsion());
}
- chainVector.push_back(chain);
- //chainVector.push_back(chain2);
- //delete chain;
+ _generators[j].push_back(chain);
}
if(j == _cellComplex->getDim() && _cellComplex->getNumOmitted() > 0){
for(int i = 0; i < _cellComplex->getNumOmitted(); i++){
@@ -149,25 +139,15 @@ void Homology::findGenerators(std::string fileName){
std::vector<int> coeffs (_cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, _model, name, 1);
if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
- //delete chain;
- chainVector.push_back(chain);
+ _generators[j].push_back(chain);
}
}
}
- _cellComplex->writeComplexMSH(fileName);
- //writeHeaderMSH(fileName);
- for(int i = 0; i < chainVector.size(); i++){
- Chain* chain = chainVector.at(i);
- chain->writeChainMSH(fileName);
- chain->createPView();
- chainVector.at(i) = NULL;
- delete chain;
- }
- chainVector.clear();
-
+ createPViews();
+ writeGeneratorsMSH(fileName);
Msg::Info("Ranks of homology groups for primal cell complex:");
Msg::Info("H0 = %d", HRank[0]);
@@ -176,9 +156,6 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("The computation of generators in %d highest dimensions was omitted.", omitted);
- Msg::Info("Wrote results to %s.", fileName.c_str());
- printf("Wrote results to %s. \n", fileName.c_str());
-
delete chains;
printf("H0 = %d \n", HRank[0]);
@@ -252,14 +229,11 @@ void Homology::findDualGenerators(std::string fileName){
std::string name = "H" + dimension + "*" + getDomainString() + generator;
Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, _model, name, chains->getTorsion(j,i));
- chain->writeChainMSH(fileName);
- chain->createPView();
+ _generators[dim-j].push_back(chain);
if(chain->getSize() != 0){
HRank[dim-j] = HRank[dim-j] + 1;
if(chain->getTorsion() != 1) Msg::Warning("H%d* %d has torsion coefficient %d!", dim-j, i, chain->getTorsion());
- }
- delete chain;
-
+ }
}
@@ -270,16 +244,17 @@ void Homology::findDualGenerators(std::string fileName){
std::string name = "H" + dimension + "*" + getDomainString() + generator;
std::vector<int> coeffs (_cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, _model, name, 1);
- chain->writeChainMSH(fileName);
+ _generators[dim-j].push_back(chain);
if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1;
- delete chain;
-
}
}
}
+ createPViews();
+ writeGeneratorsMSH(fileName);
+
Msg::Info("Ranks of homology groups for dual cell complex:");
Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]);
@@ -287,9 +262,6 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("H3* = %d", HRank[3]);
if(omitted != 0) Msg::Info("The computation of %d highest dimension dual generators was omitted.", omitted);
- Msg::Info("Wrote results to %s.", fileName.c_str());
- printf("Wrote results to %s. \n", fileName.c_str());
-
delete chains;
printf("H0* = %d \n", HRank[0]);
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 844af0b60836f711cd87b5859ecb5cf3e5dbc6d7..3dd6e2b7e6066f0fc4936eb93fcce7a7449b283b 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -10,6 +10,7 @@
#include <sstream>
#include "CellComplex.h"
+#include "ChainComplex.h"
#if defined(HAVE_KBIPACK)
@@ -25,16 +26,30 @@ bool convert(const TTypeA& input, TTypeB& output ){
class Homology
{
private:
+
+ // base cell complex and the model of the homology computation
CellComplex* _cellComplex;
GModel* _model;
+ // domain and the relative subdomain of the homology computation
std::vector<int> _domain;
std::vector<int> _subdomain;
+
+ // generator chains
+ std::vector<Chain*> _generators[4];
public:
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
- ~Homology(){ delete _cellComplex; }
+ ~Homology(){
+ delete _cellComplex;
+ for(int i = 0; i < 4; i++) {
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ delete chain;
+ }
+ }
+ }
// Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
void findGenerators(std::string fileName);
@@ -45,7 +60,10 @@ class Homology
//void swapSubdomain() { _cellComplex->swapSubdomain(); }
// Restore the cell complex to its original state before cell reductions
- void restoreHomology() { _cellComplex->restoreComplex(); }
+ void restoreHomology() {
+ _cellComplex->restoreComplex();
+ for(int i = 0; i < 4; i++) _generators[i].clear();
+ }
// Create a string describing the generator
std::string getDomainString() {
@@ -78,17 +96,31 @@ class Homology
return domainString;
}
- int writeHeaderMSH(const std::string &name){
- FILE *fp = fopen(name.c_str(), "w");
- if(!fp){
- Msg::Error("Unable to open file '%s'", name.c_str());
- printf("Unable to open file.");
- return 0;
+ // create PViews of the generators
+ void createPViews(){
+ for(int i = 0; i < 4; i++){
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ chain->createPView();
}
- fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
- fclose(fp);
- return 1;
+ }
}
+
+ // write the generators to a file
+ bool writeGeneratorsMSH(std::string fileName, bool binary=false){
+ if(!_model->writeMSH(fileName, 2.0, binary, true, false, 1.0, false)) return false;
+ for(int i = 0; i < 4; i++){
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ if(!chain->writeChainMSH(fileName)) return false;
+ }
+ }
+ Msg::Info("Wrote results to %s.", fileName.c_str());
+ printf("Wrote results to %s. \n", fileName.c_str());
+
+ return true;
+ }
+
};
#endif
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 60bfa2c41ec216e1a33e4d2aa19b63dbad8c59c0..0cad126200fab4f2c63636077d35377d09a7a6dc 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -117,7 +117,7 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
GModel *model = _steps[0]->getModel();
- if(!model->writeMSH(fileName, 2.0, binary, true)) return false;
+ if(!model->writeMSH(fileName, 2.0, binary, true, false, 1.0, false)) return false;
// append data
FILE *fp = fopen(fileName.c_str(), binary ? "ab" : "a");
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index f15f751a7da5ab99b06586407bfcc174473bb290..0335ebc0b39b46e5922317f39c16e71961106262 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -123,7 +123,6 @@ Physical Surface(74) = {46, 18, 20, 52, 22, 50, 24, 48, 66, 63, 60, 58, 56, 54};
// Complement of the domain surface respect to the four terminals
Physical Surface(75) = {46, 63, 66, 52, 50, 48, 54, 60, 58, 56};
-
// Create a mesh of the model
Mesh 3;