diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index ce80a05b9b53d6fe7390286a4dec9bf3b7174ad1..3befc1ef002de272d24932a7998af410d0a5bb40 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -7,8 +7,6 @@
#include "CellComplex.h"
-#if defined(HAVE_KBIPACK)
-
CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain ){
_domain = domain;
@@ -65,6 +63,10 @@ void CellComplex::insertCells(bool subdomain){
}
+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;
@@ -184,9 +186,9 @@ std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau
}
-void CellComplex::removeCell(Cell* cell){
+void CellComplex::removeCell(Cell* cell, bool deleteCell){
_cells[cell->getDim()].erase(cell);
- delete cell;
+ if(deleteCell) delete cell;
}
void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
Cell* c = *cell;
@@ -197,7 +199,7 @@ void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
Qset.erase(cell);
- delete cell;
+ delete cell;
}
void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
@@ -272,7 +274,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
}
-int CellComplex::coreduction(int dim){
+int CellComplex::coreduction(int dim, bool deleteCells){
if(dim < 0 || dim > 2) return 0;
@@ -287,8 +289,8 @@ int CellComplex::coreduction(int dim){
bd_c = bd(cell);
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
- removeCell(cell);
- removeCell(bd_c.at(0));
+ removeCell(cell, deleteCells);
+ removeCell(bd_c.at(0), deleteCells);
count++;
coreduced = true;
}
@@ -300,8 +302,39 @@ int CellComplex::coreduction(int dim){
}
+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, false);
+ removeCell(bd_c.at(0), false);
+ count++;
+ coreduced = true;
+ }
+
+ }
+ }
+
+ return count;
+
+}
-int CellComplex::reduction(int dim){
+
+
+int CellComplex::reduction(int dim, bool deleteCells){
if(dim < 1 || dim > 3) return 0;
std::vector<Cell*> cbd_c;
@@ -314,8 +347,8 @@ int CellComplex::reduction(int dim){
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));
+ removeCell(cell, deleteCells);
+ removeCell(cbd_c.at(0), deleteCells);
count++;
reduced = true;
}
@@ -334,36 +367,79 @@ void CellComplex::reduceComplex(){
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, false);
+ 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(){
- std::set<Cell*, Less_Cell> removedCells[4];
+ 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));
- for(int i = 0; i < 3; i++){
+ for(int i = 0; i < 4; i++){
while (getSize(i) != 0){
citer cit = firstCell(i);
Cell* cell = *cit;
- removedCells[i].insert(cell);
- _cells[i].erase(cell);
- //removeCell(cell);
- //complex.coreductionMrozek(kaka);
+ while(!cell->inSubdomain() && cit != lastCell(i)){
+ cell = *cit;
+ cit++;
+ }
+ generatorCells[i].insert(cell);
+ removeCell(cell, false);
coreduction(0);
coreduction(1);
coreduction(2);
}
}
- for(int i = 0; i < 3; i++){
- for(citer cit = removedCells[i].begin(); cit != removedCells[i].end(); cit++){
+ for(int i = 0; i < 4; i++){
+ for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
Cell* cell = *cit;
- _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;
}
+
+#if defined(HAVE_KBIPACK)
+
std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
// h_dim: C_dim -> C_(dim-1)
@@ -418,7 +494,7 @@ std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
}
return HMatrix;
}
-
+#endif
void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
@@ -541,6 +617,8 @@ void CellComplex::printComplex(int dim, bool subcomplex){
}
}
+#if defined(HAVE_KBIPACK)
+
void ChainComplex::KerCod(int dim){
if(dim < 1 || dim > 3 || _HMatrix[dim] == NULL) return;
@@ -580,11 +658,16 @@ void ChainComplex::KerCod(int dim){
return;
}
-/*
-//j:B->Z
+
+//j:B_(k+1)->Z_k
void ChainComplex::Inclusion(int dim){
- if(dim < 1 || dim > 3 || _kerH[dim] == NULL || _codH[dim] == NULL) return;
+ if(dim < 1 || dim > 2 || _kerH[dim] == NULL || _codH[dim+1] == NULL) return;
+
+ gmp_matrix* Zbasis = copy_gmp_matrix(_kerH[dim], 1, 1,
+ gmp_matrix_rows(_kerH[dim]), gmp_matrix_cols(_kerH[dim]));
+ gmp_matrix* Bbasis = copy_gmp_matrix(_codH[dim+1], 1, 1,
+ gmp_matrix_rows(_codH[dim+1]), gmp_matrix_cols(_codH[dim+1]));
int rows = gmp_matrix_rows(Zbasis);
int cols = gmp_matrix_cols(Zbasis);
@@ -594,18 +677,13 @@ void ChainComplex::Inclusion(int dim){
cols = gmp_matrix_cols(Bbasis);
if(rows < cols) return;
- gmp_matrix* Zbasis = copy_gmp_matrix(_ker[dim], 1, 1,
- gmp_matrix_rows(_kerH[dim]), gmp_matrix_cols(_kerH[dim]));
- gmp_matrix* Bbasis = copy_gmp_matrix(_cod[dim], 1, 1,
- gmp_matrix_rows(_codH[dim]), gmp_matrix_cols(_codH[dim]));
-
// A*inv(V) = U*S
- normalForm = create_gmp_Smith_normal_form(Zbasis, NOT_INVERTED, INVERTED);
+ gmp_normal_form* normalForm = create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED);
mpz_t elem;
mpz_init(elem);
- for(int i = 1; i < cols; i++){
+ for(int i = 1; i <= cols; i++){
gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
if(mpz_cmp_si(elem,0) == 0){
destroy_gmp_normal_form(normalForm);
@@ -615,10 +693,146 @@ void ChainComplex::Inclusion(int dim){
gmp_matrix_left_mult(normalForm->left, Bbasis);
+
+
+ gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1, gmp_matrix_cols(Zbasis), gmp_matrix_cols(Bbasis));
+ destroy_gmp_matrix(Bbasis);
+
+ rows = gmp_matrix_rows(LB);
+ cols = gmp_matrix_cols(LB);
+
+ mpz_t divisor;
+ mpz_init(divisor);
+ mpz_t remainder;
+ mpz_init(remainder);
+ mpz_t result;
+ mpz_init(result);
+
+ for(int i = 1; i <= rows; i++){
+ gmp_matrix_get_elem(divisor, i, i, normalForm->canonical);
+ for(int j = 1; j <= cols; j++){
+ gmp_matrix_get_elem(elem, i, j, LB);
+ mpz_cdiv_qr(result, remainder, elem, divisor);
+ if(mpz_cmp_si(remainder, 0) == 0){
+ gmp_matrix_set_elem(result, i, j, LB);
+ }
+ else return;
+ }
+ }
+
+ gmp_matrix_left_mult(normalForm->right, LB);
+
+ _JMatrix[dim] = LB;
+
+ mpz_clear(elem);
+ mpz_clear(divisor);
+ mpz_clear(result);
+ destroy_gmp_normal_form(normalForm);
+
+ return;
+
+}
+void ChainComplex::Quotient(int dim){
+
+ if(dim < 1 || dim > 3 || _JMatrix[dim] == NULL) return;
+
+ gmp_matrix* JMatrix = copy_gmp_matrix(_JMatrix[dim], 1, 1,
+ gmp_matrix_rows(_JMatrix[dim]), gmp_matrix_cols(_JMatrix[dim]));
+ int rows = gmp_matrix_rows(JMatrix);
+ int cols = gmp_matrix_cols(JMatrix);
+
+ gmp_normal_form* normalForm = create_gmp_Smith_normal_form(JMatrix, NOT_INVERTED, NOT_INVERTED);
+
+ 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){
+ destroy_gmp_normal_form(normalForm);
+ return;
+ }
+ if(mpz_cmp_si(elem,1) > 0){
+ _torsion[dim].push_back(mpz_get_si(elem));
+ }
+ }
+
+ int rank = cols - _torsion[dim].size();
+ if(rows - rank > 0){
+ gmp_matrix* Hbasis = copy_gmp_matrix(normalForm->left, 1, rank+1, rows, rows);
+ _QMatrix[dim] = Hbasis;
+ }
+
+ mpz_clear(elem);
+ destroy_gmp_normal_form(normalForm);
+ return;
+}
+
+void ChainComplex::computeHomology(){
+
+ for(int i=0; i < 4; i++){
+
+ KerCod(i+1);
+
+ // 1) no edges, but zero cells
+ if(i == 0 && gmp_matrix_cols(getHMatrix(0)) > 0 && getHMatrix(1) == NULL) {
+ _Hbasis[i] = create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(0)));
+ }
+
+ // 2) this dimension is empty
+ else if(getHMatrix(i) == NULL && getHMatrix(i+1) == NULL){
+ _Hbasis[i] = NULL;
+ }
+ // 3) No higher dimension cells -> none of the cycles are boundaries
+ else if(getHMatrix(i+1) == NULL){
+ _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+ }
+
+ // 5) General case:
+ // 1) Find the bases of boundaries B and cycles Z
+ // 2) find j: B -> Z and
+ // 3) find quotient Z/j(B)
+ else {
+
+ // 4) No lower dimension cells -> all chains are cycles
+ if(getHMatrix(i) == NULL){
+ _kerH[i] = create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(i+1)));
+ }
+
+ Inclusion(i);
+ Quotient(i);
+
+ if(getCodHMatrix(i+1) == NULL){
+ _Hbasis[i] = getKerHMatrix(i);
+ }
+ else if(getJMatrix(i) == NULL || getQMatrix(i) == NULL){
+ _Hbasis[i] = NULL;
+ }
+ else{
+ _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+ gmp_matrix_right_mult(_Hbasis[i], getQMatrix(i));
+ }
+ }
+
+ destroy_gmp_matrix(_kerH[i]);
+ destroy_gmp_matrix(_codH[i]);
+ destroy_gmp_matrix(_JMatrix[i]);
+ destroy_gmp_matrix(_QMatrix[i]);
+
+ _kerH[i] = NULL;
+ _codH[i] = NULL;
+ _JMatrix[i] = NULL;
+ _QMatrix[i] = NULL;
+
+ }
+
+
+ return;
}
-*/
+
void ChainComplex::matrixTest(){
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 5f73d241d72811ffdca7ae9cbbb2907a490bba0b..dbcfc66cf42026a2564130f55f203c45601c8f77 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -28,6 +28,8 @@ extern "C" {
#include "gmp_normal_form.h" // perhaps make c++ headers instead?
}
+#endif
+
// Abstract class representing a cell of a cell complex.
class Cell
{
@@ -364,9 +366,12 @@ class CellComplex
virtual void insertCells(bool subdomain);
public:
- CellComplex( std::set<Cell*, Less_Cell>* cells ) {
- for(int i = 0; i < 4; i++){
- //_cells[i] = cells[i];
+ CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
+ _domain = domain;
+ _subdomain = subdomain;
+ for(citer cit = cells.begin(); cit != cells.end(); cit++){
+ Cell* cell = *cit;
+ _cells[cell->getDim()].insert(cell);
}
}
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
@@ -396,24 +401,32 @@ class CellComplex
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 removeCell(Cell* cell, bool deleteCell=true);
virtual void removeCellIt(std::set<Cell*, Less_Cell>::iterator cell);
+ virtual void insertCell(Cell* cell);
+
// check whether two cells both belong to subdomain or if neither one does
virtual bool inSameDomain(Cell* c1, Cell* c2) const { return
( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
// 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, bool deleteCells=true);
+ // 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 reduction(int dim, bool deleteCells=true);
// useful functions for (co)reduction of cell complex
virtual void reduceComplex();
virtual void coreduceComplex();
+ // 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);
@@ -429,15 +442,19 @@ class CellComplex
// get the boundary operator matrix dim->dim-1
//virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+
+#if defined(HAVE_KBIPACK)
// construct boundary operator matrices of this cell complex
// used to construct a chain complex
virtual std::vector<gmp_matrix*> constructHMatrices();
+#endif
};
+#if defined(HAVE_KBIPACK)
// A class representing a chain complex of a cell complex.
// This should only be constructed for a reduced cell complex because of
-// dense matrix reprentations and great computational complexity in its methods.
+// dense matrix representations and great computational complexity in its methods.
class ChainComplex{
private:
// boundary operator matrices for this chain complex
@@ -450,7 +467,10 @@ class ChainComplex{
gmp_matrix* _JMatrix[4];
gmp_matrix* _QMatrix[4];
- std::vector<int> _torsion[4];
+ std::vector<long int> _torsion[4];
+
+ // bases for homology groups
+ gmp_matrix* _Hbasis[4];
public:
@@ -463,6 +483,7 @@ class ChainComplex{
_codH[i] = NULL;
_JMatrix[i] = NULL;
_QMatrix[i] = NULL;
+ _Hbasis[i] = NULL;
}
}
@@ -473,25 +494,30 @@ class ChainComplex{
_codH[i] = NULL;
_JMatrix[i] = NULL;
_QMatrix[i] = NULL;
+ _Hbasis[i] = NULL;
}
}
virtual ~ChainComplex(){}
// get the boundary operator matrix dim->dim-1
- virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
- virtual gmp_matrix* getKerHMatrix(int dim) { return _kerH[dim]; }
- virtual gmp_matrix* getCodHMatrix(int dim) { return _codH[dim]; }
- virtual gmp_matrix* getJMatrix(int dim) { return _JMatrix[dim]; }
- virtual gmp_matrix* getQMatrix(int dim) { return _QMatrix[dim]; }
+ 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;}
// 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);
// Compute matrix representation J for inclusion relation from dim-cells who are boundary of dim+1-cells
// to cycles of dim-cells (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
- virtual void Inclusion(int dim){}
+ virtual void Inclusion(int dim);
// Compute quotient problem for given inclusion relation j to find representatives of homology groups
// and possible torsion coeffcients
- virtual void Quotient(int dim){}
+ virtual void Quotient(int dim);
+
+ // Compute bases for the homology groups of this chain complex
+ virtual void computeHomology();
virtual int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index 2d07724b0804ac8f598a68b0ca9fc37842758892..c6dcb11fac8ecec18015404d31fddd310f20c005 100755
--- a/contrib/kbipack/gmp_matrix.c
+++ b/contrib/kbipack/gmp_matrix.c
@@ -22,7 +22,7 @@
P.O.Box 692, FIN-33101 Tampere, Finland
saku.suuriniemi@tut.fi
- $Id: gmp_matrix.c,v 1.2 2009-04-03 11:06:12 matti Exp $
+ $Id: gmp_matrix.c,v 1.3 2009-04-14 10:02:22 matti Exp $
*/
@@ -281,6 +281,24 @@ gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col,
return EXIT_SUCCESS;
}
+/* (matrix(row, col)) <- elem */
+int
+gmp_matrix_set_elem(mpz_t elem, size_t row, size_t col,
+ const gmp_matrix * m)
+{
+#ifdef DEBUG
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+#endif
+
+ mpz_set(m -> storage[(col-1)*(m -> rows)+row-1], elem);
+ return EXIT_SUCCESS;
+}
+
int
gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix * m)
{
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
index 784a87c20609f436bacda294851d6887107914b3..5d0e7c57dc853238729492548f1da24eb634b837 100755
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -22,7 +22,7 @@
P.O.Box 692, FIN-33101 Tampere, Finland
saku.suuriniemi@tut.fi
- $Id: gmp_matrix.h,v 1.2 2009-04-03 11:06:13 matti Exp $
+ $Id: gmp_matrix.h,v 1.3 2009-04-14 10:02:22 matti Exp $
*/
#ifndef __GMP_MATRIX_H__
@@ -72,6 +72,11 @@ int
gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col,
const gmp_matrix *);
+/* (matrix(row, col)) <- elem */
+int
+gmp_matrix_set_elem(mpz_t elem, size_t row, size_t col,
+ const gmp_matrix *);
+
int
gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix *);