diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index cb8818c06f6735e92453f889b50c7e098daa911a..09876bbe70a13170881d0abd81f58b7c7bcb82c4 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -506,7 +506,7 @@ int CellComplex::reduction(int dim){
}
*/
-int CellComplex::reduceComplex(int omit){
+int CellComplex::reduceComplex(){
double t1 = Cpu();
@@ -518,7 +518,7 @@ int CellComplex::reduceComplex(int omit){
int omitted = 0;
_store.clear();
- if(omit > getDim()) omit = getDim();
+ //if(omit > getDim()) omit = getDim();
CellComplex::removeSubdomain();
@@ -592,7 +592,7 @@ void CellComplex::removeSubdomain(){
}
-int CellComplex::coreduceComplex(int omit){
+int CellComplex::coreduceComplex(){
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -669,6 +669,7 @@ void CellComplex::computeBettiNumbers(){
int CellComplex::cocombine(int dim){
+ double t1 = Cpu();
printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -726,9 +727,9 @@ int CellComplex::cocombine(int dim){
}
}
-
- printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ double t2 = Cpu();
+ printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
+ getSize(3), getSize(2), getSize(1), getSize(0), t2-t1);
return count;
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 8a647f1c85290c8eff91c868f8a79c2f190920bd..d8ed040f9a8eb7fed3ceeb7f46e577127feab4ea 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -99,6 +99,17 @@ class Cell
virtual int getSortedVertex(int vertex) const = 0;
virtual std::vector<MVertex*> getVertexVector() const = 0;
+ virtual void restoreCell(){
+ _boundary = _obd;
+ _coboundary = _ocbd;
+ _bdSize = _boundary.size();
+ _cbdSize = _coboundary.size();
+ _combined = false;
+ _index = 0;
+ _immune = false;
+
+ }
+
// returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
// otherwise returns 0
// implementation will vary depending on cell type
@@ -986,6 +997,7 @@ class CellComplex
public:
+ /*
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
_domain = domain;
_subdomain = subdomain;
@@ -994,41 +1006,29 @@ class CellComplex
_cells[cell->getDim()].insert(cell);
}
}
- /*
+
+
CellComplex(CellComplex* cellComplex){
- _domain = cellComplex->_domain;
- _subdomain = cellComplex->_subdomain;
- _boundary = cellComplex->_boundary;
- _domainVertices = cellComplex->_domainVertices;
+ _domain = cellComplex->getDomain();
+ _subdomain = cellComplex->getSubdomain;
+ _boundary = cellComplex->getBoundary;
+ _domainVertices = cellComplex->getDomainVertices;
for(int i = 0; i < 4; i++){
- _betti[i] = cellComplex->_betti[i];
-
- for(citer cit = cellComplex->_cells[i].begin(); cit != cellComplex->_cells[i].end(); cit++){
- Cell* cell = *cit;
- if(i == 0) _cells[i].insert(new ZeroSimplex(*cell));
-
- }
-
- _originalCells[i] = _cells[i];
+ _betti[i] = cellComplex->getBetti(i);
+ _cells[i] = ocells[i];
+ _ocells[i] = ocells[i];
}
- _dim = cellComplex->_dim;
-
- }*/
+ _dim = cellComplex->getDim();
+ }
+ */
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
CellComplex(){}
~CellComplex(){
for(int i = 0; i < 4; i++){
- /*
- for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
- Cell* cell = *cit;
- delete cell;
- //if(cell->isCombined()) delete cell;
- }
- */
_cells[i].clear();
for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
@@ -1041,6 +1041,7 @@ class CellComplex
//emptyTrash();
}
+
void emptyTrash(){
for(std::list<Cell*>::iterator cit = _trash.begin(); cit != _trash.end(); cit++){
Cell* cell = *cit;
@@ -1049,6 +1050,20 @@ class CellComplex
_trash.clear();
}
+ void restoreComplex(){
+ for(int i = 0; i < 4; i++){
+ _betti[i] = 0;
+ _cells[i] = _ocells[i];
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ cell->restoreCell();
+ }
+ }
+ _store.clear();
+ _ecells.clear();
+ _trash.clear();
+ }
+
// get the number of certain dimensional cells
int getSize(int dim){ return _cells[dim].size(); }
int getOrgSize(int dim){ return _ocells[dim].size(); }
@@ -1091,8 +1106,8 @@ class CellComplex
// useful functions for (co)reduction of cell complex
- int reduceComplex(int omit = 0);
- int coreduceComplex(int omit = 0);
+ int reduceComplex();
+ int coreduceComplex();
// add every volume, face and edge its missing boundary cells
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 58a00afadb7bbc4a79a5aeb8fb9146871a90585d..ed1fa49b202e77ea9eea04137c437ea255c54443 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -527,145 +527,344 @@ int Chain::findOrientation(Cell* b, Cell* c){
return 0;
}
-void Chain::smoothenChain(){
+std::map<Cell*, int, Less_Cell> Chain::getBdCellsInChain(Cell* cell){
+
+ std::map<Cell*, int, Less_Cell> cells;
+ std::map<Cell*, int, Less_Cell> boundary = cell->getOrgBd();
+ for(citer cit = boundary.begin(); cit != boundary.end(); cit++){
+ Cell* BdCell = (*cit).first;
+ int BdCellO = (*cit).second;
+ if(hasCell(BdCell) || BdCell->inSubdomain() ) cells.insert( std::make_pair(BdCell, BdCellO) );
+ }
- int start = getSize();
- double t1 = Cpu();
+ return cells;
+}
+
+bool Chain::removeBoundary( std::pair<Cell*, int> cell ){
- if(getDim() == 1){
- for(citer i = _cells.begin(); i != _cells.end(); i++){
- Cell* c1 = (*i).first;
- int c1c = (*i).second;
-
- for(citer j = _cells.begin(); j != _cells.end(); j++){
- Cell* c2 = (*j).first;
- int c2c = (*j).second;
- if ( !(*c2 == *c1) && (c1c == 1 || c1c == -1) && (c2c == 1 || c2c == -1)){
-
- Cell* b = findCommonCbdCell(c1, c2);
- if(b != NULL){
-
- std::pair<Cell*, int> c3p = std::make_pair(b, 0);
- c3p = findRemainingBoundary(b, c1, c2);
- int c1o = findOrientation(b, c1);
- int c2o = findOrientation(b, c2);
-
- int temp1 = c1c - c1o;
- int temp2 = c2c - c2o;
-
- if(c3p.second != 0 && !c3p.first->inSubdomain()) {
-
- Cell* c3 = c3p.first;
- int c3o = c3p.second;
-
- int c3c = -c3o;
- if(temp1 != 0) c3c= c3c*-1;
-
- this->removeCell(c1, c1c);
- this->removeCell(c2, c2c);
- this->addCell(c3, c3c);
- c1c = (*i).second;
- c2c = (*j).second;
-
- }
- }
- }
+ Cell* c1 = cell.first;
+ int c1c = cell.second;
+ if(c1c == 0) return false;
+
+ std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
+ for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
+ Cell* c1CbdCell = (*cit).first;
+
+ std::map<Cell*, int, Less_Cell> cells = getBdCellsInChain(c1CbdCell);
+
+ if( (getDim() == 1 && cells.size() == 3) || (getDim() == 2 && cells.size() == 4)){
+ for(citer cit2 = cells.begin(); cit2 != cells.end(); cit2++){
+ Cell* cell = (*cit2).first;
+ removeCell(cell);
}
+ return true;
}
}
- /*
- else if(getDim() == 2){
- for(citer i = _cells.begin(); i != _cells.end(); i++){
- Cell* c1 = (*i).first;
- int c1c = (*i).second;
+ return false;
+}
+
+bool Chain::straightenChain( std::pair<Cell*, int> cell ){
+
+ Cell* c1 = cell.first;
+ int c1c = cell.second;
+ if(c1c == 0 || c1->getImmune()) return false;
+
+ int c1o = 0;
+
+ Cell* c2 = NULL;
+ int c2c = 0;
+ int c2o = 0;
+
+ Cell* c3 = NULL;
+ int c3c = 0;
+ int c3o = 0;
+
+ Cell* b = NULL;
+
+ std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
+ for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
+ Cell* c1CbdCell = (*cit).first;
+ c1o = (*cit).second;
+
+ /*
+ std::map<Cell*, int, Less_Cell> cells = getBdCellsInChain(c1CbdCell);
+ if((getDim() == 1 && cells.size() != 2) || (getDim() == 2 && cells.size() != 3) ) break;
+ else if( getDim() == 1 && cells.size() == 2){
- /*
- int c1o = 0;
+ }
+ else if( getDim() == 2 && cells.size() == 3){
- Cell* c2 = NULL;
- int c2c = 0;
- int c2o = 0;
+ }
+ */
+
+
+ std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
+ int count = 0;
+ for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
+ Cell* c1CbdBdCell = (*cit2).first;
+ int c1CbdBdCellO = (*cit2).second;
+ int coeff = getCoeff(c1CbdBdCell);
+ if( (coeff != 0 || c1CbdBdCell->inSubdomain() ) && !(*c1CbdBdCell == *c1) && !c1CbdBdCell->getImmune()){
+ if(c1->getDim() == 1){
+ count++;
+ c2 = c1CbdBdCell; c2c = coeff; c2o = c1CbdBdCellO;
+ b = c1CbdCell; break;
+ }
+ else if(c1->getDim() == 2){
+ count++;
+ if(count == 1) { c2 = c1CbdBdCell; c2c = coeff; c2o = c1CbdBdCellO; }
+ else if(count == 2) { c3 = c1CbdBdCell; c3c = coeff; c3o = c1CbdBdCellO; b = c1CbdCell; break;}
+ }
+ }
+ }
+
+ if (b != NULL) break;
+ }
+
+ if(c1->getDim() == 1 &&
+ b != NULL && c2 != NULL && !(*c2 == *c1)){
+
+ int temp1 = c1c - c1o;
+
+ std::pair<Cell*, int> c4p = std::make_pair(b, 0);
+ c4p = findRemainingBoundary(b, c1, c2);
+ Cell* c4 = c4p.first;
+ int c4o = c4p.second;
+
+ if(c4o != 0 && !c2->getImmune() && !c4->getImmune()
+ && ( hasCell(c1) || c1->inSubdomain() ) && (hasCell(c2) || c2->inSubdomain() ) && !hasCell(c4) ){
+
+ int c4c = -c4o;
+ if(temp1 != 0) c4c= c4c*-1;
- Cell* c3 = NULL;
- int c3c = 0;
- int c3o = 0;
+ this->removeCell(c1);
+ this->removeCell(c2);
+ c1->setImmune(false);
+ c2->setImmune(false);
+ c4->setImmune(false);
+ if(!c4->inSubdomain()) this->addCell(c4, c4c);
+ return true;
+ }
+ }
+
+ else if(c1->getDim() == 2 &&
+ b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3)){
+
+ int temp1 = c1c - c1o;
+
+ std::pair<Cell*, int> c4p = std::make_pair(b, 0);
+ c4p = findRemainingBoundary(b, c1, c2, c3);
+ Cell* c4 = c4p.first;
+ int c4o = c4p.second;
+
+ if(c4o != 0 && !c2->getImmune() && !c3->getImmune() && !c4->getImmune()
+ && (hasCell(c1) || c1->inSubdomain()) && (hasCell(c2) || c2->inSubdomain())
+ && (hasCell(c3) || c3->inSubdomain()) && !hasCell(c4)) {
- std::map<Cell*, int, Less_Cell> c1Cbd c1->getOrgCbd();
- for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
- Cell* c1CbdCell = (*cit).first;
- std::map<Cell*, int, Less_Cell> c1CbdBd c1CbdCell->getOrgBd();
-
- for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
- Cell* c1CbdBdCell = (*cit2).first;
- int c1CbdBdCellc = (*cit2).first;
- int count = 0;
- int coeff = getCoeff(c1CbdBdCell);
- if(coeff != 0 && !(*c1CbdBdCell == *c1)){
- count++;
- if(count == 0) { c2 = coeff; c2o=c1CbdBdCellc; }
- if(count == 1) { c3 = coeff; c3o=c1CbdBdCellc; }
- if(count == 2) break;
- }
- }
+ int c4c = -c4o;
+ if(temp1 != 0) c4c= c4c*-1;
+ this->removeCell(c1);
+ this->removeCell(c2);
+ this->removeCell(c3);
+ c1->setImmune(false);
+ c2->setImmune(false);
+ c3->setImmune(false);
+ c4->setImmune(false);
+ if(!c4->inSubdomain()) this->addCell(c4, c4c);
+ return true;
+ }
+
+ }
+ return false;
+}
+
+bool Chain::bendChain( std::pair<Cell*, int> cell ){
+
+ Cell* c1 = cell.first;
+ int c1c = cell.second;
+ if(c1c == 0 || c1->getImmune()) return false;
+ int c1o = 0;
+
+ Cell* c2 = NULL;
+ int c2c = 0;
+ int c2o = 0;
+
+ Cell* c3 = NULL;
+ int c3c = 0;
+ int c3o = 0;
+
+ Cell* b = NULL;
+
+ std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
+ for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
+ Cell* c1CbdCell = (*cit).first;
+ c1o = (*cit).second;
+ std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
+
+ int count = 0;
+ for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
+ Cell* c1CbdBdCell = (*cit2).first;
+ int c1CbdBdCellO = (*cit2).second;
+ if(!hasCell(c1CbdBdCell) && !c1CbdBdCell->getImmune() ){
+ count++;
+ if(count == 1) { c2 = c1CbdBdCell; c2o = c1CbdBdCellO; }
+ else if(count == 2) { c3 = c1CbdBdCell; c3o = c1CbdBdCellO; b = c1CbdCell; break;}
}
+ }
+
+ if (b != NULL) break;
+ else c2 = NULL;
+ }
+
+ if(c1->getDim() == 2 &&
+ b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3) &&
+ (c1c == 1 || c1c == -1) && !c2->getImmune() && !c3->getImmune() ){
+
+ std::pair<Cell*, int> c4p = std::make_pair(b, 0);
+ c4p = findRemainingBoundary(b, c1, c2, c3);
+ int c4c = getCoeff(c4p.first);
+ Cell* c4 = c4p.first;
+
+ int temp1 = c1c - c1o;
+
+ if(c4p.second != 0 && c4c != 0 && !c2->inSubdomain() && !c3->inSubdomain()
+ && hasCell(c1) && hasCell(c4) && !hasCell(c2) && !hasCell(c3)) {
+ c2c = -c2o;
+ c3c = -c3o;
+ if(temp1 != 0) c2c= c2c*-1;
+ if(temp1 != 0) c3c= c3c*-1;
+ this->removeCell(c1);
+ this->removeCell(c4);
+ this->addCell(c2, c2c);
+ this->addCell(c3, c3c);
+ c1->setImmune(false);
+ c2->setImmune(true);
+ c3->setImmune(false);
+ c4->setImmune(false);
+ return true;
+ }
+
+ }
+
+ else if(c1->getDim() == 1 &&
+ b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3) &&
+ (c1c == 1 || c1c == -1) && !c2->getImmune() && !c3->getImmune() ){
+
+ int temp1 = c1c - c1o;
+
+ if(!c2->inSubdomain() && !c3->inSubdomain() && hasCell(c1) && !hasCell(c2) && !hasCell(c3)) {
- for(citer j = _cells.begin(); j != _cells.end(); j++){
- Cell* c2 = (*j).first;
- int c2c = (*j).second;
-
- for(citer k = _cells.begin(); k != _cells.end(); k++){
- Cell* c3 = (*k).first;
- int c3c = (*k).second;
-
- if ( !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3) && (c1c == 1 || c1c == -1) && (c2c == 1 || c2c == -1) && (c3c == 1 || c3c == -1)){
-
- Cell* b = findCommonCbdCell(c1, c2, c3);
- if(b != NULL){
-
- std::pair<Cell*, int> c4p = std::make_pair(b, 0);
- c4p = findRemainingBoundary(b, c1, c2, c3);
- int c1o = findOrientation(b, c1);
- int c2o = findOrientation(b, c2);
- int c3o = findOrientation(b, c3);
-
- int temp1 = c1c - c1o;
- int temp2 = c2c - c2o;
- int temp3 = c3c - c3o;
-
- if(c4p.second != 0 && !c4p.first->inSubdomain()) {
-
- Cell* c4 = c4p.first;
- int c4o = c4p.second;
-
- int c4c = -c4o;
- if(temp1 != 0) c4c= c4c*-1;
-
- this->removeCell(c1, c1c);
- this->removeCell(c2, c2c);
- this->removeCell(c3, c3c);
- this->addCell(c4, c4c);
- c1c = (*i).second;
- c2c = (*j).second;
- c3c = (*k).second;
- }
- }
- }
+ //printf("c2: %d, c3; %d \n", getCoeff(c2), getCoeff(c3));
+
+ c2c = -c2o;
+ c3c = -c3o;
+ if(temp1 != 0) c2c= c2c*-1;
+ if(temp1 != 0) c3c= c3c*-1;
+
+ this->removeCell(c1);
+ this->addCell(c2, c2c);
+ this->addCell(c3, c3c);
+ c1->setImmune(false);
+ c2->setImmune(true);
+ c3->setImmune(false);
+ return true;
+ }
+
+ }
+
+
+ return false;
+}
+
+void Chain::smoothenChain(){
+
+ int start = getSize();
+ double t1 = Cpu();
+ const int MAXROUNDS = 2;
+
+ bool smoothened = true;
+ int round = 0;
+ while(smoothened){
+ round++;
+ smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(straightenChain(*cit)) { smoothened = true; }
+ if(removeBoundary(*cit)) { smoothened = true; }
+ }
+ if(round > MAXROUNDS) smoothened = false;
+ }
+ eraseNullCells();
+
+ if(getDim() == 2){
+ Chain* bestChain = new Chain(this);
+ int before = getSize();
+ deImmuneCells();
+ srand ( time(NULL) );
+ int n = 0;
+ while( n < 10){
+ double t = 1;
+ while(t>0){
+ double tt = 1;
+
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ int random = rand() % 100 + 1;
+ double r = random*t;
+ //printf("random: %d, t: %d \n", random, t);
+ if(r > 80) {
+ bendChain(*cit);
+ //printf("random: %d, t: %g, random*t: %g.\n", random, t, r);
+ //cit = _cells.begin();
+ //tt = tt - 0.1;
+ //if(tt <= 0) tt = 0;
}
}
-
+ eraseNullCells();
+
+ smoothened = true;
+ round = 0;
+ while(smoothened){
+ round++;
+ smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(straightenChain(*cit)){smoothened = true;}
+ if(removeBoundary(*cit)) { smoothened = true; }
+ }
+ if(round > MAXROUNDS) smoothened = false;
+ }
+ eraseNullCells();
+ if(this->getSize() < bestChain->getSize()) bestChain = this;
+ t = t - 0.1;
}
- }*/
+ deImmuneCells();
+ n=n+1;
+ }
+
+ deImmuneCells();
+ smoothened = true;
+ round = 0;
+ while(smoothened){
+ round++;
+ smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(straightenChain(*cit)){smoothened = true;}
+ if(removeBoundary(*cit)) { smoothened = true; }
+ }
+ if(round > MAXROUNDS) smoothened = false;
+ }
eraseNullCells();
+ if(this->getSize() < bestChain->getSize()) bestChain = this;
+ printf("%d-chain simulated annealing removed %d cells. \n", getDim(), before - getSize());
+ }
double t2 = Cpu();
- printf("Smoothened a %d-chain from %d cells to %d cells (%g s). \n", getDim(), start, getSize(), t2-t1);
+ printf("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
return;
}
-
+
int Chain::writeChainMSH(const std::string &name){
//_cellComplex->writeComplexMSH(name);
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index d4527bbb98445a4691d20a5663464d398bee1ec6..e18aec776e8604d7081956ad43c7923079f7dcdd 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -146,6 +146,10 @@ class Chain{
std::pair<Cell*, int> findRemainingBoundary(Cell* b, Cell* c1, Cell* c2, Cell* c3=NULL);
Cell* findCommonCbdCell(Cell* c1, Cell* c2, Cell* c3=NULL);
int findOrientation(Cell* b, Cell* c);
+ bool straightenChain( std::pair<Cell*, int> cell);
+ bool bendChain( std::pair<Cell*, int> cell);
+ bool removeBoundary( std::pair<Cell*, int> cell);
+ std::map<Cell*, int, Less_Cell> getBdCellsInChain(Cell* cell);
public:
Chain(){}
@@ -168,38 +172,21 @@ class Chain{
- void removeCell(Cell* cell, int coeff) {
+ void removeCell(Cell* cell) {
citer it = _cells.find(cell);
if(it != _cells.end()){
- //int coeff2 = (*it).second;
- //(*it).second = coeff2 - coeff;
(*it).second = 0;
- //printf("removed %d, %d :", cell->getNum(), coeff); cell->printCell();
- //_cells.erase(it);
- //if((*it).second == 0) _cells.erase(it);
}
return;
}
void addCell(Cell* cell, int coeff) {
- citer it = _cells.find(cell);
- /*
- if(it != _cells.end()){
- int coeff2 = (*it).second;
- (*it).second = coeff2 + coeff;
- //if((*it).second == 0) _cells.erase(it);
- }
- else{*/
- //printf("added %d, %d :", cell->getNum(), coeff); cell->printCell();
- _cells.insert( std::make_pair( cell, coeff));
- //cell->setImmune(true);
- //cell->clearBoundary();
- //cell->clearCoboundary();
- //_cellComplex->insertCell(cell);
+ std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
+ if(!insert.second && (*insert.first).second == 0) (*insert.first).second = coeff;
+ else if (!insert.second && (*insert.first).second != 0) printf("Error: invalid chain smoothening add! \n");
+
if(!_cellComplex->hasCell(cell)){
_cellComplex->insertCell(cell);
}
-
- //}
return;
}
/*
@@ -221,6 +208,7 @@ class Chain{
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;
@@ -243,6 +231,8 @@ class Chain{
void eraseNullCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
+ //cit++;
+ //_cells.erase(--cit);
_cells.erase(cit);
++cit;
}
@@ -252,12 +242,17 @@ class Chain{
_cells.erase(cit);
cit = _cells.begin();
}
- }
-
-
+ }
return;
}
+ void deImmuneCells(){
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ cell->setImmune(false);
+ }
+ }
+
// number of cells in this chain
int getSize() {
//eraseNullCells();
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 3be2e4a497c810b90c08d07e0dac6eecf32bc3c3..3005ad3f51f500647f13acd5d7268a8998b9c233 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -14,8 +14,6 @@
Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
_model = model;
- _combine = true;
- _omit = 1;
_domain = physicalDomain;
_subdomain = physicalSubdomain;
@@ -81,18 +79,18 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
//_cellComplex->printEuler();
- int omitted = _cellComplex->reduceComplex(_omit);
+ int omitted = _cellComplex->reduceComplex();
//_cellComplex->printEuler();
//_cellComplex->emptyTrash();
- if(getCombine()){
- _cellComplex->combine(3);
- _cellComplex->reduction(2);
- _cellComplex->combine(2);
- _cellComplex->reduction(1);
- _cellComplex->combine(1);
- }
+
+ _cellComplex->combine(3);
+ _cellComplex->reduction(2);
+ _cellComplex->combine(2);
+ _cellComplex->reduction(1);
+ _cellComplex->combine(1);
+
_cellComplex->checkCoherence();
//_cellComplex->printEuler();
@@ -128,7 +126,7 @@ 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, name, chains->getTorsion(j,i));
- //Chain* chain2 = new Chain(chain);
+ Chain* chain2 = new Chain(chain);
//printf("chain %d \n", i);
t1 = Cpu();
int start = chain->getSize();
@@ -140,7 +138,7 @@ void Homology::findGenerators(std::string fileName){
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);
+ chainVector.push_back(chain2);
//delete chain;
}
if(j == _cellComplex->getDim() && _cellComplex->getNumOmitted() > 0){
@@ -174,7 +172,7 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]);
Msg::Info("H3 = %d", HRank[3]);
- if(omitted != 0) Msg::Info("The computation of generators in %d highest dimensions was omitted.", _omit);
+ 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());
@@ -195,7 +193,7 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
- int omitted = _cellComplex->coreduceComplex(_omit);
+ int omitted = _cellComplex->coreduceComplex();
//_cellComplex->emptyTrash();
/*
@@ -209,11 +207,11 @@ void Homology::findDualGenerators(std::string fileName){
_cellComplex->makeDualComplex();
*/
- if(getCombine()){
- _cellComplex->cocombine(0);
- _cellComplex->cocombine(1);
- _cellComplex->cocombine(2);
- }
+
+ _cellComplex->cocombine(0);
+ _cellComplex->cocombine(1);
+ _cellComplex->cocombine(2);
+
//_cellComplex->emptyTrash();
@@ -284,7 +282,7 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("H1* = %d", HRank[1]);
Msg::Info("H2* = %d", HRank[2]);
Msg::Info("H3* = %d", HRank[3]);
- if(omitted != 0) Msg::Info("The computation of %d highest dimension dual generators was omitted.", _omit);
+ 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());
diff --git a/Geo/Homology.h b/Geo/Homology.h
index c6e3a28ffe48029a7a0f7f14cb7f69b9fe78eda5..5308ce105ad5492da59950b21402db645ddf8dc3 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -27,9 +27,6 @@ class Homology
private:
CellComplex* _cellComplex;
GModel* _model;
-
- bool _combine;
- int _omit;
std::vector<int> _domain;
std::vector<int> _subdomain;
@@ -39,9 +36,6 @@ class Homology
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
~Homology(){ delete _cellComplex; }
- bool getCombine() { return _combine; }
- bool setCombine(bool combine) { _combine = combine; return _combine; }
-
void findGenerators(std::string fileName);
void findDualGenerators(std::string fileName);
@@ -49,22 +43,6 @@ class Homology
void swapSubdomain() { _cellComplex->swapSubdomain(); }
- int getOmit() {return _omit; }
- void setOmit(int omit) {
- if(omit == 0) _omit = 0;
- else _omit = 1;
-
- /*
- if(omit > _cellComplex->getDim() || omit < 0) {
- Msg::Error("Invalid number of dimensions to omit. Must be between 0 - %d.", _cellComplex->getDim());
- Msg::Warning("Set to omit 1 dimension.");
- _omit = 1;
- }
- else _omit = omit;
- */
-
- }
-
std::string getDomainString() {
std::string domainString = "({";
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 337b2a79c83a3782e2ca36c4712bcf66988d0329..eb5734d831ffa9472f8a667389e3791730251c75 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -94,7 +94,7 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
Homology* homology = new Homology(m, domain, subdomain);
//if(combine == 0) homology->setCombine(false);
- homology->setOmit(1);
+ //homology->setOmit(1);
//if(swap == 1) homology->swapSubdomain();