diff --git a/Geo/Cell.h b/Geo/Cell.h
index ccf6213b31d3937356ff2a2993b83016c29bf304..c429f775ddec80ea422ddf76baee9fdd407cccd4 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -84,6 +84,23 @@ class Cell
virtual MElement* getImageMElement() const { return _image; };
+
+ // get the number of vertices this cell has
+ virtual int getNumVertices() const { return _image->getNumVertices(); }
+ // get the number of facets of this cell
+ virtual int getNumFacets() const;
+ // get the vertices on a facet of this cell
+ virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
+ // get boundary cell orientation
+ virtual int getFacetOri(std::vector<MVertex*> &v);
+ virtual int getFacetOri(Cell* cell) {
+ std::vector<MVertex*> v;
+ for(int i = 0; i < cell->getNumVertices(); i++) {
+ v.push_back(cell->getVertex(i));
+ }
+ return getFacetOri(v);
+ }
+
virtual int getDim() const { return _dim; };
virtual int getIndex() const { return _index; };
virtual void setIndex(int index) { _index = index; };
@@ -94,9 +111,6 @@ class Cell
virtual void setDeleteImage(bool deleteImage) {
_deleteImage = deleteImage; };
virtual bool getDeleteImage() const { return _deleteImage; };
-
- // get the number of vertices this cell has
- virtual int getNumVertices() const { return _image->getNumVertices(); }
virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
// restores the cell information to its original state before reduction
@@ -145,22 +159,7 @@ class Cell
// print cell debug info
virtual void printCell();
virtual void printBoundary(bool org=false);
- virtual void printCoboundary(bool org=false);
-
- // get the number of facets of this cell
- virtual int getNumFacets() const;
- // get the vertices on a facet of this cell
- virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
-
- // get boundary cell orientation
- virtual int getFacetOri(std::vector<MVertex*> &v);
- virtual int getFacetOri(Cell* cell) {
- std::vector<MVertex*> v;
- for(int i = 0; i < cell->getNumVertices(); i++) {
- v.push_back(cell->getVertex(i));
- }
- return getFacetOri(v);
- }
+ virtual void printCoboundary(bool org=false);
// tools for combined cells
virtual bool isCombined() { return _combined; }
@@ -183,6 +182,14 @@ class Cell
}
return true;
}
+ /*
+ Cell operator=(const Cell& c2) {
+ Cell cell;
+ cell._ocbd = c2._ocbd;
+
+ return cell;
+ }
+ Cell(const Cell& c2){ *this = c2; }*/
};
// A cell that is a combination of cells of same dimension
@@ -194,18 +201,21 @@ class CombinedCell : public Cell{
// list of cells this cell is a combination of
std::list< std::pair<int, Cell*> > _cells;
- MVertex* getVertex(int vertex) const { return NULL; }
+ MVertex* getVertex(int vertex) const {
+ printf("ERROR: No mesh vertices for combined cell."); }
public:
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
~CombinedCell() {}
- int getNum() const { return 0; }
- int getType() const { return 0; }
- int getTypeForMSH() const { return 0; }
- int getPartition() const { return 0; }
-
+ MElement* getImageMElement() const {
+ printf("ERROR: No image mesh element for combined cell."); }
+ int getNumFacets() const { return 0; }
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {}
+ int getFacetOri(std::vector<MVertex*> &v) { return 0; }
+ int getFacetOri(Cell* cell) { return 0; }
+
int getNumVertices() const { return _vs.size(); }
int getSortedVertex(int vertex) const { return _vs.at(vertex); }
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 49e5b626001d2cfa10883af0279c0881245608ba..e0754a5bd484eadf455e57632f7b7611b26e851a 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -621,4 +621,23 @@ void CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
}
}
}
+
+void CellComplex::restoreComplex()
+{
+ for(int i = 0; i < 4; i++){
+ _cells[i] = _ocells[i];
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ cell->restoreCell();
+ }
+ }
+ for(unsigned int i = 0; i < _newcells.size(); i++){
+ Cell* cell = _newcells.at(i);
+ delete cell;
+ }
+ _newcells.clear();
+ _store.clear();
+}
+
+
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 30df99cfeae572e2fa2cbd8cda364667374560c9..d533fb5fffab6afc155a084e6bb6ca13dba0cf94 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -130,6 +130,8 @@ class CellComplex
int getNumOmitted() { return _store.size(); }
std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
+
+ void restoreComplex();
};
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 7b0cd2bab4cb4d5c95be6d34f7ab99c39c27e83a..33b67e6bf187d2f422700677d9e6f6646bdd814a 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -108,10 +108,14 @@ Homology::~Homology()
}
-void Homology::findGenerators()
+void Homology::findGenerators(CellComplex* cellComplex)
{
- CellComplex* cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
+ bool ownComplex = false;
+ if(cellComplex==NULL){
+ CellComplex* cellComplex = createCellComplex(_domainEntities,
+ _subdomainEntities);
+ ownComplex = true;
+ }
std::string domainString = getDomainString(_domain, _subdomain);
Msg::Info("Reducing the Cell Complex...");
@@ -206,7 +210,7 @@ void Homology::findGenerators()
if(_fileName != "") writeGeneratorsMSH();
- delete cellComplex;
+ if(ownComplex) delete cellComplex;
delete chains;
Msg::Info("Ranks of homology spaces for primal cell complex:");
@@ -221,10 +225,14 @@ void Homology::findGenerators()
HRank[0], HRank[1], HRank[2], HRank[3]);
}
-void Homology::findDualGenerators()
+void Homology::findDualGenerators(CellComplex* cellComplex)
{
- CellComplex* cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
+ bool ownComplex = false;
+ if(cellComplex==NULL){
+ CellComplex* cellComplex = createCellComplex(_domainEntities,
+ _subdomainEntities);
+ ownComplex = true;
+ }
Msg::Info("Reducing Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
@@ -321,7 +329,7 @@ void Homology::findDualGenerators()
if(_fileName != "") writeGeneratorsMSH();
- delete cellComplex;
+ if(ownComplex) delete cellComplex;
delete chains;
Msg::Info("Ranks of homology spaces for the dual cell complex:");
@@ -445,7 +453,7 @@ void Homology::findHomSequence(){
}
if(_fileName != "") writeGeneratorsMSH();
- delete cellComplex;
+ delete cellComplex;
}
std::string Homology::getDomainString(const std::vector<int>& domain,
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 705a27bfff223d442d52a2aba219363fef3ea40e..1412d3afb9d0f6284f6da395a78ae52b61174cef 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -42,8 +42,7 @@ class Homology
std::vector<GEntity*> _domainEntities;
std::vector<GEntity*> _subdomainEntities;
- // generator chains, and their physical group number
- //std::map<int, Chain*> _generators;
+ // physical group numbers of generator chains in model
std::vector<int> _generators;
std::string _fileName;
@@ -54,17 +53,20 @@ class Homology
std::vector<int> physicalSubdomain);
~Homology();
+
CellComplex* createCellComplex(std::vector<GEntity*>& domainEntities,
std::vector<GEntity*>& subdomainEntities);
+ CellComplex* createCellComplex() {
+ return createCellComplex(_domainEntities, _subdomainEntities); }
void setFileName(std::string fileName) { _fileName = fileName; }
// Find the generators/duals of homology spaces,
// or just compute the ranks of homology spaces
- void findGenerators();
- void findDualGenerators();
- void computeRanks() {}
-
+ void findGenerators(CellComplex* cellComplex=NULL);
+ void findDualGenerators(CellComplex* cellComplex=NULL);
+
+ void computeRanks() {}
void findHomSequence();
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 33739fc361737f6d5767478a4c93252d8dac8a81..683305092ed18504e23e16fef78abd7234bed072 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -21,7 +21,7 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "PhysicalGroupForSubdomain2", NULL, 0.},
{GMSH_FULLRC, "ComputeGenerators", NULL, 1.},
{GMSH_FULLRC, "ComputeCuts", NULL, 0.},
- {GMSH_FULLRC, "ComputeRanks", NULL, 0.},
+ //{GMSH_FULLRC, "ComputeRanks", NULL, 0.},
};
StringXString HomologyComputationOptions_String[] = {
@@ -89,29 +89,27 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
int gens = (int)HomologyComputationOptions_Number[4].def;
int cuts = (int)HomologyComputationOptions_Number[5].def;
- int rank = (int)HomologyComputationOptions_Number[6].def;
+ //int rank = (int)HomologyComputationOptions_Number[6].def;
GModel *m = GModel::current();
+ Homology* homology = new Homology(m, domain, subdomain);
+ homology->setFileName(fileName);
+ CellComplex* cc = homology->createCellComplex();
if(gens != 0){
- Homology* homology = new Homology(m, domain, subdomain);
- homology->setFileName(fileName);
- homology->findGenerators();
- delete homology;
+ homology->findGenerators(cc);
}
if(cuts != 0){
- Homology* homology = new Homology(m, domain, subdomain);
- homology->setFileName(fileName);
- homology->findDualGenerators();
- delete homology;
+ cc->restoreComplex();
+ homology->findDualGenerators(cc);
}
- if(rank != 0){
- Homology* homology = new Homology(m, domain, subdomain);
- homology->setFileName(fileName);
+ /*if(rank != 0){
homology->computeRanks();
- delete homology;
- }
+ }*/
+ delete cc;
+ delete homology;
+
return 0;
}
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index 6b6aecc25229581eee025548bdbd4faf0ccbfa3c..00d7639c5c1845b5cd1ce8842d9d29c81772dd63 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -12,7 +12,7 @@
#include <gmsh/GModel.h>
#include <gmsh/MElement.h>
#include <gmsh/CellComplex.h>
-#include <gmsh/ChainComplex.h>
+//#include <gmsh/ChainComplex.h>
#include <gmsh/Homology.h>
int main(int argc, char **argv)
@@ -25,21 +25,31 @@ int main(int argc, char **argv)
// OR
// m->readMSH("model.msh");
- // List of physical regions as domain for homology computation (relative to subdomain).
+ // List of physical regions as domain for homology computation
+ // (relative to subdomain).
std::vector<int> domain;
std::vector<int> subdomain;
- Homology* homology1 = new Homology(m, domain, subdomain);
- std::string fileName = "homology.msh";
- homology->findGenerators();
- homology->setFileName(fileName);
- delete homology1;
-
- Homology* homology2 = new Homology(m, domain, subdomain);
- homology->findDualGenerators();
- homology->setFileName(fileName);
- delete homology2;
-
+ // initialize
+ Homology* homology = new Homology(m, domain, subdomain);
+
+ // save all resulting chains to a file as physical groups
+ homology->setFileName("homology.msh");
+
+ // construct cell complex of the meshed domain
+ CellComplex cc = homology->createCellComplex();
+
+ // find homology generator chains
+ homology->findGenerators(cc);
+
+ // restore cell complex for a new run
+ cc->restoreComplex();
+
+ // find thick cuts
+ homology->findDualGenerators(cc);
+
+ delete cc;
+ delete homology;
delete m;
GmshFinalize();