1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle 2 // 3 // See the LICENSE.txt file in the Gmsh root directory for license information. 4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues. 5 // 6 // Contributed by Matti Pellikka <matti.pellikka@gmail.com>. 7 8 #ifndef HOMOLOGY_H 9 #define HOMOLOGY_H 10 11 #include <sstream> 12 #include "CellComplex.h" 13 #include "ChainComplex.h" 14 #include "Chain.h" 15 #include "OS.h" 16 #include "GModel.h" 17 #include "Options.h" 18 19 #if defined(HAVE_KBIPACK) 20 21 std::vector<int> vecN0(int n); 22 23 // Interface class for homology computation in Gmsh 24 class Homology { 25 private: 26 // the Gmsh model for homology computation 27 GModel *_model; 28 29 // domain and the relative subdomain of the homology computation 30 // physical group IDs 31 std::vector<int> _domain; 32 std::vector<int> _subdomain; 33 std::vector<int> _nondomain; 34 std::vector<int> _nonsubdomain; 35 std::vector<int> _imdomain; 36 // corresponding geometrical entities 37 std::vector<GEntity *> _domainEntities; 38 std::vector<GEntity *> _subdomainEntities; 39 std::vector<GEntity *> _nondomainEntities; 40 std::vector<GEntity *> _nonsubdomainEntities; 41 std::vector<GEntity *> _immuneEntities; 42 43 // save original cell complex 44 bool _saveOrig; 45 46 // use cell combining 47 int _combine; 48 // use cell omit 49 bool _omit; 50 // use chain smoothning 51 bool _smoothen; 52 // corecution heuristic 53 int _heuristic; 54 55 // file name to store the results 56 std::string _fileName; 57 58 // cell complex of the domain 59 CellComplex *_cellComplex; 60 61 // whether representatives of (co)homology bases are available 62 bool _homologyComputed[4]; 63 bool _cohomologyComputed[4]; 64 65 // resulting betti numbers and chains 66 int _betti[4]; 67 std::vector<Chain<int> *> _chains[4]; 68 std::vector<Chain<int> *> _cochains[4]; 69 70 typedef std::map<Cell *, int, CellPtrLessThan>::iterator citer; 71 72 void _getEntities(const std::vector<int> &physicalGroups, 73 std::vector<GEntity *> &entities); 74 void _getElements(const std::vector<GEntity *> &entities, 75 std::vector<MElement *> &elements); 76 77 // create a string describing the generator 78 std::string _getDomainString(const std::vector<int> &domain, 79 const std::vector<int> &subdomain) const; 80 81 // construct the cell complex 82 void _createCellComplex(); 83 84 // create and store a chain from homology solver result 85 void _createChain(std::map<Cell *, int, CellPtrLessThan> &preChain, 86 const std::string &name, bool co); 87 88 void _deleteChains(std::vector<int> dim = vecN0(4)); 89 void _deleteCochains(std::vector<int> dim = vecN0(4)); 90 91 std::vector<int> _addToModel(int dim, bool co, bool post, 92 int physicalNumRequest) const; 93 94 public: 95 // Determine domain and relative subdomain of (co)homology computation 96 Homology(GModel *model, const std::vector<int> &physicalDomain, 97 const std::vector<int> &physicalSubdomain, 98 const std::vector<int> &physicalIm, bool saveOrig = true, 99 int combine = 3, bool omit = true, bool smoothen = true, 100 int heuristic = 1); 101 ~Homology(); 102 getModel()103 GModel *getModel() const { return _model; } setFileName(const std::string & fileName)104 void setFileName(const std::string &fileName) { _fileName = fileName; } 105 getDomain(std::vector<int> & domain)106 void getDomain(std::vector<int> &domain) const { domain = _domain; } getSubdomain(std::vector<int> & subdomain)107 void getSubdomain(std::vector<int> &subdomain) const 108 { 109 subdomain = _subdomain; 110 } 111 112 // find the bases of (co)homology spaces 113 // if dim is not provided,, find 0-,1-,2-,3-(co)homology spaces bases 114 // otherwise only find those indicated in dim 115 void findHomologyBasis(std::vector<int> dim = vecN0(4)); 116 void findCohomologyBasis(std::vector<int> dim = vecN0(4)); 117 118 // find a homology and cohomology basis pair such that 119 // the incidence matrix of the bases is an identity matrix 120 // if master==0, homology basis determines the cohomology basis 121 // if dim is not provided, find 0-,1-,2-,3-(co)homology spaces bases 122 // otherwise only find those indicated in dim 123 void findCompatibleBasisPair(int master = 0, std::vector<int> dim = vecN0(4)); 124 125 // is the (co)homology in given dimensions already compited 126 // if dim is not provided, return true only if computed in all dimensions 127 bool isHomologyComputed(std::vector<int> dim = vecN0(4)) const; 128 bool isCohomologyComputed(std::vector<int> dim = vecN0(4)) const; 129 130 // add chains to Gmsh model 131 // dim: only add dim-chains if dim != -1 132 // post: create a post-processing view 133 // physicalNumRequest: number the chains starting from this if available 134 // returns physical group numbers 135 std::vector<int> addChainsToModel(int dim = -1, bool post = true, 136 int physicalNumRequest = -1) const; 137 std::vector<int> addCochainsToModel(int dim = -1, bool post = true, 138 int physicalNumRequest = -1) const; 139 140 // get representative chains of a (co)homology space basis 141 void getHomologyBasis(int dim, std::vector<Chain<int> > &hom); 142 void getCohomologyBasis(int dim, std::vector<Chain<int> > &coh); 143 144 // find Betti numbers 145 // faster than finding bases for (co)homology spaces 146 void findBettiNumbers(); 147 148 // are Betti numbers computed 149 bool isBettiComputed() const; 150 151 // get a Betti number 152 int betti(int dim); 153 154 // get the Euler characteristic 155 int eulerCharacteristic(); 156 157 // write the generators to a file 158 bool writeBasisMSH(bool binary = false); 159 160 // store dim-dimensional cells of cellComplex as a physical group 161 // in _model, for debugging 162 void storeCells(CellComplex *cellComplex, int dim); 163 }; 164 165 #endif 166 167 #endif 168