diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 01164cc1e492994adadc26f1f2da3c02e1351982..0494ade2e52c493d60f2abdbbf24c2ecacc89b26 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -17,6 +17,58 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_dim = 0;
+ // find boundary elements
+ find_boundary(domain, subdomain);
+
+ // insert cells into cell complex
+ // subdomain need to be inserted first!
+ insert_cells(true, true);
+ insert_cells(false, true);
+ insert_cells(false, false);
+
+ // find vertices in domain
+ find_vertices(domain, subdomain);
+
+
+ //int tag = 1;
+ for(int i = 0; i < 4; i++){
+ /*
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ cell->setTag(tag);
+ tag++;
+ }*/
+ _ocells[i] = _cells[i];
+ _betti[i] = 0;
+ if(getSize(i) > _dim) _dim = i;
+ }
+
+
+}
+
+void CellComplex::find_vertices(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
+ // find mesh vertices in the domain
+ for(unsigned int j=0; j < domain.size(); j++) {
+ for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
+ for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
+ MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
+ _domainVertices.insert(vertex);
+ }
+ }
+ }
+
+ for(unsigned int j=0; j < subdomain.size(); j++) {
+ for(unsigned int i=0; i < subdomain.at(j)->getNumMeshElements(); i++){
+ for(int k=0; k < subdomain.at(j)->getMeshElement(i)->getNumVertices(); k++){
+ MVertex* vertex = subdomain.at(j)->getMeshElement(i)->getVertex(k);
+ _domainVertices.insert(vertex);
+ }
+ }
+ }
+}
+
+void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
+
// determine mesh entities on boundary of the domain
bool duplicate = false;
for(unsigned int j=0; j < _domain.size(); j++){
@@ -72,46 +124,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
}
}
- // insert cells into cell complex
- // subdomain need to be inserted first!
- insert_cells(true, true);
- insert_cells(false, true);
- insert_cells(false, false);
-
- // find mesh vertices in the domain
- for(unsigned int j=0; j < domain.size(); j++) {
- for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
- MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
- _domainVertices.insert(vertex);
- }
- }
- }
-
- for(unsigned int j=0; j < subdomain.size(); j++) {
- for(unsigned int i=0; i < subdomain.at(j)->getNumMeshElements(); i++){
- for(int k=0; k < subdomain.at(j)->getMeshElement(i)->getNumVertices(); k++){
- MVertex* vertex = subdomain.at(j)->getMeshElement(i)->getVertex(k);
- _domainVertices.insert(vertex);
- }
- }
- }
-
-
- // attach individual tags to cells
- int tag = 1;
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- //Cell* cell = *cit;
- //cell->setTag(tag);
- tag++;
- }
- _ocells[i] = _cells[i];
- _betti[i] = 0;
- if(getSize(i) > _dim) _dim = i;
- }
-
-
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
@@ -852,6 +864,7 @@ int CellComplex::combine(int dim){
}
*/
+/*
void CellComplex::swapSubdomain(){
for(int i = 0; i < 4; i++){
@@ -877,7 +890,7 @@ void CellComplex::swapSubdomain(){
}
return;
-}
+}*/
int CellComplex::writeComplexMSH(const std::string &name){
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index d8ed040f9a8eb7fed3ceeb7f46e577127feab4ea..742d85a49d1ec20193975a845618e149d9532ca5 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#ifndef _CELLCOMPLEX_H_
#define _CELLCOMPLEX_H_
@@ -48,7 +48,7 @@ class Cell
// cell dimension
int _dim;
- // whether this cell belongs to a subdomain
+ // whether this cell belongs to a subdomain, immutable
// used in relative homology computation
bool _inSubdomain;
@@ -58,19 +58,22 @@ class Cell
// whether this cell a combinded cell of elemetary cells
bool _combined;
- // unique tag for each cell
+ // unused
int _tag;
+
+ // mutable index for each cell
int _index;
+ // for some algorithms to omit this cell
bool _immune;
- // cells on the boundary and on the coboundary of this cell
+ // mutable list of cells on the boundary and on the coboundary of this cell
std::map< Cell*, int, Less_Cell > _boundary;
std::map< Cell*, int, Less_Cell > _coboundary;
int _bdSize;
int _cbdSize;
- // original boundary and coboundary before the reduction of the cell complex
+ // immutable original boundary and coboundary before the reduction of the cell complex
std::map<Cell*, int, Less_Cell > _obd;
std::map<Cell*, int, Less_Cell > _ocbd;
@@ -99,6 +102,7 @@ class Cell
virtual int getSortedVertex(int vertex) const = 0;
virtual std::vector<MVertex*> getVertexVector() const = 0;
+ // restores the cell information to its original state before reduction
virtual void restoreCell(){
_boundary = _obd;
_coboundary = _ocbd;
@@ -110,19 +114,16 @@ class Cell
}
- // returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
- // otherwise returns 0
- // implementation will vary depending on cell type
- //virtual int kappa(Cell* tau) const = 0;
-
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const = 0;
// (co)boundary cell iterator
typedef std::map<Cell*, int, Less_Cell>::iterator biter;
- virtual int getBoundarySize() { return _bdSize; }
- virtual int getCoboundarySize() { return _cbdSize; }
+ //virtual int getBoundarySize() { return _bdSize; }
+ //virtual int getCoboundarySize() { return _cbdSize; }
+ virtual int getBoundarySize() { return _boundary.size(); }
+ virtual int getCoboundarySize() { return _coboundary.size(); }
virtual std::map<Cell*, int, Less_Cell > getOrientedBoundary() { return _boundary; }
virtual std::list< Cell* > getBoundary() {
@@ -156,21 +157,6 @@ class Cell
return true;
}
_boundary.insert( std::make_pair(cell, orientation ) );
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- Cell* cell2 = (*it).second;
- if(*cell2 == *cell) {
- (*it).first = (*it).first + orientation;
- if((*it).first == 0) {
- _boundary.erase(it); _bdSize--;
- cell2->removeCoboundaryCell(this,false);
- return false;
- }
- return true;
- }
- }
- _boundary.push_back( std::make_pair(orientation, cell ) );
- */
return true;
}
@@ -198,22 +184,11 @@ class Cell
_bdSize--;
return (*it).second;
}
-
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- Cell* cell2 = (*it).second;
- int ori = (*it).first;
- if(*cell2 == *cell) {
- _boundary.erase(it);
- if(other) cell2->removeCoboundaryCell(this, false);
- _bdSize--;
- return ori;
- }
- }
- */
return 0;
}
+
+
virtual int removeCoboundaryCell(Cell* cell, bool other=true) {
biter it = _coboundary.find(cell);
if(it != _coboundary.end()){
@@ -225,15 +200,9 @@ class Cell
return 0;
}
+ // true if has given cell on (original) boundary
virtual bool hasBoundary(Cell* cell, bool org=false){
if(!org){
- /*
- for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
- Cell* cell2 = (*it).second;
- if(*cell2 == *cell) return true;
- }
- return false;
- */
biter it = _boundary.find(cell);
if(it != _boundary.end()) return true;
return false;
@@ -244,6 +213,8 @@ class Cell
return false;
}
}
+
+ // true if has given cell on (original) boundary
virtual bool hasCoboundary(Cell* cell, bool org=false){
if(!org){
biter it = _coboundary.find(cell);
@@ -261,6 +232,7 @@ class Cell
virtual void clearBoundary() { _boundary.clear(); }
virtual void clearCoboundary() { _coboundary.clear(); }
+ // algebraic dual of the cell
virtual void makeDualCell(){
std::map<Cell*, int, Less_Cell > temp = _boundary;
_boundary = _coboundary;
@@ -305,13 +277,11 @@ class Cell
cell2->printCell();
if(_ocbd.empty()) printf("Cell coboundary is empty. \n");
}
- }
-
-
+ }
virtual bool inSubdomain() const { return _inSubdomain; }
- virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
+ //virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
virtual bool onDomainBoundary() const { return _onDomainBoundary; }
virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
@@ -326,8 +296,8 @@ class Cell
virtual std::list< std::pair<int, Cell*> > getCells() { std::list< std::pair<int, Cell*> >cells; cells.push_back( std::make_pair(1, this)); return cells; }
virtual int getNumCells() {return 1;}
- bool operator==(const Cell& c2){
-
+
+ bool operator==(const Cell& c2){
if(this->getNumVertices() != c2.getNumVertices()){
return false;
}
@@ -337,10 +307,6 @@ class Cell
}
}
return true;
-
- //return (this->getTag() == c2.getTag());
-
-
}
bool operator<(const Cell& c2) const {
@@ -353,10 +319,9 @@ class Cell
else if (this->getSortedVertex(i) > c2.getSortedVertex(i)) return false;
}
return false;
-
- //return (this->getTag() < c2.getTag());
}
+ // checks whether lower dimensional simplex tau is on the boundary of this cell
virtual int kappa(Cell* tau) const;
};
@@ -369,10 +334,6 @@ class Simplex : public Cell
public:
Simplex() : Cell() {}
~Simplex(){}
-
- // kappa for simplices
- // checks whether lower dimensional simplex tau is on the boundary of this cell
- //virtual int kappa(Cell* tau) const;
};
@@ -409,6 +370,7 @@ class ZeroSimplex : public Simplex, public MPoint
class OneSimplex : public Simplex, public MLine
{
private:
+ // sorted list of vertices
int _vs[2];
public:
@@ -453,7 +415,7 @@ class OneSimplex : public Simplex, public MLine
class TwoSimplex : public Simplex, public MTriangle
{
private:
-
+ // sorted list of vertices
int _vs[3];
public:
@@ -497,7 +459,7 @@ class TwoSimplex : public Simplex, public MTriangle
class CQuadrangle : public Cell, public MQuadrangle
{
private:
-
+ // sorted list of vertices
int _vs[4];
public:
@@ -545,6 +507,7 @@ class CQuadrangle : public Cell, public MQuadrangle
class ThreeSimplex : public Simplex, public MTetrahedron
{
private:
+ // sorted list of vertices
int _vs[4];
public:
@@ -586,6 +549,7 @@ class ThreeSimplex : public Simplex, public MTetrahedron
class CHexahedron : public Cell, public MHexahedron
{
private:
+ // sorted list of vertices
int _vs[8];
public:
@@ -633,6 +597,7 @@ class CHexahedron : public Cell, public MHexahedron
class CPrism : public Cell, public MPrism
{
private:
+ // sorted list of vertices
int _vs[6];
public:
@@ -678,6 +643,7 @@ class CPrism : public Cell, public MPrism
class CPyramid : public Cell, public MPyramid
{
private:
+ // sorted list of vertices
int _vs[5];
public:
@@ -716,57 +682,6 @@ class CPyramid : public Cell, public MPyramid
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _inSubdomain); }
};
-/*
-// Ordering for the cells.
-class Less_Cell{
- public:
- bool operator()(const Cell* c1, const Cell* c2) const {
-
- // cells with fever vertices first
-
- if(c1->getNumVertices() != c2->getNumVertices()){
- return (c1->getNumVertices() < c2->getNumVertices());
- }
-
-
-
- // "natural number" -like ordering (the number of a vertice corresponds a digit)
-
- for(int i=0; i < c1->getNumVertices();i++){
- if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
- else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
- }
-
- return false;
-
-
- //return (c1->getTag() < c2->getTag());
-
- }
-};
-*/
-
-
-class Equal_Cell{
- public:
- bool operator ()(const Cell* c1, const Cell* c2){
-
-
- if(c1->getNumVertices() != c2->getNumVertices()){
- return false;
- }
- for(int i=0; i < c1->getNumVertices();i++){
- if(c1->getSortedVertex(i) != c2->getSortedVertex(i)){
- return false;
- }
- }
- return true;
-
- //return (c1->getTag() == c2->getTag());
-
- }
-};
-
// Ordering for the finite element mesh vertices.
class Less_MVertex{
@@ -780,8 +695,11 @@ class Less_MVertex{
class CombinedCell : public Cell{
private:
+ // vertices
std::vector<MVertex*> _v;
+ // sorted list of all vertices
std::vector<int> _vs;
+ // list of cells this cell is a combination of
std::list< std::pair<int, Cell*> > _cells;
int _num;
@@ -813,9 +731,11 @@ class CombinedCell : public Cell{
if(!this->hasVertex(c2->getVertex(i)->getNum())) _v.push_back(c2->getVertex(i));
}
+ // sorted vertices
for(unsigned int i = 0; i < _v.size(); i++) _vs.push_back(_v.at(i)->getNum());
std::sort(_vs.begin(), _vs.end());
+ // cells
_cells = c1->getCells();
std::list< std::pair<int, Cell*> > c2Cells = c2->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = c2Cells.begin(); it != c2Cells.end(); it++){
@@ -823,12 +743,10 @@ class CombinedCell : public Cell{
_cells.push_back(*it);
}
-
- //_boundary = c1->getOrientedBoundary();
+ // boundary cells
std::map< Cell*, int, Less_Cell > c1Boundary = c1->getOrientedBoundary();
std::map< Cell*, int, Less_Cell > c2Boundary = c2->getOrientedBoundary();
-
for(std::map<Cell*, int, Less_Cell>::iterator it = c1Boundary.begin(); it != c1Boundary.end(); it++){
Cell* cell = (*it).first;
int ori = (*it).second;
@@ -839,18 +757,11 @@ class CombinedCell : public Cell{
Cell* cell = (*it).first;
if(!orMatch) (*it).second = -1*(*it).second;
int ori = (*it).second;
- cell->removeCoboundaryCell(c2);
- //if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
+ cell->removeCoboundaryCell(c2);
if(co){
std::map<Cell*, int, Less_Cell >::iterator it2 = c1Boundary.find(cell);
bool old = false;
if(it2 != c1Boundary.end()) old = true;
- /*
- for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Boundary.begin(); it2 != c1Boundary.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cell) old = true;
- }
- */
if(!old){ if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this); }
}
else{
@@ -858,7 +769,7 @@ class CombinedCell : public Cell{
}
}
- //_coboundary = c1->getOrientedCoboundary();
+ // coboundary cells
std::map<Cell*, int, Less_Cell > c1Coboundary = c1->getOrientedCoboundary();
std::map<Cell*, int, Less_Cell > c2Coboundary = c2->getOrientedCoboundary();
@@ -872,18 +783,11 @@ class CombinedCell : public Cell{
Cell* cell = (*it).first;
if(!orMatch) (*it).second = -1*(*it).second;
int ori = (*it).second;
- cell->removeBoundaryCell(c2);
- //if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
+ cell->removeBoundaryCell(c2);
if(!co){
std::map<Cell*, int, Less_Cell >::iterator it2 = c1Coboundary.find(cell);
bool old = false;
if(it2 != c1Coboundary.end()) old = true;
- /*
- for(std::list< std::pair<int, Cell* > >::iterator it2 = c1Coboundary.begin(); it2 != c1Coboundary.end(); it2++){
- Cell* cell2 = (*it2).second;
- if(*cell2 == *cell) old = true;
- }
- */
if(!old) { if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this); }
}
else {
@@ -947,6 +851,7 @@ class CellComplex
std::vector<GEntity*> _boundary;
+ // mesh vertices in this domain
std::set<MVertex*, Less_MVertex> _domainVertices;
// sorted containers of unique cells in this cell complex
@@ -954,13 +859,13 @@ class CellComplex
std::set<Cell*, Less_Cell> _cells[4];
// storage for cell pointers to delete them
- //std::set<Cell*, Less_Cell> _cells2[4];
std::list<Cell*> _trash;
+
std::vector< std::set<Cell*, Less_Cell> > _store;
std::set<Cell*, Less_Cell> _ecells;
-
+ // original cells of this cell complex
std::set<Cell*, Less_Cell> _ocells[4];
// Betti numbers of this cell complex (ranks of homology groups)
@@ -979,9 +884,10 @@ class CellComplex
// remove cell from the queue set
void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
- // insert cells into this cell complex
- //virtual void insert_cells(bool subdomain, bool boundary);
+ // for constructor
void insert_cells(bool subdomain, bool boundary);
+ void find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain);
+ void find_vertices(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain);
// remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
@@ -997,17 +903,8 @@ class CellComplex
public:
- /*
- CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
- _domain = domain;
- _subdomain = subdomain;
- for(citer cit = cells.begin(); cit != cells.end(); cit++){
- Cell* cell = *cit;
- _cells[cell->getDim()].insert(cell);
- }
- }
-
+ /*
CellComplex(CellComplex* cellComplex){
_domain = cellComplex->getDomain();
@@ -1050,6 +947,7 @@ class CellComplex
_trash.clear();
}
+ // restore this cell complex to its original state
void restoreComplex(){
for(int i = 0; i < 4; i++){
_betti[i] = 0;
@@ -1109,18 +1007,14 @@ class CellComplex
int reduceComplex();
int coreduceComplex();
-
- // add every volume, face and edge its missing boundary cells
- // void repairComplex(int i=3);
- // change non-subdomain cells to be in subdomain, subdomain cells to not to be in subdomain
- void swapSubdomain();
+ // remove cells in subdomain from this cell complex
void removeSubdomain();
// print the vertices of cells of certain dimension
void printComplex(int dim);
- // write this cell complex in 2.0 MSH ASCII format
+ // write this cell complex in 2.1 MSH ASCII format
int writeComplexMSH(const std::string &name);
// Cell combining for reduction and coreduction
@@ -1142,16 +1036,6 @@ class CellComplex
int getNumOmitted() { return _store.size(); };
std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
-
- void setImmuneCell(int num){
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- if(cell->getNum() == num) cell->setImmune(true);
- }
- }
-
- }
// change roles of boundaries and coboundaries of the cells in this cell complex
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index f25ec4c13d076a19e3292d6e30230fc5380b6663..85929ec46949a3f8c05391d5e6a67d2a11ff98a0 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "ChainComplex.h"
#include "OS.h"
@@ -567,10 +567,15 @@ bool Chain::straightenChain( std::pair<Cell*, int> cell ){
Cell* c1 = cell.first;
int c1c = cell.second;
+ int c1o = 0;
if(c1c == 0 || c1->getImmune()) return false;
+ if(cell.second == 0 || cell.first->getImmune()) return false;
+ /*
+ Cell* c1 = NULL;
+ int c1c = 0;
int c1o = 0;
-
+ */
Cell* c2 = NULL;
int c2c = 0;
int c2o = 0;
@@ -581,16 +586,29 @@ bool Chain::straightenChain( std::pair<Cell*, int> cell ){
Cell* b = NULL;
- std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
+ std::map<Cell*, int, Less_Cell> c1Cbd = cell.first->getOrgCbd();
for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
Cell* c1CbdCell = (*cit).first;
c1o = (*cit).second;
/*
std::map<Cell*, int, Less_Cell> cells = getBdCellsInChain(c1CbdCell);
- if((getDim() == 1 && cells.size() != 2) || (getDim() == 2 && cells.size() != 3) ) break;
- else if( getDim() == 1 && cells.size() == 2){
+ //if((getDim() == 1 && cells.size() != 2) || (getDim() == 2 && cells.size() != 3) ) continue;
+ if( getDim() == 1 && cells.size() == 2){
+ citer cit2 = cells.end();
+ c1 = (*cit2).first;
+ //c1->printCell();
+
+ c1c = getCoeff(c1);
+ c1o = (*cit2).second;
+
+ cit2 = cells.begin();
+ c2 = (*cit2).first;
+ c2c = getCoeff(c2);
+ c2o = (*cit2).second;
+ b = c1CbdCell;
+ printf("kakak4 \n");
}
else if( getDim() == 2 && cells.size() == 3){
@@ -619,6 +637,7 @@ bool Chain::straightenChain( std::pair<Cell*, int> cell ){
}
if (b != NULL) break;
+
}
if(c1->getDim() == 1 &&
@@ -798,86 +817,11 @@ void Chain::smoothenChain(){
if (useless > 5) break;
}
-
- /*
- const int MAXROUNDS = 2;
- bool smoothened = true;
- int round = 0;
- while(smoothened){
- round++;
- smoothened = false;
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(straightenChain(*cit)) { smoothened = true; }
- if(removeBoundary(*cit)) { smoothened = true; }
- }
- if(round > MAXROUNDS) smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ straightenChain(*cit);
}
- eraseNullCells();
- if(getDim() == 2){
- Chain* bestChain = new Chain(this);
- int before = getSize();
- deImmuneCells();
- srand ( time(NULL) );
- int n = 0;
- while( n < 10){
- double t = 1;
- while(t>0){
- double tt = 1;
-
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- int random = rand() % 100 + 1;
- double r = random*t;
- //printf("random: %d, t: %d \n", random, t);
- if(r > 80) {
- bendChain(*cit);
- //printf("random: %d, t: %g, random*t: %g.\n", random, t, r);
- //cit = _cells.begin();
- //tt = tt - 0.1;
- //if(tt <= 0) tt = 0;
- }
- }
- eraseNullCells();
-
- smoothened = true;
- round = 0;
- while(smoothened){
- round++;
- smoothened = false;
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(straightenChain(*cit)){smoothened = true;}
- if(removeBoundary(*cit)) { smoothened = true; }
- }
- if(round > MAXROUNDS) smoothened = false;
- }
- eraseNullCells();
- if(this->getSize() < bestChain->getSize()) bestChain = this;
- t = t - 0.1;
- }
- deImmuneCells();
- n=n+1;
- }
-
-
- deImmuneCells();
- smoothened = true;
- round = 0;
- while(smoothened){
- round++;
- smoothened = false;
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(straightenChain(*cit)){smoothened = true;}
- if(removeBoundary(*cit)) { smoothened = true; }
- }
- if(round > MAXROUNDS) smoothened = false;
- }
-
-
- eraseNullCells();
- if(this->getSize() < bestChain->getSize()) bestChain = this;
- printf("%d-chain simulated annealing removed %d cells. \n", getDim(), before - getSize());
- }
- */
+
double t2 = Cpu();
printf("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
return;
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 6a7dc870913e89718a69027f4c02138cd73686c2..7739b36c37bc9a2f81f4a1f4ebeea8bb7a73da7c 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -115,8 +115,11 @@ 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
std::vector<int> getCoeffVector(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; }
int printMatrix(gmp_matrix* matrix){
@@ -139,6 +142,7 @@ class Chain{
// cell complex this chain belongs to
CellComplex* _cellComplex;
+ // torsion coefficient
int _torsion;
int _dim;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 2c8e69a5bc2831d85ba41522e82105c7335ad472..8175fc95f22d1fb6bd95325ae78417ac60a833f2 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -42,7 +42,7 @@ class Homology
void computeBettiNumbers();
- void swapSubdomain() { _cellComplex->swapSubdomain(); }
+ //void swapSubdomain() { _cellComplex->swapSubdomain(); }
// Restore the cell complex to its original state before cell reductions
void restoreHomology() { _cellComplex->restoreComplex(); }
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index eb5734d831ffa9472f8a667389e3791730251c75..4f6193f5a72b8341555fb44c0e49da90fc04bb08 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include <stdlib.h>
#include "Gmsh.h"
@@ -23,8 +23,6 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "ComputeGenerators", NULL, 1.},
{GMSH_FULLRC, "ComputeDualGenerators", NULL, 0.},
{GMSH_FULLRC, "ComputeBettiNumbers", NULL, 0.},
- //{GMSH_FULLRC, "OmitDimensions", NULL, 1.},
- //{GMSH_FULLRC, "CombineCells", NULL, 1.},
};
StringXString HomologyComputationOptions_String[] = {
@@ -86,15 +84,11 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
int gens = (int)HomologyComputationOptions_Number[4].def;
int cuts = (int)HomologyComputationOptions_Number[5].def;
int betti = (int)HomologyComputationOptions_Number[6].def;
- //int omit = (int)HomologyComputationOptions_Number[6].def;
- //int combine = (int)HomologyComputationOptions_Number[7].def;
GModel *m = GModel::current();
Homology* homology = new Homology(m, domain, subdomain);
- //if(combine == 0) homology->setCombine(false);
- //homology->setOmit(1);
//if(swap == 1) homology->swapSubdomain();
diff --git a/Plugin/HomologyComputation.h b/Plugin/HomologyComputation.h
index 26f510307372227019cc8c4669fc0120ada92ad7..a792885670f32488adcdbe78d0b9c1ce714645e3 100644
--- a/Plugin/HomologyComputation.h
+++ b/Plugin/HomologyComputation.h
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#ifndef _HOMOLOGY_COMPUTATION_H_
#define _HOMOLOGY_COMPUTATION_H_