Skip to content
Snippets Groups Projects
Commit 69a8c6ad authored by Matti Pellika's avatar Matti Pellika
Browse files

Added files ChainComplex.cpp/.h, bugfixes.

parent 33792858
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/CellComplex.cpp
+ 4
− 298
View file @ 69a8c6ad
......@@ -403,7 +403,7 @@ 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];
......@@ -423,6 +423,7 @@ void CellComplex::coreduceComplex(){
coreduction(1);
coreduction(2);
}
if(generatorDim == i) break;
}
for(int i = 0; i < 4; i++){
......@@ -437,65 +438,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){
std::vector<GEntity*> domain;
......@@ -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
Geo/CellComplex.h
+ 5
− 105
View file @ 69a8c6ad
......@@ -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,10 +331,6 @@ 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:
// 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
Geo/ChainComplex.cpp 0 → 100644
+ 307
− 0
View file @ 69a8c6ad
// 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
Geo/ChainComplex.h 0 → 100644
+ 99
− 0
View file @ 69a8c6ad
// 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
Geo/Makefile
+ 3
− 1
View file @ 69a8c6ad
......@@ -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
Makefile
+ 1
− 0
View file @ 69a8c6ad
......@@ -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\
......
contrib/kbipack/gmp_matrix.c
+ 3
− 1
View file @ 69a8c6ad
......@@ -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)
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment