diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 02a7c12836dfee280ec786c7b23009521823e175..892fc24c667f1169656bd9f764d670802775495e 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 9ab40aadf53569055f1cfae0350fc2344a4f847f..5eec690c5ac830d2c988829a3f7a4257f16a0800 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 88d02b75d2b2d90fd9c8d92a75e978465d114229..365599a63012258e51505a4e7c7f1d9fab807292 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 adc5450660cee15d46c820ecc57e583997d06164..b8a5e2ec577dde443fd530dd2c68492dd33e94fa 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 5bc4920a9ad710ffbfdcad3786a084cd8ef9e92e..4fe82d739b793405b77af7a7f45865dcbde804a7 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;
}