Bugfixes in CellComplex and in the examples.

Begun writing ChainComplex class. Will have own file some day.

Geo/CellComplex.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, 16.3.2009.
+// Contributed by Matti Pellikka.
 
 #include "CellComplex.h"
 
@@ -18,7 +18,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
   insertCells(true);
   insertCells(false);
   
-  
 }
 void CellComplex::insertCells(bool subdomain){  
   
@@ -79,7 +78,7 @@ int Simplex::kappa(Cell* tau) const{
     value = value*-1;
   }
   
-  return 1;  
+  return value;  
 }
 std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices){
   if(dim == 0) return _cells[dim].find(new ZeroSimplex(vertices.at(0)));
@@ -289,7 +288,6 @@ int CellComplex::coreduction(int dim){
       bd_c = bd(cell);
       if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
         removeCell(cell);
-        //cit = firstCell(dim-1);
         removeCell(bd_c.at(0));
         count++;
         coreduced =true;
@@ -314,11 +312,9 @@ int CellComplex::reduction(int dim){
     reduced = false;
     for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
       Cell* cell = *cit;
-      //if(reductionMrozek(cell) !=0) cit = firstCell(dim-1);   
       cbd_c = cbd(cell);
       if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
         removeCell(cell);
-        //cit = firstCell(dim-1);
         removeCell(cbd_c.at(0));
         count++;
         reduced =true;
@@ -355,49 +351,117 @@ void CellComplex::coreduceComplex(){
   printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
          getSize(3), getSize(2), getSize(1), getSize(0));
 }
-
+/*
 void CellComplex::constructHMatrix(int dim){
   
-  // h: C_dim -> C_(dim-1)
+  // h_dim: C_dim -> C_(dim-1)
+
   
   if(dim > 3 || dim < 1){
     return;
   }
   
+  destroy_gmp_matrix(_HMatrix[dim]);
+  
   if( _cells[dim].size() == 0 ){
     _HMatrix[dim] = create_gmp_matrix_zero(1, 1);
+    //gmp_matrix_printf(_HMatrix[dim]);
     return;
   }
   unsigned int cols = _cells[dim].size(); 
   
   if( _cells[dim-1].size() == 0){ 
     _HMatrix[dim] = create_gmp_matrix_zero(1, cols);
+    //gmp_matrix_printf(_HMatrix[dim]);
     return;
   }
   unsigned int rows = _cells[dim-1].size();
   
   mpz_t elems[rows*cols];
-  mpz_array_init( elems[0], rows*cols, 512 );  
+  mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );  
     
   citer high = firstCell(dim);
   citer low = firstCell(dim-1);
   
-  //printf("laa %d %d %d \n", rows, cols, rows*cols);
+  /*
+  printf("laa %d %d %d \n", rows, cols, rows*cols);
   for(unsigned int i=0; i < rows*cols; i++){
-    //printf("%d, ", i);
+    printf(" %d, ", i);
     if(low == lastCell(dim-1)){
-      //printf("rowfull %d", i);
+      printf("rowfull %d ", i);
       high++;
       low = firstCell(dim-1);
     }
     mpz_set_ui(elems[i], kappa(*high, *low));
     low++;
   }
+  */
+  /*
+  unsigned int i = 0;
+  while(i < rows*cols){
+    while(low != lastCell(dim-1)){
+      mpz_set_si(elems[i], kappa(*high, *low));
+      //printf(" %d %d, ",i,  kappa(*high, *low));
+      i++;
+      low++;
+    }
+    low = firstCell(dim-1);
+    high++;
+  }
+  
+  _HMatrix[dim] = create_gmp_matrix(rows, cols,(const mpz_t *) elems);
    
-  _HMatrix[dim] = create_gmp_matrix(rows, cols, elems);
+  //gmp_matrix_printf(_HMatrix[dim]);
   
   return; 
 }
+*/
+std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
+  
+  // h_dim: C_dim -> C_(dim-1)
+  
+  std::vector<gmp_matrix*> HMatrix;
+  
+  HMatrix.push_back(create_gmp_matrix_zero(1, 1));
+  
+  for(int dim = 1; dim < 4; dim++){
+    if( _cells[dim].size() == 0 ){
+      HMatrix.push_back( create_gmp_matrix_zero(1, 1));
+      //gmp_matrix_printf(HMatrix[dim]);
+      break;
+    }
+    unsigned int cols = _cells[dim].size(); 
+    
+    if( _cells[dim-1].size() == 0){ 
+      HMatrix.push_back(create_gmp_matrix_zero(1, cols));
+      //gmp_matrix_printf(HMatrix[dim]);
+      break;
+    }
+    unsigned int rows = _cells[dim-1].size();
+    mpz_t elems[rows*cols];
+    mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );  
+    
+    citer high = firstCell(dim);
+    citer low = firstCell(dim-1);
+    
+
+    unsigned int i = 0;
+    while(i < rows*cols){
+      while(low != lastCell(dim-1)){
+        mpz_set_si(elems[i], kappa(*high, *low));
+        //printf(" %d %d, ",i,  kappa(*high, *low));
+        i++;
+        low++;
+      }
+      low = firstCell(dim-1);
+      high++;
+    }
+    HMatrix.push_back(create_gmp_matrix(rows, cols,(const mpz_t *) elems));
+    
+    //gmp_matrix_printf(_HMatrix[dim]);
+  }
+  return HMatrix; 
+}
 
 
 void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
@@ -521,4 +585,21 @@ void CellComplex::printComplex(int dim){
   }
 }
 
+gmp_matrix* ChainComplex::ker(gmp_matrix* HMatrix){
+  
+  // H = USV -> WH = US, W = inv(V)
+  gmp_normal_form* W_USform = create_gmp_Smith_normal_form(HMatrix, NOT_INVERTED, INVERTED);
+ 
+  
+  
+  return W_USform->canonical;
+  
+}
+  
+
+
+
+
 #endif
+
+

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, 16.3.2009.
+// Contributed by Matti Pellikka.
 
 #ifndef _CELLCOMPLEX_H_
 #define _CELLCOMPLEX_H_
@@ -155,6 +155,7 @@ class ZeroSimplex : public Simplex
    
    virtual ~ZeroSimplex(){}
    
+   virtual int getDim() const { return 0; }
    virtual int getNumVertices() const { return 1; }
    virtual int getVertex(int vertex) const {return _v; }
    virtual bool hasVertex(int vertex) const {return (_v == vertex); }
@@ -198,6 +199,7 @@ class OneSimplex : public Simplex
    
    virtual ~OneSimplex(){}
    
+   virtual int getDim() const { return 1; }
    virtual int getNumVertices() const { return 2; }
    virtual int getVertex(int vertex) const {return _v[vertex]; }
    virtual bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
@@ -236,6 +238,7 @@ class TwoSimplex : public Simplex
    
    virtual ~TwoSimplex(){}
    
+   virtual int getDim() const { return 2; }
    virtual int getNumVertices() const { return 3; }
    virtual int getVertex(int vertex) const {return _v[vertex]; }
    virtual bool hasVertex(int vertex) const {return 
@@ -276,6 +279,7 @@ class ThreeSimplex : public Simplex
    
    virtual ~ThreeSimplex(){}
    
+   virtual int getDim() const { return 3; }
    virtual int getNumVertices() const { return 4; }
    virtual int getVertex(int vertex) const {return _v[vertex]; }
    virtual bool hasVertex(int vertex) const {return 
@@ -294,18 +298,14 @@ class Less_Cell{
        return (c1->getNumVertices() < c2->getNumVertices());
      }
      for(int i=0; i < c1->getNumVertices();i++){
-       if(c1->getVertex(i) < c2->getVertex(i)){
-         return true;
-       }
-       else if (c1->getVertex(i) > c2->getVertex(i)){
-         return false;
-       }
+       if(c1->getVertex(i) < c2->getVertex(i)) return true;
+       else if (c1->getVertex(i) > c2->getVertex(i)) return false;
      }
      return false;
    }
 };
 
-// Ordering for the finite element mesh vertices
+// Ordering for the finite element mesh vertices.
 class Less_MVertex{
   public:
    bool operator()(const MVertex* v1, const MVertex* v2) const {
@@ -324,14 +324,13 @@ class CellComplex
    // used in relative homology computation, may be empty
    std::vector<GEntity*> _subdomain;
    
-   // sorted container of unique cells in this cell complex 
-   // mapped according to their dimension
-   std::map< int, std::set<Cell*, Less_Cell> >  _cells;
+   // sorted containers of unique cells in this cell complex 
+   // one for each dimension
+   std::set<Cell*, Less_Cell>  _cells[4];
     
-   // boundary operator matrices for this cell complex
+   // boundary operator matrices for this chain complex
    // h_k: C_k -> C_(k-1)
-   std::map<int, gmp_matrix*>  _HMatrix; 
-   
+   //gmp_matrix* _HMatrix[4];
    
  public:
    // iterator for the cells of same dimension
@@ -355,6 +354,11 @@ class CellComplex
    virtual void insertCells(bool subdomain);
    
   public: 
+   CellComplex( std::set<Cell*, Less_Cell>* cells ) {
+     for(int i = 0; i < 4; i++){
+     _cells[i] = cells[i]; 
+     }
+   }
    CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
    virtual ~CellComplex(){}
    
@@ -404,20 +408,53 @@ class CellComplex
    // slower, but produces cleaner result
    virtual int coreductionMrozek(Cell* generator);
    
+   // print the vertices of cells of certain dimension
+   virtual void printComplex(int dim);
+   
+   // write this cell complex in legacy .msh format
+   virtual int writeComplexMSH(const std::string &name); 
    
     // construct boundary operator matrix dim->dim-1
-   virtual void constructHMatrix(int dim);
+   //virtual void constructHMatrix(int dim);
+   // get the boundary operator matrix dim->dim-1
+   //virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+  
+   virtual std::vector<gmp_matrix*> constructHMatrices();
+   
+};
+
+class ChainComplex{
+  private:
+   // boundary operator matrices for this chain complex
+   // h_k: C_k -> C_(k-1)
+   gmp_matrix* _HMatrix[4];
+   
+  public:
+   
+   ChainComplex( std::vector<gmp_matrix*> HMatrix ){
+     for(int i = 0; i < HMatrix.size(); i++){
+       _HMatrix[i] = HMatrix.at(i);
+     }  
+   }
+   ChainComplex(){
+     for(int i = 0; i < 4; i++){
+       _HMatrix[i] = create_gmp_matrix_zero(1,1);
+     }
+   }
+   virtual ~ChainComplex(){}
+   
    // get the boundary operator matrix dim->dim-1
    virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
    
-   // print the vertices of cells of certain dimension
-   virtual void printComplex(int dim);
+   virtual gmp_matrix* ker(gmp_matrix* HMatrix);
    
-   // write this cell complex in legacy .msh format
-   virtual int writeComplexMSH(const std::string &name); 
+   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); } 
      
 };
 
+
 #endif
 
 #endif

utils/misc/celldriver.cpp

 // configure, compile and install Gmsh as a library with
 //
-//   ./configure --disable-gui --disable-netgen --disable-chaco 
-//        --disable-metis --disable-tetgen --prefix=/usr/local/
+//   ./configure --disable-gui --disable-netgen --disable-chaco --prefix=/usr/local/
 //   make install-lib
 //
-// Then compile this driver with "g++ kaka.cpp -lGmsh -llapack -lblas -lgmp"
+// 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.
+
 
 #include <stdio.h>
 #include <gmsh/Gmsh.h>

utils/misc/transformer.geo

@@ -89,11 +89,11 @@ Line Loop(57) = {34, -39, -38, -44};
 Ruled Surface(58) = {57};
 Line Loop(59) = {33, 44, -37, -42};
 Ruled Surface(60) = {59};
-Line Loop(61) = {28, -1, 25, 16, 29, -11, -32, 8};
-Line Loop(62) = {33, 34, 35, 36};
-Ruled Surface(63) = {61, 62};
-Line Loop(64) = {3, -26, -14, -30, 9, 31, -6, -27};
-Line Loop(65) = {37, 38, -40, 41};
-Ruled Surface(66) = {64, 65};
-Surface Loop(67) = {63, 46, 66, 18, 48, 24, 50, 22, 52, 20, 54, 60, 58, 56};
+Line Loop(61) = {27, 6, -31, -9, 30, 14, 26, -3};
+Line Loop(62) = {37, 38, -40, 41};
+Plane Surface(63) = {61, 62};
+Line Loop(64) = {25, 16, 29, -11, -32, 8, 28, -1};
+Line Loop(65) = {34, 35, 36, 33};
+Plane Surface(66) = {64, 65};
+Surface Loop(67) = {66, 48, 24, 63, 46, 20, 52, 22, 50, 18, 54, 60, 58, 56};
 Volume(68) = {67};