From b6489b1535389dcc797343ccc7c0438da25c9850 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Wed, 27 May 2009 07:11:22 +0000
Subject: [PATCH] Small update on reduceComplex and coreduceComplex methods to
correspond current insight of reduction and coreduction: coreduction is the
same thing to the dual cell complex as reduction is to the primal cell
complex, and hence has no "priviledges" over that.
When coreduction is applied to primal complex, we obtain the thick cuts
of generators of homology groups obtained via reduction of the primal
complex.
---
Geo/CellComplex.cpp | 303 +++++++++++++++++++++++++++++++++++---
Geo/CellComplex.h | 113 +++++++++++---
Geo/ChainComplex.cpp | 80 +++++-----
Geo/ChainComplex.h | 47 +++---
utils/misc/celldriver.cpp | 9 +-
5 files changed, 454 insertions(+), 98 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 02a7c12836..892fc24c66 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -12,6 +12,8 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_domain = domain;
_subdomain = subdomain;
+ _dim = 0;
+
// determine mesh entities on boundary of the domain
bool duplicate = false;
for(unsigned int j=0; j < _domain.size(); j++){
@@ -101,7 +103,9 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
cell->setTag(tag);
tag++;
}
- _originalCells[i] = _cells[i];
+ //_originalCells[i] = _cells[i];
+ _betti[i] = 0;
+ if(getSize(i) > _dim) _dim = i;
}
}
@@ -270,7 +274,7 @@ int Quadrangle::kappa(Cell* tau) const{
return value;
}
-
+/*
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
Cell* cell;
@@ -299,7 +303,7 @@ std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex,
delete cell;
return cit;
}
-
+*/
void CellComplex::removeCell(Cell* cell){
@@ -335,7 +339,7 @@ void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, st
}
}
-int CellComplex::coreductionMrozek(Cell* generator){
+int CellComplex::coreduction(Cell* generator){
int coreductions = 0;
@@ -344,10 +348,9 @@ int CellComplex::coreductionMrozek(Cell* generator){
Q.push(generator);
Qset.insert(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;
@@ -377,7 +380,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
-
+/*
int CellComplex::coreduction(int dim){
if(dim < 0 || dim > 2) return 0;
@@ -405,6 +408,51 @@ int CellComplex::coreduction(int dim){
return count;
}
+*/
+/*
+int CellComplex::reduction(Cell* generator){
+
+ int count = 0;
+
+ std::queue<Cell*> Q;
+ std::set<Cell*, Less_Cell> Qset;
+
+ Q.push(generator);
+ Qset.insert(generator);
+
+
+ std::list<Cell*> cbd_s;
+ std::list<Cell*> bd_c;
+ Cell* s;
+ int round = 0;
+ while( !Q.empty() ){
+ round++;
+ //printf("%d ", round);
+
+ s = Q.front();
+ Q.pop();
+ removeCellQset(s, Qset);
+
+ cbd_s = s->getCoboundary();
+ if( cbd_s.size() == 1 && inSameDomain(s, cbd_s.front()) ){
+ removeCell(s);
+ bd_c = cbd_s.front()->getBoundary();
+ enqueueCells(bd_c, Q, Qset);
+ removeCell(cbd_s.front());
+ count++;
+ }
+ else if(cbd_s.empty()){
+ bd_c = s->getBoundary();
+ enqueueCells(bd_c, Q, Qset);
+ }
+
+
+ }
+
+ return count;
+}
+*/
+
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
@@ -431,18 +479,58 @@ int CellComplex::reduction(int dim){
}
-void CellComplex::reduceComplex(){
+void CellComplex::reduceComplex(bool omitHighdim){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
- reduction(3);
- reduction(2);
- reduction(1);
+
+
+ //reduction(3);
+ //reduction(2);
+ //reduction(1);
+
+
+
+ int count = 0;
+ for(int i = 3; i > 0; i--) count = count + reduction(i);
+
+
+
+ if(count == 0 && omitHighdim){
+
+ removeSubdomain();
+ std::set<Cell*, Less_Cell> generatorCells;
+
+ while (getSize(getDim()) != 0){
+ citer cit = firstCell(getDim());
+ Cell* cell = *cit;
+ generatorCells.insert(cell);
+ removeCell(cell);
+ reduction(3);
+ reduction(2);
+ reduction(1);
+
+ }
+
+
+ for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
+ Cell* cell = *cit;
+ cell->clearBoundary();
+ cell->clearCoboundary();
+ _cells[getDim()].insert(cell);
+ }
+
+
+ }
+
+
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> generatorCells[4];
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -456,31 +544,194 @@ void CellComplex::coreduceComplex(int generatorDim){
cell = *cit;
cit++;
}
- generatorCells[i].insert(cell);
+ if(!cell->inSubdomain()) generatorCells[i].insert(cell);
removeCell(cell);
- //coreduction(0);
- //coreduction(1);
- //coreduction(2);
- coreductionMrozek(cell);
+ coreduction(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);
+ _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::removeSubdomain(){
+
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()) {
+ removeCell(cell);
+ cit = firstCell(i);
+ }
+ }
+
+ }
+ return;
+}
+
+/*
+int CellComplex::coreduction(int dim){
+
+ int count = 0;
+
+ std::queue<Cell*> Q;
+ std::set<Cell*, Less_Cell> Qset;
+
+ std::list<Cell*> bd_s;
+ std::list<Cell*> cbd_c;
+
+ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+
+ Cell* cell = *cit;
+ cbd_c = cell->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+
+ int round = 0;
+ while( !Q.empty() ){
+ round++;
+ //printf("%d ", round);
+
+ Cell* s = Q.front();
+ Q.pop();
+
+ 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, Qset);
+ removeCell(bd_s.front());
+ cit = firstCell(dim);
+ count++;
+ }
+ else if(bd_s.empty()){
+ cbd_c = s->getCoboundary();
+ enqueueCells(cbd_c, Q, Qset);
+ }
+
+ removeCellQset(s, Qset);
+
+ }
+
+
+ }
+ return count;
+}
+*/
+void CellComplex::coreduceComplex(bool omitLowdim){
+
+ printf("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;
+
+ removeSubdomain();
+
+
+ for(int dim = 0; dim < 4; dim++){
+ citer cit = firstCell(dim);
+ if(cit != lastCell(dim)){
+ Cell* cell = *cit;
+ count = count + coreduction(cell);
+ }
+ }
+
+ if(count == 0 && omitLowdim){
+ std::set<Cell*, Less_Cell> generatorCells;
+
+ while (getSize(0) != 0){
+ citer cit = firstCell(0);
+ Cell* cell = *cit;
+ generatorCells.insert(cell);
+ removeCell(cell);
+ coreduction(cell);
+ }
+
+ for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
+ Cell* cell = *cit;
+ cell->clearBoundary();
+ cell->clearCoboundary();
+ _cells[0].insert(cell);
+ }
+
+
+ }
+ /*
+ std::set<Cell*, Less_Cell> generatorCells;
+
+ while (getSize(0) != 0){
+ citer cit = firstCell(0);
+ Cell* cell = *cit;
+ while(!cell->inSubdomain() && cit != lastCell(0)){
+ cell = *cit;
+ cit++;
+ }
+ if(!cell->inSubdomain()) generatorCells.insert(cell);
+ removeCell(cell);
+
+ coreduction(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::computeBettiNumbers(){
+
+ /*
+ removeSubdomain();
+
+ int count = 0;
+ for(int dim = 0; dim < 4; dim++){
+ citer cit = firstCell(dim);
+ if(cit != lastCell(dim)){
+ Cell* cell = *cit;
+ count = count + coreduction(cell);
+ }
+ }
+ */
+ for(int i = 0; i < 4; i++){
+ printf("Betti number computation process: step %d of 4 \n", i+1);
+
+
+ while (getSize(i) != 0){
+
+ citer cit = firstCell(i);
+ Cell* cell = *cit;
+ while(!cell->inSubdomain() && cit != lastCell(i)){
+ cell = *cit;
+ cit++;
+ }
+ if(!cell->inSubdomain()) _betti[i] = _betti[i] + 1;
+ removeCell(cell);
+
+ coreduction(cell);
+ }
+
+ }
+ printf("Cell complex Betti numbers: \n b0 = %d \n b1 = %d \n b2 = %d \n b3 = %d \n",
+ getBettiNumber(0), getBettiNumber(1), getBettiNumber(2), getBettiNumber(3));
+
+ return;
+}
+
+
void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co){
@@ -551,7 +802,7 @@ 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;
+ if(dim < 0 || dim > 2) return 0;
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
@@ -662,7 +913,7 @@ int CellComplex::combine(int dim){
return count;
}
-
+/*
void CellComplex::repairComplex(int i){
printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -695,7 +946,8 @@ void CellComplex::repairComplex(int i){
return;
}
-
+*/
+
void CellComplex::swapSubdomain(){
for(int i = 0; i < 4; i++){
@@ -828,3 +1080,12 @@ void CellComplex::printComplex(int dim){
+DualCellComplex::DualCellComplex(CellComplex* cellComplex){
+
+ for(int i = 0; i < 4; i++){
+ _cells[i] = cellComplex->getCells(i);
+ }
+
+
+
+}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 9ab40aadf5..5eec690c5a 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -116,7 +116,18 @@ class Cell
if(*cell2 == *cell) { _coboundary.erase(it); break; }
}
}
-
+
+ virtual void clearBoundary() { _boundary.clear(); }
+ virtual void clearCoboundary() { _coboundary.clear(); }
+
+ virtual void makeDualCell(){
+ std::list< std::pair<int, Cell*> > temp = _boundary;
+ _boundary = _coboundary;
+ _coboundary = temp;
+ _dim = 3-_dim;
+
+ }
+
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);
@@ -215,7 +226,7 @@ class ZeroSimplex : public Simplex
}
~ZeroSimplex(){}
- int getDim() const { return 0; }
+ //int getDim() const { return 0; }
int getNumVertices() const { return 1; }
int getVertex(int vertex) const {return _v; }
int getSortedVertex(int vertex) const {return _v; }
@@ -268,7 +279,7 @@ class OneSimplex : public Simplex
~OneSimplex(){}
- int getDim() const { return 1; }
+ //int getDim() const { return 1; }
int getNumVertices() const { return 2; }
int getNumFacets() const { return 2; }
int getVertex(int vertex) const {return _v[vertex]; }
@@ -335,7 +346,7 @@ class TwoSimplex : public Simplex
~TwoSimplex(){}
- int getDim() const { return 2; }
+ //int getDim() const { return 2; }
int getNumVertices() const { return 3; }
int getNumFacets() const { return 3; }
int getVertex(int vertex) const {return _v[vertex]; }
@@ -405,7 +416,7 @@ class ThreeSimplex : public Simplex
~ThreeSimplex(){}
- int getDim() const { return 3; }
+ //int getDim() const { return 3; }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
int getVertex(int vertex) const {return _v[vertex]; }
@@ -466,7 +477,7 @@ class Quadrangle : public Cell
}
~Quadrangle(){}
- int getDim() const { return 2; }
+ //int getDim() const { return 2; }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
int getVertex(int vertex) const {return _v[vertex]; }
@@ -650,7 +661,7 @@ class CombinedCell : public Cell{
// A class representing a cell complex made out of a finite element mesh.
class CellComplex
{
- private:
+ protected:
// the domain in the model which this cell complex covers
std::vector<GEntity*> _domain;
@@ -666,7 +677,12 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
- std::set<Cell*, Less_Cell> _originalCells[4];
+ //std::set<Cell*, Less_Cell> _originalCells[4];
+
+ // Betti numbers of this cell complex (ranks of homology groups)
+ int _betti[4];
+
+ int _dim;
public:
// iterator for the cells of same dimension
@@ -692,13 +708,40 @@ class CellComplex
_cells[cell->getDim()].insert(cell);
}
}
+ /*
+ CellComplex(CellComplex* cellComplex){
+
+ _domain = cellComplex->_domain;
+ _subdomain = cellComplex->_subdomain;
+ _boundary = cellComplex->_boundary;
+ _domainVertices = cellComplex->_domainVertices;
+
+ for(int i = 0; i < 4; i++){
+ _betti[i] = cellComplex->_betti[i];
+
+ for(citer cit = cellComplex->_cells[i].begin(); cit != cellComplex->_cells[i].end(); cit++){
+ Cell* cell = *cit;
+ if(i == 0) _cells[i].insert(new ZeroSimplex(*cell));
+
+ }
+
+ _originalCells[i] = _cells[i];
+ }
+
+ _dim = cellComplex->_dim;
+
+ }*/
+
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
+ CellComplex(){}
virtual ~CellComplex(){}
-
+
// get the number of certain dimensional cells
virtual int getSize(int dim){ return _cells[dim].size(); }
+ virtual int getDim() {return _dim; }
+
virtual std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
// iterators to the first and last cells of certain dimension
@@ -706,8 +749,8 @@ class CellComplex
virtual citer lastCell(int dim) {return _cells[dim].end(); }
// find a cell in this cell complex
- virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices, bool original=false);
- virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
+ //virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices, bool original=false);
+ //virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
@@ -726,27 +769,33 @@ 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 coreduction(int dim);
+ //virtual int coreduction();
// stores removed cells
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
virtual int reduction(int dim);
+ //virtual int reduction(Cell* generator);
// useful functions for (co)reduction of cell complex
- virtual void reduceComplex();
- // coreduction up to generators of dimension generatorDim
- virtual void coreduceComplex(int generatorDim=3);
+ virtual void reduceComplex(bool omitHighdim = false);
+ virtual void coreduceComplex(bool omitlowDim = false);
+
+ //virtual void coreduceComplex(int generatorDim);
+ //virtual void coreduceComplex();
// queued coreduction presented in Mrozek's paper
// slower, but produces cleaner result
- virtual int coreductionMrozek(Cell* generator);
+ virtual int coreduction(Cell* generator);
// add every volume, face and edge its missing boundary cells
- virtual void repairComplex(int i=3);
+ //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();
+ virtual void removeSubdomain();
+
// print the vertices of cells of certain dimension
virtual void printComplex(int dim);
@@ -757,6 +806,36 @@ class CellComplex
virtual int combine(int dim);
virtual int cocombine(int dim);
+ virtual void computeBettiNumbers();
+ virtual int getBettiNumber(int i) { if(i > -1 && i < 4) return _betti[i]; else return 0; }
+
+ virtual 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();
+ }
+ }
+ }
+
+};
+
+
+class DualCellComplex : public CellComplex
+{
+
+ public:
+
+ DualCellComplex(CellComplex* cellComplex);
+ ~DualCellComplex(){}
+
};
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 88d02b75d2..365599a630 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -11,17 +11,21 @@
ChainComplex::ChainComplex(CellComplex* cellComplex){
- _dim = 0;
- _HMatrix[0] = NULL;
- _kerH[0] = NULL;
- _codH[0] = NULL;
- _JMatrix[0] = NULL;
- _QMatrix[0] = NULL;
- _Hbasis[0] = NULL;
-
- for(int dim = 1; dim < 4; dim++){
+ _dim = cellComplex->getDim();
+ _cellComplex = cellComplex;
+
+ _HMatrix[4] = NULL;
+ _kerH[4] = NULL;
+ _codH[4] = NULL;
+ _JMatrix[4] = NULL;
+ _QMatrix[4] = NULL;
+ _Hbasis[4] = NULL;
+
+
+ for(int dim = 0; dim < 4; dim++){
unsigned int cols = cellComplex->getSize(dim);
- unsigned int rows = cellComplex->getSize(dim-1);
+ unsigned int rows = 0;
+ if(dim > 0) rows = cellComplex->getSize(dim-1);
int index = 1;
@@ -36,17 +40,18 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
}
index = 1;
- for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){
- Cell* cell = *cit;
- cell->setIndex(index);
- index++;
- if(cell->inSubdomain()){
- index--;
- rows--;
+ if(dim > 0){
+ for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){
+ Cell* cell = *cit;
+ cell->setIndex(index);
+ index++;
+ if(cell->inSubdomain()){
+ index--;
+ rows--;
+ }
}
}
-
if( cols == 0 ){
//_HMatrix[dim] = create_gmp_matrix_zero(rows, 1);
_HMatrix[dim] = NULL;
@@ -62,6 +67,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
mpz_init(elem);
_HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
+ //printMatrix(_HMatrix[dim]);
for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if(!cell->inSubdomain()){
@@ -70,6 +76,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
Cell* bdCell = (*it).second;
if(!bdCell->inSubdomain()){
int old_elem = 0;
+ //printf("kaka1: %d, kaka2: %d \n", bdCell->getIndex(), cell->getIndex());
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);
@@ -126,8 +133,6 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
_QMatrix[dim] = NULL;
_Hbasis[dim] = NULL;
- if(cellComplex->getSize(dim) != 0) _dim = dim;
-
}
return;
@@ -254,7 +259,7 @@ void ChainComplex::Inclusion(int lowDim, int highDim){
void ChainComplex::Quotient(int dim){
- if(dim < 0 || dim > 3 || _JMatrix[dim] == NULL) return;
+ if(dim < 0 || dim > 4 || _JMatrix[dim] == NULL) return;
gmp_matrix* JMatrix = copy_gmp_matrix(_JMatrix[dim], 1, 1,
gmp_matrix_rows(_JMatrix[dim]), gmp_matrix_cols(_JMatrix[dim]));
@@ -297,35 +302,43 @@ void ChainComplex::computeHomology(bool dual){
int lowDim = 0;
int highDim = 0;
+ int setDim = 0;
+
- for(int i=0; i < 4; i++){
+ for(int i=-1; i < 4; i++){
if(dual){
lowDim = getDim()+1-i;
highDim = getDim()+1-(i+1);
+ setDim = highDim;
//KerCod(lowDim);
}
else{
lowDim = i;
highDim = i+1;
+ setDim = lowDim;
//KerCod(highDim);
}
- printf("Homology computation process: step %d of 4 \n", i+1 );
+ printf("Homology computation process: step %d of 4 \n", i+1);
KerCod(highDim);
// 1) no edges, but zero cells
- if(lowDim == 0 && gmp_matrix_cols(getHMatrix(lowDim)) > 0 && getHMatrix(highDim) == NULL) {
- setHbasis( lowDim, create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) );
- }
+ if(lowDim == 0 && !dual && gmp_matrix_cols(getHMatrix(lowDim)) > 0 && getHMatrix(highDim) == NULL) {
+ setHbasis( setDim, create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) );
+ }
+ else if(highDim == 0 && dual && gmp_matrix_rows(getHMatrix(highDim)) > 0 && getHMatrix(lowDim) == NULL) {
+ setHbasis( setDim, create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
+ }
+
// 2) this dimension is empty
else if(getHMatrix(lowDim) == NULL && getHMatrix(highDim) == NULL){
- setHbasis(lowDim, NULL);
+ setHbasis(setDim, NULL);
}
// 3) No higher dimension cells -> none of the cycles are boundaries
else if(getHMatrix(highDim) == NULL){
- setHbasis( lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ setHbasis( setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
gmp_matrix_rows(getKerHMatrix(lowDim)),
gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
@@ -344,19 +357,19 @@ void ChainComplex::computeHomology(bool dual){
Quotient(lowDim);
if(getCodHMatrix(highDim) == NULL){
- setHbasis(lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ setHbasis(setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
gmp_matrix_rows(getKerHMatrix(lowDim)),
gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
else if(getJMatrix(lowDim) == NULL || getQMatrix(lowDim) == NULL){
- setHbasis(lowDim, NULL);
+ setHbasis(setDim, NULL);
}
else{
- setHbasis(lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ setHbasis(setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
gmp_matrix_rows(getKerHMatrix(lowDim)),
gmp_matrix_cols(getKerHMatrix(lowDim))) );
- gmp_matrix_right_mult(getHbasis(lowDim), getQMatrix(lowDim));
+ gmp_matrix_right_mult(getHbasis(setDim), getQMatrix(lowDim));
}
}
@@ -401,7 +414,7 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){
std::vector<int> coeffVector;
- if(dim < 0 || dim > 3) return coeffVector;
+ if(dim < 0 || dim > 4) return coeffVector;
if(_Hbasis[dim] == NULL || gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return coeffVector;
int rows = gmp_matrix_rows(_Hbasis[dim]);
@@ -435,6 +448,7 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
if(coeffs.at(i) != 0) _cells.push_back( std::make_pair(cell, coeffs.at(i)) );
i++;
}
+
}
_name = name;
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index adc5450660..b8a5e2ec57 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -36,28 +36,29 @@ class ChainComplex{
private:
// boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
- gmp_matrix* _HMatrix[4];
+ gmp_matrix* _HMatrix[5];
// Basis matrices for the kernels and codomains of the boundary operator matrices
- gmp_matrix* _kerH[4];
- gmp_matrix* _codH[4];
+ gmp_matrix* _kerH[5];
+ gmp_matrix* _codH[5];
- gmp_matrix* _JMatrix[4];
- gmp_matrix* _QMatrix[4];
- std::vector<long int> _torsion[4];
+ gmp_matrix* _JMatrix[5];
+ gmp_matrix* _QMatrix[5];
+ std::vector<long int> _torsion[5];
// bases for homology groups
- gmp_matrix* _Hbasis[4];
+ gmp_matrix* _Hbasis[5];
int _dim;
+ CellComplex* _cellComplex;
// set the matrices
- virtual void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _HMatrix[dim] = matrix;}
- virtual void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _kerH[dim] = matrix;}
- virtual void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _codH[dim] = matrix;}
- virtual void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _JMatrix[dim] = matrix;}
- virtual void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _QMatrix[dim] = matrix;}
- virtual void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _Hbasis[dim] = matrix;}
+ virtual void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
+ virtual void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
+ virtual void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _codH[dim] = matrix;}
+ virtual void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
+ virtual void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
+ virtual void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
@@ -66,7 +67,7 @@ class ChainComplex{
ChainComplex(CellComplex* cellComplex);
ChainComplex(){
- for(int i = 0; i < 4; i++){
+ for(int i = 0; i < 5; i++){
_HMatrix[i] = create_gmp_matrix_zero(1,1);
_kerH[i] = NULL;
_codH[i] = NULL;
@@ -81,12 +82,12 @@ class ChainComplex{
virtual int getDim() { return _dim; }
// get the boundary operator matrix dim->dim-1
- virtual gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 4) return _HMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 4) return _kerH[dim]; else return NULL;}
- virtual gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 4) return _codH[dim]; else return NULL;}
- virtual gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 4) return _JMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 4) return _QMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 4) return _Hbasis[dim]; else return NULL;}
+ virtual gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
+ virtual gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
+ virtual gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
// Compute basis for kernel and codomain of boundary operator matrix of dimension dim (ie. ker(h_dim) and cod(h_dim) )
virtual void KerCod(int dim);
@@ -98,14 +99,14 @@ class ChainComplex{
virtual void Quotient(int dim);
// transpose the boundary operator matrices, these are boundary operator matrices for the dual mesh
- virtual void transposeHMatrices() { for(int i = 0; i < 4; i++) gmp_matrix_transp(_HMatrix[i]); }
- virtual void transposeHMatrix(int dim) { if(dim > -1 && dim < 4) gmp_matrix_transp(_HMatrix[dim]); }
+ virtual void transposeHMatrices() { for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
+ virtual void transposeHMatrix(int dim) { if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
// Compute bases for the homology groups of this chain complex
virtual void computeHomology(bool dual=false);
virtual std::vector<int> getCoeffVector(int dim, int chainNumber);
- virtual int getBasisSize(int dim) { if(dim > -1 && dim < 4) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
+ virtual int getBasisSize(int dim) { if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
virtual int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index 5bc4920a9a..4fe82d739b 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -66,7 +66,7 @@ int main(int argc, char **argv)
complex->getSize(3), complex->getSize(2), complex->getSize(1), complex->getSize(0));
// reduce the complex in order to ease homology computation
- complex->reduceComplex();
+ complex->reduceComplex(true);
// perform series of cell combinatios in order to further ease homology computation
complex->combine(3);
@@ -104,9 +104,10 @@ int main(int argc, char **argv)
CellComplex* complex2 = new CellComplex(domain, subdomain);
// reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
- complex2->coreduceComplex(0);
+ complex2->coreduceComplex(true);
// perform series of "dual complex" cell combinations in order to ease homology computation
+ complex2->cocombine(0);
complex2->cocombine(1);
complex2->cocombine(2);
@@ -127,10 +128,10 @@ int main(int argc, char **argv)
std::string generator;
std::string dimension;
convert(i, generator);
- convert(3-(j-1), dimension);
+ convert(3-j, dimension);
std::string name = dimension + "D Dual cut " + generator;
- Chain* chain = new Chain(complex2->getCells(j-1), dualChains->getCoeffVector(j,i), complex2, name);
+ Chain* chain = new Chain(complex2->getCells(j), dualChains->getCoeffVector(j,i), complex2, name);
chain->writeChainMSH("chains.msh");
delete chain;
}
--
GitLab