From 61a9dafb349ae154ec5c22d474ea4ea2c6e843ca Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Fri, 15 Jan 2010 11:12:55 +0000
Subject: [PATCH] Code cleanup, increased verbosity.
---
Geo/CMakeLists.txt | 2 +-
Geo/CellComplex.cpp | 252 +++++++++++-------------
Geo/CellComplex.h | 447 +++++--------------------------------------
Geo/ChainComplex.cpp | 79 +++++++-
Geo/ChainComplex.h | 80 +-------
Geo/Homology.cpp | 124 ++++++++++--
Geo/Homology.h | 73 +------
tutorial/t10.geo | 5 +-
8 files changed, 354 insertions(+), 708 deletions(-)
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 89e0c22ee7..a004d0e60a 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -31,7 +31,7 @@ set(SRC
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp
- CellComplex.cpp ChainComplex.cpp Homology.cpp
+ Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Geo *.h)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 7162aefca0..6e67fe2e6a 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -18,7 +18,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_dim = 0;
_simplicial = true;
- // find boundary elements
+ // find boundary entities
find_boundary(domain, subdomain);
// insert cells into cell complex
@@ -182,7 +182,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
cell = new CPyramid(vertices, tag, subdomain, boundary);
_simplicial = false;
}*/
- else printf("Error: mesh element %d not implemented yet! \n", type);
+ else Msg::Error("Error: mesh element %d not implemented yet! \n", type);
tag++;
//if(domain.at(j)->getMeshElement(i)->getNum() == 8196) cell->setImmune(true);
//else
@@ -203,7 +203,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
if(dim == 3){
if(vertices.size() == 3) newCell = new TwoSimplex(vertices, tag, subdomain, boundary);
else if(vertices.size() == 4) newCell = new CQuadrangle(vertices, tag, subdomain, boundary);
- else printf("Error: invalid face! \n");
+ else Msg::Debug("Error: invalid face! \n");
}
else if(dim == 2) newCell = new OneSimplex(vertices, tag, subdomain, boundary);
else if(dim == 1) newCell = new ZeroSimplex(vertices, tag, subdomain, boundary);
@@ -236,98 +236,54 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
+CellComplex::~CellComplex(){
+ for(int i = 0; i < 4; i++){
+ _cells[i].clear();
+
+ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
+ _ocells[i].clear();
+
+ }
+ //emptyTrash();
+ }
+
+ /*
+ 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();
+ }
+ */
+
void CellComplex::insertCell(Cell* cell){
_ecells.insert(cell);
}
-/*
-int Simplex::kappa(Cell* tau) const{
- for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
- }
-
- if(this->getDim() - tau->getDim() != 1) return 0;
-
- int value=1;
- for(int i=0; i < tau->getNumVertices(); i++){
- if(this->getSortedVertex(i) != tau->getSortedVertex(i)) return value;
- value = value*-1;
- }
-
- return value;
-}
-*/
-/*
-int Cell::kappa(Cell* tau) const{
- for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
- }
-
- if(this->getDim() - tau->getDim() != 1) return 0;
-
- int value=1;
-
- for(int i = 0; i < this->getNumFacets(); i++){
- std::vector<MVertex*> vTau = tau->getVertexVector();
- std::vector<MVertex*> v;
- this->getFacetVertices(i, v);
- value = -1;
-
- if(v.size() != vTau.size()) printf("Error: invalid facet!");
-
- do {
- value = value*-1;
- if(v == vTau) return value;
- }
- while (std::next_permutation(vTau.begin(), vTau.end()) );
-
- vTau = tau->getVertexVector();
- value = -1;
- do {
- value = value*-1;
- if(v == vTau) return value;
- }
- while (std::prev_permutation(vTau.begin(), vTau.end()) );
-
-
- }
-
- return 0;
-}
-
-
-int OneSimplex::kappa(Cell* tau) const{
-
- for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
- }
-
- if(tau->getDim() != 0) return 0;
-
- if(tau->getVertex(0) == this->getVertex(0)) return -1;
- else return 1;
-
-}*/
-
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::list<Cell*> boundary = cell->getBoundary();
- //std::list<Cell*> coboundary = cell->getCoboundary();
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).first;
- //Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell, other);
- }
-
+ }
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).first;
- //Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell, other);
}
@@ -527,7 +483,7 @@ int CellComplex::reduceComplex(){
double t1 = Cpu();
- printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+ Msg::Debug("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
@@ -535,23 +491,15 @@ int CellComplex::reduceComplex(){
int omitted = 0;
_store.clear();
- //if(omit > getDim()) omit = getDim();
-
- CellComplex::removeSubdomain();
- //std::set<Cell*, Less_Cell> generatorCells;
+ removeSubdomain();
while (getSize(getDim()) != 0){
citer cit = firstCell(getDim());
-
Cell* cell = *cit;
- //generatorCells.insert(cell);
removeCell(cell);
- //makeDualComplex();
- //coreduction(cell);
- //makeDualComplex();
omitted++;
std::set< Cell*, Less_Cell > omittedCells;
_store.push_back(omittedCells);
@@ -561,17 +509,9 @@ int CellComplex::reduceComplex(){
}
- /*
- for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
- Cell* cell = *cit;
- cell->clearBoundary();
- cell->clearCoboundary();
- _cells[cell->getDim()].insert(cell);
- }*/
-
double t2 = Cpu();
- printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
+ Msg::Debug("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
/*
@@ -611,7 +551,7 @@ void CellComplex::removeSubdomain(){
int CellComplex::coreduceComplex(){
- printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+ Msg::Debug("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
@@ -631,18 +571,18 @@ int CellComplex::coreduceComplex(){
int omitted = 0;
_store.clear();
- //std::set<Cell*, Less_Cell> generatorCells;
while (getSize(0) != 0){
citer cit = firstCell(0);
Cell* cell = *cit;
- //generatorCells.insert(cell);
+
removeCell(cell, false);
std::set< Cell*, Less_Cell > omittedCells;
omitted++;
_store.push_back(omittedCells);
_store.at(omitted-1).insert(cell);
coreduction(cell, omitted);
- }
+ }
+
/*
for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
Cell* cell = *cit;
@@ -653,7 +593,7 @@ int CellComplex::coreduceComplex(){
*/
- printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+ Msg::Debug("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
return omitted;
@@ -663,7 +603,7 @@ void CellComplex::computeBettiNumbers(){
for(int i = 0; i < 4; i++){
_betti[i] = 0;
- printf("Betti number computation process: step %d of 4 \n", i+1);
+ Msg::Debug("Betti number computation process: step %d of 4 \n", i+1);
while (getSize(i) != 0){
citer cit = firstCell(i);
@@ -677,7 +617,7 @@ void CellComplex::computeBettiNumbers(){
coreduction(cell);
}
}
- printf("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n",
+ Msg::Debug("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n",
getBettiNumber(0), getBettiNumber(1), getBettiNumber(2), getBettiNumber(3));
return;
@@ -687,7 +627,7 @@ void CellComplex::computeBettiNumbers(){
int CellComplex::cocombine(int dim){
double t1 = Cpu();
- printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
+ Msg::Debug("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
if(dim < 0 || dim > 2) return 0;
@@ -735,7 +675,7 @@ int CellComplex::cocombine(int dim){
removeCell(c1);
removeCell(c2);
std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
- if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n");
+ if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
cit = firstCell(dim);
count++;
@@ -746,7 +686,7 @@ int CellComplex::cocombine(int dim){
}
}
double t2 = Cpu();
- printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
+ Msg::Debug("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
getSize(3), getSize(2), getSize(1), getSize(0), t2-t1);
return count;
@@ -754,7 +694,7 @@ int CellComplex::cocombine(int dim){
int CellComplex::combine(int dim){
double t1 = Cpu();
- printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
+ Msg::Debug("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
if(dim < 1 || dim > 3) return 0;
@@ -802,7 +742,7 @@ int CellComplex::combine(int dim){
removeCell(c1);
removeCell(c2);
std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
- if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n");
+ if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
cit = firstCell(dim);
count++;
}
@@ -812,7 +752,7 @@ int CellComplex::combine(int dim){
}
}
double t2 = Cpu();
- printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
+ Msg::Debug("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
@@ -820,7 +760,7 @@ int CellComplex::combine(int dim){
}
/*
-int CellComplex::combine(int dim){
+int CellComplex::combine(int dim){ // version of combine that doesn't have a queue
double t1 = Cpu();
printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -900,22 +840,17 @@ void CellComplex::swapSubdomain(){
int CellComplex::writeComplexMSH(const std::string &name){
-
FILE *fp = fopen(name.c_str(), "w");
if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str());
- printf("Unable to open file.");
+ Msg::Debug("Unable to open file.");
return 0;
- }
-
-
+ }
fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
-
fprintf(fp, "$Nodes\n");
-
fprintf(fp, "%d\n", (int)_domainVertices.size());
for(std::set<MVertex*, Less_MVertex>::iterator vit = _domainVertices.begin(); vit != _domainVertices.end(); vit++){
@@ -1073,14 +1008,14 @@ bool CellComplex::checkCoherence(){
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");
+ Msg::Debug("Warning! Boundary cell not in cell complex! Boundary removed. \n");
//printf(" "); cell->printCell();
//printf(" "); bdCell->printCell();
cell->removeBoundaryCell(bdCell);
coherent = false;
}
if(!bdCell->hasCoboundary(cell)){
- printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
+ Msg::Debug("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
//printf(" "); cell->printCell();
//printf(" "); cell->printBoundary();
//printf(" "); bdCell->printCell();
@@ -1096,14 +1031,14 @@ bool CellComplex::checkCoherence(){
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");
+ Msg::Debug("Warning! Coboundary cell not in cell complex! Coboundary removed. \n");
//printf(" "); cell->printCell();
//printf(" "); cbdCell->printCell();
cell->removeCoboundaryCell(cbdCell);
coherent = false;
}
if(!cbdCell->hasBoundary(cell)){
- printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
+ Msg::Debug("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
//printf(" "); cell->printCell();
//printf(" "); cell->printCoboundary();
//printf(" "); cbdCell->printCell();
@@ -1119,24 +1054,55 @@ bool CellComplex::checkCoherence(){
return coherent;
}
+ void CellComplex::emptyTrash(){
+ for(std::list<Cell*>::iterator cit = _trash.begin(); cit != _trash.end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
+ _trash.clear();
+ }
+
+ void CellComplex::restoreComplex(){
+ for(int i = 0; i < 4; i++){
+ _betti[i] = 0;
+ _cells[i] = _ocells[i];
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ cell->restoreCell();
+ }
+ }
+ _store.clear();
+ _ecells.clear();
+ _trash.clear();
+ }
+
+ bool CellComplex::hasCell(Cell* cell, bool org){
+ 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;
+ }
+ }
+
+ void CellComplex::makeDualComplex(){
+ std::set<Cell*, Less_Cell> temp = _cells[0];
+ _cells[0] = _cells[3];
+ _cells[3] = temp;
+ temp = _cells[1];
+ _cells[1] = _cells[2];
+ _cells[2] = temp;
+
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ cell->makeDualCell();
+ }
+ }
+ }
-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 92ac46b194..3ae5943c43 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -109,16 +109,7 @@ class Cell
virtual std::vector<MVertex*> getVertexVector() const = 0;
// restores the cell information to its original state before reduction
- virtual void restoreCell(){
- _boundary = _obd;
- _coboundary = _ocbd;
- _bdSize = _boundary.size();
- _cbdSize = _coboundary.size();
- _combined = false;
- _index = 0;
- _immune = false;
-
- }
+ virtual void restoreCell();
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const = 0;
@@ -126,165 +117,43 @@ class Cell
// (co)boundary cell iterator
typedef std::map<Cell*, int, Less_Cell>::iterator biter;
- //virtual int getBoundarySize() { return _bdSize; }
- //virtual int getCoboundarySize() { return _cbdSize; }
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() {
- std::list<Cell*> boundary;
- for( biter it= _boundary.begin(); it != _boundary.end();it++){
- Cell* cell = (*it).first;
- boundary.push_back(cell);
- }
- return boundary;
- }
+ virtual std::list< Cell* > getBoundary();
virtual std::map<Cell*, int, Less_Cell > getOrientedCoboundary() { return _coboundary; }
- virtual std::list< 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;
- }
-
- 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 ) );
- return true;
- }
-
- virtual bool addCoboundaryCell(int orientation, Cell* cell) {
- _cbdSize++;
- 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.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;
- }
-
- return 0;
- }
+ virtual std::list< Cell* > getCoboundary();
+ // add (co)boundary cell
+ virtual bool addBoundaryCell(int orientation, Cell* cell);
+ virtual bool addCoboundaryCell(int orientation, Cell* cell);
- virtual int removeCoboundaryCell(Cell* cell, bool other=true) {
- 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;
- }
-
- // true if has given cell on (original) boundary
- virtual bool hasBoundary(Cell* cell, bool org=false){
- if(!org){
- 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;
- }
- }
-
- // true if has given cell on (original) boundary
- virtual bool hasCoboundary(Cell* cell, bool org=false){
- if(!org){
- biter it = _coboundary.find(cell);
- if(it != _coboundary.end()) return true;
- return false;
- }
- else{
- biter it = _ocbd.find(cell);
- if(it != _ocbd.end()) return true;
- return false;
- }
- }
+ // 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(){
- std::map<Cell*, int, Less_Cell > temp = _boundary;
- _boundary = _coboundary;
- _coboundary = temp;
- _dim = 3-_dim;
-
- }
+ virtual void makeDualCell();
- virtual void printBoundary() {
- 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(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 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() {
- 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 void printOrgBd();
+ virtual void printOrgCbd();
virtual bool inSubdomain() const { return _inSubdomain; }
//virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
@@ -297,15 +166,16 @@ class Cell
virtual int getNumFacets() const { return 0; }
virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const {};
+
// get boundary cell orientation
virtual int getFacetOri(Cell* cell) { std::vector<MVertex*> v = cell->getVertexVector(); return getFacetOri(v); }
virtual int getFacetOri(std::vector<MVertex*> &v) { return 0; }
+ // 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;}
-
-
+
bool operator==(const Cell& c2){
if(this->getNumVertices() != c2.getNumVertices()){
return false;
@@ -330,21 +200,16 @@ class Cell
return false;
}
- // checks whether lower dimensional simplex tau is on the boundary of this cell
- //virtual int kappa(Cell* tau) const;
-
- virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, int num=0, int elementary=1, int physical=1, int parentNum=0) {}
+ //virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, int num=0, int elementary=1, int physical=1, int parentNum=0) {}
};
// Simplicial cell type.
class Simplex : public Cell
-{
-
+{
public:
Simplex() : Cell() {}
- ~Simplex(){}
-
+ ~Simplex(){}
};
// Zero simplex cell type.
@@ -385,9 +250,7 @@ class OneSimplex : public Simplex, public MLine
private:
// sorted list of vertices
int _vs[2];
- public:
-
-
+ public:
OneSimplex(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
: MLine(v, num, part){
_tag = tag;
@@ -427,8 +290,6 @@ class OneSimplex : public Simplex, public MLine
else if(v[0] == _v[1]) return 1;
else return 0;
}
-
- //int kappa(Cell* tau) const;
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
@@ -798,105 +659,8 @@ class CombinedCell : public Cell{
public:
- CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
-
- // use "smaller" cell as c2
- if(c1->getNumVertices() < c2->getNumVertices()){
- Cell* temp = c1;
- c1 = c2;
- c2 = temp;
- }
-
- _index = c1->getIndex();
- _tag = c1->getTag();
- _dim = c1->getDim();
- _num = c1->getNum();
- _inSubdomain = c1->inSubdomain();
- _onDomainBoundary = c1->onDomainBoundary();
- _combined = true;
-
- _v.reserve(c1->getNumVertices() + c2->getNumVertices());
- _vs.reserve(c1->getNumVertices() + c2->getNumVertices());
-
-
- _v = c1->getVertexVector();
- for(int i = 0; i < c2->getNumVertices(); 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());
- 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++){
- 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();
-
- 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++){
- Cell* cell = (*it).first;
- if(!orMatch) (*it).second = -1*(*it).second;
- int ori = (*it).second;
- cell->removeCoboundaryCell(c2);
- if(co){
- 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); }
- }
- else{
- if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
- }
- }
-
- // coboundary cells
- std::map<Cell*, int, Less_Cell > c1Coboundary = c1->getOrientedCoboundary();
- std::map<Cell*, int, Less_Cell > c2Coboundary = c2->getOrientedCoboundary();
-
- 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++){
- Cell* cell = (*it).first;
- if(!orMatch) (*it).second = -1*(*it).second;
- int ori = (*it).second;
- cell->removeBoundaryCell(c2);
- if(!co){
- 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); }
- }
- else {
- if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
- }
- }
-
- }
-
- ~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;
- }
- }
+ CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
+ ~CombinedCell();
int getNum() { return _num; }
int getNumVertices() const { return _v.size(); }
@@ -904,27 +668,10 @@ class CombinedCell : public Cell{
int getSortedVertex(int vertex) const { return _vs.at(vertex); }
std::vector<MVertex*> getVertexVector() const { return _v; }
- //int kappa(Cell* tau) const { return 0; }
-
// true if this cell has given vertex
- bool hasVertex(int vertex) const {
- //std::vector<int>::iterator it = std::find(_vs.begin(), _vs.end(), vertex);
- //if (it == _vs.end()) return false;
- //else return true;
- for(unsigned int i = 0; i < _v.size(); i++){
- if(_v.at(i)->getNum() == vertex) return true;
- }
- return false;
- }
+ bool hasVertex(int vertex) const;
- virtual void printCell() const {
- printf("Cell dimension: %d, ", getDim() );
- printf("Vertices: ");
- for(int i = 0; i < this->getNumVertices(); i++){
- printf("%d ", this->getSortedVertex(i));
- }
- printf(", in subdomain: %d\n", _inSubdomain);
- }
+ virtual void printCell() const;
std::list< std::pair<int, Cell*> > getCells() { return _cells; }
int getNumCells() {return _cells.size();}
@@ -952,8 +699,7 @@ class CellComplex
std::set<Cell*, Less_Cell> _cells[4];
// storage for cell pointers to delete them
- std::list<Cell*> _trash;
-
+ std::list<Cell*> _trash;
std::vector< std::set<Cell*, Less_Cell> > _store;
std::set<Cell*, Less_Cell> _ecells;
@@ -969,11 +715,6 @@ class CellComplex
// Is the cell complex simplicial
bool _simplicial;
- public:
- // iterator for the cells of same dimension
- typedef std::set<Cell*, Less_Cell>::iterator citer;
-
- private:
// enqueue cells in queue if they are not there already
void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
@@ -987,76 +728,21 @@ class CellComplex
// remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
- public:
- void insertCell(Cell* cell);
- // reduction of this cell complex
- // removes reduction pairs of cell of dimension dim and dim-1
- int reduction(int dim, int omitted=0);
- private:
// queued coreduction presented in Mrozek's paper
int coreduction(Cell* generator, int omitted=0);
-
public:
- /*
- 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( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
+ //CellComplex(CellComplex* cellComplex)
CellComplex(){}
- ~CellComplex(){
- for(int i = 0; i < 4; i++){
- _cells[i].clear();
-
- for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
- Cell* cell = *cit;
- delete cell;
- }
- _ocells[i].clear();
-
- }
- //emptyTrash();
- }
+ ~CellComplex();
-
- void emptyTrash(){
- for(std::list<Cell*>::iterator cit = _trash.begin(); cit != _trash.end(); cit++){
- Cell* cell = *cit;
- delete cell;
- }
- _trash.clear();
- }
+ void emptyTrash();
// restore this cell complex to its original state
- void restoreComplex(){
- for(int i = 0; i < 4; i++){
- _betti[i] = 0;
- _cells[i] = _ocells[i];
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->restoreCell();
- }
- }
- _store.clear();
- _ecells.clear();
- _trash.clear();
- }
+ void restoreComplex();
// get the number of certain dimensional cells
int getSize(int dim){ return _cells[dim].size(); }
@@ -1069,6 +755,9 @@ class CellComplex
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(); }
@@ -1077,49 +766,36 @@ class CellComplex
// true if cell complex has given cell
- 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;
- }
- }
-
-
- // kappa for two cells of this cell complex
- // implementation will vary depending on cell type
- //inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
-
+ 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()) ); }
-
+ void insertCell(Cell* cell);
+
+ // 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);
- // write this cell complex in 2.1 MSH ASCII format
+ // write this cell complex in 2.0 MSH ASCII format
+ // for debugging only
int writeComplexMSH(const std::string &name);
// Cell combining for reduction and coreduction
int combine(int dim);
int cocombine(int dim);
-
void computeBettiNumbers();
int getBettiNumber(int i) { if(i > -1 && i < 4) return _betti[i]; else return 0; }
@@ -1127,36 +803,15 @@ class CellComplex
// also check whether all (co)boundary cells of a cell belong to this cell complex
bool checkCoherence();
- int eulerCharacteristic(){
- return getSize(0) - getSize(1) + getSize(2) - getSize(3);
- }
+ 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); }
-
+ 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(){
- std::set<Cell*, Less_Cell> temp = _cells[0];
- _cells[0] = _cells[3];
- _cells[3] = temp;
- temp = _cells[1];
- _cells[1] = _cells[2];
- _cells[2] = temp;
-
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->makeDualCell();
- }
- }
- }
-
-
-
-
+ void makeDualComplex();
};
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 172fe594e7..7509908513 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -82,17 +82,15 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
//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){
- printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
+ Msg::Debug("Warning: Index out of bound! HMatrix: %d. \n", dim);
}
else{
gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).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);
- }*/
+ if( abs((old_elem + (*it).second)) > 1){
+ Msg::Debug("Incidence index: %d! HMatrix: %d.", (old_elem + (*it).second), dim);
+ }
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
}
@@ -291,7 +289,7 @@ void ChainComplex::computeHomology(bool dual){
//KerCod(highDim);
}
- printf("Homology computation process: step %d of 4 \n", i+1);
+ Msg::Debug("Homology computation process: step %d of 4 \n", i+1);
KerCod(highDim);
@@ -560,7 +558,7 @@ void Chain::smoothenChain(){
eraseNullCells();
double t2 = Cpu();
- printf("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
+ Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
return;
}
@@ -574,8 +572,8 @@ int Chain::writeChainMSH(const std::string &name){
FILE *fp = fopen(name.c_str(), "a");
if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str());
- printf("Unable to open file.");
- return 0;
+ Msg::Debug("Unable to open file.");
+ return 0;
}
fprintf(fp, "\n$ElementData\n");
@@ -674,4 +672,65 @@ void Chain::createPView(){
}
+void Chain::removeCell(Cell* cell) {
+ citer it = _cells.find(cell);
+ if(it != _cells.end()){
+ (*it).second = 0;
+ }
+ return;
+}
+
+void Chain::addCell(Cell* cell, int coeff) {
+ std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
+ if(!insert.second && (*insert.first).second == 0) (*insert.first).second = coeff;
+ else if (!insert.second && (*insert.first).second != 0) Msg::Debug("Error: invalid chain smoothening add! \n");
+
+ if(!_cellComplex->hasCell(cell)){
+ _cellComplex->insertCell(cell);
+ }
+ return;
+}
+
+bool Chain::hasCell(Cell* c){
+ citer it = _cells.find(c);
+ if(it != _cells.end() && (*it).second != 0) return true;
+ return false;
+}
+Cell* Chain::findCell(Cell* c){
+ citer it = _cells.find(c);
+ if(it != _cells.end() && (*it).second != 0) return (*it).first;
+ return NULL;
+}
+int Chain::getCoeff(Cell* c){
+ citer it = _cells.find(c);
+ if(it != _cells.end()) return (*it).second;
+ return 0;
+}
+
+void Chain::eraseNullCells(){
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if( (*cit).second == 0){
+ //cit++;
+ //_cells.erase(--cit);
+ _cells.erase(cit);
+ ++cit;
+ }
+ }
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if( (*cit).second == 0){
+ _cells.erase(cit);
+ cit = _cells.begin();
+ }
+ }
+ return;
+}
+
+
+void Chain::deImmuneCells(){
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ cell->setImmune(false);
+ }
+}
+
#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 3a474b3a2c..3897b6cc82 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -115,7 +115,7 @@ class ChainComplex{
// Compute bases for the homology groups of this chain complex
void computeHomology(bool dual=false);
- // get coefficient vector for dim-dimensional chain chainNumber
+ // 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);
@@ -171,48 +171,15 @@ class 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; }
- // get coeffcient of i:th cell of this chain
- //int getCoeff(int i) { return _cells.at(i).second; }
-
-
// remove a cell from this chain
- void removeCell(Cell* cell) {
- citer it = _cells.find(cell);
- if(it != _cells.end()){
- (*it).second = 0;
- }
- return;
- }
+ void removeCell(Cell* cell);
// add a cell to this chain
- void addCell(Cell* cell, int coeff) {
- std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
- if(!insert.second && (*insert.first).second == 0) (*insert.first).second = coeff;
- else if (!insert.second && (*insert.first).second != 0) printf("Error: invalid chain smoothening add! \n");
-
- if(!_cellComplex->hasCell(cell)){
- _cellComplex->insertCell(cell);
- }
- return;
- }
+ void addCell(Cell* cell, int coeff);
- 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;
- }
+ bool hasCell(Cell* c);
+ Cell* findCell(Cell* c);
+ int getCoeff(Cell* c);
int getTorsion() {return _torsion;}
@@ -222,40 +189,12 @@ class Chain{
std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
// erase cells from the chain with zero coefficient
- void eraseNullCells(){
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if( (*cit).second == 0){
- //cit++;
- //_cells.erase(--cit);
- _cells.erase(cit);
- ++cit;
- }
- }
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if( (*cit).second == 0){
- _cells.erase(cit);
- cit = _cells.begin();
- }
- }
- return;
- }
+ void eraseNullCells();
-
- void deImmuneCells(){
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- Cell* cell = (*cit).first;
- cell->setImmune(false);
- }
- }
+ void deImmuneCells();
// number of cells in this chain
- int getSize() {
- //eraseNullCells();
- return _cells.size();
- }
- int getNumCells() {
- return _cells.size();
- }
+ int getSize() { return _cells.size();}
// get/set chain name
std::string getName() { return _name; }
@@ -269,6 +208,7 @@ class Chain{
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.
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 9d83783a41..7b0a96c0ab 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -19,6 +19,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
_subdomain = physicalSubdomain;
Msg::Info("Creating a Cell Complex...");
+ Msg::StatusBar(1, false, "Cell Complex...");
+ Msg::StatusBar(2, false, "");
double t1 = Cpu();
std::map<int, std::vector<GEntity*> > groups[4];
@@ -68,13 +70,25 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
Msg::Info("Cell Complex complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
+ Msg::StatusBar(2, false, "%d V, %d F, %d E, %d V.",
+ _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;
+ }
+ }
+ }
void Homology::findGenerators(std::string fileName){
- Msg::Info("Reducing Cell Complex...");
+ Msg::Info("Reducing the Cell Complex...");
+ Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
//_cellComplex->printEuler();
int omitted = _cellComplex->reduceComplex();
@@ -98,10 +112,12 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
-
+ Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
+ _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
- Msg::Info("Computing homology groups...");
+ Msg::Info("Computing homology spaces...");
+ Msg::StatusBar(1, false, "Computing...");
t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex);
chains->computeHomology();
@@ -149,19 +165,22 @@ void Homology::findGenerators(std::string fileName){
createPViews();
if(fileName != "") writeGeneratorsMSH(fileName);
- Msg::Info("Ranks of homology groups for primal cell complex:");
+ Msg::Info("Ranks of homology spaces for primal cell complex:");
Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]);
Msg::Info("H3 = %d", HRank[3]);
- if(omitted != 0) Msg::Info("The computation of generators in %d highest dimensions was omitted.", omitted);
+ if(omitted != 0) Msg::Info("The computation of generators in the highest dimension was omitted.");
delete chains;
- printf("H0 = %d \n", HRank[0]);
- printf("H1 = %d \n", HRank[1]);
- printf("H2 = %d \n", HRank[2]);
- printf("H3 = %d \n", HRank[3]);
+ Msg::Debug("H0 = %d \n", HRank[0]);
+ Msg::Debug("H1 = %d \n", HRank[1]);
+ Msg::Debug("H2 = %d \n", HRank[2]);
+ Msg::Debug("H3 = %d \n", HRank[3]);
+
+ Msg::StatusBar(1, false, "Homology");
+ Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", HRank[0], HRank[1], HRank[2], HRank[3]);
return;
}
@@ -171,6 +190,7 @@ void Homology::findDualGenerators(std::string fileName){
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
Msg::Info("Reducing Cell Complex...");
+ Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex();
//_cellComplex->emptyTrash();
@@ -199,12 +219,14 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
-
+ Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
+ _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
- _cellComplex->writeComplexMSH(fileName);
+ //_cellComplex->writeComplexMSH(fileName);
- Msg::Info("Computing homology groups...");
+ Msg::Info("Computing homology spaces...");
+ Msg::StatusBar(1, false, "Computing...");
t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex);
chains->transposeHMatrices();
@@ -255,7 +277,7 @@ void Homology::findDualGenerators(std::string fileName){
createPViews();
if(fileName != "") writeGeneratorsMSH(fileName);
- Msg::Info("Ranks of homology groups for dual cell complex:");
+ Msg::Info("Ranks of homology spaces for dual cell complex:");
Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]);
Msg::Info("H2* = %d", HRank[2]);
@@ -264,10 +286,13 @@ void Homology::findDualGenerators(std::string fileName){
delete chains;
- printf("H0* = %d \n", HRank[0]);
- printf("H1* = %d \n", HRank[1]);
- printf("H2* = %d \n", HRank[2]);
- printf("H3* = %d \n", HRank[3]);
+ Msg::Debug("H0* = %d \n", HRank[0]);
+ Msg::Debug("H1* = %d \n", HRank[1]);
+ Msg::Debug("H2* = %d \n", HRank[2]);
+ Msg::Debug("H3* = %d \n", HRank[3]);
+
+ Msg::StatusBar(1, false, "Homology");
+ Msg::StatusBar(2, false, "H0*: %d, H1*: %d, H2*: %d, H3*: %d.", HRank[0], HRank[1], HRank[2], HRank[3]);
return;
}
@@ -275,6 +300,7 @@ void Homology::findDualGenerators(std::string fileName){
void Homology::computeBettiNumbers(){
Msg::Info("Running coreduction...");
+ Msg::StatusBar(1, false, "Computing...");
double t1 = Cpu();
_cellComplex->computeBettiNumbers();
double t2 = Cpu();
@@ -285,7 +311,69 @@ void Homology::computeBettiNumbers(){
Msg::Info("H2 = %d", _cellComplex->getBettiNumber(2));
Msg::Info("H3 = %d", _cellComplex->getBettiNumber(3));
+ Msg::StatusBar(1, false, "Homology");
+ Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", _cellComplex->getBettiNumber(0), _cellComplex->getBettiNumber(1), _cellComplex->getBettiNumber(2), _cellComplex->getBettiNumber(3));
+
return;
}
+
+ void Homology::restoreHomology() {
+ _cellComplex->restoreComplex();
+ for(int i = 0; i < 4; i++) _generators[i].clear();
+ }
+
+ std::string Homology::getDomainString() {
+ std::string domainString = "({";
+ for(unsigned int i = 0; i < _domain.size(); i++){
+ std::string temp = "";
+ convert(_domain.at(i),temp);
+ domainString += temp;
+ if (_domain.size()-1 > i){
+ domainString += ", ";
+ }
+ }
+ domainString += "}";
+
+ if(!_subdomain.empty()){
+ domainString += ", {";
+
+ for(unsigned int i = 0; i < _subdomain.size(); i++){
+ std::string temp = "";
+ convert(_subdomain.at(i),temp);
+ domainString += temp;
+ if (_subdomain.size()-1 > i){
+ domainString += ", ";
+ }
+ }
+ domainString += "}";
+
+ }
+ domainString += ") ";
+ return domainString;
+ }
+
+ void Homology::createPViews(){
+ for(int i = 0; i < 4; i++){
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ chain->createPView();
+ }
+ }
+ }
+
+ bool Homology::writeGeneratorsMSH(std::string fileName, bool binary){
+ if(!_model->writeMSH(fileName, 2.0, binary)) return false;
+ /*
+ for(int i = 0; i < 4; i++){
+ for(int j = 0; j < _generators[i].size(); j++){
+ Chain* chain = _generators[i].at(j);
+ if(!chain->writeChainMSH(fileName)) return false;
+ }
+ }*/
+ Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
+ Msg::Debug("Wrote homology computation results to %s. \n", fileName.c_str());
+
+ return true;
+ }
#endif
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 99af3629f1..fbc47be58c 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -41,16 +41,7 @@ class Homology
public:
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
- ~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();
// Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
void findGenerators(std::string fileName);
@@ -61,67 +52,15 @@ class Homology
//void swapSubdomain() { _cellComplex->swapSubdomain(); }
// Restore the cell complex to its original state before cell reductions
- void restoreHomology() {
- _cellComplex->restoreComplex();
- for(int i = 0; i < 4; i++) _generators[i].clear();
- }
+ void restoreHomology();
// Create a string describing the generator
- std::string getDomainString() {
- std::string domainString = "({";
- for(unsigned int i = 0; i < _domain.size(); i++){
- std::string temp = "";
- convert(_domain.at(i),temp);
- domainString += temp;
- if (_domain.size()-1 > i){
- domainString += ", ";
- }
- }
- domainString += "}";
-
- if(!_subdomain.empty()){
- domainString += ", {";
-
- for(unsigned int i = 0; i < _subdomain.size(); i++){
- std::string temp = "";
- convert(_subdomain.at(i),temp);
- domainString += temp;
- if (_subdomain.size()-1 > i){
- domainString += ", ";
- }
- }
- domainString += "}";
-
- }
- domainString += ") ";
- return domainString;
- }
-
- // create PViews of the generators
- void createPViews(){
- for(int i = 0; i < 4; i++){
- for(int j = 0; j < _generators[i].size(); j++){
- Chain* chain = _generators[i].at(j);
- chain->createPView();
- }
- }
- }
+ 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){
- if(!_model->writeMSH(fileName, 2.0, binary)) return false;
- /*
- for(int i = 0; i < 4; i++){
- for(int j = 0; j < _generators[i].size(); j++){
- Chain* chain = _generators[i].at(j);
- if(!chain->writeChainMSH(fileName)) return false;
- }
- }*/
- Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
- printf("Wrote homology computation results to %s. \n", fileName.c_str());
-
- return true;
- }
+ bool writeGeneratorsMSH(std::string fileName, bool binary=false);
};
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index 464d688417..a97eb68816 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -133,9 +133,9 @@ Mesh 3;
// Save the generator chains to t10_hom.msh.
HomGen("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
-// Find the corresponding cuts.
+// Find the corresponding thin cuts.
// Save the cut chains to t10_hom.msh.
-HomCut("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
+HomGen("t10_hom.msh") = {{69}, {75}};
// Only find and print the ranks of the relative homology spaces (Betti numbers).
HomRank {{69},{70, 71, 72, 73}};
@@ -143,6 +143,5 @@ HomRank {{69},{70, 71, 72, 73}};
// More examples (uncomment):
// HomGen("t10_homgen.msh_1") = {{69}, {}};
// HomGen("t10_homgen.msh_2") = {{69}, {74}};
-// HomGen("t10_homgen.msh_3") = {{69}, {75}};
--
GitLab