Added Chain class in order to visualize the homology generator chains.

Modified utils/misc/celldriver.cpp example to demonstrate this.

Geo/CellComplex.cpp

@@ -16,6 +16,15 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
   insertCells(true);
   insertCells(false);
 
+  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++;
+    }
+  }
+  
 }
 void CellComplex::insertCells(bool subdomain){  
   
@@ -493,13 +502,18 @@ int CellComplex::writeComplexMSH(const std::string &name){
   
     
   FILE *fp = fopen(name.c_str(), "w");
+  
   if(!fp){
     Msg::Error("Unable to open file '%s'", name.c_str());
     printf("Unable to open file.");
     return 0;
   }
   
-  fprintf(fp, "$NOD\n");
+  
+  
+  fprintf(fp, "$MeshFormat\n2.0 0 8\n$EndMeshFormat\n");
+  
+  fprintf(fp, "$Nodes\n");
   
   std::set<MVertex*, Less_MVertex> domainVertices;
   getDomainVertices(domainVertices, true);
@@ -513,38 +527,34 @@ int CellComplex::writeComplexMSH(const std::string &name){
   }
   
       
-  fprintf(fp, "$ENDNOD\n");
-  fprintf(fp, "$ELM\n");
+  fprintf(fp, "$EndNodes\n");
+  fprintf(fp, "$Elements\n");
 
   fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size());
   
-  int index = 1;
-  
   for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
     Cell* vertex = *cit;
-    fprintf(fp, "%d %d %d %d %d %d\n", index, 15, 0, 1, 1, vertex->getVertex(0));
-    index++;
+    fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, 0, vertex->getVertex(0));
   }
   
   
   for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
     Cell* edge = *cit;
-    fprintf(fp, "%d %d %d %d %d %d %d\n", index, 1, 0, 1, 2, edge->getVertex(0), edge->getVertex(1));
-    index++;
+    fprintf(fp, "%d %d %d %d %d %d %d %d\n", edge->getTag(), 1, 3, 0, 0, 0, edge->getVertex(0), edge->getVertex(1));
   }
   
   for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
     Cell* face = *cit;
-    fprintf(fp, "%d %d %d %d %d %d %d %d\n", index, 2, 0, 1, 3, face->getVertex(0), face->getVertex(1), face->getVertex(2));
-    index++;
+    fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", face->getTag(), 2, 3, 0, 0, 0, face->getVertex(0), face->getVertex(1), face->getVertex(2));
   }
   for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
     Cell* volume = *cit;
-    fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", index, 4, 0, 1, 4, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
-    index++;
+    fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", volume->getTag(), 4, 3, 0, 0, 0, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
   }
     
-  fprintf(fp, "$ENDELM\n");
+  fprintf(fp, "$EndElements\n");
+  
+  fclose(fp);
   
   return 1;
 }

Geo/CellComplex.h

@@ -40,12 +40,15 @@ class Cell
    // used in relative homology computation
    bool _inSubdomain;
    
+   int _tag;
+   
   public:
    Cell(){}
    virtual ~Cell(){}
    
    virtual int getDim() const { return _dim; };
-   //virtual int getTag() const { return _tag; };
+   virtual int getTag() const { return _tag; };
+   virtual void setTag(int tag) { _tag = tag; };
    
    // get the number of vertices this cell has
    virtual int getNumVertices() const = 0; //{return _vertices.size();}
@@ -112,7 +115,7 @@ class Simplex : public Cell
 { 
    
  public:
-   Simplex(){
+   Simplex() : Cell() {
      _cbdSize = 1000; // big number
      _bdSize = 1000;
    }
@@ -345,10 +348,6 @@ class CellComplex
    // remove cell from the queue set
    virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
       
-   // get all the finite element mesh vertices in the domain of this cell complex
-   virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, 
-                                  bool subdomain);
-   
    // insert cells into this cell complex
    virtual void insertCells(bool subdomain);
    
@@ -369,6 +368,9 @@ class CellComplex
    
    virtual std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
    
+   // get all the finite element mesh vertices in the domain of this cell complex
+   virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain);
+   
    // iterators to the first and last cells of certain dimension
    virtual citer firstCell(int dim) {return _cells[dim].begin(); }
    virtual citer lastCell(int dim) {return _cells[dim].end(); }
@@ -377,7 +379,6 @@ class CellComplex
    virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices);
    virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
    
-   
    // kappa for two cells of this cell complex
    // implementation will vary depending on cell type
    virtual int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
@@ -424,7 +425,7 @@ class CellComplex
    // print the vertices of cells of certain dimension
    virtual void printComplex(int dim, bool subcomplex);
    
-   // write this cell complex in legacy .msh format
+   // write this cell complex in .msh format
    virtual int writeComplexMSH(const std::string &name); 
    
    

Geo/ChainComplex.cpp

@@ -12,6 +12,11 @@
 ChainComplex::ChainComplex(CellComplex* cellComplex){
   
   _HMatrix[0] = NULL;
+  _kerH[0] = NULL;
+  _codH[0] = NULL;
+  _JMatrix[0] = NULL;
+  _QMatrix[0] = NULL;
+  _Hbasis[0] = NULL;
   
   for(int dim = 1; dim < 4; dim++){
     unsigned int cols = cellComplex->getSize(dim);
@@ -263,8 +268,10 @@ void ChainComplex::computeHomology(){
       else{
         _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1, 
                                      gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+        
         gmp_matrix_right_mult(_Hbasis[i], getQMatrix(i));
       } 
+      
     } 
     
     
@@ -303,5 +310,84 @@ void ChainComplex::matrixTest(){
   return; 
 }
 
+std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){
+  
+  std::vector<int> coeffVector;
+  
+  if(dim < 0 || dim > 3) return coeffVector;
+  if(_Hbasis[dim] == NULL || gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return coeffVector;
+  
+  int rows = gmp_matrix_rows(_Hbasis[dim]);
+  
+  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]);
+    elemli = mpz_get_si(elem);
+    elemi = elemli;
+    coeffVector.push_back(elemi);
+    //printf("coeff: %d \n", coeffVector.at(i-1));
+  }
+  
+  mpz_clear(elem);
+  return coeffVector;
+  
+}
 
 #endif
+
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name){
+  
+  int i = 0;
+  for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
+    Cell* cell = *cit;
+    if(!cell->inSubdomain() && coeffs.size() > i){
+      if(coeffs.at(i) != 0) _cells.push_back( std::make_pair(cell, coeffs.at(i)) );
+      i++;
+    }
+  }
+  
+  _name = name;
+  _cellComplex = cellComplex;
+  
+}
+
+int Chain::writeChainMSH(const std::string &name){
+
+  //_cellComplex->writeComplexMSH(name);
+  
+  FILE *fp = fopen(name.c_str(), "a");
+  if(!fp){
+    Msg::Error("Unable to open file '%s'", name.c_str());
+        printf("Unable to open file.");
+        return 0;
+  }
+
+  fprintf(fp, "\n$ElementData\n");
+  
+  fprintf(fp, "1 \n");
+  fprintf(fp, "\"%s\" \n", getName().c_str());
+  fprintf(fp, "1 \n");
+  fprintf(fp, "0.0 \n");
+  fprintf(fp, "4 \n");
+  fprintf(fp, "0 \n");
+  fprintf(fp, "1 \n");
+  fprintf(fp, "%d \n", getSize());
+  fprintf(fp, "0 \n");
+  
+  for(int i = 0; i < _cells.size(); i++){
+    
+    fprintf(fp, "%d %d \n", getCell(i)->getTag(), getCoeff(i));
+    
+  }
+  
+  fprintf(fp, "$EndElementData\n");
+  
+  fclose(fp);
+  
+  return 1;
+  
+}

Geo/ChainComplex.h

@@ -49,6 +49,8 @@ class ChainComplex{
    // bases for homology groups
    gmp_matrix* _Hbasis[4];
    
+   
+   
   public:
    
    ChainComplex(CellComplex* cellComplex);
@@ -85,6 +87,9 @@ class ChainComplex{
    // Compute bases for the homology groups of this chain complex 
    virtual void computeHomology();
    
+   virtual std::vector<int> getCoeffVector(int dim, int chainNumber);
+   virtual int getBasisSize(int dim) {  if(dim > -1 || dim < 4) return gmp_matrix_cols(_Hbasis[dim]); else return 0; } 
+   
    virtual int printMatrix(gmp_matrix* matrix){ 
      printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix)); 
      return gmp_matrix_printf(matrix); } 
@@ -93,7 +98,40 @@ class ChainComplex{
    virtual void matrixTest();
 };
 
-
 #endif
 
+// A class representing a chain.
+// Used to visualize generators of the homology groups.
+class Chain{
+   
+  private:
+   // cells and their coefficients in this chain
+   std::vector< std::pair <Cell*, int> > _cells;
+   // name of the chain (optional)
+   std::string _name;
+   // cell complex this chain belongs to
+   CellComplex* _cellComplex;
+   
+  public:
+   Chain(){}
+   Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain");
+   ~Chain(){}
+   
+   // get i:th cell of this chain
+   virtual 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; }
+   
+   // number of cells in this chain 
+   virtual int getSize() { return _cells.size(); }
+   
+   // get/set chain name
+   virtual std::string getName() { return _name; }
+   virtual void setName(std::string name) { _name=name; }
+
+   // append this chain to a 2.0 .msh file as $ElementData
+   virtual int writeChainMSH(const std::string &name);
+   
+};
+
 #endif

utils/misc/celldriver.cpp

@@ -5,15 +5,26 @@
 //
 // Then compile this driver with "g++ celldriver.cpp -lGmsh -llapack -lblas -lgmp"
 //
-// This program creates .msh files reduced_complex.msh and coreduced_complex.msh
-// of an two port transformer model to represent its relative homology groups.
+//
+// This program creates .msh file chains.msh which represents the generators of 
+// an two port transformer model's homology groups.
 
 
 #include <stdio.h>
+#include <sstream>
 #include <gmsh/Gmsh.h>
 #include <gmsh/GModel.h>
 #include <gmsh/MElement.h>
 #include <gmsh/CellComplex.h>
+#include <gmsh/ChainComplex.h>
+
+template <class TTypeA, class TTypeB>
+bool convert(const TTypeA& input, TTypeB& output ){
+  std::stringstream stream;
+  stream << input;
+  stream >> output;
+  return stream.good();
+}
 
 
 int main(int argc, char **argv)
@@ -46,16 +57,45 @@ int main(int argc, char **argv)
   s= *(++sub);
   subdomain.push_back(s);
   
-  CellComplex complex = CellComplex(domain, subdomain);
+  // create a cell complex
+  CellComplex* complex = new CellComplex(domain, subdomain);
 
   printf("Cell complex of this model has: %d volumes, %d faces, %d edges and %d vertices\n",
-         complex.getSize(3), complex.getSize(2), complex.getSize(1), complex.getSize(0));  
+         complex->getSize(3), complex->getSize(2), complex->getSize(1), complex->getSize(0));  
+  
+  // reduce the complex in order to ease homology computation
+  complex->reduceComplex();
+  
+  // create a chain complex of a cell complex (construct boundary operator matrices)
+  ChainComplex* chains = new ChainComplex(complex);
+  // compute the homology using the boundary operator matrices
+  chains->computeHomology();
+  
+  // write the reduced cell complex to a .msh file
+  complex->writeComplexMSH("chains.msh");
+  
+  // Append the homology generators to the .msh file as post-processing views. 
+  // To visualize the result, open chains.msh with Gmsh GUI and deselect all mesh
+  // entities from Tools->Options->Mesh->Visibility.
+  for(int j = 0; j < 4; j++){
+    for(int i = 1; i <= chains->getBasisSize(j); i++){
+    
+      std::string generator;
+      std::string dimension;
+      convert(i, generator);
+      convert(j, dimension);
+      std::string name = dimension + "D Generator " + generator;
+      
+      Chain* chain = new Chain(complex->getCells(j), chains->getCoeffVector(j,i), complex, name);
+      chain->writeChainMSH("chains.msh");
+      delete chain;
+      
+    }
+  }
   
-  complex.reduceComplex();
-  complex.writeComplexMSH("reduced_complex.msh");
-  complex.coreduceComplex();
-  complex.writeComplexMSH("coreduced_complex.msh");
     
+  delete chains;
+  delete complex;  
   delete m;
   GmshFinalize();