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 #include "Homology.h"
9 #include "fullMatrix.h"
10
11 #if defined(HAVE_KBIPACK)
12
Homology(GModel * model,const std::vector<int> & physicalDomain,const std::vector<int> & physicalSubdomain,const std::vector<int> & physicalImdomain,bool saveOrig,int combine,bool omit,bool smoothen,int heuristic)13 Homology::Homology(GModel *model, const std::vector<int> &physicalDomain,
14 const std::vector<int> &physicalSubdomain,
15 const std::vector<int> &physicalImdomain, bool saveOrig,
16 int combine, bool omit, bool smoothen, int heuristic)
17 : _model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
18 _imdomain(physicalImdomain), _saveOrig(saveOrig), _combine(combine),
19 _omit(omit), _smoothen(smoothen), _heuristic(heuristic),
20 _cellComplex(nullptr)
21 {
22 _fileName = "";
23
24 // default to the whole model
25 if(_domain.empty()) {
26 int dim = _model->getDim();
27 std::vector<GEntity *> entities;
28 _model->getEntities(entities);
29 for(auto it = entities.begin(); it != entities.end(); it++) {
30 if((*it)->dim() == dim) _domainEntities.push_back(*it);
31 }
32 }
33 else {
34 _getEntities(_domain, _domainEntities);
35 _getEntities(_subdomain, _subdomainEntities);
36 _getEntities(_nondomain, _nondomainEntities);
37 _getEntities(_nonsubdomain, _nonsubdomainEntities);
38 _getEntities(_imdomain, _immuneEntities);
39 }
40
41 for(std::size_t i = 0; i < 4; i++) {
42 _homologyComputed[i] = false;
43 _cohomologyComputed[i] = false;
44 _betti[i] = -1;
45 }
46
47 if(abs(_heuristic) > 1) _heuristic = 0;
48 }
49
vecN0(int n)50 std::vector<int> vecN0(int n)
51 {
52 std::vector<int> v;
53 v.reserve(n);
54 for(int i = 0; i < n; i++) v.push_back(i);
55 return v;
56 }
57
_getEntities(const std::vector<int> & physicalGroups,std::vector<GEntity * > & entities)58 void Homology::_getEntities(const std::vector<int> &physicalGroups,
59 std::vector<GEntity *> &entities)
60 {
61 entities.clear();
62 std::map<int, std::vector<GEntity *> > groups[4];
63 _model->getPhysicalGroups(groups);
64
65 for(std::size_t i = 0; i < physicalGroups.size(); i++) {
66 for(int j = 0; j < 4; j++) {
67 auto it = groups[j].find(physicalGroups.at(i));
68 if(it != groups[j].end()) {
69 std::vector<GEntity *> physicalGroup = (*it).second;
70 for(std::size_t k = 0; k < physicalGroup.size(); k++) {
71 entities.push_back(physicalGroup.at(k));
72 }
73 }
74 }
75 }
76 }
77
_getElements(const std::vector<GEntity * > & entities,std::vector<MElement * > & elements)78 void Homology::_getElements(const std::vector<GEntity *> &entities,
79 std::vector<MElement *> &elements)
80 {
81 elements.clear();
82 for(std::size_t j = 0; j < entities.size(); j++) {
83 for(std::size_t i = 0; i < entities.at(j)->getNumMeshElements(); i++) {
84 MElement *element = entities.at(j)->getMeshElement(i);
85 elements.push_back(element);
86 }
87 }
88 }
89
_createCellComplex()90 void Homology::_createCellComplex()
91 {
92 Msg::StatusBar(true, "Creating cell complex...");
93 double t1 = Cpu(), w1 = TimeOfDay();
94
95 if(_domainEntities.empty()) Msg::Error("Domain is empty");
96 if(_subdomainEntities.empty()) Msg::Info("Subdomain is empty");
97
98 std::vector<MElement *> domainElements;
99 std::vector<MElement *> subdomainElements;
100 std::vector<MElement *> nondomainElements;
101 std::vector<MElement *> nonsubdomainElements;
102 std::vector<MElement *> immuneElements;
103
104 _getElements(_domainEntities, domainElements);
105 _getElements(_subdomainEntities, subdomainElements);
106 _getElements(_nondomainEntities, nondomainElements);
107 _getElements(_nonsubdomainEntities, nonsubdomainElements);
108 _getElements(_immuneEntities, immuneElements);
109
110 if(_cellComplex != nullptr) delete _cellComplex;
111 _cellComplex = new CellComplex(_model, domainElements, subdomainElements,
112 nondomainElements, nonsubdomainElements,
113 immuneElements, _saveOrig);
114
115 if(_cellComplex->getSize(0) == 0) {
116 Msg::Error("Cell Complex is empty: check the domain and the mesh");
117 }
118 double t2 = Cpu(), w2 = TimeOfDay();
119 Msg::StatusBar(true, "Done creating cell complex (Wall %gs, CPU %gs)",
120 w2 - w1, t2 - t1);
121 Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
122 _cellComplex->getSize(3), _cellComplex->getSize(2),
123 _cellComplex->getSize(1), _cellComplex->getSize(0));
124 }
125
_deleteChains(std::vector<int> dim)126 void Homology::_deleteChains(std::vector<int> dim)
127 {
128 for(std::size_t j = 0; j < dim.size(); j++) {
129 int d = dim.at(j);
130 if(d < 0 || d > 3) continue;
131 for(std::size_t i = 0; i < _chains[d].size(); i++) {
132 delete _chains[d].at(i);
133 }
134 _chains[d].clear();
135 _homologyComputed[d] = false;
136 }
137 }
138
_deleteCochains(std::vector<int> dim)139 void Homology::_deleteCochains(std::vector<int> dim)
140 {
141 for(std::size_t j = 0; j < dim.size(); j++) {
142 int d = dim.at(j);
143 if(d < 0 || d > 3) continue;
144 for(std::size_t i = 0; i < _cochains[d].size(); i++) {
145 delete _cochains[d].at(i);
146 }
147 _cochains[d].clear();
148 _cohomologyComputed[d] = false;
149 }
150 }
151
~Homology()152 Homology::~Homology()
153 {
154 if(_cellComplex != nullptr) delete _cellComplex;
155 _deleteChains();
156 _deleteCochains();
157 }
158
findHomologyBasis(std::vector<int> dim)159 void Homology::findHomologyBasis(std::vector<int> dim)
160 {
161 double t0 = Cpu(), w0 = TimeOfDay();
162 std::string domain = _getDomainString(_domain, _subdomain);
163 Msg::Info("");
164 Msg::Info("To compute domain (%s) homology spaces", domain.c_str());
165
166 if(dim.empty()) {
167 findBettiNumbers();
168 return;
169 }
170
171 if(_cellComplex == nullptr) _createCellComplex();
172 if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
173
174 Msg::StatusBar(true, "Reducing cell complex...");
175
176 double t1 = Cpu(), w1 = TimeOfDay();
177 double size1 = _cellComplex->getSize(-1);
178 _cellComplex->reduceComplex(_combine, _omit);
179
180 if(_combine > 1 && !_smoothen) {
181 for(int i = 1; i <= 3; i++) {
182 if(!std::binary_search(dim.begin(), dim.end(), i)) {
183 _cellComplex->cocombine(i - 1);
184 }
185 }
186 }
187
188 double t2 = Cpu(), w2 = TimeOfDay();
189 double size2 = _cellComplex->getSize(-1);
190 Msg::StatusBar(true, "Done reducing cell complex (Wall %gs, CPU %gs, %g%%)",
191 w2 - w1, t2 - t1, (1. - size2 / size1) * 100.);
192 Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
193 _cellComplex->getSize(3), _cellComplex->getSize(2),
194 _cellComplex->getSize(1), _cellComplex->getSize(0));
195
196 Msg::StatusBar(true, "Computing homology space bases...");
197 t1 = Cpu();
198 w1 = TimeOfDay();
199 ChainComplex chainComplex = ChainComplex(_cellComplex);
200 chainComplex.computeHomology();
201 t2 = Cpu();
202 w2 = TimeOfDay();
203 Msg::StatusBar(true,
204 "Done computing homology space bases (Wall %gs, CPU %gs)",
205 w2 - w1, t2 - t1);
206
207 _deleteChains(dim);
208 for(int j = 0; j < 4; j++) {
209 _betti[j] = 0;
210 std::string dimension = convertInt(j);
211 for(int i = 1; i <= chainComplex.getBasisSize(j, 3); i++) {
212 std::string generator = convertInt(i);
213
214 std::string name = "H_" + dimension + domain + generator;
215 std::map<Cell *, int, CellPtrLessThan> chain;
216 chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
217 int torsion = chainComplex.getTorsion(j, i);
218 if(!chain.empty()) {
219 _createChain(chain, name, false);
220 _betti[j] = _betti[j] + 1;
221 if(torsion != 1) {
222 Msg::Warning("H_%d %d has torsion coefficient %d!", j, i, torsion);
223 }
224 }
225 }
226 }
227
228 if(_fileName != "") writeBasisMSH();
229
230 Msg::Info("Ranks of domain (%s) homology spaces:", domain.c_str());
231 Msg::Info("H_0 = %d", _betti[0]);
232 Msg::Info("H_1 = %d", _betti[1]);
233 Msg::Info("H_2 = %d", _betti[2]);
234 Msg::Info("H_3 = %d", _betti[3]);
235
236 double t3 = Cpu(), w3 = TimeOfDay();
237 Msg::Info("Done computing (%s) homology spaces (Wall %gs, CPU %gs)",
238 domain.c_str(), w3 - w0, t3 - t0);
239 Msg::StatusBar(false, "H_0: %d, H_1: %d, H_2: %d, H_3: %d", _betti[0],
240 _betti[1], _betti[2], _betti[3]);
241
242 for(std::size_t i = 0; i < dim.size(); i++) {
243 int d = dim.at(i);
244 if(d >= 0 && d < 4) _homologyComputed[d] = true;
245 }
246 }
247
findCohomologyBasis(std::vector<int> dim)248 void Homology::findCohomologyBasis(std::vector<int> dim)
249 {
250 double t0 = Cpu(), w0 = TimeOfDay();
251 std::string domain = _getDomainString(_domain, _subdomain);
252 Msg::Info("");
253 Msg::Info("To compute domain (%s) cohomology spaces", domain.c_str());
254
255 if(dim.empty()) {
256 findBettiNumbers();
257 return;
258 }
259
260 if(_cellComplex == nullptr) _createCellComplex();
261 if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
262
263 Msg::StatusBar(true, "Reducing cell complex...");
264
265 double t1 = Cpu(), w1 = TimeOfDay();
266 double size1 = _cellComplex->getSize(-1);
267
268 _cellComplex->coreduceComplex(_combine, _omit, _heuristic);
269
270 std::sort(dim.begin(), dim.end());
271 if(_combine > 1) {
272 for(int i = 2; i >= 0; i--) {
273 if(!std::binary_search(dim.begin(), dim.end(), i)) {
274 _cellComplex->combine(i + 1);
275 }
276 }
277 }
278
279 double t2 = Cpu(), w2 = TimeOfDay();
280 double size2 = _cellComplex->getSize(-1);
281
282 Msg::StatusBar(true, "Done reducing cell complex (Wall %gs, CPU %gs, %g %%)",
283 w2 - w1, t2 - t1, (1. - size2 / size1) * 100.);
284 Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
285 _cellComplex->getSize(3), _cellComplex->getSize(2),
286 _cellComplex->getSize(1), _cellComplex->getSize(0));
287
288 Msg::StatusBar(true, "Computing cohomology space bases ...");
289 t1 = Cpu();
290 w1 = TimeOfDay();
291 ChainComplex chainComplex = ChainComplex(_cellComplex);
292 chainComplex.computeHomology(true);
293 t2 = Cpu();
294 w2 = TimeOfDay();
295 Msg::StatusBar(true,
296 "Done computing cohomology space bases (Wall %gs, CPU %gs)",
297 w2 - w1, t2 - t1);
298
299 _deleteCochains(dim);
300 for(int i = 0; i < 4; i++) _betti[i] = 0;
301 for(int j = 3; j > -1; j--) {
302 std::string dimension = convertInt(j);
303
304 for(int i = 1; i <= chainComplex.getBasisSize(j, 3); i++) {
305 std::string generator = convertInt(i);
306
307 std::string name = "H^" + dimension + domain + generator;
308 std::map<Cell *, int, CellPtrLessThan> chain;
309 chainComplex.getBasisChain(chain, i, j, 3, false);
310 int torsion = chainComplex.getTorsion(j, i);
311 if(!chain.empty()) {
312 _createChain(chain, name, true);
313 _betti[j] = _betti[j] + 1;
314 if(torsion != 1) {
315 Msg::Warning("H^%d %d has torsion coefficient %d!", j, i, torsion);
316 }
317 }
318 }
319 }
320
321 if(_fileName != "") writeBasisMSH();
322
323 Msg::Info("Ranks of domain (%s) cohomology spaces:", domain.c_str());
324 Msg::Info("H^0 = %d", _betti[0]);
325 Msg::Info("H^1 = %d", _betti[1]);
326 Msg::Info("H^2 = %d", _betti[2]);
327 Msg::Info("H^3 = %d", _betti[3]);
328
329 double t3 = Cpu(), w3 = TimeOfDay();
330 Msg::Info("Done computing (%s) cohomology spaces (Wall %gs, CPU %gs)",
331 domain.c_str(), w3 - w0, t3 - t0);
332 Msg::StatusBar(false, "H^0: %d, H^1: %d, H^2: %d, H^3: %d", _betti[0],
333 _betti[1], _betti[2], _betti[3]);
334
335 for(std::size_t i = 0; i < dim.size(); i++) {
336 int d = dim.at(i);
337 if(d >= 0 && d < 4) _cohomologyComputed[d] = true;
338 }
339 }
340
isBettiComputed() const341 bool Homology::isBettiComputed() const
342 {
343 bool computed = true;
344 for(int i = 0; i < 4; i++) {
345 if(_betti[i] == -1) computed = false;
346 }
347 return computed;
348 }
349
isHomologyComputed(std::vector<int> dim) const350 bool Homology::isHomologyComputed(std::vector<int> dim) const
351 {
352 bool computed = true;
353 for(std::size_t i = 0; i < dim.size(); i++) {
354 int d = dim.at(i);
355 if(d < 0 || d > 3) continue;
356 computed = computed && _homologyComputed[d];
357 }
358 return computed;
359 }
360
isCohomologyComputed(std::vector<int> dim) const361 bool Homology::isCohomologyComputed(std::vector<int> dim) const
362 {
363 bool computed = true;
364 for(std::size_t i = 0; i < dim.size(); i++) {
365 int d = dim.at(i);
366 if(d < 0 || d > 3) continue;
367 computed = computed && _cohomologyComputed[d];
368 }
369 return computed;
370 }
371
findCompatibleBasisPair(int master,std::vector<int> dim)372 void Homology::findCompatibleBasisPair(int master, std::vector<int> dim)
373 {
374 if(!this->isHomologyComputed(dim)) this->findHomologyBasis(dim);
375 if(!this->isCohomologyComputed(dim)) this->findCohomologyBasis(dim);
376 for(std::size_t idim = 0; idim < dim.size(); idim++) {
377 int d = dim.at(idim);
378 if(d < 1 || d > 2) continue;
379 int n = this->betti(d);
380 if(n < 2) continue;
381 if((int)_chains[d].size() != n || (int)_cochains[d].size() != n) {
382 Msg::Warning("Cannot produce compatible %d-(co)homology bases.", d);
383 Msg::Debug("%d basis %d-chains and %d basis %d-cochains.",
384 (int)_chains[d].size(), d, (int)_cochains[d].size(), d);
385 continue;
386 }
387 fullMatrix<double> m(n, n);
388 for(int i = 0; i < n; i++) {
389 for(int j = 0; j < n; j++) {
390 if(master == 0)
391 m(i, j) = incidence(*_cochains[d].at(i), *_chains[d].at(j));
392 else
393 m(i, j) = incidence(*_chains[d].at(i), *_cochains[d].at(j));
394 }
395 }
396
397 int det = m.determinant();
398 if(abs(det) != 1 || !m.invertInPlace()) {
399 Msg::Warning("Cannot produce compatible %d-(co)homology bases.", d);
400 Msg::Debug("Incidence matrix: ");
401 for(int i = 0; i < n; i++)
402 for(int j = 0; j < n; j++) Msg::Debug("(%d, %d) = %d", i, j, m(i, j));
403 continue;
404 }
405
406 std::vector<Chain<int> *> newBasis(n);
407
408 if(master == 0) {
409 for(int i = 0; i < n; i++) {
410 newBasis.at(i) = new Chain<int>();
411 for(int j = 0; j < n; j++) {
412 *newBasis.at(i) += (int)m(i, j) * (*_cochains[d].at(j));
413 }
414 }
415 for(int i = 0; i < n; i++) {
416 newBasis.at(i)->setName(_cochains[d].at(i)->getName());
417 delete _cochains[d].at(i);
418 _cochains[d].at(i) = newBasis.at(i);
419 }
420 }
421 else {
422 for(int i = 0; i < n; i++) {
423 newBasis.at(i) = new Chain<int>();
424 for(int j = 0; j < n; j++) {
425 *newBasis.at(i) += (int)m(i, j) * (*_chains[d].at(j));
426 }
427 }
428 for(int i = 0; i < n; i++) {
429 newBasis.at(i)->setName(_chains[d].at(i)->getName());
430 delete _chains[d].at(i);
431 _chains[d].at(i) = newBasis.at(i);
432 }
433 }
434 }
435 }
436
_addToModel(int dim,bool co,bool post,int physicalNumRequest) const437 std::vector<int> Homology::_addToModel(int dim, bool co, bool post,
438 int physicalNumRequest) const
439 {
440 std::vector<int> physicals;
441 if(dim < 0 || dim > 3) return physicals;
442 int pgnum = -1;
443 if(!co) {
444 for(std::size_t i = 0; i < _chains[dim].size(); i++) {
445 if(physicalNumRequest != -1)
446 pgnum = physicalNumRequest + i;
447 else
448 pgnum = -1;
449 physicals.push_back(
450 _chains[dim].at(i)->addToModel(this->getModel(), post, pgnum));
451 }
452 }
453 else {
454 for(std::size_t i = 0; i < _cochains[dim].size(); i++) {
455 if(physicalNumRequest != -1)
456 pgnum = physicalNumRequest + i;
457 else
458 pgnum = -1;
459 physicals.push_back(
460 _cochains[dim].at(i)->addToModel(this->getModel(), post, pgnum));
461 }
462 }
463 if(!physicals.empty()) {
464 std::vector<int> empty;
465 std::string span = _getDomainString(physicals, empty);
466 std::string domain = _getDomainString(_domain, _subdomain);
467 if(!co)
468 Msg::Info("Span H_%d(%s) = %s", dim, domain.c_str(), span.c_str());
469 else
470 Msg::Info("Span H^%d(%s) = %s", dim, domain.c_str(), span.c_str());
471 }
472 return physicals;
473 }
474
addChainsToModel(int dim,bool post,int physicalNumRequest) const475 std::vector<int> Homology::addChainsToModel(int dim, bool post,
476 int physicalNumRequest) const
477 {
478 std::vector<int> physicals;
479 if(dim > -1 && !_homologyComputed[dim])
480 Msg::Warning("%d-Homology is not computed", dim);
481 if(dim == -1) {
482 for(int j = 0; j < 4; j++) {
483 std::vector<int> p = _addToModel(j, false, post, physicalNumRequest);
484 physicals.insert(physicals.end(), p.begin(), p.end());
485 }
486 }
487 else if(dim > -1 && dim < 4) {
488 physicals = _addToModel(dim, false, post, physicalNumRequest);
489 }
490 return physicals;
491 }
492
addCochainsToModel(int dim,bool post,int physicalNumRequest) const493 std::vector<int> Homology::addCochainsToModel(int dim, bool post,
494 int physicalNumRequest) const
495 {
496 std::vector<int> physicals;
497 if(dim > -1 && !_cohomologyComputed[dim])
498 Msg::Warning("%d-Cohomology is not computed", dim);
499 if(dim == -1) {
500 for(int j = 0; j < 4; j++) {
501 std::vector<int> p = _addToModel(j, true, post, physicalNumRequest);
502 physicals.insert(physicals.end(), p.begin(), p.end());
503 }
504 }
505 else if(dim > -1 && dim < 4) {
506 physicals = _addToModel(dim, true, post, physicalNumRequest);
507 }
508 return physicals;
509 }
510
getHomologyBasis(int dim,std::vector<Chain<int>> & hom)511 void Homology::getHomologyBasis(int dim, std::vector<Chain<int> > &hom)
512 {
513 if(dim < 0 || dim > 3) return;
514 if(!_homologyComputed[dim]) findHomologyBasis();
515
516 hom.resize(_chains[dim].size());
517 for(std::size_t i = 0; i < _chains[dim].size(); i++)
518 hom[i] = *_chains[dim].at(i);
519 }
520
getCohomologyBasis(int dim,std::vector<Chain<int>> & coh)521 void Homology::getCohomologyBasis(int dim, std::vector<Chain<int> > &coh)
522 {
523 if(dim < 0 || dim > 3) return;
524 if(!_cohomologyComputed[dim]) findCohomologyBasis();
525
526 coh.resize(_cochains[dim].size());
527 for(std::size_t i = 0; i < _cochains[dim].size(); i++)
528 coh[i] = *_cochains[dim].at(i);
529 }
530
findBettiNumbers()531 void Homology::findBettiNumbers()
532 {
533 if(!isBettiComputed()) {
534 if(_cellComplex == nullptr) _createCellComplex();
535 if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
536
537 Msg::StatusBar(true, "Reducing cell complex...");
538 double t1 = Cpu(), w1 = TimeOfDay();
539 double size1 = _cellComplex->getSize(-1);
540
541 _cellComplex->bettiReduceComplex();
542
543 double t2 = Cpu(), w2 = TimeOfDay();
544 double size2 = _cellComplex->getSize(-1);
545
546 Msg::StatusBar(true,
547 "Done reducing cell complex (Wall %gs, CPU %gs, %g %%)",
548 w2 - w1, t2 - t1, (1. - size2 / size1) * 100.);
549 Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
550 _cellComplex->getSize(3), _cellComplex->getSize(2),
551 _cellComplex->getSize(1), _cellComplex->getSize(0));
552
553 Msg::StatusBar(true, "Computing betti numbers...");
554 t1 = Cpu();
555 w1 = TimeOfDay();
556 ChainComplex chainComplex = ChainComplex(_cellComplex);
557 chainComplex.computeHomology();
558
559 for(int i = 0; i < 4; i++) _betti[i] = chainComplex.getBasisSize(i, 3);
560
561 t2 = Cpu();
562 w2 = TimeOfDay();
563 Msg::StatusBar(true, "Betti numbers computed (Wall %gs, CPU %gs)", w2 - w1,
564 t2 - t1);
565 }
566
567 std::string domain = _getDomainString(_domain, _subdomain);
568 Msg::Info("Domain %s Betti numbers:", domain.c_str());
569 Msg::Info("b0 = %d", _betti[0]);
570 Msg::Info("b1 = %d", _betti[1]);
571 Msg::Info("b2 = %d", _betti[2]);
572 Msg::Info("b3 = %d", _betti[3]);
573
574 Msg::StatusBar(false, "b0: %d, b1: %d, b2: %d, b3: %d", _betti[0], _betti[1],
575 _betti[2], _betti[3]);
576 }
577
betti(int dim)578 int Homology::betti(int dim)
579 {
580 if(dim < 0 || dim > 3) return 0;
581 if(_betti[dim] != -1) return _betti[dim];
582
583 findBettiNumbers();
584 return _betti[dim];
585 }
586
eulerCharacteristic()587 int Homology::eulerCharacteristic()
588 {
589 if(_cellComplex == nullptr) _createCellComplex();
590 return _cellComplex->eulerCharacteristic();
591 }
592
_createChain(std::map<Cell *,int,CellPtrLessThan> & preChain,const std::string & name,bool co)593 void Homology::_createChain(std::map<Cell *, int, CellPtrLessThan> &preChain,
594 const std::string &name, bool co)
595 {
596 Chain<int> *chain = new Chain<int>();
597 chain->setName(name);
598
599 for(auto cit = preChain.begin(); cit != preChain.end(); cit++) {
600 Cell *cell = cit->first;
601 int coeff = cit->second;
602 if(coeff == 0) continue;
603
604 std::vector<MVertex *> v;
605 cell->getMeshVertices(v);
606 chain->addElemChain(ElemChain(cell->getDim(), v), coeff);
607 }
608 if(co)
609 _cochains[chain->getDim()].push_back(chain);
610 else
611 _chains[chain->getDim()].push_back(chain);
612 }
613
_getDomainString(const std::vector<int> & domain,const std::vector<int> & subdomain) const614 std::string Homology::_getDomainString(const std::vector<int> &domain,
615 const std::vector<int> &subdomain) const
616 {
617 std::string domainString = "{";
618 if(domain.empty())
619 domainString += "0";
620 else {
621 for(std::size_t i = 0; i < domain.size(); i++) {
622 std::string temp = convertInt(domain.at(i));
623 domainString += temp;
624 if(domain.size() - 1 > i) { domainString += ","; }
625 }
626 }
627 domainString += "}";
628
629 if(!subdomain.empty()) {
630 domainString += ",{";
631 for(std::size_t i = 0; i < subdomain.size(); i++) {
632 std::string temp = convertInt(subdomain.at(i));
633 domainString += temp;
634 if(subdomain.size() - 1 > i) { domainString += ","; }
635 }
636 domainString += "}";
637 }
638
639 return domainString;
640 }
641
writeBasisMSH(bool binary)642 bool Homology::writeBasisMSH(bool binary)
643 {
644 if(_fileName.empty()) return false;
645 if(!_model->writeMSH(_fileName, 2.0, binary)) return false;
646 Msg::Info("Wrote homology computation results to %s", _fileName.c_str());
647 return true;
648 }
649
storeCells(CellComplex * cellComplex,int dim)650 void Homology::storeCells(CellComplex *cellComplex, int dim)
651 {
652 std::vector<MElement *> elements;
653 MElementFactory factory;
654
655 for(auto cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim);
656 cit++) {
657 Cell *cell = *cit;
658 std::map<Cell *, int, CellPtrLessThan> cells;
659 cell->getCells(cells);
660 for(auto it = cells.begin(); it != cells.end(); it++) {
661 Cell *subCell = it->first;
662 std::vector<MVertex *> v;
663 subCell->getMeshVertices(v);
664
665 MElement *e = factory.create(subCell->getTypeMSH(), v);
666 elements.push_back(e);
667 }
668 }
669
670 int max[4];
671 for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
672 int entityNum = *std::max_element(max, max + 4) + 1;
673 for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
674 int physicalNum = *std::max_element(max, max + 4) + 1;
675
676 std::map<int, std::vector<MElement *> > entityMap;
677 entityMap[entityNum] = elements;
678 std::map<int, std::map<int, std::string> > physicalMap;
679 std::map<int, std::string> physicalInfo;
680 physicalInfo[physicalNum] = "Cell Complex";
681 physicalMap[entityNum] = physicalInfo;
682
683 _model->storeChain(dim, entityMap, physicalMap);
684 _model->setPhysicalName("Cell Complex", dim, physicalNum);
685 }
686
687 #endif
688