diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 58553fc608f45bd42e1c352814ab6cb00bd4b36c..71b4415876ba5d5e02bb631ddd3c6bc151c48703 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -280,15 +280,10 @@ CombinedCell::~CombinedCell(){
}
}
-// true if this cell has given vertex
bool CombinedCell::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;
+ std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), vertex);
+ if (it != _vs.end()) return true;
+ else return false;
}
void CombinedCell::printCell() const {
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 6e67fe2e6a7a92ee7157f629235589069b89200e..8b49f25cee53fdaa0bf5285a8e4d8e4188807a4b 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -268,10 +268,6 @@ CellComplex::~CellComplex(){
}
*/
-void CellComplex::insertCell(Cell* cell){
- _ecells.insert(cell);
-}
-
void CellComplex::removeCell(Cell* cell, bool other){
std::map<Cell*, int, Less_Cell > coboundary = cell->getOrientedCoboundary();
@@ -337,10 +333,8 @@ int CellComplex::coreduction(Cell* generator, int omitted){
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());
- else _trash.push_back(bd_s.front());
coreductions++;
- _trash.push_back(s);
}
else if(bd_s.empty()){
@@ -375,10 +369,7 @@ int CellComplex::reduction(int dim, int omitted){
++cit;
removeCell(cbd_c.front());
removeCell(cell);
- _trash.push_back(cell);
- if(dim == getDim() && omitted > 0) _store.at(omitted-1).insert(cbd_c.front());
- else _trash.push_back(cbd_c.front());
-
+ if(dim == getDim() && omitted > 0) _store.at(omitted-1).insert(cbd_c.front());
count++;
reduced = true;
@@ -619,7 +610,6 @@ void CellComplex::computeBettiNumbers(){
}
Msg::Debug("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n",
getBettiNumber(0), getBettiNumber(1), getBettiNumber(2), getBettiNumber(3));
-
return;
}
@@ -664,7 +654,6 @@ int CellComplex::cocombine(int dim){
&& c2->getNumVertices() < getSize(dim)){
removeCell(s);
- _trash.push_back(s);
cbd_c = c1->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
@@ -731,7 +720,6 @@ int CellComplex::combine(int dim){
&& c2->getNumVertices() < getSize(dim)){
removeCell(s);
- _trash.push_back(s);
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
@@ -872,7 +860,6 @@ int CellComplex::writeComplexMSH(const std::string &name){
for(int i = 0; i < _store.size(); i++){
count = count + _store.at(i).size();
}
- count = count + _ecells.size();
fprintf(fp, "%d\n", count);
@@ -957,6 +944,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
}
}
+ /*
for(citer cit = _ecells.begin(); cit != _ecells.end(); cit++){
Cell* cell = *cit;
type = cell->getTypeForMSH();
@@ -975,9 +963,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
}
- }
-
-
+ }*/
fprintf(fp, "$EndElements\n");
@@ -1054,55 +1040,45 @@ 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();
- }
- }
- }
+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();
+}
+
+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();
+ }
+ }
+}
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 3ae5943c43567bc32fe739585bb0b73ce1e556d2..ff2bd3b08cc81892497804f68e15cb19341dd649 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -66,8 +66,6 @@ class Cell
// 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;
@@ -689,6 +687,7 @@ class CellComplex
// 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;
// mesh vertices in this domain
@@ -698,11 +697,8 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
- // storage for cell pointers to delete them
- std::list<Cell*> _trash;
-
+ // temporary store for omitted cells (generators of the highest dimension)
std::vector< std::set<Cell*, Less_Cell> > _store;
- std::set<Cell*, Less_Cell> _ecells;
// original cells of this cell complex
std::set<Cell*, Less_Cell> _ocells[4];
@@ -712,7 +708,7 @@ class CellComplex
int _dim;
- // Is the cell complex simplicial
+ // is the cell complex simplicial
bool _simplicial;
// enqueue cells in queue if they are not there already
@@ -739,7 +735,7 @@ class CellComplex
CellComplex(){}
~CellComplex();
- void emptyTrash();
+ //void emptyTrash();
// restore this cell complex to its original state
void restoreComplex();
@@ -771,8 +767,6 @@ class CellComplex
// 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
@@ -796,6 +790,7 @@ class CellComplex
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; }
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 750990851363d9b6899012372c64e01e41b6dd9c..73c36f04bacfe481ab406900ef330798a0e8624d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -55,11 +55,11 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
}
- if( cols == 0 ){
+ if( cols == 0 ){ // no dim-cells, no map
//_HMatrix[dim] = create_gmp_matrix_zero(rows, 1);
_HMatrix[dim] = NULL;
}
- else if( rows == 0){
+ else if( rows == 0){ // no dim-1-cells, maps everything to zero
_HMatrix[dim] = create_gmp_matrix_zero(1, cols);
//_HMatrix[dim] = NULL;
}
@@ -127,6 +127,7 @@ void ChainComplex::KerCod(int dim){
mpz_t elem;
mpz_init(elem);
+ // find the rank
while(rank < minRowCol){
gmp_matrix_get_elem(elem, rank+1, rank+1, normalForm->canonical);
if(mpz_cmp_si(elem,0) == 0) break;
@@ -684,10 +685,6 @@ 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;
}
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 7b0a96c0ab073252219ad3c94fa1dfb6327f5ae6..0a8d7bc04cd617d86c78c779add13b08f33219c7 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -317,63 +317,63 @@ void Homology::computeBettiNumbers(){
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 += "}";
+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 += ", {";
- }
- 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;
+ 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