Tuning.

Geo/Cell.h

@@ -84,6 +84,23 @@ class Cell
   
   virtual MElement* getImageMElement() const { return _image; };
   
+
+  // get the number of vertices this cell has
+  virtual int getNumVertices() const { return _image->getNumVertices(); }
+  // get the number of facets of this cell
+  virtual int getNumFacets() const;
+  // get the vertices on a facet of this cell
+  virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
+  // get boundary cell orientation
+  virtual int getFacetOri(std::vector<MVertex*> &v); 
+  virtual int getFacetOri(Cell* cell) { 
+    std::vector<MVertex*> v; 
+    for(int i = 0; i < cell->getNumVertices(); i++) {
+      v.push_back(cell->getVertex(i));
+    }
+    return getFacetOri(v);
+  }
+
   virtual int getDim() const { return _dim; };
   virtual int getIndex() const { return _index; };
   virtual void setIndex(int index) { _index = index; };
@@ -94,9 +111,6 @@ class Cell
   virtual void setDeleteImage(bool deleteImage) { 
     _deleteImage = deleteImage; };
   virtual bool getDeleteImage() const { return _deleteImage; };
-  
-  // get the number of vertices this cell has
-  virtual int getNumVertices() const { return _image->getNumVertices(); }
   virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
   
   // restores the cell information to its original state before reduction
@@ -147,21 +161,6 @@ class Cell
   virtual void printBoundary(bool org=false);
   virtual void printCoboundary(bool org=false);   
   
-  // get the number of facets of this cell
-  virtual int getNumFacets() const;
-  // get the vertices on a facet of this cell
-  virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
-  
-  // get boundary cell orientation
-  virtual int getFacetOri(std::vector<MVertex*> &v); 
-  virtual int getFacetOri(Cell* cell) { 
-    std::vector<MVertex*> v; 
-    for(int i = 0; i < cell->getNumVertices(); i++) {
-      v.push_back(cell->getVertex(i));
-    }
-    return getFacetOri(v);
-  }   
-  
   // tools for combined cells
   virtual bool isCombined() { return _combined; }
   virtual void getCells( std::list< std::pair<int, Cell*> >& cells) {  
@@ -183,6 +182,14 @@ class Cell
     }
     return true;
   }
+  /*
+  Cell operator=(const Cell& c2) {
+    Cell cell;
+    cell._ocbd = c2._ocbd;
+
+    return cell;
+  }
+  Cell(const Cell& c2){ *this = c2; }*/
 };
 
 // A cell that is a combination of cells of same dimension
@@ -194,17 +201,20 @@ class CombinedCell : public Cell{
   // list of cells this cell is a combination of
   std::list< std::pair<int, Cell*> > _cells;
   
-  MVertex* getVertex(int vertex) const { return NULL; }
+  MVertex* getVertex(int vertex) const {
+    printf("ERROR: No mesh vertices for combined cell."); } 
   
  public:
   
   CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
   ~CombinedCell() {}
   
-  int getNum() const { return 0; }
-  int getType() const { return 0; }
-  int getTypeForMSH() const { return 0; }
-  int getPartition() const { return 0; }
+  MElement* getImageMElement() const { 
+    printf("ERROR: No image mesh element for combined cell."); }
+  int getNumFacets() const { return 0; }
+  void getFacetVertices(const int num, std::vector<MVertex*> &v) const {}
+  int getFacetOri(std::vector<MVertex*> &v) { return 0; }
+  int getFacetOri(Cell* cell) { return 0; }
 
   int getNumVertices() const { return _vs.size(); } 
   int getSortedVertex(int vertex) const { return _vs.at(vertex); }

Geo/CellComplex.cpp

@@ -621,4 +621,23 @@ void  CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
     }
   }
 }
+
+void CellComplex::restoreComplex()
+{
+  for(int i = 0; i < 4; i++){
+    _cells[i] = _ocells[i];
+    for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+      Cell* cell = *cit;
+      cell->restoreCell();
+    }
+  }
+  for(unsigned int i = 0; i < _newcells.size(); i++){
+    Cell* cell = _newcells.at(i);
+    delete cell;
+  }
+  _newcells.clear();
+  _store.clear();
+}
+
+
 #endif

Geo/CellComplex.h

@@ -131,6 +131,8 @@ class CellComplex
   int getNumOmitted() { return _store.size(); }
   std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }  
 
+  void restoreComplex();
+  
 };
 
 #endif

Geo/Homology.cpp

@@ -108,10 +108,14 @@ Homology::~Homology()
   
 }
 
-void Homology::findGenerators()
+void Homology::findGenerators(CellComplex* cellComplex)
 {
+  bool ownComplex = false;
+  if(cellComplex==NULL){
     CellComplex* cellComplex = createCellComplex(_domainEntities, 
 						 _subdomainEntities);
+    ownComplex = true;
+  }
   std::string domainString = getDomainString(_domain, _subdomain);
 
   Msg::Info("Reducing the Cell Complex...");
@@ -206,7 +210,7 @@ void Homology::findGenerators()
   
   if(_fileName != "") writeGeneratorsMSH();
   
-  delete cellComplex;
+  if(ownComplex) delete cellComplex;
   delete chains;
   
   Msg::Info("Ranks of homology spaces for primal cell complex:");
@@ -221,10 +225,14 @@ void Homology::findGenerators()
 		 HRank[0], HRank[1], HRank[2], HRank[3]);
 }
 
-void Homology::findDualGenerators()
+void Homology::findDualGenerators(CellComplex* cellComplex)
 { 
+  bool ownComplex = false;
+  if(cellComplex==NULL){
     CellComplex* cellComplex = createCellComplex(_domainEntities, 
 						 _subdomainEntities);
+    ownComplex = true;
+  }
 
   Msg::Info("Reducing Cell Complex...");
   Msg::StatusBar(1, false, "Reducing...");
@@ -321,7 +329,7 @@ void Homology::findDualGenerators()
 
   if(_fileName != "") writeGeneratorsMSH();
 
-  delete cellComplex;
+  if(ownComplex) delete cellComplex;
   delete chains;
   
   Msg::Info("Ranks of homology spaces for the dual cell complex:");

Geo/Homology.h

@@ -42,8 +42,7 @@ class Homology
   std::vector<GEntity*> _domainEntities;
   std::vector<GEntity*> _subdomainEntities;  
 
-  // generator chains, and their physical group number
-  //std::map<int, Chain*> _generators;
+  // physical group numbers of generator chains in model
   std::vector<int> _generators;
 
   std::string _fileName;
@@ -54,17 +53,20 @@ class Homology
 	   std::vector<int> physicalSubdomain);
   ~Homology();
   
+  
   CellComplex* createCellComplex(std::vector<GEntity*>& domainEntities,
 				 std::vector<GEntity*>& subdomainEntities);
+  CellComplex* createCellComplex() { 
+    return createCellComplex(_domainEntities, _subdomainEntities); }
 
   void setFileName(std::string fileName) { _fileName = fileName; }
 
   // Find the generators/duals of homology spaces,
   // or just compute the ranks of homology spaces
-  void findGenerators();
-  void findDualGenerators();
-  void computeRanks() {}
+  void findGenerators(CellComplex* cellComplex=NULL);
+  void findDualGenerators(CellComplex* cellComplex=NULL);
 
+  void computeRanks() {}  
   void findHomSequence();
 
    

Plugin/HomologyComputation.cpp

@@ -21,7 +21,7 @@ StringXNumber HomologyComputationOptions_Number[] = {
   {GMSH_FULLRC, "PhysicalGroupForSubdomain2", NULL, 0.},
   {GMSH_FULLRC, "ComputeGenerators", NULL, 1.},
   {GMSH_FULLRC, "ComputeCuts", NULL, 0.},
-  {GMSH_FULLRC, "ComputeRanks", NULL, 0.},
+  //{GMSH_FULLRC, "ComputeRanks", NULL, 0.},
 };
 
 StringXString HomologyComputationOptions_String[] = {
@@ -89,28 +89,26 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
 
   int gens = (int)HomologyComputationOptions_Number[4].def;
   int cuts = (int)HomologyComputationOptions_Number[5].def;
-  int rank = (int)HomologyComputationOptions_Number[6].def;  
+  //int rank = (int)HomologyComputationOptions_Number[6].def;  
 
   GModel *m = GModel::current();
   
-  if(gens != 0){
   Homology* homology = new Homology(m, domain, subdomain);
   homology->setFileName(fileName);
-    homology->findGenerators();
-    delete homology; 
+  CellComplex* cc = homology->createCellComplex();
+  if(gens != 0){
+    homology->findGenerators(cc);
   }
   if(cuts != 0){
-    Homology* homology = new Homology(m, domain, subdomain);
-    homology->setFileName(fileName);
-    homology->findDualGenerators();
-    delete homology; 
+    cc->restoreComplex();
+    homology->findDualGenerators(cc);
   }
-  if(rank != 0){
-    Homology* homology = new Homology(m, domain, subdomain);
-    homology->setFileName(fileName);
+  /*if(rank != 0){
     homology->computeRanks();
+    }*/
+  
+  delete cc;
   delete homology;
-  }
 
   return 0;
 }

utils/api_demos/mainHomology.cpp

@@ -12,7 +12,7 @@
 #include <gmsh/GModel.h>
 #include <gmsh/MElement.h>
 #include <gmsh/CellComplex.h>
-#include <gmsh/ChainComplex.h>
+//#include <gmsh/ChainComplex.h>
 #include <gmsh/Homology.h>
 
 int main(int argc, char **argv)
@@ -25,21 +25,31 @@ int main(int argc, char **argv)
   // OR
   // m->readMSH("model.msh");
   
-  // List of physical regions as domain for homology computation (relative to subdomain).
+  // List of physical regions as domain for homology computation
+  // (relative to subdomain).
   std::vector<int> domain;
   std::vector<int> subdomain;
   
-  Homology* homology1 = new Homology(m, domain, subdomain);
-  std::string fileName = "homology.msh";
-  homology->findGenerators();
-  homology->setFileName(fileName);
-  delete homology1;
+  // initialize
+  Homology* homology = new Homology(m, domain, subdomain);
 
-  Homology* homology2 = new Homology(m, domain, subdomain);
-  homology->findDualGenerators();
-  homology->setFileName(fileName);
-  delete homology2;
+  // save all resulting chains to a file as physical groups
+  homology->setFileName("homology.msh");
 
+  // construct cell complex of the meshed domain
+  CellComplex cc = homology->createCellComplex();
+
+  // find homology generator chains
+  homology->findGenerators(cc);
+
+  // restore cell complex for a new run
+  cc->restoreComplex(); 
+
+  // find thick cuts
+  homology->findDualGenerators(cc);
+
+  delete cc;
+  delete homology;
   delete m;
   GmshFinalize();