From 4e37dafc555df6658bc12877e8163b5dc7ab0fd5 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Tue, 20 Mar 2012 17:14:38 +0000
Subject: [PATCH] Cleaner homology solver API
---
CMakeLists.txt | 2 +-
Geo/Cell.cpp | 5 +-
Geo/CellComplex.cpp | 54 +++++---
Geo/CellComplex.h | 15 ++-
Geo/GModel.cpp | 37 +++---
Geo/Homology.cpp | 207 +++++++++++++++++++++----------
Geo/Homology.h | 47 +++++--
Plugin/HomologyComputation.cpp | 16 +--
demos/homology.geo | 19 ---
utils/api_demos/mainHomology.cpp | 13 +-
10 files changed, 256 insertions(+), 159 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 63c7263b35..1f4de45753 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -93,7 +93,7 @@ set(GMSH_API
Geo/discreteFace.h Geo/discreteRegion.h Geo/SPoint2.h Geo/SPoint3.h
Geo/SVector3.h Geo/STensor3.h Geo/SBoundingBox3d.h Geo/Pair.h Geo/Range.h
Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h
- Geo/Homology.h Geo/partitionEdge.h Geo/CGNSOptions.h Geo/gmshLevelset.h
+ Geo/Homology.h Geo/Chain.h Geo/partitionEdge.h Geo/CGNSOptions.h Geo/gmshLevelset.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h Mesh/meshPartition.h
Mesh/meshGFaceDelaunayInsertion.h Mesh/simple3D.h Mesh/meshPartitionOptions.h
Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h
diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index cec524dac0..165425e4a3 100644
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -351,8 +351,6 @@ void Cell::restoreCell(){
if(it->second.get() == 0) toRemove.push_back(it->first);
}
for(unsigned int i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]);
- _combined = false;
- _immune = false;
}
void Cell::addBoundaryCell(int orientation, Cell* cell, bool other)
@@ -475,7 +473,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co)
_num = ++_globalNum;
_domain = c1->getDomain();
_combined = true;
- _immune = false;
+ _immune = (c1->getImmune() || c2->getImmune());
// cells
c1->getCells(_cells);
@@ -538,6 +536,7 @@ CombinedCell::CombinedCell(std::vector<Cell*>& cells)
// cells
for(unsigned int i = 0; i < cells.size(); i++){
Cell* c = cells.at(i);
+ if(c->getImmune()) _immune = true;
_cells[c] = 1;
}
}
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 4df6ce7a65..cf4021e3fa 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -11,12 +11,14 @@
CellComplex::CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
std::vector<MElement*>& subdomainElements,
+ std::vector<MElement*>& immuneElements,
bool saveOriginalComplex) :
_model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex)
{
_deleteCount = 0;
_insertCells(subdomainElements, 1);
_insertCells(domainElements, 0);
+ _immunizeCells(immuneElements);
int num = 0;
for(int dim = 0; dim < 4; dim++){
@@ -29,6 +31,7 @@ CellComplex::CellComplex(GModel* model,
cell->saveOriginalBd();
}
}
+ _reduced = false;
}
bool CellComplex::_insertCells(std::vector<MElement*>& elements,
@@ -70,6 +73,18 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
return true;
}
+bool CellComplex::_immunizeCells(std::vector<MElement*>& elements)
+{
+ for(unsigned int i=0; i < elements.size(); i++){
+ MElement* element = elements.at(i);
+ Cell* cell = new Cell(element, 0);
+ int dim = cell->getDim();
+ citer cit = _cells[dim].find(cell);
+ if(cit != lastCell(dim)) (*cit)->setImmune(true);
+ delete cell;
+ }
+}
+
CellComplex::~CellComplex()
{
for(int i = 0; i < 4; i++){
@@ -169,8 +184,8 @@ int CellComplex::coreduction(Cell* startCell, bool omit,
s = Q.front();
Q.pop();
Qset.erase(s);
- if(s->getBoundarySize() == 1
- && inSameDomain(s, s->firstBoundary()->first) ){
+ if(s->getBoundarySize() == 1 &&
+ inSameDomain(s, s->firstBoundary()->first)){
s->getBoundary(bd_s);
removeCell(s);
bd_s.begin()->first->getCoboundary(cbd_c);
@@ -186,6 +201,7 @@ int CellComplex::coreduction(Cell* startCell, bool omit,
enqueueCells(cbd_c, Q, Qset);
}
}
+ _reduced = true;
return coreductions;
}
@@ -204,7 +220,9 @@ int CellComplex::reduction(int dim, bool omit,
while(cit != lastCell(dim-1)){
Cell* cell = *cit;
if(cell->getCoboundarySize() == 1 &&
- inSameDomain(cell, cell->firstCoboundary()->first)){
+ inSameDomain(cell, cell->firstCoboundary()->first) &&
+ !cell->getImmune() &&
+ !cell->firstCoboundary()->first->getImmune()){
cit++;
if(dim == getDim() && omit){
omittedCells.push_back(cell->firstCoboundary()->first);
@@ -219,6 +237,7 @@ int CellComplex::reduction(int dim, bool omit,
cit++;
}
}
+ _reduced = true;
return count;
}
@@ -236,8 +255,8 @@ int CellComplex::coreduction(int dim, bool omit,
citer cit = firstCell(dim);
while(cit != lastCell(dim)){
Cell* cell = *cit;
- if( cell->getBoundarySize() == 1
- && inSameDomain(cell, cell->firstBoundary()->first)){
+ if(cell->getBoundarySize() == 1 &&
+ inSameDomain(cell, cell->firstBoundary()->first)) {
++cit;
if(dim-1 == 0 && omit){
omittedCells.push_back(cell->firstBoundary()->first);
@@ -252,6 +271,7 @@ int CellComplex::coreduction(int dim, bool omit,
cit++;
}
}
+ _reduced = true;
return count;
}
@@ -289,7 +309,7 @@ int CellComplex::reduceComplex(bool docombine, bool omit)
}
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
if(docombine) combine(3);
reduction(2, false, empty);
@@ -298,8 +318,8 @@ int CellComplex::reduceComplex(bool docombine, bool omit)
if(docombine) combine(1);
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ _reduced = true;
return 0;
}
@@ -313,6 +333,7 @@ void CellComplex::removeSubdomain()
}
}
for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
+ _reduced = true;
}
int CellComplex::coreduceComplex(bool docombine, bool omit)
@@ -354,7 +375,7 @@ int CellComplex::coreduceComplex(bool docombine, bool omit)
}
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
if(docombine) cocombine(0);
coreduction(1, false, empty);
@@ -364,8 +385,8 @@ int CellComplex::coreduceComplex(bool docombine, bool omit)
coreduction(3, false, empty);
coherent();
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ _reduced = true;
return 0;
}
@@ -389,7 +410,7 @@ int CellComplex::combine(int dim)
Cell* s = Q.front();
Q.pop();
- if(s->getCoboundarySize() == 2){
+ if(s->getCoboundarySize() == 2 && !s->getImmune()){
Cell::biter it = s->firstCoboundary();
int or1 = it->second.get();
Cell* c1 = it->first;
@@ -399,7 +420,7 @@ int CellComplex::combine(int dim)
Cell* c2 = it->first;
if(!(*c1 == *c2) && abs(or1) == abs(or2)
- && inSameDomain(s, c1) && inSameDomain(s, c2)){
+ && inSameDomain(s, c1) && inSameDomain(s, c2)) {
removeCell(s);
c1->getBoundary(bd_c);
@@ -423,7 +444,7 @@ int CellComplex::combine(int dim)
Msg::Debug("Cell complex after combining:");
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
getSize(3), getSize(2), getSize(1), getSize(0));
-
+ _reduced = true;
return count;
}
@@ -485,7 +506,7 @@ int CellComplex::cocombine(int dim)
Msg::Debug("Cell complex after cocombining:");
Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
getSize(3), getSize(2), getSize(1), getSize(0));
-
+ _reduced = true;
return count;
}
@@ -577,8 +598,9 @@ bool CellComplex::restoreComplex()
}
_newcells.clear();
Msg::Info("Restored Cell Complex:");
- Msg::Info(" %d volumes, %d faces, %d edges and %d vertices",
+ Msg::Info("%d volumes, %d faces, %d edges and %d vertices",
getSize(3), getSize(2), getSize(1), getSize(0));
+ _reduced = false;
return true;
}
else {
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 71b362050a..3f149d4195 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -44,8 +44,12 @@ class CellComplex
int _deleteCount;
+ bool _reduced;
+
// for constructor
- bool _insertCells(std::vector<MElement*>& elements, int domain);
+ bool _insertCells(std::vector<MElement*>& elements, int domain);
+
+ bool _immunizeCells(std::vector<MElement*>& elements);
// enqueue cells in queue if they are not there already
void enqueueCells(std::map<Cell*, short int, Less_Cell>& cells,
@@ -63,6 +67,7 @@ class CellComplex
CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
std::vector<MElement*>& subdomainElements,
+ std::vector<MElement*>& immuneElements,
bool saveOriginalComplex=true);
~CellComplex();
@@ -122,9 +127,11 @@ class CellComplex
int reduceComplex(bool docombine=true, bool omit=true);
int coreduceComplex(bool docombine=true, bool omit=true);
- int eulerCharacteristic(){
- return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
- void printEuler(){
+ bool isReduced() const { return _reduced; }
+
+ int eulerCharacteristic() {
+ return getSize(0) - getSize(1) + getSize(2) - getSize(3); }
+ void printEuler() {
printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
// restore the cell complex to its original state before (co)reduction
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 07df48e583..9d96bc5800 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2849,27 +2849,26 @@ void GModel::computeHomology()
Homology* homology = new Homology(this, itp.first->first.first,
itp.first->first.second,
prepareToRestore);
- CellComplex *cellcomplex = homology->createCellComplex();
- if(cellcomplex->getSize(0)){
- for(std::multimap<dpair, std::string>::iterator itt = itp.first;
- itt != itp.second; itt++){
- // restore cell complex to non-reduced state if we are reusing it
- if(itt != itp.first) cellcomplex->restoreComplex();
- std::string type = itt->second;
- if(type == "Homology")
- homology->findHomologyBasis(cellcomplex);
- else if(type == "Cohomology")
- homology->findCohomologyBasis(cellcomplex);
- else
- Msg::Error("Unknown type of homology computation: %s", type.c_str());
+
+ // do not save 0-, and n-chains, where n is the dimension of the model
+ // (usually not needed for anything, available through the plugin)
+ for(std::multimap<dpair, std::string>::iterator itt = itp.first;
+ itt != itp.second; itt++){
+ std::string type = itt->second;
+ if(type == "Homology") {
+ homology->findHomologyBasis();
+ if(this->getDim() != 1) homology->addChainsToModel(1);
+ if(this->getDim() != 2) homology->addChainsToModel(2);
}
- // do not save 0-, and n-chains, where n is the dimension of the model
- // (usually not needed for anything, available through the plugin)
- if(this->getDim() != 1) homology->addChainsToModel(1);
- if(this->getDim() != 2) homology->addChainsToModel(2);
- _pruneMeshVertexAssociations();
+ else if(type == "Cohomology") {
+ homology->findCohomologyBasis();
+ if(this->getDim() != 1) homology->addCochainsToModel(1);
+ if(this->getDim() != 2) homology->addCochainsToModel(2);
+ }
+ else
+ Msg::Error("Unknown type of homology computation: %s", type.c_str());
}
- delete cellcomplex;
+ _pruneMeshVertexAssociations();
delete homology;
}
#else
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 8b37fc422c..9f9b8ccfd8 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -20,7 +20,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
bool combine, bool omit, bool smoothen) :
_model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
_saveOrig(saveOrig),
- _combine(combine), _omit(omit), _smoothen(smoothen)
+ _combine(combine), _omit(omit), _smoothen(smoothen), _cellComplex(NULL),
+ _homologyComputed(false), _cohomologyComputed(false)
{
_fileName = "";
@@ -61,89 +62,121 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
}
}
}
+ for(unsigned int i = 0; i < _imdomain.size(); i++){
+ for(int j = 0; j < 4; j++){
+ it = groups[j].find(_imdomain.at(i));
+ if(it != groups[j].end()){
+ std::vector<GEntity*> physicalGroup = (*it).second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ _immuneEntities.push_back(physicalGroup.at(k));
+ }
+ }
+ }
+ }
}
-CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
- std::vector<GEntity*>& subdomainEntities)
+void Homology::_createCellComplex()
{
Msg::StatusBar(2, true, "Creating cell complex...");
double t1 = Cpu();
- if(domainEntities.empty()) Msg::Error("Domain is empty");
- if(subdomainEntities.empty()) Msg::Info("Subdomain is empty");
+ if(_domainEntities.empty()) Msg::Error("Domain is empty");
+ if(_subdomainEntities.empty()) Msg::Info("Subdomain is empty");
std::vector<MElement*> domainElements;
std::vector<MElement*> subdomainElements;
- for(unsigned int j=0; j < domainEntities.size(); j++) {
- for(unsigned int i=0; i < domainEntities.at(j)->getNumMeshElements(); i++){
- MElement* element = domainEntities.at(j)->getMeshElement(i);
+ std::vector<MElement*> immuneElements;
+ for(unsigned int j=0; j < _domainEntities.size(); j++) {
+ for(unsigned int i=0; i < _domainEntities.at(j)->getNumMeshElements(); i++){
+ MElement* element = _domainEntities.at(j)->getMeshElement(i);
domainElements.push_back(element);
}
}
- for(unsigned int j=0; j < subdomainEntities.size(); j++) {
- for(unsigned int i=0; i < subdomainEntities.at(j)->getNumMeshElements();
+ for(unsigned int j=0; j < _subdomainEntities.size(); j++) {
+ for(unsigned int i=0; i < _subdomainEntities.at(j)->getNumMeshElements();
i++){
- MElement* element = subdomainEntities.at(j)->getMeshElement(i);
+ MElement* element = _subdomainEntities.at(j)->getMeshElement(i);
subdomainElements.push_back(element);
}
}
+ for(unsigned int j=0; j < _immuneEntities.size(); j++) {
+ for(unsigned int i=0; i < _immuneEntities.at(j)->getNumMeshElements();
+ i++){
+ MElement* element = _immuneEntities.at(j)->getMeshElement(i);
+ immuneElements.push_back(element);
+ }
+ }
- CellComplex* cellComplex = new CellComplex(_model,
- domainElements,
- subdomainElements,
- _saveOrig);
+ if(_cellComplex != NULL) delete _cellComplex;
+ _cellComplex = new CellComplex(_model,
+ domainElements,
+ subdomainElements,
+ immuneElements,
+ _saveOrig);
- if(cellComplex->getSize(0) == 0){
+ if(_cellComplex->getSize(0) == 0) {
Msg::Error("Cell Complex is empty: check the domain and the mesh");
}
double t2 = Cpu();
Msg::StatusBar(2, true, "Done creating cell complex (%g s)", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
- return cellComplex;
+ _cellComplex->getSize(3), _cellComplex->getSize(2),
+ _cellComplex->getSize(1), _cellComplex->getSize(0));
}
-CellComplex* Homology::createCellComplex()
+void Homology::_deleteChains()
{
- return createCellComplex(_domainEntities, _subdomainEntities);
+ for(int j = 0; j < 4; j ++) {
+ for(unsigned int i = 0; i < _chains[j].size(); i++) {
+ delete _chains[j].at(i);
+ _chains[j].clear();
+ }
+ }
+ _homologyComputed = false;
}
-Homology::~Homology()
+void Homology::_deleteCochains()
{
- for(unsigned int i = 0; i < _chains.size(); i++)
- delete _chains.at(i);
+ for(int j = 0; j < 4; j ++) {
+ for(unsigned int i = 0; i < _cochains[j].size(); i++) {
+ delete _cochains[j].at(i);
+ _cochains[j].clear();
+ }
+ }
+ _cohomologyComputed = false;
}
-void Homology::findHomologyBasis(CellComplex* cellComplex)
+Homology::~Homology()
{
- bool ownComplex = false;
- if(cellComplex==NULL){
- cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
- ownComplex = true;
- }
+ if(_cellComplex != NULL) delete _cellComplex;
+ _deleteChains();
+ _deleteCochains();
+}
+void Homology::findHomologyBasis()
+{
+ if(_cellComplex == NULL) _createCellComplex();
+ if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
Msg::StatusBar(2, true, "Reducing cell complex...");
double t1 = Cpu();
- int omitted = cellComplex->reduceComplex(_combine, _omit);
+ int omitted = _cellComplex->reduceComplex(_combine, _omit);
double t2 = Cpu();
Msg::StatusBar(2, true, "Done reducing cell complex (%g s)", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ _cellComplex->getSize(3), _cellComplex->getSize(2),
+ _cellComplex->getSize(1), _cellComplex->getSize(0));
Msg::StatusBar(2, true, "Computing homology space basis ...");
t1 = Cpu();
- ChainComplex chainComplex = ChainComplex(cellComplex);
+ ChainComplex chainComplex = ChainComplex(_cellComplex);
chainComplex.computeHomology();
t2 = Cpu();
Msg::StatusBar(2, true, "Done computing homology space basis (%g s)", t2 - t1);
- int dim = cellComplex->getDim();
std::string domain = _getDomainString(_domain, _subdomain);
+ _deleteChains();
int HRank[4];
for(int j = 0; j < 4; j++){
HRank[j] = 0;
@@ -159,7 +192,7 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
if(!chain.empty()) {
- _createChain(chain, name);
+ _createChain(chain, name, false);
HRank[j] = HRank[j] + 1;
if(torsion != 1){
Msg::Warning("H%d %d has torsion coefficient %d!",
@@ -171,8 +204,6 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
if(_fileName != "") writeBasisMSH();
- if(ownComplex) delete cellComplex;
-
Msg::Info("Ranks of domain %shomology spaces:", domain.c_str());
Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]);
@@ -182,38 +213,35 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d",
HRank[0], HRank[1], HRank[2], HRank[3]);
+ _homologyComputed = true;
}
-void Homology::findCohomologyBasis(CellComplex* cellComplex)
+void Homology::findCohomologyBasis()
{
- bool ownComplex = false;
- if(cellComplex==NULL){
- cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
- ownComplex = true;
- }
+ if(_cellComplex == NULL) _createCellComplex();
+ if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
Msg::StatusBar(2, true, "Reducing cell complex...");
double t1 = Cpu();
- int omitted = cellComplex->coreduceComplex(_combine, _omit);
+ int omitted = _cellComplex->coreduceComplex(_combine, _omit);
double t2 = Cpu();
Msg::StatusBar(2, true, "Done reducing cell complex (%g s)", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ _cellComplex->getSize(3), _cellComplex->getSize(2),
+ _cellComplex->getSize(1), _cellComplex->getSize(0));
Msg::StatusBar(2, true, "Computing cohomology space basis ...");
t1 = Cpu();
- ChainComplex chainComplex = ChainComplex(cellComplex);
+ ChainComplex chainComplex = ChainComplex(_cellComplex);
chainComplex.transposeHMatrices();
chainComplex.computeHomology(true);
t2 = Cpu();
Msg::StatusBar(2, true, "Done computing cohomology space basis (%g s)", t2- t1);
- int dim = cellComplex->getDim();
std::string domain = _getDomainString(_domain, _subdomain);
+ _deleteCochains();
int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0;
for(int j = 3; j > -1; j--){
@@ -231,7 +259,7 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
if(!chain.empty()) {
- _createChain(chain, name);
+ _createChain(chain, name, true);
HRank[j] = HRank[j] + 1;
if(torsion != 1){
Msg::Warning("H%d* %d has torsion coefficient %d!", j, i, torsion);
@@ -242,8 +270,6 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
if(_fileName != "") writeBasisMSH();
- if(ownComplex) delete cellComplex;
-
Msg::Info("Ranks of domain %scohomology spaces:", domain.c_str());
Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]);
@@ -253,19 +279,71 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
Msg::StatusBar(2, false, "H0*: %d, H1*: %d, H2*: %d, H3*: %d",
HRank[0], HRank[1], HRank[2], HRank[3]);
+ _cohomologyComputed = true;
}
void Homology::addChainsToModel(int dim, bool post)
{
- for(unsigned int i = 0; i < _chains.size(); i++) {
- if(dim != -1 && _chains.at(i)->getDim() == dim)
- _chains.at(i)->addToModel(this->getModel(), post);
- else if(dim == -1) _chains.at(i)->addToModel(this->getModel(), post);
- }
+ if(!_homologyComputed) Msg::Warning("Homology is not computed");
+ if(dim == -1)
+ for(int j = 0; j < 4; j++)
+ for(unsigned int i = 0; i < _chains[j].size(); i++)
+ _chains[j].at(i)->addToModel(this->getModel(), post);
+ else if (dim > -1 && dim < 4)
+ for(unsigned int i = 0; i < _chains[dim].size(); i++)
+ _chains[dim].at(i)->addToModel(this->getModel(), post);
+}
+
+void Homology::addCochainsToModel(int dim, bool post)
+{
+ if(!_cohomologyComputed) Msg::Warning("Cohomology is not computed");
+ if(dim == -1)
+ for(int j = 0; j < 4; j++)
+ for(unsigned int i = 0; i < _cochains[j].size(); i++)
+ _cochains[j].at(i)->addToModel(this->getModel(), post);
+ else if (dim > -1 && dim < 4)
+ for(unsigned int i = 0; i < _cochains[dim].size(); i++)
+ _cochains[dim].at(i)->addToModel(this->getModel(), post);
+}
+
+void Homology::getHomologyBasis(int dim, std::vector<Chain<int> >& hom)
+{
+ if(dim < 0 || dim > 3) return;
+ if(!_homologyComputed) findHomologyBasis();
+
+ hom.resize(_chains[dim].size());
+ for(unsigned int i = 0; i < _chains[dim].size(); i++)
+ hom[i] = *_chains[dim].at(i);
+}
+
+void Homology::getCohomologyBasis(int dim, std::vector<Chain<int> >& coh)
+{
+ if(dim < 0 || dim > 3) return;
+ if(!_cohomologyComputed) findCohomologyBasis();
+
+ coh.resize(_cochains[dim].size());
+ for(unsigned int i = 0; i < _cochains[dim].size(); i++)
+ coh[i] = *_cochains[dim].at(i);
+}
+
+int Homology::betti(int dim)
+{
+ if(dim < 0 || dim > 3) return 0;
+ if(_homologyComputed) return _chains[dim].size();
+ if(_cohomologyComputed) return _cochains[dim].size();
+
+ findHomologyBasis();
+ return _chains[dim].size();
+}
+
+int Homology::eulerCharacteristic()
+{
+ if(_cellComplex == NULL) _createCellComplex();
+ return _cellComplex->eulerCharacteristic();
}
void Homology::_createChain(std::map<Cell*, int, Less_Cell>& preChain,
- std::string name)
+ std::string name, bool co)
{
Chain<int>* chain = new Chain<int>();
chain->setName(name);
@@ -279,7 +357,8 @@ void Homology::_createChain(std::map<Cell*, int, Less_Cell>& preChain,
cell->getMeshVertices(v);
chain->addElemChain(ElemChain(cell->getDim(), v), coeff);
}
- _chains.push_back(chain);
+ if(co) _cochains[chain->getDim()].push_back(chain);
+ else _chains[chain->getDim()].push_back(chain);
}
std::string Homology::_getDomainString(const std::vector<int>& domain,
@@ -331,16 +410,14 @@ void Homology::storeCells(CellComplex* cellComplex, int dim)
for(CellComplex::citer cit = cellComplex->firstCell(dim);
cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
-
std::map<Cell*, int, Less_Cell > cells;
cell->getCells(cells);
for(Cell::citer it = cells.begin(); it != cells.end(); it++){
Cell* subCell = it->first;
-
std::vector<MVertex*> v;
- cell->getMeshVertices(v);
+ subCell->getMeshVertices(v);
- MElement* e = factory.create(cell->getTypeMSH(), v);
+ MElement* e = factory.create(subCell->getTypeMSH(), v);
elements.push_back(e);
}
}
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 6f396fe20f..16bf789646 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -30,9 +30,14 @@ class Homology
// physical group IDs
std::vector<int> _domain;
std::vector<int> _subdomain;
+ std::vector<int> _imdomain;
// corresponding geometrical entities
std::vector<GEntity*> _domainEntities;
std::vector<GEntity*> _subdomainEntities;
+ std::vector<GEntity*> _immuneEntities;
+
+ // save original cell complex
+ bool _saveOrig;
// use cell combining
bool _combine;
@@ -41,14 +46,18 @@ class Homology
// use chain smoothning
bool _smoothen;
- // save original cell complex
- bool _saveOrig;
-
// file name to store the results
std::string _fileName;
+ // cell complex of the domain
+ CellComplex* _cellComplex;
+
+ bool _homologyComputed;
+ bool _cohomologyComputed;
+
// resulting chains
- std::vector<Chain<int>*> _chains;
+ std::vector<Chain<int>*> _chains[4];
+ std::vector<Chain<int>*> _cochains[4];
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
@@ -56,35 +65,47 @@ class Homology
std::string _getDomainString(const std::vector<int>& domain,
const std::vector<int>& subdomain);
+ // construct the cell complex
+ void _createCellComplex();
+
// create and store a chain from homology solver result
void _createChain(std::map<Cell*, int, Less_Cell>& preChain,
- std::string name);
+ std::string name, bool co);
+ void _deleteChains();
+ void _deleteCochains();
public:
+ // Determine domain and relative subdomain of (co)homology computation
Homology(GModel* model, std::vector<int> physicalDomain,
std::vector<int> physicalSubdomain,
bool saveOrig=true,
bool combine=true, bool omit=true, bool smoothen=true);
~Homology();
- // create a cell complex from a mesh in geometrical entities of Gmsh
- CellComplex* createCellComplex(std::vector<GEntity*>& domainEntities,
- std::vector<GEntity*>& subdomainEntities);
- CellComplex* createCellComplex();
-
GModel* getModel() const { return _model; }
void setFileName(std::string fileName) { _fileName = fileName; }
- // find the generators/dual generarators of homology spaces,
- void findHomologyBasis(CellComplex* cellComplex=NULL);
- void findCohomologyBasis(CellComplex* cellComplex=NULL);
+ // find the bases of (co)homology spaces
+ void findHomologyBasis();
+ void findCohomologyBasis();
// add chains to Gmsh model
// dim: only add dim-chains if dim != -1
// post: create a post-processing view
void addChainsToModel(int dim=-1, bool post=true);
+ void addCochainsToModel(int dim=-1, bool post=true);
+
+ // get representative chains of a (co)homology space basis
+ void getHomologyBasis(int dim, std::vector<Chain<int> >& hom);
+ void getCohomologyBasis(int dim, std::vector<Chain<int> >& coh);
+
+ // get Betti number
+ int betti(int dim);
+
+ // get Euler characteristic
+ int eulerCharacteristic();
// write the generators to a file
bool writeBasisMSH(bool binary=false);
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 5c38a6e77f..ca7ecbe2e2 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -110,20 +110,16 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
Homology* homology = new Homology(m, domain, subdomain);
homology->setFileName(fileName);
- CellComplex* cc = homology->createCellComplex();
- if(hom != 0){
- homology->findHomologyBasis(cc);
- }
- if(coh != 0){
- cc->restoreComplex();
- homology->findCohomologyBasis(cc);
- }
+
+ if(hom != 0) homology->findHomologyBasis();
+ if(coh != 0) homology->findCohomologyBasis();
+
for(unsigned int i = 0; i < dimsave.size(); i++) {
int dim = dimsave.at(i);
- if(dim > -1 && dim < 4) homology->addChainsToModel(dim);
+ if(dim > -1 && dim < 4 && hom != 0) homology->addChainsToModel(dim);
+ if(dim > -1 && dim < 4 && coh != 0) homology->addCochainsToModel(dim);
}
- delete cc;
delete homology;
return 0;
diff --git a/demos/homology.geo b/demos/homology.geo
index de0500b9a4..a8031ed426 100644
--- a/demos/homology.geo
+++ b/demos/homology.geo
@@ -147,22 +147,3 @@ Cohomology {{69}, {70, 71, 72, 73}};
// Homology = {{69}, {}};
// Homology {{69}, {74}};
// Homology {{74}, {70, 71, 72, 73}};
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index 1697a079f7..b721163118 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -32,26 +32,21 @@ int main(int argc, char **argv)
// initialize
Homology* homology = new Homology(m, domain, subdomain);
- // construct cell complex of the meshed domain
- CellComplex* cc = homology->createCellComplex();
-
// find homology basis elements
- homology->findHomologyBasis(cc);
-
- // restore cell complex for a new run
- cc->restoreComplex();
+ homology->findHomologyBasis();
// find cohomology basis elements
- homology->findCohomologyBasis(cc);
+ homology->findCohomologyBasis();
// add 1 and 2 dimensional result chains to model
homology->addChainsToModel(1);
homology->addChainsToModel(2);
+ homology->addCochainsToModel(1);
+ homology->addCochainsToModel(2);
// write mesh with (co)homology computation result chains
m->writeMSH("model_hom.msh");
- delete cc;
delete homology;
delete m;
GmshFinalize();
--
GitLab