diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index f59cc141332e9154eba9df04d172ffa75fbd0eba..c19622eb2463c7641e2ba1f109e612eaac27e6df 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -57,8 +57,8 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
if(type == TYPE_TET) newtype = MSH_TRI_3;
else if(type == TYPE_HEX) newtype = MSH_QUA_4;
else if(type == TYPE_PRI) {
- if(vertices.size() == 3) newtype = MSH_TRI_3;
- else if(vertices.size() == 4) newtype = MSH_QUA_4;
+ if(vertices.size() == 3) newtype = MSH_TRI_3;
+ else if(vertices.size() == 4) newtype = MSH_QUA_4;
}
}
else if(_dim == 2) newtype = MSH_LIN_2;
@@ -68,7 +68,7 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
return false;
}
MElement* element = factory.create(newtype, vertices, 0,
- _image->getPartition());
+ _image->getPartition());
Cell* cell = new Cell(element);
bdCells.push_back(cell);
}
@@ -101,7 +101,7 @@ void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const
return;
}
-int Cell::getFacetOri(Cell* cell)
+int Cell::findBoundaryCellOrientation(Cell* cell)
{
if(_image == NULL){
printf("ERROR: No image mesh element for cell. \n");
@@ -141,7 +141,7 @@ int Cell::getFacetOri(Cell* cell)
bool Cell::hasVertex(int vertex) const
{
std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(),
- vertex);
+ vertex);
if (it != _vs.end()) return true;
else return false;
}
@@ -166,7 +166,7 @@ void Cell::restoreCell(){
}
void Cell::addBoundaryCell(int orientation, Cell* cell,
- bool orig, bool other)
+ bool orig, bool other)
{
if(orig) _obd.insert( std::make_pair(cell, orientation ) );
biter it = _bd.find(cell);
@@ -184,7 +184,7 @@ void Cell::addBoundaryCell(int orientation, Cell* cell,
}
void Cell::addCoboundaryCell(int orientation, Cell* cell,
- bool orig, bool other)
+ bool orig, bool other)
{
if(orig) _ocbd.insert( std::make_pair(cell, orientation ) );
biter it = _cbd.find(cell);
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 46b516af3b0cafbcdbe3ab70880692b691ea9fd2..0d0dd3e6caf3312c295e3852176ea06a334d8a9f 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -59,15 +59,17 @@ class Cell
std::map<Cell*, int, Less_Cell > _obd;
std::map<Cell*, int, Less_Cell > _ocbd;
+ // sorted vertices of this cell (used for ordering of the cells)
+ std::vector<int> _vs;
+
+
+
// the mesh element that is the image of this cell
MElement* _image;
// whether to delete the mesh element when done
- // (created for homology computation only)
+ // (set to true if created for homology computation only)
bool _delimage;
- // sorted vertices of this cell (used for ordering of the cells)
- std::vector<int> _vs;
-
// get vertex
MVertex* getVertex(int vertex) const {
if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
@@ -99,10 +101,11 @@ class Cell
void setDeleteImage(bool delimage) { _delimage = delimage; }
bool getDeleteImage() const { return _delimage; }
- // find the cells on the boundary of this cell
+ // find the cells on the boundary of this cell
bool findBoundaryCells(std::vector<Cell*>& bdCells);
// get boundary cell orientation
- int getFacetOri(Cell* cell);
+ int findBoundaryCellOrientation(Cell* cell);
+
int getDim() const { return _dim; };
int getIndex() const { return _index; };
@@ -123,6 +126,8 @@ class Cell
// (co)boundary cell iterator
typedef std::map<Cell*, int, Less_Cell>::iterator biter;
+ // iterators to (first/last (co)boundary cells of this cell
+ // (orig: to original (co)boundary cells of this cell)
biter firstBoundary(bool orig=false){
return orig ? _obd.begin() : _bd.begin(); }
biter lastBoundary(bool orig=false){
@@ -135,7 +140,7 @@ class Cell
int getBoundarySize() { return _bd.size(); }
int getCoboundarySize() { return _cbd.size(); }
- // get the cell boundary
+ // get the (orig: original) cell boundary
void getBoundary(std::map<Cell*, int, Less_Cell >& boundary,
bool orig=false){
orig ? boundary = _obd : boundary = _bd; }
@@ -143,17 +148,19 @@ class Cell
bool orig=false){
orig ? coboundary = _ocbd : coboundary = _cbd; }
- // add (co)boundary cell
+ // add (co)boundary cell (orig: as original)
+ // (other: reciprocally also add this cell from the other cell's (co)boundary)
void addBoundaryCell(int orientation, Cell* cell,
bool orig, bool other);
void addCoboundaryCell(int orientation, Cell* cell,
bool orig, bool other);
- // remove (co)boundary cell
+ // remove (co)boundary cell
+ // (other: reciprocally also revove this cell from the other cell's (co)boundary)
void removeBoundaryCell(Cell* cell, bool other);
void removeCoboundaryCell(Cell* cell, bool other);
- // true if has given cell on (original) (co)boundary
+ // true if has given cell on (orig: original) (co)boundary
bool hasBoundary(Cell* cell, bool orig=false);
bool hasCoboundary(Cell* cell, bool orig=false);
@@ -184,7 +191,7 @@ class Cell
}
for(int i=0; i < this->getNumSortedVertices();i++){
if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
- return false;
+ return false;
}
}
return true;
@@ -212,6 +219,7 @@ class CombinedCell : public Cell{
int getNum() const { return _num; }
bool operator==(const Cell& c2) const {
+ if(this->isCombined() != c2.isCombined()) return false;
if(this->getNum() != c2.getNum()) return false;
else return true;
}
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 4a87e357193fd655be500ee0ff0bc6cf78cbac6e..26162fee1f55ddda8c94d3d70f90a764fc4ad5e0 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -10,7 +10,7 @@
#if defined(HAVE_KBIPACK)
CellComplex::CellComplex( std::vector<MElement*>& domainElements,
- std::vector<MElement*>& subdomainElements)
+ std::vector<MElement*>& subdomainElements)
{
_dim = 0;
_simplicial = true;
@@ -41,7 +41,7 @@ void CellComplex::panic_exit(){
}
bool CellComplex::insert_cells(std::vector<MElement*>& elements,
- bool subdomain)
+ bool subdomain)
{
// add highest dimensional cells
for(unsigned int i=0; i < elements.size(); i++){
@@ -49,7 +49,7 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
int type = element->getType();
if(_simplicial && !(type == TYPE_PNT || type == TYPE_LIN
- || type == TYPE_TRI || type == TYPE_TET) ){
+ || type == TYPE_TRI || type == TYPE_TET) ){
_simplicial = false;
}
//FIXME: no getFaceInfo methods for these MElements
@@ -70,19 +70,19 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
std::vector<Cell*> bdCells;
if(!cell->findBoundaryCells(bdCells)) return false;
for(unsigned int i = 0; i < bdCells.size(); i++){
- Cell* newCell = bdCells.at(i);
- newCell->setInSubdomain(subdomain);
- newCell->setDeleteImage(true);
- std::pair<citer, bool> insertInfo =
- _cells[newCell->getDim()].insert(newCell);
- if(!insertInfo.second){ // the cell was already in the cell complex
- delete newCell;
- newCell = *(insertInfo.first);
- }
- if(!subdomain) {
- int ori = cell->getFacetOri(newCell);
- cell->addBoundaryCell( ori, newCell, true, true);
- }
+ Cell* newCell = bdCells.at(i);
+ newCell->setInSubdomain(subdomain);
+ newCell->setDeleteImage(true);
+ std::pair<citer, bool> insertInfo =
+ _cells[newCell->getDim()].insert(newCell);
+ if(!insertInfo.second){ // the cell was already in the cell complex
+ delete newCell;
+ newCell = *(insertInfo.first);
+ }
+ if(!subdomain) {
+ int ori = cell->findBoundaryCellOrientation(newCell);
+ cell->addBoundaryCell( ori, newCell, true, true);
+ }
}
}
}
@@ -138,14 +138,14 @@ void CellComplex::removeCell(Cell* cell, bool other)
}
void CellComplex::removeCellQset(Cell* cell,
- std::set<Cell*, Less_Cell>& Qset)
+ std::set<Cell*, Less_Cell>& Qset)
{
Qset.erase(cell);
}
void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
- std::queue<Cell*>& Q,
- std::set<Cell*, Less_Cell>& Qset)
+ std::queue<Cell*>& Q,
+ std::set<Cell*, Less_Cell>& Qset)
{
for(std::map<Cell*, int, Less_Cell>::iterator cit = cells.begin();
cit != cells.end(); cit++){
@@ -185,7 +185,7 @@ int CellComplex::coreduction(Cell* startCell, int omitted)
enqueueCells(cbd_c, Q, Qset);
removeCell(bd_s.begin()->first);
if(bd_s.begin()->first->getDim() == 0 && omitted > 0){
- _store.at(omitted-1).insert(bd_s.begin()->first);
+ _store.at(omitted-1).insert(bd_s.begin()->first);
}
coreductions++;
}
@@ -212,15 +212,15 @@ int CellComplex::reduction(int dim, int omitted)
while(cit != lastCell(dim-1)){
Cell* cell = *cit;
if( cell->getCoboundarySize() == 1
- && inSameDomain(cell, cell->firstCoboundary()->first)){
- ++cit;
- if(dim == getDim() && omitted > 0){
- _store.at(omitted-1).insert(cell->firstCoboundary()->first);
- }
- removeCell(cell->firstCoboundary()->first);
- removeCell(cell);
- count++;
- reduced = true;
+ && inSameDomain(cell, cell->firstCoboundary()->first)){
+ ++cit;
+ if(dim == getDim() && omitted > 0){
+ _store.at(omitted-1).insert(cell->firstCoboundary()->first);
+ }
+ removeCell(cell->firstCoboundary()->first);
+ removeCell(cell);
+ count++;
+ reduced = true;
}
if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
@@ -247,9 +247,9 @@ int CellComplex::coreduction(int dim, int omitted)
if( cell->getBoundarySize() == 1
&& inSameDomain(cell, cell->firstBoundary()->first)){
++cit;
- if(dim-1 == 0 && omitted > 0){
- _store.at(omitted-1).insert(cell->firstBoundary()->first);
- }
+ if(dim-1 == 0 && omitted > 0){
+ _store.at(omitted-1).insert(cell->firstBoundary()->first);
+ }
removeCell(cell->firstBoundary()->first);
removeCell(cell);
count++;
@@ -266,7 +266,7 @@ int CellComplex::coreduction(int dim, int omitted)
int CellComplex::reduceComplex(bool omit)
{
printf("Cell Complex: \n %d volumes, %d faces, %d edges and %d vertices. \n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
@@ -301,7 +301,7 @@ int CellComplex::reduceComplex(bool omit)
combine(1);
printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
return 0;
}
@@ -321,7 +321,7 @@ void CellComplex::removeSubdomain()
int CellComplex::coreduceComplex(bool omit)
{
printf("Cell Complex: \n %d volumes, %d faces, %d edges and %d vertices. \n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
removeSubdomain();
@@ -354,7 +354,7 @@ int CellComplex::coreduceComplex(bool omit)
}
printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
cocombine(0);
coreduction(1);
@@ -365,7 +365,7 @@ int CellComplex::coreduceComplex(bool omit)
coherent();
printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ getSize(3), getSize(2), getSize(1), getSize(0));
return omitted;
}
@@ -390,7 +390,7 @@ int CellComplex::cocombine(int dim)
Cell* s = Q.front();
Q.pop();
if(s->getBoundarySize() == 2){
- Cell::biter it = s->firstBoundary();
+ Cell::biter it = s->firstBoundary();
int or1 = (*it).second;
Cell* c1 = (*it).first;
it++;
@@ -399,7 +399,7 @@ int CellComplex::cocombine(int dim)
if(!(*c1 == *c2) && abs(or1) == abs(or2)
&& inSameDomain(s, c1) && inSameDomain(s, c2)){
- removeCell(s);
+ removeCell(s);
c1->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
@@ -407,7 +407,7 @@ int CellComplex::cocombine(int dim)
enqueueCells(cbd_c, Q, Qset);
CombinedCell* newCell = new CombinedCell(c1, c2,
- (or1 != or2), true );
+ (or1 != or2), true );
removeCell(c1);
removeCell(c2);
insertCell(newCell);
@@ -448,7 +448,7 @@ int CellComplex::combine(int dim)
Q.pop();
if(s->getCoboundarySize() == 2){
- Cell::biter it = s->firstCoboundary();
+ Cell::biter it = s->firstCoboundary();
int or1 = (*it).second;
Cell* c1 = (*it).first;
it++;
@@ -468,7 +468,7 @@ int CellComplex::combine(int dim)
removeCell(c1);
removeCell(c2);
insertCell(newCell);
-
+
cit = firstCell(dim);
count++;
}
@@ -492,7 +492,7 @@ bool CellComplex::coherent()
std::map<Cell*, int, Less_Cell> boundary;
cell->getBoundary(boundary);
for(Cell::biter it = boundary.begin();
- it != boundary.end(); it++){
+ it != boundary.end(); it++){
Cell* bdCell = (*it).first;
int ori = (*it).second;
citer cit = _cells[bdCell->getDim()].find(bdCell);
@@ -503,7 +503,7 @@ bool CellComplex::coherent()
}
if(!bdCell->hasCoboundary(cell)){
printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
- bdCell->addCoboundaryCell(ori, cell, false, false);
+ bdCell->addCoboundaryCell(ori, cell, false, false);
coherent = false;
}
@@ -511,7 +511,7 @@ bool CellComplex::coherent()
std::map<Cell*, int, Less_Cell> coboundary;
cell->getCoboundary(coboundary);
for(Cell::biter it = coboundary.begin();
- it != coboundary.end(); it++){
+ it != coboundary.end(); it++){
Cell* cbdCell = (*it).first;
int ori = (*it).second;
citer cit = _cells[cbdCell->getDim()].find(cbdCell);
@@ -522,7 +522,7 @@ bool CellComplex::coherent()
}
if(!cbdCell->hasBoundary(cell)){
printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
- cbdCell->addBoundaryCell(ori, cell, false, false);
+ cbdCell->addBoundaryCell(ori, cell, false, false);
coherent = false;
}
@@ -543,12 +543,12 @@ bool CellComplex::hasCell(Cell* cell, bool orig)
}
void CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
- int dim, int domain){
+ int dim, int domain){
cells.clear();
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
if( (domain == 0 && !cell->inSubdomain()) || domain == 1
- || (domain == 2 && cell->inSubdomain()) ){
+ || (domain == 2 && cell->inSubdomain()) ){
cells.insert(cell);
}
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index ae6ced80e378880da21b28b0da12ecaec9a2021e..a801d1f9281fa219dcf7892442a169bf76e2ccb5 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -48,7 +48,7 @@ class CellComplex
// enqueue cells in queue if they are not there already
void enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
- std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
+ std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
// remove cell from the queue set
void removeCellQset(Cell* cell, std::set<Cell*, Less_Cell>& Qset);
@@ -66,7 +66,7 @@ class CellComplex
public:
CellComplex( std::vector<MElement*>& domainElements,
- std::vector<MElement*>& subdomainElements);
+ std::vector<MElement*>& subdomainElements);
~CellComplex();
// get the number of certain dimensional cells
@@ -99,7 +99,7 @@ class CellComplex
// check whether two cells both belong to subdomain or if neither one does
bool inSameDomain(Cell* c1, Cell* c2) const {
return ( (!c1->inSubdomain() && !c2->inSubdomain())
- || (c1->inSubdomain() && c2->inSubdomain() )); }
+ || (c1->inSubdomain() && c2->inSubdomain() )); }
// remove cells in subdomain from this cell complex
void removeSubdomain();
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 3bffbd8d68504da11a652c86d860641e1dd7e314..2a0c624583cf60e26c53da26973fbd1ab6a4fab9 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -15,55 +15,43 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
_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 i = 0; i < 5; i++){
+ _HMatrix[i] = NULL;
+ _kerH[i] = NULL;
+ _codH[i] = NULL;
+ _JMatrix[i] = NULL;
+ _QMatrix[i] = NULL;
+ _Hbasis[i] = NULL;
+ }
+
+ int lastCols = 0;
for(int dim = 0; dim < 4; dim++){
unsigned int cols = cellComplex->getSize(dim);
unsigned int rows = 0;
- if(dim > 0) rows = cellComplex->getSize(dim-1);
int index = 1;
// ignore cells depending on domain
for(CellComplex::citer cit = cellComplex->firstCell(dim);
- cit != cellComplex->lastCell(dim); cit++){
+ cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
cell->setIndex(0);
cols--;
- if((domain == 0 && !cell->inSubdomain())
- || (domain == 2 && cell->inSubdomain()) ){
+ if((domain == 0 && !cell->inSubdomain()) || domain == 1
+ || (domain == 2 && cell->inSubdomain()) ){
cols++;
- cell->setIndex(index);
- index++;
- _cellIndices[dim][cell->getIndex()] = cell;
- }
- }
- index = 1;
- if(dim > 0){
- for(CellComplex::citer cit = cellComplex->firstCell(dim-1);
- cit != cellComplex->lastCell(dim-1); cit++){
- Cell* cell = *cit;
- cell->setIndex(0);
- rows--;
- if((domain == 0 && !cell->inSubdomain())
- || (domain == 2 && cell->inSubdomain()) ){
- rows++;
- cell->setIndex(index);
- index++;
- _cellIndices[dim-1][cell->getIndex()] = cell;
- }
+ cell->setIndex(index);
+ index++;
+ _cellIndices[dim][cell] = cell->getIndex();
}
}
+ if(dim > 0) rows = lastCols;
+ lastCols = cols;
- if( cols == 0 ){ // no dim-cells, no map
+ if(cols == 0){ // no dim-cells, no map
//_HMatrix[dim] = create_gmp_matrix_zero(rows, 1);
_HMatrix[dim] = NULL;
}
- else if( rows == 0){ // no dim-1-cells, maps everything to zero
+ else if(rows == 0){ // no dim-1-cells, maps everything to zero
_HMatrix[dim] = create_gmp_matrix_zero(1, cols);
//_HMatrix[dim] = NULL;
}
@@ -73,34 +61,34 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
mpz_init(elem);
_HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
for( std::set<Cell*, Less_Cell>::iterator cit =
- cellComplex->firstCell(dim);
- cit != cellComplex->lastCell(dim); cit++){
+ cellComplex->firstCell(dim);
+ cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if( (domain == 0 && !cell->inSubdomain()) || domain == 1
- || (domain == 2 && cell->inSubdomain()) ){
+ || (domain == 2 && cell->inSubdomain()) ){
for(Cell::biter it = cell->firstBoundary();
- it != cell->lastBoundary(); it++){
+ it != cell->lastBoundary(); it++){
Cell* bdCell = (*it).first;
if((domain == 0 && !bdCell->inSubdomain()) || domain == 1
- || (domain == 2 && cell->inSubdomain()) ){
+ || (domain == 2 && cell->inSubdomain()) ){
int old_elem = 0;
if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim])
- || bdCell->getIndex() < 1
+ || bdCell->getIndex() < 1
|| cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim])
- || cell->getIndex() < 1){
- //printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
+ || cell->getIndex() < 1){
+ printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
}
else{
gmp_matrix_get_elem(elem, bdCell->getIndex(),
- cell->getIndex(), _HMatrix[dim]);
+ cell->getIndex(), _HMatrix[dim]);
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).second);
if( abs((old_elem + (*it).second)) > 1){
- //printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
+ //printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(),
- cell->getIndex(), _HMatrix[dim]);
+ cell->getIndex(), _HMatrix[dim]);
}
}
}
@@ -115,9 +103,19 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
_QMatrix[dim] = NULL;
_Hbasis[dim] = NULL;
}
- return;
}
+ChainComplex::~ChainComplex()
+{
+ for(int i = 0; i < 5; i++){
+ destroy_gmp_matrix(_HMatrix[i]);
+ destroy_gmp_matrix(_kerH[i]);
+ destroy_gmp_matrix(_codH[i]);
+ destroy_gmp_matrix(_JMatrix[i]);
+ destroy_gmp_matrix(_QMatrix[i]);
+ destroy_gmp_matrix(_Hbasis[i]);
+ }
+}
void ChainComplex::KerCod(int dim)
{
@@ -125,8 +123,8 @@ void ChainComplex::KerCod(int dim)
gmp_matrix* HMatrix
= copy_gmp_matrix(_HMatrix[dim], 1, 1,
- gmp_matrix_rows(_HMatrix[dim]),
- gmp_matrix_cols(_HMatrix[dim]));
+ gmp_matrix_rows(_HMatrix[dim]),
+ gmp_matrix_cols(_HMatrix[dim]));
gmp_normal_form* normalForm
= create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED);
@@ -135,7 +133,7 @@ void ChainComplex::KerCod(int dim)
//printMatrix(normalForm->right);
int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical),
- gmp_matrix_cols(normalForm->canonical));
+ gmp_matrix_cols(normalForm->canonical));
int rank = 0;
mpz_t elem;
mpz_init(elem);
@@ -150,14 +148,14 @@ void ChainComplex::KerCod(int dim)
if(rank != (int)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));
+ 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_rows(normalForm->canonical), rank);
gmp_matrix_left_mult(normalForm->left, _codH[dim]);
}
@@ -167,7 +165,7 @@ void ChainComplex::KerCod(int dim)
return;
}
-//j:B_(k+1)->Z_k
+//j:B_k->Z_k
void ChainComplex::Inclusion(int lowDim, int highDim)
{
if(getKerHMatrix(lowDim) == NULL
@@ -176,12 +174,12 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
gmp_matrix* Zbasis =
copy_gmp_matrix(_kerH[lowDim], 1, 1,
- gmp_matrix_rows(_kerH[lowDim]),
- gmp_matrix_cols(_kerH[lowDim]));
+ gmp_matrix_rows(_kerH[lowDim]),
+ gmp_matrix_cols(_kerH[lowDim]));
gmp_matrix* Bbasis
= copy_gmp_matrix(_codH[highDim], 1, 1,
- gmp_matrix_rows(_codH[highDim]),
- gmp_matrix_cols(_codH[highDim]));
+ gmp_matrix_rows(_codH[highDim]),
+ gmp_matrix_cols(_codH[highDim]));
int rows = gmp_matrix_rows(Bbasis);
@@ -192,7 +190,7 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
cols = gmp_matrix_cols(Zbasis);
if(rows < cols) return;
- // A*inv(V) = U*S
+ // inv(U)*A*inv(V) = S
gmp_normal_form* normalForm
= create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED);
@@ -211,8 +209,8 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
gmp_matrix_left_mult(normalForm->left, Bbasis);
gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1,
- gmp_matrix_cols(Zbasis),
- gmp_matrix_cols(Bbasis));
+ gmp_matrix_cols(Zbasis),
+ gmp_matrix_cols(Bbasis));
destroy_gmp_matrix(Bbasis);
rows = gmp_matrix_rows(LB);
@@ -245,8 +243,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
mpz_clear(divisor);
mpz_clear(result);
destroy_gmp_normal_form(normalForm);
-
- return;
}
void ChainComplex::Quotient(int dim)
@@ -255,8 +251,8 @@ void ChainComplex::Quotient(int dim)
gmp_matrix* JMatrix =
copy_gmp_matrix(_JMatrix[dim], 1, 1,
- gmp_matrix_rows(_JMatrix[dim]),
- gmp_matrix_cols(_JMatrix[dim]));
+ gmp_matrix_rows(_JMatrix[dim]),
+ gmp_matrix_cols(_JMatrix[dim]));
int rows = gmp_matrix_rows(JMatrix);
int cols = gmp_matrix_cols(JMatrix);
@@ -321,13 +317,13 @@ void ChainComplex::computeHomology(bool dual)
&& gmp_matrix_cols(getHMatrix(lowDim)) > 0
&& getHMatrix(highDim) == NULL) {
setHbasis( setDim,
- create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) );
+ create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) );
}
else if(highDim == 0 && dual
- && gmp_matrix_rows(getHMatrix(highDim)) > 0
- && getHMatrix(lowDim) == NULL) {
+ && gmp_matrix_rows(getHMatrix(highDim)) > 0
+ && getHMatrix(lowDim) == NULL) {
setHbasis( setDim,
- create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
+ create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
}
// 2) this dimension is empty
@@ -337,9 +333,9 @@ void ChainComplex::computeHomology(bool dual)
// 3) No higher dimension cells -> none of the cycles are boundaries
else if(getHMatrix(highDim) == NULL){
setHbasis( setDim,
- copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
- gmp_matrix_rows(getKerHMatrix(lowDim)),
- gmp_matrix_cols(getKerHMatrix(lowDim))) );
+ copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(lowDim)),
+ gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
@@ -352,37 +348,37 @@ void ChainComplex::computeHomology(bool dual)
// 4) No lower dimension cells -> all chains are cycles
if(getHMatrix(lowDim) == NULL){
setKerHMatrix(lowDim,
- create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
+ create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
}
Inclusion(lowDim, highDim);
Quotient(lowDim);
if(getCodHMatrix(highDim) == NULL){
setHbasis(setDim,
- copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
- gmp_matrix_rows(getKerHMatrix(lowDim)),
- gmp_matrix_cols(getKerHMatrix(lowDim))) );
+ 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(setDim, NULL);
}
else{
setHbasis(setDim,
- copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
- gmp_matrix_rows(getKerHMatrix(lowDim)),
- gmp_matrix_cols(getKerHMatrix(lowDim))) );
+ copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(lowDim)),
+ gmp_matrix_cols(getKerHMatrix(lowDim))) );
gmp_matrix_right_mult(getHbasis(setDim), getQMatrix(lowDim));
}
}
-
- destroy_gmp_matrix(getKerHMatrix(lowDim));
- destroy_gmp_matrix(getCodHMatrix(lowDim));
+
+ //destroy_gmp_matrix(getKerHMatrix(lowDim));
+ //destroy_gmp_matrix(getCodHMatrix(lowDim));
destroy_gmp_matrix(getJMatrix(lowDim));
destroy_gmp_matrix(getQMatrix(lowDim));
- setKerHMatrix(lowDim, NULL);
- setCodHMatrix(lowDim, NULL);
+ //setKerHMatrix(lowDim, NULL);
+ //setCodHMatrix(lowDim, NULL);
setJMatrix(lowDim, NULL);
setQMatrix(lowDim, NULL);
}
@@ -433,50 +429,406 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber)
return coeffVector;
}
-void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain,
- int dim, int chainNumber)
+gmp_matrix* ChainComplex::getBasis(int dim, int basis)
{
+ if(dim > -2 && dim < 5 && basis == 2) return _codH[dim+1];
+ if(dim < 0 || dim > 4) return NULL;
+ if(basis == 0) return create_gmp_matrix_identity(getBasisSize(dim, 0));
+ else if(basis == 1) return _kerH[dim];
+ else if(basis == 3) return _Hbasis[dim];
+ else return NULL;
+}
+
+void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
+ int num, int dim, int basis)
+{
+ gmp_matrix* basisMatrix;
+ if(basis == 0) basisMatrix = getBasis(dim, 0);
+ else if(basis == 1) basisMatrix = getBasis(dim, 1);
+ else if(basis == 2) basisMatrix = getBasis(dim, 2);
+ else if(basis == 3) basisMatrix = getBasis(dim, 3);
+ else return;
+
chain.clear();
if(dim < 0 || dim > 4) return;
- if(_Hbasis[dim] == NULL
- || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return;
+ if(basisMatrix == NULL
+ || (int)gmp_matrix_cols(basisMatrix) < num) return;
- int rows = gmp_matrix_rows(_Hbasis[dim]);
+ int rows = gmp_matrix_rows(basisMatrix);
int elemi;
long int elemli;
mpz_t elem;
mpz_init(elem);
- for(int i = 1; i <= rows; i++){
- gmp_matrix_get_elem(elem, i, chainNumber, _Hbasis[dim]);
+ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+ Cell* cell = cit->first;
+ int index = cit->second;
+ gmp_matrix_get_elem(elem, index, num, basisMatrix);
elemli = mpz_get_si(elem);
elemi = elemli;
- Cell* cell = _cellIndices[dim][i];
- if(cell == NULL){
- //printf("Warning: Missing cell in %d-chain %d. \n", dim, chainNumber);
- continue;
- }
-
- std::map<Cell*, int, Less_Cell > subCells;
- cell->getCells(subCells);
- for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
- Cell* subCell = (*it).first;
- int coeff = (*it).second;
- chain[subCell] = coeff*elemi;
+ if(elemli != 0){
+ std::map<Cell*, int, Less_Cell > subCells;
+ cell->getCells(subCells);
+ for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
+ Cell* subCell = (*it).first;
+ int coeff = (*it).second;
+ chain[subCell] = coeff*elemi;
+ }
}
}
mpz_clear(elem);
}
-int ChainComplex::getTorsion(int dim, int chainNumber)
+int ChainComplex::getBasisSize(int dim, int basis)
+{
+ gmp_matrix* basisMatrix;
+ if(basis == 0 && _HMatrix[dim] != NULL){
+ return gmp_matrix_cols(_HMatrix[dim]);
+ }
+ else if(basis == 1) basisMatrix = getBasis(dim, 1);
+ else if(basis == 2) basisMatrix = getBasis(dim, 2);
+ else if(basis == 3) basisMatrix = getBasis(dim, 3);
+ else return 0;
+
+ if(basisMatrix != NULL) return gmp_matrix_cols(basisMatrix);
+ else return 0;
+}
+
+
+int ChainComplex::getTorsion(int dim, int num)
{
if(dim < 0 || dim > 4) return 0;
if(_Hbasis[dim] == NULL
- || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0;
+ || (int)gmp_matrix_cols(_Hbasis[dim]) < num) return 0;
if(_torsion[dim].empty()
- || (int)_torsion[dim].size() < chainNumber) return 1;
- else return _torsion[dim].at(chainNumber-1);
+ || (int)_torsion[dim].size() < num) return 1;
+ else return _torsion[dim].at(num-1);
+}
+
+
+HomologySequence::HomologySequence(ChainComplex* subcomplex,
+ ChainComplex* complex,
+ ChainComplex* relcomplex)
+{
+ _subcomplex = subcomplex;
+ _complex = complex;
+ _relcomplex = relcomplex;
+
+ mpz_t elem;
+ mpz_init_set_si(elem, -1);
+
+ for(int i = 0; i < 4; i++){
+
+ //printf("Dimension %d. \n", i);
+ _Ic_sub[i] = NULL;
+ _Ic_rel[i] = NULL;
+
+ _Dh[i] = NULL;
+ _invDh[i] = NULL;
+
+ _Jh[i] = NULL;
+ _Ih[i] = NULL;
+ _invJh[i] = NULL;
+ _invIh[i] = NULL;
+
+ mpz_t one;
+ mpz_init_set_si(one, 1);
+ if(_complex->getBasisSize(i, 0) > 0
+ && _subcomplex->getBasisSize(i, 0) > 0){
+ _Ic_sub[i] = create_gmp_matrix_zero(_complex->getBasisSize(i, 0),
+ _subcomplex->getBasisSize(i, 0));
+ //printf("rows %d, cols %d. \n", _complex->getBasisSize(i, 0),
+ // _subcomplex->getBasisSize(i, 0));
+ for(ChainComplex::citer cit = _complex->firstCell(i);
+ cit != _complex->lastCell(i); cit++){
+ Cell* cell = cit->first;
+ int row = cit->second;
+ int col = _subcomplex->cellIndex(cell);
+ //printf("row %d, col %d. \n", row, col);
+ if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_sub[i]);
+ }
+ }
+
+ if(_complex->getBasisSize(i, 0) > 0
+ && _relcomplex->getBasisSize(i, 0) > 0){
+ _Ic_rel[i] = create_gmp_matrix_zero(_complex->getBasisSize(i, 0),
+ _relcomplex->getBasisSize(i, 0));
+ //printf("rows %d, cols %d. \n", _complex->getBasisSize(i, 0),
+ for(ChainComplex::citer cit = _complex->firstCell(i);
+ cit != _complex->lastCell(i); cit++){
+ Cell* cell = cit->first;
+ int row = cit->second;
+ int col = _relcomplex->cellIndex(cell);
+ //printf("row %d, col %d. \n", row, col);
+ if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_rel[i]);
+ }
+ }
+ mpz_clear(one);
+
+ if(_Ic_sub[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i]),
+ gmp_matrix_cols(_Ic_sub[i]));
+ gmp_matrix_right_mult(IH, _subcomplex->getBasis(i, 3));
+ _Ih[i] = createIncMap(IH, _complex->getBasis(i, 3));
+ }
+ if(_Ic_sub[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i]),
+ gmp_matrix_cols(_Ic_sub[i]));
+ gmp_matrix_transp(IH);
+ gmp_matrix_right_mult(IH, _complex->getBasis(i, 3));
+ _invIh[i] = createIncMap(IH, _subcomplex->getBasis(i, 3));
+ }
+
+ if(_Ic_rel[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _relcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_transp(JH);
+ gmp_matrix_right_mult(JH, _complex->getBasis(i, 3));
+ _Jh[i] = createIncMap(JH, _relcomplex->getBasis(i, 3));
+ }
+ if(_Ic_rel[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _relcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_right_mult(JH, _relcomplex->getBasis(i, 3));
+ _invJh[i] = createIncMap(JH, _complex->getBasis(i, 3));
+ }
+
+ //printMatrix(_Ih[i]);
+ //printMatrix(_invIh[i]);
+ //printMatrix(_Jh[i]);
+ //printMatrix(_invJh[i]);
+
+ if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i-1, 3) > 0
+ && _complex->getBoundaryOp(i) != NULL){
+ gmp_matrix* JDIH = copy_gmp_matrix(_Ic_sub[i-1], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i-1]),
+ gmp_matrix_cols(_Ic_sub[i-1]));
+ gmp_matrix_transp(JDIH);
+ gmp_matrix_right_mult(JDIH, _complex->getBoundaryOp(i));
+ gmp_matrix_right_mult(JDIH, _Ic_rel[i]);
+ gmp_matrix_right_mult(JDIH, _relcomplex->getBasis(i, 3));
+ _Dh[i] = createIncMap(JDIH, _subcomplex->getBasis(i-1, 3));
+ }
+
+ if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i-1, 3) > 0
+ && _complex->getBoundaryOp(i) != NULL){
+ gmp_matrix* JDIH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_transp(JDIH);
+ gmp_matrix* bd = _complex->getBoundaryOp(i);
+ gmp_matrix_transp(bd);
+ gmp_matrix_right_mult(JDIH, bd);
+ gmp_matrix_transp(bd);
+ gmp_matrix_right_mult(JDIH, _Ic_sub[i-1]);
+ gmp_matrix_right_mult(JDIH, _subcomplex->getBasis(i-1, 3));
+ _invDh[i] = createIncMap(JDIH, _relcomplex->getBasis(i, 3));
+ }
+
+ //printMatrix(_Dh[i]);
+ //printMatrix(_invDh[i]);
+
+ }
+ for(int i = 0; i < 3; i++){
+ blockHBasis(_Dh[i+1], _invIh[i], _subcomplex, i);
+ blockHBasis(_Ih[i], _invJh[i], _complex, i);
+ blockHBasis(_Jh[i], _invDh[i], _relcomplex, i);
+ }
+
+}
+
+HomologySequence::~HomologySequence()
+{
+ for(int i = 0; i < 4; i++){
+ destroy_gmp_matrix(_Ic_sub[i]);
+ destroy_gmp_matrix(_Ic_rel[i]);
+ destroy_gmp_matrix(_Ih[i]);
+ destroy_gmp_matrix(_Jh[i]);
+ destroy_gmp_matrix(_invIh[i]);
+ destroy_gmp_matrix(_invJh[i]);
+ destroy_gmp_matrix(_Dh[i]);
+ destroy_gmp_matrix(_invDh[i]);
+ }
}
+//i: a->b : aBasis = bBasis*incMap
+gmp_matrix* HomologySequence::createIncMap(gmp_matrix* domBasis,
+ gmp_matrix* codBasis)
+{
+ if(domBasis == NULL || codBasis == NULL){
+ printf("ERROR: null matrix given. \n");
+ return NULL;
+ }
+
+ int rows = gmp_matrix_rows(domBasis);
+ int cols = gmp_matrix_cols(domBasis);
+ if(rows < cols || rows == 0 || cols == 0) return NULL;
+
+ rows = gmp_matrix_rows(codBasis);
+ cols = gmp_matrix_cols(codBasis);
+ if(rows < cols || rows == 0 || cols == 0) return NULL;
+
+ gmp_matrix* temp = codBasis;
+ codBasis = copy_gmp_matrix(temp, 1, 1,
+ gmp_matrix_rows(temp),
+ gmp_matrix_cols(temp));
+ // inv(U)*A*inv(V) = S
+ gmp_normal_form* normalForm
+ = create_gmp_Smith_normal_form(codBasis, 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 NULL;
+ }
+ }
+
+ gmp_matrix_left_mult(normalForm->left, domBasis);
+ gmp_matrix* LB = copy_gmp_matrix(domBasis, 1, 1,
+ gmp_matrix_cols(codBasis),
+ gmp_matrix_cols(domBasis));
+ destroy_gmp_matrix(domBasis);
+
+ 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 NULL;
+ }
+ }
+
+ gmp_matrix_left_mult(normalForm->right, LB);
+
+ mpz_clear(elem);
+ mpz_clear(divisor);
+ mpz_clear(result);
+ destroy_gmp_normal_form(normalForm);
+ return LB;
+}
+
+
+gmp_matrix* HomologySequence::removeZeroCols(gmp_matrix* matrix) // FIXME
+{
+ mpz_t elem;
+ mpz_init(elem);
+
+ int rows = gmp_matrix_rows(matrix);
+ int cols = gmp_matrix_cols(matrix);
+
+ std::vector<int> zcols;
+
+ for(int j = 1; j <= cols; j++){
+ bool zcol = true;
+ for(int i = 1; i <= rows; i++){
+ gmp_matrix_get_elem(elem, i, j, matrix);
+ if(mpz_cmp_si(elem, 0) != 0){
+ zcol = false;
+ break;
+ }
+ }
+ if(zcol) zcols.push_back(j);
+ }
+ if(zcols.empty()) return matrix;
+
+ gmp_matrix* newMatrix = create_gmp_matrix_zero(rows, cols-zcols.size());
+ if(cols-zcols.size() == 0) return newMatrix;
+
+ int k = 0;
+ for(int j = 1; j <= cols; j++){
+ if(zcols.size()-1 < k) break;
+ if(zcols.at(k) == j) { k++; continue; }
+ for(int i = 1; i <= rows; i++){
+ gmp_matrix_get_elem(elem, i, j, matrix);
+ gmp_matrix_set_elem(elem, i, j-k, newMatrix);
+ }
+ }
+
+ destroy_gmp_matrix(matrix);
+ mpz_clear(elem);
+ return newMatrix;
+}
+
+void HomologySequence::blockHBasis(gmp_matrix* block1T,
+ gmp_matrix* block2T,
+ ChainComplex* complex, int dim)
+{
+ if(block1T == NULL && block2T == NULL) return;
+
+ gmp_matrix* Hbasis = complex->getBasis(dim, 3);
+
+ if(block1T == NULL && block2T != NULL){
+ gmp_matrix_right_mult(Hbasis, block2T);
+ return;
+ }
+ if(block1T != NULL && block2T == NULL){
+ gmp_matrix_right_mult(Hbasis, block1T);
+ return;
+ }
+
+ int rows = gmp_matrix_rows(Hbasis);
+ int cols = gmp_matrix_cols(Hbasis);
+ gmp_matrix* temp1 = copy_gmp_matrix(Hbasis, 1, 1, rows, cols);
+ gmp_matrix* temp2 = copy_gmp_matrix(Hbasis, 1, 1, rows, cols);
+
+ gmp_matrix_right_mult(temp1, block1T);
+ gmp_matrix_right_mult(temp2, block2T);
+
+
+ temp1 = removeZeroCols(temp1);
+ temp2 = removeZeroCols(temp2);
+ //printMatrix(temp1);
+ //printMatrix(temp2);
+
+ int bcol = gmp_matrix_cols(temp1);
+
+ mpz_t elem;
+ mpz_init(elem);
+ for(int i = 1; i <= rows; i++){
+ for(int j = 1; j <= cols; j++){
+ if(j <= bcol) gmp_matrix_get_elem(elem, i, j, temp1);
+ else gmp_matrix_get_elem(elem, i, j-bcol, temp2);
+ gmp_matrix_set_elem(elem, i, j, Hbasis);
+ }
+ }
+
+ //printMatrix(Hbasis);
+ mpz_clear(elem);
+ destroy_gmp_matrix(temp1);
+ destroy_gmp_matrix(temp2);
+}
+
+
#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index ef2308128e8b82c9285c2e803fd45b76c67debb9..dcf9b2fe3923202566cb2eb076864be3166aa186 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -33,7 +33,8 @@ class CellComplex;
// This should only be constructed for a reduced cell complex because of
// dense matrix representations and great computational complexity in
// its methods.
-class ChainComplex{
+class ChainComplex
+{
private:
// boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
@@ -44,17 +45,23 @@ class ChainComplex{
gmp_matrix* _kerH[5];
gmp_matrix* _codH[5];
+ // matrix of the mapping B_k -> Z_k
gmp_matrix* _JMatrix[5];
+ // matrix of the mapping H_k -> Z_k
gmp_matrix* _QMatrix[5];
- std::vector<long int> _torsion[5];
-
+
// bases for homology groups
gmp_matrix* _Hbasis[5];
+ // torsion coefficients of homology generators
+ // corresponding the columns of _Hbasis
+ std::vector<long int> _torsion[5];
int _dim;
CellComplex* _cellComplex;
- std::map<int, Cell*> _cellIndices[4];
+ // index to cell map
+ // matrix indices correspond to these cells in _cellComplex
+ std::map<Cell*, int, Less_Cell> _cellIndices[4];
// set the matrices
void setHMatrix(int dim, gmp_matrix* matrix) {
@@ -69,49 +76,51 @@ class ChainComplex{
if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
void setHbasis(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
-
- public:
-
- // domain = 0 : relative chain space
- // domain = 1 : absolute chain space of all cells in cellComplex
- // domain = 2 : absolute chain space of cells in subdomain
- ChainComplex(CellComplex* cellComplex, int domain=0);
- ~ChainComplex(){
- for(int i = 0; i < 5; i++){
- destroy_gmp_matrix(_HMatrix[i]);
- destroy_gmp_matrix(_kerH[i]);
- destroy_gmp_matrix(_codH[i]);
- destroy_gmp_matrix(_JMatrix[i]);
- destroy_gmp_matrix(_QMatrix[i]);
- destroy_gmp_matrix(_Hbasis[i]);
- }
- }
-
- int getDim() { return _dim; }
-
+
+
// get the boundary operator matrix dim->dim-1
- gmp_matrix* getHMatrix(int dim) {
+ gmp_matrix* getHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
- gmp_matrix* getKerHMatrix(int dim) {
+ gmp_matrix* getKerHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
- gmp_matrix* getCodHMatrix(int dim) {
+ gmp_matrix* getCodHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
- gmp_matrix* getJMatrix(int dim) {
+ gmp_matrix* getJMatrix(int dim) const {
if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
- gmp_matrix* getQMatrix(int dim) {
+ gmp_matrix* getQMatrix(int dim) const {
if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
- gmp_matrix* getHbasis(int dim) {
+ gmp_matrix* getHbasis(int dim) const {
if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
+
+
+ public:
+ // domain = 0 : relative chain space
+ // domain = 1 : absolute chain space of all cells in cellComplex
+ // domain = 2 : absolute chain space of cells in subdomain
+ ChainComplex(CellComplex* cellComplex, int domain=0);
+ ~ChainComplex();
+
+ int getDim() const { return _dim; }
- // Compute basis for kernel and codomain of boundary operator matrix
+ // 0 : C basis (chains)
+ // 1 : Z basis (cycles)
+ // 2 : B basis (boundaries)
+ // 3 : H basis (homology)
+ // get the bases for various spaces
+ gmp_matrix* getBasis(int dim, int basis);
+ gmp_matrix* getBoundaryOp(int dim) {
+ if(dim > -1 && dim < 5) return _HMatrix[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) )
void KerCod(int dim);
- // Compute matrix representation J for inclusion relation from dim-cells
+ // 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) )
void Inclusion(int lowDim, int highDim);
- // Compute quotient problem for given inclusion relation j to find
- //representatives of homology groups and possible torsion coeffcients
+ // compute quotient problem for given inclusion relation j to find
+ // representatives of homology group generators and possible
+ // torsion coeffcients
void Quotient(int dim);
// transpose the boundary operator matrices, these are boundary operator
@@ -124,26 +133,89 @@ class ChainComplex{
// Compute bases for the homology groups of this chain complex
void computeHomology(bool dual=false);
- // get coefficient vector for dim-dimensional chain chainNumber
- // (columns of _Hbasis[dim])
+
+ typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+ citer firstCell(int dim){ return _cellIndices[dim].begin(); }
+ citer lastCell(int dim){ return _cellIndices[dim].end(); }
+ // get the cell index
+ int cellIndex(Cell* cell){
+ citer cit = _cellIndices[cell->getDim()].find(cell);
+ if(cit != lastCell(cell->getDim())) return cit->second;
+ else return 0;
+ }
+
+ // get coefficient vector for dim-dimensional Hbasis chain chainNumber
std::vector<int> getCoeffVector(int dim, int chainNumber);
- void getChain(std::map<Cell*, int, Less_Cell>& chain, int dim, int chainNumber);
- // torsion coefficient for dim-dimensional chain chainNumber
- int getTorsion(int dim, int chainNumber);
- // get rank of homology group of dimension dim
- int getBasisSize(int dim) {
- if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]);
- else return 0; }
+ // get basis chain from a basis matrix
+ void getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
+ int num, int dim, int basis);
+ // get rank of a basis
+ int getBasisSize(int dim, int basis);
+ // homology torsion coefficient for dim-dimensional chain num
+ int getTorsion(int dim, int num);
+
+ // apply a transformation T to a basis (T should be unimodular)
+ void transformBasis(gmp_matrix* T, int dim, int basis){
+ if(basis == 3 && _Hbasis[dim] != NULL) {
+ gmp_matrix_right_mult(_Hbasis[dim], T);
+ }
+ }
+ // debugging aid
int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n",
- (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
+ (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
return gmp_matrix_printf(matrix); }
-
- // debugging aid
void matrixTest();
};
+
+// An experimental class to modify computed bases for homnology spaces
+// so that the basis chains are decomposed according to the long
+// exact homology sequence.
+class HomologySequence
+{
+ private:
+ ChainComplex* _subcomplex;
+ ChainComplex* _complex;
+ ChainComplex* _relcomplex;
+
+ gmp_matrix* _Ic_sub[4];
+ gmp_matrix* _Ic_rel[4];
+
+ gmp_matrix* _Ih[4];
+ gmp_matrix* _Jh[4];
+ gmp_matrix* _invIh[4];
+ gmp_matrix* _invJh[4];
+
+
+ gmp_matrix* _Dh[4];
+ gmp_matrix* _invDh[4];
+
+ public:
+
+ HomologySequence(ChainComplex* subcomplex, ChainComplex* complex,
+ ChainComplex* relcomplex);
+ ~HomologySequence();
+
+ // create an inclusion map from domBasis to codBasis
+ // (deletes domBasis, leaves codBasis unaffected)
+ gmp_matrix* createIncMap(gmp_matrix* domBasis,
+ gmp_matrix* codBasis);
+
+
+ gmp_matrix* removeZeroCols(gmp_matrix* matrix);
+ void blockHBasis(gmp_matrix* block1T, gmp_matrix* block2T,
+ ChainComplex* complex, int dim);
+
+ int printMatrix(gmp_matrix* matrix){
+ if(matrix == NULL){ printf("NULL matrix. \n"); return 0; }
+ printf("%d rows and %d columns\n",
+ (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
+ return gmp_matrix_printf(matrix); }
+
+};
+
#endif
#endif
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index c0e571a8f8252882fcdfe86c7d80531bf746340b..eef7aca4cb1be41eb8711ba1416128cab67b88e4 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -11,7 +11,7 @@
#if defined(HAVE_KBIPACK)
Homology::Homology(GModel* model, std::vector<int> physicalDomain,
- std::vector<int> physicalSubdomain)
+ std::vector<int> physicalSubdomain)
{
_model = model;
_domain = physicalDomain;
@@ -37,10 +37,10 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
for(int j = 0; j < 4; j++){
it = groups[j].find(_domain.at(i));
if(it != groups[j].end()){
- std::vector<GEntity*> physicalGroup = (*it).second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- _domainEntities.push_back(physicalGroup.at(k));
- }
+ std::vector<GEntity*> physicalGroup = (*it).second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ _domainEntities.push_back(physicalGroup.at(k));
+ }
}
}
}
@@ -48,10 +48,10 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
for(int j = 0; j < 4; j++){
it = groups[j].find(_subdomain.at(i));
if(it != groups[j].end()){
- std::vector<GEntity*> physicalGroup = (*it).second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- _subdomainEntities.push_back(physicalGroup.at(k));
- }
+ std::vector<GEntity*> physicalGroup = (*it).second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ _subdomainEntities.push_back(physicalGroup.at(k));
+ }
}
}
}
@@ -59,7 +59,7 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
}
CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
- std::vector<GEntity*>& subdomainEntities){
+ std::vector<GEntity*>& subdomainEntities){
Msg::Info("Creating a Cell Complex...");
Msg::StatusBar(1, false, "Cell Complex...");
Msg::StatusBar(2, false, "");
@@ -79,14 +79,14 @@ CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
}
for(unsigned int j=0; j < subdomainEntities.size(); j++) {
for(unsigned int i=0; i < subdomainEntities.at(j)->getNumMeshElements();
- i++){
+ i++){
MElement* element = subdomainEntities.at(j)->getMeshElement(i);
subdomainElements.push_back(element);
}
}
CellComplex* cellComplex = new CellComplex(domainElements,
- subdomainElements);
+ subdomainElements);
if(cellComplex->getSize(0) == 0){
Msg::Error("Cell Complex is empty!");
@@ -96,10 +96,10 @@ CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
Msg::Info("Cell Complex complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::StatusBar(2, false, "%d V, %d F, %d E, %d V.",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
return cellComplex;
}
@@ -113,7 +113,7 @@ void Homology::findGenerators(CellComplex* cellComplex)
bool ownComplex = false;
if(cellComplex==NULL){
cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
+ _subdomainEntities);
ownComplex = true;
}
std::string domainString = getDomainString(_domain, _subdomain);
@@ -127,10 +127,10 @@ void Homology::findGenerators(CellComplex* cellComplex)
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
@@ -145,29 +145,33 @@ void Homology::findGenerators(CellComplex* cellComplex)
HRank[j] = 0;
std::string dimension = "";
convert(j, dimension);
- for(int i = 1; i <= chains->getBasisSize(j); i++){
+ for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
std::string generator = "";
convert(i, generator);
std::string name = "H" + dimension + domainString + generator;
- std::set<Cell*, Less_Cell> cells;
+ std::map<Cell*, int, Less_Cell> protoChain;
+ chains->getBasisChain(protoChain, i, j, 3);
+ Chain* chain = new Chain(protoChain, cellComplex, _model, name,
+ chains->getTorsion(j,i));
+ /* std::set<Cell*, Less_Cell> cells;
cellComplex->getCells(cells, j);
Chain* chain = new Chain(cells,
chains->getCoeffVector(j,i), cellComplex,
- _model, name, chains->getTorsion(j,i));
+ _model, name, chains->getTorsion(j,i));*/
t1 = Cpu();
int start = chain->getSize();
chain->smoothenChain();
t2 = Cpu();
Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
- j, i, start, chain->getSize(), t2 - t1);
+ j, i, start, chain->getSize(), t2 - t1);
if(chain->getSize() != 0) {
HRank[j] = HRank[j] + 1;
if(chain->getTorsion() != 1){
- Msg::Warning("H%d %d has torsion coefficient %d!",
- j, i, chain->getTorsion());
- }
+ Msg::Warning("H%d %d has torsion coefficient %d!",
+ j, i, chain->getTorsion());
+ }
}
_generators.push_back(chain->createPGroup());
delete chain;
@@ -179,10 +183,10 @@ void Homology::findGenerators(CellComplex* cellComplex)
std::string name = "H" + dimension + domainString + generator;
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
- cellComplex, _model, name, 1);
+ cellComplex, _model, name, 1);
if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
- _generators.push_back(chain->createPGroup());
- delete chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
}
@@ -201,7 +205,7 @@ void Homology::findGenerators(CellComplex* cellComplex)
Msg::StatusBar(1, false, "Homology");
Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.",
- HRank[0], HRank[1], HRank[2], HRank[3]);
+ HRank[0], HRank[1], HRank[2], HRank[3]);
}
void Homology::findDualGenerators(CellComplex* cellComplex)
@@ -209,7 +213,7 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
bool ownComplex = false;
if(cellComplex==NULL){
cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
+ _subdomainEntities);
ownComplex = true;
}
@@ -223,10 +227,10 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
@@ -238,33 +242,37 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
Msg::Info("Homology Computation complete (%g s).", t2- t1);
int dim = cellComplex->getDim();
-
+
int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0;
for(int j = 3; j > -1; j--){
std::string dimension = "";
convert(dim-j, dimension);
- for(int i = 1; i <= chains->getBasisSize(j); i++){
+ for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
std::string generator = "";
convert(i, generator);
std::string name = "H" + dimension + "*" +
- getDomainString(_domain, _subdomain) + generator;
- std::set<Cell*, Less_Cell> cells;
+ getDomainString(_domain, _subdomain) + generator;
+ std::map<Cell*, int, Less_Cell> protoChain;
+ chains->getBasisChain(protoChain, i, j, 3);
+ Chain* chain = new Chain(protoChain, cellComplex, _model, name,
+ chains->getTorsion(j,i));
+ /*std::set<Cell*, Less_Cell> cells;
cellComplex->getCells(cells, j);
Chain* chain = new Chain(cells,
- chains->getCoeffVector(j,i), cellComplex,
- _model, name, chains->getTorsion(j,i));
+ chains->getCoeffVector(j,i), cellComplex,
+ _model, name, chains->getTorsion(j,i));*/
_generators.push_back(chain->createPGroup());
delete chain;
if(chain->getSize() != 0){
HRank[dim-j] = HRank[dim-j] + 1;
if(chain->getTorsion() != 1){
- Msg::Warning("H%d* %d has torsion coefficient %d!",
- dim-j, i, chain->getTorsion());
- }
+ Msg::Warning("H%d* %d has torsion coefficient %d!",
+ dim-j, i, chain->getTorsion());
+ }
}
}
@@ -274,14 +282,14 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
std::string generator;
convert(i+1, generator);
std::string name
- = "H" + dimension + "*" +
- getDomainString(_domain, _subdomain) + generator;
+ = "H" + dimension + "*" +
+ getDomainString(_domain, _subdomain) + generator;
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
- cellComplex, _model, name, 1);
+ cellComplex, _model, name, 1);
if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1;
- _generators.push_back(chain->createPGroup());
- delete chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
}
@@ -300,12 +308,12 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
Msg::StatusBar(1, false, "Homology");
Msg::StatusBar(2, false, "H0*: %d, H1*: %d, H2*: %d, H3*: %d.",
- HRank[0], HRank[1], HRank[2], HRank[3]);
+ HRank[0], HRank[1], HRank[2], HRank[3]);
}
void Homology::findHomSequence(){
CellComplex* cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
+ _subdomainEntities);
Msg::Info("Reducing the Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
@@ -316,58 +324,79 @@ void Homology::findHomSequence(){
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
- cellComplex->getSize(3), cellComplex->getSize(2),
- cellComplex->getSize(1), cellComplex->getSize(0));
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
t1 = Cpu();
+ ChainComplex* subcomplex = new ChainComplex(cellComplex, 2);
+ ChainComplex* complex = new ChainComplex(cellComplex, 1);
+ ChainComplex* relcomplex = new ChainComplex(cellComplex, 0);
+
+ subcomplex->computeHomology();
+ complex->computeHomology();
+ relcomplex->computeHomology();
+
+ t2 = Cpu();
+ Msg::Info("Homology computation complete (%g s).", t2 - t1);
+
+ Msg::Info("Computing homology sequence...");
+ HomologySequence* seq = new HomologySequence(subcomplex,
+ complex, relcomplex);
+ t1 = Cpu();
+ Msg::Info("Homology sequence computation complete (%g s).", t1 - t2);
+
for(int task = 0; task < 3; task++){
- ChainComplex* chains = new ChainComplex(cellComplex, task);
+ ChainComplex* chains;
+
std::string domainString = "";
std::vector<int> empty;
- if(task == 0) domainString = getDomainString(_domain, _subdomain);
- else if(task == 1) domainString = getDomainString(_domain, empty);
- else if(task == 2) domainString = getDomainString(_subdomain, empty);
- chains->computeHomology();
- t2 = Cpu();
- Msg::Info("Homology Computation complete (%g s).", t2 - t1);
+ if(task == 0) {
+ chains = subcomplex;
+ domainString = getDomainString(_subdomain, empty);
+ }
+ else if(task == 1) {
+ chains = complex;
+ domainString = getDomainString(_domain, empty);
+ }
+ else{
+ chains = relcomplex;
+ domainString = getDomainString(_domain, _subdomain);
+ }
int HRank[4];
for(int j = 0; j < 4; j++){
HRank[j] = 0;
std::string dimension = "";
convert(j, dimension);
- for(int i = 1; i <= chains->getBasisSize(j); i++){
-
- std::string generator = "";
- convert(i, generator);
+ for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
+
+ std::string generator = "";
+ convert(i, generator);
- std::string name = "H" + dimension + domainString + generator;
- std::set<Cell*, Less_Cell> cells;
- cellComplex->getCells(cells, j, task);
- Chain* chain = new Chain(cells,
- chains->getCoeffVector(j,i),
- cellComplex, _model, name,
- chains->getTorsion(j,i));
- t1 = Cpu();
- int start = chain->getSize();
- //chain->smoothenChain();
- t2 = Cpu();
- //Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
- // j, i, start, chain->getSize(), t2 - t1);
- if(chain->getSize() != 0) {
- HRank[j] = HRank[j] + 1;
- if(chain->getTorsion() != 1){
- Msg::Warning("H%d %d has torsion coefficient %d!",
- j, i, chain->getTorsion());
- }
- }
- _generators.push_back(chain->createPGroup());
- delete chain;
+ std::string name = "H" + dimension + domainString + generator;
+ std::map<Cell*, int, Less_Cell> protoChain;
+ chains->getBasisChain(protoChain, i, j, 3);
+ Chain* chain = new Chain(protoChain, cellComplex, _model, name,
+ chains->getTorsion(j,i));
+ if(chain->getSize() == 0) continue;
+ t1 = Cpu();
+ int start = chain->getSize();
+ //chain->smoothenChain();
+ t2 = Cpu();
+ //Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
+ // j, i, start, chain->getSize(), t2 - t1);
+ HRank[j] = HRank[j] + 1;
+ if(chain->getTorsion() != 1){
+ Msg::Warning("H%d %d has torsion coefficient %d!",
+ j, i, chain->getTorsion());
+ }
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
@@ -387,16 +416,19 @@ void Homology::findHomSequence(){
Msg::StatusBar(1, false, "Homology");
Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.",
- HRank[0], HRank[1], HRank[2], HRank[3]);
- delete chains;
+ HRank[0], HRank[1], HRank[2], HRank[3]);
}
+ delete subcomplex;
+ delete complex;
+ delete relcomplex;
+
if(_fileName != "") writeGeneratorsMSH();
delete cellComplex;
}
std::string Homology::getDomainString(const std::vector<int>& domain,
- const std::vector<int>& subdomain)
+ const std::vector<int>& subdomain)
{
std::string domainString = "({";
if(domain.empty()) domainString += "0";
@@ -406,7 +438,7 @@ std::string Homology::getDomainString(const std::vector<int>& domain,
convert(domain.at(i),temp);
domainString += temp;
if (domain.size()-1 > i){
- domainString += ", ";
+ domainString += ", ";
}
}
}
@@ -436,8 +468,8 @@ bool Homology::writeGeneratorsMSH(bool binary)
return true;
}
Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
- CellComplex* cellComplex, GModel* model,
- std::string name, int torsion)
+ CellComplex* cellComplex, GModel* model,
+ std::string name, int torsion)
{
int i = 0;
for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin();
@@ -447,7 +479,7 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
if((int)coeffs.size() > i){
if(coeffs.at(i) != 0){
std::map<Cell*, int, Less_Cell > subCells;
- cell->getCells(subCells);
+ cell->getCells(subCells);
for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
Cell* subCell = (*it).first;
int coeff = (*it).second;
@@ -502,8 +534,8 @@ void Homology::storeCells(CellComplex* cellComplex, int dim)
}
Chain::Chain(std::map<Cell*, int, Less_Cell>& chain,
- CellComplex* cellComplex, GModel* model,
- std::string name, int torsion)
+ CellComplex* cellComplex, GModel* model,
+ std::string name, int torsion)
{
_cells = chain;
if(!_cells.empty()) _dim = firstCell()->first->getDim();
@@ -516,7 +548,7 @@ Chain::Chain(std::map<Cell*, int, Less_Cell>& chain,
}
bool Chain::deform(std::map<Cell*, int, Less_Cell>& cellsInChain,
- std::map<Cell*, int, Less_Cell>& cellsNotInChain)
+ std::map<Cell*, int, Less_Cell>& cellsNotInChain)
{
std::vector<int> cc;
std::vector<int> bc;
@@ -565,12 +597,12 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend)
Cell* c1CbdCell = (*cit).first;
for(citer cit2 = c1CbdCell->firstBoundary(true);
- cit2 != c1CbdCell->lastBoundary(true); cit2++){
+ cit2 != c1CbdCell->lastBoundary(true); cit2++){
Cell* c1CbdBdCell = (*cit2).first;
int coeff = (*cit2).second;
if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0)
- || c1CbdBdCell->inSubdomain()){
- cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ || c1CbdBdCell->inSubdomain()){
+ cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
}
else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
}
@@ -586,30 +618,30 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend)
}
for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
- cit2++){
+ cit2++){
Cell* c = (*cit2).first;
if(c->inSubdomain()) next = true;
}
if(next) continue;
if( (getDim() == 1 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 1) ||
- (getDim() == 2 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 1)){
+ && cellsNotInChain.size() == 1) ||
+ (getDim() == 2 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 1)){
//printf("straighten \n");
return deform(cellsInChain, cellsNotInChain);
}
else if ( (getDim() == 1 && cellsInChain.size() == 1
- && cellsNotInChain.size() == 2 && bend) ||
+ && cellsNotInChain.size() == 2 && bend) ||
(getDim() == 2 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 2 && bend)){
+ && cellsNotInChain.size() == 2 && bend)){
//printf("bend \n");
return deform(cellsInChain, cellsNotInChain);
}
else if ((getDim() == 1 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 0) ||
+ && cellsNotInChain.size() == 0) ||
(getDim() == 2 && cellsInChain.size() == 4
- && cellsNotInChain.size() == 0)){
+ && cellsNotInChain.size() == 0)){
//printf("remove boundary \n");
return deform(cellsInChain, cellsNotInChain);
}
@@ -648,7 +680,7 @@ void Chain::smoothenChain()
eraseNullCells();
double t2 = Cpu();
Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n",
- getDim(), start, getSize(), t2-t1);
+ getDim(), start, getSize(), t2-t1);
return;
}
@@ -747,7 +779,7 @@ int Chain::createPGroup()
// create PView for visualization
PView* chain = new PView(getName(), "ElementData", getGModel(),
- data, 0, 1);
+ data, 0, 1);
}
return physicalNum;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 9bd285f598362c063d9e5a637b934548e096957b..2f7298ea0cef021627cee8d4f0dd3c89e1045288 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -33,57 +33,60 @@ class Homology
{
private:
- // the model of the homology computation
+ // the Gmsh model for homology computation
GModel* _model;
// domain and the relative subdomain of the homology computation
+ // physical group IDs
std::vector<int> _domain;
std::vector<int> _subdomain;
+ // corresponding geometrical entities
std::vector<GEntity*> _domainEntities;
std::vector<GEntity*> _subdomainEntities;
- // physical group numbers of generator chains in model
+ // physical group numbers of resulted generator chains in model
std::vector<int> _generators;
+ // file name to store the results
std::string _fileName;
public:
Homology(GModel* model, std::vector<int> physicalDomain,
- std::vector<int> physicalSubdomain);
+ std::vector<int> physicalSubdomain);
~Homology();
-
+ // create a cell complex from a mesh in geometrical entities of Gmsh
CellComplex* createCellComplex(std::vector<GEntity*>& domainEntities,
- std::vector<GEntity*>& subdomainEntities);
+ std::vector<GEntity*>& subdomainEntities);
CellComplex* createCellComplex() {
return createCellComplex(_domainEntities, _subdomainEntities); }
void setFileName(std::string fileName) { _fileName = fileName; }
- // Find the generators/duals of homology spaces,
- // or just compute the ranks of homology spaces
+ // find the generators/dual generarators of homology spaces,
void findGenerators(CellComplex* cellComplex=NULL);
void findDualGenerators(CellComplex* cellComplex=NULL);
- void computeRanks() {}
- void findHomSequence();
+ // experimental
+ void findHomSequence();
+ void computeRanks() {}
- // Create a string describing the generator
+ // create a string describing the generator
std::string getDomainString(const std::vector<int>& domain,
- const std::vector<int>& subdomain);
+ const std::vector<int>& subdomain);
// write the generators to a file
bool writeGeneratorsMSH(bool binary=false);
- // store dim-dimensional cells of cellComplex as a pshysical group
+ // store dim-dimensional cells of cellComplex as a physical group
// in _model, for debugging
void storeCells(CellComplex* cellComplex, int dim);
};
// A class representing a chain.
-// Used to visualize generators of the homology groups.
+// Used to visualize and post-process generators of the homology spaces in Gmsh.
class Chain{
private:
@@ -103,20 +106,20 @@ class Chain{
int _dim;
bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
- std::map<Cell*, int, Less_Cell> &cellsNotInChain);
+ std::map<Cell*, int, Less_Cell> &cellsNotInChain);
bool deformChain(std::pair<Cell*, int> cell, bool bend);
public:
Chain() {}
Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
- CellComplex* cellComplex, GModel* model,
- std::string name="Chain", int torsion=0);
+ CellComplex* cellComplex, GModel* model,
+ std::string name="Chain", int torsion=0);
Chain(std::map<Cell*, int, Less_Cell>& chain,
CellComplex* cellComplex, GModel* model,
std::string name="Chain", int torsion=0);
Chain(Chain* chain){
- _cells = chain->getCells();
+ chain->getCells(_cells);
_torsion = chain->getTorsion();
_name = chain->getName();
_cellComplex = chain->getCellComplex();
@@ -141,25 +144,26 @@ class Chain{
int getCoeff(Cell* c);
- int getTorsion() { return _torsion; }
- int getDim() { return _dim; }
- CellComplex* getCellComplex() { return _cellComplex; }
- GModel* getGModel() { return _model; }
- std::map<Cell*, int, Less_Cell> getCells() { return _cells; }
+ int getTorsion() const { return _torsion; }
+ int getDim() const { return _dim; }
+ CellComplex* getCellComplex() const { return _cellComplex; }
+ GModel* getGModel() const { return _model; }
+ void getCells(std::map<Cell*, int, Less_Cell> cells) const{ cells = _cells; }
// erase cells from the chain with zero coefficient
void eraseNullCells();
-
+
+ // set all cell in chain not immune
void deImmuneCells();
// number of cells in this chain
- int getSize() { return _cells.size();}
+ int getSize() const { return _cells.size();}
// get/set chain name
- std::string getName() { return _name; }
+ std::string getName() const { return _name; }
void setName(std::string name) { _name=name; }
// get/set physical group number
- int getNum() { return _num; }
+ int getNum() const { return _num; }
void setNum(int num) { _num=num; }
// make local deformations to the chain to make it smoother and smaller
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index c690292c6c70a4c07757ea8a974761aae76f4938..a9a4f9eeee77aef68f069eb842ed68f63895c769 100644
--- a/contrib/kbipack/gmp_matrix.c
+++ b/contrib/kbipack/gmp_matrix.c
@@ -619,6 +619,36 @@ gmp_matrix_transp(gmp_matrix * M)
return EXIT_SUCCESS;
}
+int
+gmp_matrix_sum(gmp_matrix * A, const gmp_matrix * B)
+{
+ size_t i,j;
+ size_t rows_A, cols_A, rows_B, cols_B;
+
+ if((A == NULL) || (B == NULL))
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows_A = A->rows;
+ cols_A = A->cols;
+ rows_B = B->rows;
+ cols_B = B->cols;
+
+ if(cols_A != cols_B || rows_A != rows_B)
+ {
+ return EXIT_FAILURE;
+ }
+ mpz_t a;
+ mpz_init_set_si(a, 1);
+ /* Compute the sum */
+ for(i = 1; i <= rows_A; i++)
+ {
+
+ }
+
+ return EXIT_SUCCESS;
+}
int
gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B)
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
index 931f3828036a15ed73221d9d54a43cdddabeebe5..318dd85657c06d91b0bea17a99be6563d220af2d 100644
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -127,6 +127,10 @@ gmp_matrix_row_inz (size_t r, size_t c1, size_t c2,
int
gmp_matrix_transp(gmp_matrix * M);
+/* A <- A+B */
+int
+gmp_matrix_sum(gmp_matrix * A, const gmp_matrix * B);
+
/* A <- A*B */
int
gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B);