Bugfixes. Added notion about possible torsion of generators of homology
groups.

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);

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;
   
 }
 

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);
    
 };
 

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;
             
     }