diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 03031c0b017f46a935717c83626833004a0c069a..176260c60b7c6b70ca4959c0c2e6c2e50359c804 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -114,7 +114,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
-
+
// works only for simplcial meshes at the moment!
std::vector<GEntity*> domain;
@@ -215,8 +215,8 @@ int Simplex::kappa(Cell* tau) const{
int value=1;
for(int i = 0; i < this->getNumFacets(); i++){
- std::vector<int> vTau = tau->getVertexVector();
- std::vector<int> v;
+ std::vector<MVertex*> vTau = tau->getVertexVector();
+ std::vector<MVertex*> v;
this->getFacetVertices(i, v);
value = -1;
@@ -239,7 +239,8 @@ int Simplex::kappa(Cell* tau) const{
return 0;
}
-
+*/
+/*
int OneSimplex::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
@@ -255,25 +256,28 @@ int OneSimplex::kappa(Cell* tau) const{
void CellComplex::removeCell(Cell* cell){
- //_trash.insert(cell);
+ std::list< std::pair< int, Cell*> > coboundary = cell->getOrientedCoboundary();
+ std::list< std::pair< int, Cell*> > boundary = cell->getOrientedBoundary();
+ //std::list<Cell*> boundary = cell->getBoundary();
+ //std::list<Cell*> coboundary = cell->getCoboundary();
- _cells[cell->getDim()].erase(cell);
-
- std::list<Cell*> coboundary = cell->getCoboundary();
- std::list<Cell*> boundary = cell->getBoundary();
-
- for(std::list<Cell*>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
- Cell* cbdCell = *it;
+ for(std::list< std::pair< int, Cell*> >::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ //for(std::list<Cell*>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ Cell* cbdCell = (*it).second;
+ //Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell);
}
- for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
- Cell* bdCell = *it;
+
+ for(std::list< std::pair< int, Cell*> >::iterator it = boundary.begin(); it != boundary.end(); it++){
+ //for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
+ Cell* bdCell = (*it).second;
+ //Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell);
}
-
- //cell->clearBoundary();
- //cell->clearCoboundary();
+
+ _cells[cell->getDim()].erase(cell);
+
}
@@ -346,11 +350,11 @@ int CellComplex::reduction(int dim){
for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
Cell* cell = *cit;
cbd_c = cell->getCoboundary();
- if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front()) ){
-
- removeCell(cell);
+ if( cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front()) ){
removeCell(cbd_c.front());
+ removeCell(cell);
+ //removeCell(cbd_c.front());
count++;
reduced = true;
}
@@ -385,11 +389,8 @@ int CellComplex::reduceComplex(bool omitHighdim){
removeCell(cell);
reduction(3);
-
reduction(2);
-
reduction(1);
-
omitted++;
}
@@ -451,13 +452,13 @@ int CellComplex::coreduceComplex(bool omitLowdim){
int omitted = 0;
if(count == 0 && omitLowdim){
std::set<Cell*, Less_Cell> generatorCells;
-
while (getSize(0) != 0){
citer cit = firstCell(0);
Cell* cell = *cit;
generatorCells.insert(cell);
removeCell(cell);
- coreduction(cell);
+ coreduceComplex(false);
+ //coreduction(cell);
omitted++;
}
@@ -505,111 +506,6 @@ 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();
- std::list< std::pair<int, Cell*> > boundary2 = c2->getOrientedBoundary();
-
-
- for(std::list< std::pair<int, Cell*> >::iterator it = coboundary1.begin(); it != coboundary1.end(); it++){
- Cell* cbdCell = (*it).second;
- int ori = (*it).first;
- cbdCell->removeBoundaryCell(c1);
- cbdCell->addBoundaryCell(ori, newCell, true);
- }
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = coboundary1.begin(); it != coboundary1.end(); it++){
- Cell* cbdCell = (*it).second;
- int ori = (*it).first;
- cbdCell->removeBoundaryCell(c1);
- if(!co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = coboundary2.begin(); it2 != coboundary2.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cbdCell) old = true;
- }
- if(!old) cbdCell->addBoundaryCell(ori, newCell, true);
- }
- else cbdCell->addBoundaryCell(ori, newCell, true);
- }
- */
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
- Cell* cbdCell = (*it).second;
- int ori = (*it).first;
- if(!orMatch) ori = -1*ori;
- cbdCell->removeBoundaryCell(c2);
- if(!co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = coboundary1.begin(); it2 != coboundary1.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cbdCell) old = true;
- }
- if(!old) cbdCell->addBoundaryCell(ori, newCell, true);
- }
- 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, 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);
- if(co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = boundary2.begin(); it2 != boundary2.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *bdCell) old = true;
- }
- if(!old) bdCell->addCoboundaryCell(ori, newCell, true);
- }
- else bdCell->addCoboundaryCell(ori, newCell, true);
-
- }
- */
- /*
- for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
- Cell* bdCell = (*it).second;
- int ori = (*it).first;
- if(!orMatch) ori = -1*ori;
- bdCell->removeCoboundaryCell(c2);
- if(co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = boundary1.begin(); it2 != boundary1.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *bdCell) old = true;
- }
- if(!old) bdCell->addCoboundaryCell(ori, newCell, true);
- }
- else bdCell->addCoboundaryCell(ori, newCell, true);
-
- }
- */
- _cells[dim].erase(c1);
- _cells[dim].erase(c2);
- //removeCell(c1);
- //removeCell(c2);
- //_trash.insert(c1);
- //_trash.insert(c2);
- _cells[dim].insert(newCell);
-
-
-
- return;
-}
-
int CellComplex::cocombine(int dim){
printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -635,14 +531,15 @@ int CellComplex::cocombine(int dim){
bd_c = s->getOrientedBoundary();
- 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) ){
+ if(s->getBoundarySize() == 2 && !(*(bd_c.front().second) == *(bd_c.back().second))
+ && inSameDomain(s, bd_c.front().second) && inSameDomain(s, bd_c.back().second)
+ && bd_c.front().second->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
+ && bd_c.back().second->getNumVertices() < getSize(dim)){
int or1 = bd_c.front().first;
int or2 = bd_c.back().first;
Cell* c1 = bd_c.front().second;
Cell* c2 = bd_c.back().second;
- if( c1->getNumVertices() < getSize(dim) && c2->getNumVertices() < getSize(dim) ){
removeCell(s);
@@ -652,14 +549,12 @@ int CellComplex::cocombine(int dim){
enqueueCells(cbd_c, Q, Qset);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
- //replaceCells(c1, c2, newCell, (or1 != or2), true);
- _cells[dim].erase(c1);
- _cells[dim].erase(c2);
+ removeCell(c1);
+ removeCell(c2);
_cells[dim].insert(newCell);
cit = firstCell(dim);
count++;
- }
}
removeCellQset(s, Qset);
@@ -687,44 +582,39 @@ int CellComplex::combine(int dim){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
-
Cell* cell = *cit;
bd_c = cell->getBoundary();
enqueueCells(bd_c, Q, Qset);
while(Q.size() != 0){
Cell* s = Q.front();
- Q.pop();
+ Q.pop();
cbd_c = s->getOrientedCoboundary();
- 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)){
+ if(s->getCoboundarySize() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
+ && inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second)
+ && cbd_c.front().second->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
+ && cbd_c.back().second->getNumVertices() < getSize(dim)){
int or1 = cbd_c.front().first;
int or2 = cbd_c.back().first;
Cell* c1 = cbd_c.front().second;
Cell* c2 = cbd_c.back().second;
- if( c1->getNumVertices() < getSize(dim) && c2->getNumVertices() < getSize(dim)){
removeCell(s);
+
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
bd_c = c2->getBoundary();
enqueueCells(bd_c, Q, Qset);
-
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
- _cells[dim].erase(c1);
- _cells[dim].erase(c2);
- _cells[dim].insert(newCell);
- //replaceCells(c1, c2, newCell, (or1 != or2));
- //removeCell(s);
+ removeCell(c1);
+ removeCell(c2);
+ std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
+ if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n");
cit = firstCell(dim);
- //cit++;
count++;
-
- }
-
}
removeCellQset(s, Qset);
@@ -868,4 +758,61 @@ void CellComplex::printComplex(int dim){
}
}
+bool CellComplex::checkCoherence(){
+ bool coherent = true;
+ 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 = cell->getOrientedBoundary();
+ for(std::list< std::pair<int, Cell* > >::iterator it = boundary.begin(); it != boundary.end(); it++){
+ Cell* bdCell = (*it).second;
+ int ori = (*it).first;
+ citer cit = _cells[bdCell->getDim()].find(bdCell);
+ if(cit == lastCell(bdCell->getDim())){
+ printf("Warning! Boundary cell not in cell complex! Boundary removed. \n", cell->getDim());
+ //printf(" "); cell->printCell();
+ //printf(" "); bdCell->printCell();
+ cell->removeBoundaryCell(bdCell);
+ coherent = false;
+ }
+ if(!bdCell->hasCoboundary(cell)){
+ printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
+ //printf(" "); cell->printCell();
+ //printf(" "); cell->printBoundary();
+ //printf(" "); bdCell->printCell();
+ //printf(" "); bdCell->printCoboundary();
+ bdCell->addCoboundaryCell(ori, cell);
+ coherent = false;
+ }
+
+ }
+ std::list< std::pair<int, Cell*> > coboundary = cell->getOrientedCoboundary();
+ for(std::list< std::pair<int, Cell* > >::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ Cell* cbdCell = (*it).second;
+ int ori = (*it).first;
+ citer cit = _cells[cbdCell->getDim()].find(cbdCell);
+ if(cit == lastCell(cbdCell->getDim())){
+ printf("Warning! Coboundary cell not in cell complex! Coboundary removed. \n", cell->getDim());
+ //printf(" "); cell->printCell();
+ //printf(" "); cbdCell->printCell();
+ cell->removeCoboundaryCell(cbdCell);
+ coherent = false;
+ }
+ if(!cbdCell->hasBoundary(cell)){
+ printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
+ //printf(" "); cell->printCell();
+ //printf(" "); cell->printCoboundary();
+ //printf(" "); cbdCell->printCell();
+ //printf(" "); cbdCell->printBoundary();
+ cbdCell->addBoundaryCell(ori, cell);
+ coherent = false;
+ }
+
+ }
+
+ }
+ }
+ return coherent;
+}
+
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index e387628f0b34e95b7335c6e72dbb89c638bd9dee..021437788976d355e60784612fc9f889c2bb3b7e 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -52,10 +52,14 @@ class Cell
// cells on the boundary and on the coboundary of thhis cell
std::list< std::pair<int, Cell*> > _boundary;
std::list< std::pair<int, Cell*> > _coboundary;
-
+ int _bdSize;
+ int _cbdSize;
public:
- Cell(){}
+ Cell(){
+ _bdSize = 0;
+ _cbdSize = 0;
+ }
virtual ~Cell(){}
virtual int getDim() const { return _dim; };
@@ -79,8 +83,8 @@ class Cell
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const = 0;
- virtual int getBoundarySize() { return _boundary.size(); }
- virtual int getCoboundarySize() { return _coboundary.size(); }
+ virtual int getBoundarySize() { return _bdSize; }
+ virtual int getCoboundarySize() { return _cbdSize; }
virtual std::list< std::pair<int, Cell*> > getOrientedBoundary() { return _boundary; }
virtual std::list< Cell* > getBoundary() {
@@ -101,43 +105,65 @@ class Cell
return coboundary;
}
- 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;
+ virtual bool addBoundaryCell(int orientation, Cell* cell) {
+ _bdSize++;
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) {
+ (*it).first = (*it).first + orientation;
+ if((*it).first == 0) {
+ _boundary.erase(it); _bdSize--;
+ cell2->removeCoboundaryCell(this,false);
+ return false;
+ }
+ return true;
}
}
- _boundary.push_back( std::make_pair(orientation, cell));
-
+ _boundary.push_back( std::make_pair(orientation, cell ) );
+ return true;
}
- 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;
+
+ virtual bool addCoboundaryCell(int orientation, Cell* cell) {
+ _cbdSize++;
+ for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) {
+ (*it).first = (*it).first + orientation;
+ if((*it).first == 0){
+ _coboundary.erase(it); _cbdSize--;
+ cell2->removeBoundaryCell(this,false);
+ return false;
+ }
+ return true;
}
}
- _coboundary.push_back( std::make_pair(orientation, cell));
+ _coboundary.push_back( std::make_pair(orientation, cell) );
+ return true;
}
- virtual int removeBoundaryCell(Cell* cell) {
- int count = 0;
+ virtual int removeBoundaryCell(Cell* cell, bool other=true) {
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); count++; it = _boundary.begin(); }
+ int ori = (*it).first;
+ if(*cell2 == *cell) {
+ _boundary.erase(it);
+ if(other) cell2->removeCoboundaryCell(this, false);
+ _bdSize--;
+ return ori;
+ }
}
- return count;
}
- virtual int removeCoboundaryCell(Cell* cell) {
- int count = 0;
+ virtual int removeCoboundaryCell(Cell* cell, bool other=true) {
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); count++; it = _coboundary.begin(); }
+ int ori = (*it).first;
+ if(*cell2 == *cell) {
+ _coboundary.erase(it);
+ if(other) cell2->removeBoundaryCell(this, false);
+ _cbdSize--;
+ return ori;
+ }
}
- return count;
}
virtual bool hasBoundary(Cell* cell){
@@ -147,6 +173,14 @@ class Cell
}
return false;
}
+ virtual bool hasCoboundary(Cell* cell){
+ for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) return true;
+ }
+ return false;
+ }
+
virtual void clearBoundary() { _boundary.clear(); }
virtual void clearCoboundary() { _coboundary.clear(); }
@@ -165,12 +199,14 @@ class Cell
Cell* cell2 = (*it).second;
cell2->printCell();
}
+ if(_boundary.empty()) printf("Cell boundary is empty. \n");
}
virtual void printCoboundary() {
for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
- printf("Coboundary cell orientation: %d ", (*it).first);
+ printf("Coboundary cell orientation: %d, ", (*it).first);
Cell* cell2 = (*it).second;
cell2->printCell();
+ if(_coboundary.empty()) printf("Cell coboundary is empty. \n");
}
}
@@ -188,11 +224,12 @@ class Cell
virtual int getNumFacets() const { return 0; }
virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const {};
- virtual bool combined() { return _combined; }
+ virtual bool isCombined() { return _combined; }
virtual std::list< std::pair<int, Cell*> > getCells() { std::list< std::pair<int, Cell*> >cells; cells.push_back( std::make_pair(1, this)); return cells; }
virtual int getNumCells() {return 1;}
bool operator==(const Cell& c2){
+
if(this->getNumVertices() != c2.getNumVertices()){
return false;
}
@@ -202,6 +239,9 @@ class Cell
}
}
return true;
+
+
+
}
bool operator<(const Cell& c2) const {
@@ -261,7 +301,7 @@ class ZeroSimplex : public Simplex, public MPoint
std::vector<MVertex*> getVertexVector() const { std::vector<MVertex*> v; v.push_back(_v[0]); return v; }
std::vector<int> getSortedVertexVector() const { return std::vector<int>(1,_v[0]->getNum()); }
- void printCell() const { printf("Vertices: %d, in subdomain: %d \n", _v[0]->getNum(), _inSubdomain); }
+ void printCell() const { printf("Cell dimension: %d, Vertices: %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _inSubdomain); }
};
@@ -306,7 +346,7 @@ class OneSimplex : public Simplex, public MLine
//int kappa(Cell* tau) const;
- void printCell() const { printf("Vertices: %d %d, in subdomain: %d \n", _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
+ void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
};
// Two simplex cell type.
@@ -350,7 +390,7 @@ class TwoSimplex : public Simplex, public MTriangle
MTriangle::getEdgeVertices(num, v);
}
- void printCell() const { printf("Vertices: %d %d %d, in subdomain: %d\n", _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
+ void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d, in subdomain: %d\n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
};
// Three simplex cell type.
@@ -391,7 +431,7 @@ class ThreeSimplex : public Simplex, public MTetrahedron
MTetrahedron::getFaceVertices(num, v);
}
- virtual void printCell() const { printf("Vertices: %d %d %d %d, in subdomain: %d \n", _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
+ virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
};
// Ordering for the cells.
@@ -454,6 +494,14 @@ class CombinedCell : public Cell{
public:
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
+
+ // use "smaller" cell as c2
+ if(c1->getNumVertices() < c2->getNumVertices()){
+ Cell* temp = c1;
+ c1 = c2;
+ c2 = temp;
+ }
+
_index = c1->getIndex();
_tag = c1->getTag();
_dim = c1->getDim();
@@ -481,96 +529,66 @@ class CombinedCell : public Cell{
_cells.push_back(*it);
}
- _boundary = c1->getOrientedBoundary();
+
+ //_boundary = c1->getOrientedBoundary();
std::list< std::pair<int, Cell*> > c1Boundary = c1->getOrientedBoundary();
std::list< std::pair<int, Cell*> > c2Boundary = c2->getOrientedBoundary();
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
- Cell* cell = (*it).second;
- if(co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = c2Boundary.begin(); it2 != c2Boundary.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cell) old = true;
- }
- if(!old) _boundary.push_back(*it);
- }
- else _boundary.push_back(*it);
- }
- */
-
for(std::list< std::pair<int, Cell*> >::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
Cell* cell = (*it).second;
int ori = (*it).first;
- cell->removeCoboundaryCell(c1);
- cell->addCoboundaryCell(ori, this, true);
+ cell->removeCoboundaryCell(c1);
+ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
}
-
for(std::list< std::pair<int, Cell*> >::iterator it = c2Boundary.begin(); it != c2Boundary.end(); it++){
- int ori2 = (*it).first;
- if(!orMatch) (*it).first = -1*(*it).first;
Cell* cell = (*it).second;
+ if(!orMatch) (*it).first = -1*(*it).first;
int ori = (*it).first;
cell->removeCoboundaryCell(c2);
+ //if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
if(co){
bool old = false;
for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Boundary.begin(); it2 != c1Boundary.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *cell) old = true;
}
- if(!old){ _boundary.push_back(*it); cell->addCoboundaryCell(ori, this, true); }
+ if(!old){ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this); }
+ }
+ else{
+ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
}
- else { _boundary.push_back(*it); cell->addCoboundaryCell(ori, this, true); }
}
- _coboundary = c1->getOrientedCoboundary();
+ //_coboundary = c1->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > c1Coboundary = c1->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > c2Coboundary = c2->getOrientedCoboundary();
- /*
- for(std::list< std::pair<int, Cell* > >::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
-
- Cell* cell = (*it).second;
- if(!co){
- bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = c2Coboundary.begin(); it2 != c2Coboundary.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cell) old = true;
- }
- if(!old) _coboundary.push_back(*it);
- }
- else _coboundary.push_back(*it);
- }
- */
-
for(std::list< std::pair<int, Cell*> >::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
Cell* cell = (*it).second;
int ori = (*it).first;
- cell->removeBoundaryCell(c1);
- cell->addBoundaryCell(ori, this,true);
+ cell->removeBoundaryCell(c1);
+ if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
}
-
-
- for(std::list< std::pair<int, Cell* > >::iterator it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
- int ori2 = (*it).first;
- if(!orMatch) (*it).first = -1*(*it).first;
+ for(std::list< std::pair<int, Cell*> >::iterator it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
Cell* cell = (*it).second;
+ if(!orMatch) (*it).first = -1*(*it).first;
int ori = (*it).first;
cell->removeBoundaryCell(c2);
+ //if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
if(!co){
bool old = false;
for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Coboundary.begin(); it2 != c1Coboundary.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *cell) old = true;
}
- if(!old) { _coboundary.push_back(*it); cell->addBoundaryCell(ori, this, true); }
-
+ if(!old) { if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this); }
+ }
+ else {
+ if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
}
- else { _coboundary.push_back(*it); cell->addBoundaryCell(ori, this, true); }
}
-
+
}
~CombinedCell(){}
@@ -595,9 +613,10 @@ class CombinedCell : public Cell{
}
virtual void printCell() const {
+ printf("Cell dimension: %d, ", getDim() );
printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){
- printf("%d ", this->getVertex(i)->getNum());
+ printf("%d ", this->getSortedVertex(i));
}
printf(", in subdomain: %d\n", _inSubdomain);
}
@@ -626,9 +645,8 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
- // trashbin for cell pointers removed from cell complex
+ // storage for cell pointers to delete them
std::set<Cell*, Less_Cell> _cells2[4];
- std::set<Cell*, Less_Cell> _trash;
//std::set<Cell*, Less_Cell> _originalCells[4];
@@ -654,13 +672,13 @@ class CellComplex
// 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);
-
+
+ public:
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
int reduction(int dim);
-
+ private:
// queued coreduction presented in Mrozek's paper
int coreduction(Cell* generator);
@@ -701,16 +719,11 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
CellComplex(){}
~CellComplex(){
- //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++){
for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
Cell* cell = *cit;
- if(cell->combined()) delete cell;
+ if(cell->isCombined()) delete cell;
}
_cells[i].clear();
for(citer cit = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
@@ -734,6 +747,14 @@ class CellComplex
citer firstCell(int dim) {return _cells[dim].begin(); }
citer lastCell(int dim) {return _cells[dim].end(); }
+ // true if cell complex has given cell
+ bool hasCell(Cell* cell){
+ citer cit = _cells[cell->getDim()].find(cell);
+ if( cit == lastCell(cell->getDim()) ) return false;
+ else return true;
+ }
+
+
// 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); }
@@ -762,19 +783,20 @@ class CellComplex
// write this cell complex in 2.0 MSH ASCII format
int writeComplexMSH(const std::string &name);
+ // Cell combining for reduction and coreduction
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; }
+ // check whether all boundary cells of a cell has this cell as coboundary cell and vice versa
+ // also check whether all (co)boundary cells of a cell belong to this cell complex
+ bool checkCoherence();
+
+ // change roles of boundaries and coboundaries of the cells in this cell complex
+ // equivalent to transposing boundary operator matrices
void makeDualComplex(){
std::set<Cell*, Less_Cell> temp = _cells[0];
_cells[0] = _cells[3];
@@ -791,26 +813,8 @@ 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!");
- }
-
- }
- }
- }
+
+
};
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 22e44eb8342d09814a528e56877fc4a0badc4bb5..56091c0688b0328cbe06be61d675d6bec9de4919 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -86,8 +86,9 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
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);
+ printf("Warning: Invalid incidence index: %d! HMatrix: %d.", (old_elem + (*it).first), dim);
+ printf(" Set to %d. \n", (old_elem + (*it).first) % 2);
+ mpz_set_si(elem, (old_elem + (*it).first) % 2);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
@@ -351,6 +352,7 @@ void ChainComplex::computeHomology(bool dual){
gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
+
// 5) General case:
// 1) Find the bases of boundaries B and cycles Z
// 2) find j: B -> Z and
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index f6edfd5a40c4ba246dfc193461c1722fcc4f9ae1..39b7eba797efbc2a5017d84590e77e7961fdaf6b 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -78,23 +78,22 @@ void Homology::findGenerators(std::string fileName){
double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(true);
- //_cellComplex->checkCoherence();
if(getCombine()){
_cellComplex->combine(3);
- _cellComplex->reduceComplex(false);
+ _cellComplex->reduction(2);
_cellComplex->combine(2);
- _cellComplex->reduceComplex(false);
+ _cellComplex->reduction(1);
_cellComplex->combine(1);
- _cellComplex->reduceComplex(false);
}
- //_cellComplex->checkCoherence();
+ _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));
- //for(int i = 0; i < 4; i++) _cellComplex->printComplex(i);
-
+ //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
+
_cellComplex->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
@@ -146,16 +145,20 @@ void Homology::findGenerators(std::string fileName){
void Homology::findThickCuts(std::string fileName){
+ //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
+
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(true);
+ //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
+ //_cellComplex->coreduceComplex(false);
//_cellComplex->checkCoherence();
if(getCombine()){
_cellComplex->cocombine(0);
_cellComplex->cocombine(1);
_cellComplex->cocombine(2);
}
- //_cellComplex->checkCoherence();
+ _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.",