diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index d7823ca2fbf12f0d142f47655ef9648f105c467e..f5f2c7b446aa365da1b8fd55de11cd0ba0d1bdcc 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -272,6 +272,8 @@ void CellComplex::removeCell(Cell* cell){
bdCell->removeCoboundaryCell(cell);
}
+ //cell->clearBoundary();
+ //cell->clearCoboundary();
}
@@ -283,7 +285,8 @@ void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, st
for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
Cell* cell = *cit;
citer it = Qset.find(cell);
- if(it == Qset.end()){
+ //citer it2 = _cells[cell->getDim()].find(cell);
+ if(it == Qset.end()){// && it2 != _cells[cell->getDim()].end()){
Qset.insert(cell);
Q.push(cell);
}
@@ -519,6 +522,23 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
cbdCell->removeBoundaryCell(c1);
cbdCell->addBoundaryCell(ori, newCell, true);
}
+ /*
+ for(std::list< std::pair<int, Cell*> >::iterator it = coboundary1.begin(); it != coboundary1.end(); it++){
+ Cell* cbdCell = (*it).second;
+ int ori = (*it).first;
+ cbdCell->removeBoundaryCell(c1);
+ if(!co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = coboundary2.begin(); it2 != coboundary2.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cbdCell) old = true;
+ }
+ if(!old) cbdCell->addBoundaryCell(ori, newCell, true);
+ }
+ else cbdCell->addBoundaryCell(ori, newCell, true);
+ }
+ */
+
for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
Cell* cbdCell = (*it).second;
int ori = (*it).first;
@@ -535,12 +555,31 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
else cbdCell->addBoundaryCell(ori, newCell, true);
}
+
for(std::list< std::pair<int, Cell* > >::iterator it = boundary1.begin(); it != boundary1.end(); it++){
Cell* bdCell = (*it).second;
int ori = (*it).first;
bdCell->removeCoboundaryCell(c1);
bdCell->addCoboundaryCell(ori, newCell, true);
}
+ /*
+ for(std::list< std::pair<int, Cell* > >::iterator it = boundary1.begin(); it != boundary1.end(); it++){
+ Cell* bdCell = (*it).second;
+ int ori = (*it).first;
+ bdCell->removeCoboundaryCell(c1);
+ if(co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = boundary2.begin(); it2 != boundary2.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *bdCell) old = true;
+ }
+ if(!old) bdCell->addCoboundaryCell(ori, newCell, true);
+ }
+ else bdCell->addCoboundaryCell(ori, newCell, true);
+
+ }
+ */
+
for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
Cell* bdCell = (*it).second;
int ori = (*it).first;
@@ -560,6 +599,8 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
_cells[dim].erase(c1);
_cells[dim].erase(c2);
+ //removeCell(c1);
+ //removeCell(c2);
//_trash.insert(c1);
//_trash.insert(c2);
_cells[dim].insert(newCell);
@@ -601,19 +642,21 @@ int CellComplex::cocombine(int dim){
int or2 = bd_c.back().first;
Cell* c1 = bd_c.front().second;
Cell* c2 = bd_c.back().second;
+ if( c1->getNumVertices() < getSize(dim) && c2->getNumVertices() < getSize(dim) ){
+
+ removeCell(s);
cbd_c = c1->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
cbd_c = c2->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
- removeCell(s);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
replaceCells(c1, c2, newCell, (or1 != or2), true);
cit = firstCell(dim);
count++;
-
+ }
}
removeCellQset(s, Qset);
@@ -653,23 +696,28 @@ int CellComplex::combine(int dim){
cbd_c = s->getOrientedCoboundary();
if(cbd_c.size() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
- && inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second) ){
+ && inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second)){
int or1 = cbd_c.front().first;
int or2 = cbd_c.back().first;
Cell* c1 = cbd_c.front().second;
Cell* c2 = cbd_c.back().second;
-
+ if( c1->getNumVertices() < getSize(dim) && c2->getNumVertices() < getSize(dim)){
+
+ removeCell(s);
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
bd_c = c2->getBoundary();
enqueueCells(bd_c, Q, Qset);
- removeCell(s);
+
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
replaceCells(c1, c2, newCell, (or1 != or2));
+ //removeCell(s);
cit = firstCell(dim);
//cit++;
count++;
+
+ }
}
removeCellQset(s, Qset);
@@ -807,10 +855,10 @@ int CellComplex::writeComplexMSH(const std::string &name){
void CellComplex::printComplex(int dim){
for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
- for(int i = 0; i < cell->getNumVertices(); i++){
- printf("%d ", cell->getVertex(i)->getNum());
- }
- printf("\n");
+ cell->printCell();
+ cell->printBoundary();
+ cell->printCoboundary();
+ printf("--- \n");
}
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 4a3e8d587785e219fdb77df1dadd35cf6ee2b3de..16f264c6a03fe3e717928d53352869551929c718 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -141,7 +141,15 @@ class Cell
}
return count;
}
-
+
+ virtual bool hasBoundary(Cell* cell){
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) return true;
+ }
+ return false;
+ }
+
virtual void clearBoundary() { _boundary.clear(); }
virtual void clearCoboundary() { _coboundary.clear(); }
@@ -160,6 +168,14 @@ class Cell
cell2->printCell();
}
}
+ virtual void printCoboundary() {
+ for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
+ printf("Coboundary cell orientation: %d ", (*it).first);
+ Cell* cell2 = (*it).second;
+ cell2->printCell();
+ }
+ }
+
virtual bool inSubdomain() const { return _inSubdomain; }
@@ -435,6 +451,7 @@ class CombinedCell : public Cell{
std::vector<MVertex*> _v;
std::vector<int> _vs;
std::list< std::pair<int, Cell*> > _cells;
+ int _num;
public:
@@ -442,10 +459,15 @@ class CombinedCell : public Cell{
_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));
@@ -461,14 +483,31 @@ class CombinedCell : public Cell{
_cells.push_back(*it);
}
- _boundary = c1->getOrientedBoundary();
+ if(!co) _boundary = c1->getOrientedBoundary();
+
+ std::list< std::pair<int, Cell*> > c1Boundary = c1->getOrientedBoundary();
std::list< std::pair<int, Cell*> > c2Boundary = c2->getOrientedBoundary();
+ /*
+ for(std::list< std::pair<int, Cell*> >::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
+ Cell* cell = (*it).second;
+ if(co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = c2Boundary.begin(); it2 != c2Boundary.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cell) old = true;
+ }
+ if(!old) _boundary.push_back(*it);
+ }
+ else _boundary.push_back(*it);
+ }
+ */
+
for(std::list< std::pair<int, Cell*> >::iterator it = c2Boundary.begin(); it != c2Boundary.end(); it++){
if(!orMatch) (*it).first = -1*(*it).first;
Cell* cell = (*it).second;
if(co){
bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = _boundary.begin(); it2 != _boundary.end(); it2++){
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Boundary.begin(); it2 != c1Boundary.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *cell) old = true;
}
@@ -477,14 +516,32 @@ class CombinedCell : public Cell{
else _boundary.push_back(*it);
}
- _coboundary = c1->getOrientedCoboundary();
+ if(co) _coboundary = c1->getOrientedCoboundary();
+ std::list< std::pair<int, Cell*> > c1Coboundary = c1->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > c2Coboundary = c2->getOrientedCoboundary();
+
+ /*
+ for(std::list< std::pair<int, Cell* > >::iterator it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
+
+ Cell* cell = (*it).second;
+ if(!co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = c2Coboundary.begin(); it2 != c2Coboundary.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cell) old = true;
+ }
+ if(!old) _coboundary.push_back(*it);
+ }
+ else _coboundary.push_back(*it);
+ }
+ */
+
for(std::list< std::pair<int, Cell* > >::iterator it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
if(!orMatch) (*it).first = -1*(*it).first;
Cell* cell = (*it).second;
if(!co){
bool old = false;
- for(std::list< std::pair<int, Cell* > >::iterator it2 = _coboundary.begin(); it2 != _coboundary.end(); it2++){
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Coboundary.begin(); it2 != c1Coboundary.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *cell) old = true;
}
@@ -497,6 +554,7 @@ class CombinedCell : public Cell{
~CombinedCell(){}
+ int getNum() { return _num; }
int getNumVertices() const { return _v.size(); }
MVertex* getVertex(int vertex) const { return _v.at(vertex); }
int getSortedVertex(int vertex) const { return _vs.at(vertex); }
@@ -506,6 +564,9 @@ class CombinedCell : public Cell{
// 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(int i = 0; i < _v.size(); i++){
if(_v.at(i)->getNum() == vertex) return true;
}
@@ -522,7 +583,6 @@ class CombinedCell : public Cell{
std::list< std::pair<int, Cell*> > getCells() { return _cells; }
int getNumCells() {return _cells.size();}
-
};
@@ -626,11 +686,17 @@ class CellComplex
//}
for(int i = 0; i < 4; i++){
+
+ for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
+ Cell* cell = *cit;
+ if(cell->combined()) delete cell;
+ }
_cells[i].clear();
for(citer cit = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
Cell* cell = *cit;
delete cell;
}
+ _cells2[i].clear();
}
}
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index fefc4e10063f08eb952b04073f1276866282f7df..22e44eb8342d09814a528e56877fc4a0badc4bb5 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -87,7 +87,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
mpz_set_si(elem, old_elem + (*it).first);
if( (old_elem + (*it).first) > 1 || (old_elem + (*it).first) < -1 ){
printf("Warning: Invalid incidence index: %d! HMatrix: %d. \n", (old_elem + (*it).first), dim);
- mpz_set_si(elem, 0);
+ //mpz_set_si(elem, 0);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 18cdc81db9e0968e08602d949da4584eaf92929d..f6edfd5a40c4ba246dfc193461c1722fcc4f9ae1 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -72,9 +72,12 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
void Homology::findGenerators(std::string fileName){
+
+
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(true);
+
//_cellComplex->checkCoherence();
if(getCombine()){
_cellComplex->combine(3);
@@ -90,6 +93,7 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
+ //for(int i = 0; i < 4; i++) _cellComplex->printComplex(i);
_cellComplex->writeComplexMSH(fileName);