From 6091a8d35e46ff1c78b0559da4f7526147cac418 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Tue, 19 May 2009 11:24:33 +0000
Subject: [PATCH] Code cleanup.
Added cell combining methods (about O(n^2)) to further reduce the size of
the cell complex before applying Smith normal form computation.
(Don't know why there seems to be some benchmarks/ files modified too,
Somebody else committed between my update and commit?)
---
Geo/CellComplex.cpp | 834 +++++++++++++++-----------------------
Geo/CellComplex.h | 423 ++++++++++++++-----
Geo/ChainComplex.cpp | 72 +++-
Geo/ChainComplex.h | 10 +
utils/misc/celldriver.cpp | 9 +
5 files changed, 727 insertions(+), 621 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 1d6337ad85..02a7c12836 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -70,9 +70,9 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
// insert cells into cell complex
// subdomain need to be inserted first!
- insert_cells2(true, true);
- insert_cells2(false, true);
- insert_cells2(false, false);
+ insert_cells(true, true);
+ insert_cells(false, true);
+ insert_cells(false, false);
// find mesh vertices in the domain
for(unsigned int j=0; j < domain.size(); j++) {
@@ -93,94 +93,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
- // attach induvivial tags to cells
+ // attach individual tags to cells
int tag = 1;
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
cell->setTag(tag);
tag++;
-
- // make sure that boundary and coboundary information is coherent
- std::list<Cell*> boundary = cell->getBoundary();
- boundary.sort(Less_Cell());
- boundary.unique(Equal_Cell());
- cell->setBoundary(boundary);
-
- std::list<Cell*> coboundary = cell->getCoboundary();
- coboundary.sort(Less_Cell());
- coboundary.unique(Equal_Cell());
- cell->setCoboundary(coboundary);
-
-
}
-
-
_originalCells[i] = _cells[i];
}
-
- _simplicial = true;
-
-}
-void CellComplex::insert_cells(bool subdomain, bool boundary){
-
- std::vector<GEntity*> domain;
-
- if(subdomain) domain = _subdomain;
- else if(boundary) domain = _boundary;
- else domain = _domain;
-
- std::vector<int> vertices;
- int vertex = 0;
-
- std::pair<citer, bool> insertInfo;
-
- for(unsigned int j=0; j < domain.size(); j++) {
- for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
- vertices.clear();
-
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
- MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
- vertices.push_back(vertex->getNum());
- //_cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain, vertex->x(), vertex->y(), vertex->z() )); // Add vertices
- Cell* cell = new ZeroSimplex(vertex->getNum(), subdomain, boundary);
- insertInfo = _cells[0].insert(cell);
- if(!insertInfo.second) delete cell;
- }
- if(domain.at(j)->getMeshElement(i)->getDim() == 3){
- Cell* cell = new ThreeSimplex(vertices, subdomain, boundary); // Add volumes
- insertInfo = _cells[3].insert(cell);
- if(!insertInfo.second) delete cell;
- }
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumFaces(); k++){
- vertices.clear();
- for(int l=0; l < domain.at(j)->getMeshElement(i)->getFace(k).getNumVertices(); l++){
- vertex = domain.at(j)->getMeshElement(i)->getFace(k).getVertex(l)->getNum();
- vertices.push_back(vertex);
- }
- Cell* cell = new TwoSimplex(vertices, subdomain, boundary); // Add faces
- insertInfo = _cells[2].insert(cell);
- if(!insertInfo.second) delete cell;
- }
-
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumEdges(); k++){
- vertices.clear();
- for(int l=0; l < domain.at(j)->getMeshElement(i)->getEdge(k).getNumVertices(); l++){
- vertex = domain.at(j)->getMeshElement(i)->getEdge(k).getVertex(l)->getNum();
- vertices.push_back(vertex);
- }
- Cell* cell = new OneSimplex(vertices, subdomain, boundary); // Add edges
- insertInfo = _cells[1].insert(cell);
- if(!insertInfo.second) delete cell;
- }
-
- }
- }
-
}
-void CellComplex::insert_cells2(bool subdomain, bool boundary){
+void CellComplex::insert_cells(bool subdomain, bool boundary){
// works only for simplcial meshes at the moment!
@@ -208,10 +134,10 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){
int dim = domain.at(j)->getMeshElement(i)->getDim();
int type = domain.at(j)->getMeshElement(i)->getTypeForMSH();
Cell* cell;
- if(dim == 3) cell = new ThreeSimplex(vertices, subdomain, boundary);
- else if(dim == 2) cell = new TwoSimplex(vertices, subdomain, boundary);
- else if(dim == 1) cell = new OneSimplex(vertices, subdomain, boundary);
- else cell = new ZeroSimplex(vertices.at(0), subdomain, boundary);
+ if(dim == 3) cell = new ThreeSimplex(vertices, 0, subdomain, boundary);
+ else if(dim == 2) cell = new TwoSimplex(vertices, 0, subdomain, boundary);
+ else if(dim == 1) cell = new OneSimplex(vertices, 0, subdomain, boundary);
+ else cell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
}
@@ -225,19 +151,23 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
Cell* newCell;
- if(dim == 3) newCell = new TwoSimplex(vertices, subdomain, boundary);
- else if(dim == 2) newCell = new OneSimplex(vertices, subdomain, boundary);
- else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), subdomain, boundary);
+ if(dim == 3) newCell = new TwoSimplex(vertices, 0, subdomain, boundary);
+ else if(dim == 2) newCell = new OneSimplex(vertices, 0, subdomain, boundary);
+ else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
delete newCell;
Cell* oldCell = *(insertInfo.first);
- oldCell->addCoboundaryCell(cell);
- cell->addBoundaryCell(oldCell);
+ if(!subdomain && !boundary){
+ int ori = cell->kappa(oldCell);
+ oldCell->addCoboundaryCell( ori, cell );
+ cell->addBoundaryCell( ori, oldCell);
+ }
}
- else{
- cell->addBoundaryCell(newCell);
- newCell->addCoboundaryCell(cell);
+ else if(!subdomain && !boundary) {
+ int ori = cell->kappa(newCell);
+ cell->addBoundaryCell( ori, newCell );
+ newCell->addCoboundaryCell( ori, cell);
}
}
}
@@ -249,6 +179,7 @@ void CellComplex::insertCell(Cell* cell){
_cells[cell->getDim()].insert(cell);
}
+
int Simplex::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
@@ -258,12 +189,88 @@ int Simplex::kappa(Cell* tau) const{
int value=1;
for(int i=0; i < tau->getNumVertices(); i++){
- if(this->getVertex(i) != tau->getVertex(i)) return value;
+ if(this->getSortedVertex(i) != tau->getSortedVertex(i)) return value;
value = value*-1;
}
return value;
}
+
+/*
+int Simplex::kappa(Cell* tau) const{
+ for(int i=0; i < tau->getNumVertices(); i++){
+ if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+ }
+
+ if(this->getDim() - tau->getDim() != 1) return 0;
+
+ int value=1;
+
+ for(int i = 0; i < this->getNumFacets(); i++){
+ std::vector<int> vTau = tau->getVertexVector();
+ std::vector<int> v;
+ this->getFacetVertices(i, v);
+ value = -1;
+
+ 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))) ) return 0;
+ }
+
+ if(tau->getDim() != 0) return 0;
+
+ if(tau->getVertex(0) == this->getVertex(0)) return -1;
+ else return 1;
+
+}
+*/
+
+int Quadrangle::kappa(Cell* tau) const{
+ for(int i=0; i < tau->getNumVertices(); i++){
+ if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+ }
+
+ if(tau->getDim() != 1) return 0;
+
+ int value=1;
+
+ std::vector<int> vTau = tau->getVertexVector();
+
+ for(int i = 0; i < this->getNumFacets(); i++){
+ std::vector<int> v;
+ this->getFacetVertices(i, v);
+ value = -1;
+ do {
+ value = value*-1;
+ if(v == vTau) return value;
+ }
+ while (std::next_permutation(v.begin(), v.end()) );
+ }
+
+ return value;
+}
+
+
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
Cell* cell;
@@ -294,122 +301,8 @@ std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex,
}
-std::vector<Cell*> CellComplex::bd(Cell* sigma){
- std::vector<Cell*> boundary;
- int dim = sigma->getDim();
- if(dim < 1) return boundary;
-
- if(simplicial()){ // faster
- std::vector<int> vertices = sigma->getVertexVector();
-
- for(int j=0; j < vertices.size(); j++){
- std::vector<int> bVertices;
-
- // simplicial mesh assumed here!
- for(int k=0; k < vertices.size(); k++){
- if (k !=j ) bVertices.push_back(vertices.at(k));
- }
- citer cit2 = findCell(dim-1, bVertices);
- if(cit2 != lastCell(dim-1)){
- boundary.push_back(*cit2);
- if(boundary.size() == sigma->getBdMaxSize()){
- return boundary;
- }
- }
- }
-
- }
- else{
- citer start = findCell(dim-1, sigma->getVertex(0), 0);
- if(start == lastCell(dim-1)) return boundary;
-
- citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1));
- if(end != lastCell(dim-1)) end++;
-
- //for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
- for(citer cit = start; cit != end; cit++){
- //if(kappa(sigma, *cit) != 0){
- if(sigma->kappa(*cit) != 0){
- boundary.push_back(*cit);
- if(boundary.size() == sigma->getBdMaxSize()){
- return boundary;
- }
- }
- }
- }
- return boundary;
-}
-std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::bdIt(Cell* sigma){
-
- int dim = sigma->getDim();
- std::vector< std::set<Cell*, Less_Cell>::iterator > boundary;
- if(dim < 1){
- return boundary;
- }
-
- citer start = findCell(dim-1, sigma->getVertex(0), 0);
- if(start == lastCell(dim-1)) return boundary;
-
- citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1));
- if(end != lastCell(dim-1)) end++;
-
- for(citer cit = start; cit != end; cit++){
- if(kappa(sigma, *cit) != 0){
- boundary.push_back(cit);
- if(boundary.size() == sigma->getBdMaxSize()){
- return boundary;
- }
- }
- }
-
- return boundary;
-}
-
-
-
-std::vector<Cell*> CellComplex::cbd(Cell* tau){
- std::vector<Cell*> coboundary;
- int dim = tau->getDim();
- if(dim > 2){
- return coboundary;
- }
-
- for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
- //if(kappa(*cit, tau) != 0){
- Cell* cell = *cit;
- if(cell->kappa(tau) != 0){
- coboundary.push_back(*cit);
- if(coboundary.size() == tau->getCbdMaxSize()){
- return coboundary;
- }
- }
- }
-
- return coboundary;
-}
-
-std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau){
-
- std::vector< std::set<Cell*, Less_Cell>::iterator > coboundary;
- int dim = tau->getDim();
- if(dim > 2){
- return coboundary;
- }
-
- for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
- if(kappa(*cit, tau) != 0){
- coboundary.push_back(cit);
- if(coboundary.size() == tau->getCbdMaxSize()){
- return coboundary;
- }
- }
- }
-
- return coboundary;
-}
-
-
void CellComplex::removeCell(Cell* cell){
+
_cells[cell->getDim()].erase(cell);
std::list<Cell*> coboundary = cell->getCoboundary();
@@ -419,34 +312,18 @@ void CellComplex::removeCell(Cell* cell){
Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell);
}
+
for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell);
}
-}
-void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
- Cell* c = *cell;
- int dim = c->getDim();
- _cells[dim].erase(cell);
-
}
void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
Qset.erase(cell);
}
-void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
- std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
- for(unsigned int i=0; i < cells.size(); i++){
- Cell* cell = *cells.at(i);
- citer cit = Qset.find(cell);
- if(cit == Qset.end()){
- Qset.insert(cell);
- Q.push(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;
@@ -458,18 +335,6 @@ void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, st
}
}
-void CellComplex::enqueueCells(std::list<Cell*>& cells, std::list<Cell*>& Q){
- for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
- Cell* cell = *cit;
- std::list<Cell*>::iterator it = std::find(Q.begin(), Q.end(), cell);
- if(it == Q.end()){
- Q.push_back(cell);
- }
- }
-}
-
-
-
int CellComplex::coreductionMrozek(Cell* generator){
int coreductions = 0;
@@ -482,6 +347,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
//removeCell(generator);
std::list<Cell*> bd_s;
+ //std::list< std::pair<int, Cell*> >bd_s;
std::list<Cell*> cbd_c;
Cell* s;
int round = 0;
@@ -512,145 +378,11 @@ int CellComplex::coreductionMrozek(Cell* generator){
return coreductions;
}
-int CellComplex::coreductionMrozek2(Cell* generator){
-
- int coreductions = 0;
-
- std::list<Cell*> Q;
-
- Q.push_back(generator);
-
- std::list<Cell*> bd_s;
- std::list<Cell*> cbd_c;
- Cell* s;
- int round = 0;
- while( !Q.empty() ){
- round++;
- //printf("%d ", round);
- s = Q.front();
- Q.pop_front();
- bd_s = s->getBoundary();
- if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
- removeCell(s);
- cbd_c = bd_s.front()->getCoboundary();
- enqueueCells(cbd_c, Q);
- removeCell(bd_s.front());
- coreductions++;
- }
- else if(bd_s.empty()){
- cbd_c = s->getCoboundary();
- enqueueCells(cbd_c, Q);
- }
-
- //Q.unique(Equal_Cell());
-
- }
- //printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
- return coreductions;
-}
-
-int CellComplex::coreductionMrozek3(Cell* generator){
-
- int coreductions = 0;
-
- std::queue<Cell*> Q;
- std::set<Cell*, Less_Cell> Qset;
-
- Q.push(generator);
- Qset.insert(generator);
-
- std::vector< std::set<Cell*, Less_Cell>::iterator > bd_s;
- std::vector< std::set<Cell*, Less_Cell>::iterator > cbd_c;
-
- Cell* s;
- int round = 0;
- while( !Q.empty() ){
- round++;
- //printf("%d ", round);
- s = Q.front();
- Q.pop();
- removeCellQset(s, Qset);
-
- bd_s = bdIt(s);
- if( bd_s.size() == 1 && inSameDomain(s, *bd_s.at(0)) ){
- removeCell(s);
- cbd_c = cbdIt(*bd_s.at(0));
- enqueueCellsIt(cbd_c, Q, Qset);
- removeCellIt(bd_s.at(0));
- coreductions++;
- }
- else if(bd_s.empty()){
- cbd_c = cbdIt(s);
- enqueueCellsIt(cbd_c, Q, Qset);
- }
- }
- //printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
- return coreductions;
-}
int CellComplex::coreduction(int dim){
if(dim < 0 || dim > 2) return 0;
- std::vector<Cell*> bd_c;
- int count = 0;
-
- bool coreduced = true;
- while (coreduced){
- coreduced = false;
- for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
- Cell* cell = *cit;
-
- bd_c = bd(cell);
- if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
- removeCell(cell);
- removeCell(bd_c.at(0));
- count++;
- coreduced = true;
- }
-
- }
- }
-
- return count;
-
-}
-
-int CellComplex::coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells){
-
- if(dim < 0 || dim > 2) return 0;
-
- std::vector<Cell*> bd_c;
- int count = 0;
-
- bool coreduced = true;
- while (coreduced){
- coreduced = false;
- for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
- Cell* cell = *cit;
-
- bd_c = bd(cell);
- if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
- removedCells.insert(cell);
- removedCells.insert(bd_c.at(0));
- removeCell(cell);
- removeCell(bd_c.at(0));
- count++;
- coreduced = true;
- }
-
- }
- }
-
- return count;
-
-}
-
-int CellComplex::coreduction2(int dim){
-
- if(dim < 0 || dim > 2) return 0;
-
std::list<Cell*> bd_c;
-
int count = 0;
bool coreduced = true;
@@ -669,40 +401,14 @@ int CellComplex::coreduction2(int dim){
}
}
+
return count;
+
}
-
-
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
- std::vector<Cell*> cbd_c;
- std::vector<Cell*> bd_c;
- int count = 0;
-
- bool reduced = true;
- while (reduced){
- reduced = false;
- for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
- Cell* cell = *cit;
- cbd_c = cbd(cell);
- if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
- removeCell(cell);
- removeCell(cbd_c.at(0));
- count++;
- reduced = true;
- }
- }
- }
- return count;
-}
-
-
-int CellComplex::reduction2(int dim){
- if(dim < 1 || dim > 3) return 0;
-
- std::list<Cell*> bd_c;
std::list<Cell*> cbd_c;
int count = 0;
@@ -728,48 +434,13 @@ int CellComplex::reduction2(int dim){
void CellComplex::reduceComplex(){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
- reduction2(3);
- reduction2(2);
- reduction2(1);
+ reduction(3);
+ reduction(2);
+ reduction(1);
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
}
-void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>& removedCells){
-
- std::set<Cell*, Less_Cell> generatorCells[4];
-
- printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- int i = generatorDim;
- while (getSize(i) != 0){
- citer cit = firstCell(i);
- Cell* cell = *cit;
- while(!cell->inSubdomain() && cit != lastCell(i)){
- cell = *cit;
- cit++;
- }
- generatorCells[i].insert(cell);
- removedCells.insert(cell);
- removeCell(cell);
- coreduction(0, removedCells);
- coreduction(1, removedCells);
- coreduction(2, removedCells);
- }
-
-
- for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
- Cell* cell = *cit;
- if(!cell->inSubdomain()) _cells[i].insert(cell);
- if(!cell->inSubdomain()) removedCells.insert(cell);
- }
- printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- return;
-
-}
void CellComplex::coreduceComplex(int generatorDim){
std::set<Cell*, Less_Cell> generatorCells[4];
@@ -788,26 +459,210 @@ void CellComplex::coreduceComplex(int generatorDim){
generatorCells[i].insert(cell);
removeCell(cell);
- //coreduction2(0);
- //coreduction2(1);
- //coreduction2(2);
+ //coreduction(0);
+ //coreduction(1);
+ //coreduction(2);
coreductionMrozek(cell);
}
if(generatorDim == i) break;
}
+ /*
for(int i = 0; i < 4; i++){
for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
Cell* cell = *cit;
- //if(!cell->inSubdomain()) _cells[i].insert(cell);
+ if(!cell->inSubdomain()) _cells[i].insert(cell);
}
}
+ */
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
return;
}
+
+void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co){
+
+ int dim = c1->getDim();
+
+ std::list< std::pair<int, Cell*> > coboundary1 = c1->getOrientedCoboundary();
+ std::list< std::pair<int, Cell*> > coboundary2 = c2->getOrientedCoboundary();
+ std::list< std::pair<int, Cell*> > boundary1 = c1->getOrientedBoundary();
+ std::list< std::pair<int, Cell*> > boundary2 = c2->getOrientedBoundary();
+
+
+ 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);
+ cbdCell->addBoundaryCell(ori, newCell);
+ }
+ for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
+ Cell* cbdCell = (*it).second;
+ int ori = (*it).first;
+ if(!orMatch) ori = -1*ori;
+ cbdCell->removeBoundaryCell(c2);
+ if(!co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = coboundary1.begin(); it2 != coboundary1.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cbdCell) old = true;
+ }
+ if(!old) cbdCell->addBoundaryCell(ori, newCell);
+ }
+ else cbdCell->addBoundaryCell(ori, newCell);
+ }
+
+ 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);
+ }
+ for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
+ Cell* bdCell = (*it).second;
+ int ori = (*it).first;
+ if(!orMatch) ori = -1*ori;
+ bdCell->removeCoboundaryCell(c2);
+ if(co){
+ bool old = false;
+ for(std::list< std::pair<int, Cell* > >::iterator it2 = boundary1.begin(); it2 != boundary1.end(); it2++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *bdCell) old = true;
+ }
+ if(!old) bdCell->addCoboundaryCell(ori, newCell);
+ }
+ else bdCell->addCoboundaryCell(ori, newCell);
+
+ }
+
+ _cells[dim].erase(c1);
+ _cells[dim].erase(c2);
+ _cells[dim].insert(newCell);
+
+
+
+ return;
+}
+
+int CellComplex::cocombine(int dim){
+
+ printf("Cell complex before cocombining: %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;
+
+ std::queue<Cell*> Q;
+ std::set<Cell*, Less_Cell> Qset;
+ std::list<Cell*> cbd_c;
+ std::list< std::pair< int, Cell*> > bd_c;
+ int count = 0;
+
+ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+
+ Cell* cell = *cit;
+ cbd_c = cell->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+ while(Q.size() != 0){
+
+ Cell* s = Q.front();
+ Q.pop();
+
+ bd_c = s->getOrientedBoundary();
+
+ if(bd_c.size() == 2 && !(*(bd_c.front().second) == *(bd_c.back().second))
+ && inSameDomain(s, bd_c.front().second) && inSameDomain(s, bd_c.back().second) ){
+
+ Cell* c1 = bd_c.front().second;
+ Cell* c2 = bd_c.back().second;
+
+ cbd_c = c1->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+ cbd_c = c2->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+
+ int or1 = bd_c.front().first;
+ int or2 = bd_c.back().first;
+
+ 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);
+
+ }
+ }
+
+ printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+
+ return count;
+}
+
+int CellComplex::combine(int dim){
+
+ printf("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;
+
+ std::queue<Cell*> Q;
+ std::set<Cell*, Less_Cell> Qset;
+ std::list< std::pair<int, Cell*> > cbd_c;
+ std::list<Cell*> bd_c;
+ int count = 0;
+
+
+ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+
+ Cell* cell = *cit;
+ bd_c = cell->getBoundary();
+ enqueueCells(bd_c, Q, Qset);
+ while(Q.size() != 0){
+
+ Cell* s = Q.front();
+ Q.pop();
+
+ 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) ){
+
+ Cell* c1 = cbd_c.front().second;
+ Cell* c2 = cbd_c.back().second;
+
+ bd_c = c1->getBoundary();
+ enqueueCells(bd_c, Q, Qset);
+ bd_c = c2->getBoundary();
+ enqueueCells(bd_c, Q, Qset);
+
+ int or1 = cbd_c.front().first;
+ int or2 = cbd_c.back().first;
+
+ removeCell(s);
+ CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
+ replaceCells(c1, c2, newCell, (or1 != or2));
+
+ cit = firstCell(dim);
+ count++;
+ }
+ removeCellQset(s, Qset);
+
+ }
+ }
+
+ printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices.\n",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+
+
+ return count;
+}
+
+
void CellComplex::repairComplex(int i){
printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -851,61 +706,24 @@ void CellComplex::swapSubdomain(){
}
}
- return;
-}
-
-
-void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
-
- std::vector<GEntity*> domain;
- if(subdomain) domain = _subdomain;
- else domain = _domain;
-
- for(unsigned int j=0; j < domain.size(); j++) {
- for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
- MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
- domainVertices.insert(vertex);
- }
- }
- }
-
-/* for(unsigned int j=0; j < _domain.size(); j++) {
- for(unsigned int i=0; i < _domain.at(j)->getNumMeshVertices(); i++){
- MVertex* vertex = _domain.at(j)->getMeshVertex(i);
- domainVertices.insert(vertex);
- }
-
- std::list<GFace*> faces = _domain.at(j)->faces();
- for(std::list<GFace*>::iterator fit = faces.begin(); fit != faces.end(); fit++){
- GFace* face = *fit;
- for(unsigned int i=0; i < face->getNumMeshVertices(); i++){
- MVertex* vertex = face->getMeshVertex(i);
- domainVertices.insert(vertex);
- }
- }
- std::list<GEdge*> edges = _domain.at(j)->edges();
- for(std::list<GEdge*>::iterator eit = edges.begin(); eit != edges.end(); eit++){
- GEdge* edge = *eit;
- for(unsigned int i=0; i < edge->getNumMeshVertices(); i++){
- MVertex* vertex = edge->getMeshVertex(i);
- domainVertices.insert(vertex);
- }
- }
- std::list<GVertex*> vertices = _domain.at(j)->vertices();
- for(std::list<GVertex*>::iterator vit = vertices.begin(); vit != vertices.end(); vit++){
- GVertex* vertex = *vit;
- for(unsigned int i=0; i < vertex->getNumMeshVertices(); i++){
- MVertex* mvertex = vertex->getMeshVertex(i);
- domainVertices.insert(mvertex);
+ // make subdomain closed
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()){
+ std::list<Cell*> boundary = cell->getBoundary();
+ for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
+ Cell* bdCell = *it;
+ bdCell->setInSubdomain(true);
+ }
}
}
-
}
- */
+
return;
}
+
int CellComplex::writeComplexMSH(const std::string &name){
@@ -923,10 +741,6 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$Nodes\n");
- //std::set<MVertex*, Less_MVertex> domainVertices;
- //getDomainVertices(domainVertices, true);
- //getDomainVertices(domainVertices, false);
-
fprintf(fp, "%d\n", _domainVertices.size());
for(std::set<MVertex*, Less_MVertex>::iterator vit = _domainVertices.begin(); vit != _domainVertices.end(); vit++){
@@ -938,7 +752,15 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$EndNodes\n");
fprintf(fp, "$Elements\n");
- fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size());
+ int count = 0;
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ count = count + cell->getNumCells();
+ }
+ }
+
+ fprintf(fp, "%d\n", count);
int physical = 0;
@@ -950,13 +772,17 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0));
}
-
+ std::list< std::pair<int, Cell*> > cells;
for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
Cell* edge = *cit;
if(edge->inSubdomain()) physical = 3;
else if(edge->onDomainBoundary()) physical = 2;
else physical = 1;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", edge->getTag(), 1, 3, 0, 0, physical, edge->getVertex(0), edge->getVertex(1));
+ cells = edge->getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+ Cell* cell = (*it).second;
+ fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getTag(), 1, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1));
+ }
}
for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
@@ -964,14 +790,22 @@ int CellComplex::writeComplexMSH(const std::string &name){
if(face->inSubdomain()) physical = 3;
else if(face->onDomainBoundary()) physical = 2;
else physical = 1;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", face->getTag(), 2, 3, 0, 0, physical, face->getVertex(0), face->getVertex(1), face->getVertex(2));
+ cells = face->getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+ Cell* cell = (*it).second;
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getTag(), 2, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2));
+ }
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
Cell* volume = *cit;
if(volume->inSubdomain()) physical = 3;
else if(volume->onDomainBoundary()) physical = 2;
else physical = 1;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", volume->getTag(), 4, 3, 0, 0, physical, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
+ cells = volume->getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+ Cell* cell = (*it).second;
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getTag(), 4, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2), cell->getVertex(3));
+ }
}
fprintf(fp, "$EndElements\n");
@@ -982,11 +816,11 @@ int CellComplex::writeComplexMSH(const std::string &name){
}
-void CellComplex::printComplex(int dim, bool subcomplex){
+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++){
- if(!subcomplex && !cell->inSubdomain()) printf("%d ", cell->getVertex(i));
+ printf("%d ", cell->getVertex(i));
}
printf("\n");
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index ebeafaedec..9ab40aadf5 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -21,7 +21,7 @@
#include "GFace.h"
#include "GVertex.h"
-// Abstract class representing a cell of a cell complex.
+// Abstract class representing an elemtary cell of a cell complex.
class Cell
{
protected:
@@ -40,12 +40,16 @@ class Cell
// whether this cell belongs to the boundary of a cell complex
bool _onDomainBoundary;
+ bool _combined;
+
// unique tag for each cell
int _tag;
+ int _index;
// cells on the boundary and on the coboundary of thhis cell
- std::list<Cell*> _boundary;
- std::list<Cell*> _coboundary;
+ std::list< std::pair<int, Cell*> > _boundary;
+ std::list< std::pair<int, Cell*> > _coboundary;
+
public:
Cell(){}
@@ -54,11 +58,14 @@ class Cell
virtual int getDim() const { return _dim; };
virtual int getTag() const { return _tag; };
virtual void setTag(int tag) { _tag = tag; };
+ virtual int getIndex() const { return _index; };
+ virtual void setIndex(int index) { _index = index; };
+
// get the number of vertices this cell has
virtual int getNumVertices() const = 0; //{return _vertices.size();}
-
virtual int getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
+ virtual int getSortedVertex(int vertex) const = 0;
virtual std::vector<int> getVertexVector() const = 0;
// returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
@@ -67,37 +74,58 @@ class Cell
virtual int kappa(Cell* tau) const = 0;
// true if this cell has given vertex
- virtual bool hasVertex(int vertex) const = 0;
-
+ virtual bool hasVertex(int vertex) const =0;
+
virtual unsigned int getBdMaxSize() const { return _bdMaxSize; };
virtual unsigned int getCbdMaxSize() const { return _cbdMaxSize; };
virtual int getBoundarySize() { return _boundary.size(); }
virtual int getCoboundarySize() { return _coboundary.size(); }
-
- virtual void setBoundary(std::list<Cell*> boundary){ _boundary = boundary; }
- virtual void setBoundary(std::vector<Cell*> boundary){
- for(int i = 0; i < boundary.size(); i++) _boundary.push_back(boundary.at(i)); }
- virtual void setCoboundary(std::list<Cell*> coboundary){ _coboundary = coboundary; }
- virtual void setCoboundary(std::vector<Cell*> coboundary){
- for(int i = 0; i < coboundary.size(); i++) _coboundary.push_back(coboundary.at(i)); }
- virtual std::list<Cell*> getBoundary() { return _boundary; }
- virtual std::list<Cell*> getCoboundary() { return _coboundary; }
- virtual void addBoundaryCell(Cell* cell) { _boundary.push_back(cell); }
- virtual void addCoboundaryCell(Cell* cell) { _coboundary.push_back(cell); }
+
+ virtual std::list< std::pair<int, Cell*> > getOrientedBoundary() { return _boundary; }
+ virtual std::list< Cell* > getBoundary() {
+ std::list<Cell*> boundary;
+ for( std::list< std::pair<int, Cell*> >::iterator it= _boundary.begin();it!= _boundary.end();it++){
+ Cell* cell = (*it).second;
+ boundary.push_back(cell);
+ }
+ return boundary;
+ }
+ virtual std::list<std::pair<int, Cell*> >getOrientedCoboundary() { return _coboundary; }
+ virtual std::list< Cell* > getCoboundary() {
+ std::list<Cell*> coboundary;
+ for( std::list< std::pair<int, Cell*> >::iterator it= _coboundary.begin();it!= _coboundary.end();it++){
+ Cell* cell = (*it).second;
+ coboundary.push_back(cell);
+ }
+ return coboundary;
+ }
+
+ virtual void addBoundaryCell(int orientation, Cell* cell) { _boundary.push_back( std::make_pair(orientation, cell)); }
+ virtual void addCoboundaryCell(int orientation, Cell* cell) { _coboundary.push_back( std::make_pair(orientation, cell)); }
+
virtual void removeBoundaryCell(Cell* cell) {
- for(std::list<Cell*>::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- Cell* cell2 = *it;
- if(cell2->getTag() == cell->getTag()) { _boundary.erase(it); break; }
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) { _boundary.erase(it); break; }
}
}
virtual void removeCoboundaryCell(Cell* cell) {
- for(std::list<Cell*>::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
- Cell* cell2 = *it;
- if(cell2->getTag() == cell->getTag()) { _coboundary.erase(it); break; }
+ for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
+ Cell* cell2 = (*it).second;
+ if(*cell2 == *cell) { _coboundary.erase(it); break; }
}
}
+ virtual void printBoundary() {
+ for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
+ printf("Boundary cell orientation: %d ", (*it).first);
+ Cell* cell2 = (*it).second;
+ cell2->printCell();
+ }
+ }
+
+
virtual bool inSubdomain() const { return _inSubdomain; }
virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
@@ -114,26 +142,30 @@ class Cell
virtual int getNumFacets() const { return 0; }
virtual void getFacetVertices(const int num, std::vector<int> &v) const {};
-
- bool operator==(const Cell* c2){
- if(this->getNumVertices() != c2->getNumVertices()){
+
+ virtual bool combined() { 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;
}
for(int i=0; i < this->getNumVertices();i++){
- if(this->getVertex(i) != c2->getVertex(i)){
+ if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
return false;
}
}
return true;
}
- bool operator<(const Cell* c2) const {
- if(this->getNumVertices() != c2->getNumVertices()){
- return (this->getNumVertices() < c2->getNumVertices());
+ bool operator<(const Cell& c2) const {
+ if(this->getNumVertices() != c2.getNumVertices()){
+ return (this->getNumVertices() < c2.getNumVertices());
}
for(int i=0; i < this->getNumVertices();i++){
- if(this->getVertex(i) < c2->getVertex(i)) return true;
- else if (this->getVertex(i) > c2->getVertex(i)) return false;
+ if(this->getSortedVertex(i) < c2.getSortedVertex(i)) return true;
+ else if (this->getSortedVertex(i) > c2.getSortedVertex(i)) return false;
}
return false;
}
@@ -141,15 +173,16 @@ class Cell
};
+
// Simplicial cell type.
class Simplex : public Cell
{
public:
Simplex() : Cell() {
-
+ _combined = false;
}
- virtual ~Simplex(){}
+ ~Simplex(){}
// kappa for simplices
// checks whether lower dimensional simplex tau is on the boundary of this cell
@@ -168,8 +201,9 @@ class ZeroSimplex : public Simplex
double _x, _y, _z;
public:
- ZeroSimplex(int vertex, bool subdomain=false, bool boundary=false, double x=0, double y=0, double z=0) : Simplex(){
+ ZeroSimplex(int vertex, int tag=0, bool subdomain=false, bool boundary=false, double x=0, double y=0, double z=0) : Simplex(){
_v = vertex;
+ _tag = tag;
_dim = 0;
_bdMaxSize = 0;
_cbdMaxSize = 1000; // big number
@@ -179,20 +213,22 @@ class ZeroSimplex : public Simplex
_inSubdomain = subdomain;
_onDomainBoundary = boundary;
}
- virtual ~ZeroSimplex(){}
+ ~ZeroSimplex(){}
- virtual int getDim() const { return 0; }
- virtual int getNumVertices() const { return 1; }
- virtual int getVertex(int vertex) const {return _v; }
- virtual bool hasVertex(int vertex) const {return (_v == vertex); }
+ int getDim() const { return 0; }
+ int getNumVertices() const { return 1; }
+ int getVertex(int vertex) const {return _v; }
+ int getSortedVertex(int vertex) const {return _v; }
+ bool hasVertex(int vertex) const {return (_v == vertex); }
- virtual std::vector<int> getVertexVector() const { return std::vector<int>(1,_v); }
-
- virtual inline double x() const { return _x; }
- virtual inline double y() const { return _y; }
- virtual inline double z() const { return _z; }
+ std::vector<int> getVertexVector() const { return std::vector<int>(1,_v); }
+ std::vector<int> getSortedVertexVector() const { return std::vector<int>(1,_v); }
- virtual void printCell()const { printf("Vertices: %d\n", _v); }
+ inline double x() const { return _x; }
+ inline double y() const { return _y; }
+ inline double z() const { return _z; }
+
+ void printCell() const { printf("Vertices: %d, in subdomain: %d \n", _v, _inSubdomain); }
};
@@ -203,12 +239,16 @@ class OneSimplex : public Simplex
// numbers of the vertices of this simplex
// same as the corresponding vertices in the finite element mesh
int _v[2];
+ int _vs[2];
public:
- OneSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
- sort(vertices.begin(), vertices.end());
+ OneSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
+ sort(vertices.begin(), vertices.end());
+ _vs[0] = vertices.at(0);
+ _vs[1] = vertices.at(1);
+ _tag = tag;
_dim = 1;
_bdMaxSize = 2;
_cbdMaxSize = 1000;
@@ -222,8 +262,8 @@ class OneSimplex : public Simplex
// used to find another definite one simplex from the cell complex
// first vertex gives the lower bound from where to look
OneSimplex(int vertex, int dummy){
- _v[0] = vertex;
- _v[1] = dummy;
+ _vs[0] = vertex;
+ _vs[1] = dummy;
}
~OneSimplex(){}
@@ -232,18 +272,23 @@ class OneSimplex : public Simplex
int getNumVertices() const { return 2; }
int getNumFacets() const { return 2; }
int getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) );
+ }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(1);
v[0] = _v[num];
}
+ //int kappa(Cell* tau) const;
- void printCell() const { printf("Vertices: %d %d\n", _v[0], _v[1]); }
+ void printCell() const { printf("Vertices: %d %d, in subdomain: %d \n", _v[0], _v[1], _inSubdomain); }
};
// Two simplex cell type.
@@ -253,6 +298,7 @@ class TwoSimplex : public Simplex
// numbers of the vertices of this simplex
// same as the corresponding vertices in the finite element mesh
int _v[3];
+ int _vs[3];
int edges_tri(const int edge, const int vert) const{
static const int e[3][2] = {
@@ -265,11 +311,15 @@ class TwoSimplex : public Simplex
public:
- TwoSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
- sort(vertices.begin(), vertices.end());
+ TwoSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
_v[2] = vertices.at(2);
+ sort(vertices.begin(), vertices.end());
+ _vs[0] = vertices.at(0);
+ _vs[1] = vertices.at(1);
+ _vs[2] = vertices.at(2);
+ _tag = tag;
_dim = 2;
_bdMaxSize = 3;
_cbdMaxSize = 2;
@@ -289,10 +339,14 @@ class TwoSimplex : public Simplex
int getNumVertices() const { return 3; }
int getNumFacets() const { return 3; }
int getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
bool hasVertex(int vertex) const {return
(_v[0] == vertex || _v[1] == vertex || _v[2] == vertex); }
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) );
+ }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(2);
@@ -300,7 +354,7 @@ class TwoSimplex : public Simplex
v[1] = _v[edges_tri(num, 1)];
}
- void printCell() const { printf("Vertices: %d %d %d\n", _v[0], _v[1], _v[2]); }
+ void printCell() const { printf("Vertices: %d %d %d, in subdomain: %d\n", _v[0], _v[1], _v[2], _inSubdomain); }
};
// Three simplex cell type.
@@ -310,6 +364,7 @@ class ThreeSimplex : public Simplex
// numbers of the vertices of this simplex
// same as the corresponding vertices in the finite element mesh
int _v[4];
+ int _vs[4];
int faces_tetra(const int face, const int vert) const{
static const int f[4][3] = {
@@ -323,12 +378,17 @@ class ThreeSimplex : public Simplex
public:
- ThreeSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
- sort(vertices.begin(), vertices.end());
+ ThreeSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
_v[2] = vertices.at(2);
_v[3] = vertices.at(3);
+ sort(vertices.begin(), vertices.end());
+ _vs[0] = vertices.at(0);
+ _vs[1] = vertices.at(1);
+ _vs[2] = vertices.at(2);
+ _vs[3] = vertices.at(3);
+ _tag = tag;
_dim = 3;
_bdMaxSize = 4;
_cbdMaxSize = 0;
@@ -349,10 +409,13 @@ class ThreeSimplex : public Simplex
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
int getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
bool hasVertex(int vertex) const {return
(_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(3);
@@ -361,34 +424,107 @@ class ThreeSimplex : public Simplex
v[2] = _v[faces_tetra(num, 2)];
}
- virtual void printCell() const { printf("Vertices: %d %d %d %d\n", _v[0], _v[1], _v[2], _v[3]); }
+ virtual void printCell() const { printf("Vertices: %d %d %d %d, in subdomain: %d \n", _v[0], _v[1], _v[2], _v[3], _inSubdomain); }
};
+
+class Quadrangle : public Cell
+{
+ private:
+
+ int _v[4];
+ int _vs[4];
+
+ int edges_quad(const int edge, const int vert) const{
+ static const int e[4][2] = {
+ {0, 1},
+ {1, 2},
+ {2, 3},
+ {3, 0}
+ };
+ return e[edge][vert];
+ }
+
+ public:
+
+ Quadrangle(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Cell(){
+ _v[0] = vertices.at(0);
+ _v[1] = vertices.at(1);
+ _v[2] = vertices.at(2);
+ _v[3] = vertices.at(3);
+ sort(vertices.begin(), vertices.end());
+ _vs[0] = vertices.at(0);
+ _vs[1] = vertices.at(1);
+ _vs[2] = vertices.at(2);
+ _vs[3] = vertices.at(3);
+ _tag = tag;
+ _dim = 2;
+ _bdMaxSize = 4;
+ _cbdMaxSize = 2;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ }
+ ~Quadrangle(){}
+
+ int getDim() const { return 2; }
+ int getNumVertices() const { return 4; }
+ int getNumFacets() const { return 4; }
+ int getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const { return _vs[vertex]; }
+ int kappa(Cell* tau) const;
+
+ bool hasVertex(int vertex) const {return
+ (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
+ std::vector<int> getVertexVector() const {
+ return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
+
+ void getFacetVertices(const int num, std::vector<int> &v) const {
+ v.resize(2);
+ v[0] = _v[edges_quad(num, 0)];
+ v[1] = _v[edges_quad(num, 1)];
+ }
+
+};
+
+
// Ordering for the cells.
class Less_Cell{
public:
bool operator()(const Cell* c1, const Cell* c2) const {
- // subcomplex cells in the end
- //if(c1->inSubdomain() != c2->inSubdomain()){
- // if(c1->inSubdomain()) return false;
- // else return true;
- //}
-
// 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->getVertex(i) < c2->getVertex(i)) return true;
- else if (c1->getVertex(i) > c2->getVertex(i)) return false;
+ if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
+ else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
+ }
+
+ /*
+ std::vector<int> c1v = c1->getVertexVector();
+ std::vector<int> c2v = c2->getVertexVector();
+ std::sort(c1v.begin(), c1v.end());
+ std::sort(c2v.begin(), c2v.end());
+
+ for(int i=0; i < c1v.size();i++){
+ if(c1v.at(i) < c2v.at(i)) return true;
+ else if (c1v.at(i) > c2v.at(i)) return false;
}
+ */
return false;
}
};
+
class Equal_Cell{
public:
bool operator ()(const Cell* c1, const Cell* c2){
@@ -396,7 +532,7 @@ class Equal_Cell{
return false;
}
for(int i=0; i < c1->getNumVertices();i++){
- if(c1->getVertex(i) != c2->getVertex(i)){
+ if(c1->getSortedVertex(i) != c2->getSortedVertex(i)){
return false;
}
}
@@ -413,6 +549,104 @@ class Less_MVertex{
}
};
+class CombinedCell : public Cell{
+
+ private:
+ std::vector<int> _v;
+ std::vector<int> _vs;
+ std::list< std::pair<int, Cell*> > _cells;
+
+ public:
+
+ CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
+ _tag = c1->getTag();
+ _dim = c1->getDim();
+ _bdMaxSize = 1000000;
+ _cbdMaxSize = 1000000;
+ _inSubdomain = c1->inSubdomain();
+ _onDomainBoundary = c1->onDomainBoundary();
+ _combined = true;
+
+ _v = c1->getVertexVector();
+ for(int i = 0; i < c2->getNumVertices(); i++){
+ if(!this->hasVertex(c2->getVertex(i))) _v.push_back(c2->getVertex(i));
+ }
+
+ _vs = _v;
+ std::sort(_vs.begin(), _vs.end());
+
+ _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 = c1->getOrientedBoundary();
+ std::list< std::pair<int, Cell*> > c2Boundary = c2->getOrientedBoundary();
+ 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++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cell) old = true;
+ }
+ if(!old) _boundary.push_back(*it);
+ }
+ else _boundary.push_back(*it);
+ }
+
+ _coboundary = c1->getOrientedCoboundary();
+ std::list< std::pair<int, Cell*> > c2Coboundary = c2->getOrientedCoboundary();
+ 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++){
+ Cell* cell2 = (*it2).second;
+ if(*cell2 == *cell) old = true;
+ }
+ if(!old) _coboundary.push_back(*it);
+ }
+ else _coboundary.push_back(*it);
+ }
+
+ }
+
+ ~CombinedCell(){}
+
+ int getNumVertices() const { return _v.size(); }
+ int getVertex(int vertex) const { return _v.at(vertex); }
+ int getSortedVertex(int vertex) const { return _vs.at(vertex); }
+ std::vector<int> getVertexVector() const { return _v; }
+
+ int kappa(Cell* tau) const { return 0; }
+
+ // true if this cell has given vertex
+ bool hasVertex(int vertex) const {
+ for(int i = 0; i < _v.size(); i++){
+ if(_v.at(i) == vertex) return true;
+ }
+ return false;
+ }
+
+ virtual void printCell() const {
+ printf("Vertices: ");
+ for(int i = 0; i < this->getNumVertices(); i++){
+ printf("%d ", this->getVertex(i));
+ }
+ printf(", in subdomain: %d\n", _inSubdomain);
+ }
+
+ std::list< std::pair<int, Cell*> > getCells() { return _cells; }
+ int getNumCells() {return _cells.size();}
+
+
+};
+
// A class representing a cell complex made out of a finite element mesh.
class CellComplex
{
@@ -434,8 +668,6 @@ class CellComplex
std::set<Cell*, Less_Cell> _originalCells[4];
- bool _simplicial;
-
public:
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
@@ -444,16 +676,12 @@ class CellComplex
// enqueue cells in queue if they are not there already
virtual void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
- virtual void enqueueCells(std::list<Cell*>& cells, std::list<Cell*>& Q);
- virtual void enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
- std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
-
// remove cell from the queue set
virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
// insert cells into this cell complex
+ //virtual void insert_cells(bool subdomain, bool boundary);
virtual void insert_cells(bool subdomain, bool boundary);
- virtual void insert_cells2(bool subdomain, bool boundary);
public:
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
@@ -467,16 +695,12 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
virtual ~CellComplex(){}
- virtual bool simplicial() { return _simplicial; }
// get the number of certain dimensional cells
virtual int getSize(int dim){ return _cells[dim].size(); }
virtual std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
-
- // get all the finite element mesh vertices in the domain of this cell complex
- virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain);
-
+
// iterators to the first and last cells of certain dimension
virtual citer firstCell(int dim) {return _cells[dim].begin(); }
virtual citer lastCell(int dim) {return _cells[dim].end(); }
@@ -489,16 +713,10 @@ class CellComplex
// implementation will vary depending on cell type
virtual inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
- // cells on the boundary of a cell
- virtual std::vector<Cell*> bd(Cell* sigma);
- virtual std::vector< std::set<Cell*, Less_Cell>::iterator > bdIt(Cell* sigma);
- // cells on the coboundary of a cell
- virtual std::vector<Cell*> cbd(Cell* tau);
- virtual std::vector< std::set<Cell*, Less_Cell>::iterator > cbdIt(Cell* tau);
-
// remove a cell from this cell complex
virtual void removeCell(Cell* cell);
- virtual void removeCellIt(std::set<Cell*, Less_Cell>::iterator cell);
+ virtual void replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co=false);
+
virtual void insertCell(Cell* cell);
@@ -509,54 +727,35 @@ class CellComplex
// coreduction of this cell complex
// removes corection pairs of cells of dimension dim and dim+1
virtual int coreduction(int dim);
- virtual int coreduction2(int dim);
+
// stores removed cells
- virtual int coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells);
+
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
virtual int reduction(int dim);
- virtual int reduction2(int dim);
// useful functions for (co)reduction of cell complex
virtual void reduceComplex();
// coreduction up to generators of dimension generatorDim
virtual void coreduceComplex(int generatorDim=3);
- // stores cells removed after removal of generatos of dimension generatorDim
- virtual void coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>& removedCells);
-
-
// queued coreduction presented in Mrozek's paper
// slower, but produces cleaner result
virtual int coreductionMrozek(Cell* generator);
- virtual int coreductionMrozek2(Cell* generator);
- virtual int coreductionMrozek3(Cell* generator);
-
- // never use, O(n^2) complexity
- virtual void restoreComplex(){
- for(int i = 0; i < 4; i++) _cells[i] = _originalCells[i];
-
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->setBoundary(bd(cell));
- cell->setCoboundary(cbd(cell));
- }
- }
-
- }
-
+
// add every volume, face and edge its missing boundary cells
virtual void repairComplex(int i=3);
// change non-subdomain cells to be in subdomain, subdomain cells to not to be in subdomain
virtual void swapSubdomain();
// print the vertices of cells of certain dimension
- virtual void printComplex(int dim, bool subcomplex);
+ virtual void printComplex(int dim);
// write this cell complex in 2.0 MSH ASCII format
virtual int writeComplexMSH(const std::string &name);
+ virtual int combine(int dim);
+ virtual int cocombine(int dim);
};
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index c0eae87ce5..88d02b75d2 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -24,15 +24,26 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
unsigned int rows = cellComplex->getSize(dim-1);
-
+ int index = 1;
// ignore subdomain cells
for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
- if(cell->inSubdomain()) cols--;
+ cell->setIndex(index);
+ index++;
+ if(cell->inSubdomain()){
+ index--;
+ cols--;
+ }
}
+ index = 1;
for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){
Cell* cell = *cit;
- if(cell->inSubdomain()) rows--;
+ cell->setIndex(index);
+ index++;
+ if(cell->inSubdomain()){
+ index--;
+ rows--;
+ }
}
@@ -45,7 +56,34 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
//_HMatrix[dim] = NULL;
}
+ else{
+
+ mpz_t elem;
+ mpz_init(elem);
+
+ _HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
+ for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
+ Cell* cell = *cit;
+ if(!cell->inSubdomain()){
+ std::list< std::pair<int,Cell*> >bdCell = cell->getOrientedBoundary();
+ for(std::list< std::pair<int, Cell*> >::iterator it = bdCell.begin(); it != bdCell.end(); it++){
+ Cell* bdCell = (*it).second;
+ if(!bdCell->inSubdomain()){
+ int old_elem = 0;
+ gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ old_elem = mpz_get_si(elem);
+ mpz_set_si(elem, old_elem + (*it).first);
+ gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
+ }
+ }
+ }
+ }
+
+ mpz_clear(elem);
+
+ }
+ /*
else{
long int elems[rows*cols];
@@ -59,6 +97,16 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
Cell* lowcell = *low;
Cell* highcell = *high;
if(!(highcell->inSubdomain() || lowcell->inSubdomain())){
+
+
+ std::list< std::pair<int, Cell*> >bdHigh = highcell->getBoundary();
+ for(std::list< std::pair<int, Cell*> >::iterator it = bdHigh.begin(); it != bdHigh.end(); it++){
+ Cell* bdCell = (*it).second;
+ if(bdCell->getTag() == lowcell->getTag()) elems[i] = (*it).first;
+ else elems[i] = 0;
+ }
+
+
elems[i] = cellComplex->kappa(*high, *low);
i++;
}
@@ -69,7 +117,9 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
_HMatrix[dim] = create_gmp_matrix_int(rows, cols, elems);
}
-
+ */
+
+
_kerH[dim] = NULL;
_codH[dim] = NULL;
_JMatrix[dim] = NULL;
@@ -220,7 +270,7 @@ void ChainComplex::Quotient(int dim){
mpz_t elem;
mpz_init(elem);
-
+
for(int i = 1; i <= cols; i++){
gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
if(mpz_cmp_si(elem,0) == 0){
@@ -414,13 +464,17 @@ int Chain::writeChainMSH(const std::string &name){
fprintf(fp, "4 \n");
fprintf(fp, "0 \n");
fprintf(fp, "1 \n");
- fprintf(fp, "%d \n", getSize());
+ fprintf(fp, "%d \n", getNumCells());
fprintf(fp, "0 \n");
+ std::list< std::pair<int, Cell*> > cells;
for(int i = 0; i < getSize(); i++){
-
- fprintf(fp, "%d %d \n", getCell(i)->getTag(), getCoeff(i));
-
+ Cell* cell = getCell(i);
+ cells = cell->getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+ Cell* cell2 = (*it).second;
+ fprintf(fp, "%d %d \n", cell2->getTag(), getCoeff(i)*(*it).first );
+ }
}
fprintf(fp, "$EndElementData\n");
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 0ae7756866..adc5450660 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -142,6 +142,16 @@ class Chain{
// number of cells in this chain
virtual int getSize() { return _cells.size(); }
+ virtual int getNumCells() {
+ int count = 0;
+ for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
+ Cell* cell = (*it).first;
+ count = count + cell->getNumCells();
+ }
+ return count;
+ }
+
+
// get/set chain name
virtual std::string getName() { return _name; }
virtual void setName(std::string name) { _name=name; }
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index 14040cae5b..5bc4920a9a 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -68,6 +68,11 @@ int main(int argc, char **argv)
// reduce the complex in order to ease homology computation
complex->reduceComplex();
+ // perform series of cell combinatios in order to further ease homology computation
+ complex->combine(3);
+ complex->combine(2);
+ complex->combine(1);
+
// create a chain complex of a cell complex (construct boundary operator matrices)
ChainComplex* chains = new ChainComplex(complex);
// compute the homology using the boundary operator matrices
@@ -101,6 +106,10 @@ int main(int argc, char **argv)
// reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
complex2->coreduceComplex(0);
+ // perform series of "dual complex" cell combinations in order to ease homology computation
+ complex2->cocombine(1);
+ complex2->cocombine(2);
+
// create a chain complex of a cell complex (construct boundary operator matrices)
ChainComplex* dualChains = new ChainComplex(complex2);
--
GitLab