diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index fdfb1a6fd44c1d9e0c0f1c42b16647a5e4105ef8..a102d9dbc76ca1bfa371dabc72b73c8b81543a8f 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -18,7 +18,8 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const {
return (c1->getNumVertices() < c2->getNumVertices());
}
- //"natural number" -like ordering (the number of a vertice corresponds a digit)
+ // "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;
@@ -28,8 +29,26 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const {
return false;
}
+
+Cell::Cell(MElement* image, bool subdomain, bool boundary) :
+ _combined(false), _index(0), _immune(false), _image(NULL),
+ _deleteImage(false) {
+ _onDomainBoundary = boundary;
+ _inSubdomain = subdomain;
+ _dim = image->getDim();
+ _image = image;
+ for(int i = 0; i < image->getNumVertices(); i++)
+ _vs.push_back(image->getVertex(i)->getNum());
+ std::sort(_vs.begin(), _vs.end());
+}
+
+Cell::~Cell() {
+ if(_deleteImage) delete _image;
+}
+
bool Cell::hasVertex(int vertex) const {
- std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), vertex);
+ std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(),
+ vertex);
if (it != _vs.end()) return true;
else return false;
}
@@ -96,25 +115,6 @@ void Cell::restoreCell(){
_immune = false;
}
-std::list< Cell* > Cell::getBoundary() {
- std::list<Cell*> boundary;
- for( biter it= _boundary.begin(); it != _boundary.end();it++){
- Cell* cell = (*it).first;
- boundary.push_back(cell);
- }
- return boundary;
-}
-
-std::list< Cell* > Cell::getCoboundary() {
- std::list<Cell*> coboundary;
- for( biter it= _coboundary.begin();it!= _coboundary.end();it++){
- Cell* cell = (*it).first;
- coboundary.push_back(cell);
- }
- return coboundary;
-}
-
-
bool Cell::addBoundaryCell(int orientation, Cell* cell) {
biter it = _boundary.find(cell);
if(it != _boundary.end()){
@@ -256,32 +256,42 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
for(int i = 0; i < c1->getNumVertices(); i++) _v.push_back(c1->getVertex(i));
for(int i = 0; i < c2->getNumVertices(); i++){
- if(!this->hasVertex(c2->getVertex(i)->getNum())) _v.push_back(c2->getVertex(i));
+ if(!this->hasVertex(c2->getVertex(i)->getNum())){
+ _v.push_back(c2->getVertex(i));
+ }
}
// sorted vertices
- for(unsigned int i = 0; i < _v.size(); i++) _vs.push_back(_v.at(i)->getNum());
+ for(unsigned int i = 0; i < _v.size(); i++){
+ _vs.push_back(_v.at(i)->getNum());
+ }
std::sort(_vs.begin(), _vs.end());
// cells
- _cells = c1->getCells();
- std::list< std::pair<int, Cell*> > c2Cells = c2->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = c2Cells.begin(); it != c2Cells.end(); it++){
+ c1->getCells(_cells);
+ std::list< std::pair<int, Cell*> > c2Cells;
+ c2->getCells(c2Cells);
+ for(std::list< std::pair<int, Cell*> >::iterator it = c2Cells.begin();
+ it != c2Cells.end(); it++){
if(!orMatch) (*it).first = -1*(*it).first;
_cells.push_back(*it);
}
// boundary cells
- std::map< Cell*, int, Less_Cell > c1Boundary = c1->getOrientedBoundary();
- std::map< Cell*, int, Less_Cell > c2Boundary = c2->getOrientedBoundary();
+ std::map< Cell*, int, Less_Cell > c1Boundary;
+ c1->getBoundary(c1Boundary);
+ std::map< Cell*, int, Less_Cell > c2Boundary;
+ c2->getBoundary(c2Boundary);
- for(std::map<Cell*, int, Less_Cell>::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
+ 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::map<Cell*, int, Less_Cell >::iterator it = c2Boundary.begin(); it != c2Boundary.end(); it++){
+ 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;
@@ -290,7 +300,11 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
std::map<Cell*, int, Less_Cell >::iterator it2 = c1Boundary.find(cell);
bool old = false;
if(it2 != c1Boundary.end()) old = true;
- if(!old){ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this); }
+ if(!old){
+ if(this->addBoundaryCell(ori, cell)) {
+ cell->addCoboundaryCell(ori, this);
+ }
+ }
}
else{
if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
@@ -298,16 +312,20 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
}
// coboundary cells
- std::map<Cell*, int, Less_Cell > c1Coboundary = c1->getOrientedCoboundary();
- std::map<Cell*, int, Less_Cell > c2Coboundary = c2->getOrientedCoboundary();
+ std::map<Cell*, int, Less_Cell > c1Coboundary;
+ c1->getCoboundary(c1Coboundary);
+ std::map<Cell*, int, Less_Cell > c2Coboundary;
+ c2->getCoboundary(c2Coboundary);
- for(std::map<Cell*, int, Less_Cell>::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
+ 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::map<Cell*, int, Less_Cell>::iterator it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
+ 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;
@@ -316,7 +334,11 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
std::map<Cell*, int, Less_Cell >::iterator it2 = c1Coboundary.find(cell);
bool old = false;
if(it2 != c1Coboundary.end()) old = true;
- if(!old) { if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this); }
+ if(!old) {
+ if(this->addCoboundaryCell(ori, cell)){
+ cell->addBoundaryCell(ori, this);
+ }
+ }
}
else {
if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
@@ -325,13 +347,4 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
}
-
-CombinedCell::~CombinedCell(){
- std::list< std::pair<int, Cell*> > cells = getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell2 = (*it).second;
- delete cell2;
- }
-}
-
#endif
diff --git a/Geo/Cell.h b/Geo/Cell.h
index dd5b41bb3b2f5595df8f44e71062c32741c9f34e..f42d376cddcf2000bce15e510c44d8684a254c82 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -40,208 +40,203 @@ class Less_Cell{
// Class representing an elementary cell of a cell complex.
class Cell
{
- protected:
-
- // cell dimension
- int _dim;
-
- // whether this cell belongs to a subdomain, immutable
- // used in relative homology computation
- bool _inSubdomain;
-
- // whether this cell belongs to the boundary of a cell complex
- bool _onDomainBoundary;
-
- // whether this cell a combinded cell of elemetary cells
- bool _combined;
-
- // mutable index for each cell (used to create boundary operator matrices)
- int _index;
-
- // for some algorithms to omit this cell
- bool _immune;
+ protected:
- // mutable list of 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;
-
- // immutable 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;
-
- // The mesh element that is the image of this cell
- MElement* _image;
- // Whether to delete the mesh element when done (created for homology computation only)
- bool _deleteImage;
+ // cell dimension
+ int _dim;
- // sorted vertices of this cell (used for ordering of the cells)
- std::vector<int> _vs;
-
- public:
- Cell(){
- _combined = false;
- _index = 0;
- _immune = false;
- _image = NULL;
- _deleteImage = false;
- }
-
- Cell(MElement* image, bool subdomain, bool boundary){
- _combined = false;
- _index = 0;
- _immune = false;
- _image = image;
- _onDomainBoundary = boundary;
- _inSubdomain = subdomain;
- _dim = image->getDim();
- _deleteImage = false;
- for(int i = 0; i < image->getNumVertices(); i++) _vs.push_back(image->getVertex(i)->getNum());
- std::sort(_vs.begin(), _vs.end());
- }
- virtual ~Cell() {if(_deleteImage) delete _image; }
-
- virtual MElement* getImageMElement() const { return _image; };
-
- virtual int getDim() const { return _dim; };
- virtual int getIndex() const { return _index; };
- virtual void setIndex(int index) { _index = index; };
- virtual int getNum() { return _image->getNum(); }
- virtual int getType() { return _image->getType(); }
- virtual int getTypeForMSH() { return _image->getTypeForMSH(); }
- virtual int getPartition() { return _image->getPartition(); }
- virtual void setImmune(bool immune) { _immune = immune; };
- virtual bool getImmune() const { return _immune; };
- virtual void setDeleteImage(bool deleteImage) { _deleteImage = deleteImage; };
- virtual bool getDeleteImage() const { return _deleteImage; };
-
- // get the number of vertices this cell has
- virtual int getNumVertices() const { return _image->getNumVertices(); }
- virtual MVertex* getVertex(int vertex) const { return _image->getVertex(vertex); }
- virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
-
- // restores the cell information to its original state before reduction
- virtual void restoreCell();
-
- // true if this cell has given vertex
- virtual bool hasVertex(int vertex) const;
+ // whether this cell belongs to a subdomain, immutable
+ // used in relative homology computation
+ bool _inSubdomain;
+
+ // whether this cell belongs to the boundary of a cell complex
+ bool _onDomainBoundary;
+
+ // whether this cell a combinded cell of elemetary cells
+ bool _combined;
+
+ // mutable index for each cell (used to create boundary operator matrices)
+ int _index;
+
+ // for some algorithms to omit this cell
+ bool _immune;
+
+ // mutable list of 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;
+
+ // immutable 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;
+
+ // The mesh element that is the image of this cell
+ MElement* _image;
+ // Whether to delete the mesh element when done
+ // (created for homology computation only)
+ bool _deleteImage;
+
+ // sorted vertices of this cell (used for ordering of the cells)
+ std::vector<int> _vs;
+
+ public:
- // (co)boundary cell iterator
- typedef std::map<Cell*, int, Less_Cell>::iterator biter;
-
- virtual int getBoundarySize() { return _boundary.size(); }
- virtual int getCoboundarySize() { return _coboundary.size(); }
-
- // get the cell boundary
- virtual std::map<Cell*, int, Less_Cell > getOrientedBoundary() { return _boundary; }
- virtual std::list< Cell* > getBoundary();
- virtual std::map<Cell*, int, Less_Cell > getOrientedCoboundary() { return _coboundary; }
- virtual std::list< Cell* > getCoboundary();
-
- // add (co)boundary cell
- virtual bool addBoundaryCell(int orientation, Cell* cell);
- virtual bool addCoboundaryCell(int orientation, Cell* cell);
-
- // remove (co)boundary cell
- virtual int removeBoundaryCell(Cell* cell, bool other=true);
- virtual int removeCoboundaryCell(Cell* cell, bool other=true);
-
- // true if has given cell on (original) (co)boundary
- virtual bool hasBoundary(Cell* cell, bool org=false);
- virtual bool hasCoboundary(Cell* cell, bool org=false);
-
- virtual void clearBoundary() { _boundary.clear(); }
- virtual void clearCoboundary() { _coboundary.clear(); }
-
+ Cell() : _combined(false), _index(0), _immune(false), _image(NULL),
+ _deleteImage(false) {}
+ Cell(MElement* image, bool subdomain, bool boundary);
+
+ virtual ~Cell();
+
+ virtual MElement* getImageMElement() const { return _image; };
+
+ virtual int getDim() const { return _dim; };
+ virtual int getIndex() const { return _index; };
+ virtual void setIndex(int index) { _index = index; };
+ virtual int getNum() { return _image->getNum(); }
+ virtual int getType() { return _image->getType(); }
+ virtual int getTypeForMSH() { return _image->getTypeForMSH(); }
+ virtual int getPartition() { return _image->getPartition(); }
+ virtual void setImmune(bool immune) { _immune = immune; };
+ virtual bool getImmune() const { return _immune; };
+ virtual void setDeleteImage(bool deleteImage) {
+ _deleteImage = deleteImage; };
+ virtual bool getDeleteImage() const { return _deleteImage; };
+
+ // get the number of vertices this cell has
+ virtual int getNumVertices() const { return _image->getNumVertices(); }
+ virtual MVertex* getVertex(int vertex) const {
+ return _image->getVertex(vertex); }
+ virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
+
+ // restores the cell information to its original state before reduction
+ virtual void restoreCell();
+
+ // true if this cell has given vertex
+ virtual bool hasVertex(int vertex) const;
+
+ // (co)boundary cell iterator
+ typedef std::map<Cell*, int, Less_Cell>::iterator biter;
+
+ virtual int getBoundarySize() { return _boundary.size(); }
+ virtual int getCoboundarySize() { return _coboundary.size(); }
+
+ // get the cell boundary
+ virtual void getBoundary(std::map<Cell*, int, Less_Cell >& boundary){
+ boundary = _boundary; }
+ virtual void getCoboundary( std::map<Cell*, int, Less_Cell >& coboundary ){
+ coboundary = _coboundary; }
+
+ // add (co)boundary cell
+ virtual bool addBoundaryCell(int orientation, Cell* cell);
+ virtual bool addCoboundaryCell(int orientation, Cell* cell);
+
+ // remove (co)boundary cell
+ virtual int removeBoundaryCell(Cell* cell, bool other=true);
+ virtual int removeCoboundaryCell(Cell* cell, bool other=true);
+
+ // true if has given cell on (original) (co)boundary
+ virtual bool hasBoundary(Cell* cell, bool org=false);
+ virtual bool hasCoboundary(Cell* cell, bool org=false);
+
+ virtual void clearBoundary() { _boundary.clear(); }
+ virtual void clearCoboundary() { _coboundary.clear(); }
+
// algebraic dual of the cell
- virtual void makeDualCell();
-
- // print cell info
- virtual void printCell();
- virtual void printBoundary();
- virtual void printCoboundary();
-
- // original (co)boundary
- 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();
- virtual void printOrgCbd();
-
- virtual bool inSubdomain() const { return _inSubdomain; }
- virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
-
- virtual bool onDomainBoundary() const { return _onDomainBoundary; }
- virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
-
- // get the number of facets of this cell
- virtual int getNumFacets() const;
- // get the vertices on a facet of this cell
- virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
-
- // get boundary cell orientation
- virtual int getFacetOri(Cell* cell) {
+ virtual void makeDualCell();
+
+ // print cell info
+ virtual void printCell();
+ virtual void printBoundary();
+ virtual void printCoboundary();
+
+ // original (co)boundary
+ 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 void getOrgBd( std::map<Cell*, int, Less_Cell >& obd) {
+ obd = _obd; }
+ virtual void getOrgCbd( std::map<Cell*, int, Less_Cell >& obdc) {
+ obdc = _ocbd; }
+ virtual void printOrgBd();
+ virtual void printOrgCbd();
+
+ virtual bool inSubdomain() const { return _inSubdomain; }
+ virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
+
+ virtual bool onDomainBoundary() const { return _onDomainBoundary; }
+ virtual void setOnDomainBoundary(bool domainboundary) {
+ _onDomainBoundary = domainboundary; }
+
+ // get the number of facets of this cell
+ virtual int getNumFacets() const;
+ // get the vertices on a facet of this cell
+ virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
+
+ // get boundary cell orientation
+ virtual int getFacetOri(Cell* cell) {
std::vector<MVertex*> v;
- for(int i = 0; i < cell->getNumVertices(); i++) v.push_back(cell->getVertex(i));
+ for(int i = 0; i < cell->getNumVertices(); i++) {
+ v.push_back(cell->getVertex(i));
+ }
return getFacetOri(v);
- }
- virtual int getFacetOri(std::vector<MVertex*> &v);
-
- // tools for combined cells
- 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;}
-
- // equivalence
- bool operator==(const Cell& c2) const {
- if(this->getNumVertices() != c2.getNumVertices()){
- return false;
- }
- for(int i=0; i < this->getNumVertices();i++){
- if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
- return false;
- }
- }
- return true;
- }
+ }
+ virtual int getFacetOri(std::vector<MVertex*> &v);
+
+ // tools for combined cells
+ virtual bool isCombined() { return _combined; }
+ virtual void getCells( std::list< std::pair<int, Cell*> >& cells) {
+ cells.clear();
+ cells.push_back( std::make_pair(1, this));
+ return;
+ }
+ virtual int getNumCells() {return 1;}
+
+ // equivalence
+ bool operator==(const Cell& c2) const {
+ if(this->getNumVertices() != c2.getNumVertices()){
+ return false;
+ }
+ for(int i=0; i < this->getNumVertices();i++){
+ if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
+ return false;
+ }
+ }
+ return true;
+ }
};
// A cell that is a combination of cells of same dimension
class CombinedCell : public Cell{
-
- private:
- // vertices
- std::vector<MVertex*> _v;
- // sorted list of all vertices
- std::vector<int> _vs;
- // list of cells this cell is a combination of
- std::list< std::pair<int, Cell*> > _cells;
- int _num;
-
- public:
-
- CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
- ~CombinedCell();
-
- int getNum() { return _num; }
- int getType() { return 0; }
- int getTypeForMSH() { return 0; }
- int getPartition() { return 0; }
-
- int getNumVertices() const { return _v.size(); }
- MVertex* getVertex(int vertex) const { return _v.at(vertex); }
- int getSortedVertex(int vertex) const { return _vs.at(vertex); }
-
- std::list< std::pair<int, Cell*> > getCells() { return _cells; }
- int getNumCells() {return _cells.size();}
-
+
+ private:
+ // vertices
+ std::vector<MVertex*> _v;
+ // sorted list of all vertices
+ std::vector<int> _vs;
+ // list of cells this cell is a combination of
+ std::list< std::pair<int, Cell*> > _cells;
+ int _num;
+
+ public:
+
+ CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
+ ~CombinedCell() {}
+
+ int getNum() { return _num; }
+ int getType() { return 0; }
+ int getTypeForMSH() { return 0; }
+ int getPartition() { return 0; }
+
+ int getNumVertices() const { return _v.size(); }
+ MVertex* getVertex(int vertex) const { return _v.at(vertex); }
+ int getSortedVertex(int vertex) const { return _vs.at(vertex); }
+
+ void getCells(std::list< std::pair<int, Cell*> >& cells) { cells = _cells; }
+ int getNumCells() {return _cells.size();}
+
};
-
#endif
#endif
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index cd88192b041548c497c53806307439547d7d9018..e2112bf73831fd780e3e27f52658fe3497333f76 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -47,17 +47,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
-void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
-
+void CellComplex::find_boundary(std::vector<GEntity*>& domain,
+ std::vector<GEntity*>& subdomain){
+
// determine mesh entities on boundary of the domain
bool duplicate = false;
for(unsigned int j=0; j < _domain.size(); j++){
if(_domain.at(j)->dim() == 3){
std::list<GFace*> faces = _domain.at(j)->faces();
- for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+ for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end();
+ it++){
GFace* face = *it;
duplicate = false;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin();
+ itb != _boundary.end(); itb++){
GEntity* entity = *itb;
if(face->tag() == entity->tag()){
_boundary.erase(itb);
@@ -71,11 +74,13 @@ void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEnti
else if(_domain.at(j)->dim() == 2){
std::list<GEdge*> edges = _domain.at(j)->edges();
- for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++){
+ for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end();
+ it++){
GEdge* edge = *it;
duplicate = false;
std::vector<GEntity*>::iterator erase;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin();
+ itb != _boundary.end(); itb++){
GEntity* entity = *itb;
if(edge->tag() == entity->tag()){
_boundary.erase(itb);
@@ -88,10 +93,12 @@ void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEnti
}
else if(_domain.at(j)->dim() == 1){
std::list<GVertex*> vertices = _domain.at(j)->vertices();
- for(std::list<GVertex*>::iterator it = vertices.begin(); it != vertices.end(); it++){
+ for(std::list<GVertex*>::iterator it = vertices.begin();
+ it != vertices.end(); it++){
GVertex* vertex = *it;
duplicate = false;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin();
+ itb != _boundary.end(); itb++){
GEntity* entity = *itb;
if(vertex->tag() == entity->tag()){
_boundary.erase(itb);
@@ -132,10 +139,13 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
int type = domain.at(j)->getMeshElement(i)->getTypeForMSH();
Cell* cell;
// simplex types
- if(type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4 || type == MSH_LIN_3 || type == MSH_TRI_6 ||
- type == MSH_TET_10 || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10 || type == MSH_TRI_12 ||
- type == MSH_TRI_15 || type == MSH_TRI_15I || type == MSH_TRI_21 || type == MSH_LIN_4 ||
- type == MSH_LIN_5 || type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35 || type == MSH_TET_56
+ if(type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4
+ || type == MSH_LIN_3 || type == MSH_TRI_6 || type == MSH_TET_10
+ || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10
+ || type == MSH_TRI_12 || type == MSH_TRI_15 || type == MSH_TRI_15I
+ || type == MSH_TRI_21 || type == MSH_LIN_4 || type == MSH_LIN_5
+ || type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35
+ || type == MSH_TET_56
|| type == MSH_TET_34 || type == MSH_TET_52 ){
cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
}
@@ -193,10 +203,11 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
/*else if(type == MSH_LIN_3 || type == MSH_LIN_4 ||
type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
}
- if(newtype == 0) Msg::Error("Error: mesh element %d not implemented yet! \n", type);
-
- //printf("dim: %d type: %d, newtype: %d, vertices: %d \n", dim, type, newtype, vertices.size());
- MElement* element = factory.create(newtype, vertices, 0, cell->getPartition());
+ if(newtype == 0){
+ Msg::Error("Error: mesh element %d not implemented yet! \n", type); }
+
+ MElement* element = factory.create(newtype, vertices, 0,
+ cell->getPartition());
Cell* newCell = new Cell(element, subdomain, boundary);
newCell->setImmune(cell->getImmune());
newCell->setDeleteImage(true);
@@ -235,37 +246,45 @@ CellComplex::~CellComplex(){
_ocells[i].clear();
}
+ for(int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
+ _newcells.clear();
}
/*
- CellComplex::CellComplex(CellComplex* cellComplex){
- *
- _domain = cellComplex->getDomain();
- _subdomain = cellComplex->getSubdomain;
- _boundary = cellComplex->getBoundary;
- _domainVertices = cellComplex->getDomainVertices;
- *
- for(int i = 0; i < 4; i++){
- _betti[i] = cellComplex->getBetti(i);
- _cells[i] = ocells[i];
- _ocells[i] = ocells[i];
- }
- *
- _dim = cellComplex->getDim();
- }
+CellComplex::CellComplex(CellComplex* cellComplex){
+ _domain = cellComplex->getDomain();
+ _subdomain = cellComplex->getSubdomain();
+ _boundary = cellComplex->getBoundary();
+
+ for(int i = 0; i < 4; i++){
+ _betti[i] = cellComplex->getBetti(i);
+ _ocells[i] = cellComplex->getOcells(i);
+ }
+ _dim = cellComplex->getDim();
+ restoreComplex();
+}
*/
+void CellComplex::insertCell(Cell* cell){
+ _newcells.push_back(cell);
+ std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
+ if(!insertInfo.second) Msg::Debug("Warning: Cell not inserted! \n");
+}
void CellComplex::removeCell(Cell* cell, bool other){
- std::map<Cell*, int, Less_Cell > coboundary = cell->getOrientedCoboundary();
- std::map<Cell*, int, Less_Cell > boundary = cell->getOrientedBoundary();
+ std::map<Cell*, int, Less_Cell > coboundary;
+ cell->getCoboundary(coboundary);
+ std::map<Cell*, int, Less_Cell > boundary;
+ cell->getBoundary(boundary);
- for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin();
+ it != coboundary.end(); it++){
Cell* cbdCell = (*it).first;
cbdCell->removeBoundaryCell(cell, other);
}
- for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){
+ for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin();
+ it != boundary.end(); it++){
Cell* bdCell = (*it).first;
bdCell->removeCoboundaryCell(cell, other);
}
@@ -274,13 +293,17 @@ void CellComplex::removeCell(Cell* cell, bool other){
}
-void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
- Qset.erase(cell);
+void CellComplex::removeCellQset(Cell*& cell,
+ std::set<Cell*, Less_Cell>& Qset){
+ Qset.erase(cell);
}
-void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
- for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
- Cell* cell = *cit;
+void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
+ std::queue<Cell*>& Q,
+ std::set<Cell*, Less_Cell>& Qset){
+ for(std::map<Cell*, int, Less_Cell>::iterator cit = cells.begin();
+ cit != cells.end(); cit++){
+ Cell* cell = (*cit).first;
citer it = Qset.find(cell);
//citer it2 = _cells[cell->getDim()].find(cell);
if(it == Qset.end()){// && it2 != _cells[cell->getDim()].end()){
@@ -300,9 +323,10 @@ int CellComplex::coreduction(Cell* generator, int omitted){
Q.push(generator);
Qset.insert(generator);
-
- std::list<Cell*> bd_s;
- std::list<Cell*> cbd_c;
+ std::map<Cell*, int, Less_Cell > bd_s;
+ std::map<Cell*, int, Less_Cell > cbd_c;
+ //std::list<Cell*> bd_s;
+ //std::list<Cell*> cbd_c;
Cell* s;
int round = 0;
while( !Q.empty() ){
@@ -312,20 +336,22 @@ int CellComplex::coreduction(Cell* generator, int omitted){
s = Q.front();
Q.pop();
removeCellQset(s, Qset);
- bd_s = s->getBoundary();
+ s->getBoundary(bd_s);
- if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
+ if( bd_s.size() == 1 && inSameDomain(s, bd_s.begin()->first) ){
removeCell(s);
- cbd_c = bd_s.front()->getCoboundary();
+ bd_s.begin()->first->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
- removeCell(bd_s.front());
- if(bd_s.front()->getDim() == 0 && omitted > 0) _store.at(omitted-1).insert(bd_s.front());
+ removeCell(bd_s.begin()->first);
+ if(bd_s.begin()->first->getDim() == 0 && omitted > 0){
+ _store.at(omitted-1).insert(bd_s.begin()->first);
+ }
coreductions++;
}
else if(bd_s.empty()){
- cbd_c = s->getCoboundary();
+ s->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
}
@@ -337,7 +363,8 @@ int CellComplex::coreduction(Cell* generator, int omitted){
int CellComplex::reduction(int dim, int omitted){
if(dim < 1 || dim > 3) return 0;
- std::list<Cell*> cbd_c;
+ //std::list<Cell*> cbd_c;
+ std::map<Cell*, int, Less_Cell > cbd_c;
int count = 0;
bool reduced = true;
@@ -350,13 +377,15 @@ int CellComplex::reduction(int dim, int omitted){
Cell* cell = *cit;
- cbd_c = cell->getCoboundary();
- if( cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front())){
+ cell->getCoboundary(cbd_c);
+ if( cbd_c.size() == 1 && inSameDomain(cell, cbd_c.begin()->first)){
//&& ( (!cell->getImmune() && !cbd_c.front()->getImmune() ) )){
++cit;
- removeCell(cbd_c.front());
+ removeCell(cbd_c.begin()->first);
removeCell(cell);
- if(dim == getDim() && omitted > 0) _store.at(omitted-1).insert(cbd_c.front());
+ if(dim == getDim() && omitted > 0){
+ _store.at(omitted-1).insert(cbd_c.begin()->first);
+ }
count++;
reduced = true;
@@ -367,7 +396,7 @@ int CellComplex::reduction(int dim, int omitted){
}
//if(!reduced && ignoreCells) { ignoreCells = false; reduced = true;}
-
+
}
return count;
@@ -611,21 +640,21 @@ int CellComplex::cocombine(int dim){
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
- std::list<Cell*> cbd_c;
+ std::map<Cell*, int, Less_Cell> cbd_c;
std::map<Cell*, int, Less_Cell > bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
- cbd_c = cell->getCoboundary();
+ cell->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
while(Q.size() != 0){
Cell* s = Q.front();
Q.pop();
- bd_c = s->getOrientedBoundary();
+ s->getBoundary(bd_c);
if(s->getBoundarySize() == 2){
std::map<Cell*, int, Less_Cell>::iterator it = bd_c.begin();
@@ -639,19 +668,19 @@ int CellComplex::cocombine(int dim){
&& inSameDomain(s, c1) && inSameDomain(s, c2)
&& c1->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
&& c2->getNumVertices() < getSize(dim)){
-
+
removeCell(s);
- cbd_c = c1->getCoboundary();
+ c1->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
- cbd_c = c2->getCoboundary();
+ c2->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
- CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
+ CombinedCell* newCell = new CombinedCell(c1, c2,
+ (or1 != or2), true );
removeCell(c1);
removeCell(c2);
- std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
- if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
+ insertCell(newCell);
cit = firstCell(dim);
count++;
@@ -678,20 +707,20 @@ int CellComplex::combine(int dim){
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::map<Cell*, int, Less_Cell> cbd_c;
- std::list<Cell*> bd_c;
+ std::map<Cell*, int, Less_Cell> bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
- bd_c = cell->getBoundary();
+ cell->getBoundary(bd_c);
enqueueCells(bd_c, Q, Qset);
while(Q.size() != 0){
Cell* s = Q.front();
Q.pop();
- cbd_c = s->getOrientedCoboundary();
+ s->getCoboundary(cbd_c);
if(s->getCoboundarySize() == 2){
@@ -709,17 +738,16 @@ int CellComplex::combine(int dim){
removeCell(s);
- bd_c = c1->getBoundary();
+ c1->getBoundary(bd_c);
enqueueCells(bd_c, Q, Qset);
- bd_c = c2->getBoundary();
+ c2->getBoundary(bd_c);
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) Msg::Debug("Warning: Combined cell not inserted! \n");
+ insertCell(newCell);
+
cit = firstCell(dim);
count++;
}
@@ -752,7 +780,7 @@ int CellComplex::combine(int dim){ // version of combine that doesn't have a que
combined = false;
for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
Cell* s = *cit;
- cbd_c = s->getOrientedCoboundary();
+ cbd_c = s->getCoboundary();
if(s->getCoboundarySize() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
&& inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second)
@@ -831,8 +859,10 @@ bool CellComplex::checkCoherence(){
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
- std::map<Cell*, int, Less_Cell> boundary = cell->getOrientedBoundary();
- for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){
+ std::map<Cell*, int, Less_Cell> boundary;
+ cell->getBoundary(boundary);
+ 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);
@@ -854,8 +884,10 @@ bool CellComplex::checkCoherence(){
}
}
- std::map<Cell*, int, Less_Cell> coboundary = cell->getOrientedCoboundary();
- for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+ std::map<Cell*, int, Less_Cell> coboundary;
+ cell->getCoboundary(coboundary);
+ 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);
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 017925888365673ca31c0b2c4782cfd4c7201318..e0cb5d8592df6c8b88980347958d2b7bff6a99fd 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -31,131 +31,141 @@ class Cell;
class CellComplex
{
private:
-
- // the domain in the model which this cell complex covers
- std::vector<GEntity*> _domain;
- // a subdomain of the given domain
- // used in relative homology computation, may be empty
- std::vector<GEntity*> _subdomain;
-
- // entities on the boundary of the homology computation domain
- std::vector<GEntity*> _boundary;
-
- // sorted containers of unique cells in this cell complex
- // one for each dimension
- std::set<Cell*, Less_Cell> _cells[4];
-
- // temporary store for omitted cells (generators of the highest dimension)
- std::vector< std::set<Cell*, Less_Cell> > _store;
-
- // original cells of this cell complex
- std::set<Cell*, Less_Cell> _ocells[4];
-
- // Betti numbers of this cell complex (ranks of homology groups)
- int _betti[4];
-
- int _dim;
-
- // is the cell complex simplicial
- bool _simplicial;
-
- // does the domain contain entities of different dimensions
- bool _multidim;
-
- // enqueue cells in queue if they are not there already
- void enqueueCells(std::list<Cell*>& cells,
- std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
- // remove cell from the queue set
- void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
-
- // for constructor
- void insert_cells(bool subdomain, bool boundary);
- void find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain);
-
- // remove a cell from this cell complex
- void removeCell(Cell* cell, bool other=true);
-
- // queued coreduction presented in Mrozek's paper
- int coreduction(Cell* generator, int omitted=0);
-
- public:
-
- CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
- //CellComplex(CellComplex* cellComplex)
- CellComplex(){}
- ~CellComplex();
-
- //void emptyTrash();
-
- // restore this cell complex to its original state
- void restoreComplex();
-
- // 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; }
-
- bool simplicial() { return _simplicial; }
-
- std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
- std::set<Cell*, Less_Cell> getOrgCells(int dim){ return _ocells[dim]; }
-
- // iterator for the cells of same dimension
- typedef std::set<Cell*, Less_Cell>::iterator citer;
-
- // 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, bool org=false);
-
- // check whether two cells both belong to subdomain or if neither one does
- bool inSameDomain(Cell* c1, Cell* c2) const { return
- ( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
-
- // reduction of this cell complex
- // removes reduction pairs of cell of dimension dim and dim-1
- int reduction(int dim, int omitted=0);
+
+ // the domain in the model which this cell complex covers
+ std::vector<GEntity*> _domain;
+ // a subdomain of the given domain
+ // used in relative homology computation, may be empty
+ std::vector<GEntity*> _subdomain;
+
+ // entities on the boundary of the homology computation domain
+ std::vector<GEntity*> _boundary;
+
+ // sorted containers of unique cells in this cell complex
+ // one for each dimension
+ std::set<Cell*, Less_Cell> _cells[4];
+
+ // temporary store for omitted cells (generators of the highest dimension)
+ std::vector< std::set<Cell*, Less_Cell> > _store;
+
+ // original cells of this cell complex
+ std::set<Cell*, Less_Cell> _ocells[4];
+
+ // new cells created during reductions
+ std::vector<Cell*> _newcells;
+
+ // Betti numbers of this cell complex (ranks of homology groups)
+ int _betti[4];
+
+ int _dim;
+
+ // is the cell complex simplicial
+ bool _simplicial;
+
+ // does the domain contain entities of different dimensions
+ bool _multidim;
+
+ // enqueue cells in queue if they are not there already
+ void enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
+ std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
+ // remove cell from the queue set
+ void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
+
+ // for constructor
+ void insert_cells(bool subdomain, bool boundary);
+ void find_boundary(std::vector<GEntity*>& domain,
+ std::vector<GEntity*>& subdomain);
+
+ // insert/remove a cell from this cell complex
+ void removeCell(Cell* cell, bool other=true);
+ void insertCell(Cell* cell);
+
+ // queued coreduction presented in Mrozek's paper
+ int coreduction(Cell* generator, int omitted=0);
+
+ public:
+
+ CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
+ //CellComplex(CellComplex* cellComplex)
+ ~CellComplex();
+
+ // restore this cell complex to its original state
+ void restoreComplex();
- // useful functions for (co)reduction of cell complex
- int reduceComplex();
- int coreduceComplex();
-
- // remove cells in subdomain from this cell complex
- void removeSubdomain();
-
- // print the vertices of cells of certain dimension
- void printComplex(int dim);
-
- // Cell combining for reduction and coreduction
- int combine(int dim);
- int cocombine(int dim);
-
- // Compute betti numbers of this cell complex
- 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();
-
- // make cells on the domain boundary that are not in the subdomain subdomain cells and vice versa
- bool swapSubdomain();
-
- int eulerCharacteristic(){ return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
- void printEuler(){ printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
-
- int getNumOmitted() { return _store.size(); };
- std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
-
- // change roles of boundaries and coboundaries of the cells in this cell complex
- // equivalent to transposing boundary operator matrices
- void makeDualComplex();
+ // 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; }
+
+ bool simplicial() { return _simplicial; }
+
+ std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
+ std::set<Cell*, Less_Cell> getOrgCells(int dim){ return _ocells[dim]; }
+
+ // iterator for the cells of same dimension
+ typedef std::set<Cell*, Less_Cell>::iterator citer;
+
+ // 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, bool org=false);
+
+ // check whether two cells both belong to subdomain or if neither one does
+ bool inSameDomain(Cell* c1, Cell* c2) const {
+ return ( (!c1->inSubdomain() && !c2->inSubdomain())
+ || (c1->inSubdomain() && c2->inSubdomain()) ); }
+
+ // reduction of this cell complex
+ // removes reduction pairs of cell of dimension dim and dim-1
+ int reduction(int dim, int omitted=0);
+
+ // useful functions for (co)reduction of cell complex
+ int reduceComplex();
+ int coreduceComplex();
+
+ // remove cells in subdomain from this cell complex
+ void removeSubdomain();
+
+ // print the vertices of cells of certain dimension
+ void printComplex(int dim);
+
+ // Cell combining for reduction and coreduction
+ int combine(int dim);
+ int cocombine(int dim);
+
+ // Compute betti numbers of this cell complex
+ 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();
+
+ // make cells on the domain boundary that are not in the subdomain
+ // subdomain cells and vice versa
+ bool swapSubdomain();
+
+ int eulerCharacteristic(){
+ return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
+ void printEuler(){
+ printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
+
+ int getNumOmitted() { return _store.size(); };
+ std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
+
+ // change roles of boundaries and coboundaries of the cells in this
+ // cell complex
+ // equivalent to transposing boundary operator matrices
+ void makeDualComplex();
};
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index cbfd41bb1897c5d54cb0f17563f4473222007dcd..dde548d1dd9d9fc15c22ac5e55763e8c9e017d12 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -71,27 +71,36 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
_HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
//printMatrix(_HMatrix[dim]);
- for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
+ for( std::set<Cell*, Less_Cell>::iterator cit =
+ cellComplex->firstCell(dim);
+ cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if(!cell->inSubdomain()){
- std::map<Cell*, int, Less_Cell> bdCell = cell->getOrientedBoundary();
- for(std::map<Cell*, int, Less_Cell>::iterator it = bdCell.begin(); it != bdCell.end(); it++){
+ std::map<Cell*, int, Less_Cell> bdCell;
+ cell->getBoundary(bdCell);
+ 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());
- if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1
- || cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){
+ if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim])
+ || bdCell->getIndex() < 1
+ || cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim])
+ || cell->getIndex() < 1){
Msg::Debug("Warning: Index out of bound! HMatrix: %d. \n", dim);
}
else{
- gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ gmp_matrix_get_elem(elem, bdCell->getIndex(),
+ cell->getIndex(), _HMatrix[dim]);
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).second);
if( abs((old_elem + (*it).second)) > 1){
- Msg::Debug("Incidence index: %d! HMatrix: %d.", (old_elem + (*it).second), dim);
+ Msg::Debug("Incidence index: %d! HMatrix: %d.",
+ (old_elem + (*it).second), dim);
}
- gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ gmp_matrix_set_elem(elem, bdCell->getIndex(),
+ cell->getIndex(), _HMatrix[dim]);
}
}
}
@@ -415,15 +424,19 @@ int ChainComplex::getTorsion(int dim, int chainNumber){
}
-Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, GModel* model, std::string name, int torsion){
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
+ CellComplex* cellComplex, GModel* model,
+ std::string name, int torsion){
int i = 0;
- for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
+ 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){
- std::list< std::pair<int, Cell*> > subCells = cell->getCells();
+ std::list< std::pair<int, Cell*> > subCells;
+ cell->getCells(subCells);
for(std::list< std::pair<int, Cell*> >::iterator it = subCells.begin(); it != subCells.end(); it++){
Cell* subCell = (*it).second;
int coeff = (*it).first;
@@ -441,7 +454,8 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
}
-bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*, int, Less_Cell> &cellsNotInChain){
+bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
+ std::map<Cell*, int, Less_Cell> &cellsNotInChain){
std::vector<int> cc;
std::vector<int> bc;
@@ -461,7 +475,9 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*
if(cc[i]*bc[i] != inout) return false;
}
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++) removeCell((*cit).first);
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ removeCell((*cit).first);
+ }
int n = 1;
for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end(); cit++){
@@ -479,18 +495,23 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*
bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){
Cell* c1 = cell.first;
- std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
+ std::map<Cell*, int, Less_Cell> c1Cbd;
+ c1->getOrgCbd(c1Cbd);
for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
std::map<Cell*, int, Less_Cell> cellsInChain;
std::map<Cell*, int, Less_Cell> cellsNotInChain;
Cell* c1CbdCell = (*cit).first;
- std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
+ std::map<Cell*, int, Less_Cell> c1CbdBd;
+ c1CbdCell->getOrgBd(c1CbdBd);
for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
Cell* c1CbdBdCell = (*cit2).first;
int coeff = (*cit2).second;
- if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0) || c1CbdBdCell->inSubdomain()) cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0)
+ || c1CbdBdCell->inSubdomain()){
+ cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ }
else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
}
@@ -504,7 +525,8 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){
if(coeff > 1 || coeff < -1) next = true;
}
- for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end(); cit2++){
+ for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
+ cit2++){
Cell* c = (*cit2).first;
if(c->inSubdomain()) next = true;
}
@@ -513,18 +535,24 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){
//printf("dim: %d, in chain: %d, not in chain: %d \n", getDim(), cellsInChain.size(), cellsNotInChain.size());
- if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1) ||
- (getDim() == 2 && cellsInChain.size() == 3 && cellsNotInChain.size() == 1)){
+ if( (getDim() == 1 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 1) ||
+ (getDim() == 2 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 1)){
//printf("straighten \n");
return deform(cellsInChain, cellsNotInChain);
}
- else if ( (getDim() == 1 && cellsInChain.size() == 1 && cellsNotInChain.size() == 2 && bend) ||
- (getDim() == 2 && cellsInChain.size() == 2 && cellsNotInChain.size() == 2 && bend)){
+ else if ( (getDim() == 1 && cellsInChain.size() == 1
+ && cellsNotInChain.size() == 2 && bend) ||
+ (getDim() == 2 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 2 && bend)){
//printf("bend \n");
return deform(cellsInChain, cellsNotInChain);
}
- else if ((getDim() == 1 && cellsInChain.size() == 3 && cellsNotInChain.size() == 0) ||
- (getDim() == 2 && cellsInChain.size() == 4 && cellsNotInChain.size() == 0)){
+ else if ((getDim() == 1 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 0) ||
+ (getDim() == 2 && cellsInChain.size() == 4
+ && cellsNotInChain.size() == 0)){
//printf("remove boundary \n");
return deform(cellsInChain, cellsNotInChain);
}
@@ -652,9 +680,9 @@ void Chain::createPView(){
_model->setPhysicalName(getName(), getDim(), physicalNum);
// only for visualization
- PView *chain = new PView(getName(), "ElementData", getGModel(), data, 0, 1);
+ PView* chain = new PView(getName(), "ElementData", getGModel(), data, 0, 1);
}
-
+
return;
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 3897b6cc824d2737d05e6d88373992e0fdf575e4..5b85d8a1d867c64bbf055297ebaf5d7d88704a10 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -32,190 +32,213 @@ extern "C" {
// A class representing a chain complex of a cell complex.
// This should only be constructed for a reduced cell complex because of
-// dense matrix representations and great computational complexity in its methods.
+// dense matrix representations and great computational complexity in
+// its methods.
class ChainComplex{
- private:
- // boundary operator matrices for this chain complex
- // h_k: C_k -> C_(k-1)
- gmp_matrix* _HMatrix[5];
-
- // Basis matrices for the kernels and codomains of the boundary operator matrices
- gmp_matrix* _kerH[5];
- gmp_matrix* _codH[5];
-
- gmp_matrix* _JMatrix[5];
- gmp_matrix* _QMatrix[5];
- std::vector<long int> _torsion[5];
-
- // bases for homology groups
- gmp_matrix* _Hbasis[5];
-
- int _dim;
- CellComplex* _cellComplex;
-
- // set the matrices
- void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
- void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
- void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _codH[dim] = matrix;}
- void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
- void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
- void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
-
-
-
- public:
-
- ChainComplex(CellComplex* cellComplex);
-
- ChainComplex(){
- for(int i = 0; i < 5; i++){
- _HMatrix[i] = create_gmp_matrix_zero(1,1);
- _kerH[i] = NULL;
- _codH[i] = NULL;
- _JMatrix[i] = NULL;
- _QMatrix[i] = NULL;
- _Hbasis[i] = NULL;
- }
- _dim = 0;
- }
- ~ChainComplex(){
- for(int i = 0; i < 5; i++){
- destroy_gmp_matrix(_HMatrix[i]);
- destroy_gmp_matrix(_kerH[i]);
- destroy_gmp_matrix(_codH[i]);
- destroy_gmp_matrix(_JMatrix[i]);
- destroy_gmp_matrix(_QMatrix[i]);
- destroy_gmp_matrix(_Hbasis[i]);
- }
- }
-
- int getDim() { return _dim; }
-
- // get the boundary operator matrix dim->dim-1
- gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
- gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
- gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
- gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
- gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
- gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
-
- // Compute basis for kernel and codomain of boundary operator matrix of dimension dim (ie. ker(h_dim) and cod(h_dim) )
- void KerCod(int dim);
- // Compute matrix representation J for inclusion relation from dim-cells who are boundary of dim+1-cells
- // to cycles of dim-cells (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
- void Inclusion(int lowDim, int highDim);
- // Compute quotient problem for given inclusion relation j to find representatives of homology groups
- // and possible torsion coeffcients
- void Quotient(int dim);
-
- // transpose the boundary operator matrices, these are boundary operator matrices for the dual mesh
- void transposeHMatrices() { for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
- void transposeHMatrix(int dim) { if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
-
- // Compute bases for the homology groups of this chain complex
- void computeHomology(bool dual=false);
-
- // get coefficient vector for dim-dimensional chain chainNumber (columns of _Hbasis[dim])
- std::vector<int> getCoeffVector(int dim, int chainNumber);
- // torsion coefficient for dim-dimensional chain chainNumber
- int getTorsion(int dim, int chainNumber);
- // get rank of homology group of dimension dim
- int getBasisSize(int dim) { if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
-
- int printMatrix(gmp_matrix* matrix){
- printf("%d rows and %d columns\n", (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
- return gmp_matrix_printf(matrix); }
-
- // debugging aid
- void matrixTest();
+ private:
+ // boundary operator matrices for this chain complex
+ // h_k: C_k -> C_(k-1)
+ gmp_matrix* _HMatrix[5];
+
+ // Basis matrices for the kernels and codomains of the boundary operator
+ // matrices
+ gmp_matrix* _kerH[5];
+ gmp_matrix* _codH[5];
+
+ gmp_matrix* _JMatrix[5];
+ gmp_matrix* _QMatrix[5];
+ std::vector<long int> _torsion[5];
+
+ // bases for homology groups
+ gmp_matrix* _Hbasis[5];
+
+ int _dim;
+ CellComplex* _cellComplex;
+
+ // set the matrices
+ void setHMatrix(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
+ void setKerHMatrix(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
+ void setCodHMatrix(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _codH[dim] = matrix;}
+ void setJMatrix(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
+ void setQMatrix(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
+ void setHbasis(int dim, gmp_matrix* matrix) {
+ if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
+
+ public:
+
+ ChainComplex(CellComplex* cellComplex);
+
+ ChainComplex(){
+ for(int i = 0; i < 5; i++){
+ _HMatrix[i] = create_gmp_matrix_zero(1,1);
+ _kerH[i] = NULL;
+ _codH[i] = NULL;
+ _JMatrix[i] = NULL;
+ _QMatrix[i] = NULL;
+ _Hbasis[i] = NULL;
+ }
+ _dim = 0;
+ }
+ ~ChainComplex(){
+ for(int i = 0; i < 5; i++){
+ destroy_gmp_matrix(_HMatrix[i]);
+ destroy_gmp_matrix(_kerH[i]);
+ destroy_gmp_matrix(_codH[i]);
+ destroy_gmp_matrix(_JMatrix[i]);
+ destroy_gmp_matrix(_QMatrix[i]);
+ destroy_gmp_matrix(_Hbasis[i]);
+ }
+ }
+
+ int getDim() { return _dim; }
+
+ // get the boundary operator matrix dim->dim-1
+ gmp_matrix* getHMatrix(int dim) {
+ if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
+ gmp_matrix* getKerHMatrix(int dim) {
+ if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
+ gmp_matrix* getCodHMatrix(int dim) {
+ if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
+ gmp_matrix* getJMatrix(int dim) {
+ if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
+ gmp_matrix* getQMatrix(int dim) {
+ if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
+ gmp_matrix* getHbasis(int dim) {
+ if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
+
+ // Compute basis for kernel and codomain of boundary operator matrix
+ // of dimension dim (ie. ker(h_dim) and cod(h_dim) )
+ void KerCod(int dim);
+ // Compute matrix representation J for inclusion relation from dim-cells
+ // who are boundary of dim+1-cells to cycles of dim-cells
+ // (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
+ void Inclusion(int lowDim, int highDim);
+ // Compute quotient problem for given inclusion relation j to find
+ //representatives of homology groups and possible torsion coeffcients
+ void Quotient(int dim);
+
+ // transpose the boundary operator matrices, these are boundary operator
+ // matrices for the dual mesh
+ void transposeHMatrices() {
+ for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
+ void transposeHMatrix(int dim) {
+ if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
+
+ // Compute bases for the homology groups of this chain complex
+ void computeHomology(bool dual=false);
+
+ // get coefficient vector for dim-dimensional chain chainNumber
+ // (columns of _Hbasis[dim])
+ std::vector<int> getCoeffVector(int dim, int chainNumber);
+ // torsion coefficient for dim-dimensional chain chainNumber
+ int getTorsion(int dim, int chainNumber);
+ // get rank of homology group of dimension dim
+ int getBasisSize(int dim) {
+ if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]);
+ else return 0; }
+
+ int printMatrix(gmp_matrix* matrix){
+ printf("%d rows and %d columns\n",
+ (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
+ return gmp_matrix_printf(matrix); }
+
+ // debugging aid
+ void matrixTest();
};
// A class representing a chain.
// Used to visualize generators of the homology groups.
class Chain{
-
- private:
- // cells and their coefficients in this chain
- std::map< Cell*, int, Less_Cell > _cells;
- // name of the chain (optional)
- std::string _name;
- // physical group number of the chain
- int _num;
- // cell complex this chain belongs to
- CellComplex* _cellComplex;
- GModel* _model;
-
- // torsion coefficient
- int _torsion;
-
- int _dim;
-
-
- bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*, int, Less_Cell> &cellsNotInChain);
- bool deformChain(std::pair<Cell*, int> cell, bool bend);
-
-
- public:
- Chain(){}
- Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, GModel* model,
- 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();
- _model = chain->getGModel();
- }
- ~Chain(){}
-
- typedef std::map<Cell*, int, Less_Cell>::iterator citer;
-
- // remove a cell from this chain
- void removeCell(Cell* cell);
-
- // add a cell to this chain
- void addCell(Cell* cell, int coeff);
+
+ private:
+ // cells and their coefficients in this chain
+ std::map< Cell*, int, Less_Cell > _cells;
+ // name of the chain (optional)
+ std::string _name;
+ // physical group number of the chain
+ int _num;
+ // cell complex this chain belongs to
+ CellComplex* _cellComplex;
+ GModel* _model;
+
+ // torsion coefficient
+ int _torsion;
+
+ int _dim;
+
+ bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
+ std::map<Cell*, int, Less_Cell> &cellsNotInChain);
+ bool deformChain(std::pair<Cell*, int> cell, bool bend);
+
+
+ public:
+ Chain(){}
+ Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
+ CellComplex* cellComplex, GModel* model,
+ 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();
+ _model = chain->getGModel();
+ }
+ ~Chain(){}
+
+ typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+
+ // remove a cell from this chain
+ void removeCell(Cell* cell);
+
+ // add a cell to this chain
+ void addCell(Cell* cell, int coeff);
+
+ bool hasCell(Cell* c);
+ Cell* findCell(Cell* c);
+ int getCoeff(Cell* c);
+
+
+ int getTorsion() {return _torsion;}
+ int getDim() {return _dim;}
+ CellComplex* getCellComplex() {return _cellComplex;}
+ GModel* getGModel() {return _model;}
+ std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
+
+ // erase cells from the chain with zero coefficient
+ void eraseNullCells();
+
+ void deImmuneCells();
+
+ // number of cells in this chain
+ int getSize() { return _cells.size();}
+
+ // get/set chain name
+ std::string getName() { return _name; }
+ void setName(std::string name) { _name=name; }
+ // get/set physical group number
+ int getNum() { return _num; }
+ void setNum(int num) { _num=num; }
+
+ // make local deformations to the chain to make it smoother and smaller
+ // (for primary complex chains only, not for dual chains represented
+ // by primary cells (yet).)
+ void smoothenChain();
+
+ // append this chain to a MSH ASCII file as $ElementData
+ // for debugging only
+ int writeChainMSH(const std::string &name);
- bool hasCell(Cell* c);
- Cell* findCell(Cell* c);
- int getCoeff(Cell* c);
-
-
- int getTorsion() {return _torsion;}
- int getDim() {return _dim;}
- CellComplex* getCellComplex() {return _cellComplex;}
- GModel* getGModel() {return _model;}
- std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
-
- // erase cells from the chain with zero coefficient
- void eraseNullCells();
-
- void deImmuneCells();
-
- // number of cells in this chain
- int getSize() { return _cells.size();}
-
- // get/set chain name
- std::string getName() { return _name; }
- void setName(std::string name) { _name=name; }
- // get/set physical group number
- int getNum() { return _num; }
- void setNum(int num) { _num=num; }
-
- // make local deformations to the chain to make it smoother and smaller
- // (for primary complex chains only, not for dual chains represented by primary cells (yet).)
- void smoothenChain();
-
- // append this chain to a MSH ASCII file as $ElementData
- // for debugging only
- int writeChainMSH(const std::string &name);
-
- // create a PView of this chain.
- void createPView();
-
+ // create a PView of this chain.
+ void createPView();
+
};
-
+
#endif
#endif
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 1a5344f3ea95c33466c064feb59943c8873bf325..7736b9ee8ce0751a21db3a37e36df46aa4372f82 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -74,16 +74,16 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
}
- Homology::~Homology(){
- delete _cellComplex;
- for(int i = 0; i < 4; i++) {
- for(int j = 0; j < _generators[i].size(); j++){
- Chain* chain = _generators[i].at(j);
- //_model->deletePhysicalGroup(chain->getDim(), chain->getNum());
- delete chain;
- }
- }
- }
+Homology::~Homology(){
+ delete _cellComplex;
+ for(int i = 0; i < 4; i++) {
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ //_model->deletePhysicalGroup(chain->getDim(), chain->getNum());
+ delete chain;
+ }
+ }
+}
void Homology::findGenerators(std::string fileName){
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 0443ba40fdcf4f486572dc36440350953100b8b0..07e1bcd8428301452681a31120b3f8280a4c3613 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -15,52 +15,55 @@
#if defined(HAVE_KBIPACK)
template <class TTypeA, class TTypeB>
-bool convert(const TTypeA& input, TTypeB& output ){
- std::stringstream stream;
- stream << input;
- stream >> output;
- return stream.good();
+ bool convert(const TTypeA& input, TTypeB& output ){
+ std::stringstream stream;
+ stream << input;
+ stream >> output;
+ return stream.good();
}
class Homology
{
- private:
-
- // base cell complex and the model of the homology computation
- CellComplex* _cellComplex;
- GModel* _model;
-
- // domain and the relative subdomain of the homology computation
- std::vector<int> _domain;
- std::vector<int> _subdomain;
+ private:
+
+ // base cell complex and the model of the homology computation
+ CellComplex* _cellComplex;
+ GModel* _model;
- // generator chains
- std::vector<Chain*> _generators[4];
-
- public:
-
- Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
- ~Homology();
-
- // Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
- void findGenerators(std::string fileName);
- void findDualGenerators(std::string fileName);
- void computeBettiNumbers();
-
- bool swapSubdomain() { return _cellComplex->swapSubdomain(); }
-
- // Restore the cell complex to its original state before cell reductions
- void restoreHomology();
-
- // Create a string describing the generator
- std::string getDomainString();
+ // domain and the relative subdomain of the homology computation
+ std::vector<int> _domain;
+ std::vector<int> _subdomain;
+
+ // generator chains
+ std::vector<Chain*> _generators[4];
- // create PViews of the generators and save the chain mesh elements to the mesh of the model
- void createPViews();
- // write the generators to a file
- bool writeGeneratorsMSH(std::string fileName, bool binary=false);
+ public:
+
+ Homology(GModel* model, std::vector<int> physicalDomain,
+ std::vector<int> physicalSubdomain);
+ ~Homology();
+
+ // Find the generators/duals of homology spaces,
+ // or just compute the ranks of homology spaces
+ void findGenerators(std::string fileName);
+ void findDualGenerators(std::string fileName);
+ void computeBettiNumbers();
+
+ bool swapSubdomain() { return _cellComplex->swapSubdomain(); }
+
+ // Restore the cell complex to its original state before cell reductions
+ void restoreHomology();
+ // Create a string describing the generator
+ std::string getDomainString();
+
+ // create PViews of the generators and save the chain mesh elements
+ // to the mesh of the model
+ void createPViews();
+ // write the generators to a file
+ bool writeGeneratorsMSH(std::string fileName, bool binary=false);
+
};
#endif
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index acf2812edc22ede0319feb38dbe26fd36df83e89..cba560a26e69485ebc62ed51518f14738f788c5a 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -29,19 +29,19 @@ int main(int argc, char **argv)
std::vector<int> domain;
std::vector<int> subdomain;
- Homology* homology = new Homology(m, domain, subdomain);
- std::string fileName = "homology_generators.msh";
+ Homology* homology1 = new Homology(m, domain, subdomain);
+ std::string fileName = "homology.msh";
homology->findGenerators(fileName);
- homology->restoreHomology();
+ delete homology1;
- Homology* homology = new Homology(m, domain, subdomain);
- fileName = "homology_dualgenerators.msh";
+ Homology* homology2 = new Homology(m, domain, subdomain);
+ fileName = "homology.msh";
homology->findDualGenerators(fileName);
- homology->restoreHomology();
+ delete homology2;
- Homology* homology = new Homology(m, domain, subdomain);
+ Homology* homology3 = new Homology(m, domain, subdomain);
homology->computeBettiNumbers();
- delete homology;
+ delete homology3;
delete m;
GmshFinalize();