diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 742d85a49d1ec20193975a845618e149d9532ca5..79667eee33a96b5ccefa5ba02c78481b161f1a56 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -40,7 +40,7 @@ class Less_Cell{
bool operator()(const Cell* c1, const Cell* c2) const;
};
-// Abstract class representing an elemtary cell of a cell complex.
+// Abstract class representing an elementary cell of a cell complex.
class Cell
{
protected:
@@ -93,6 +93,9 @@ class Cell
virtual int getIndex() const { return _index; };
virtual void setIndex(int index) { _index = index; };
virtual int getNum() { return -1; }
+ virtual int getType() { return -1; }
+ virtual int getTypeForMSH() { return -1; }
+ virtual int getPartition() { return -1; }
virtual void setImmune(bool immune) { _immune = immune; };
virtual bool getImmune() const { return _immune; };
@@ -354,6 +357,9 @@ class ZeroSimplex : public Simplex, public MPoint
int getDim() const { return 0; }
int getNum() { return MPoint::getNum(); }
+ int getType() { return MPoint::getType(); }
+ int getTypeForMSH() { return MPoint::getTypeForMSH(); }
+ int getPartition() { return MPoint::getPartition(); }
int getNumVertices() const { return 1; }
MVertex* getVertex(int vertex) const {return _v[0]; }
int getSortedVertex(int vertex) const {return _v[0]->getNum(); }
@@ -389,6 +395,9 @@ class OneSimplex : public Simplex, public MLine
int getDim() const { return 1; }
int getNum() { return MLine::getNum(); }
+ int getType() { return MLine::getType(); }
+ int getTypeForMSH() { return MLine::getTypeForMSH(); }
+ int getPartition() { return MLine::getPartition(); }
int getNumVertices() const { return 2; }
int getNumFacets() const { return 2; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -435,6 +444,9 @@ class TwoSimplex : public Simplex, public MTriangle
int getDim() const { return 2; }
int getNum() { return MTriangle::getNum(); }
+ int getType() { return MTriangle::getType(); }
+ int getTypeForMSH() { return MTriangle::getTypeForMSH(); }
+ int getPartition() { return MTriangle::getPartition(); }
int getNumVertices() const { return 3; }
int getNumFacets() const { return 3; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -480,6 +492,9 @@ class CQuadrangle : public Cell, public MQuadrangle
int getDim() const { return 2; }
int getNum() { return MQuadrangle::getNum(); }
+ int getType() { return MQuadrangle::getType(); }
+ int getTypeForMSH() { return MQuadrangle::getTypeForMSH(); }
+ int getPartition() { return MQuadrangle::getPartition(); }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -527,6 +542,9 @@ class ThreeSimplex : public Simplex, public MTetrahedron
int getDim() const { return 3; }
int getNum() { return MTetrahedron::getNum(); }
+ int getType() { return MTetrahedron::getType(); }
+ int getTypeForMSH() { return MTetrahedron::getTypeForMSH(); }
+ int getPartition() { return MTetrahedron::getPartition(); }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -573,6 +591,9 @@ class CHexahedron : public Cell, public MHexahedron
int getDim() const { return 3; }
int getNum() { return MHexahedron::getNum(); }
+ int getType() { return MHexahedron::getType(); }
+ int getTypeForMSH() { return MHexahedron::getTypeForMSH(); }
+ int getPartition() { return MHexahedron::getPartition(); }
int getNumVertices() const { return 8; }
int getNumFacets() const { return 6; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -619,6 +640,9 @@ class CPrism : public Cell, public MPrism
int getDim() const { return 3; }
int getNum() { return MPrism::getNum(); }
+ int getType() { return MPrism::getType(); }
+ int getTypeForMSH() { return MPrism::getTypeForMSH(); }
+ int getPartition() { return MPrism::getPartition(); }
int getNumVertices() const { return 6; }
int getNumFacets() const { return 5; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -664,6 +688,9 @@ class CPyramid : public Cell, public MPyramid
int getDim() const { return 3; }
int getNum() { return MPyramid::getNum(); }
+ int getType() { return MPyramid::getType(); }
+ int getTypeForMSH() { return MPyramid::getTypeForMSH(); }
+ int getPartition() { return MPyramid::getPartition(); }
int getNumVertices() const { return 5; }
int getNumFacets() const { return 5; }
MVertex* getVertex(int vertex) const {return _v[vertex]; }
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 85929ec46949a3f8c05391d5e6a67d2a11ff98a0..5f78674b7773fcc301320f2c9ffb698bf8520d7d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -7,6 +7,7 @@
#include "ChainComplex.h"
#include "OS.h"
+#include "PView.h"
#if defined(HAVE_KBIPACK)
@@ -455,7 +456,7 @@ int ChainComplex::getTorsion(int dim, int chainNumber){
}
-Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name, int torsion){
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, GModel* model, std::string name, int torsion){
int i = 0;
for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
@@ -478,6 +479,7 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
_name = name;
_cellComplex = cellComplex;
_torsion = torsion;
+ _model = model;
}
@@ -821,7 +823,7 @@ void Chain::smoothenChain(){
straightenChain(*cit);
}
-
+ eraseNullCells();
double t2 = Cpu();
printf("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
return;
@@ -875,39 +877,32 @@ int Chain::writeChainMSH(const std::string &name){
}
-/*
-void Chain::getData(std::map<int, std::vector<double> > & data){
+void Chain::createPView(){
- if(getSize() == 0) return;
+ std::vector<MElement*> elements;
+ MElementFactory factory;
+
+ std::map<int, std::vector<double> > data;
- //std::list< std::pair<int, Cell*> > cells;
- for(int i = 0; i < getSize(); i++){
- Cell* cell = getCell(i);
- /*
- cells = cell->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell2 = (*it).second;
- std::vector<double> coeff;
- coeff.push_back(getCoeff(i)*(*it).first);
- std::pair<int, std::vector<double> > dataPair = std::make_pair(cell2->getTag(), coeff);
- data.insert(dataPair);
- //printf("%d, %d, \n", cell2->getNum(), (int)coeff.at(0));
-
- }
-
- std::vector<double> coeff;
- coeff.push_back(getCoeff(i));
- std::pair<int, std::vector<double> > dataPair = std::make_pair(cell->getTag(), coeff);
- data.insert(dataPair);
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ int coeff = (*cit).second;
+ std::vector<MVertex*> v = cell->getVertexVector();
+ MElement *e = factory.create(cell->getTypeForMSH(), v, cell->getNum(), cell->getPartition());
+ elements.push_back(e);
+ std::vector<double> coeffs (1,coeff);
+ data.insert(std::make_pair(e->getNum(), coeffs));
}
- //for(std::map<int, std::vector<double> >::iterator it = data.begin(); it != data.end(); it++){
- // printf("%d, %d, \n", (*it).first, (int)(*it).second.at(0));
- //}
+ std::map<int, std::vector<MElement*> > map;
+ map.insert(std::make_pair(0, elements));
+ _model->storeElementsInEntities(map);
+
+ if(!data.empty()) PView *chain = new PView(getName(), "ElementData", getGModel(), data, 0, 1);
- return;
+ return;
}
-*/
+
#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 7739b36c37bc9a2f81f4a1f4ebeea8bb7a73da7c..06e85c8569077d5e08fa9cab280adbf8c176979c 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -141,6 +141,7 @@ class Chain{
std::string _name;
// cell complex this chain belongs to
CellComplex* _cellComplex;
+ GModel* _model;
// torsion coefficient
int _torsion;
@@ -157,13 +158,15 @@ class Chain{
public:
Chain(){}
- Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain", int torsion=0);
+ Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, GModel* model,
+ std::string name="Chain", int torsion=0);
Chain(Chain* chain){
_cells = chain->getCells();
_torsion = chain->getTorsion();
_name = chain->getName();
_cellComplex = chain->getCellComplex();
_dim = chain->getDim();
+ _model = chain->getGModel();
}
~Chain(){}
@@ -175,7 +178,7 @@ class Chain{
//int getCoeff(int i) { return _cells.at(i).second; }
-
+ // remove a cell from this chain
void removeCell(Cell* cell) {
citer it = _cells.find(cell);
if(it != _cells.end()){
@@ -183,6 +186,8 @@ class Chain{
}
return;
}
+
+ // add a cell to this chain
void addCell(Cell* cell, int coeff) {
std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
if(!insert.second && (*insert.first).second == 0) (*insert.first).second = coeff;
@@ -193,26 +198,12 @@ class Chain{
}
return;
}
- /*
- void insertCells(){
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- Cell* cell = (*cit).first;
- std::vector<Cell*> cells;
- if(!_cellComplex->hasCell(cell)){
- _cellComplex->insertCell(cell);
- }
-
- }
- return;
- }
- */
-
+
bool hasCell(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end() && (*it).second != 0) return true;
return false;
- }
-
+ }
Cell* findCell(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end() && (*it).second != 0) return (*it).first;
@@ -225,13 +216,13 @@ class Chain{
}
-
-
int getTorsion() {return _torsion;}
int getDim() {return _dim;}
CellComplex* getCellComplex() {return _cellComplex;}
+ GModel* getGModel() {return _model;}
std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
+ // erase cells from the chain with zero coefficient
void eraseNullCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
@@ -250,6 +241,7 @@ class Chain{
return;
}
+
void deImmuneCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
@@ -262,29 +254,26 @@ class Chain{
//eraseNullCells();
return _cells.size();
}
-
int getNumCells() {
- //int count = 0;
- //for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
- // Cell* cell = (*it).first;
- // count = count + cell->getNumCells();
- //}
- //return count;
return _cells.size();
}
-
// get/set chain name
std::string getName() { return _name; }
void setName(std::string name) { _name=name; }
+ // make local deformations to the chain to make it smoother and smaller
+ // (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
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();
- //void getData(std::map<int, std::vector<double> >& data);
-
};
#endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 7aab300314ef13a9ea3afd663bb68fabf2e9a0f7..fc5f44e668045db87a548bfed469fcbe4da250d2 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -442,6 +442,10 @@ class GModel
void load(std::string fileName);
static void registerBindings(binding *b);
+
+ void storeElementsInEntities(std::map<int, std::vector<MElement*> > &map){
+ _storeElementsInEntities(map);
+ }
};
#endif
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index a4447364402d4521382bea6c46ea0b045a97920b..36c986821104f3f81b6df2df419b7615615c7286 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -9,8 +9,6 @@
#include "Homology.h"
#include "ChainComplex.h"
#include "OS.h"
-#include "PView.h"
-#include "PViewDataGModel.h"
#if defined(HAVE_KBIPACK)
Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
@@ -127,7 +125,7 @@ void Homology::findGenerators(std::string fileName){
convert(i, generator);
std::string name = "H" + dimension + getDomainString() + generator;
- Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name, chains->getTorsion(j,i));
+ 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();
@@ -149,7 +147,7 @@ void Homology::findGenerators(std::string fileName){
convert(i+1, generator);
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, name, 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);
@@ -160,9 +158,11 @@ void Homology::findGenerators(std::string fileName){
}
_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;
}
@@ -251,8 +251,9 @@ void Homology::findDualGenerators(std::string fileName){
convert(i, generator);
std::string name = "H" + dimension + "*" + getDomainString() + generator;
- Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name, chains->getTorsion(j,i));
+ Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, _model, name, chains->getTorsion(j,i));
chain->writeChainMSH(fileName);
+ chain->createPView();
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());
@@ -268,7 +269,7 @@ void Homology::findDualGenerators(std::string fileName){
convert(i+1, generator);
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, name, 1);
+ Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, _model, name, 1);
chain->writeChainMSH(fileName);
if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1;
delete chain;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 8175fc95f22d1fb6bd95325ae78417ac60a833f2..844af0b60836f711cd87b5859ecb5cf3e5dbc6d7 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -78,6 +78,17 @@ 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;
+ }
+ fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
+ fclose(fp);
+ return 1;
+ }
};
#endif
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index af4b553952e783637331c3d18d349bf6302742e2..f15f751a7da5ab99b06586407bfcc174473bb290 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -128,11 +128,11 @@ Physical Surface(75) = {46, 63, 66, 52, 50, 48, 54, 60, 58, 56};
Mesh 3;
// Find generators of relative homology spaces of the domain modulo the four terminals.
-// Save the generator chains to t10_homgen_1.msh.
+// Save the generator chains to t10_homgen.msh.
HomGen("t10_homgen.msh") = {{69}, {70, 71, 72, 73}};
// Find the corresponding cuts.
-// Save the cut chains to t10_homcut_1.msh.
+// Save the cut chains to t10_homcut.msh.
HomCut("t10_homcut.msh") = {{69}, {70, 71, 72, 73}};
// Only find and print the ranks of the relative homology spaces (Betti numbers).
@@ -144,14 +144,10 @@ Mesh.Lines = 0;
Mesh.Triangles = 0;
Mesh.Tetrahedra = 0;
-// Show the results as post-processing views
-Merge "t10_homgen.msh";
-Merge "t10_homcut.msh";
-
// More examples (uncomment):
-// HomGen("t10_homgen.msh") = {{69}, {}}; Merge "t10_homgen.msh";
-// HomGen("t10_homgen.msh") = {{69}, {74}}; Merge "t10_homgen.msh";
-// HomGen("t10_homgen.msh") = {{69}, {75}}; Merge "t10_homgen.msh";
+// HomGen("t10_homgen.msh_1") = {{69}, {}};
+// HomGen("t10_homgen.msh_2") = {{69}, {74}};
+// HomGen("t10_homgen.msh_3") = {{69}, {75}};