From ec180db5ca03db44118415a946003c4750748d58 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Thu, 11 Jun 2009 12:41:55 +0000
Subject: [PATCH] Bugfixes to cell combining. I can still find examples where
it fails.
---
Geo/CellComplex.cpp | 65 ++++++++++--------
Geo/CellComplex.h | 122 +++++++++++++++++++++++++--------
Geo/ChainComplex.cpp | 20 ++++--
Geo/ChainComplex.h | 11 ++-
Geo/Homology.cpp | 46 ++++++++++---
Geo/Homology.h | 5 ++
Plugin/HomologyComputation.cpp | 24 ++++---
7 files changed, 213 insertions(+), 80 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index cee51bcecc..d7823ca2fb 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -105,11 +105,12 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
cell->setTag(tag);
tag++;
}
- //_originalCells[i] = _cells[i];
+ _cells2[i] = _cells[i];
_betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
}
+
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
@@ -254,6 +255,8 @@ int OneSimplex::kappa(Cell* tau) const{
void CellComplex::removeCell(Cell* cell){
+ //_trash.insert(cell);
+
_cells[cell->getDim()].erase(cell);
std::list<Cell*> coboundary = cell->getCoboundary();
@@ -269,7 +272,6 @@ void CellComplex::removeCell(Cell* cell){
bdCell->removeCoboundaryCell(cell);
}
- //_trash.push_back(cell);
}
@@ -332,7 +334,6 @@ int CellComplex::coreduction(Cell* generator){
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
-
std::list<Cell*> cbd_c;
int count = 0;
@@ -343,7 +344,9 @@ int CellComplex::reduction(int dim){
Cell* cell = *cit;
cbd_c = cell->getCoboundary();
if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front()) ){
+
removeCell(cell);
+
removeCell(cbd_c.front());
count++;
reduced = true;
@@ -362,24 +365,33 @@ int CellComplex::reduceComplex(bool omitHighdim){
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
+
+
int omitted = 0;
if(count == 0 && omitHighdim){
- removeSubdomain();
+ CellComplex::removeSubdomain();
+ CellComplex::removeSubdomain();
std::set<Cell*, Less_Cell> generatorCells;
-
+
while (getSize(getDim()) != 0){
citer cit = firstCell(getDim());
+
Cell* cell = *cit;
generatorCells.insert(cell);
removeCell(cell);
+
reduction(3);
+
reduction(2);
+
reduction(1);
+
omitted++;
+
}
-
+
for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
Cell* cell = *cit;
cell->clearBoundary();
@@ -390,7 +402,6 @@ int CellComplex::reduceComplex(bool omitHighdim){
}
-
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -422,7 +433,7 @@ int CellComplex::coreduceComplex(bool omitLowdim){
int count = 0;
CellComplex::removeSubdomain();
-
+ CellComplex::removeSubdomain();
for(int dim = 0; dim < 4; dim++){
citer cit = firstCell(dim);
@@ -495,7 +506,7 @@ void CellComplex::computeBettiNumbers(){
void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co){
int dim = c1->getDim();
-
+
std::list< std::pair<int, Cell*> > coboundary1 = c1->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > coboundary2 = c2->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > boundary1 = c1->getOrientedBoundary();
@@ -506,7 +517,7 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cbdCell = (*it).second;
int ori = (*it).first;
cbdCell->removeBoundaryCell(c1);
- cbdCell->addBoundaryCell(ori, newCell);
+ cbdCell->addBoundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
Cell* cbdCell = (*it).second;
@@ -519,16 +530,16 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cell2 = (*it2).second;
if(*cell2 == *cbdCell) old = true;
}
- if(!old) cbdCell->addBoundaryCell(ori, newCell);
+ if(!old) cbdCell->addBoundaryCell(ori, newCell, true);
}
- else cbdCell->addBoundaryCell(ori, newCell);
+ else cbdCell->addBoundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell* > >::iterator it = boundary1.begin(); it != boundary1.end(); it++){
Cell* bdCell = (*it).second;
int ori = (*it).first;
bdCell->removeCoboundaryCell(c1);
- bdCell->addCoboundaryCell(ori, newCell);
+ bdCell->addCoboundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
Cell* bdCell = (*it).second;
@@ -541,14 +552,16 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cell2 = (*it2).second;
if(*cell2 == *bdCell) old = true;
}
- if(!old) bdCell->addCoboundaryCell(ori, newCell);
+ if(!old) bdCell->addCoboundaryCell(ori, newCell, true);
}
- else bdCell->addCoboundaryCell(ori, newCell);
+ else bdCell->addCoboundaryCell(ori, newCell, true);
}
_cells[dim].erase(c1);
_cells[dim].erase(c2);
+ //_trash.insert(c1);
+ //_trash.insert(c2);
_cells[dim].insert(newCell);
@@ -584,6 +597,8 @@ int CellComplex::cocombine(int dim){
if(bd_c.size() == 2 && !(*(bd_c.front().second) == *(bd_c.back().second))
&& inSameDomain(s, bd_c.front().second) && inSameDomain(s, bd_c.back().second) ){
+ int or1 = bd_c.front().first;
+ int or2 = bd_c.back().first;
Cell* c1 = bd_c.front().second;
Cell* c2 = bd_c.back().second;
@@ -591,16 +606,14 @@ int CellComplex::cocombine(int dim){
enqueueCells(cbd_c, Q, Qset);
cbd_c = c2->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
-
- int or1 = bd_c.front().first;
- int or2 = bd_c.back().first;
-
+
removeCell(s);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
replaceCells(c1, c2, newCell, (or1 != or2), true);
cit = firstCell(dim);
count++;
+
}
removeCellQset(s, Qset);
@@ -641,25 +654,23 @@ int CellComplex::combine(int dim){
if(cbd_c.size() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
&& inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second) ){
-
+ int or1 = cbd_c.front().first;
+ int or2 = cbd_c.back().first;
Cell* c1 = cbd_c.front().second;
Cell* c2 = cbd_c.back().second;
-
+
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
bd_c = c2->getBoundary();
enqueueCells(bd_c, Q, Qset);
-
- int or1 = cbd_c.front().first;
- int or2 = cbd_c.back().first;
-
+
removeCell(s);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
replaceCells(c1, c2, newCell, (or1 != or2));
-
cit = firstCell(dim);
//cit++;
count++;
+
}
removeCellQset(s, Qset);
@@ -667,7 +678,7 @@ int CellComplex::combine(int dim){
}
printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
return count;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 7c217e15ff..4a3e8d5877 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -88,6 +88,7 @@ class Cell
for( std::list< std::pair<int, Cell*> >::iterator it= _boundary.begin();it!= _boundary.end();it++){
Cell* cell = (*it).second;
boundary.push_back(cell);
+ if((*it).first == 0) boundary.push_back(cell);
}
return boundary;
}
@@ -97,26 +98,50 @@ class Cell
for( std::list< std::pair<int, Cell*> >::iterator it= _coboundary.begin();it!= _coboundary.end();it++){
Cell* cell = (*it).second;
coboundary.push_back(cell);
+ if((*it).first == 0) coboundary.push_back(cell);
}
return coboundary;
}
- virtual void addBoundaryCell(int orientation, Cell* cell) { _boundary.push_back( std::make_pair(orientation, cell)); }
- virtual void addCoboundaryCell(int orientation, Cell* cell) { _coboundary.push_back( std::make_pair(orientation, cell)); }
+ virtual void addBoundaryCell(int orientation, Cell* cell, bool duplicates=false) {
+ if(!duplicates){
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ int ori2 = (*it).first;
+ if(*cell2 == *cell && !duplicates) return;
+ }
+ }
+ _boundary.push_back( std::make_pair(orientation, cell));
+
+ }
+ virtual void addCoboundaryCell(int orientation, Cell* cell, bool duplicates=false) {
+ if(!duplicates){
+ for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ int ori2 = (*it).first;
+ if(*cell2 == *cell && !duplicates) return;
+ }
+ }
+ _coboundary.push_back( std::make_pair(orientation, cell));
+ }
- virtual void removeBoundaryCell(Cell* cell) {
+ virtual int removeBoundaryCell(Cell* cell) {
+ int count = 0;
for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = (*it).second;
- if(*cell2 == *cell) { _boundary.erase(it); break; }
+ if(*cell2 == *cell) { _boundary.erase(it); count++; it = _boundary.begin(); }
}
+ return count;
}
- virtual void removeCoboundaryCell(Cell* cell) {
+ virtual int removeCoboundaryCell(Cell* cell) {
+ int count = 0;
for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
Cell* cell2 = (*it).second;
- if(*cell2 == *cell) { _coboundary.erase(it); break; }
+ if(*cell2 == *cell) { _coboundary.erase(it); count++; it = _coboundary.begin(); }
}
+ return count;
}
-
+
virtual void clearBoundary() { _boundary.clear(); }
virtual void clearCoboundary() { _coboundary.clear(); }
@@ -187,6 +212,7 @@ class Simplex : public Cell
public:
Simplex() : Cell() {
_combined = false;
+ _index = 0;
}
~Simplex(){}
@@ -413,6 +439,7 @@ class CombinedCell : public Cell{
public:
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
+ _index = c1->getIndex();
_tag = c1->getTag();
_dim = c1->getDim();
_inSubdomain = c1->inSubdomain();
@@ -518,7 +545,9 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
- std::vector<Cell*> _trash;
+ // trashbin for cell pointers removed from cell complex
+ std::set<Cell*, Less_Cell> _cells2[4];
+ std::set<Cell*, Less_Cell> _trash;
//std::set<Cell*, Less_Cell> _originalCells[4];
@@ -531,7 +560,7 @@ class CellComplex
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
- protected:
+ private:
// enqueue cells in queue if they are not there already
void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
@@ -542,6 +571,19 @@ class CellComplex
//virtual void insert_cells(bool subdomain, bool boundary);
void insert_cells(bool subdomain, bool boundary);
+ // remove a cell from this cell complex
+ void removeCell(Cell* cell);
+ void replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co=false);
+ void insertCell(Cell* cell);
+
+ // reduction of this cell complex
+ // removes reduction pairs of cell of dimension dim and dim-1
+ int reduction(int dim);
+
+ // queued coreduction presented in Mrozek's paper
+ int coreduction(Cell* generator);
+
+
public:
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
_domain = domain;
@@ -578,10 +620,19 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
CellComplex(){}
~CellComplex(){
- for(int i = 0; i < _trash.size(); i++){
- Cell* cell = _trash.at(i);
- delete cell;
+ //for(std::set<Cell*, Less_Cell>::iterator it = _trash.begin(); it != _trash.end(); it++){
+ // Cell* cell = *it;
+ // delete cell;
+ //}
+
+ for(int i = 0; i < 4; i++){
+ _cells[i].clear();
+ for(citer cit = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
}
+
}
@@ -593,36 +644,23 @@ class CellComplex
std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
// iterators to the first and last cells of certain dimension
- citer firstCell(int dim) {return _cells[dim].begin(); }
- citer lastCell(int dim) {return _cells[dim].end(); }
+ citer firstCell(int dim) {return _cells[dim].begin(); }
+ citer lastCell(int dim) {return _cells[dim].end(); }
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
- // remove a cell from this cell complex
- void removeCell(Cell* cell);
- void replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co=false);
-
-
- void insertCell(Cell* cell);
// check whether two cells both belong to subdomain or if neither one does
bool inSameDomain(Cell* c1, Cell* c2) const { return
( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
-
-
- // reduction of this cell complex
- // removes reduction pairs of cell of dimension dim and dim-1
- int reduction(int dim);
-
+
// useful functions for (co)reduction of cell complex
int reduceComplex(bool omitHighdim = false);
int coreduceComplex(bool omitlowDim = false);
- // queued coreduction presented in Mrozek's paper
- int coreduction(Cell* generator);
// add every volume, face and edge its missing boundary cells
// void repairComplex(int i=3);
@@ -640,6 +678,13 @@ class CellComplex
int combine(int dim);
int cocombine(int dim);
+ void emptyTrash() {
+ for(std::set<Cell*, Less_Cell>::iterator it = _trash.begin(); it != _trash.end(); it++){
+ Cell* cell = *it;
+ delete cell;
+ }
+ }
+
void computeBettiNumbers();
int getBettiNumber(int i) { if(i > -1 && i < 4) return _betti[i]; else return 0; }
@@ -659,6 +704,27 @@ class CellComplex
}
}
+ void checkCoherence(){
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ std::list< std::pair<int, Cell*> > boundary;
+ for(std::list< std::pair<int, Cell* > >::iterator it = boundary.begin(); it != boundary.end(); it++){
+ Cell* bdCell = (*it).second;
+ citer cit = _cells[bdCell->getDim()].find(bdCell);
+ if(cit == lastCell(bdCell->getDim())) printf("Warning! Incoherent boundary!");
+ }
+ std::list< std::pair<int, Cell*> > coboundary;
+ for(std::list< std::pair<int, Cell* > >::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ Cell* cbdCell = (*it).second;
+ citer cit = _cells[cbdCell->getDim()].find(cbdCell);
+ if(cit == lastCell(cbdCell->getDim())) printf("Warning! Incoherent coboundary!");
+ }
+
+ }
+ }
+ }
+
};
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 769e610b4d..fefc4e1006 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -76,11 +76,21 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
Cell* bdCell = (*it).second;
if(!bdCell->inSubdomain()){
int old_elem = 0;
- //printf("kaka1: %d, kaka2: %d \n", bdCell->getIndex(), cell->getIndex());
- gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
- old_elem = mpz_get_si(elem);
- mpz_set_si(elem, old_elem + (*it).first);
- gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ //printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex());
+ if(bdCell->getIndex() > gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1
+ || cell->getIndex() > gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){
+ printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
+ }
+ else{
+ gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ old_elem = mpz_get_si(elem);
+ mpz_set_si(elem, old_elem + (*it).first);
+ if( (old_elem + (*it).first) > 1 || (old_elem + (*it).first) < -1 ){
+ printf("Warning: Invalid incidence index: %d! HMatrix: %d. \n", (old_elem + (*it).first), dim);
+ mpz_set_si(elem, 0);
+ }
+ gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ }
}
}
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index f5569a52d5..1b41eb7d02 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -77,7 +77,16 @@ class ChainComplex{
}
_dim = 0;
}
- ~ChainComplex(){}
+ ~ChainComplex(){
+ for(int i = 0; i < 5; i++){
+ destroy_gmp_matrix(_HMatrix[i]);
+ destroy_gmp_matrix(_kerH[i]);
+ destroy_gmp_matrix(_codH[i]);
+ destroy_gmp_matrix(_JMatrix[i]);
+ destroy_gmp_matrix(_QMatrix[i]);
+ destroy_gmp_matrix(_Hbasis[i]);
+ }
+ }
int getDim() { return _dim; }
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 8a7b7b1543..18cdc81db9 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -14,6 +14,7 @@
Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
_model = model;
+ _combine = true;
Msg::Info("Creating a Cell Complex...");
double t1 = Cpu();
@@ -74,14 +75,22 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(true);
- _cellComplex->combine(3);
- _cellComplex->combine(2);
- _cellComplex->combine(1);
+ //_cellComplex->checkCoherence();
+ if(getCombine()){
+ _cellComplex->combine(3);
+ _cellComplex->reduceComplex(false);
+ _cellComplex->combine(2);
+ _cellComplex->reduceComplex(false);
+ _cellComplex->combine(1);
+ _cellComplex->reduceComplex(false);
+ }
+ //_cellComplex->checkCoherence();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
@@ -119,29 +128,37 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("Computation of %dD generators was omitted.", _cellComplex->getDim());
-
Msg::Info("Wrote results to %s.", fileName.c_str());
delete chains;
+ printf("H0 = %d \n", HRank[0]);
+ printf("H1 = %d \n", HRank[1]);
+ printf("H2 = %d \n", HRank[2]);
+ printf("H3 = %d \n", HRank[3]);
+
return;
}
void Homology::findThickCuts(std::string fileName){
- _cellComplex->removeSubdomain();
-
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(true);
- _cellComplex->cocombine(0);
- _cellComplex->cocombine(1);
- _cellComplex->cocombine(2);
+ //_cellComplex->checkCoherence();
+ if(getCombine()){
+ _cellComplex->cocombine(0);
+ _cellComplex->cocombine(1);
+ _cellComplex->cocombine(2);
+ }
+ //_cellComplex->checkCoherence();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
@@ -153,6 +170,8 @@ void Homology::findThickCuts(std::string fileName){
Msg::Info("Homology Computation complete (%g s).", t2- t1);
int dim = _cellComplex->getDim();
+
+
int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0;
@@ -175,7 +194,7 @@ void Homology::findThickCuts(std::string fileName){
}
}
-
+
Msg::Info("Ranks of homology groups for dual cell complex:");
Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]);
@@ -183,10 +202,15 @@ void Homology::findThickCuts(std::string fileName){
Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("Computation of %dD thick cuts was omitted.", dim);
-
Msg::Info("Wrote results to %s.", fileName.c_str());
delete chains;
+ printf("H0 = %d \n", HRank[0]);
+ printf("H1 = %d \n", HRank[1]);
+ printf("H2 = %d \n", HRank[2]);
+ printf("H3 = %d \n", HRank[3]);
+
+ return;
}
#endif
diff --git a/Geo/Homology.h b/Geo/Homology.h
index b5f013653b..e53ca8bd8d 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -28,11 +28,16 @@ class Homology
CellComplex* _cellComplex;
GModel* _model;
+ bool _combine;
+
public:
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
~Homology(){ delete _cellComplex; }
+ bool getCombine() { return _combine; }
+ bool setCombine(bool combine) { _combine = combine; }
+
void findGenerators(std::string fileName);
void findThickCuts(std::string fileName);
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index b7d4825630..41d6e8cbfb 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -15,11 +15,14 @@
#if defined(HAVE_KBIPACK)
StringXNumber HomologyComputationOptions_Number[] = {
- {GMSH_FULLRC, "Physical group for domain", NULL, 0.},
- {GMSH_FULLRC, "Physical group for subdomain", NULL, 0.},
+ {GMSH_FULLRC, "1. Physical group for domain", NULL, 0.},
+ {GMSH_FULLRC, "2. Physical group for domain", NULL, 0.},
+ {GMSH_FULLRC, "1. Physical group for subdomain", NULL, 0.},
+ {GMSH_FULLRC, "2. Physical group for subdomain", NULL, 0.},
{GMSH_FULLRC, "Compute generators", NULL, 1.},
{GMSH_FULLRC, "Compute thick cuts", NULL, 0.},
{GMSH_FULLRC, "Swap subdomain", NULL, 0.},
+ {GMSH_FULLRC, "Combine cells", NULL, 1.}
};
StringXString HomologyComputationOptions_String[] = {
@@ -89,22 +92,27 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
std::vector<int> subdomain;
domain.push_back((int)HomologyComputationOptions_Number[0].def);
- subdomain.push_back((int)HomologyComputationOptions_Number[1].def);
+ domain.push_back((int)HomologyComputationOptions_Number[1].def);
+ subdomain.push_back((int)HomologyComputationOptions_Number[2].def);
+ subdomain.push_back((int)HomologyComputationOptions_Number[3].def);
- int gen = (int)HomologyComputationOptions_Number[2].def;
- int cuts = (int)HomologyComputationOptions_Number[3].def;
- int swap = (int)HomologyComputationOptions_Number[4].def;
+ int gens = (int)HomologyComputationOptions_Number[4].def;
+ int cuts = (int)HomologyComputationOptions_Number[5].def;
+ int swap = (int)HomologyComputationOptions_Number[6].def;
+ int combine = (int)HomologyComputationOptions_Number[7].def;
GModel *m = GModel::current();
Homology* homology = new Homology(m, domain, subdomain);
+ if(combine == 0) homology->setCombine(false);
+
if(swap == 1) homology->swapSubdomain();
- if(gen == 1 && cuts != 1) {
+ if(gens == 1 && cuts != 1) {
homology->findGenerators(fileName);
GmshMergeFile(fileName);
}
- else if(cuts == 1 && gen != 1) {
+ else if(cuts == 1 && gens != 1) {
homology->findThickCuts(fileName);
GmshMergeFile(fileName);
}
--
GitLab