diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 3befc1ef002de272d24932a7998af410d0a5bb40..5d17436d71154629fd7a492d10a75a2dddd3888f 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -403,13 +403,13 @@ void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>&
return;
}
-void CellComplex::coreduceComplex(){
+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",
getSize(3), getSize(2), getSize(1), getSize(0));
- for(int i = 0; i < 4; i++){
+ for(int i = 0; i < 4; i++){
while (getSize(i) != 0){
citer cit = firstCell(i);
Cell* cell = *cit;
@@ -423,6 +423,7 @@ void CellComplex::coreduceComplex(){
coreduction(1);
coreduction(2);
}
+ if(generatorDim == i) break;
}
for(int i = 0; i < 4; i++){
@@ -436,65 +437,6 @@ void CellComplex::coreduceComplex(){
return;
}
-
-
-#if defined(HAVE_KBIPACK)
-
-std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
-
- // h_dim: C_dim -> C_(dim-1)
-
- std::vector<gmp_matrix*> HMatrix;
-
- HMatrix.push_back(NULL);
-
- for(int dim = 1; dim < 4; dim++){
- unsigned int cols = _cells[dim].size();
- unsigned int rows = _cells[dim-1].size();
-
- for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
- Cell* cell = *cit;
- if(cell->inSubdomain()) cols--;
- }
- for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
- Cell* cell = *cit;
- if(cell->inSubdomain()) rows--;
- }
-
-
- if( cols == 0 ){
- HMatrix.push_back(NULL);
- }
- else if( rows == 0){
- HMatrix.push_back(create_gmp_matrix_zero(1, cols));
- }
- else{
- long int elems[rows*cols];
-
- citer high = firstCell(dim);
- citer low = firstCell(dim-1);
-
- unsigned int i = 0;
- while(i < rows*cols){
- while(low != lastCell(dim-1)){
- Cell* lowcell = *low;
- Cell* highcell = *high;
- if(!(highcell->inSubdomain() || lowcell->inSubdomain())){
- elems[i] = kappa(*high, *low);
- i++;
- }
- low++;
- }
- low = firstCell(dim-1);
- high++;
- }
- HMatrix.push_back(create_gmp_matrix_int(rows, cols, elems));
-
- }
- }
- return HMatrix;
-}
-#endif
void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
@@ -553,6 +495,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
FILE *fp = fopen(name.c_str(), "w");
if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str());
+ printf("Unable to open file.");
return 0;
}
@@ -617,242 +560,5 @@ 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;
-
- gmp_matrix* HMatrix = copy_gmp_matrix(_HMatrix[dim], 1, 1,
- gmp_matrix_rows(_HMatrix[dim]), gmp_matrix_cols(_HMatrix[dim]));
-
- gmp_normal_form* normalForm = create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED);
- //printMatrix(normalForm->left);
- //printMatrix(normalForm->canonical);
- //printMatrix(normalForm->right);
-
- int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical), gmp_matrix_cols(normalForm->canonical));
- int rank = 0;
- mpz_t elem;
- mpz_init(elem);
-
- while(rank < minRowCol){
- gmp_matrix_get_elem(elem, rank+1, rank+1, normalForm->canonical);
- if(mpz_cmp_si(elem,0) == 0) break;
- rank++;
- }
-
- if(rank != gmp_matrix_cols(normalForm->canonical)){
- _kerH[dim] = copy_gmp_matrix(normalForm->right, 1, rank+1,
- gmp_matrix_rows(normalForm->right), gmp_matrix_cols(normalForm->right));
- }
-
- if(rank > 0){
- _codH[dim] = copy_gmp_matrix(normalForm->canonical, 1, 1,
- gmp_matrix_rows(normalForm->canonical), rank);
- gmp_matrix_left_mult(normalForm->left, _codH[dim]);
- }
-
- mpz_clear(elem);
- destroy_gmp_normal_form(normalForm);
-
- return;
-}
-
-//j:B_(k+1)->Z_k
-void ChainComplex::Inclusion(int dim){
-
- 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);
- if(rows < cols) return;
-
- rows = gmp_matrix_rows(Bbasis);
- cols = gmp_matrix_cols(Bbasis);
- if(rows < cols) return;
-
- // A*inv(V) = U*S
- 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++){
- gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
- if(mpz_cmp_si(elem,0) == 0){
- destroy_gmp_normal_form(normalForm);
- return;
- }
- }
-
- 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(){
-
- int rows = 3;
- int cols = 6;
-
- long int elems[rows*cols];
- for(int i = 1; i<=rows*cols; i++) elems[i-1] = i;
-
- gmp_matrix* matrix = create_gmp_matrix_int(rows, cols, elems);
-
- gmp_matrix* copymatrix = copy_gmp_matrix(matrix, 3, 2, 3, 5);
-
- printMatrix(matrix);
- printMatrix(copymatrix);
-
- return;
-}
-
-
-#endif
-
-
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index dbcfc66cf42026a2564130f55f203c45601c8f77..9cf34a2a59f32cb51129f360572d3596efdd0a93 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -21,15 +21,6 @@
#include "GFace.h"
#include "GVertex.h"
-#if defined(HAVE_KBIPACK)
-
-#include "gmp.h"
-extern "C" {
- #include "gmp_normal_form.h" // perhaps make c++ headers instead?
-}
-
-#endif
-
// Abstract class representing a cell of a cell complex.
class Cell
{
@@ -340,11 +331,7 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
- // boundary operator matrices for this chain complex
- // h_k: C_k -> C_(k-1)
- //gmp_matrix* _HMatrix[4];
-
- public:
+ public:
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
@@ -380,6 +367,8 @@ class CellComplex
// 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]; }
+
// 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(); }
@@ -421,7 +410,8 @@ class CellComplex
// useful functions for (co)reduction of cell complex
virtual void reduceComplex();
- virtual void coreduceComplex();
+ // 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);
@@ -437,97 +427,7 @@ class CellComplex
// write this cell complex in legacy .msh format
virtual int writeComplexMSH(const std::string &name);
- // construct boundary operator matrix dim->dim-1
- //virtual void constructHMatrix(int dim);
- // 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 representations and great computational complexity in its methods.
-class ChainComplex{
- private:
- // boundary operator matrices for this chain complex
- // h_k: C_k -> C_(k-1)
- gmp_matrix* _HMatrix[4];
-
- // Basis matrices for the kernels and codomains of the boundary operator matrices
- gmp_matrix* _kerH[4];
- gmp_matrix* _codH[4];
-
- gmp_matrix* _JMatrix[4];
- gmp_matrix* _QMatrix[4];
- std::vector<long int> _torsion[4];
-
- // bases for homology groups
- gmp_matrix* _Hbasis[4];
-
- public:
- ChainComplex( std::vector<gmp_matrix*> HMatrix ){
- for(int i = 0; i < HMatrix.size(); i++){
- _HMatrix[i] = HMatrix.at(i);
- }
- for(int i = 0; i < 4; i++){
- _kerH[i] = NULL;
- _codH[i] = NULL;
- _JMatrix[i] = NULL;
- _QMatrix[i] = NULL;
- _Hbasis[i] = NULL;
- }
-
- }
- ChainComplex(){
- for(int i = 0; i < 4; i++){
- _HMatrix[i] = create_gmp_matrix_zero(1,1);
- _kerH[i] = NULL;
- _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) { 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);
- // Compute quotient problem for given inclusion relation j to find representatives of homology groups
- // and possible torsion coeffcients
- 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));
- return gmp_matrix_printf(matrix); }
-
- // debugging aid
- virtual void matrixTest();
};
-
-#endif
-
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5bdcd8f5170c35654bbbdf9f15146b66f3731dae
--- /dev/null
+++ b/Geo/ChainComplex.cpp
@@ -0,0 +1,307 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka.
+
+#include "ChainComplex.h"
+
+#if defined(HAVE_KBIPACK)
+
+ChainComplex::ChainComplex(CellComplex* cellComplex){
+
+ _HMatrix[0] = NULL;
+
+ for(int dim = 1; dim < 4; dim++){
+ unsigned int cols = cellComplex->getSize(dim);
+ unsigned int rows = cellComplex->getSize(dim-1);
+
+ for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()) cols--;
+ }
+ 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--;
+ }
+
+
+ if( cols == 0 ){
+ _HMatrix[dim] = NULL;
+ }
+ else if( rows == 0){
+ _HMatrix[dim] = create_gmp_matrix_zero(1, cols);
+ }
+ else{
+
+ long int elems[rows*cols];
+
+ std::set<Cell*, Less_Cell>::iterator high = cellComplex->firstCell(dim);
+ std::set<Cell*, Less_Cell>::iterator low = cellComplex->firstCell(dim-1);
+
+ unsigned int i = 0;
+ while(i < rows*cols){
+ while(low != cellComplex->lastCell(dim-1)){
+ Cell* lowcell = *low;
+ Cell* highcell = *high;
+ if(!(highcell->inSubdomain() || lowcell->inSubdomain())){
+ elems[i] = cellComplex->kappa(*high, *low);
+ i++;
+ }
+ low++;
+ }
+ low = cellComplex->firstCell(dim-1);
+ high++;
+ }
+ _HMatrix[dim] = create_gmp_matrix_int(rows, cols, elems);
+ }
+
+ _kerH[dim] = NULL;
+ _codH[dim] = NULL;
+ _JMatrix[dim] = NULL;
+ _QMatrix[dim] = NULL;
+ _Hbasis[dim] = NULL;
+ }
+ return;
+}
+
+
+void ChainComplex::KerCod(int dim){
+
+ if(dim < 1 || dim > 3 || _HMatrix[dim] == NULL) return;
+
+ gmp_matrix* HMatrix = copy_gmp_matrix(_HMatrix[dim], 1, 1,
+ gmp_matrix_rows(_HMatrix[dim]), gmp_matrix_cols(_HMatrix[dim]));
+
+ gmp_normal_form* normalForm = create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED);
+ //printMatrix(normalForm->left);
+ //printMatrix(normalForm->canonical);
+ //printMatrix(normalForm->right);
+
+ int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical), gmp_matrix_cols(normalForm->canonical));
+ int rank = 0;
+ mpz_t elem;
+ mpz_init(elem);
+
+ while(rank < minRowCol){
+ gmp_matrix_get_elem(elem, rank+1, rank+1, normalForm->canonical);
+ if(mpz_cmp_si(elem,0) == 0) break;
+ rank++;
+ }
+
+ if(rank != gmp_matrix_cols(normalForm->canonical)){
+ _kerH[dim] = copy_gmp_matrix(normalForm->right, 1, rank+1,
+ gmp_matrix_rows(normalForm->right), gmp_matrix_cols(normalForm->right));
+ }
+
+ if(rank > 0){
+ _codH[dim] = copy_gmp_matrix(normalForm->canonical, 1, 1,
+ gmp_matrix_rows(normalForm->canonical), rank);
+ gmp_matrix_left_mult(normalForm->left, _codH[dim]);
+ }
+
+ mpz_clear(elem);
+ destroy_gmp_normal_form(normalForm);
+
+ return;
+}
+
+//j:B_(k+1)->Z_k
+void ChainComplex::Inclusion(int dim){
+
+ 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);
+ if(rows < cols) return;
+
+ rows = gmp_matrix_rows(Bbasis);
+ cols = gmp_matrix_cols(Bbasis);
+ if(rows < cols) return;
+
+ // A*inv(V) = U*S
+ 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++){
+ gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
+ if(mpz_cmp_si(elem,0) == 0){
+ destroy_gmp_normal_form(normalForm);
+ return;
+ }
+ }
+
+ 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] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(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(){
+
+ int rows = 3;
+ int cols = 6;
+
+ long int elems[rows*cols];
+ for(int i = 1; i<=rows*cols; i++) elems[i-1] = i;
+
+ gmp_matrix* matrix = create_gmp_matrix_int(rows, cols, elems);
+
+ gmp_matrix* copymatrix = copy_gmp_matrix(matrix, 3, 2, 3, 5);
+
+ printMatrix(matrix);
+ printMatrix(copymatrix);
+
+ return;
+}
+
+
+#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
new file mode 100644
index 0000000000000000000000000000000000000000..89b505a900b7722385bdcd131a310c47dd597f34
--- /dev/null
+++ b/Geo/ChainComplex.h
@@ -0,0 +1,99 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka.
+
+#ifndef _CHAINCOMPLEX_H_
+#define _CHAINCOMPLEX_H_
+
+#include <stdio.h>
+#include <string>
+#include <algorithm>
+#include <set>
+#include <queue>
+#include "GmshConfig.h"
+#include "MElement.h"
+#include "GModel.h"
+#include "GEntity.h"
+#include "GRegion.h"
+#include "GFace.h"
+#include "GVertex.h"
+#include "CellComplex.h"
+
+#if defined(HAVE_KBIPACK)
+
+#include "gmp.h"
+extern "C" {
+ #include "gmp_normal_form.h" // perhaps make c++ headers instead?
+}
+
+// A class representing a chain complex of a cell complex.
+// This should only be constructed for a reduced cell complex because of
+// dense matrix representations and great computational complexity in its methods.
+class ChainComplex{
+ private:
+ // boundary operator matrices for this chain complex
+ // h_k: C_k -> C_(k-1)
+ gmp_matrix* _HMatrix[4];
+
+ // Basis matrices for the kernels and codomains of the boundary operator matrices
+ gmp_matrix* _kerH[4];
+ gmp_matrix* _codH[4];
+
+ gmp_matrix* _JMatrix[4];
+ gmp_matrix* _QMatrix[4];
+ std::vector<long int> _torsion[4];
+
+ // bases for homology groups
+ gmp_matrix* _Hbasis[4];
+
+ public:
+
+ ChainComplex(CellComplex* cellComplex);
+
+ ChainComplex(){
+ for(int i = 0; i < 4; i++){
+ _HMatrix[i] = create_gmp_matrix_zero(1,1);
+ _kerH[i] = NULL;
+ _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) { 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);
+ // Compute quotient problem for given inclusion relation j to find representatives of homology groups
+ // and possible torsion coeffcients
+ 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));
+ return gmp_matrix_printf(matrix); }
+
+ // debugging aid
+ virtual void matrixTest();
+};
+
+
+#endif
+
+#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index 33438837597c590354d4661b209464cd5c8947b9..fcb2a22ac12dbb1f4c49ebae596c16a6e712d9ac 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -39,7 +39,8 @@ SRC = GEntity.cpp\
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp\
MHexahedron.cpp MPrism.cpp MPyramid.cpp\
MZone.cpp MZoneBoundary.cpp\
- CellComplex.cpp
+ CellComplex.cpp\
+ ChainComplex.cpp
OBJ = ${SRC:.cpp=${OBJEXT}}
@@ -387,3 +388,4 @@ CellComplex${OBJEXT}: CellComplex.cpp CellComplex.h ../Common/GmshConfig.h \
../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
GModel.h GVertex.h GEntity.h Range.h SBoundingBox3d.h GPoint.h GEdge.h \
GFace.h GEdgeLoop.h Pair.h GRegion.h
+ChainComplex${OBJEXT}: ChainComplex.cpp ChainComplex.h
diff --git a/Makefile b/Makefile
index 36c7086ca76157afb02554fd1617a1bc8abb62e5..1be159addb4e14ba67d8659fe59fcae9744e55c2 100644
--- a/Makefile
+++ b/Makefile
@@ -27,6 +27,7 @@ GMSH_API =\
Geo/discreteVertex.h Geo/discreteEdge.h Geo/discreteFace.h Geo/discreteRegion.h\
Geo/SPoint2.h Geo/SPoint3.h Geo/SVector3.h Geo/SBoundingBox3d.h Geo/Pair.h Geo/Range.h\
Geo/CellComplex.h\
+ Geo/ChainComplex.h\
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h contrib/kbipack/gmp_blas.h\
Numeric/Gauss.h Numeric/FunctionSpace.h Numeric/GmshMatrix.h\
Numeric/gmshAssembler.h Numeric/gmshTermOfFormulation.h Numeric/gmshLaplace.h\
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index c6dcb11fac8ecec18015404d31fddd310f20c005..47bf32dd827b313831a69612947b097d5479e49b 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.3 2009-04-14 10:02:22 matti Exp $
+ $Id: gmp_matrix.c,v 1.4 2009-04-21 07:06:22 matti Exp $
*/
@@ -182,6 +182,8 @@ copy_gmp_matrix(const gmp_matrix * matrix,
r = end_row-start_row+1;
c = end_col-start_col+1;
+ if(r < 1 || c < 1) return NULL;
+
new_matrix -> storage = (mpz_t *) calloc(r*c, sizeof(mpz_t));
if(new_matrix -> storage == NULL)
{