diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 9c004e0e0eee30b84b658ecdf7758e09d12a5976..cb8818c06f6735e92453f889b50c7e098daa911a 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -106,7 +106,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
//cell->setTag(tag);
tag++;
}
- //_cells2[i] = _cells[i];
+ _ocells[i] = _cells[i];
_betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
}
@@ -116,7 +116,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!
+ // works only for simplicial meshes at the moment!
std::vector<GEntity*> domain;
@@ -200,13 +200,17 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
if(!subdomain && !boundary){
int ori = cell->kappa(oldCell);
oldCell->addCoboundaryCell( ori, cell );
+ oldCell->addOrgCbdCell( ori, cell );
cell->addBoundaryCell( ori, oldCell);
+ cell->addOrgBdCell( ori, oldCell);
}
}
else if(!subdomain && !boundary) {
int ori = cell->kappa(newCell);
cell->addBoundaryCell( ori, newCell );
+ cell->addOrgBdCell( ori, newCell );
newCell->addCoboundaryCell( ori, cell);
+ newCell->addOrgCbdCell( ori, cell);
}
}
}
@@ -215,7 +219,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
void CellComplex::insertCell(Cell* cell){
- _cells[cell->getDim()].insert(cell);
+ _ecells.insert(cell);
}
/*
@@ -290,22 +294,22 @@ int OneSimplex::kappa(Cell* tau) const{
void CellComplex::removeCell(Cell* cell, bool other){
- std::list< std::pair< int, Cell*> > coboundary = cell->getOrientedCoboundary();
- std::list< std::pair< int, Cell*> > boundary = cell->getOrientedBoundary();
+ std::map<Cell*, int, Less_Cell > coboundary = cell->getOrientedCoboundary();
+ std::map<Cell*, int, Less_Cell > boundary = cell->getOrientedBoundary();
//std::list<Cell*> boundary = cell->getBoundary();
//std::list<Cell*> coboundary = cell->getCoboundary();
- for(std::list< std::pair< int, Cell*> >::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ for(std::map<Cell*, int, Less_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).first;
//Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell, other);
}
- for(std::list< std::pair< int, Cell*> >::iterator it = boundary.begin(); it != boundary.end(); it++){
+ for(std::map<Cell*, int, Less_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).first;
//Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell, other);
}
@@ -673,7 +677,7 @@ int CellComplex::cocombine(int dim){
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::list<Cell*> cbd_c;
- std::list< std::pair< int, Cell*> > bd_c;
+ std::map<Cell*, int, Less_Cell > bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
@@ -687,32 +691,36 @@ int CellComplex::cocombine(int dim){
Q.pop();
bd_c = s->getOrientedBoundary();
+ if(s->getBoundarySize() == 2){
- 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)){
+ std::map<Cell*, int, Less_Cell>::iterator it = bd_c.begin();
+ int or1 = (*it).second;
+ Cell* c1 = (*it).first;
+ it++;
+ int or2 = (*it).second;
+ Cell* c2 = (*it).first;
- int or1 = bd_c.front().first;
- int or2 = bd_c.back().first;
- Cell* c1 = bd_c.front().second;
- Cell* c2 = bd_c.back().second;
-
- removeCell(s);
- _trash.push_back(s);
-
- cbd_c = c1->getCoboundary();
- enqueueCells(cbd_c, Q, Qset);
- cbd_c = c2->getCoboundary();
- enqueueCells(cbd_c, Q, Qset);
+ if(!(*c1 == *c2)
+ && inSameDomain(s, c1) && inSameDomain(s, c2)
+ && c1->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
+ && c2->getNumVertices() < getSize(dim)){
+
+ removeCell(s);
+ _trash.push_back(s);
- CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
- removeCell(c1);
- removeCell(c2);
- _cells[dim].insert(newCell);
+ cbd_c = c1->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+ cbd_c = c2->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
- cit = firstCell(dim);
- count++;
+ CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
+ removeCell(c1);
+ removeCell(c2);
+ _cells[dim].insert(newCell);
+
+ cit = firstCell(dim);
+ count++;
+ }
}
removeCellQset(s, Qset);
@@ -735,7 +743,7 @@ int CellComplex::combine(int dim){
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
- std::list< std::pair<int, Cell*> > cbd_c;
+ std::map<Cell*, int, Less_Cell> cbd_c;
std::list<Cell*> bd_c;
int count = 0;
@@ -748,33 +756,37 @@ int CellComplex::combine(int dim){
Cell* s = Q.front();
Q.pop();
-
cbd_c = s->getOrientedCoboundary();
-
- 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(s->getCoboundarySize() == 2){
+ std::map<Cell*, int, Less_Cell>::iterator it = cbd_c.begin();
+ int or1 = (*it).second;
+ Cell* c1 = (*it).first;
+ it++;
+ int or2 = (*it).second;
+ Cell* c2 = (*it).first;
+
+ if(!(*c1 == *c2)
+ && inSameDomain(s, c1) && inSameDomain(s, c2)
+ && c1->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
+ && c2->getNumVertices() < getSize(dim)){
+
+ removeCell(s);
+ _trash.push_back(s);
- removeCell(s);
- _trash.push_back(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) );
- 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);
- count++;
+ 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) );
+ 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);
+ count++;
+ }
}
removeCellQset(s, Qset);
@@ -905,7 +917,8 @@ int CellComplex::writeComplexMSH(const std::string &name){
for(int i = 0; i < _store.size(); i++){
count = count + _store.at(i).size();
}
-
+ count = count + _ecells.size();
+
fprintf(fp, "%d\n", count);
int partition = 0;
@@ -979,6 +992,26 @@ int CellComplex::writeComplexMSH(const std::string &name){
}
}
+ for(citer cit = _ecells.begin(); cit != _ecells.end(); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()) partition = 3;
+ else if(cell->onDomainBoundary()) partition = 2;
+ else partition = 1;
+ if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), 15, 3, 0, 0, partition, cell->getVertex(0)->getNum());
+ if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+ if(cell->getDim() == 2){
+ if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 3, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ }
+ if(cell->getDim() == 3){
+ if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 12, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
+ if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 13, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
+ if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 14, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
+ }
+ }
+
+
fprintf(fp, "$EndElements\n");
@@ -1003,10 +1036,10 @@ bool CellComplex::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 = 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;
+ std::map<Cell*, int, Less_Cell> boundary = cell->getOrientedBoundary();
+ for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){
+ Cell* bdCell = (*it).first;
+ int ori = (*it).second;
citer cit = _cells[bdCell->getDim()].find(bdCell);
if(cit == lastCell(bdCell->getDim())){
printf("Warning! Boundary cell not in cell complex! Boundary removed. \n");
@@ -1026,10 +1059,10 @@ bool CellComplex::checkCoherence(){
}
}
- 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;
+ std::map<Cell*, int, Less_Cell> coboundary = cell->getOrientedCoboundary();
+ for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ Cell* cbdCell = (*it).first;
+ int ori = (*it).second;
citer cit = _cells[cbdCell->getDim()].find(cbdCell);
if(cit == lastCell(cbdCell->getDim())){
printf("Warning! Coboundary cell not in cell complex! Coboundary removed. \n");
@@ -1055,4 +1088,24 @@ bool CellComplex::checkCoherence(){
return coherent;
}
+
+bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const {
+
+ //cells with fever vertices first
+
+ if(c1->getNumVertices() != c2->getNumVertices()){
+ return (c1->getNumVertices() < c2->getNumVertices());
+ }
+
+
+
+ //"natural number" -like ordering (the number of a vertice corresponds a digit)
+
+ for(int i=0; i < c1->getNumVertices();i++){
+ if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
+ else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
+ }
+
+ return false;
+}
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 8d561d861c2843898e51b2b0a2086dd2ad721706..8a647f1c85290c8eff91c868f8a79c2f190920bd 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -32,6 +32,14 @@
#include "GFace.h"
#include "GVertex.h"
+class Cell;
+
+// Ordering of the cells
+class Less_Cell{
+ public:
+ bool operator()(const Cell* c1, const Cell* c2) const;
+};
+
// Abstract class representing an elemtary cell of a cell complex.
class Cell
{
@@ -56,12 +64,16 @@ class Cell
bool _immune;
- // 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;
+ // cells on the boundary and on the coboundary of this cell
+ std::map< Cell*, int, Less_Cell > _boundary;
+ std::map< Cell*, int, Less_Cell > _coboundary;
int _bdSize;
int _cbdSize;
+ // original boundary and coboundary before the reduction of the cell complex
+ std::map<Cell*, int, Less_Cell > _obd;
+ std::map<Cell*, int, Less_Cell > _ocbd;
+
public:
Cell(){
_bdSize = 0;
@@ -94,24 +106,27 @@ class Cell
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const = 0;
-
+
+ // (co)boundary cell iterator
+ typedef std::map<Cell*, int, Less_Cell>::iterator biter;
+
virtual int getBoundarySize() { return _bdSize; }
virtual int getCoboundarySize() { return _cbdSize; }
- virtual std::list< std::pair<int, Cell*> > getOrientedBoundary() { return _boundary; }
+ virtual std::map<Cell*, int, Less_Cell > getOrientedBoundary() { return _boundary; }
virtual std::list< Cell* > getBoundary() {
std::list<Cell*> boundary;
- for( std::list< std::pair<int, Cell*> >::iterator it= _boundary.begin();it!= _boundary.end();it++){
- Cell* cell = (*it).second;
+ for( biter it= _boundary.begin(); it != _boundary.end();it++){
+ Cell* cell = (*it).first;
boundary.push_back(cell);
}
return boundary;
}
- virtual std::list<std::pair<int, Cell*> >getOrientedCoboundary() { return _coboundary; }
+ virtual std::map<Cell*, int, Less_Cell > getOrientedCoboundary() { return _coboundary; }
virtual std::list< Cell* > getCoboundary() {
std::list<Cell*> coboundary;
- for( std::list< std::pair<int, Cell*> >::iterator it= _coboundary.begin();it!= _coboundary.end();it++){
- Cell* cell = (*it).second;
+ for( biter it= _coboundary.begin();it!= _coboundary.end();it++){
+ Cell* cell = (*it).first;
coboundary.push_back(cell);
}
return coboundary;
@@ -119,6 +134,18 @@ class Cell
virtual bool addBoundaryCell(int orientation, Cell* cell) {
_bdSize++;
+ biter it = _boundary.find(cell);
+ if(it != _boundary.end()){
+ (*it).second = (*it).second + orientation;
+ if((*it).second == 0) {
+ _boundary.erase(it); _bdSize--;
+ (*it).first->removeCoboundaryCell(this,false);
+ return false;
+ }
+ return true;
+ }
+ _boundary.insert( std::make_pair(cell, orientation ) );
+ /*
for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = (*it).second;
if(*cell2 == *cell) {
@@ -132,28 +159,36 @@ class Cell
}
}
_boundary.push_back( std::make_pair(orientation, cell ) );
+ */
return true;
}
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;
+ biter it = _coboundary.find(cell);
+ if(it != _coboundary.end()){
+ (*it).second = (*it).second + orientation;
+ if((*it).second == 0) {
+ _coboundary.erase(it); _cbdSize--;
+ (*it).first->removeBoundaryCell(this,false);
+ return false;
}
+ return true;
}
- _coboundary.push_back( std::make_pair(orientation, cell) );
+ _coboundary.insert( std::make_pair(cell, orientation ) );
return true;
}
virtual int removeBoundaryCell(Cell* cell, bool other=true) {
+ biter it = _boundary.find(cell);
+ if(it != _boundary.end()){
+ _boundary.erase(it);
+ if(other) (*it).first->removeCoboundaryCell(this, false);
+ _bdSize--;
+ return (*it).second;
+ }
+
+ /*
for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = (*it).second;
int ori = (*it).first;
@@ -164,35 +199,51 @@ class Cell
return ori;
}
}
+ */
+
return 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;
- int ori = (*it).first;
- if(*cell2 == *cell) {
- _coboundary.erase(it);
- if(other) cell2->removeBoundaryCell(this, false);
- _cbdSize--;
- return ori;
- }
+ biter it = _coboundary.find(cell);
+ if(it != _coboundary.end()){
+ _coboundary.erase(it);
+ if(other) (*it).first->removeBoundaryCell(this, false);
+ _cbdSize--;
+ return (*it).second;
}
return 0;
}
- virtual bool hasBoundary(Cell* cell){
- for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- Cell* cell2 = (*it).second;
- if(*cell2 == *cell) return true;
+ virtual bool hasBoundary(Cell* cell, bool org=false){
+ if(!org){
+ /*
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) return true;
+ }
+ return false;
+ */
+ biter it = _boundary.find(cell);
+ if(it != _boundary.end()) return true;
+ return false;
+ }
+ else{
+ biter it = _obd.find(cell);
+ if(it != _obd.end()) return true;
+ return false;
}
- 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;
+ virtual bool hasCoboundary(Cell* cell, bool org=false){
+ if(!org){
+ biter it = _coboundary.find(cell);
+ if(it != _coboundary.end()) return true;
+ return false;
}
- return false;
+ else{
+ biter it = _ocbd.find(cell);
+ if(it != _ocbd.end()) return true;
+ return false;
+ }
}
@@ -200,7 +251,7 @@ class Cell
virtual void clearCoboundary() { _coboundary.clear(); }
virtual void makeDualCell(){
- std::list< std::pair<int, Cell*> > temp = _boundary;
+ std::map<Cell*, int, Less_Cell > temp = _boundary;
_boundary = _coboundary;
_coboundary = temp;
_dim = 3-_dim;
@@ -208,22 +259,44 @@ class Cell
}
virtual void printBoundary() {
- for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- printf("Boundary cell orientation: %d ", (*it).first);
- Cell* cell2 = (*it).second;
+ for(biter it = _boundary.begin(); it != _boundary.end(); it++){
+ printf("Boundary cell orientation: %d ", (*it).second);
+ Cell* cell2 = (*it).first;
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);
- Cell* cell2 = (*it).second;
+ for(biter it = _coboundary.begin(); it != _coboundary.end(); it++){
+ printf("Coboundary cell orientation: %d, ", (*it).second);
+ Cell* cell2 = (*it).first;
cell2->printCell();
if(_coboundary.empty()) printf("Cell coboundary is empty. \n");
}
}
+ virtual void addOrgBdCell(int orientation, Cell* cell) { _obd.insert( std::make_pair(cell, orientation ) ); };
+ virtual void addOrgCbdCell(int orientation, Cell* cell) { _ocbd.insert( std::make_pair(cell, orientation ) ); };
+ virtual std::map<Cell*, int, Less_Cell > getOrgBd() { return _obd; }
+ virtual std::map<Cell*, int, Less_Cell > getOrgCbd() { return _ocbd; }
+ virtual void printOrgBd() {
+ for(biter it = _obd.begin(); it != _obd.end(); it++){
+ printf("Boundary cell orientation: %d ", (*it).second);
+ Cell* cell2 = (*it).first;
+ cell2->printCell();
+ }
+ if(_obd.empty()) printf("Cell boundary is empty. \n");
+ }
+ virtual void printOrgCbd() {
+ for(biter it = _ocbd.begin(); it != _ocbd.end(); it++){
+ printf("Coboundary cell orientation: %d, ", (*it).second);
+ Cell* cell2 = (*it).first;
+ cell2->printCell();
+ if(_ocbd.empty()) printf("Cell coboundary is empty. \n");
+ }
+ }
+
+
virtual bool inSubdomain() const { return _inSubdomain; }
@@ -632,8 +705,7 @@ class CPyramid : public Cell, public MPyramid
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _inSubdomain); }
};
-
-
+/*
// Ordering for the cells.
class Less_Cell{
public:
@@ -661,6 +733,7 @@ class Less_Cell{
}
};
+*/
class Equal_Cell{
@@ -692,6 +765,7 @@ class Less_MVertex{
}
};
+// A cell that is a combination of cells of same dimension
class CombinedCell : public Cell{
private:
@@ -740,28 +814,32 @@ class CombinedCell : public Cell{
//_boundary = c1->getOrientedBoundary();
- std::list< std::pair<int, Cell*> > c1Boundary = c1->getOrientedBoundary();
- std::list< std::pair<int, Cell*> > c2Boundary = c2->getOrientedBoundary();
+ std::map< Cell*, int, Less_Cell > c1Boundary = c1->getOrientedBoundary();
+ std::map< Cell*, int, Less_Cell > c2Boundary = c2->getOrientedBoundary();
- for(std::list< std::pair<int, Cell*> >::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
- Cell* cell = (*it).second;
- int ori = (*it).first;
+ for(std::map<Cell*, int, Less_Cell>::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
+ Cell* cell = (*it).first;
+ int ori = (*it).second;
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++){
- Cell* cell = (*it).second;
- if(!orMatch) (*it).first = -1*(*it).first;
- int ori = (*it).first;
+ for(std::map<Cell*, int, Less_Cell >::iterator it = c2Boundary.begin(); it != c2Boundary.end(); it++){
+ Cell* cell = (*it).first;
+ if(!orMatch) (*it).second = -1*(*it).second;
+ int ori = (*it).second;
cell->removeCoboundaryCell(c2);
//if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
if(co){
+ std::map<Cell*, int, Less_Cell >::iterator it2 = c1Boundary.find(cell);
bool old = false;
+ if(it2 != c1Boundary.end()) old = true;
+ /*
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){ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this); }
}
else{
@@ -770,27 +848,31 @@ class CombinedCell : public Cell{
}
//_coboundary = c1->getOrientedCoboundary();
- std::list< std::pair<int, Cell*> > c1Coboundary = c1->getOrientedCoboundary();
- std::list< std::pair<int, Cell*> > c2Coboundary = c2->getOrientedCoboundary();
+ std::map<Cell*, int, Less_Cell > c1Coboundary = c1->getOrientedCoboundary();
+ std::map<Cell*, int, Less_Cell > c2Coboundary = c2->getOrientedCoboundary();
- for(std::list< std::pair<int, Cell*> >::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
- Cell* cell = (*it).second;
- int ori = (*it).first;
+ for(std::map<Cell*, int, Less_Cell>::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
+ Cell* cell = (*it).first;
+ int ori = (*it).second;
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++){
- Cell* cell = (*it).second;
- if(!orMatch) (*it).first = -1*(*it).first;
- int ori = (*it).first;
+ for(std::map<Cell*, int, Less_Cell>::iterator it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
+ Cell* cell = (*it).first;
+ if(!orMatch) (*it).second = -1*(*it).second;
+ int ori = (*it).second;
cell->removeBoundaryCell(c2);
//if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
if(!co){
+ std::map<Cell*, int, Less_Cell >::iterator it2 = c1Coboundary.find(cell);
bool old = false;
+ if(it2 != c1Coboundary.end()) old = true;
+ /*
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) { if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this); }
}
else {
@@ -865,8 +947,10 @@ class CellComplex
std::list<Cell*> _trash;
std::vector< std::set<Cell*, Less_Cell> > _store;
+ std::set<Cell*, Less_Cell> _ecells;
+
- //std::set<Cell*, Less_Cell> _originalCells[4];
+ std::set<Cell*, Less_Cell> _ocells[4];
// Betti numbers of this cell complex (ranks of homology groups)
int _betti[4];
@@ -890,9 +974,9 @@ class CellComplex
// remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
+ public:
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, int omitted=0);
@@ -938,22 +1022,23 @@ class CellComplex
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 = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
+
+ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
Cell* cell = *cit;
delete cell;
}
- _cells2[i].clear();
- */
+ _ocells[i].clear();
+
}
- emptyTrash();
+ //emptyTrash();
}
void emptyTrash(){
@@ -966,20 +1051,32 @@ class CellComplex
// get the number of certain dimensional cells
int getSize(int dim){ return _cells[dim].size(); }
+ int getOrgSize(int dim){ return _ocells[dim].size(); }
int getDim() {return _dim; }
std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
-
+ std::set<Cell*, Less_Cell> getOrgCells(int dim){ return _ocells[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 firstOrgCell(int dim) {return _ocells[dim].begin(); }
+ citer lastOrgCell(int dim) {return _ocells[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;
+ bool hasCell(Cell* cell, bool org=false){
+ if(!org){
+ citer cit = _cells[cell->getDim()].find(cell);
+ if( cit == lastCell(cell->getDim()) ) return false;
+ else return true;
+ }
+ else{
+ citer cit = _ocells[cell->getDim()].find(cell);
+ if( cit == lastOrgCell(cell->getDim()) ) return false;
+ else return true;
+ }
}
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 01245f710c510a5b4ceff3391aaafb74316d5815..58a00afadb7bbc4a79a5aeb8fb9146871a90585d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -6,6 +6,7 @@
// Contributed by Matti Pellikka.
#include "ChainComplex.h"
+#include "OS.h"
#if defined(HAVE_KBIPACK)
@@ -71,9 +72,9 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if(!cell->inSubdomain()){
- std::list< std::pair<int,Cell*> >bdCell = cell->getOrientedBoundary();
- for(std::list< std::pair<int, Cell*> >::iterator it = bdCell.begin(); it != bdCell.end(); it++){
- Cell* bdCell = (*it).second;
+ std::map<Cell*, int, Less_Cell> bdCell = cell->getOrientedBoundary();
+ for(std::map<Cell*, int, Less_Cell>::iterator it = bdCell.begin(); it != bdCell.end(); it++){
+ Cell* bdCell = (*it).first;
if(!bdCell->inSubdomain()){
int old_elem = 0;
//printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex());
@@ -84,11 +85,11 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
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.", (old_elem + (*it).first), dim);
- printf(" Set to %d. \n", (old_elem + (*it).first) % 2);
- mpz_set_si(elem, (old_elem + (*it).first) % 2);
+ mpz_set_si(elem, old_elem + (*it).second);
+ if( (old_elem + (*it).second) > 1 || (old_elem + (*it).second) < -1 ){
+ printf("Warning: Invalid incidence index: %d! HMatrix: %d.", (old_elem + (*it).second), dim);
+ printf(" Set to %d. \n", (old_elem + (*it).second) % 2);
+ mpz_set_si(elem, (old_elem + (*it).second) % 2);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
@@ -459,8 +460,17 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
int i = 0;
for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
Cell* cell = *cit;
+ _dim = cell->getDim();
if(!cell->inSubdomain() && (int)coeffs.size() > i){
- if(coeffs.at(i) != 0) _cells.push_back( std::make_pair(cell, coeffs.at(i)) );
+ if(coeffs.at(i) != 0){
+ std::list< std::pair<int, Cell*> > subCells = cell->getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = subCells.begin(); it != subCells.end(); it++){
+ Cell* subCell = (*it).second;
+ int coeff = (*it).first;
+ _cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
+ }
+ //_cells.push_back( std::make_pair(cell, coeffs.at(i)) );
+ }
i++;
}
@@ -471,8 +481,193 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
}
-int Chain::writeChainMSH(const std::string &name){
+Cell* Chain::findCommonCbdCell(Cell* c1, Cell* c2, Cell* c3){
+ /*
+ c1->printCell();
+ c1->printOrgCbd();
+ c2->printOrgCbd();
+ c2->printCell();
+ printf("------ \n");
+ */
+
+ std::map< Cell*, int, Less_Cell > c1Cbd = c1->getOrgCbd();
+ for(std::map< Cell*, int, Less_Cell >::iterator it = c1Cbd.begin(); it != c1Cbd.end(); it++){
+ Cell* c1CbdCell = (*it).first;
+ if(c3 == NULL){
+ if(c2->hasCoboundary(c1CbdCell, true)) return c1CbdCell;
+ }
+ else{
+ if(c2->hasCoboundary(c1CbdCell, true) && c3->hasCoboundary(c1CbdCell, true)) return c1CbdCell;
+ }
+ }
+
+ return NULL;
+
+}
+
+std::pair<Cell*, int> Chain::findRemainingBoundary(Cell* b, Cell* c1, Cell* c2, Cell* c3){
+ std::map<Cell*, int, Less_Cell> bBd = b->getOrgBd();
+ for(std::map<Cell*, int, Less_Cell >::iterator it = bBd.begin(); it != bBd.end(); it++){
+ Cell* bBdCell = (*it).first;
+ if(c3 == NULL) {
+ if( !(*bBdCell == *c1) && !(*bBdCell == *c2)) return *it;
+ }
+ else{
+ if( !(*bBdCell == *c1) && !(*bBdCell == *c2) && !(*bBdCell == *c3) ) return *it;
+ }
+ }
+
+ return std::make_pair(b, 0);
+}
+
+int Chain::findOrientation(Cell* b, Cell* c){
+ std::map< Cell*, int, Less_Cell > bBd = b->getOrgBd();
+ std::map< Cell*, int, Less_Cell >::iterator it = bBd.find(c);
+ if(it != bBd.end()) return (*it).second;
+ return 0;
+}
+
+void Chain::smoothenChain(){
+
+ int start = getSize();
+ double t1 = Cpu();
+
+ 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;
+
+ }
+ }
+ }
+ }
+ }
+ }
+ /*
+ else if(getDim() == 2){
+ for(citer i = _cells.begin(); i != _cells.end(); i++){
+ Cell* c1 = (*i).first;
+ int c1c = (*i).second;
+
+ /*
+ int c1o = 0;
+
+ Cell* c2 = NULL;
+ int c2c = 0;
+ int c2o = 0;
+
+ Cell* c3 = NULL;
+ int c3c = 0;
+ int c3o = 0;
+
+ 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;
+ }
+ }
+
+ }
+
+
+
+ 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;
+ }
+ }
+ }
+ }
+ }
+
+ }
+ }*/
+
+ eraseNullCells();
+ double t2 = Cpu();
+ 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);
if(getSize() == 0) return 0;
@@ -483,7 +678,7 @@ int Chain::writeChainMSH(const std::string &name){
printf("Unable to open file.");
return 0;
}
-
+
fprintf(fp, "\n$ElementData\n");
fprintf(fp, "1 \n");
@@ -493,17 +688,21 @@ int Chain::writeChainMSH(const std::string &name){
fprintf(fp, "4 \n");
fprintf(fp, "0 \n");
fprintf(fp, "1 \n");
- fprintf(fp, "%d \n", getNumCells());
+ fprintf(fp, "%d \n", getSize());
fprintf(fp, "0 \n");
- std::list< std::pair<int, Cell*> > cells;
- for(int i = 0; i < getSize(); i++){
- Cell* cell = getCell(i);
- cells = cell->getCells();
+ //std::list< std::pair<int, Cell*> > cells;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ int coeff = (*cit).second;
+ fprintf(fp, "%d %d \n", cell->getNum(), coeff );
+ /*
+ cells = cell->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell2 = (*it).second;
- fprintf(fp, "%d %d \n", cell2->getNum(), getCoeff(i)*(*it).first );
+ fprintf(fp, "%d %d \n", cell2->getNum(), getCoeff(i)*(*it).first );
}
+ */
}
fprintf(fp, "$EndElementData\n");
@@ -514,13 +713,15 @@ int Chain::writeChainMSH(const std::string &name){
}
+/*
void Chain::getData(std::map<int, std::vector<double> > & data){
if(getSize() == 0) return;
- std::list< std::pair<int, Cell*> > cells;
+ //std::list< std::pair<int, Cell*> > cells;
for(int i = 0; i < getSize(); i++){
Cell* cell = getCell(i);
+ /*
cells = cell->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell2 = (*it).second;
@@ -531,6 +732,12 @@ void Chain::getData(std::map<int, std::vector<double> > & data){
//printf("%d, %d, \n", cell2->getNum(), (int)coeff.at(0));
}
+
+ std::vector<double> coeff;
+ coeff.push_back(getCoeff(i));
+ std::pair<int, std::vector<double> > dataPair = std::make_pair(cell->getTag(), coeff);
+ data.insert(dataPair);
+
}
//for(std::map<int, std::vector<double> >::iterator it = data.begin(); it != data.end(); it++){
@@ -539,5 +746,6 @@ void Chain::getData(std::map<int, std::vector<double> > & data){
return;
}
+*/
#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index c54ed874db0810061003e57343653a80f747914d..d4527bbb98445a4691d20a5663464d398bee1ec6 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -133,7 +133,7 @@ class Chain{
private:
// cells and their coefficients in this chain
- std::vector< std::pair <Cell*, int> > _cells;
+ std::map< Cell*, int, Less_Cell > _cells;
// name of the chain (optional)
std::string _name;
// cell complex this chain belongs to
@@ -141,28 +141,137 @@ class Chain{
int _torsion;
+ int _dim;
+
+ 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);
+
public:
Chain(){}
Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain", int torsion=0);
+ Chain(Chain* chain){
+ _cells = chain->getCells();
+ _torsion = chain->getTorsion();
+ _name = chain->getName();
+ _cellComplex = chain->getCellComplex();
+ _dim = chain->getDim();
+ }
~Chain(){}
+ typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+
// get i:th cell of this chain
- Cell* getCell(int i) { return _cells.at(i).first; }
+ //Cell* getCell(int i) { return _cells.at(i).first; }
// get coeffcient of i:th cell of this chain
- int getCoeff(int i) { return _cells.at(i).second; }
+ //int getCoeff(int i) { return _cells.at(i).second; }
+
+
+
+ void removeCell(Cell* cell, int coeff) {
+ 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);
+ if(!_cellComplex->hasCell(cell)){
+ _cellComplex->insertCell(cell);
+ }
+
+ //}
+ return;
+ }
+ /*
+ void insertCells(){
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ std::vector<Cell*> cells;
+ if(!_cellComplex->hasCell(cell)){
+ _cellComplex->insertCell(cell);
+ }
+
+ }
+ return;
+ }
+ */
+
+ bool hasCell(Cell* c){
+ citer it = _cells.find(c);
+ 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;
+ return NULL;
+ }
+ int getCoeff(Cell* c){
+ citer it = _cells.find(c);
+ if(it != _cells.end()) return (*it).second;
+ return 0;
+ }
+
+
+
int getTorsion() {return _torsion;}
+ int getDim() {return _dim;}
+ CellComplex* getCellComplex() {return _cellComplex;}
+ std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
+
+ void eraseNullCells(){
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if( (*cit).second == 0){
+ _cells.erase(cit);
+ ++cit;
+ }
+ }
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if( (*cit).second == 0){
+ _cells.erase(cit);
+ cit = _cells.begin();
+ }
+ }
+
+
+ return;
+ }
// number of cells in this chain
- int getSize() { return _cells.size(); }
+ int getSize() {
+ //eraseNullCells();
+ return _cells.size();
+ }
int getNumCells() {
- int count = 0;
- for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
- Cell* cell = (*it).first;
- count = count + cell->getNumCells();
- }
- return count;
+ //int count = 0;
+ //for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
+ // Cell* cell = (*it).first;
+ // count = count + cell->getNumCells();
+ //}
+ //return count;
+ return _cells.size();
}
@@ -170,13 +279,16 @@ class Chain{
std::string getName() { return _name; }
void setName(std::string name) { _name=name; }
- // append this chain to a 2.0 MSH ASCII file as $ElementData
- int writeChainMSH(const std::string &name);
+ void smoothenChain();
- void getData(std::map<int, std::vector<double> >& data);
+ // append this chain to a 2.1 MSH ASCII file as $ElementData
+ int writeChainMSH(const std::string &name);
+ //void getData(std::map<int, std::vector<double> >& data);
+
};
-
+
#endif
-
+
#endif
+
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index a81a6b82ca1c012aea9b8dbbe19eb67d62c8a5fe..4a044c4b9938ac18ea1d4f5fdf0e4794a5d7609e 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -84,7 +84,7 @@ void Homology::findGenerators(std::string fileName){
int omitted = _cellComplex->reduceComplex(_omit);
//_cellComplex->printEuler();
- _cellComplex->emptyTrash();
+ //_cellComplex->emptyTrash();
if(getCombine()){
_cellComplex->combine(3);
@@ -96,7 +96,7 @@ void Homology::findGenerators(std::string fileName){
_cellComplex->checkCoherence();
//_cellComplex->printEuler();
- _cellComplex->emptyTrash();
+ //_cellComplex->emptyTrash();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
@@ -105,7 +105,7 @@ void Homology::findGenerators(std::string fileName){
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
- _cellComplex->writeComplexMSH(fileName);
+ //_cellComplex->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
t1 = Cpu();
@@ -113,7 +113,9 @@ void Homology::findGenerators(std::string fileName){
chains->computeHomology();
t2 = Cpu();
Msg::Info("Homology Computation complete (%g s).", t2 - t1);
-
+
+ std::vector<Chain*> chainVector;
+
int HRank[4];
for(int j = 0; j < 4; j++){
HRank[j] = 0;
@@ -126,12 +128,20 @@ 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->writeChainMSH(fileName);
+ Chain* chain2 = new Chain(chain);
+ //printf("chain %d \n", i);
+ t1 = Cpu();
+ int start = chain->getSize();
+ chain->smoothenChain();
+ t2 = Cpu();
+ Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).", j, i, start, chain->getSize(), t2 - t1);
if(chain->getSize() != 0) {
HRank[j] = HRank[j] + 1;
if(chain->getTorsion() != 1) Msg::Warning("H%d %d has torsion coefficient %d!", j, i, chain->getTorsion());
}
- delete chain;
+ chainVector.push_back(chain);
+ chainVector.push_back(chain2);
+ //delete chain;
}
if(j == _cellComplex->getDim() && _cellComplex->getNumOmitted() > 0){
for(int i = 0; i < _cellComplex->getNumOmitted(); i++){
@@ -140,16 +150,25 @@ void Homology::findGenerators(std::string fileName){
std::string name = "H" + dimension + getDomainString() + generator;
std::vector<int> coeffs (_cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, name, 1);
- chain->writeChainMSH(fileName);
if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
- delete chain;
-
+ //delete chain;
+ chainVector.push_back(chain);
}
}
}
+ _cellComplex->writeComplexMSH(fileName);
+ for(int i = 0; i < chainVector.size(); i++){
+ Chain* chain = chainVector.at(i);
+ chain->writeChainMSH(fileName);
+ chainVector.at(i) = NULL;
+ delete chain;
+ }
+ chainVector.clear();
+
+
Msg::Info("Ranks of homology groups for primal cell complex:");
Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]);
@@ -177,7 +196,7 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(_omit);
- _cellComplex->emptyTrash();
+ //_cellComplex->emptyTrash();
/*
_cellComplex->makeDualComplex();
@@ -196,7 +215,7 @@ void Homology::findDualGenerators(std::string fileName){
_cellComplex->cocombine(2);
}
- _cellComplex->emptyTrash();
+ //_cellComplex->emptyTrash();
_cellComplex->checkCoherence();
double t2 = Cpu();
@@ -288,10 +307,10 @@ void Homology::computeBettiNumbers(){
double t2 = Cpu();
Msg::Info("Betti number computation complete (%g s).", t2- t1);
- Msg::Info("H0 = %d \n", _cellComplex->getBettiNumber(0));
- Msg::Info("H1 = %d \n", _cellComplex->getBettiNumber(1));
- Msg::Info("H2 = %d \n", _cellComplex->getBettiNumber(2));
- Msg::Info("H3 = %d \n", _cellComplex->getBettiNumber(3));
+ Msg::Info("H0 = %d", _cellComplex->getBettiNumber(0));
+ Msg::Info("H1 = %d", _cellComplex->getBettiNumber(1));
+ Msg::Info("H2 = %d", _cellComplex->getBettiNumber(2));
+ Msg::Info("H3 = %d", _cellComplex->getBettiNumber(3));
return;
}
diff --git a/Geo/Homology.h b/Geo/Homology.h
index c1e2e7b077dba11e074ce7d1aebeebe9ef11a825..c6e3a28ffe48029a7a0f7f14cb7f69b9fe78eda5 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -30,7 +30,7 @@ class Homology
bool _combine;
int _omit;
-
+
std::vector<int> _domain;
std::vector<int> _subdomain;