Added Geo/Homology.h and Geo/Homology.cpp as mid-interface for homology
computation (the role of utils/misc/celldriver.cpp so far). From there on
it's easy to either make plugin interface or parser interface, or both.

Made Homology Computation plugin usable.

The IO & graphics for post-processing views of the results of the homology
computation are crude at the moment: it saves the results to a separate
.msh file and reads them back by merging that to the current Gmsh session
(it uses File->Merge method).

Geo/CellComplex.cpp

@@ -268,6 +268,8 @@ void CellComplex::removeCell(Cell* cell){
     bdCell->removeCoboundaryCell(cell);
   }
 
+  //_trash.push_back(cell);
+  
 }
 
 void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
@@ -360,7 +362,6 @@ void CellComplex::reduceComplex(bool omitHighdim){
   for(int i = 3; i > 0; i--) count = count + reduction(i);
   
   
-  
   if(count == 0 && omitHighdim){
     
     removeSubdomain();
@@ -737,7 +738,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
     if(vertex->inSubdomain()) physical = 3;
     else if(vertex->onDomainBoundary()) physical = 2;
     else physical = 1;
-     fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0)->getNum());
+     fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getNum(), 15, 3, 0, 0, physical, vertex->getVertex(0)->getNum());
   }
   
   std::list< std::pair<int, Cell*> > cells;
@@ -749,7 +750,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
     cells = edge->getCells();
     for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
       Cell* cell = (*it).second;
-      fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getTag(), 1, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+      fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
     }
   }
   
@@ -761,7 +762,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
     cells = face->getCells();
     for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
       Cell* cell = (*it).second;
-      fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getTag(), 2, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+      fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
     }
   }
   for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
@@ -772,7 +773,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
     cells = volume->getCells();
     for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
       Cell* cell = (*it).second;
-      fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getTag(), 4, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+      fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
     }
   }
     

Geo/CellComplex.h

@@ -62,7 +62,7 @@ class Cell
    virtual void setTag(int tag) { _tag = tag; };
    virtual int getIndex() const { return _index; };
    virtual void setIndex(int index) { _index = index; };
-   
+   virtual int getNum() { return -1; }
    
    // get the number of vertices this cell has
    virtual int getNumVertices() const = 0;
@@ -211,6 +211,7 @@ class ZeroSimplex : public Simplex, public MPoint
    ~ZeroSimplex(){}
    
    int getDim() const { return 0; }
+   int getNum() { return MPoint::getNum(); }
    int getNumVertices() const { return 1; }
    MVertex* getVertex(int vertex) const {return _v[0]; }
    int getSortedVertex(int vertex) const {return _v[0]->getNum(); }
@@ -244,6 +245,7 @@ class OneSimplex : public Simplex, public MLine
    ~OneSimplex(){}
    
    int getDim() const { return 1; }
+   int getNum() { return MLine::getNum(); }
    int getNumVertices() const { return 2; }
    int getNumFacets() const {  return 2; }
    MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -289,6 +291,7 @@ class TwoSimplex : public Simplex, public MTriangle
    ~TwoSimplex(){}
    
    int getDim() const { return 2; }
+   int getNum() { return MTriangle::getNum(); }
    int getNumVertices() const { return 3; }
    int getNumFacets() const { return 3; }
    MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -331,6 +334,7 @@ class ThreeSimplex : public Simplex, public MTetrahedron
    ~ThreeSimplex(){}
    
    int getDim() const { return 3; }
+   int getNum() { return MTetrahedron::getNum(); }
    int getNumVertices() const { return 4; }
    int getNumFacets() const { return 4; }
    MVertex* getVertex(int vertex) const {return _v[vertex]; }
@@ -513,6 +517,8 @@ class CellComplex
    // one for each dimension
    std::set<Cell*, Less_Cell>  _cells[4];
    
+   std::vector<Cell*> _trash;
+   
    //std::set<Cell*, Less_Cell>  _originalCells[4];
    
    // Betti numbers of this cell complex (ranks of homology groups)
@@ -570,7 +576,12 @@ class CellComplex
    
    CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
    CellComplex(){}
-   ~CellComplex(){}
+   ~CellComplex(){ 
+     for(int i = 0; i < _trash.size(); i++){
+       Cell* cell = _trash.at(i);
+       delete cell;
+     }
+   }
 
    
    // get the number of certain dimensional cells

Geo/ChainComplex.cpp

@@ -487,7 +487,7 @@ int Chain::writeChainMSH(const std::string &name){
     cells = cell->getCells();
     for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
       Cell* cell2 = (*it).second;
-      fprintf(fp, "%d %d \n", cell2->getTag(), getCoeff(i)*(*it).first );
+      fprintf(fp, "%d %d \n", cell2->getNum(), getCoeff(i)*(*it).first );
     }
   }
   
@@ -498,3 +498,29 @@ int Chain::writeChainMSH(const std::string &name){
   return 1;
   
 }
+
+void Chain::getData(std::map<int, std::vector<double> > & data){
+  
+  if(getSize() == 0) return;
+  
+  std::list< std::pair<int, Cell*> > cells;
+  for(int i = 0; i < getSize(); i++){
+    Cell* cell = getCell(i);
+    cells = cell->getCells();
+    for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+      Cell* cell2 = (*it).second;
+      std::vector<double> coeff;
+      coeff.push_back(getCoeff(i)*(*it).first);
+      std::pair<int, std::vector<double> >  dataPair = std::make_pair(cell2->getTag(), coeff);
+      data.insert(dataPair);
+      //printf("%d, %d, \n", cell2->getNum(), (int)coeff.at(0));
+      
+    }
+  }
+  
+  //for(std::map<int, std::vector<double> >::iterator it = data.begin(); it != data.end(); it++){
+  //  printf("%d, %d, \n", (*it).first, (int)(*it).second.at(0));  
+  //}
+  
+  return; 
+}

Geo/ChainComplex.h

@@ -160,6 +160,8 @@ class Chain{
    // append this chain to a 2.0 MSH ASCII file as $ElementData
    virtual int writeChainMSH(const std::string &name);
    
+   virtual void getData(std::map<int, std::vector<double> >& data);
+   
 };
 
 #endif

Geo/Homology.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka.
+ 
+
+#include "Homology.h"
+#include "ChainComplex.h"
+
+
+Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
+  
+  _model = model;
+  
+  Msg::Info("Creating a Cell Complex...");
+  
+  
+  std::map<int, std::vector<GEntity*> > groups[4];
+  model->getPhysicalGroups(groups);
+  
+  
+  std::map<int, std::vector<GEntity*> >::iterator it;
+  std::vector<GEntity*> domainEntities;
+  std::vector<GEntity*> subdomainEntities;
+    
+  for(int i = 0; i < physicalDomain.size(); i++){
+    for(int j = 0; j < 4; j++){
+      it = groups[j].find(physicalDomain.at(i));
+      if(it != groups[j].end()){
+        std::vector<GEntity*> physicalGroup = (*it).second;
+        for(int k = 0; k < physicalGroup.size(); k++){
+          domainEntities.push_back(physicalGroup.at(k));
+        }
+      }
+    }
+  }
+  for(int i = 0; i < physicalSubdomain.size(); i++){           
+    for(int j = 0; j < 4; j++){
+      it = groups[j].find(physicalSubdomain.at(i));
+      if(it != groups[j].end()){
+        std::vector<GEntity*> physicalGroup = (*it).second;
+        for(int k = 0; k < physicalGroup.size(); k++){
+          subdomainEntities.push_back(physicalGroup.at(k));
+        }
+        
+      }
+    }
+  }
+  
+  if(domainEntities.empty()) Msg::Warning("Domain is empty.");
+  if(subdomainEntities.empty()) Msg::Info("Subdomain is empty.");
+  
+  
+  _cellComplex =  new CellComplex(domainEntities, subdomainEntities);
+  
+  if(_cellComplex->getSize(0) == 0){ 
+    Msg::Error("Cell Complex is empty!");
+    Msg::Error("Check the domain & the mesh.");
+    return;
+  }
+  
+  Msg::Info("Cell Complex complete.");
+  Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
+            _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
+  
+}
+
+
+void Homology::findGenerators(std::string fileName){
+  
+  Msg::Info("Reducing Cell Complex...");
+  _cellComplex->reduceComplex(true);
+  _cellComplex->combine(3);
+  _cellComplex->combine(2);
+  _cellComplex->combine(1);
+  Msg::Info("Cell Complex reduction complete.");
+  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...");
+  ChainComplex* chains = new ChainComplex(_cellComplex);
+  chains->computeHomology();
+  Msg::Info("Homology Computation complete.");
+  
+  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(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name);
+      chain->writeChainMSH(fileName);
+      delete chain;
+    }
+  }
+  Msg::Info("Wrote results to %s.", fileName.c_str());
+  
+  delete chains;
+  
+  return;
+}
+
+void Homology::findThickCuts(std::string fileName){
+  
+  Msg::Info("Reducing Cell Complex...");
+  _cellComplex->coreduceComplex(true);
+  _cellComplex->cocombine(0);
+  _cellComplex->cocombine(1);
+  _cellComplex->cocombine(2);
+  Msg::Info("Cell Complex reduction complete.");
+  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...");
+  ChainComplex* chains = new ChainComplex(_cellComplex);
+  chains->transposeHMatrices();
+  chains->computeHomology(true);
+  Msg::Info("Homology Computation complete.");
+  
+  for(int j = 3; j > -1; j--){
+    for(int i = 1; i <= chains->getBasisSize(j); i++){
+      
+      std::string generator;
+      std::string dimension;
+      convert(i, generator);
+      convert(3-j, dimension);
+      
+      std::string name = dimension + "D Dual cut " + generator;
+      Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name);
+      chain->writeChainMSH(fileName);
+      delete chain;
+            
+    }
+  }
+  
+  Msg::Info("Wrote results to %s.", fileName.c_str());
+  
+  delete chains;
+  
+}

Geo/Homology.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+// 
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+// 
+// Contributed by Matti Pellikka.
+
+#ifndef _HOMOLOGY_H_
+#define _HOMOLOGY_H_
+
+#include <sstream>
+
+#include "CellComplex.h"
+
+template <class TTypeA, class TTypeB>
+bool convert(const TTypeA& input, TTypeB& output ){
+   std::stringstream stream;
+   stream << input;
+   stream >> output;
+   return stream.good();
+}
+
+
+class Homology
+{
+  private:
+   CellComplex* _cellComplex;
+   GModel* _model;
+   
+  public:
+   
+   Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
+   ~Homology(){ delete _cellComplex; }
+   
+   void findGenerators(std::string fileName);
+   void findThickCuts(std::string fileName);
+      
+   void swapSubdomain() { _cellComplex->swapSubdomain(); }
+      
+};
+
+#endif

Geo/Makefile

@@ -40,8 +40,7 @@ SRC = GEntity.cpp STensor3.cpp\
         MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp\
         MHexahedron.cpp MPrism.cpp MPyramid.cpp\
       MZone.cpp MZoneBoundary.cpp\
-      CellComplex.cpp\
-      ChainComplex.cpp\
+      CellComplex.cpp ChainComplex.cpp Homology.cpp\
       STensor3.cpp
 
 
@@ -413,3 +412,4 @@ ChainComplex${OBJEXT}: ChainComplex.cpp ChainComplex.h ../Common/GmshConfig.h \
 STensor3${OBJEXT}: STensor3.cpp STensor3.h SVector3.h SPoint3.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
   ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
+Homology${OBJEXT}: Homology.cpp Homology.h CellComplex.h ChainComplex.h

Plugin/HomologyComputation.cpp

@@ -6,18 +6,24 @@
 // Contributed by Matti Pellikka
 
 #include <stdlib.h>
+#include "Gmsh.h"
 #include "GmshConfig.h"
 #include "GModel.h"
-#include "CellComplex.h"
+#include "Homology.h"
+#include "PViewDataGModel.h"
 #include "HomologyComputation.h"
 
 
 StringXNumber HomologyComputationOptions_Number[] = {
-  {GMSH_FULLRC, "test number", NULL, 1.},
+  {GMSH_FULLRC, "Physical group for domain", NULL, 0.},
+  {GMSH_FULLRC, "Physical group for subdomain", NULL, 0.},
+  {GMSH_FULLRC, "Compute generators", NULL, 1.},
+  {GMSH_FULLRC, "Compute thick cuts", NULL, 0.},
+  {GMSH_FULLRC, "Swap subdomain", NULL, 0.},
 };
 
 StringXString HomologyComputationOptions_String[] = {
-  {GMSH_FULLRC, "test string", NULL, "test string"},
+  {GMSH_FULLRC, "Filename", NULL, "homology.msh"}
 };
 
 extern "C"
@@ -36,12 +42,13 @@ void GMSH_HomologyComputationPlugin::getName(char *name) const
 void GMSH_HomologyComputationPlugin::getInfos(char *author, char *copyright,
                                    char *help_text) const
 {
-  strcpy(author, "Matti Pellikka");
+  strcpy(author, "Matti Pellikka (matti.pellikka@tut.fi)");
   strcpy(copyright, "C. Geuzaine, J.-F. Remacle");
   strcpy(help_text,
-         "Plugin(HomologyComputation) is here!\n"
+         "Plugin(HomologyComputation) computes generators \n"
+         "for (relative) homology groups and their thick cuts. \n"
          "\n"
-         "Plugin(HomologyComputation) creates a new view.\n");
+         "Plugin(HomologyComputation) creates new views.\n");
 }
 
 int GMSH_HomologyComputationPlugin::getNbOptions() const
@@ -72,28 +79,39 @@ void GMSH_HomologyComputationPlugin::catchErrorMessage(char *errorMessage) const
 
 PView *GMSH_HomologyComputationPlugin::execute(PView *v)
 {
-  std::string testString = HomologyComputationOptions_String[0].def;
-  int testInt = (int)HomologyComputationOptions_Number[0].def;
+  std::string fileName = HomologyComputationOptions_String[0].def;
+  if(fileName.empty()) {
+    Msg::Error("Filename not given!");
+    return 0;
+  }
   
-  GModel *m = GModel::current();
-  std::map<int, std::vector<GEntity*> > groups[4];
-  m->getPhysicalGroups(groups);
+  std::vector<int> domain;
+  std::vector<int> subdomain;
   
-  std::map<int, std::vector<double> > data;
+  domain.push_back((int)HomologyComputationOptions_Number[0].def);
+  subdomain.push_back((int)HomologyComputationOptions_Number[1].def);
 
-  std::vector<GEntity*> domain;
-  std::vector<GEntity*> subdomain;
+  int gen = (int)HomologyComputationOptions_Number[2].def;
+  int cuts = (int)HomologyComputationOptions_Number[3].def;
+  int swap = (int)HomologyComputationOptions_Number[4].def;
 
-  // Here's the problem, linker cannot find the method for the constructor of 
-  // CellComplex class, or any other methods coded in Geo/CellComplex.cpp
+  GModel *m = GModel::current();
   
-  //CellComplex* testComplex = new CellComplex(domain, subdomain);
+  Homology* homology = new Homology(m, domain, subdomain);
   
-  // insert homology fun here
+  if(swap == 1) homology->swapSubdomain();
+  if(gen == 1 && cuts != 1) {
+    homology->findGenerators(fileName);
+    GmshMergeFile(fileName);
+  }
+  else if(cuts == 1 && gen != 1) {
+    homology->findThickCuts(fileName);
+    GmshMergeFile(fileName);
+  }
+  else Msg::Error("Choose either generators or thick cuts to compute.");
   
-  PView *pv = new PView("homology", "ElementData", m, data, 0.);
+  delete homology; 
   
-  Msg::Error("test.");
   
   return 0;
 }

Plugin/Makefile

@@ -342,4 +342,4 @@ FiniteElement${OBJEXT}: FiniteElement.cpp ../Common/GmshConfig.h ../Geo/GModel.h
   ../Numeric/gmshLinearSystem.h
 HomologyComputation${OBJEXT}: HomologyComputation.cpp Plugin.h \
   ../Geo/GModel.h ../Geo/GEntity.h ../Geo/CellComplex.h \
-  ../Post/PView.h HomologyComputation.h
+  ../Post/PView.h HomologyComputation.h ../Geo/Homology.h