From 6f549b7a5957e5f024a6c063ebb4bff05660e960 Mon Sep 17 00:00:00 2001 From: Matti Pellika <matti.pellikka@tut.fi> Date: Wed, 3 Jun 2009 14:05:57 +0000 Subject: [PATCH] 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 | 13 +-- Geo/CellComplex.h | 15 +++- Geo/ChainComplex.cpp | 28 ++++++- Geo/ChainComplex.h | 2 + Geo/Homology.cpp | 149 +++++++++++++++++++++++++++++++++ Geo/Homology.h | 42 ++++++++++ Geo/Makefile | 4 +- Plugin/HomologyComputation.cpp | 62 +++++++++----- Plugin/Makefile | 2 +- 9 files changed, 283 insertions(+), 34 deletions(-) create mode 100644 Geo/Homology.cpp create mode 100644 Geo/Homology.h diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp index 6827a1ae31..e7b608f1fc 100644 --- a/Geo/CellComplex.cpp +++ b/Geo/CellComplex.cpp @@ -140,7 +140,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ int dim = domain.at(j)->getMeshElement(i)->getDim(); int type = domain.at(j)->getMeshElement(i)->getTypeForMSH(); Cell* cell; - if(dim == 3) cell = new ThreeSimplex(vertices, 0, subdomain, boundary); + if(dim == 3) cell = new ThreeSimplex(vertices, 0, subdomain, boundary); else if(dim == 2) cell = new TwoSimplex(vertices, 0, subdomain, boundary); else if(dim == 1) cell = new OneSimplex(vertices, 0, subdomain, boundary); else cell = new ZeroSimplex(vertices, 0, subdomain, boundary); @@ -267,6 +267,8 @@ void CellComplex::removeCell(Cell* cell){ Cell* bdCell = *it; bdCell->removeCoboundaryCell(cell); } + + //_trash.push_back(cell); } @@ -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()); } } diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h index a1f4948054..c4a415654f 100644 --- a/Geo/CellComplex.h +++ b/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 diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp index 365599a630..c5e3fca8c8 100644 --- a/Geo/ChainComplex.cpp +++ b/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; +} diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h index b8a5e2ec57..0ec42e3d6a 100644 --- a/Geo/ChainComplex.h +++ b/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 diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp new file mode 100644 index 0000000000..6ff1b406f7 --- /dev/null +++ b/Geo/Homology.cpp @@ -0,0 +1,149 @@ +// 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; + +} diff --git a/Geo/Homology.h b/Geo/Homology.h new file mode 100644 index 0000000000..1e805b0db7 --- /dev/null +++ b/Geo/Homology.h @@ -0,0 +1,42 @@ +// 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 diff --git a/Geo/Makefile b/Geo/Makefile index f3b4d07faf..dbb7036a82 100644 --- a/Geo/Makefile +++ b/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 diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp index 3f6b9ceacc..0587ec553d 100644 --- a/Plugin/HomologyComputation.cpp +++ b/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; - - GModel *m = GModel::current(); - std::map<int, std::vector<GEntity*> > groups[4]; - m->getPhysicalGroups(groups); - - std::map<int, std::vector<double> > data; + std::string fileName = HomologyComputationOptions_String[0].def; + if(fileName.empty()) { + Msg::Error("Filename not given!"); + return 0; + } + + std::vector<int> domain; + std::vector<int> subdomain; + + 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; } diff --git a/Plugin/Makefile b/Plugin/Makefile index a8b9869ea2..9d013f8352 100644 --- a/Plugin/Makefile +++ b/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 -- GitLab