diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index ecf672d2f9e7737156915d9983cf4e5e6290954c..cee51bcecc45f1bbd6d78d0a21ec2d79ca314de7 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -658,6 +658,7 @@ int CellComplex::combine(int dim){
replaceCells(c1, c2, newCell, (or1 != or2));
cit = firstCell(dim);
+ //cit++;
count++;
}
removeCellQset(s, Qset);
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 876f9ae76aff53451da2f5872abf42ee7a2a694b..769e610b4d50eaeb334f063c1f2670dcd1569638 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -275,16 +275,14 @@ void ChainComplex::Quotient(int dim){
mpz_t elem;
mpz_init(elem);
-
+
for(int i = 1; i <= cols; i++){
gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
if(mpz_cmp_si(elem,0) == 0){
destroy_gmp_normal_form(normalForm);
return;
}
- if(mpz_cmp_si(elem,1) > 0){
- _torsion[dim].push_back(mpz_get_si(elem));
- }
+ if(mpz_cmp_si(elem,1) > 0) _torsion[dim].push_back(mpz_get_si(elem));
}
int rank = cols - _torsion[dim].size();
@@ -437,7 +435,15 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){
}
-Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name){
+int ChainComplex::getTorsion(int dim, int chainNumber){
+ if(dim < 0 || dim > 4) return 0;
+ if(_Hbasis[dim] == NULL || gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0;
+ if(_torsion[dim].empty() || _torsion[dim].size() < chainNumber) return 1;
+ else return _torsion[dim].at(chainNumber-1);
+
+}
+
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name, int torsion){
int i = 0;
for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
@@ -448,9 +454,9 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
}
}
-
_name = name;
_cellComplex = cellComplex;
+ _torsion = torsion;
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 558cc50cd5f41253c8b86c5d18bc6ee27c5ee136..f5569a52d502d473d622d59d1c32b4c513eb8250 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -53,19 +53,19 @@ class ChainComplex{
CellComplex* _cellComplex;
// set the matrices
- virtual void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
- virtual void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
- virtual void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _codH[dim] = matrix;}
- virtual void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
- virtual void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
- virtual void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
+ void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
+ void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
+ void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _codH[dim] = matrix;}
+ void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
+ void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
+ void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
public:
ChainComplex(CellComplex* cellComplex);
-
+
ChainComplex(){
for(int i = 0; i < 5; i++){
_HMatrix[i] = create_gmp_matrix_zero(1,1);
@@ -77,43 +77,44 @@ class ChainComplex{
}
_dim = 0;
}
- virtual ~ChainComplex(){}
+ ~ChainComplex(){}
- virtual int getDim() { return _dim; }
+ int getDim() { return _dim; }
// get the boundary operator matrix dim->dim-1
- virtual gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
- virtual gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
- virtual gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
+ gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
+ gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
+ gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
+ gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
+ gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
+ gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
// Compute basis for kernel and codomain of boundary operator matrix of dimension dim (ie. ker(h_dim) and cod(h_dim) )
- virtual void KerCod(int dim);
+ void KerCod(int dim);
// Compute matrix representation J for inclusion relation from dim-cells who are boundary of dim+1-cells
// to cycles of dim-cells (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
- virtual void Inclusion(int lowDim, int highDim);
+ void Inclusion(int lowDim, int highDim);
// Compute quotient problem for given inclusion relation j to find representatives of homology groups
// and possible torsion coeffcients
- virtual void Quotient(int dim);
+ void Quotient(int dim);
// transpose the boundary operator matrices, these are boundary operator matrices for the dual mesh
- virtual void transposeHMatrices() { for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
- virtual void transposeHMatrix(int dim) { if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
+ void transposeHMatrices() { for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
+ void transposeHMatrix(int dim) { if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
// Compute bases for the homology groups of this chain complex
- virtual void computeHomology(bool dual=false);
+ void computeHomology(bool dual=false);
- virtual std::vector<int> getCoeffVector(int dim, int chainNumber);
- virtual int getBasisSize(int dim) { if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
+ std::vector<int> getCoeffVector(int dim, int chainNumber);
+ int getTorsion(int dim, int chainNumber);
+ int getBasisSize(int dim) { if(dim > -1 && dim < 5) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
- virtual int printMatrix(gmp_matrix* matrix){
+ int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
return gmp_matrix_printf(matrix); }
-
+
// debugging aid
- virtual void matrixTest();
+ void matrixTest();
};
// A class representing a chain.
@@ -128,20 +129,24 @@ class Chain{
// cell complex this chain belongs to
CellComplex* _cellComplex;
+ int _torsion;
+
public:
Chain(){}
- Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain");
+ Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain", int torsion=0);
~Chain(){}
// get i:th cell of this chain
- virtual Cell* getCell(int i) { return _cells.at(i).first; }
+ Cell* getCell(int i) { return _cells.at(i).first; }
// get coeffcient of i:th cell of this chain
- virtual int getCoeff(int i) { return _cells.at(i).second; }
+ int getCoeff(int i) { return _cells.at(i).second; }
+
+ int getTorsion() {return _torsion;}
// number of cells in this chain
- virtual int getSize() { return _cells.size(); }
+ int getSize() { return _cells.size(); }
- virtual int getNumCells() {
+ int getNumCells() {
int count = 0;
for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
Cell* cell = (*it).first;
@@ -152,13 +157,13 @@ class Chain{
// get/set chain name
- virtual std::string getName() { return _name; }
- virtual void setName(std::string name) { _name=name; }
+ std::string getName() { return _name; }
+ void setName(std::string name) { _name=name; }
// append this chain to a 2.0 MSH ASCII file as $ElementData
- virtual int writeChainMSH(const std::string &name);
+ int writeChainMSH(const std::string &name);
- virtual void getData(std::map<int, std::vector<double> >& data);
+ void getData(std::map<int, std::vector<double> >& data);
};
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 7174b65e0f2540f3bc5fa5e847c794931a45e33f..8a7b7b154353a0fe51ac53c624ab41c4775ba035 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -8,6 +8,7 @@
#include "Homology.h"
#include "ChainComplex.h"
+#include "OS.h"
#if defined(HAVE_KBIPACK)
Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
@@ -15,7 +16,7 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
_model = model;
Msg::Info("Creating a Cell Complex...");
-
+ double t1 = Cpu();
std::map<int, std::vector<GEntity*> > groups[4];
model->getPhysicalGroups(groups);
@@ -60,8 +61,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
Msg::Error("Check the domain & the mesh.");
return;
}
-
- Msg::Info("Cell Complex complete.");
+ double t2 = Cpu();
+ 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));
@@ -71,20 +72,24 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
+ double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(true);
_cellComplex->combine(3);
_cellComplex->combine(2);
_cellComplex->combine(1);
- Msg::Info("Cell Complex reduction complete.");
+ double t2 = Cpu();
+ 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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
+ t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex);
chains->computeHomology();
- Msg::Info("Homology Computation complete.");
+ t2 = Cpu();
+ Msg::Info("Homology Computation complete (%g s).", t2 - t1);
int HRank[4];
for(int j = 0; j < 4; j++){
@@ -97,9 +102,12 @@ void Homology::findGenerators(std::string fileName){
convert(j, dimension);
std::string name = dimension + "D Generator " + generator;
- Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name);
+ Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name, chains->getTorsion(j,i));
chain->writeChainMSH(fileName);
- if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
+ if(chain->getSize() != 0) {
+ HRank[j] = HRank[j] + 1;
+ if(chain->getTorsion() != 1) Msg::Warning("%dD Generator %d has torsion coefficient %d!", j, i, chain->getTorsion());
+ }
delete chain;
}
}
@@ -124,27 +132,31 @@ void Homology::findThickCuts(std::string fileName){
_cellComplex->removeSubdomain();
Msg::Info("Reducing Cell Complex...");
+ double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(true);
_cellComplex->cocombine(0);
_cellComplex->cocombine(1);
_cellComplex->cocombine(2);
- Msg::Info("Cell Complex reduction complete.");
+ double t2 = Cpu();
+ 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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
+ t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex);
chains->transposeHMatrices();
chains->computeHomology(true);
- Msg::Info("Homology Computation complete.");
+ t2 = Cpu();
+ 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--){
- HRank[j] = 0;
for(int i = 1; i <= chains->getBasisSize(j); i++){
std::string generator;
@@ -153,9 +165,12 @@ void Homology::findThickCuts(std::string fileName){
convert(dim-j, dimension);
std::string name = dimension + "D Thick cut " + generator;
- Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name);
+ Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name, chains->getTorsion(j,i));
chain->writeChainMSH(fileName);
- if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
+ if(chain->getSize() != 0){
+ HRank[dim-j] = HRank[dim-j] + 1;
+ if(chain->getTorsion() != 1) Msg::Warning("%dD Thick cut %d has torsion coefficient %d!", j, i, chain->getTorsion());
+ }
delete chain;
}