Faster 2-chain smoothening. Updated api demo.

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.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
 
 #include "CellComplex.h"
 #include "OS.h"

Geo/ChainComplex.cpp

@@ -783,8 +783,24 @@ void Chain::smoothenChain(){
   
   int start = getSize();
   double t1 = Cpu();
-  const int MAXROUNDS = 2;
   
+  int useless = 0;
+  for(int i = 0; i < 20; i++){
+    int size = getSize();
+    for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+      if(!straightenChain(*cit) && getDim() == 2) bendChain(*cit);
+      removeBoundary(*cit);
+    }
+    deImmuneCells();
+    eraseNullCells();
+    if (size >= getSize()) useless++;
+    else useless = 0;
+    if (useless > 5) break;
+  }
+  
+  
+  /*
+  const int MAXROUNDS = 2;
   bool smoothened = true;
   int round = 0;
   while(smoothened){
@@ -855,10 +871,13 @@ void Chain::smoothenChain(){
       }
       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;

Geo/ChainComplex.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 _CHAINCOMPLEX_H_
 #define _CHAINCOMPLEX_H_

Geo/Homology.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 "Homology.h"

Geo/Homology.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_H_
 #define _HOMOLOGY_H_
@@ -36,15 +36,19 @@ class Homology
    Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
    ~Homology(){ delete _cellComplex; }
    
+   // Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
    void findGenerators(std::string fileName);
    void findDualGenerators(std::string fileName);
-
    void computeBettiNumbers();
    
+   
    void swapSubdomain() { _cellComplex->swapSubdomain(); }
    
-   std::string getDomainString() {
+   // Restore the cell complex to its original state before cell reductions
+   void restoreHomology() { _cellComplex->restoreComplex(); }
    
+   // Create a string describing the generator
+   std::string getDomainString() {
      std::string domainString = "({";
      for(unsigned int i = 0; i < _domain.size(); i++){
        std::string temp = "";

utils/api_demos/mainHomology.cpp

-// configure, compile and install Gmsh as a library with
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
 // 
-//   ./configure --disable-gui --disable-netgen --disable-chaco --prefix=/usr/local/
-//   make install-lib
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
 // 
-// Then compile this driver with "g++ celldriver.cpp -lGmsh -llapack -lblas -lgmp"
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
 // 
-//
-// This program creates .msh file chains.msh which represents the generators and
-// the cuts of an two port transformer model's homology groups.
-
 
 #include <stdio.h>
 #include <sstream>
@@ -17,130 +13,36 @@
 #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();
-}
-
+#include <gmsh/Homology.h>
 
 int main(int argc, char **argv)
 {
   GmshInitialize(argc, argv);
   GModel *m = new GModel();
-  m->readGEO("transformer.geo");
-  m->mesh(3);
-  
-  printf("This model has: %d GRegions, %d GFaces, %d GEdges and %d GVertices. \n" , m->getNumRegions(), m->getNumFaces(), m->getNumEdges(), m->getNumVertices());
-  
-  std::vector<GEntity*> domain;
-  std::vector<GEntity*> subdomain;
-  
-  // whole model the domain
-  for(GModel::riter rit = m->firstRegion(); rit != m->lastRegion(); rit++){
-    GEntity *r = *rit;
-    domain.push_back(r);
-  }
-
-  // the ports as subdomain
-  GModel::fiter sub = m->firstFace();
-  GEntity *s = *sub;
-  subdomain.push_back(s);
-  s = *(++sub);
-  subdomain.push_back(s);
-  s =*(++sub);
-  subdomain.push_back(s);
-  s= *(++sub);
-  subdomain.push_back(s);
-  
-  // create a cell complex
-  CellComplex* complex = new CellComplex(domain, subdomain);
-
-  // write the cell complex to a .msh file
-  complex->writeComplexMSH("chains.msh");
-
-  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));  
-  
-  // reduce the complex in order to ease homology computation
-  complex->reduceComplex(true);
-  
-  // perform series of cell combinatios in order to further ease homology computation
-  complex->combine(3);
-  complex->combine(2);
-  complex->combine(1);
-  
-  // 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();
   
-  // 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;
-      
-    }
-  }
-  
-  delete chains;
-  delete complex;
-  
-  // create new complex to compute the cuts
-  CellComplex* complex2 = new CellComplex(domain, subdomain);
-  
-  // reduce the complex by coreduction, this removes all vertices, hence 0 as an argument 
-  complex2->coreduceComplex(true);
-  
-  // perform series of "dual complex" cell combinations in order to ease homology computation
-  complex2->cocombine(0);
-  complex2->cocombine(1);
-  complex2->cocombine(2);
-  
-  // create a chain complex of a cell complex (construct boundary operator matrices)
-  ChainComplex* dualChains = new ChainComplex(complex2);
-  
-  // transpose the boundary operator matrices, 
-  // these are the boundary operator matrices for the dual chain complex
-  dualChains->transposeHMatrices();
-  
-  // compute the homology of the dual complex
-  dualChains->computeHomology(true);  
-  
-  // Append the duals of the cuts of the homology generators to the .msh file.
-  for(int j = 3; j > -1; j--){
-    for(int i = 1; i <= dualChains->getBasisSize(j); i++){
+  m->readGEO("model.geo");
+  m->mesh(3);
+  // OR
+  // m->readMSH("model.msh");
   
-      std::string generator;
-      std::string dimension;
-      convert(i, generator);
-      convert(3-j, dimension);
+  // List of physical regions as domain for homology computation (relative to subdomain).
+  std::vector<int> domain;
+  std::vector<int> subdomain;
   
-      std::string name = dimension + "D Dual cut " + generator;
-      Chain* chain = new Chain(complex2->getCells(j), dualChains->getCoeffVector(j,i), complex2, name);
-                chain->writeChainMSH("chains.msh");
-      delete chain;
-    }
-  }
+  Homology* homology = new Homology(m, domain, subdomain);
+  std::string fileName = "homology_generators.msh";
+  homology->findGenerators(fileName);
+  homology->restoreHomology();
   
-  delete dualChains;
+  Homology* homology = new Homology(m, domain, subdomain);
+  fileName = "homology_dualgenerators.msh";
+  homology->findDualGenerators(fileName);
+  homology->restoreHomology();
   
+  Homology* homology = new Homology(m, domain, subdomain);
+  homology->computeBettiNumbers();
+  delete homology;
   
-  delete complex2;  
   delete m;
   GmshFinalize();