From 8dfbe5dae6b64f9be44c86057a247d40a73443a6 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Fri, 2 Mar 2012 09:10:24 +0000
Subject: [PATCH] Chain class for homology computation post-processing
---
Geo/CMakeLists.txt | 2 +-
Geo/Cell.cpp | 27 +++
Geo/Chain.cpp | 231 ++++++++++++++++++++++++++
Geo/Chain.h | 382 +++++++++++++++++++++++++++++++++++++++++++
Geo/ChainComplex.cpp | 5 +-
Geo/GModel.cpp | 7 +-
Geo/Homology.cpp | 129 ++++-----------
Geo/Homology.h | 20 ++-
8 files changed, 698 insertions(+), 105 deletions(-)
create mode 100644 Geo/Chain.cpp
create mode 100644 Geo/Chain.h
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index aa2971af71..32bbeeb13d 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -40,7 +40,7 @@ set(SRC
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp
- Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp
+ Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp Chain.cpp
Curvature.cpp
MVertexBoundaryLayerData.cpp
)
diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 0c38a43c68..cec524dac0 100644
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -29,6 +29,14 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
return false;
}
+bool equalVertices(const std::vector<MVertex*>& v1,
+ const std::vector<MVertex*>& v2) {
+ if(v1.size() != v2.size()) return false;
+ for(unsigned int i = 0; i < v1.size(); i++)
+ if(v1[i]->getNum() != v2[i]->getNum()) return false;
+ return true;
+}
+
int Cell::_globalNum = 0;
Cell::Cell(MElement* element, int domain)
@@ -154,6 +162,25 @@ int Cell::getNumBdElements() const
int Cell::findBdCellOrientation(Cell* cell, int i) const
{
+ /*std::vector<MVertex*> v1;
+ std::vector<MVertex*> v2;
+ this->findBdElement(i, v1);
+ cell->getMeshVertices(v2);
+
+ int perm = 1;
+ if(equalVertices(v1, v2)) return perm;
+ while(std::next_permutation(v2.begin(), v2.end(), MVertexLessThanNum())) {
+ perm *= -1;
+ if(equalVertices(v1, v2)) return perm;
+ }
+ cell->getMeshVertices(v2);
+ perm = 1;
+ while(std::prev_permutation(v2.begin(), v2.end(), MVertexLessThanNum())) {
+ perm *= -1;
+ if(equalVertices(v1, v2)) return perm;
+ }
+ return 0;*/
+
std::vector<MVertex*> v;
cell->getMeshVertices(v);
switch (_dim) {
diff --git a/Geo/Chain.cpp b/Geo/Chain.cpp
new file mode 100644
index 0000000000..2751451020
--- /dev/null
+++ b/Geo/Chain.cpp
@@ -0,0 +1,231 @@
+// Gmsh - Copyright (C) 1997-2012 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 <matti.pellikka@tut.fi>.
+
+#include "Chain.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MPyramid.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+
+#if defined(HAVE_KBIPACK)
+
+std::map<GEntity*, std::set<MVertex*, MVertexLessThanNum>,
+ GEntityLessThan> ElemChain::_vertexCache;
+
+inline void ElemChain::_sortVertexIndices()
+{
+ std::map<MVertex*, int, MVertexLessThanNum> si;
+
+ for(unsigned int i = 0; i < _v.size(); i++)
+ si[_v[i]] = i;
+
+ std::map<MVertex*, int, MVertexLessThanNum>::iterator it;
+ for(it = si.begin(); it != si.end(); it++)
+ _si.push_back(it->second);
+}
+
+bool ElemChain::_equalVertices(const std::vector<MVertex*>& v2) const {
+ if(_v.size() != v2.size()) return false;
+ for(unsigned int i = 0; i < _v.size(); i++)
+ if(_v[i]->getNum() != v2[i]->getNum()) return false;
+ return true;
+}
+
+ElemChain::ElemChain(MElement* e)
+{
+ _dim = e->getDim();
+ for(int i = 0; i < e->getNumPrimaryVertices(); i++)
+ _v.push_back(e->getVertex(i));
+ _sortVertexIndices();
+}
+
+ElemChain::ElemChain(int dim, std::vector<MVertex*>& v) : _dim(dim), _v(v)
+{
+ _sortVertexIndices();
+}
+
+inline int ElemChain::getSortedVertex(int i) const
+{
+ return _v[(int)_si[i]]->getNum();
+}
+
+int ElemChain::getTypeMSH(int dim, int numVertices)
+{
+ switch (dim) {
+ case 0: return MSH_PNT;
+ case 1: return MSH_LIN_2;
+ case 2:
+ switch (numVertices) {
+ case 3: return MSH_TRI_3;
+ case 4: return MSH_QUA_4;
+ default: return 0;
+ }
+ case 3:
+ switch (numVertices) {
+ case 4: return MSH_TET_4;
+ case 5: return MSH_PYR_5;
+ case 6: return MSH_PRI_6;
+ case 8: return MSH_HEX_8;
+ default: return 0;
+ }
+ default: return 0;
+ }
+}
+
+int ElemChain::getTypeMSH() const
+{
+ return ElemChain::getTypeMSH(_dim, this->getNumVertices());
+}
+
+MElement* ElemChain::createMeshElement() const
+{
+ MElementFactory factory;
+ std::vector<MVertex*> v(_v);
+ return factory.create(this->getTypeMSH(), v);
+}
+
+int ElemChain::compareOrientation(const ElemChain& c2) const
+{
+ std::vector<MVertex*> v2;
+ c2.getMeshVertices(v2);
+
+ int perm = 1;
+ if(this->_equalVertices(v2)) return perm;
+ while(std::next_permutation(v2.begin(), v2.end(), MVertexLessThanNum())) {
+ perm *= -1;
+ if(this->_equalVertices(v2)) return perm;
+ }
+ c2.getMeshVertices(v2);
+ perm = 1;
+ while(std::prev_permutation(v2.begin(), v2.end(), MVertexLessThanNum())) {
+ perm *= -1;
+ if(this->_equalVertices(v2)) return perm;
+ }
+ return 0;
+}
+
+bool ElemChain::lessThan(const ElemChain& c2) const
+{
+ if(this->getNumSortedVertices() != c2.getNumSortedVertices())
+ return (this->getNumSortedVertices() < c2.getNumSortedVertices());
+ for(int i = 0; i < this->getNumSortedVertices(); i++){
+ if(this->getSortedVertex(i) < c2.getSortedVertex(i)) return true;
+ else if (this->getSortedVertex(i) > c2.getSortedVertex(i)) return false;
+ }
+ return false;
+}
+
+int ElemChain::getNumBoundaries(int dim, int numVertices)
+{
+ switch (dim) {
+ case 0: return 0;
+ case 1: return 2;
+ case 2:
+ switch (numVertices) {
+ case 3: return 3;
+ case 4: return 4;
+ default: return 0;
+ }
+ case 3:
+ switch (numVertices) {
+ case 4: return 4;
+ case 5: return 5;
+ case 6: return 5;
+ case 8: return 6;
+ default: return 0;
+ }
+ default: return 0;
+ }
+}
+
+int ElemChain::getNumBoundaryElemChains() const
+{
+ return ElemChain::getNumBoundaries(_dim, this->getNumVertices());
+}
+
+void ElemChain::getBoundaryVertices(int i, int dim, int numVertices,
+ const std::vector<MVertex*>& v,
+ std::vector<MVertex*>& vertices)
+{
+ vertices.clear();
+ switch (dim) {
+ case 1:
+ vertices.push_back(v[i]);
+ return;
+ case 2:
+ switch (numVertices) {
+ case 3:
+ for(int j = 0; j < 2; j++)
+ vertices.push_back(v[MTriangle::edges_tri(i, j)]);
+ return;
+ case 4:
+ for(int j = 0; j < 2; j++)
+ vertices.push_back(v[MQuadrangle::edges_quad(i, j)]);
+ return;
+ default: return;
+ }
+ case 3:
+ switch (numVertices) {
+ case 4:
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(v[MTetrahedron::faces_tetra(i, j)]);
+ return;
+ case 5:
+ if(i < 3) {
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
+ }
+ else {
+ vertices.push_back(v[0]);
+ vertices.push_back(v[3]);
+ vertices.push_back(v[2]);
+ vertices.push_back(v[1]);
+ }
+ return;
+ case 6:
+ if(i < 2)
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(v[MPrism::faces_prism(i, j)]);
+ else
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(v[MPrism::faces_prism(i, j)]);
+ return;
+ case 8:
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(v[MHexahedron::faces_hexa(i, j)]);
+ return;
+ default: return;
+ }
+ default: return;
+ }
+}
+
+ElemChain ElemChain::getBoundaryElemChain(int i) const
+{
+ std::vector<MVertex*> vertices;
+ ElemChain::getBoundaryVertices(i, _dim, this->getNumVertices(),
+ _v, vertices);
+ return ElemChain(_dim-1, vertices);
+}
+
+bool ElemChain::inEntity(GEntity* e) const
+{
+ if(_vertexCache[e].empty()) {
+ for(int i = 0; i < e->getNumMeshElements(); i++)
+ for(int j = 0; j < e->getMeshElement(i)->getNumVertices(); j++)
+ _vertexCache[e].insert(e->getMeshElement(i)->getVertex(j));
+ }
+
+ for(int i = 0; i < this->getNumVertices(); i++)
+ if(!_vertexCache[e].count(this->getMeshVertex(i))) return false;
+ return true;
+}
+
+#endif
diff --git a/Geo/Chain.h b/Geo/Chain.h
new file mode 100644
index 0000000000..88fac7d4ba
--- /dev/null
+++ b/Geo/Chain.h
@@ -0,0 +1,382 @@
+// Gmsh - Copyright (C) 1997-2012 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 <matti.pellikka@tut.fi>.
+
+#ifndef _CHAIN_H_
+#define _CHAIN_H_
+
+#include <sstream>
+#include "GModel.h"
+#include "MElement.h"
+
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "PViewOptions.h"
+#endif
+
+#if defined(HAVE_KBIPACK)
+
+template <class TTypeA, class TTypeB>
+ bool convert(const TTypeA& input, TTypeB& output ){
+ std::stringstream stream;
+ stream << input;
+ stream >> output;
+ return stream.good();
+}
+
+// Class whose derivative classes are to have partial or total order
+template <class Type>
+class PosetCat
+{
+public:
+
+ /// instantiated in derived classes
+ virtual bool lessThan(const Type& t2) const = 0;
+
+ friend bool operator<(const Type& t1, const Type& t2)
+ {
+ return t1.lessThan(t2);
+ }
+ friend bool operator>(const Type& t1, const Type& t2)
+ {
+ return !t1.lessThan(t2);
+ }
+ friend bool operator==(const Type& t1, const Type& t2) {
+ if(t1.lessThan(t2) && t2.lessThan(t1)) return true;
+ return false;
+ }
+ friend bool operator!=(const Type& t1, const Type& t2) {
+ if(t1.lessThan(t2) && t2.lessThan(t1)) return false;
+ return true;
+ }
+ friend bool operator<=(const Type& t1, const Type& t2)
+ {
+ if(t1.lessThan(t2) || t1==t2) return true;
+ return false;
+ }
+ friend bool operator>=(const Type& t1, const Type& t2)
+ {
+ if(!t1.lessThan(t2) || t1==t2) return true;
+ return false;
+ }
+
+};
+
+// Class whose derivative classes are to have vector space structure
+template <class V, class S>
+class VectorSpaceCat
+{
+public:
+
+ /// instantiated in derived classes
+ virtual V& operator+=(const V& v) = 0;
+ virtual V& operator*=(const S& s) = 0;
+ //virtual V& zero() = 0;
+ /// ---------------------
+
+ friend V operator+(const V& v1, const V& v2) {
+ V temp(v1);
+ temp += v2;
+ return temp;
+ }
+ friend V operator-(const V& v1, const V& v2) {
+ V temp(v1);
+ temp -= v2;
+ return temp;
+ }
+ friend V operator*(const S& s, const V& v) {
+ V temp(v);
+ temp*=s;
+ return temp;
+ }
+ friend V operator*(const V& v, const S& s) {
+ return s*v;
+ }
+ friend V operator/(const V& v, const S& s) {
+ S invs = 1./s;
+ return invs*v;
+ }
+
+ // Warning: assummes that the multiplicative
+ // identity element of field S can be converted from double "1."
+ // and that additive inverse of v in Abelian group V equals v*-1.
+ // (true e.g. for double and std::complex<double>),
+ // otherwise these need to be overridden by the user
+
+ virtual V& operator-() {
+ return (*this)*=(-1.);
+ }
+ virtual V& operator/=(const S& s) {
+ S temp = 1./s;
+ return (*this*=temp);
+ }
+ virtual V& operator-=(const V& v) {
+ V temp(v);
+ temp = -temp;
+ return (*this+=temp);
+ }
+};
+
+// Class to represent an 'elementary chain', a mesh cell
+class ElemChain : public PosetCat<ElemChain>
+{
+ private:
+
+ char _dim;
+ std::vector<MVertex*> _v;
+ std::vector<char> _si;
+ inline void _sortVertexIndices();
+ bool _equalVertices(const std::vector<MVertex*>& v2) const;
+
+ static std::map<GEntity*, std::set<MVertex*, MVertexLessThanNum>,
+ GEntityLessThan> _vertexCache;
+
+ public:
+
+ ElemChain(MElement* e);
+ ElemChain(int dim, std::vector<MVertex*>& v);
+
+ int getDim() const { return _dim; }
+ int getNumVertices() const { return _v.size(); }
+ MVertex* getMeshVertex(int i) const { return _v.at(i); }
+ void getMeshVertices(std::vector<MVertex*>& v) const { v = _v; }
+ int getNumSortedVertices() const { return _v.size(); }
+ inline int getSortedVertex(int i) const;
+
+ int getTypeMSH() const;
+ MElement* createMeshElement() const;
+
+ int compareOrientation(const ElemChain& c2) const;
+ bool lessThan(const ElemChain& c2) const;
+
+ int getNumBoundaryElemChains() const;
+ ElemChain getBoundaryElemChain(int i) const;
+
+ bool inEntity(GEntity* e) const;
+
+ static void clearVertexCache() { _vertexCache.clear(); }
+ static int getTypeMSH(int dim, int numVertices);
+ static int getNumBoundaries(int dim, int numVertices);
+ static void getBoundaryVertices(int i, int dim, int numVertices,
+ const std::vector<MVertex*>& v,
+ std::vector<MVertex*>& vertices);
+
+};
+
+// Class to represent a chain, formal sum of elementary chains
+template <class C>
+class Chain : public VectorSpaceCat<Chain<C>, C>
+{
+private:
+ int _dim;
+ std::map<ElemChain, C> _elemChains;
+ std::string _name;
+
+public:
+ typedef typename std::map<ElemChain, C>::iterator ecit;
+ typedef typename std::map<ElemChain, C>::const_iterator cecit;
+
+ Chain() : _dim(-1), _name("") {}
+ std::string getName() const { return _name; }
+ void setName(std::string name) { _name = name; }
+ int getDim() const { return _dim; }
+ bool isZero() const { return _elemChains.empty(); }
+ cecit firstElemChain() const { return _elemChains.begin(); }
+ cecit lastElemChain() const { return _elemChains.end(); }
+
+ void addMeshElement(MElement* e, C coeff=1);
+ void addElemChain(const ElemChain& c, C coeff=1);
+
+ Chain<C>& operator+=(const Chain<C>& chain);
+ Chain<C>& operator*=(const C& coeff);
+
+ C getCoefficient(const ElemChain& c2) const;
+ Chain<C> getBoundary() const;
+ Chain<C> getTrace(const std::vector<GEntity*>& subdomain) const;
+
+ void addToModel(GModel* m, bool post=true) const;
+};
+
+template <class C>
+C Chain<C>::getCoefficient(const ElemChain& c2) const
+{
+ cecit it = _elemChains.find(c2);
+ if(it != _elemChains.end())
+ return it->second*c2.compareOrientation(it->first);
+ else return 0;
+}
+
+template <class C>
+C incidence(const Chain<C>& c1, const Chain<C>& c2)
+{
+ C incidence = 0;
+ if(c1.getDim() != c2.getDim()) return incidence;
+ for(typename Chain<C>::cecit it = c1.firstElemChain();
+ it != c1.lastElemChain(); it++) {
+ incidence += it->second*c2.getCoefficient(it->first);
+ }
+ if(incidence != 0) {
+ Msg::Debug("%d-chains \'%s\' and \'%s\' have incidence %d", c1.getDim(),
+ c1.getName().c_str(), c2.getName().c_str(), incidence);
+
+ }
+ return incidence;
+}
+
+template <class C>
+Chain<C> boundary(const ElemChain& c)
+{
+ Chain<C> result;
+ for(int i = 0; i < c.getNumBoundaryElemChains(); i++) {
+ C coeff = 1;
+ if(c.getDim() == 1 && i == 0) coeff = -1;
+ result.addElemChain(c.getBoundaryElemChain(i), coeff);
+ }
+ return result;
+}
+
+template <class C>
+Chain<C> Chain<C>::getBoundary() const
+{
+ Chain<C> result;
+ for(cecit it = _elemChains.begin(); it != _elemChains.end(); it++) {
+ C coeff = it->second;
+ result += boundary<C>(it->first)*coeff;
+ }
+ if(result.isZero())
+ Msg::Info("The boundary chain is zero element in C%d", result.getDim());
+ return result;
+}
+
+template <class C>
+Chain<C> Chain<C>::getTrace(const std::vector<GEntity*>& subdomain) const
+{
+ Chain<C> result;
+ for(cecit it = _elemChains.begin(); it != _elemChains.end(); it++) {
+ bool inDomain = false;
+ for(unsigned int i = 0; i < subdomain.size(); i++) {
+ if(it->first.inEntity(subdomain.at(i))) {
+ inDomain = true;
+ break;
+ }
+ }
+ if(inDomain) result.addElemChain(it->first, it->second);
+ }
+ return result;
+}
+
+template <class C>
+void Chain<C>::addMeshElement(MElement* e, C coeff)
+{
+ ElemChain ce(e);
+ this->addElemChain(ce, coeff);
+}
+
+template <class C>
+void Chain<C>::addElemChain(const ElemChain& c, C coeff)
+{
+ if(_dim != -1 && _dim != c.getDim()) {
+ Msg::Error("Cannot add %d-elementrary chain to %d-chain",
+ c.getDim(), _dim);
+ return;
+ }
+ if(_dim == -1) _dim = c.getDim();
+ if(coeff == 0) return;
+ std::pair<ecit, bool> ii = _elemChains.insert( std::make_pair(c, coeff) );
+ if(!ii.second) {
+ ii.first->second += coeff*c.compareOrientation(ii.first->first);
+ if(ii.first->second == 0) _elemChains.erase(ii.first);
+ }
+}
+
+template <class C>
+Chain<C>& Chain<C>::operator+=(const Chain<C>& chain)
+{
+ for(cecit it = chain.firstElemChain();
+ it != chain.lastElemChain(); it++)
+ this->addElemChain(it->first, it->second);
+ return *this;
+}
+
+template <class C>
+Chain<C>& Chain<C>::operator*=(const C& coeff)
+{
+ if(coeff == 0)
+ _elemChains.clear();
+ else
+ for(ecit it = _elemChains.begin(); it != _elemChains.end(); it++)
+ it->second *= coeff;
+ return *this;
+ }
+
+template <class C>
+void Chain<C>::addToModel(GModel* m, bool post) const
+{
+ if(this->isZero()) {
+ Msg::Info("A chain is zero element of C%d, not added to the model",
+ this->getDim());
+ return;
+ }
+ std::vector<MElement*> elements;
+ std::map<int, std::vector<double> > data;
+ int dim = this->getDim();
+
+ for(cecit it = this->firstElemChain(); it != this->lastElemChain(); it++) {
+ MElement* e = it->first.createMeshElement();
+ C coeff = it->second;
+ elements.push_back(e);
+ if(dim > 0 && coeff < 0) e->revert();
+
+ // if elementary chain coefficient is other than -1 or 1,
+ // add multiple identical MElements to the physical group
+ for(int i = 1; i < abs(coeff); i++) {
+ MElement* ecopy = it->first.createMeshElement();
+ if(dim > 0 && coeff < 0) ecopy->revert();
+ elements.push_back(ecopy);
+ }
+
+ std::vector<double> coeffs;
+ if(dim > 0) coeffs.push_back(abs(coeff));
+ else coeffs.push_back(coeff);
+ data[e->getNum()] = coeffs;
+ }
+ int max[4];
+ for(int i = 0; i < 4; i++)
+ max[i] = m->getMaxElementaryNumber(i);
+ int entityNum = *std::max_element(max,max+4) + 1;
+ for(int i = 0; i < 4; i++)
+ max[i] = m->getMaxPhysicalNumber(i);
+ int physicalNum = *std::max_element(max,max+4) + 1;
+
+ std::map<int, std::vector<MElement*> > entityMap;
+ entityMap[entityNum] = elements;
+ std::map<int, std::map<int, std::string> > physicalMap;
+ std::map<int, std::string> physicalInfo;
+ physicalInfo[physicalNum] = _name;
+ physicalMap[entityNum] = physicalInfo;
+ m->storeChain(dim, entityMap, physicalMap);
+ m->setPhysicalName(_name, dim, physicalNum);
+
+#if defined(HAVE_POST)
+ if(post) {
+ // create PView for instant visualization
+ std::string pnum = "";
+ convert(physicalNum, pnum);
+ std::string postname = pnum + ": " + _name;
+ PView* view = new PView(postname, "ElementData", m, data, 0, 1);
+ // the user should be interested about the orientations
+ int size = 30;
+ PViewOptions* opt = view->getOptions();
+ if(opt->tangents == 0) opt->tangents = size;
+ if(opt->normals == 0) opt->normals = size;
+ view->setOptions(opt);
+ }
+#endif
+}
+
+#endif
+
+#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index a77267522e..4506f124dc 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -470,8 +470,9 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
cell->getCells(subCells);
for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
Cell* subCell = it->first;
- int coeff = it->second;
- chain[subCell] = coeff*elemi*torsion;
+ int coeff = it->second*elemi*torsion;
+ if(coeff == 0) continue;
+ chain[subCell] = coeff;
}
}
}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e21bdffaa5..ce1a2e9d45 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2844,7 +2844,8 @@ void GModel::computeHomology()
bool prepareToRestore = (itp.first != --itp.second);
itp.second++;
Homology* homology = new Homology(this, itp.first->first.first,
- itp.first->first.second, false, prepareToRestore);
+ itp.first->first.second,
+ prepareToRestore);
CellComplex *cellcomplex = homology->createCellComplex();
if(cellcomplex->getSize(0)){
for(std::multimap<dpair, std::string>::iterator itt = itp.first;
@@ -2859,6 +2860,10 @@ void GModel::computeHomology()
else
Msg::Error("Unknown type of homology computation: %s", type.c_str());
}
+ // do not save 0-, and n-chains, where n is the dimension of the model
+ // (usually not needed for anything, available through the plugin)
+ if(this->getDim() != 1) homology->addChainsToModel(1);
+ if(this->getDim() != 2) homology->addChainsToModel(2);
_pruneMeshVertexAssociations();
}
delete cellcomplex;
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index fe390b4e07..8b37fc422c 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -14,20 +14,12 @@
#if defined(HAVE_KBIPACK)
-template <class TTypeA, class TTypeB>
- bool convert(const TTypeA& input, TTypeB& output ){
- std::stringstream stream;
- stream << input;
- stream >> output;
- return stream.good();
-}
-
Homology::Homology(GModel* model, std::vector<int> physicalDomain,
std::vector<int> physicalSubdomain,
- bool save0N, bool saveOrig,
+ bool saveOrig,
bool combine, bool omit, bool smoothen) :
_model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
- _save0N(save0N), _saveOrig(saveOrig),
+ _saveOrig(saveOrig),
_combine(combine), _omit(omit), _smoothen(smoothen)
{
_fileName = "";
@@ -119,7 +111,8 @@ CellComplex* Homology::createCellComplex()
Homology::~Homology()
{
-
+ for(unsigned int i = 0; i < _chains.size(); i++)
+ delete _chains.at(i);
}
void Homology::findHomologyBasis(CellComplex* cellComplex)
@@ -165,8 +158,8 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
std::map<Cell*, int, Less_Cell> chain;
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
- if(_eraseNullCells(chain)) {
- if(!_save0N && (j != 0 && j != dim)) _createPhysicalGroup(chain, name);
+ if(!chain.empty()) {
+ _createChain(chain, name);
HRank[j] = HRank[j] + 1;
if(torsion != 1){
Msg::Warning("H%d %d has torsion coefficient %d!",
@@ -237,8 +230,8 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
std::map<Cell*, int, Less_Cell> chain;
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
- if(_eraseNullCells(chain)) {
- if(!_save0N && (j != 0 && j != dim)) _createPhysicalGroup(chain, name);
+ if(!chain.empty()) {
+ _createChain(chain, name);
HRank[j] = HRank[j] + 1;
if(torsion != 1){
Msg::Warning("H%d* %d has torsion coefficient %d!", j, i, torsion);
@@ -262,6 +255,33 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
HRank[0], HRank[1], HRank[2], HRank[3]);
}
+void Homology::addChainsToModel(int dim, bool post)
+{
+ for(unsigned int i = 0; i < _chains.size(); i++) {
+ if(dim != -1 && _chains.at(i)->getDim() == dim)
+ _chains.at(i)->addToModel(this->getModel(), post);
+ else if(dim == -1) _chains.at(i)->addToModel(this->getModel(), post);
+ }
+}
+
+void Homology::_createChain(std::map<Cell*, int, Less_Cell>& preChain,
+ std::string name)
+{
+ Chain<int>* chain = new Chain<int>();
+ chain->setName(name);
+
+ for(citer cit = preChain.begin(); cit != preChain.end(); cit++){
+ Cell* cell = cit->first;
+ int coeff = cit->second;
+ if(coeff == 0) continue;
+
+ std::vector<MVertex*> v;
+ cell->getMeshVertices(v);
+ chain->addElemChain(ElemChain(cell->getDim(), v), coeff);
+ }
+ _chains.push_back(chain);
+}
+
std::string Homology::_getDomainString(const std::vector<int>& domain,
const std::vector<int>& subdomain)
{
@@ -295,85 +315,6 @@ std::string Homology::_getDomainString(const std::vector<int>& domain,
return domainString;
}
-int Homology::_eraseNullCells(std::map<Cell*, int, Less_Cell>& chain)
-{
- std::vector<Cell*> toRemove;
- for(citer cit = chain.begin(); cit != chain.end(); cit++){
- if(cit->second == 0) toRemove.push_back(cit->first);
- }
- for(unsigned int i = 0; i < toRemove.size(); i++) chain.erase(toRemove[i]);
- return chain.size();
-}
-
-void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std::string name)
-{
- std::vector<MElement*> elements;
- std::map<int, std::vector<double> > data;
- MElementFactory factory;
- int dim = 0;
-
- typedef std::map<Cell*, int, Less_Cell>::iterator citer;
- for(citer cit = chain.begin(); cit != chain.end(); cit++){
- Cell* cell = cit->first;
- int coeff = cit->second;
- if(coeff == 0) continue;
-
- std::vector<MVertex*> v;
- cell->getMeshVertices(v);
- dim = cell->getDim();
-
- MElement* e = factory.create(cell->getTypeMSH(), v);
- if(cell->getDim() > 0 && coeff < 0) e->revert(); // flip orientation
- elements.push_back(e);
-
- // if cell coefficient is other than -1 or 1, add multiple
- // identical MElements to the physical group
- for(int i = 1; i < abs(coeff); i++) {
- MElement* ecopy = factory.create(cell->getTypeMSH(), v);
- if(cell->getDim() > 0 && coeff < 0) ecopy->revert();
- elements.push_back(ecopy);
- }
-
- std::vector<double> coeffs (1,abs(coeff));
- data[e->getNum()] = coeffs;
- }
-
- GModel* m = this->getModel();
- int max[4];
- for(int i = 0; i < 4; i++)
- max[i] = m->getMaxElementaryNumber(i);
- int entityNum = *std::max_element(max,max+4) + 1;
- for(int i = 0; i < 4; i++)
- max[i] = m->getMaxPhysicalNumber(i);
- int physicalNum = *std::max_element(max,max+4) + 1;
-
- std::map<int, std::vector<MElement*> > entityMap;
- entityMap[entityNum] = elements;
- std::map<int, std::map<int, std::string> > physicalMap;
- std::map<int, std::string> physicalInfo;
- physicalInfo[physicalNum] = name;
- physicalMap[entityNum] = physicalInfo;
-
- if(!data.empty()){
- m->storeChain(dim, entityMap, physicalMap);
- m->setPhysicalName(name, dim, physicalNum);
-
-#if defined(HAVE_POST)
- // create PView for instant visualization
- std::string pnum = "";
- convert(physicalNum, pnum);
- std::string postname = pnum + ": " + name;
- PView* view = new PView(postname, "ElementData", m, data, 0, 1);
- // the user should be interested about the orientations
- int size = 30;
- PViewOptions* opt = view->getOptions();
- if(opt->tangents == 0) opt->tangents = size;
- if(opt->normals == 0) opt->normals = size;
- view->setOptions(opt);
-#endif
- }
-}
-
bool Homology::writeBasisMSH(bool binary)
{
if(_fileName.empty()) return false;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index c5dbdc9e51..6f396fe20f 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -11,6 +11,7 @@
#include <sstream>
#include "CellComplex.h"
#include "ChainComplex.h"
+#include "Chain.h"
#include "OS.h"
#include "GModel.h"
#include "Options.h"
@@ -40,31 +41,31 @@ class Homology
// use chain smoothning
bool _smoothen;
- // save chains of 0 and highest dimension
- bool _save0N;
// save original cell complex
bool _saveOrig;
// file name to store the results
std::string _fileName;
+ // resulting chains
+ std::vector<Chain<int>*> _chains;
+
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
// create a string describing the generator
std::string _getDomainString(const std::vector<int>& domain,
const std::vector<int>& subdomain);
- // remove cells with coefficient 0
- int _eraseNullCells(std::map<Cell*, int, Less_Cell>& chain);
+ // create and store a chain from homology solver result
+ void _createChain(std::map<Cell*, int, Less_Cell>& preChain,
+ std::string name);
- // create a Gmsh physical group from a chain
- void _createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std::string name);
public:
Homology(GModel* model, std::vector<int> physicalDomain,
std::vector<int> physicalSubdomain,
- bool save0N=false, bool saveOrig=true,
+ bool saveOrig=true,
bool combine=true, bool omit=true, bool smoothen=true);
~Homology();
@@ -80,6 +81,11 @@ class Homology
void findHomologyBasis(CellComplex* cellComplex=NULL);
void findCohomologyBasis(CellComplex* cellComplex=NULL);
+ // add chains to Gmsh model
+ // dim: only add dim-chains if dim != -1
+ // post: create a post-processing view
+ void addChainsToModel(int dim=-1, bool post=true);
+
// write the generators to a file
bool writeBasisMSH(bool binary=false);
--
GitLab