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 #include <limits>
7 #include <stdlib.h>
8 #include <sstream>
9 #include <stack>
10 #include "GmshConfig.h"
11 #include "GmshMessage.h"
12 #include "GModel.h"
13 #include "GModelIO_GEO.h"
14 #include "GModelIO_OCC.h"
15 #include "MPoint.h"
16 #include "MLine.h"
17 #include "MTriangle.h"
18 #include "MQuadrangle.h"
19 #include "MTetrahedron.h"
20 #include "MHexahedron.h"
21 #include "MPrism.h"
22 #include "MPyramid.h"
23 #include "MTrihedron.h"
24 #include "MElementCut.h"
25 #include "MElementOctree.h"
26 #include "discreteRegion.h"
27 #include "discreteFace.h"
28 #include "discreteEdge.h"
29 #include "discreteVertex.h"
30 #include "partitionRegion.h"
31 #include "partitionFace.h"
32 #include "partitionEdge.h"
33 #include "partitionVertex.h"
34 #include "gmshSurface.h"
35 #include "SmoothData.h"
36 #include "Context.h"
37 #include "OS.h"
38 #include "StringUtils.h"
39 #include "GEdgeLoop.h"
40 #include "MVertexRTree.h"
41 #include "OpenFile.h"
42 #include "CreateFile.h"
43 #include "Options.h"
44 #include "GModelParametrize.h"
45
46 #if defined(HAVE_MESH)
47 #include "meshGEdge.h"
48 #include "meshGFace.h"
49 #include "meshGRegion.h"
50 #include "Field.h"
51 #include "Generator.h"
52 #include "meshPartition.h"
53 #include "HighOrder.h"
54 #include "meshMetric.h"
55 #include "meshGRegionMMG.h"
56 #include "meshGFaceBamg.h"
57 #include "meshRefine.h"
58 #endif
59
60 #if defined(HAVE_POST)
61 #include "PView.h"
62 #include "PViewDataGModel.h"
63 #endif
64
65 #if defined(HAVE_KBIPACK)
66 #include "Homology.h"
67 #endif
68
69 std::vector<GModel *> GModel::list;
70 int GModel::_current = -1;
71
GModel(const std::string & name)72 GModel::GModel(const std::string &name)
73 : _name(name), _visible(1), _elementOctree(nullptr),
74 _geo_internals(nullptr), _occ_internals(nullptr), _acis_internals(nullptr),
75 _parasolid_internals(nullptr), _fields(nullptr),
76 _currentMeshEntity(nullptr), _numPartitions(0), normals(nullptr),
77 lcCallback(nullptr)
78 {
79 _maxVertexNum = CTX::instance()->mesh.firstNodeTag - 1;
80 _maxElementNum = CTX::instance()->mesh.firstElementTag - 1;
81 _checkPointedMaxVertexNum = _maxVertexNum;
82 _checkPointedMaxElementNum = _maxElementNum;
83
84 // hide all other models
85 for(std::size_t i = 0; i < list.size(); i++) list[i]->setVisibility(0);
86
87 // push new one into the list
88 list.push_back(this);
89
90 // we always create an internal GEO model; other CAD internals are created
91 // on-demand
92 createGEOInternals();
93
94 #if defined(HAVE_MESH)
95 _fields = new FieldManager();
96 #endif
97 }
98
~GModel()99 GModel::~GModel()
100 {
101 auto it = std::find(list.begin(), list.end(), this);
102 if(it != list.end()) list.erase(it);
103
104 if(getVisibility()) {
105 // if no other model is visible, make the last one visible
106 bool othervisible = false;
107 for(std::size_t i = 0; i < list.size(); i++) {
108 if(list[i]->getVisibility()) othervisible = true;
109 }
110 if(!othervisible && list.size()) list.back()->setVisibility(1);
111 }
112
113 destroy();
114 deleteGEOInternals();
115 deleteOCCInternals();
116 deleteACISInternals();
117 deleteParasolidInternals();
118 #if defined(HAVE_MESH)
119 delete _fields;
120 #endif
121 }
122
setFileName(const std::string & fileName)123 void GModel::setFileName(const std::string &fileName)
124 {
125 _fileName = fileName;
126 _fileNames.insert(fileName);
127
128 Msg::SetOnelabString("Gmsh/Model name", fileName,
129 Msg::GetNumOnelabClients() > 1, false, true, 0, "file");
130 Msg::SetOnelabString("Gmsh/Model absolute path",
131 SplitFileName(GetAbsolutePath(fileName))[0], false,
132 false, true, 0);
133 Msg::SetWindowTitle(fileName);
134 }
135
current(int index)136 GModel *GModel::current(int index)
137 {
138 if(list.empty()) {
139 Msg::Debug("No current model available: creating one");
140 new GModel();
141 }
142 if(index >= 0) _current = index;
143 if(_current < 0 || _current >= (int)list.size()) return list.back();
144 return list[_current];
145 }
146
setCurrent(GModel * m)147 int GModel::setCurrent(GModel *m)
148 {
149 for(std::size_t i = 0; i < list.size(); i++) {
150 if(list[i] == m) {
151 _current = i;
152 break;
153 }
154 }
155 return _current;
156 }
157
findByName(const std::string & name,const std::string & fileName)158 GModel *GModel::findByName(const std::string &name, const std::string &fileName)
159 {
160 // return last mesh with given name
161 for(int i = list.size() - 1; i >= 0; i--)
162 if(list[i]->getName() == name &&
163 (fileName.empty() || !list[i]->hasFileName(fileName)))
164 return list[i];
165 return nullptr;
166 }
167
destroy(bool keepName)168 void GModel::destroy(bool keepName)
169 {
170 Msg::Debug("Destroying model %s", getName().c_str());
171
172 if(!keepName) {
173 _name.clear();
174 _fileNames.clear();
175 }
176
177 _maxVertexNum = CTX::instance()->mesh.firstNodeTag - 1;
178 _maxElementNum = CTX::instance()->mesh.firstElementTag - 1;
179 _checkPointedMaxVertexNum = _maxVertexNum;
180 _checkPointedMaxElementNum = _maxElementNum;
181 _currentMeshEntity = nullptr;
182 _lastMeshEntityError.clear();
183 _lastMeshVertexError.clear();
184
185 for(auto it = firstRegion(); it != lastRegion(); ++it) delete *it;
186 regions.clear();
187 std::set<GRegion *, GEntityPtrLessThan>().swap(regions);
188
189 for(auto it = firstFace(); it != lastFace(); ++it) delete *it;
190 faces.clear();
191 std::set<GFace *, GEntityPtrLessThan>().swap(faces);
192
193 for(auto it = firstEdge(); it != lastEdge(); ++it) delete *it;
194 edges.clear();
195 std::set<GEdge *, GEntityPtrLessThan>().swap(edges);
196
197 for(auto it = firstVertex(); it != lastVertex(); ++it) delete *it;
198 vertices.clear();
199 std::set<GVertex *, GEntityPtrLessThan>().swap(vertices);
200
201 destroyMeshCaches();
202
203 resetOCCInternals();
204
205 if(normals) delete normals;
206 normals = nullptr;
207
208 #if defined(HAVE_MESH)
209 _fields->reset();
210 #endif
211 gmshSurface::reset();
212 }
213
destroyMeshCaches()214 void GModel::destroyMeshCaches()
215 {
216 // this is called in GEntity::deleteMesh()
217 #pragma omp critical
218 {
219 _vertexVectorCache.clear();
220 std::vector<MVertex *>().swap(_vertexVectorCache);
221 _vertexMapCache.clear();
222 std::map<int, MVertex *>().swap(_vertexMapCache);
223 _elementVectorCache.clear();
224 std::vector<std::pair<MElement *, int> >().swap(_elementVectorCache);
225 _elementMapCache.clear();
226 std::map<int, std::pair<MElement *, int> >().swap(_elementMapCache);
227 _elementIndexCache.clear();
228 std::map<int, int>().swap(_elementIndexCache);
229 if(_elementOctree) {
230 delete _elementOctree;
231 _elementOctree = nullptr;
232 }
233 }
234 }
235
deleteMesh()236 void GModel::deleteMesh()
237 {
238 for(auto it = firstRegion(); it != lastRegion(); ++it) (*it)->deleteMesh();
239 for(auto it = firstFace(); it != lastFace(); ++it) (*it)->deleteMesh();
240 for(auto it = firstEdge(); it != lastEdge(); ++it) (*it)->deleteMesh();
241 for(auto it = firstVertex(); it != lastVertex(); ++it) (*it)->deleteMesh();
242 destroyMeshCaches();
243 _currentMeshEntity = nullptr;
244 _lastMeshEntityError.clear();
245 _lastMeshVertexError.clear();
246 }
247
deleteMesh(const std::vector<GEntity * > & entities)248 void GModel::deleteMesh(const std::vector<GEntity *> &entities)
249 {
250 if(entities.empty()) {
251 deleteMesh();
252 return;
253 }
254 for(std::size_t i = 0; i < entities.size(); i++) {
255 GEntity *ge = entities[i];
256 bool ok = true;
257 switch(ge->dim()) {
258 case 0: {
259 std::vector<GEdge *> e = ge->edges();
260 for(auto it = e.begin(); it != e.end(); ++it) {
261 if((*it)->getNumMeshElements()) {
262 ok = false;
263 break;
264 }
265 }
266 } break;
267 case 1: {
268 std::vector<GFace *> f = ge->faces();
269 for(auto it = f.begin(); it != f.end(); ++it) {
270 if((*it)->getNumMeshElements()) {
271 ok = false;
272 break;
273 }
274 }
275 } break;
276 case 2: {
277 std::list<GRegion *> r = ge->regions();
278 for(auto it = r.begin(); it != r.end(); ++it) {
279 if((*it)->getNumMeshElements()) {
280 ok = false;
281 break;
282 }
283 }
284 } break;
285 }
286 if(ok) { ge->deleteMesh(); }
287 else {
288 Msg::Warning("Cannot delete mesh of entity (%d, %d), connected to mesh "
289 "of higher dimensional entity",
290 ge->dim(), ge->tag());
291 }
292 }
293 destroyMeshCaches();
294 _currentMeshEntity = nullptr;
295 _lastMeshEntityError.clear();
296 _lastMeshVertexError.clear();
297 }
298
deleteVertexArrays()299 void GModel::deleteVertexArrays()
300 {
301 for(auto it = firstRegion(); it != lastRegion(); ++it)
302 (*it)->deleteVertexArrays();
303 for(auto it = firstFace(); it != lastFace(); ++it)
304 (*it)->deleteVertexArrays();
305 for(auto it = firstEdge(); it != lastEdge(); ++it)
306 (*it)->deleteVertexArrays();
307 for(auto it = firstVertex(); it != lastVertex(); ++it)
308 (*it)->deleteVertexArrays();
309 }
310
empty() const311 bool GModel::empty() const
312 {
313 return vertices.empty() && edges.empty() && faces.empty() && regions.empty();
314 }
315
getRegionByTag(int n) const316 GRegion *GModel::getRegionByTag(int n) const
317 {
318 GEntity tmp((GModel *)this, n);
319 auto it = regions.find((GRegion *)&tmp);
320 if(it != regions.end())
321 return *it;
322 else
323 return nullptr;
324 }
325
getFaceByTag(int n) const326 GFace *GModel::getFaceByTag(int n) const
327 {
328 GEntity tmp((GModel *)this, n);
329 auto it = faces.find((GFace *)&tmp);
330 if(it != faces.end())
331 return *it;
332 else
333 return nullptr;
334 }
335
getEdgeByTag(int n) const336 GEdge *GModel::getEdgeByTag(int n) const
337 {
338 GEntity tmp((GModel *)this, n);
339 auto it = edges.find((GEdge *)&tmp);
340 if(it != edges.end())
341 return *it;
342 else
343 return nullptr;
344 }
345
getVertexByTag(int n) const346 GVertex *GModel::getVertexByTag(int n) const
347 {
348 GEntity tmp((GModel *)this, n);
349 auto it = vertices.find((GVertex *)&tmp);
350 if(it != vertices.end())
351 return *it;
352 else
353 return nullptr;
354 }
355
getEntityByTag(int dim,int n) const356 GEntity *GModel::getEntityByTag(int dim, int n) const
357 {
358 switch(dim) {
359 case 0: return getVertexByTag(n);
360 case 1: return getEdgeByTag(n);
361 case 2: return getFaceByTag(n);
362 case 3: return getRegionByTag(n);
363 }
364 return nullptr;
365 }
366
changeEntityTag(int dim,int tag,int newTag)367 bool GModel::changeEntityTag(int dim, int tag, int newTag)
368 {
369 if(dim == 0) {
370 GVertex *gv = getVertexByTag(tag);
371 if(gv) {
372 vertices.erase(gv);
373 gv->setTag(newTag);
374 vertices.insert(gv);
375 }
376 else {
377 Msg::Error("Unknown model point %d", tag);
378 return false;
379 }
380 }
381 else if(dim == 1) {
382 GEdge *ge = getEdgeByTag(tag);
383 if(ge) {
384 edges.erase(ge);
385 ge->setTag(newTag);
386 edges.insert(ge);
387 }
388 else {
389 Msg::Error("Unknown model curve %d", tag);
390 return false;
391 }
392 }
393 else if(dim == 2) {
394 GFace *gf = getFaceByTag(tag);
395 if(gf) {
396 faces.erase(gf);
397 gf->setTag(newTag);
398 faces.insert(gf);
399 }
400 else {
401 Msg::Error("Unknown model surface %d", tag);
402 return false;
403 }
404 }
405 else if(dim == 3) {
406 GRegion *gr = getRegionByTag(tag);
407 if(gr) {
408 regions.erase(gr);
409 gr->setTag(newTag);
410 regions.insert(gr);
411 }
412 else {
413 Msg::Error("Unknown model volume %d", tag);
414 return false;
415 }
416 }
417 return true;
418 }
419
getTagsForPhysicalName(int dim,const std::string & name)420 std::vector<int> GModel::getTagsForPhysicalName(int dim,
421 const std::string &name)
422 {
423 std::vector<int> tags;
424 std::map<int, std::vector<GEntity *> > physicalGroups;
425 getPhysicalGroups(dim, physicalGroups);
426 std::vector<GEntity *> entities =
427 physicalGroups[getPhysicalNumber(dim, name)];
428 for(auto it = entities.begin(); it != entities.end(); it++) {
429 GEntity *ge = *it;
430 tags.push_back(ge->tag());
431 }
432 return tags;
433 }
434
remove(GRegion * r)435 bool GModel::remove(GRegion *r)
436 {
437 auto it = std::find(firstRegion(), lastRegion(), r);
438 if(it != (riter)regions.end()) {
439 regions.erase(it);
440 std::vector<GFace *> f = r->faces();
441 for(auto it = f.begin(); it != f.end(); it++) (*it)->delRegion(r);
442 return true;
443 }
444 else {
445 return false;
446 }
447 }
448
remove(GFace * f)449 bool GModel::remove(GFace *f)
450 {
451 auto it = std::find(firstFace(), lastFace(), f);
452 if(it != faces.end()) {
453 faces.erase(it);
454 std::vector<GEdge *> const &e = f->edges();
455 for(auto it = e.begin(); it != e.end(); it++) (*it)->delFace(f);
456 return true;
457 }
458 else {
459 return false;
460 }
461 }
462
remove(GEdge * e)463 bool GModel::remove(GEdge *e)
464 {
465 auto it = std::find(firstEdge(), lastEdge(), e);
466 if(it != edges.end()) {
467 edges.erase(it);
468 if(e->getBeginVertex()) e->getBeginVertex()->delEdge(e);
469 if(e->getEndVertex()) e->getEndVertex()->delEdge(e);
470 return true;
471 }
472 else {
473 return false;
474 }
475 }
476
remove(GVertex * v)477 bool GModel::remove(GVertex *v)
478 {
479 auto it = std::find(firstVertex(), lastVertex(), v);
480 if(it != vertices.end()) {
481 vertices.erase(it);
482 return true;
483 }
484 else {
485 return false;
486 }
487 }
488
remove(int dim,int tag,std::vector<GEntity * > & removed,bool recursive)489 void GModel::remove(int dim, int tag, std::vector<GEntity*> &removed,
490 bool recursive)
491 {
492 if(dim == 3) {
493 GRegion *gr = getRegionByTag(tag);
494 if(gr) {
495 if(remove(gr)) removed.push_back(gr);
496 if(recursive) {
497 std::vector<GFace *> f = gr->faces();
498 for(auto it = f.begin(); it != f.end(); it++)
499 remove(2, (*it)->tag(), removed, recursive);
500 }
501 }
502 }
503 else if(dim == 2) {
504 GFace *gf = getFaceByTag(tag);
505 if(gf) {
506 bool skip = false;
507 for(auto it = firstRegion(); it != lastRegion(); it++) {
508 std::vector<GFace *> f = (*it)->faces();
509 if(std::find(f.begin(), f.end(), gf) != f.end()) {
510 skip = true;
511 break;
512 }
513 auto ef = (*it)->embeddedFaces();
514 if(std::find(ef.begin(), ef.end(), gf) != ef.end()) {
515 skip = true;
516 break;
517 }
518 }
519 if(!skip) {
520 if(remove(gf)) removed.push_back(gf);
521 if(recursive) {
522 std::vector<GEdge *> const &e = gf->edges();
523 for(auto it = e.begin(); it != e.end(); it++)
524 remove(1, (*it)->tag(), removed, recursive);
525 }
526 }
527 }
528 }
529 else if(dim == 1) {
530 GEdge *ge = getEdgeByTag(tag);
531 if(ge) {
532 bool skip = false;
533 for(auto it = firstRegion(); it != lastRegion(); it++) {
534 auto ee = (*it)->embeddedEdges();
535 if(std::find(ee.begin(), ee.end(), ge) != ee.end()) {
536 skip = true;
537 break;
538 }
539 }
540 if(!skip) {
541 for(auto it = firstFace(); it != lastFace(); it++) {
542 std::vector<GEdge *> const &e = (*it)->edges();
543 if(std::find(e.begin(), e.end(), ge) != e.end()) {
544 skip = true;
545 break;
546 }
547 auto ee = (*it)->embeddedEdges();
548 if(std::find(ee.begin(), ee.end(), ge) != ee.end()) {
549 skip = true;
550 break;
551 }
552 }
553 }
554 if(!skip) {
555 if(remove(ge)) removed.push_back(ge);
556 if(recursive) {
557 if(ge->getBeginVertex())
558 remove(0, ge->getBeginVertex()->tag(), removed);
559 if(ge->getEndVertex())
560 remove(0, ge->getEndVertex()->tag(), removed);
561 }
562 }
563 }
564 }
565 else if(dim == 0) {
566 GVertex *gv = getVertexByTag(tag);
567 if(gv) {
568 bool skip = false;
569 for(auto it = firstEdge(); it != lastEdge(); it++) {
570 GEdge *ge = *it;
571 if(ge->getBeginVertex() == gv || ge->getEndVertex() == gv) {
572 skip = true;
573 break;
574 }
575 }
576 if(!skip) {
577 for(auto it = firstFace(); it != lastFace(); it++) {
578 auto ev = (*it)->embeddedVertices();
579 if(std::find(ev.begin(), ev.end(), gv) != ev.end()) {
580 skip = true;
581 break;
582 }
583 }
584 }
585 if(!skip) {
586 for(auto it = firstRegion(); it != lastRegion(); it++) {
587 auto ev = (*it)->embeddedVertices();
588 if(std::find(ev.begin(), ev.end(), gv) != ev.end()) {
589 skip = true;
590 break;
591 }
592 }
593 }
594 if(!skip) {
595 if(remove(gv)) removed.push_back(gv);
596 }
597 }
598 }
599 }
600
remove(const std::vector<std::pair<int,int>> & dimTags,std::vector<GEntity * > & removed,bool recursive)601 void GModel::remove(const std::vector<std::pair<int, int> > &dimTags,
602 std::vector<GEntity*> &removed, bool recursive)
603 {
604 for(std::size_t i = 0; i < dimTags.size(); i++)
605 remove(dimTags[i].first, dimTags[i].second, removed, recursive);
606 }
607
remove()608 void GModel::remove()
609 {
610 regions.clear();
611 faces.clear();
612 edges.clear();
613 vertices.clear();
614 }
615
snapVertices()616 void GModel::snapVertices()
617 {
618 if(!CTX::instance()->geom.snapPoints) return;
619
620 auto vit = firstVertex();
621
622 double tol = CTX::instance()->geom.tolerance;
623
624 while(vit != lastVertex()) {
625 std::vector<GEdge *> const &edges = (*vit)->edges();
626 for(auto it = edges.begin(); it != edges.end(); ++it) {
627 Range<double> parb = (*it)->parBounds(0);
628 double t;
629 if((*it)->getBeginVertex() == *vit) { t = parb.low(); }
630 else if((*it)->getEndVertex() == *vit) {
631 t = parb.high();
632 }
633 else {
634 Msg::Error("Weird point: impossible to snap");
635 break;
636 }
637 GPoint gp = (*it)->point(t);
638 double d = sqrt((gp.x() - (*vit)->x()) * (gp.x() - (*vit)->x()) +
639 (gp.y() - (*vit)->y()) * (gp.y() - (*vit)->y()) +
640 (gp.z() - (*vit)->z()) * (gp.z() - (*vit)->z()));
641 if(d > tol) {
642 (*vit)->setPosition(gp);
643 Msg::Info("Snapping geometry point %d to curve (distance = %g)",
644 (*vit)->tag(), d);
645 }
646 }
647 vit++;
648 }
649 }
650
getEntities(std::vector<GEntity * > & entities,int dim) const651 void GModel::getEntities(std::vector<GEntity *> &entities, int dim) const
652 {
653 entities.clear();
654 switch(dim) {
655 case 0:
656 entities.insert(entities.end(), vertices.begin(), vertices.end());
657 break;
658 case 1: entities.insert(entities.end(), edges.begin(), edges.end()); break;
659 case 2: entities.insert(entities.end(), faces.begin(), faces.end()); break;
660 case 3:
661 entities.insert(entities.end(), regions.begin(), regions.end());
662 break;
663 default:
664 entities.insert(entities.end(), vertices.begin(), vertices.end());
665 entities.insert(entities.end(), edges.begin(), edges.end());
666 entities.insert(entities.end(), faces.begin(), faces.end());
667 entities.insert(entities.end(), regions.begin(), regions.end());
668 break;
669 }
670 }
671
getEntitiesInBox(std::vector<GEntity * > & entities,const SBoundingBox3d & box,int dim) const672 void GModel::getEntitiesInBox(std::vector<GEntity *> &entities,
673 const SBoundingBox3d &box, int dim) const
674 {
675 entities.clear();
676 std::vector<GEntity *> all;
677 getEntities(all, dim);
678 // if we use this often, create an rtree to avoid the linear search
679 for(std::size_t i = 0; i < all.size(); i++) {
680 SBoundingBox3d bbox = all[i]->bounds();
681 if(bbox.min().x() >= box.min().x() && bbox.max().x() <= box.max().x() &&
682 bbox.min().y() >= box.min().y() && bbox.max().y() <= box.max().y() &&
683 bbox.min().z() >= box.min().z() && bbox.max().z() <= box.max().z())
684 entities.push_back(all[i]);
685 }
686 }
687
688 class AbsIntLessThan {
689 public:
operator ()(const int & i1,const int & i2) const690 bool operator()(const int &i1, const int &i2) const
691 {
692 if(std::abs(i1) < std::abs(i2)) return true;
693 return false;
694 }
695 };
696
getBoundaryTags(const std::vector<std::pair<int,int>> & inDimTags,std::vector<std::pair<int,int>> & outDimTags,bool combined,bool oriented,bool recursive)697 bool GModel::getBoundaryTags(const std::vector<std::pair<int, int> > &inDimTags,
698 std::vector<std::pair<int, int> > &outDimTags,
699 bool combined, bool oriented, bool recursive)
700 {
701 bool ret = true;
702 for(std::size_t i = 0; i < inDimTags.size(); i++) {
703 int dim = inDimTags[i].first;
704 int tag = std::abs(inDimTags[i].second); // abs for backward compatibility
705 bool reverse = (inDimTags[i].second < 0);
706 if(dim == 3) {
707 GRegion *gr = getRegionByTag(tag);
708 if(gr) {
709 if(recursive) {
710 std::vector<GVertex *> const &vert = gr->vertices();
711 for(auto it = vert.begin(); it != vert.end(); it++)
712 outDimTags.push_back(std::make_pair(0, (*it)->tag()));
713 }
714 else {
715 std::vector<GFace *> faces(gr->faces());
716 std::vector<int> const &orientations = gr->faceOrientations();
717 auto ito = orientations.begin();
718 for(auto it = faces.begin(); it != faces.end(); it++) {
719 int t = (*it)->tag();
720 if(oriented && ito != orientations.end()) {
721 t *= *ito;
722 ito++;
723 }
724 outDimTags.push_back(std::make_pair(2, t));
725 }
726 }
727 }
728 else {
729 Msg::Error("Unknown model region with tag %d", tag);
730 ret = false;
731 }
732 }
733 else if(dim == 2) {
734 GFace *gf = getFaceByTag(tag);
735 if(gf) {
736 if(recursive) {
737 std::vector<GVertex *> const &vert = gf->vertices();
738 for(auto it = vert.begin(); it != vert.end(); it++)
739 outDimTags.push_back(std::make_pair(0, (*it)->tag()));
740 }
741 else {
742 std::vector<GEdge *> const &edges = gf->edges();
743 std::vector<int> orientations(gf->edgeOrientations());
744 auto ito = orientations.begin();
745 for(auto it = edges.begin(); it != edges.end(); it++) {
746 int t = (*it)->tag();
747 if(oriented && ito != orientations.end()) {
748 t *= *ito;
749 ito++;
750 }
751 outDimTags.push_back(std::make_pair(1, t));
752 }
753 }
754 }
755 else {
756 Msg::Error("Unknown model face with tag %d", tag);
757 ret = false;
758 }
759 }
760 else if(dim == 1) {
761 GEdge *ge = getEdgeByTag(tag);
762 if(ge) {
763 if(reverse) { // for backward compatibility
764 if(ge->getEndVertex())
765 outDimTags.push_back(
766 std::make_pair(0, ge->getEndVertex()->tag()));
767 if(ge->getBeginVertex())
768 outDimTags.push_back(
769 std::make_pair(0, ge->getBeginVertex()->tag()));
770 }
771 else {
772 if(ge->getBeginVertex())
773 outDimTags.push_back(
774 std::make_pair(0, ge->getBeginVertex()->tag()));
775 if(ge->getEndVertex())
776 outDimTags.push_back(
777 std::make_pair(0, ge->getEndVertex()->tag()));
778 }
779 }
780 else {
781 Msg::Error("Unknown model curve with tag %d", tag);
782 ret = false;
783 }
784 }
785 else if(dim == 0) {
786 GVertex *gv = getVertexByTag(tag);
787 if(gv && recursive) {
788 outDimTags.push_back(std::make_pair(0, gv->tag()));
789 }
790 }
791 }
792
793 if(combined) {
794 // compute boundary of the combined shapes
795 std::set<int, AbsIntLessThan> c[3];
796 for(std::size_t i = 0; i < outDimTags.size(); i++) {
797 int dim = outDimTags[i].first;
798 int tag = outDimTags[i].second;
799 if(dim >= 0 && dim < 3) {
800 auto it = c[dim].find(tag);
801 if(it == c[dim].end())
802 c[dim].insert(tag);
803 else {
804 c[dim].erase(it);
805 }
806 }
807 }
808 outDimTags.clear();
809 for(int dim = 0; dim < 3; dim++) {
810 for(auto it = c[dim].begin(); it != c[dim].end(); it++)
811 outDimTags.push_back(std::make_pair(dim, *it));
812 }
813 }
814 return ret;
815 }
816
getMaxElementaryNumber(int dim)817 int GModel::getMaxElementaryNumber(int dim)
818 {
819 std::vector<GEntity *> entities;
820 getEntities(entities);
821 int num = 0;
822 for(std::size_t i = 0; i < entities.size(); i++)
823 if(dim < 0 || entities[i]->dim() == dim)
824 num = std::max(num, std::abs(entities[i]->tag()));
825 return num;
826 }
827
noPhysicalGroups()828 bool GModel::noPhysicalGroups()
829 {
830 std::vector<GEntity *> entities;
831 getEntities(entities);
832 for(std::size_t i = 0; i < entities.size(); i++)
833 if(entities[i]->physicals.size()) return false;
834 return true;
835 }
836
getPhysicalGroups(std::map<int,std::vector<GEntity * >> groups[4]) const837 void GModel::getPhysicalGroups(
838 std::map<int, std::vector<GEntity *> > groups[4]) const
839 {
840 std::vector<GEntity *> entities;
841 getEntities(entities);
842 for(std::size_t i = 0; i < entities.size(); i++) {
843 std::map<int, std::vector<GEntity *> > &group(groups[entities[i]->dim()]);
844 for(std::size_t j = 0; j < entities[i]->physicals.size(); j++) {
845 // physicals can be stored with negative signs when the entity should be
846 // "reversed"
847 int p = std::abs(entities[i]->physicals[j]);
848 group[p].push_back(entities[i]);
849 }
850 }
851 for(int dim = 0; dim < 4; ++dim) {
852 std::map<int, std::vector<GEntity *> > &group(groups[dim]);
853 for(auto it = group.begin(); it != group.end(); ++it) {
854 std::vector<GEntity *> &v = it->second;
855 std::sort(v.begin(), v.end(), GEntityPtrLessThan());
856 std::unique(v.begin(), v.end(), GEntityPtrLessThan());
857 }
858 }
859 }
860
getPhysicalGroups(int dim,std::map<int,std::vector<GEntity * >> & groups) const861 void GModel::getPhysicalGroups(
862 int dim, std::map<int, std::vector<GEntity *> > &groups) const
863 {
864 std::vector<GEntity *> entities;
865 getEntities(entities, dim);
866 for(std::size_t i = 0; i < entities.size(); i++) {
867 for(std::size_t j = 0; j < entities[i]->physicals.size(); j++) {
868 // physicals can be stored with negative signs when the entity should be
869 // "reversed"
870 int p = std::abs(entities[i]->physicals[j]);
871 groups[p].push_back(entities[i]);
872 }
873 }
874 for(auto it = groups.begin(); it != groups.end(); ++it) {
875 std::vector<GEntity *> &v = it->second;
876 std::sort(v.begin(), v.end(), GEntityPtrLessThan());
877 std::unique(v.begin(), v.end(), GEntityPtrLessThan());
878 }
879 }
880
addPhysicalGroup(int dim,int tag,const std::vector<int> & tags)881 void GModel::addPhysicalGroup(int dim, int tag, const std::vector<int> &tags)
882 {
883 for(auto t : tags) {
884 GEntity *ge = getEntityByTag(dim, std::abs(t));
885 if(ge)
886 ge->physicals.push_back((t > 0) ? tag : -tag);
887 else
888 Msg::Warning("Unknown entity of dimension %d and tag %d in physical "
889 "group %d", dim, t, tag);
890 }
891 }
892
removePhysicalGroups()893 void GModel::removePhysicalGroups()
894 {
895 std::vector<GEntity *> entities;
896 getEntities(entities);
897 for(std::size_t i = 0; i < entities.size(); i++)
898 entities[i]->physicals.clear();
899 }
900
removePhysicalGroup(int dim,int tag)901 void GModel::removePhysicalGroup(int dim, int tag)
902 {
903 // FIXME: this is very inefficient - needs to be rewriten (and we should
904 // generalize the function by taking a list of dim, tag pairs)
905 std::vector<GEntity *> entities;
906 getEntities(entities, dim);
907 for(std::size_t i = 0; i < entities.size(); i++) {
908 std::vector<int> p;
909 for(std::size_t j = 0; j < entities[i]->physicals.size(); j++)
910 if(std::abs(entities[i]->physicals[j]) != tag)
911 p.push_back(entities[i]->physicals[j]);
912 entities[i]->physicals = p;
913 }
914 _physicalNames.erase(std::make_pair(dim, tag));
915 }
916
getMaxPhysicalNumber(int dim)917 int GModel::getMaxPhysicalNumber(int dim)
918 {
919 std::vector<GEntity *> entities;
920 getEntities(entities);
921 int num = 0;
922 for(std::size_t i = 0; i < entities.size(); i++)
923 if(dim < 0 || entities[i]->dim() == dim)
924 for(std::size_t j = 0; j < entities[i]->physicals.size(); j++)
925 num = std::max(num, std::abs(entities[i]->physicals[j]));
926 return num;
927 }
928
getInnerPhysicalNamesIterators(std::vector<piter> & iterators)929 void GModel::getInnerPhysicalNamesIterators(std::vector<piter> &iterators)
930 {
931 iterators.resize(4, firstPhysicalName());
932
933 for(auto physIt = firstPhysicalName(); physIt != lastPhysicalName(); ++physIt)
934 iterators[physIt->first.first] = physIt;
935 }
936
setPhysicalName(const std::string & name,int dim,int number)937 int GModel::setPhysicalName(const std::string &name, int dim, int number)
938 {
939 // check if the name is already used
940 int findPhy = getPhysicalNumber(dim, name);
941 if(findPhy != -1) return findPhy;
942
943 // if no number is given, find the next available one
944 if(!number) number = getMaxPhysicalNumber(dim) + 1;
945 _physicalNames.insert(std::make_pair(std::make_pair(dim, number), name));
946 return number;
947 }
948
setPhysicalName(piter pos,const std::string & name,int dim,int number)949 GModel::piter GModel::setPhysicalName(piter pos, const std::string &name,
950 int dim, int number)
951 {
952 // if no number is given, find the next available one
953 if(!number) number = getMaxPhysicalNumber(dim) + 1;
954 // Insertion complexity in O(1) if position points to the element that will
955 // FOLLOW the inserted element.
956 if(pos != lastPhysicalName()) ++pos;
957 return _physicalNames.insert(pos, std::make_pair(std::make_pair(dim, number),
958 name));
959 }
960
getPhysicalName(int dim,int number) const961 std::string GModel::getPhysicalName(int dim, int number) const
962 {
963 auto it = _physicalNames.find(std::make_pair(dim, number));
964 if(it != _physicalNames.end()) return it->second;
965 return "";
966 }
967
removePhysicalName(const std::string & name)968 void GModel::removePhysicalName(const std::string &name)
969 {
970 auto it = _physicalNames.begin();
971 while(it != _physicalNames.end()) {
972 if(it->second == name)
973 // it = _physicalNames.erase(it); // C++11 only
974 _physicalNames.erase(it++);
975 else
976 ++it;
977 }
978 }
979
getPhysicalNumber(const int & dim,const std::string & name)980 int GModel::getPhysicalNumber(const int &dim, const std::string &name)
981 {
982 for(auto physIt = firstPhysicalName(); physIt != lastPhysicalName(); ++physIt)
983 if(dim == physIt->first.first && name == physIt->second)
984 return physIt->first.second;
985
986 return -1;
987 }
988
getDim() const989 int GModel::getDim() const
990 {
991 if(getNumRegions() > 0) return 3;
992 if(getNumFaces() > 0) return 2;
993 if(getNumEdges() > 0) return 1;
994 if(getNumVertices() > 0) return 0;
995 return -1;
996 }
997
getMeshDim() const998 int GModel::getMeshDim() const
999 {
1000 if(getNumMeshElements(3)) return 3;
1001 if(getNumMeshElements(2)) return 2;
1002 if(getNumMeshElements(1)) return 1;
1003 if(getNumMeshElements(0)) return 0;
1004 return -1;
1005 }
1006
getElementaryName(int dim,int tag)1007 std::string GModel::getElementaryName(int dim, int tag)
1008 {
1009 auto it = _elementaryNames.find(std::make_pair(dim, tag));
1010 if(it != _elementaryNames.end()) return it->second;
1011 return "";
1012 }
1013
removeElementaryName(const std::string & name)1014 void GModel::removeElementaryName(const std::string &name)
1015 {
1016 auto it = _elementaryNames.begin();
1017 while(it != _elementaryNames.end()) {
1018 if(it->second == name)
1019 // it = _elementaryNames.erase(it); // C++11 only
1020 _elementaryNames.erase(it++);
1021 else
1022 ++it;
1023 }
1024 }
1025
setSelection(int val)1026 void GModel::setSelection(int val)
1027 {
1028 std::vector<GEntity *> entities;
1029 getEntities(entities);
1030
1031 for(std::size_t i = 0; i < entities.size(); i++) {
1032 entities[i]->setSelection(val);
1033 // reset selection in elements (stored in the visibility flag to
1034 // save space)
1035 if(val == 0) {
1036 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++)
1037 if(entities[i]->getMeshElement(j)->getVisibility() == 2)
1038 entities[i]->getMeshElement(j)->setVisibility(1);
1039 }
1040 }
1041 }
1042
bounds(bool aroundVisible)1043 SBoundingBox3d GModel::bounds(bool aroundVisible)
1044 {
1045 std::vector<GEntity *> entities;
1046 getEntities(entities);
1047 SBoundingBox3d bb;
1048 for(std::size_t i = 0; i < entities.size(); i++) {
1049 if(!aroundVisible || entities[i]->getVisibility()) {
1050 if(entities[i]->getNativeType() == GEntity::OpenCascadeModel) {
1051 bb += entities[i]->bounds();
1052 }
1053 else {
1054 // using the mesh vertices for now
1055 if(entities[i]->dim() == 0)
1056 bb += static_cast<GVertex *>(entities[i])->xyz();
1057 else
1058 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
1059 bb += entities[i]->mesh_vertices[j]->point();
1060 }
1061 }
1062 }
1063 return bb;
1064 }
1065
mesh(int dimension)1066 int GModel::mesh(int dimension)
1067 {
1068 #if defined(HAVE_MESH)
1069 GenerateMesh(this, dimension);
1070 if(CTX::instance()->mesh.renumber) {
1071 renumberMeshVertices();
1072 renumberMeshElements();
1073 }
1074 computeHomology(); // must be done after renumbering
1075 CTX::instance()->mesh.changed = ENT_ALL;
1076 return true;
1077 #else
1078 Msg::Error("Mesh module not compiled");
1079 return false;
1080 #endif
1081 }
1082
setAllVolumesPositive()1083 bool GModel::setAllVolumesPositive()
1084 {
1085 bool ok = true;
1086 for(auto it = regions.begin(); it != regions.end(); ++it)
1087 for(std::size_t i = 0; i < (*it)->getNumMeshElements(); ++i)
1088 if(!(*it)->getMeshElement(i)->setVolumePositive()) ok = false;
1089 return ok;
1090 }
1091
1092 static void
addToMap(std::multimap<MFace,MElement *,MFaceLessThan> & faceToElement,std::map<MElement *,std::vector<std::pair<MElement *,bool>>> & elToNeighbors,const MFace & face,MElement * el)1093 addToMap(std::multimap<MFace, MElement *, MFaceLessThan> &faceToElement,
1094 std::map<MElement *, std::vector<std::pair<MElement *, bool> > >
1095 &elToNeighbors,
1096 const MFace &face, MElement *el)
1097 {
1098 auto fit = faceToElement.find(face);
1099 if(fit == faceToElement.end()) {
1100 faceToElement.insert(std::make_pair(face, el));
1101 }
1102 else { // We found the neighbor face outFace
1103 faceToElement.insert(std::make_pair(face, el));
1104 if(faceToElement.count(face) > 2) {
1105 Msg::Error(
1106 "Topological fault: Face sharing two other faces. Element %i. "
1107 "Number of nodes %i. Count of faces: %i Three first nodes %i %i %i",
1108 el->getNum(), face.getNumVertices(), faceToElement.count(face),
1109 face.getVertex(0)->getNum(), face.getVertex(1)->getNum(),
1110 face.getVertex(2)->getNum());
1111 return;
1112 }
1113 MFace outFace = fit->first;
1114 std::vector<std::pair<MElement *, bool> > &neigh = elToNeighbors[el];
1115 for(size_t iN = 0; iN < neigh.size(); iN++)
1116 if(neigh[iN].first == fit->second) // Neigbor is already in the map
1117 return;
1118 int i0 = -1;
1119 while(face.getVertex(0) != outFace.getVertex(++i0))
1120 ;
1121 bool sameOrientation =
1122 face.getVertex(1) == outFace.getVertex((i0 + 1) % face.getNumVertices());
1123 neigh.push_back(std::make_pair(fit->second, !sameOrientation));
1124 elToNeighbors[fit->second].push_back(std::make_pair(el, !sameOrientation));
1125 }
1126 }
1127
1128 static void
checkConformity(std::multimap<MFace,MElement *,MFaceLessThan> & faceToElement,std::map<MElement *,std::vector<std::pair<MElement *,bool>>> & elToNeighbors,const MFace & face,MElement * el)1129 checkConformity(std::multimap<MFace, MElement *, MFaceLessThan> &faceToElement,
1130 std::map<MElement *, std::vector<std::pair<MElement *, bool> > >
1131 &elToNeighbors,
1132 const MFace &face, MElement *el)
1133 {
1134 int connectivity = faceToElement.count(face);
1135 if(ElementType::getParentType(el->getType()) == TYPE_TRIH) {
1136 // Each face of a trihedron should exist twice (no face on the boundary)
1137 if(connectivity != 2)
1138 Msg::Error("Non conforming trihedron %i (nb connections for a face %i)",
1139 el->getNum(), faceToElement.count(face));
1140 }
1141 else {
1142 // A face can exist twice (inside) or once (boundary)
1143 if(connectivity != 2) {
1144 for(std::size_t iV = 0; iV < face.getNumVertices(); iV++)
1145 if(face.getVertex(iV)->onWhat()->dim() == 3 || connectivity != 1) {
1146 for(std::size_t jV = 0; jV < face.getNumVertices(); jV++)
1147 Msg::Info("Node %i on entity of dim %i", face.getVertex(jV)->getNum(),
1148 face.getVertex(iV)->onWhat()->dim());
1149 Msg::Error("Non conforming element %i (%i connection(s) for a face "
1150 "located on dim %i (vertex %i))",
1151 el->getNum(), connectivity,
1152 face.getVertex(iV)->onWhat()->dim(),
1153 face.getVertex(iV)->getNum());
1154 }
1155 }
1156 }
1157 }
1158
setAllVolumesPositiveTopology()1159 void GModel::setAllVolumesPositiveTopology()
1160 {
1161 Msg::Info("Orienting volumes according to topology");
1162 std::map<MElement *, std::vector<std::pair<MElement *, bool> > >
1163 elToNeighbors;
1164 std::multimap<MFace, MElement *, MFaceLessThan> faceToElement;
1165
1166 MElement *el;
1167 for(auto it = regions.begin(); it != regions.end(); ++it) {
1168 for(std::size_t iEl = 0; iEl < (*it)->getNumMeshElements(); ++iEl) {
1169 el = (*it)->getMeshElement(iEl);
1170 for(int iFace = 0; iFace < el->getNumFaces(); iFace++) {
1171 addToMap(faceToElement, elToNeighbors, el->getFace(iFace), el);
1172 }
1173 }
1174 }
1175 for(auto it = regions.begin(); it != regions.end(); ++it) {
1176 for(std::size_t iEl = 0; iEl < (*it)->getNumMeshElements(); ++iEl) {
1177 el = (*it)->getMeshElement(iEl);
1178 for(int iFace = 0; iFace < el->getNumFaces(); iFace++) {
1179 checkConformity(faceToElement, elToNeighbors, el->getFace(iFace), el);
1180 }
1181 }
1182 }
1183 std::vector<std::pair<MElement *, bool> > queue;
1184 std::set<MElement *> queued;
1185 if((*regions.begin())->tetrahedra.size() == 0) {
1186 Msg::Error(
1187 "setAllVolumePositiveTopology needs at least one tetrahedron to start");
1188 return;
1189 }
1190 el = (*regions.begin())->tetrahedra[0];
1191 queue.push_back(std::make_pair(el, true));
1192 for(size_t i = 0; i < queue.size(); i++) {
1193 el = queue[i].first;
1194 if(!queue[i].second) {
1195 el->reverse();
1196 // Msg::Info("Reverted element %i of type %i", el->getNum(),
1197 // el->getType());
1198 }
1199 const std::vector<std::pair<MElement *, bool> > &neigh = elToNeighbors[el];
1200 for(size_t iN = 0; iN < neigh.size(); iN++)
1201 if(queued.count(neigh[iN].first) == 0) {
1202 queue.push_back(
1203 std::make_pair(neigh[iN].first, neigh[iN].second == queue[i].second));
1204 // if(!(neigh[iN].second == queue[i].second))
1205 // Msg::Info("Queuing element %i (%i) from el %i (%i)",
1206 // neigh[iN].first->getNum(), neigh[iN].second,
1207 // el->getNum(), queue[i].second);
1208 queued.insert(neigh[iN].first);
1209 }
1210 }
1211 }
1212
adaptMesh()1213 int GModel::adaptMesh()
1214 {
1215 #if defined(HAVE_MESH)
1216 AdaptMesh(this);
1217 if(CTX::instance()->mesh.renumber) {
1218 renumberMeshVertices();
1219 renumberMeshElements();
1220 }
1221 CTX::instance()->mesh.changed = ENT_ALL;
1222 return 1;
1223 #else
1224 Msg::Error("Mesh module not compiled");
1225 return 0;
1226 #endif
1227 }
1228
adaptMesh(std::vector<int> technique,std::vector<simpleFunction<double> * > f,std::vector<std::vector<double>> parameters,int niter,bool meshAll)1229 int GModel::adaptMesh(std::vector<int> technique,
1230 std::vector<simpleFunction<double> *> f,
1231 std::vector<std::vector<double> > parameters, int niter,
1232 bool meshAll)
1233 {
1234 // For all algorithms:
1235 //
1236 // parameters[1] = lcmin (default : in global gmsh options
1237 // CTX::instance()->mesh.lcMin)
1238 // parameters[2] = lcmax (default : in global gmsh options
1239 // CTX::instance()->mesh.lcMax) niter is the maximum number of iterations
1240
1241 // Available algorithms:
1242 //
1243 // 1) Assume that the function is a levelset -> adapt using Coupez
1244 // technique (technique = 1)
1245 // parameters[0] = thickness of the interface (mandatory)
1246 // 2) Assume that the function is a physical quantity -> adapt using the
1247 // Hessian (technique = 2)
1248 // parameters[0] = N, the final number of elements
1249 // 3) A variant of 1) by P. Frey (= Coupez + takes curvature function into
1250 // account)
1251 // parameters[0] = thickness of the interface (mandatory)
1252 // parameters[3] = the required minimum number of elements to
1253 // represent a circle - used for curvature-based metric (default:
1254 // = 15)
1255 // 4) A variant (3), direct implementation in the metric eigendirections,
1256 // assuming a level set (ls):
1257 // - hmin is imposed in the ls gradient,
1258 // - hmax is imposed in the two eigendirections of the ls hessian that
1259 // are
1260 // (almost ?) tangent to the iso-zero plane
1261 // + the latter eigenvalues (1/hmax^2) are modified if necessary to
1262 // capture the iso-zero curvature
1263 // parameters[0] = thickness of the interface in the positive ls
1264 // direction (mandatory) parameters[4] = thickness of the interface in
1265 // the negative ls direction
1266 // (=parameters[0] if not specified)
1267 // parameters[3] = the required minimum number of elements to represent a
1268 // circle
1269 // - used for curvature-based metric (default: = 15)
1270 // 5) Same as 4, except that the transition in band E uses linear
1271 // interpolation
1272 // of h, instead of linear interpolation of metric
1273
1274 #if defined(HAVE_MESH)
1275 // copy context (in order to allow multiple calls)
1276 CTX _backup = *(CTX::instance());
1277
1278 if(getNumMeshElements() == 0) mesh(getDim());
1279 int nbElemsOld = getNumMeshElements();
1280 int nbElems;
1281
1282 FieldManager *fields = getFields();
1283 fields->reset();
1284
1285 int ITER = 0;
1286 if(meshAll) {
1287 while(1) {
1288 Msg::Info(" - Adapt mesh (all dimensions) iter. = %d", ITER);
1289 fields->reset();
1290 meshMetric *metric = new meshMetric(this);
1291 for(std::size_t imetric = 0; imetric < technique.size(); imetric++) {
1292 metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
1293 }
1294 fields->setBackgroundField(metric);
1295
1296 opt_mesh_lc_integration_precision(0, GMSH_SET, 1.e-4);
1297 opt_mesh_algo2d(0, GMSH_SET, 7.0); // bamg
1298 opt_mesh_algo3d(0, GMSH_SET, 7.0); // mmg3d
1299 opt_mesh_lc_from_points(0, GMSH_SET, 0.0); // do not mesh lines with lc
1300
1301 std::for_each(firstRegion(), lastRegion(), deMeshGRegion());
1302 std::for_each(firstFace(), lastFace(), deMeshGFace());
1303 std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
1304
1305 GenerateMesh(this, getDim());
1306 nbElems = getNumMeshElements();
1307
1308 char name[256];
1309 sprintf(name, "meshAdapt-%d.msh", ITER);
1310 writeMSH(name);
1311 // metric->exportInfo(name);
1312
1313 if(ITER++ >= niter) break;
1314 if(ITER > 3 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld)
1315 break;
1316
1317 nbElemsOld = nbElems;
1318 }
1319 }
1320 else { // adapt only upper most dimension
1321 while(1) {
1322 Msg::Info(" - Adapt mesh iter. = %d ", ITER);
1323 std::vector<MElement *> elements;
1324
1325 if(getDim() == 2) {
1326 for(auto fit = firstFace(); fit != lastFace(); ++fit) {
1327 if((*fit)->quadrangles.size()) return -1;
1328 for(unsigned i = 0; i < (*fit)->triangles.size(); i++) {
1329 elements.push_back((*fit)->triangles[i]);
1330 }
1331 }
1332 }
1333 else if(getDim() == 3) {
1334 for(auto rit = firstRegion(); rit != lastRegion(); ++rit) {
1335 if((*rit)->hexahedra.size()) return -1;
1336 for(unsigned i = 0; i < (*rit)->tetrahedra.size(); i++) {
1337 elements.push_back((*rit)->tetrahedra[i]);
1338 }
1339 }
1340 }
1341
1342 if(elements.size() == 0) return -1;
1343
1344 fields->reset();
1345 meshMetric *metric = new meshMetric(this);
1346 for(std::size_t imetric = 0; imetric < technique.size(); imetric++) {
1347 metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
1348 }
1349 fields->setBackgroundField(metric);
1350
1351 if(getDim() == 2) {
1352 for(auto fit = firstFace(); fit != lastFace(); ++fit) {
1353 if((*fit)->geomType() != GEntity::DiscreteSurface) {
1354 meshGFaceBamg(*fit);
1355 laplaceSmoothing(*fit, CTX::instance()->mesh.nbSmoothing);
1356 }
1357 if(_elementOctree) delete _elementOctree;
1358 _elementOctree = nullptr;
1359 }
1360 }
1361 else if(getDim() == 3) {
1362 for(auto rit = firstRegion(); rit != lastRegion(); ++rit) {
1363 refineMeshMMG(*rit);
1364 if(_elementOctree) delete _elementOctree;
1365 _elementOctree = nullptr;
1366 }
1367 }
1368
1369 char name[256];
1370 sprintf(name, "meshAdapt-%d.msh", ITER);
1371 writeMSH(name);
1372
1373 nbElems = getNumMeshElements();
1374 if(++ITER >= niter) break;
1375 if(ITER > 3 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld)
1376 break;
1377
1378 nbElemsOld = nbElems;
1379 }
1380 }
1381
1382 fields->reset();
1383 // copy context (in order to allow multiple calls)
1384 *(CTX::instance()) = _backup;
1385
1386 if(CTX::instance()->mesh.renumber) {
1387 renumberMeshVertices();
1388 renumberMeshElements();
1389 }
1390 CTX::instance()->mesh.changed = ENT_ALL;
1391
1392 return 0;
1393 #else
1394 Msg::Error("Mesh module not compiled");
1395 return -1;
1396 #endif
1397 }
1398
refineMesh(int linear,bool splitIntoQuads,bool splitIntoHexas,bool barycentric)1399 int GModel::refineMesh(int linear, bool splitIntoQuads, bool splitIntoHexas,
1400 bool barycentric)
1401 {
1402 #if defined(HAVE_MESH)
1403 if(!barycentric) { RefineMesh(this, linear, splitIntoQuads, splitIntoHexas); }
1404 else {
1405 BarycentricRefineMesh(this);
1406 }
1407 if(CTX::instance()->mesh.renumber) {
1408 renumberMeshVertices();
1409 renumberMeshElements();
1410 }
1411 CTX::instance()->mesh.changed = ENT_ALL;
1412 return 1;
1413 #else
1414 Msg::Error("Mesh module not compiled");
1415 return 0;
1416 #endif
1417 }
1418
recombineMesh()1419 int GModel::recombineMesh()
1420 {
1421 #if defined(HAVE_MESH)
1422 RecombineMesh(this);
1423 if(CTX::instance()->mesh.renumber) {
1424 renumberMeshVertices();
1425 renumberMeshElements();
1426 }
1427 CTX::instance()->mesh.changed = ENT_ALL;
1428 return 1;
1429 #else
1430 Msg::Error("Mesh module not compiled");
1431 return 0;
1432 #endif
1433 }
1434
optimizeMesh(const std::string & how,const bool force,int niter)1435 int GModel::optimizeMesh(const std::string &how, const bool force, int niter)
1436 {
1437 #if defined(HAVE_MESH)
1438 OptimizeMesh(this, how, force, niter);
1439 FixPeriodicMesh(this);
1440 if(CTX::instance()->mesh.renumber) {
1441 renumberMeshVertices();
1442 renumberMeshElements();
1443 }
1444 CTX::instance()->mesh.changed = ENT_ALL;
1445 return 1;
1446 #else
1447 Msg::Error("Mesh module not compiled");
1448 return 0;
1449 #endif
1450 }
1451
setOrderN(int order,int linear,int incomplete)1452 int GModel::setOrderN(int order, int linear, int incomplete)
1453 {
1454 #if defined(HAVE_MESH)
1455 if(order > 1)
1456 SetOrderN(this, order, linear, incomplete);
1457 else
1458 SetOrder1(this);
1459 FixPeriodicMesh(this);
1460 if(CTX::instance()->mesh.renumber) {
1461 renumberMeshVertices();
1462 renumberMeshElements();
1463 }
1464 CTX::instance()->mesh.changed = ENT_ALL;
1465 return true;
1466 #else
1467 Msg::Error("Mesh module not compiled");
1468 return false;
1469 #endif
1470 }
1471
getMeshStatus(bool countDiscrete)1472 int GModel::getMeshStatus(bool countDiscrete)
1473 {
1474 std::size_t numEle3D = 0;
1475 bool toMesh3D = false;
1476 bool onlyVisible = CTX::instance()->mesh.meshOnlyVisible;
1477
1478 for(auto it = firstRegion(); it != lastRegion(); ++it) {
1479 GRegion *gr = *it;
1480 if(countDiscrete || gr->geomType() != GEntity::DiscreteVolume)
1481 numEle3D += gr->getNumMeshElements();
1482 if(countDiscrete && numEle3D) return 3;
1483 if(gr->geomType() != GEntity::DiscreteVolume &&
1484 gr->meshAttributes.method != MESH_NONE)
1485 toMesh3D = true;
1486 }
1487 if(numEle3D && toMesh3D) return 3;
1488
1489 std::size_t numEle2D = 0;
1490 bool toMesh2D = false, meshDone2D = true;
1491 for(auto it = firstFace(); it != lastFace(); ++it) {
1492 GFace *gf = *it;
1493 if(countDiscrete || gf->geomType() != GEntity::DiscreteSurface)
1494 numEle2D += gf->getNumMeshElements();
1495 if(countDiscrete && numEle2D) return 2;
1496 if(gf->geomType() != GEntity::DiscreteSurface &&
1497 gf->meshAttributes.method != MESH_NONE)
1498 toMesh2D = true;
1499 if(gf->meshStatistics.status != GEntity::DONE &&
1500 (!onlyVisible || (onlyVisible && gf->getVisibility())))
1501 meshDone2D = false;
1502 }
1503 if(numEle2D && toMesh2D && meshDone2D) return 2;
1504
1505 std::size_t numEle1D = 0;
1506 bool toMesh1D = false, meshDone1D = true;
1507 for(auto it = firstEdge(); it != lastEdge(); ++it) {
1508 GEdge *ge = *it;
1509 if(countDiscrete || ge->geomType() != GEntity::DiscreteCurve)
1510 numEle1D += ge->getNumMeshElements();
1511 if(countDiscrete && numEle1D) return 1;
1512 if(ge->geomType() != GEntity::DiscreteCurve &&
1513 ge->meshAttributes.method != MESH_NONE)
1514 toMesh1D = true;
1515 if(ge->meshStatistics.status != GEntity::DONE &&
1516 (!onlyVisible || (onlyVisible && ge->getVisibility())))
1517 meshDone1D = false;
1518 }
1519 if(numEle1D && toMesh1D && meshDone1D) return 1;
1520
1521 for(auto it = firstVertex(); it != lastVertex(); ++it)
1522 if((*it)->mesh_vertices.size()) return 0;
1523
1524 return -1;
1525 }
1526
getNumMeshVertices(int dim) const1527 std::size_t GModel::getNumMeshVertices(int dim) const
1528 {
1529 std::vector<GEntity *> entities;
1530 getEntities(entities);
1531 std::size_t n = 0;
1532 for(std::size_t i = 0; i < entities.size(); i++)
1533 if(entities[i]->dim() == dim || dim < 0)
1534 n += entities[i]->getNumMeshVertices();
1535 return n;
1536 }
1537
getNumMeshElements(int dim) const1538 std::size_t GModel::getNumMeshElements(int dim) const
1539 {
1540 std::vector<GEntity *> entities;
1541 getEntities(entities);
1542 std::size_t n = 0;
1543 for(std::size_t i = 0; i < entities.size(); i++)
1544 if(entities[i]->dim() == dim || dim < 0)
1545 n += entities[i]->getNumMeshElements();
1546 return n;
1547 }
1548
getNumMeshParentElements() const1549 std::size_t GModel::getNumMeshParentElements() const
1550 {
1551 std::vector<GEntity *> entities;
1552 getEntities(entities);
1553 std::size_t n = 0;
1554 for(std::size_t i = 0; i < entities.size(); i++)
1555 n += entities[i]->getNumMeshParentElements();
1556 return n;
1557 }
1558
addMEdge(MEdge & edge,std::size_t num)1559 std::size_t GModel::addMEdge(MEdge &edge, std::size_t num)
1560 {
1561 std::pair<MEdge, std::size_t> key(std::move(edge),
1562 num ? num : _mapEdgeNum.size() + 1);
1563 std::pair<hashmapMEdge::iterator, bool> it = _mapEdgeNum.insert(std::move(key));
1564 return it.first->second;
1565 }
1566
getMEdge(MVertex * v0,MVertex * v1,MEdge & edge)1567 std::size_t GModel::getMEdge(MVertex *v0, MVertex *v1, MEdge &edge)
1568 {
1569 MEdge e(v0, v1);
1570 auto it = _mapEdgeNum.find(e);
1571 if(it != _mapEdgeNum.end()) {
1572 edge = it->first;
1573 return it->second;
1574 }
1575 else {
1576 Msg::Error("Unknown edge %d %d", v0->getNum(), v1->getNum());
1577 return 0;
1578 }
1579 }
1580
addMFace(MFace & face,std::size_t num)1581 std::size_t GModel::addMFace(MFace &face, std::size_t num)
1582 {
1583 std::pair<MFace, std::size_t> key(std::move(face),
1584 num ? num : _mapFaceNum.size() + 1);
1585 std::pair<hashmapMFace::iterator, bool> it = _mapFaceNum.insert(std::move(key));
1586 return it.first->second;
1587 }
1588
getMFace(MVertex * v0,MVertex * v1,MVertex * v2,MVertex * v3,MFace & face)1589 std::size_t GModel::getMFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
1590 MFace &face)
1591 {
1592 MFace f(v0, v1, v2, v3);
1593 auto it = _mapFaceNum.find(f);
1594 if(it != _mapFaceNum.end()) {
1595 face = it->first;
1596 return it->second;
1597 }
1598 else {
1599 Msg::Error("Unknown face %d %d %d", v0->getNum(), v1->getNum(), v2->getNum());
1600 return 0;
1601 }
1602 }
1603
1604 #if defined(HAVE_POST)
getDependentViewData(GModel * m,PViewDataGModel::DataType type,std::vector<stepData<double> * > & data)1605 static void getDependentViewData(GModel *m, PViewDataGModel::DataType type,
1606 std::vector<stepData<double> *> &data)
1607 {
1608 for(std::size_t i = 0; i < PView::list.size(); i++) {
1609 PViewDataGModel *d =
1610 dynamic_cast<PViewDataGModel *>(PView::list[i]->getData());
1611 if(d && d->getType() == type) {
1612 for(int step = 0; step < d->getNumTimeSteps(); step++) {
1613 if(d->getStepData(step)->getModel() == m)
1614 data.push_back(d->getStepData(step));
1615 }
1616 }
1617 }
1618 }
1619 #endif
1620
renumberMeshVertices()1621 void GModel::renumberMeshVertices()
1622 {
1623 destroyMeshCaches();
1624 setMaxVertexNumber(CTX::instance()->mesh.firstNodeTag - 1);
1625 std::vector<GEntity *> entities;
1626 getEntities(entities);
1627
1628 #if defined(HAVE_POST)
1629 // check if any nodal post-processing datasets depend on the model; if so,
1630 // keep track of the old node numbering
1631 std::map<MVertex *, int> old;
1632 std::vector<stepData<double> *> data;
1633 getDependentViewData(this, PViewDataGModel::NodeData, data);
1634 if(data.size()) {
1635 for(std::size_t i = 0; i < entities.size(); i++) {
1636 GEntity *ge = entities[i];
1637 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1638 MVertex *v = ge->getMeshVertex(j);
1639 old[v] = v->getNum();
1640 }
1641 }
1642 }
1643 #endif
1644
1645 // check if we will potentially only save a subset of elements, i.e. those
1646 // belonging to physical groups
1647 bool potentiallySaveSubset = false;
1648 if(!CTX::instance()->mesh.saveAll) {
1649 for(std::size_t i = 0; i < entities.size(); i++) {
1650 if(entities[i]->physicals.size()) {
1651 potentiallySaveSubset = true;
1652 break;
1653 }
1654 }
1655 }
1656
1657 std::size_t n = CTX::instance()->mesh.firstNodeTag - 1;
1658 if(potentiallySaveSubset) {
1659 Msg::Debug("Renumbering for potentially partial mesh save");
1660 // if we potentially only save a subset of elements, make sure to first
1661 // renumber the nodes that belong to those elements (so that we end up
1662 // with a dense node numbering in the output file)
1663 std::size_t nv = CTX::instance()->mesh.firstNodeTag - 1;
1664 for(std::size_t i = 0; i < entities.size(); i++) {
1665 GEntity *ge = entities[i];
1666 nv += ge->getNumMeshVertices();
1667 }
1668 for(std::size_t i = 0; i < entities.size(); i++) {
1669 GEntity *ge = entities[i];
1670 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1671 ge->getMeshVertex(j)->forceNum(nv + 1);
1672 }
1673 }
1674 for(std::size_t i = 0; i < entities.size(); i++) {
1675 GEntity *ge = entities[i];
1676 if(ge->physicals.size()) {
1677 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1678 MElement *e = ge->getMeshElement(j);
1679 for(std::size_t k = 0; k < e->getNumVertices(); k++) {
1680 e->getVertex(k)->forceNum(0);
1681 }
1682 }
1683 }
1684 }
1685 for(std::size_t i = 0; i < entities.size(); i++) {
1686 GEntity *ge = entities[i];
1687 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1688 MVertex *v = ge->getMeshVertex(j);
1689 if(v->getNum() == 0) v->forceNum(++n);
1690 }
1691 }
1692 for(std::size_t i = 0; i < entities.size(); i++) {
1693 GEntity *ge = entities[i];
1694 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1695 MVertex *v = ge->getMeshVertex(j);
1696 if(v->getNum() == nv + 1) v->forceNum(++n);
1697 }
1698 }
1699 }
1700 else {
1701 // no physical groups
1702 for(std::size_t i = 0; i < entities.size(); i++) {
1703 GEntity *ge = entities[i];
1704 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1705 ge->getMeshVertex(j)->forceNum(++n);
1706 }
1707 }
1708 }
1709
1710 #if defined(HAVE_POST)
1711 // renumber any dependent nodal post-processing datasets
1712 if(data.size()) {
1713 int n = data.size();
1714 Msg::Info("Renumbering nodal model data (%d step%s)", n, n > 1 ? "s" : "");
1715 std::map<int, int> remap;
1716 for(std::size_t i = 0; i < entities.size(); i++) {
1717 GEntity *ge = entities[i];
1718 for(std::size_t j = 0; j < ge->getNumMeshVertices(); j++) {
1719 MVertex *v = ge->getMeshVertex(j);
1720 remap[old[v]] = v->getNum();
1721 }
1722 }
1723 for(auto d : data) { d->renumberData(remap); }
1724 }
1725 #endif
1726 }
1727
renumberMeshElements()1728 void GModel::renumberMeshElements()
1729 {
1730 destroyMeshCaches();
1731 setMaxElementNumber(CTX::instance()->mesh.firstElementTag - 1);
1732 std::vector<GEntity *> entities;
1733 getEntities(entities);
1734
1735 #if defined(HAVE_POST)
1736 // check if any element-based post-processing datasets depend on the model; if
1737 // so, keep track of the old element numbering
1738 std::map<MElement *, int> old;
1739 std::vector<stepData<double> *> data[2];
1740 getDependentViewData(this, PViewDataGModel::ElementData, data[0]);
1741 getDependentViewData(this, PViewDataGModel::ElementNodeData, data[1]);
1742 if(data[0].size() || data[1].size()) {
1743 for(std::size_t i = 0; i < entities.size(); i++) {
1744 GEntity *ge = entities[i];
1745 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1746 MElement *e = ge->getMeshElement(j);
1747 old[e] = e->getNum();
1748 }
1749 }
1750 }
1751 #endif
1752
1753 // check if we will potentially only save a subset of elements, i.e. those
1754 // belonging to physical groups
1755 bool potentiallySaveSubset = false;
1756 if(!CTX::instance()->mesh.saveAll) {
1757 for(std::size_t i = 0; i < entities.size(); i++) {
1758 if(entities[i]->physicals.size()) {
1759 potentiallySaveSubset = true;
1760 break;
1761 }
1762 }
1763 }
1764
1765 std::size_t n = CTX::instance()->mesh.firstElementTag - 1;
1766 if(potentiallySaveSubset) {
1767 for(std::size_t i = 0; i < entities.size(); i++) {
1768 GEntity *ge = entities[i];
1769 if(ge->physicals.size()) {
1770 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1771 ge->getMeshElement(j)->forceNum(++n);
1772 }
1773 }
1774 }
1775 for(std::size_t i = 0; i < entities.size(); i++) {
1776 GEntity *ge = entities[i];
1777 if(ge->physicals.empty()) {
1778 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1779 ge->getMeshElement(j)->forceNum(++n);
1780 }
1781 }
1782 }
1783 }
1784 else {
1785 for(std::size_t i = 0; i < entities.size(); i++) {
1786 GEntity *ge = entities[i];
1787 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1788 ge->getMeshElement(j)->forceNum(++n);
1789 }
1790 }
1791 }
1792
1793 #if defined(HAVE_POST)
1794 // renumber any dependent element post-processing datasets
1795 if(data[0].size() || data[1].size()) {
1796 int n = data[0].size() + data[1].size();
1797 Msg::Info("Renumbering element model data (%d step%s)", n,
1798 n > 1 ? "s" : "");
1799 std::map<int, int> remap;
1800 for(std::size_t i = 0; i < entities.size(); i++) {
1801 GEntity *ge = entities[i];
1802 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
1803 MElement *e = ge->getMeshElement(j);
1804 remap[old[e]] = e->getNum();
1805 }
1806 }
1807 for(int i = 0; i < 2; i++) {
1808 for(auto d : data[i]) { d->renumberData(remap); }
1809 }
1810 }
1811 #endif
1812 }
1813
getNumMeshElements(unsigned c[6])1814 std::size_t GModel::getNumMeshElements(unsigned c[6])
1815 {
1816 c[0] = 0;
1817 c[1] = 0;
1818 c[2] = 0;
1819 c[3] = 0;
1820 c[4] = 0;
1821 c[5] = 0;
1822 for(auto it = firstRegion(); it != lastRegion(); ++it)
1823 (*it)->getNumMeshElements(c);
1824 if(c[0] + c[1] + c[2] + c[3] + c[4] + c[5]) return 3;
1825 for(auto it = firstFace(); it != lastFace(); ++it)
1826 (*it)->getNumMeshElements(c);
1827 if(c[0] + c[1] + c[2]) return 2;
1828 for(auto it = firstEdge(); it != lastEdge(); ++it)
1829 (*it)->getNumMeshElements(c);
1830 if(c[0]) return 1;
1831 return 0;
1832 }
1833
getMeshElementByCoord(SPoint3 & p,SPoint3 & param,int dim,bool strict)1834 MElement *GModel::getMeshElementByCoord(SPoint3 &p, SPoint3 ¶m, int dim,
1835 bool strict)
1836 {
1837 if(!_elementOctree) {
1838 #pragma omp barrier
1839 #pragma omp single
1840 {
1841 Msg::Debug("Rebuilding mesh element octree");
1842 _elementOctree = new MElementOctree(this);
1843 }
1844 }
1845 MElement *e = _elementOctree->find(p.x(), p.y(), p.z(), dim, strict);
1846 if(e) {
1847 double xyz[3] = {p.x(), p.y(), p.z()}, uvw[3];
1848 e->xyz2uvw(xyz, uvw);
1849 param.setPosition(uvw[0], uvw[1], uvw[2]);
1850 }
1851 else {
1852 param.setPosition(0, 0, 0);
1853 }
1854 return e;
1855 }
1856
getMeshElementsByCoord(SPoint3 & p,int dim,bool strict)1857 std::vector<MElement *> GModel::getMeshElementsByCoord(SPoint3 &p, int dim,
1858 bool strict)
1859 {
1860 if(!_elementOctree) {
1861 #pragma omp barrier
1862 #pragma omp single
1863 {
1864 Msg::Debug("Rebuilding mesh element octree");
1865 _elementOctree = new MElementOctree(this);
1866 }
1867 }
1868 return _elementOctree->findAll(p.x(), p.y(), p.z(), dim, strict);
1869 }
1870
rebuildMeshVertexCache(bool onlyIfNecessary)1871 void GModel::rebuildMeshVertexCache(bool onlyIfNecessary)
1872 {
1873 if(!onlyIfNecessary ||
1874 (_vertexVectorCache.empty() && _vertexMapCache.empty())) {
1875 _vertexVectorCache.clear();
1876 _vertexMapCache.clear();
1877 bool dense = false;
1878 if(_maxVertexNum == getNumMeshVertices()) {
1879 Msg::Debug("We have a dense node numbering in the cache");
1880 dense = true;
1881 }
1882 else if(_maxVertexNum < 10 * getNumMeshVertices()) {
1883 Msg::Debug(
1884 "We have a fairly dense node numbering - still using cache vector");
1885 dense = true;
1886 }
1887 std::vector<GEntity *> entities;
1888 getEntities(entities);
1889 if(dense) {
1890 // numbering starts at 1
1891 _vertexVectorCache.resize(_maxVertexNum + 1, (MVertex *)nullptr);
1892 for(std::size_t i = 0; i < entities.size(); i++)
1893 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
1894 _vertexVectorCache[entities[i]->mesh_vertices[j]->getNum()] =
1895 entities[i]->mesh_vertices[j];
1896 }
1897 else {
1898 for(std::size_t i = 0; i < entities.size(); i++)
1899 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
1900 _vertexMapCache[entities[i]->mesh_vertices[j]->getNum()] =
1901 entities[i]->mesh_vertices[j];
1902 }
1903 }
1904 }
1905
rebuildMeshElementCache(bool onlyIfNecessary)1906 void GModel::rebuildMeshElementCache(bool onlyIfNecessary)
1907 {
1908 if(!onlyIfNecessary ||
1909 (_elementVectorCache.empty() && _elementMapCache.empty())) {
1910 Msg::Debug("Rebuilding mesh element cache");
1911 _elementVectorCache.clear();
1912 _elementMapCache.clear();
1913 bool dense = false;
1914 if(_maxElementNum == getNumMeshElements()) {
1915 Msg::Debug("We have a dense element numbering in the cache");
1916 dense = true;
1917 }
1918 else if(_maxElementNum < 10 * getNumMeshElements()) {
1919 Msg::Debug(
1920 "We have a fairly dense element numbering - still using cache vector");
1921 dense = true;
1922 }
1923 std::vector<GEntity *> entities;
1924 getEntities(entities);
1925 if(dense) {
1926 // numbering starts at 1
1927 _elementVectorCache.resize(_maxElementNum + 1, std::make_pair(nullptr, 0));
1928 for(std::size_t i = 0; i < entities.size(); i++)
1929 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
1930 MElement *e = entities[i]->getMeshElement(j);
1931 _elementVectorCache[e->getNum()] =
1932 std::make_pair(e, entities[i]->tag());
1933 }
1934 }
1935 else {
1936 for(std::size_t i = 0; i < entities.size(); i++)
1937 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
1938 MElement *e = entities[i]->getMeshElement(j);
1939 _elementMapCache[e->getNum()] = std::make_pair(e, entities[i]->tag());
1940 }
1941 }
1942 }
1943 }
1944
getMeshVertexByTag(int n)1945 MVertex *GModel::getMeshVertexByTag(int n)
1946 {
1947 if(_vertexVectorCache.empty() && _vertexMapCache.empty()) {
1948 #pragma omp barrier
1949 #pragma omp single
1950 {
1951 Msg::Debug("Rebuilding mesh node cache");
1952 rebuildMeshVertexCache();
1953 }
1954 }
1955
1956 if(n < (int)_vertexVectorCache.size())
1957 return _vertexVectorCache[n];
1958 else
1959 return _vertexMapCache[n];
1960 }
1961
addMVertexToVertexCache(MVertex * v)1962 void GModel::addMVertexToVertexCache(MVertex* v)
1963 {
1964 if(_vertexVectorCache.empty() && _vertexMapCache.empty()) {
1965 Msg::Debug("Rebuilding mesh node cache");
1966 rebuildMeshVertexCache();
1967 }
1968 if (_vertexVectorCache.size() > 0) {
1969 #pragma omp critical
1970 if (v->getNum() >= _vertexVectorCache.size()) {
1971 _vertexVectorCache.resize(v->getNum()+1, nullptr);
1972 }
1973 _vertexVectorCache[v->getNum()] = v;
1974 } else {
1975 _vertexMapCache[v->getNum()] = v;
1976 }
1977 }
1978
getMeshVerticesForPhysicalGroup(int dim,int num,std::vector<MVertex * > & v)1979 void GModel::getMeshVerticesForPhysicalGroup(int dim, int num,
1980 std::vector<MVertex *> &v)
1981 {
1982 v.clear();
1983 std::map<int, std::vector<GEntity *> > groups;
1984 getPhysicalGroups(dim, groups);
1985 std::map<int, std::vector<GEntity *> >::const_iterator it = groups.find(num);
1986 if(it == groups.end()) return;
1987 const std::vector<GEntity *> &entities = it->second;
1988 std::set<MVertex *, MVertexPtrLessThan> sv;
1989 for(std::size_t i = 0; i < entities.size(); i++) {
1990 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
1991 MElement *e = entities[i]->getMeshElement(j);
1992 for(std::size_t k = 0; k < e->getNumVertices(); k++)
1993 sv.insert(e->getVertex(k));
1994 }
1995 }
1996 v.insert(v.begin(), sv.begin(), sv.end());
1997 }
1998
getMeshElementByTag(int n,int & entityTag)1999 MElement *GModel::getMeshElementByTag(int n, int &entityTag)
2000 {
2001 if(_elementVectorCache.empty() && _elementMapCache.empty()) {
2002 #pragma omp barrier
2003 #pragma omp single
2004 {
2005 Msg::Debug("Rebuilding mesh element cache");
2006 rebuildMeshElementCache();
2007 }
2008 }
2009
2010 std::pair<MElement*, int> ret;
2011 if(n < (int)_elementVectorCache.size())
2012 ret = _elementVectorCache[n];
2013 else
2014 ret = _elementMapCache[n];
2015 entityTag = ret.second;
2016 return ret.first;
2017 }
2018
getMeshElementIndex(MElement * e)2019 int GModel::getMeshElementIndex(MElement *e)
2020 {
2021 if(!e) return 0;
2022 if(_elementIndexCache.empty()) return e->getNum();
2023 auto it = _elementIndexCache.find(e->getNum());
2024 if(it != _elementIndexCache.end()) return it->second;
2025 return e->getNum();
2026 }
2027
setMeshElementIndex(MElement * e,int index)2028 void GModel::setMeshElementIndex(MElement *e, int index)
2029 {
2030 _elementIndexCache[e->getNum()] = index;
2031 }
2032
2033 template <class T>
removeInvisible(std::vector<T * > & elements,bool all)2034 static std::size_t removeInvisible(std::vector<T *> &elements, bool all)
2035 {
2036 std::vector<T *> tmp;
2037 std::size_t n = elements.size();
2038 for(std::size_t i = 0; i < n; i++) {
2039 if(all || !elements[i]->getVisibility())
2040 delete elements[i];
2041 else
2042 tmp.push_back(elements[i]);
2043 }
2044 elements.clear();
2045 elements = tmp;
2046 return n - elements.size();
2047 }
2048
removeInvisibleElements()2049 std::size_t GModel::removeInvisibleElements()
2050 {
2051 std::size_t n = 0;
2052 for(auto it = firstRegion(); it != lastRegion(); ++it) {
2053 bool all = !(*it)->getVisibility();
2054 n += removeInvisible((*it)->tetrahedra, all);
2055 n += removeInvisible((*it)->hexahedra, all);
2056 n += removeInvisible((*it)->prisms, all);
2057 n += removeInvisible((*it)->pyramids, all);
2058 n += removeInvisible((*it)->trihedra, all);
2059 n += removeInvisible((*it)->polyhedra, all);
2060 (*it)->deleteVertexArrays();
2061 }
2062 for(auto it = firstFace(); it != lastFace(); ++it) {
2063 bool all = !(*it)->getVisibility();
2064 n += removeInvisible((*it)->triangles, all);
2065 n += removeInvisible((*it)->quadrangles, all);
2066 n += removeInvisible((*it)->polygons, all);
2067 (*it)->deleteVertexArrays();
2068 }
2069 for(auto it = firstEdge(); it != lastEdge(); ++it) {
2070 bool all = !(*it)->getVisibility();
2071 n += removeInvisible((*it)->lines, all);
2072 (*it)->deleteVertexArrays();
2073 }
2074 destroyMeshCaches();
2075 Msg::Info("Removed %lu elements", n);
2076 return n;
2077 }
2078
2079 template <class T>
reverseInvisible(std::vector<T * > & elements,bool all)2080 static std::size_t reverseInvisible(std::vector<T *> &elements, bool all)
2081 {
2082 std::size_t n = 0;
2083 std::vector<T *> tmp;
2084 for(std::size_t i = 0; i < elements.size(); i++) {
2085 if(all || !elements[i]->getVisibility()) {
2086 elements[i]->reverse();
2087 elements[i]->setVisibility(1);
2088 n++;
2089 }
2090 }
2091 return n;
2092 }
2093
reverseInvisibleElements()2094 std::size_t GModel::reverseInvisibleElements()
2095 {
2096 std::size_t n = 0;
2097 for(auto it = firstRegion(); it != lastRegion(); ++it) {
2098 bool all = !(*it)->getVisibility();
2099 n += reverseInvisible((*it)->tetrahedra, all);
2100 n += reverseInvisible((*it)->hexahedra, all);
2101 n += reverseInvisible((*it)->prisms, all);
2102 n += reverseInvisible((*it)->pyramids, all);
2103 n += reverseInvisible((*it)->trihedra, all);
2104 n += reverseInvisible((*it)->polyhedra, all);
2105 (*it)->deleteVertexArrays();
2106 if(all) (*it)->setVisibility(1);
2107 }
2108 for(auto it = firstFace(); it != lastFace(); ++it) {
2109 bool all = !(*it)->getVisibility();
2110 n += reverseInvisible((*it)->triangles, all);
2111 n += reverseInvisible((*it)->quadrangles, all);
2112 n += reverseInvisible((*it)->polygons, all);
2113 (*it)->deleteVertexArrays();
2114 if(all) (*it)->setVisibility(1);
2115 }
2116 for(auto it = firstEdge(); it != lastEdge(); ++it) {
2117 bool all = !(*it)->getVisibility();
2118 n += reverseInvisible((*it)->lines, all);
2119 (*it)->deleteVertexArrays();
2120 if(all) (*it)->setVisibility(1);
2121 }
2122 destroyMeshCaches();
2123 Msg::Info("Reversed %lu elements", n);
2124 return n;
2125 }
2126
indexMeshVertices(bool all,int singlePartition,bool renumber)2127 std::size_t GModel::indexMeshVertices(bool all, int singlePartition,
2128 bool renumber)
2129 {
2130 std::vector<GEntity *> entities;
2131 getEntities(entities);
2132
2133 // tag all mesh nodes with -1 (negative nodes will not be saved)
2134 for(std::size_t i = 0; i < entities.size(); i++)
2135 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
2136 entities[i]->mesh_vertices[j]->setIndex(-1);
2137
2138 // tag all mesh nodes belonging to elements that need to be saved with 0, or
2139 // with -2 if they need to be taken into account in the numbering but need not
2140 // to be saved (because we save a single partition and they are not used in
2141 // that partition)
2142 for(std::size_t i = 0; i < entities.size(); i++) {
2143 if(all || entities[i]->physicals.size() ||
2144 (entities[i]->getParentEntity() &&
2145 entities[i]->getParentEntity()->physicals.size())) {
2146 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
2147 MElement *e = entities[i]->getMeshElement(j);
2148 for(std::size_t k = 0; k < e->getNumVertices(); k++) {
2149 if(singlePartition <= 0 || e->getPartition() == singlePartition)
2150 e->getVertex(k)->setIndex(0);
2151 else if(e->getVertex(k)->getIndex() == -1)
2152 e->getVertex(k)->setIndex(-2);
2153 }
2154 }
2155 }
2156 }
2157
2158 // renumber all the mesh nodes tagged with 0
2159 std::size_t numVertices = 0;
2160 long int index = 0;
2161 for(std::size_t i = 0; i < entities.size(); i++) {
2162 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++) {
2163 MVertex *v = entities[i]->mesh_vertices[j];
2164 if(!v->getIndex()) {
2165 index++;
2166 numVertices++;
2167 if(renumber)
2168 v->setIndex(index);
2169 else
2170 v->setIndex(v->getNum());
2171 }
2172 else if(v->getIndex() == -2) {
2173 index++;
2174 }
2175 }
2176 }
2177
2178 return numVertices;
2179 }
2180
scaleMesh(double factor)2181 void GModel::scaleMesh(double factor)
2182 {
2183 std::vector<GEntity *> entities;
2184 getEntities(entities);
2185 for(std::size_t i = 0; i < entities.size(); i++)
2186 for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++) {
2187 MVertex *v = entities[i]->mesh_vertices[j];
2188 v->x() *= factor;
2189 v->y() *= factor;
2190 v->z() *= factor;
2191 }
2192 }
2193
setCurrentMeshEntity(GEntity * e)2194 void GModel::setCurrentMeshEntity(GEntity *e)
2195 {
2196 #pragma omp atomic write
2197 _currentMeshEntity = e;
2198 }
2199
partitionMesh(int numPart,std::vector<std::pair<MElement *,int>> elementPartition)2200 int GModel::partitionMesh(
2201 int numPart, std::vector<std::pair<MElement *, int> > elementPartition)
2202 {
2203 #if defined(HAVE_MESH) && (defined(HAVE_METIS))
2204 if(numPart > 0) {
2205 if(_numPartitions > 0) UnpartitionMesh(this);
2206 if(elementPartition.empty())
2207 return PartitionMesh(this, numPart);
2208 else
2209 return PartitionUsingThisSplit(this, elementPartition);
2210 }
2211 return 1;
2212 #else
2213 Msg::Error("Mesh or Metis module not compiled");
2214 return 1;
2215 #endif
2216 }
2217
unpartitionMesh()2218 int GModel::unpartitionMesh()
2219 {
2220 #if defined(HAVE_MESH)
2221 return UnpartitionMesh(this);
2222 #else
2223 Msg::Error("Mesh module not compiled");
2224 return 1;
2225 #endif
2226 }
2227
convertOldPartitioningToNewOne()2228 int GModel::convertOldPartitioningToNewOne()
2229 {
2230 #if defined(HAVE_MESH) && (defined(HAVE_METIS))
2231 int ier = ConvertOldPartitioningToNewOne(this);
2232 return ier;
2233 #else
2234 Msg::Error("Mesh or Metis module not compiled");
2235 return 1;
2236 #endif
2237 }
2238
storeChain(int dim,std::map<int,std::vector<MElement * >> & entityMap,std::map<int,std::map<int,std::string>> & physicalMap)2239 void GModel::storeChain(int dim,
2240 std::map<int, std::vector<MElement *> > &entityMap,
2241 std::map<int, std::map<int, std::string> > &physicalMap)
2242 {
2243 _storeElementsInEntities(entityMap);
2244 _storePhysicalTagsInEntities(dim, physicalMap);
2245 std::map<int, std::vector<MElement *> >::iterator it;
2246 for(it = entityMap.begin(); it != entityMap.end(); it++) {
2247 if(dim == 0)
2248 _chainVertices.insert(getVertexByTag(it->first));
2249 else if(dim == 1)
2250 _chainEdges.insert(getEdgeByTag(it->first));
2251 else if(dim == 2)
2252 _chainFaces.insert(getFaceByTag(it->first));
2253 else if(dim == 3)
2254 _chainRegions.insert(getRegionByTag(it->first));
2255 }
2256 }
2257
2258 template <class T>
_addElements(std::vector<T * > & dst,const std::vector<MElement * > & src)2259 static void _addElements(std::vector<T *> &dst,
2260 const std::vector<MElement *> &src)
2261 {
2262 for(std::size_t i = 0; i < src.size(); i++) dst.push_back((T *)src[i]);
2263 }
2264
_storeElementsInEntities(std::map<int,std::vector<MElement * >> & map)2265 void GModel::_storeElementsInEntities(
2266 std::map<int, std::vector<MElement *> > &map)
2267 {
2268 std::map<int, std::vector<MElement *> >::const_iterator it = map.begin();
2269 for(; it != map.end(); ++it) {
2270 if(!it->second.size()) continue;
2271 int type = it->second[0]->getType();
2272 switch(type) {
2273 case TYPE_PNT: {
2274 GVertex *v = getVertexByTag(it->first);
2275 if(!v) {
2276 double x = it->second[0]->getVertex(0)->x();
2277 double y = it->second[0]->getVertex(0)->y();
2278 double z = it->second[0]->getVertex(0)->z();
2279 v = new discreteVertex(this, it->first, x, y, z);
2280 add(v);
2281 }
2282 if(!v->points.empty()) { // CAD points already have one by default
2283 v->points.clear();
2284 v->mesh_vertices.clear();
2285 }
2286 _addElements(v->points, it->second);
2287 } break;
2288 case TYPE_LIN: {
2289 GEdge *e = getEdgeByTag(it->first);
2290 if(!e) {
2291 e = new discreteEdge(this, it->first, nullptr, nullptr);
2292 add(e);
2293 }
2294 _addElements(e->lines, it->second);
2295 } break;
2296 case TYPE_TRI:
2297 case TYPE_QUA:
2298 case TYPE_POLYG: {
2299 GFace *f = getFaceByTag(it->first);
2300 if(!f) {
2301 f = new discreteFace(this, it->first);
2302 add(f);
2303 }
2304 if(type == TYPE_TRI)
2305 _addElements(f->triangles, it->second);
2306 else if(type == TYPE_QUA)
2307 _addElements(f->quadrangles, it->second);
2308 else
2309 _addElements(f->polygons, it->second);
2310 } break;
2311 case TYPE_TET:
2312 case TYPE_HEX:
2313 case TYPE_PYR:
2314 case TYPE_TRIH:
2315 case TYPE_PRI:
2316 case TYPE_POLYH: {
2317 GRegion *r = getRegionByTag(it->first);
2318 if(!r) {
2319 r = new discreteRegion(this, it->first);
2320 add(r);
2321 }
2322 if(type == TYPE_TET)
2323 _addElements(r->tetrahedra, it->second);
2324 else if(type == TYPE_HEX)
2325 _addElements(r->hexahedra, it->second);
2326 else if(type == TYPE_PRI)
2327 _addElements(r->prisms, it->second);
2328 else if(type == TYPE_PYR)
2329 _addElements(r->pyramids, it->second);
2330 else if(type == TYPE_TRIH)
2331 _addElements(r->trihedra, it->second);
2332 else
2333 _addElements(r->polyhedra, it->second);
2334 } break;
2335 }
2336 }
2337 }
2338
_storeParentsInSubElements(std::map<int,std::vector<MElement * >> & map)2339 void GModel::_storeParentsInSubElements(
2340 std::map<int, std::vector<MElement *> > &map)
2341 {
2342 std::map<int, std::vector<MElement *> >::const_iterator it;
2343 for(it = map.begin(); it != map.end(); ++it)
2344 for(std::size_t i = 0; i < it->second.size(); ++i)
2345 it->second[i]->updateParent(this);
2346 }
2347
2348 template <class T>
_associateEntityWithElementVertices(GEntity * ge,std::vector<T * > & elements,bool force=false)2349 static void _associateEntityWithElementVertices(GEntity *ge,
2350 std::vector<T *> &elements,
2351 bool force = false)
2352 {
2353 for(std::size_t i = 0; i < elements.size(); i++) {
2354 for(std::size_t j = 0; j < elements[i]->getNumVertices(); j++) {
2355 if(force || !elements[i]->getVertex(j)->onWhat() ||
2356 elements[i]->getVertex(j)->onWhat()->dim() > ge->dim())
2357 elements[i]->getVertex(j)->setEntity(ge);
2358 }
2359 }
2360 }
2361
createGeometryOfDiscreteEntities(const std::vector<std::pair<int,int>> & dimTags)2362 void GModel::createGeometryOfDiscreteEntities(
2363 const std::vector<std::pair<int, int> > &dimTags)
2364 {
2365 std::vector<discreteEdge *> e;
2366 std::vector<discreteFace *> f;
2367 std::vector<discreteRegion *> r;
2368
2369 if(dimTags.empty()) {
2370 for(auto it = firstEdge(); it != lastEdge(); it++) {
2371 discreteEdge *de = dynamic_cast<discreteEdge *>(*it);
2372 if(de) e.push_back(de);
2373 }
2374 for(auto it = firstFace(); it != lastFace(); it++) {
2375 discreteFace *df = dynamic_cast<discreteFace *>(*it);
2376 if(df) f.push_back(df);
2377 }
2378 for(auto it = firstRegion(); it != lastRegion(); it++) {
2379 discreteRegion *dr = dynamic_cast<discreteRegion *>(*it);
2380 if(dr) r.push_back(dr);
2381 }
2382 }
2383 else {
2384 for(std::size_t i = 0; i < dimTags.size(); i++) {
2385 int dim = dimTags[i].first;
2386 int tag = dimTags[i].second;
2387 if(dim == 1) {
2388 discreteEdge *de = dynamic_cast<discreteEdge *>(getEdgeByTag(tag));
2389 if(de)
2390 e.push_back(de);
2391 else
2392 Msg::Error("No discrete curve with tag %d", tag);
2393 }
2394 else if(dim == 2) {
2395 discreteFace *df = dynamic_cast<discreteFace *>(getFaceByTag(tag));
2396 if(df)
2397 f.push_back(df);
2398 else
2399 Msg::Error("No discrete surface with tag %d", tag);
2400 }
2401 else if(dim == 3) {
2402 discreteRegion *dr =
2403 dynamic_cast<discreteRegion *>(getRegionByTag(tag));
2404 if(dr)
2405 r.push_back(dr);
2406 else
2407 Msg::Error("No discrete volume with tag %d", tag);
2408 }
2409 }
2410 }
2411
2412 if(e.size()) {
2413 Msg::StatusBar(true, "Creating geometry of discrete curves...");
2414 double t1 = Cpu(), w1 = TimeOfDay();
2415 for(std::size_t i = 0; i < e.size(); i++) {
2416 if(e[i]->createGeometry())
2417 Msg::Error("Could not create geometry of discrete curve %d",
2418 e[i]->tag());
2419 }
2420 double t2 = Cpu(), w2 = TimeOfDay();
2421 Msg::StatusBar(true,
2422 "Done creating geometry of discrete curves "
2423 "(Wall %gs, CPU %gs)",
2424 w2 - w1, t2 - t1);
2425 }
2426
2427 if(f.size()) {
2428 Msg::StatusBar(true, "Creating geometry of discrete surfaces...");
2429 double t1 = Cpu(), w1 = TimeOfDay();
2430 Msg::StartProgressMeter(f.size());
2431 for(std::size_t i = 0; i < f.size(); i++) {
2432 Msg::ProgressMeter(i, true, "Creating geometry");
2433 if(f[i]->createGeometry())
2434 Msg::Error("Could not create geometry of discrete surface %d",
2435 f[i]->tag());
2436 }
2437 Msg::StopProgressMeter();
2438 double t2 = Cpu();
2439 double w2 = TimeOfDay();
2440 Msg::StatusBar(true,
2441 "Done creating geometry of discrete surfaces "
2442 "(Wall %gs, CPU %gs)",
2443 w2 - w1, t2 - t1);
2444 }
2445
2446 if(r.size()) {
2447 Msg::StatusBar(true, "Creating geometry of discrete volumes...");
2448 double t1 = Cpu(), w1 = TimeOfDay();
2449 for(std::size_t i = 0; i < r.size(); i++) {
2450 if(r[i]->createGeometry())
2451 Msg::Error("Could not create geometry of discrete volume %d",
2452 r[i]->tag());
2453 }
2454 double t2 = Cpu(), w2 = TimeOfDay();
2455 Msg::StatusBar(true,
2456 "Done creating geometry of discrete volumes "
2457 "(Wall %gs, CPU %gs)",
2458 w2 - w1, t2 - t1);
2459 }
2460 }
2461
_associateEntityWithMeshVertices()2462 void GModel::_associateEntityWithMeshVertices()
2463 {
2464 // loop on regions, then on faces, edges and vertices and store the entity
2465 // pointer in the the elements' vertices (this way we associate the entity of
2466 // lowest geometrical degree with each vertex)
2467 for(auto it = firstRegion(); it != lastRegion(); ++it) {
2468 _associateEntityWithElementVertices(*it, (*it)->tetrahedra);
2469 _associateEntityWithElementVertices(*it, (*it)->hexahedra);
2470 _associateEntityWithElementVertices(*it, (*it)->prisms);
2471 _associateEntityWithElementVertices(*it, (*it)->pyramids);
2472 _associateEntityWithElementVertices(*it, (*it)->trihedra);
2473 _associateEntityWithElementVertices(*it, (*it)->polyhedra);
2474 }
2475 for(auto it = firstFace(); it != lastFace(); ++it) {
2476 _associateEntityWithElementVertices(*it, (*it)->triangles);
2477 _associateEntityWithElementVertices(*it, (*it)->quadrangles);
2478 _associateEntityWithElementVertices(*it, (*it)->polygons);
2479 }
2480 for(auto it = firstEdge(); it != lastEdge(); ++it) {
2481 _associateEntityWithElementVertices(*it, (*it)->lines);
2482 }
2483 for(auto it = firstVertex(); it != lastVertex(); ++it) {
2484 _associateEntityWithElementVertices(*it, (*it)->points);
2485 }
2486 }
2487
_storeVerticesInEntities(std::map<int,MVertex * > & vertices)2488 void GModel::_storeVerticesInEntities(std::map<int, MVertex *> &vertices)
2489 {
2490 auto it = vertices.begin();
2491 for(; it != vertices.end(); ++it) {
2492 MVertex *v = it->second;
2493 GEntity *ge = v->onWhat();
2494 if(ge)
2495 ge->mesh_vertices.push_back(v);
2496 else {
2497 delete v; // we delete all unused vertices
2498 it->second = 0;
2499 }
2500 }
2501 }
2502
_storeVerticesInEntities(std::vector<MVertex * > & vertices)2503 void GModel::_storeVerticesInEntities(std::vector<MVertex *> &vertices)
2504 {
2505 for(std::size_t i = 0; i < vertices.size(); i++) {
2506 MVertex *v = vertices[i];
2507 if(v) { // the vector is allowed to have null entries
2508 GEntity *ge = v->onWhat();
2509 if(ge)
2510 ge->mesh_vertices.push_back(v);
2511 else {
2512 delete v; // we delete all unused vertices
2513 vertices[i] = nullptr;
2514 }
2515 }
2516 }
2517 }
2518
pruneMeshVertexAssociations()2519 void GModel::pruneMeshVertexAssociations()
2520 {
2521 std::vector<GEntity *> entities;
2522 std::set<MVertex *, MVertexPtrLessThan> vertSet;
2523 getEntities(entities);
2524 for(std::size_t i = 0; i < entities.size(); i++) {
2525 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
2526 MElement *e = entities[i]->getMeshElement(j);
2527 for(std::size_t k = 0; k < e->getNumVertices(); k++) {
2528 MVertex *v = e->getVertex(k);
2529 v->setEntity(nullptr);
2530 vertSet.insert(v);
2531 }
2532 }
2533 entities[i]->mesh_vertices.clear();
2534 }
2535 std::vector<MVertex *> vertices(vertSet.begin(), vertSet.end());
2536 _associateEntityWithMeshVertices();
2537 // associate mesh nodes primarily with chain entities
2538 for(auto it = _chainRegions.begin(); it != _chainRegions.end(); ++it) {
2539 _associateEntityWithElementVertices(*it, (*it)->tetrahedra, true);
2540 _associateEntityWithElementVertices(*it, (*it)->hexahedra, true);
2541 _associateEntityWithElementVertices(*it, (*it)->prisms, true);
2542 _associateEntityWithElementVertices(*it, (*it)->pyramids, true);
2543 _associateEntityWithElementVertices(*it, (*it)->trihedra, true);
2544 _associateEntityWithElementVertices(*it, (*it)->polyhedra, true);
2545 }
2546 for(auto it = _chainFaces.begin(); it != _chainFaces.end(); ++it) {
2547 _associateEntityWithElementVertices(*it, (*it)->triangles, true);
2548 _associateEntityWithElementVertices(*it, (*it)->quadrangles, true);
2549 _associateEntityWithElementVertices(*it, (*it)->polygons, true);
2550 }
2551 for(auto it = _chainEdges.begin(); it != _chainEdges.end(); ++it) {
2552 _associateEntityWithElementVertices(*it, (*it)->lines, true);
2553 }
2554 for(auto it = _chainVertices.begin(); it != _chainVertices.end(); ++it) {
2555 _associateEntityWithElementVertices(*it, (*it)->points, true);
2556 }
2557 _storeVerticesInEntities(vertices);
2558 }
2559
_storePhysicalTagsInEntities(int dim,std::map<int,std::map<int,std::string>> & map)2560 void GModel::_storePhysicalTagsInEntities(
2561 int dim, std::map<int, std::map<int, std::string> > &map)
2562 {
2563 std::map<int, std::map<int, std::string> >::const_iterator it = map.begin();
2564 for(; it != map.end(); ++it) {
2565 GEntity *ge = nullptr;
2566 switch(dim) {
2567 case 0: ge = getVertexByTag(it->first); break;
2568 case 1: ge = getEdgeByTag(it->first); break;
2569 case 2: ge = getFaceByTag(it->first); break;
2570 case 3: ge = getRegionByTag(it->first); break;
2571 }
2572
2573 if(ge) {
2574 for(auto it2 = it->second.begin(); it2 != it->second.end(); ++it2) {
2575 if(std::find(ge->physicals.begin(), ge->physicals.end(), it2->first) ==
2576 ge->physicals.end()) {
2577 ge->physicals.push_back(it2->first);
2578 }
2579 if(it2->second.size() && it2->second != "unnamed")
2580 _physicalNames.insert(
2581 std::make_pair(std::make_pair(dim, it2->first), it2->second));
2582 }
2583 }
2584 }
2585 }
2586
checkMeshCoherence(double tolerance)2587 void GModel::checkMeshCoherence(double tolerance)
2588 {
2589 int numEle = getNumMeshElements();
2590 if(!numEle) return;
2591
2592 Msg::StatusBar(true, "Checking mesh coherence (%d elements)...", numEle);
2593
2594 SBoundingBox3d bbox = bounds();
2595 double lc = bbox.empty() ? 1. : norm(SVector3(bbox.max(), bbox.min()));
2596 double eps = lc * tolerance;
2597
2598 std::vector<GEntity *> entities;
2599 getEntities(entities);
2600
2601 // check for duplicate mesh nodes
2602 {
2603 Msg::Info("Checking for duplicate nodes...");
2604 std::vector<MVertex *> vertices;
2605 for(std::size_t i = 0; i < entities.size(); i++)
2606 vertices.insert(vertices.end(), entities[i]->mesh_vertices.begin(),
2607 entities[i]->mesh_vertices.end());
2608 MVertexRTree pos(eps);
2609 std::set<MVertex *, MVertexPtrLessThan> duplicates;
2610 int num = pos.insert(vertices, true, &duplicates);
2611 if(num) {
2612 Msg::Error("%d duplicate node%s: see `duplicate_node.pos'", num,
2613 num > 1 ? "s" : "");
2614 FILE *fp = Fopen("duplicate_nodes.pos", "w");
2615 if(fp) {
2616 fprintf(fp, "View \"duplicate vertices\"{\n");
2617 for(auto it = duplicates.begin(); it != duplicates.end(); it++) {
2618 MVertex *v = *it;
2619 fprintf(fp, "SP(%.16g,%.16g,%.16g){%lu};\n", v->x(), v->y(), v->z(),
2620 v->getNum());
2621 }
2622 fprintf(fp, "};\n");
2623 fclose(fp);
2624 }
2625 }
2626 }
2627
2628 // check for duplicate elements and inverted or zero-volume elements
2629 {
2630 Msg::Info("Checking for duplicate elements...");
2631 std::vector<MVertex *> vertices;
2632 for(std::size_t i = 0; i < entities.size(); i++) {
2633 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
2634 MElement *e = entities[i]->getMeshElement(j);
2635 double vol = e->getVolume();
2636 if(vol < 0)
2637 Msg::Warning("Element %d has negative volume", e->getNum());
2638 else if(vol < eps * eps * eps)
2639 Msg::Warning("Element %d has zero volume", e->getNum());
2640 SPoint3 p = e->barycenter();
2641 vertices.push_back(new MVertex(p.x(), p.y(), p.z()));
2642 }
2643 }
2644 MVertexRTree pos(eps);
2645 int num = pos.insert(vertices, true);
2646 for(std::size_t i = 0; i < vertices.size(); i++) delete vertices[i];
2647 if(num) Msg::Error("%d duplicate element%s", num, num > 1 ? "s" : "");
2648 }
2649
2650 // check for isolated nodes (not belonging to any elements
2651 {
2652 Msg::Info("Checking for isolated nodes...");
2653 std::vector<GEntity *> entities2;
2654 getEntities(entities2, getMeshDim());
2655 std::set<MVertex *, MVertexPtrLessThan> allv;
2656 for(std::size_t i = 0; i < entities2.size(); i++) {
2657 for(std::size_t j = 0; j < entities2[i]->getNumMeshElements(); j++) {
2658 MElement *e = entities2[i]->getMeshElement(j);
2659 for(std::size_t k = 0; k < e->getNumVertices(); k++) {
2660 allv.insert(e->getVertex(k));
2661 }
2662 }
2663 }
2664 int diff = (int)(getNumMeshVertices() - allv.size());
2665 if(diff) {
2666 Msg::Warning("%d node%s not connected to any %dD elements", diff,
2667 (diff > 1) ? "s" : "", getMeshDim());
2668 }
2669 }
2670
2671 Msg::StatusBar(true, "Done checking mesh coherence");
2672 }
2673
removeDuplicateMeshVertices(double tolerance)2674 int GModel::removeDuplicateMeshVertices(double tolerance)
2675 {
2676 Msg::StatusBar(true, "Removing duplicate mesh nodes...");
2677
2678 SBoundingBox3d bbox = bounds();
2679 double lc = bbox.empty() ? 1. : norm(SVector3(bbox.max(), bbox.min()));
2680 double eps = lc * tolerance;
2681
2682 // get entities (in order of increasing dimensions so that topological
2683 // classification of vertices remains correct)
2684 std::vector<GEntity *> entities;
2685 getEntities(entities);
2686
2687 // re-index all vertices (don't use MVertex::getNum(), as we want to be able
2688 // to remove duplicate vertices from "incorrect" meshes, where vertices with
2689 // the same number are duplicated)
2690 int n = 0;
2691 for(std::size_t i = 0; i < entities.size(); i++) {
2692 GEntity *ge = entities[i];
2693 for(std::size_t j = 0; j < ge->mesh_vertices.size(); j++) {
2694 MVertex *v = ge->mesh_vertices[j];
2695 v->setIndex(++n);
2696 }
2697 }
2698
2699 MVertexRTree pos(eps);
2700 std::map<int, MVertex *> vertices;
2701 std::map<MVertex *, MVertex *> duplicates;
2702 for(std::size_t i = 0; i < entities.size(); i++) {
2703 GEntity *ge = entities[i];
2704 for(std::size_t j = 0; j < ge->mesh_vertices.size(); j++) {
2705 MVertex *v = ge->mesh_vertices[j];
2706 MVertex *v2 = pos.insert(v);
2707 if(v2)
2708 duplicates[v] = v2; // v should be removed
2709 else
2710 vertices[v->getIndex()] = v;
2711 }
2712 }
2713
2714 int num = (int)duplicates.size();
2715 Msg::Info("Found %d duplicate nodes ", num);
2716
2717 if(!num) {
2718 Msg::Info("No duplicate nodes found");
2719 return 0;
2720 }
2721
2722 for(std::size_t i = 0; i < entities.size(); i++) {
2723 GEntity *ge = entities[i];
2724 // clear list of vertices owned by entity
2725 ge->mesh_vertices.clear();
2726 // replace vertices in element
2727 for(std::size_t j = 0; j < ge->getNumMeshElements(); j++) {
2728 MElement *e = ge->getMeshElement(j);
2729 for(std::size_t k = 0; k < e->getNumVertices(); k++) {
2730 auto it = duplicates.find(e->getVertex(k));
2731 if(it != duplicates.end()) e->setVertex(k, it->second);
2732 }
2733 }
2734 // replace vertices in periodic copies
2735 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &corrVtcs = ge->correspondingVertices;
2736 if(corrVtcs.size()) {
2737 std::map<MVertex *, MVertex *>::iterator cIter;
2738 for(cIter = duplicates.begin(); cIter != duplicates.end(); ++cIter) {
2739 MVertex *oldTgt = cIter->first;
2740 MVertex *newTgt = cIter->second;
2741 auto cvIter = corrVtcs.find(oldTgt);
2742 if(cvIter != corrVtcs.end()) {
2743 MVertex *src = cvIter->second;
2744 corrVtcs.erase(cvIter);
2745 corrVtcs[newTgt] = src;
2746 }
2747 }
2748 for(cIter = corrVtcs.begin(); cIter != corrVtcs.end(); ++cIter) {
2749 MVertex *oldSrc = cIter->second;
2750 auto nIter = duplicates.find(oldSrc);
2751 if(nIter != duplicates.end()) {
2752 MVertex *tgt = cIter->first;
2753 MVertex *newSrc = nIter->second;
2754 corrVtcs[tgt] = newSrc;
2755 }
2756 }
2757 }
2758 }
2759
2760 destroyMeshCaches();
2761 _associateEntityWithMeshVertices();
2762 _storeVerticesInEntities(vertices);
2763
2764 // delete duplicates
2765 std::vector<MVertex *> to_delete;
2766 for(auto it = duplicates.begin(); it != duplicates.end(); it++)
2767 to_delete.push_back(it->first);
2768 for(std::size_t i = 0; i < to_delete.size(); i++) delete to_delete[i];
2769
2770 if(num)
2771 Msg::Info("Removed %d duplicate mesh node%s", num, num > 1 ? "s" : "");
2772
2773 Msg::StatusBar(true, "Done removing duplicate mesh nodes");
2774 return num;
2775 }
2776
alignPeriodicBoundaries()2777 void GModel::alignPeriodicBoundaries()
2778 {
2779 // Is this still necessary/useful?
2780 // 1) It's quite horrible
2781 // 2) It's only called when reading MSH2 files
2782
2783 Msg::Debug("Aligning periodic boundaries");
2784
2785 // realigning edges
2786
2787 for(auto it = firstEdge(); it != lastEdge(); ++it) {
2788 GEdge *tgt = *it;
2789 GEdge *src = dynamic_cast<GEdge *>(tgt->getMeshMaster());
2790
2791 if(src != nullptr && src != tgt) {
2792 // compose a search list on master edge
2793
2794 std::map<MEdge, MLine *, MEdgeLessThan> srcLines;
2795 for(std::size_t i = 0; i < src->getNumMeshElements(); i++) {
2796 MLine *srcLine = dynamic_cast<MLine *>(src->getMeshElement(i));
2797 if(!srcLine) {
2798 Msg::Debug("Master element %d is not a line",
2799 src->getMeshElement(i)->getNum());
2800 return;
2801 }
2802 srcLines[MEdge(srcLine->getVertex(0), srcLine->getVertex(1))] = srcLine;
2803 }
2804
2805 // run through slave edge elements
2806 // - check whether we find a counterpart (if not, abort)
2807 // - check orientation and reorient if necessary
2808
2809 for(std::size_t i = 0; i < tgt->getNumMeshElements(); ++i) {
2810 MLine *tgtLine = dynamic_cast<MLine *>(tgt->getMeshElement(i));
2811
2812 if(!tgtLine) {
2813 Msg::Debug("Slave element %d is not a line",
2814 tgt->getMeshElement(i)->getNum());
2815 return;
2816 }
2817
2818 MVertex *tgtVtcs[2];
2819 for(int iVtx = 0; iVtx < 2; iVtx++) {
2820 MVertex *tgtVtx = tgtLine->getVertex(iVtx);
2821 GEntity *ge = tgtVtx->onWhat();
2822 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &geV2v = ge->correspondingVertices;
2823 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &v2v = tgt->correspondingVertices;
2824 auto srcIter = v2v.find(tgtVtx);
2825 if(srcIter == v2v.end() || !srcIter->second) {
2826 srcIter = geV2v.find(tgtVtx);
2827 if(srcIter == geV2v.end() || !srcIter->second) {
2828 Msg::Debug(
2829 "Could not find periodic counterpart of node %d on curve %d "
2830 "or on entity %d of dimension %d",
2831 tgtVtx->getNum(), tgt->tag(), ge->tag(), ge->dim());
2832 return;
2833 }
2834 else
2835 tgtVtcs[iVtx] = srcIter->second;
2836 }
2837 else
2838 tgtVtcs[iVtx] = srcIter->second;
2839 }
2840
2841 MEdge tgtEdge(tgtVtcs[0], tgtVtcs[1]);
2842
2843 auto sIter = srcLines.find(tgtEdge);
2844
2845 if(sIter == srcLines.end() || !sIter->second) {
2846 Msg::Debug(
2847 "Could not find periodic counterpart of mesh edge %d-%d on "
2848 "curve %d for mesh edge %d-%d on curve %d",
2849 tgtLine->getVertex(0)->getNum(), tgtLine->getVertex(1)->getNum(),
2850 tgt->tag(), tgtVtcs[0]->getNum(), tgtVtcs[1]->getNum(), src->tag());
2851 return;
2852 }
2853 else {
2854 MLine *srcLine = sIter->second;
2855 MEdge srcEdge(srcLine->getVertex(0), srcLine->getVertex(1));
2856 if(tgtEdge.computeCorrespondence(srcEdge) == -1) tgtLine->reverse();
2857 }
2858 }
2859 }
2860 }
2861
2862 // run through all model faces
2863
2864 for(auto it = firstFace(); it != lastFace(); ++it) {
2865 GFace *tgt = *it;
2866 GFace *src = dynamic_cast<GFace *>(tgt->getMeshMaster());
2867 if(src != nullptr && src != tgt) {
2868 std::map<MFace, MElement *, MFaceLessThan> srcElmts;
2869
2870 for(std::size_t i = 0; i < src->getNumMeshElements(); ++i) {
2871 MElement *srcElmt = src->getMeshElement(i);
2872 int nbVtcs = 0;
2873 if(dynamic_cast<MTriangle *>(srcElmt)) nbVtcs = 3;
2874 if(dynamic_cast<MQuadrangle *>(srcElmt)) nbVtcs = 4;
2875 std::vector<MVertex *> vtcs;
2876 vtcs.reserve(nbVtcs);
2877 for(int iVtx = 0; iVtx < nbVtcs; iVtx++) {
2878 vtcs.push_back(srcElmt->getVertex(iVtx));
2879 }
2880 srcElmts[MFace(vtcs)] = srcElmt;
2881 }
2882
2883 for(std::size_t i = 0; i < tgt->getNumMeshElements(); ++i) {
2884 MElement *tgtElmt = tgt->getMeshElement(i);
2885 MTriangle *tgtTri = dynamic_cast<MTriangle *>(tgtElmt);
2886 MQuadrangle *tgtQua = dynamic_cast<MQuadrangle *>(tgtElmt);
2887
2888 int nbVtcs = 0;
2889 if(tgtTri) nbVtcs = 3;
2890 if(tgtQua) nbVtcs = 4;
2891
2892 std::vector<MVertex *> vtcs;
2893 for(int iVtx = 0; iVtx < nbVtcs; iVtx++) {
2894 MVertex *vtx = tgtElmt->getVertex(iVtx);
2895 GEntity *ge = vtx->onWhat();
2896
2897 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &geV2v = ge->correspondingVertices;
2898 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &v2v = tgt->correspondingVertices;
2899
2900 auto vIter = v2v.find(vtx);
2901 if(vIter == v2v.end() || !vIter->second) {
2902 vIter = geV2v.find(vtx);
2903 if(vIter == geV2v.end() || !vIter->second) {
2904 Msg::Debug("Could not find periodic counterpart of node %d in "
2905 "surface %d or in entity %d of dimension %d",
2906 vtx->getNum(), tgt->tag(), ge->tag(), ge->dim());
2907 return;
2908 }
2909 else
2910 vtcs.push_back(vIter->second);
2911 }
2912 else
2913 vtcs.push_back(vIter->second);
2914 }
2915 MFace tgtFace(vtcs);
2916
2917 auto mIter = srcElmts.find(tgtFace);
2918 if(mIter == srcElmts.end()) {
2919 std::ostringstream faceDef;
2920 for(int iVtx = 0; iVtx < nbVtcs; iVtx++) {
2921 faceDef << vtcs[iVtx]->getNum() << " ";
2922 }
2923 Msg::Debug("Could not find periodic counterpart of mesh face %s in "
2924 "surface %d connected to surface %d",
2925 faceDef.str().c_str(), tgt->tag(), src->tag());
2926 return;
2927 }
2928 else {
2929 const MFace &srcFace = mIter->first;
2930 MElement *srcElmt = mIter->second;
2931 std::vector<MVertex *> srcVtcs;
2932
2933 if((tgtTri && !dynamic_cast<MTriangle *>(srcElmt)) ||
2934 (tgtQua && !dynamic_cast<MQuadrangle *>(srcElmt))) {
2935 Msg::Error("Invalid source/target elements");
2936 return;
2937 }
2938
2939 int rotation = 0;
2940 bool swap = false;
2941
2942 if(!tgtFace.computeCorrespondence(srcFace, rotation, swap)) {
2943 Msg::Debug(
2944 "Could not find correspondance between mesh face %d-%d-%d "
2945 "(slave) "
2946 "and %d-%d-%d (master)",
2947 tgtElmt->getVertex(0)->getNum(), tgtElmt->getVertex(1)->getNum(),
2948 tgtElmt->getVertex(2)->getNum(), srcElmt->getVertex(0)->getNum(),
2949 srcElmt->getVertex(1)->getNum(), srcElmt->getVertex(2)->getNum());
2950 return;
2951 }
2952
2953 if(tgtTri) tgtTri->reorient(rotation, swap);
2954 if(tgtQua) tgtQua->reorient(rotation, swap);
2955 }
2956 }
2957 }
2958 }
2959 Msg::Debug("Done aligning periodic boundaries");
2960 }
2961
connectMElementsByMFace(const MFace & f,std::multimap<MFace,MElement *,MFaceLessThan> & e2f,std::set<MElement * > & group,std::set<MFace,MFaceLessThan> & touched,int recur_level)2962 static void connectMElementsByMFace(
2963 const MFace &f, std::multimap<MFace, MElement *, MFaceLessThan> &e2f,
2964 std::set<MElement *> &group, std::set<MFace, MFaceLessThan> &touched,
2965 int recur_level)
2966 {
2967 // this is very slow...
2968 std::stack<MFace> _stack;
2969 _stack.push(f);
2970
2971 while(!_stack.empty()) {
2972 MFace ff = _stack.top();
2973 _stack.pop();
2974 if(touched.find(ff) == touched.end()) {
2975 touched.insert(ff);
2976 for(auto it = e2f.lower_bound(ff); it != e2f.upper_bound(ff); ++it) {
2977 group.insert(it->second);
2978 for(int i = 0; i < it->second->getNumFaces(); ++i) {
2979 _stack.push(it->second->getFace(i));
2980 }
2981 }
2982 }
2983 }
2984 }
2985
connectedVolumes(std::vector<MElement * > & elements,std::vector<std::vector<MElement * >> & regs)2986 static int connectedVolumes(std::vector<MElement *> &elements,
2987 std::vector<std::vector<MElement *> > ®s)
2988 {
2989 std::multimap<MFace, MElement *, MFaceLessThan> e2f;
2990 for(std::size_t i = 0; i < elements.size(); ++i) {
2991 for(int j = 0; j < elements[i]->getNumFaces(); j++) {
2992 e2f.insert(std::make_pair(elements[i]->getFace(j), elements[i]));
2993 }
2994 }
2995 while(!e2f.empty()) {
2996 std::set<MElement *> group;
2997 std::set<MFace, MFaceLessThan> touched;
2998 connectMElementsByMFace(e2f.begin()->first, e2f, group, touched, 0);
2999 std::vector<MElement *> temp;
3000 temp.insert(temp.begin(), group.begin(), group.end());
3001 regs.push_back(temp);
3002 for(auto it = touched.begin(); it != touched.end(); ++it) e2f.erase(*it);
3003 }
3004 return regs.size();
3005 }
3006
connectMElementsByMEdge(const MEdge & e,std::multimap<MEdge,MElement *,MEdgeLessThan> & e2e,std::set<MElement * > & group,std::set<MEdge,MEdgeLessThan> & touched)3007 static void connectMElementsByMEdge(
3008 const MEdge &e, std::multimap<MEdge, MElement *, MEdgeLessThan> &e2e,
3009 std::set<MElement *> &group, std::set<MEdge, MEdgeLessThan> &touched)
3010 {
3011 // this is very slow...
3012 std::stack<MEdge> _stack;
3013 _stack.push(e);
3014
3015 while(!_stack.empty()) {
3016 MEdge ee = _stack.top();
3017 _stack.pop();
3018 if(touched.find(ee) == touched.end()) {
3019 touched.insert(ee);
3020 for(auto it = e2e.lower_bound(ee); it != e2e.upper_bound(ee); ++it) {
3021 group.insert(it->second);
3022 for(int i = 0; i < it->second->getNumEdges(); ++i) {
3023 _stack.push(it->second->getEdge(i));
3024 }
3025 }
3026 }
3027 }
3028 }
3029
connectedSurfaces(std::vector<MElement * > & elements,std::vector<std::vector<MElement * >> & faces)3030 static int connectedSurfaces(std::vector<MElement *> &elements,
3031 std::vector<std::vector<MElement *> > &faces)
3032 {
3033 std::multimap<MEdge, MElement *, MEdgeLessThan> e2e;
3034 for(std::size_t i = 0; i < elements.size(); ++i) {
3035 for(int j = 0; j < elements[i]->getNumEdges(); j++) {
3036 e2e.insert(std::make_pair(elements[i]->getEdge(j), elements[i]));
3037 }
3038 }
3039 while(!e2e.empty()) {
3040 std::set<MElement *> group;
3041 std::set<MEdge, MEdgeLessThan> touched;
3042 connectMElementsByMEdge(e2e.begin()->first, e2e, group, touched);
3043 std::vector<MElement *> temp;
3044 temp.insert(temp.begin(), group.begin(), group.end());
3045 faces.push_back(temp);
3046 for(auto it = touched.begin(); it != touched.end(); ++it) e2e.erase(*it);
3047 }
3048 return faces.size();
3049 }
3050
makeDiscreteRegionsSimplyConnected()3051 void GModel::makeDiscreteRegionsSimplyConnected()
3052 {
3053 Msg::Info("Making discrete regions simply connected...");
3054
3055 std::vector<discreteRegion *> discRegions;
3056 for(auto it = firstRegion(); it != lastRegion(); it++)
3057 if((*it)->geomType() == GEntity::DiscreteVolume)
3058 discRegions.push_back((discreteRegion *)*it);
3059
3060 std::set<MVertex *> touched;
3061
3062 for(auto itR = discRegions.begin(); itR != discRegions.end(); itR++) {
3063 std::vector<MElement *> allElements((*itR)->getNumMeshElements());
3064 for(std::size_t i = 0; i < (*itR)->getNumMeshElements(); i++)
3065 allElements[i] = (*itR)->getMeshElement(i);
3066
3067 std::vector<std::vector<MElement *> > conRegions;
3068 int nbRegions = connectedVolumes(allElements, conRegions);
3069 if(nbRegions > 1) remove(*itR);
3070
3071 for(int ire = 0; ire < nbRegions; ire++) {
3072 int numR =
3073 (nbRegions == 1) ? (*itR)->tag() : getMaxElementaryNumber(3) + 1;
3074 discreteRegion *r = new discreteRegion(this, numR);
3075 add(r);
3076 std::vector<MElement *> myElements = conRegions[ire];
3077 std::set<MVertex *> myVertices;
3078 for(std::size_t i = 0; i < myElements.size(); i++) {
3079 MElement *e = myElements[i];
3080 std::vector<MVertex *> verts;
3081 e->getVertices(verts);
3082 for(std::size_t k = 0; k < verts.size(); k++) {
3083 if(verts[k]->onWhat() && verts[k]->onWhat()->dim() == 3) {
3084 if(touched.find(verts[k]) == touched.end()) {
3085 verts[k]->setEntity(r);
3086 myVertices.insert(verts[k]);
3087 touched.insert(verts[k]);
3088 }
3089 }
3090 }
3091 MElementFactory factory;
3092 MElement *e2 = factory.create(e->getTypeForMSH(), verts, e->getNum(),
3093 e->getPartition());
3094 switch(e2->getType()) {
3095 case TYPE_TET: r->tetrahedra.push_back((MTetrahedron *)e2); break;
3096 case TYPE_HEX: r->hexahedra.push_back((MHexahedron *)e2); break;
3097 case TYPE_PRI: r->prisms.push_back((MPrism *)e2); break;
3098 case TYPE_PYR: r->pyramids.push_back((MPyramid *)e2); break;
3099 case TYPE_TRIH: r->trihedra.push_back((MTrihedron *)e2); break;
3100 }
3101 }
3102 r->mesh_vertices.insert(r->mesh_vertices.begin(), myVertices.begin(),
3103 myVertices.end());
3104 }
3105 }
3106
3107 Msg::Info("Done making discrete regions simply connected");
3108 }
3109
makeDiscreteFacesSimplyConnected()3110 void GModel::makeDiscreteFacesSimplyConnected()
3111 {
3112 Msg::Info("Making discrete faces simply connected...");
3113
3114 std::vector<discreteFace *> discFaces;
3115 for(auto it = firstFace(); it != lastFace(); it++)
3116 if((*it)->geomType() == GEntity::DiscreteSurface)
3117 discFaces.push_back((discreteFace *)*it);
3118
3119 std::set<MVertex *> touched;
3120
3121 for(auto itF = discFaces.begin(); itF != discFaces.end(); itF++) {
3122 std::vector<MElement *> allElements((*itF)->getNumMeshElements());
3123 for(std::size_t i = 0; i < (*itF)->getNumMeshElements(); i++)
3124 allElements[i] = (*itF)->getMeshElement(i);
3125
3126 std::vector<std::vector<MElement *> > conFaces;
3127 int nbFaces = connectedSurfaces(allElements, conFaces);
3128 if(nbFaces > 1) remove(*itF);
3129
3130 for(int ifa = 0; ifa < nbFaces; ifa++) {
3131 int numF = (nbFaces == 1) ? (*itF)->tag() : getMaxElementaryNumber(2) + 1;
3132 discreteFace *f = new discreteFace(this, numF);
3133 add(f);
3134 std::vector<MElement *> myElements = conFaces[ifa];
3135 std::set<MVertex *> myVertices;
3136 for(std::size_t i = 0; i < myElements.size(); i++) {
3137 MElement *e = myElements[i];
3138 std::vector<MVertex *> verts;
3139 e->getVertices(verts);
3140 for(std::size_t k = 0; k < verts.size(); k++) {
3141 if(verts[k]->onWhat() && verts[k]->onWhat()->dim() == 2) {
3142 if(touched.find(verts[k]) == touched.end()) {
3143 verts[k]->setEntity(f);
3144 myVertices.insert(verts[k]);
3145 touched.insert(verts[k]);
3146 }
3147 }
3148 }
3149 MElementFactory factory;
3150 MElement *e2 = factory.create(e->getTypeForMSH(), verts, e->getNum(),
3151 e->getPartition());
3152 if(e2->getType() == TYPE_TRI)
3153 f->triangles.push_back((MTriangle *)e2);
3154 else
3155 f->quadrangles.push_back((MQuadrangle *)e2);
3156 }
3157 f->mesh_vertices.insert(f->mesh_vertices.begin(), myVertices.begin(),
3158 myVertices.end());
3159 }
3160 }
3161
3162 Msg::Info("Done making discrete faces simply connected");
3163 }
3164
3165 static void
makeSimplyConnected(std::map<int,std::vector<MElement * >> elements[11])3166 makeSimplyConnected(std::map<int, std::vector<MElement *> > elements[11])
3167 {
3168 // only for tetras and triangles
3169 Msg::Info("Make simply connected regions and surfaces");
3170 std::vector<int> regs;
3171 for(auto it = elements[4].begin(); it != elements[4].end(); it++)
3172 regs.push_back(it->first);
3173 std::multimap<MFace, MElement *, MFaceLessThan> f2e;
3174 if(regs.size() > 2) {
3175 for(std::size_t i = 0; i < regs.size(); i++) {
3176 for(std::size_t j = 0; j < elements[4][regs[i]].size(); j++) {
3177 MElement *el = elements[4][regs[i]][j];
3178 for(int k = 0; k < el->getNumFaces(); k++)
3179 f2e.insert(std::make_pair(el->getFace(k), el));
3180 }
3181 }
3182 }
3183 for(std::size_t i = 0; i < regs.size(); i++) {
3184 int ri = regs[i];
3185 std::vector<MElement *> allElements;
3186 for(std::size_t j = 0; j < elements[4][ri].size(); j++)
3187 allElements.push_back(elements[4][ri][j]);
3188 std::vector<std::vector<MElement *> > conRegions;
3189 int nbConRegions = connectedVolumes(allElements, conRegions);
3190 Msg::Info("%d connected regions (reg=%d)", nbConRegions, ri);
3191 std::size_t maxNumEl = 1;
3192 for(int j = 0; j < nbConRegions; j++)
3193 if(conRegions[j].size() > maxNumEl) maxNumEl = conRegions[j].size();
3194 for(int j = 0; j < nbConRegions; j++) {
3195 // remove conRegions containing few elements
3196 if(conRegions[j].size() < maxNumEl * 1.e-4) {
3197 // find adjacent region
3198 int r2 = ri;
3199 if(regs.size() == 2)
3200 r2 = (ri + 1) % 2;
3201 else {
3202 for(std::size_t k = 0; k < conRegions[j].size(); k++) {
3203 MElement *el = conRegions[j][k];
3204 for(int l = 0; l < el->getNumFaces(); l++) {
3205 MFace mf = el->getFace(l);
3206 auto itl = f2e.lower_bound(mf);
3207 for(; itl != f2e.upper_bound(mf); itl++) {
3208 if(itl->second != el) break;
3209 }
3210 MElement *el2 = itl->second;
3211 bool sameRegion = false;
3212 for(std::size_t m = 0; m < conRegions[j].size(); m++)
3213 if(conRegions[j][m] == el2) {
3214 sameRegion = true;
3215 break;
3216 }
3217 if(sameRegion) continue;
3218 for(std::size_t m = 0; m < regs.size(); m++) {
3219 int rm = regs[m];
3220 if(rm == ri) continue;
3221 for(std::size_t n = 0; n < elements[4][rm].size(); n++)
3222 if(elements[4][rm][n] == el2) {
3223 r2 = rm;
3224 break;
3225 }
3226 if(r2 != ri) break;
3227 }
3228 if(r2 != ri) break;
3229 }
3230 if(r2 != ri) break;
3231 }
3232 if(r2 == ri)
3233 Msg::Warning("Element not found for simply connected regions");
3234 }
3235
3236 for(std::size_t k = 0; k < conRegions[j].size(); k++) {
3237 MElement *el = conRegions[j][k];
3238 std::size_t l = 0;
3239 for(; l < elements[4][ri].size(); l++)
3240 if(elements[4][ri][l] == el) break;
3241 elements[4][ri].erase(elements[4][ri].begin() + l);
3242 elements[4][r2].push_back(el);
3243 }
3244 }
3245 }
3246 }
3247
3248 std::vector<int> faces;
3249 for(auto it = elements[2].begin(); it != elements[2].end(); it++)
3250 faces.push_back(it->first);
3251 std::multimap<MEdge, MElement *, MEdgeLessThan> e2e;
3252 if(faces.size() > 2) {
3253 for(std::size_t i = 0; i < faces.size(); i++) {
3254 for(std::size_t j = 0; j < elements[2][faces[i]].size(); j++) {
3255 MElement *el = elements[2][faces[i]][j];
3256 for(int k = 0; k < el->getNumEdges(); k++)
3257 e2e.insert(std::make_pair(el->getEdge(k), el));
3258 }
3259 }
3260 }
3261 for(std::size_t i = 0; i < faces.size(); i++) {
3262 int fi = faces[i];
3263 std::vector<MElement *> allElements;
3264 for(std::size_t j = 0; j < elements[2][fi].size(); j++)
3265 allElements.push_back(elements[2][fi][j]);
3266 std::vector<std::vector<MElement *> > conSurfaces;
3267 int nbConSurfaces = connectedSurfaces(allElements, conSurfaces);
3268 Msg::Info("%d connected surfaces (reg=%d)", nbConSurfaces, fi);
3269 std::size_t maxNumEl = 1;
3270 for(int j = 0; j < nbConSurfaces; j++)
3271 if(conSurfaces[j].size() > maxNumEl) maxNumEl = conSurfaces[j].size();
3272 for(int j = 0; j < nbConSurfaces; j++) {
3273 // remove conSurfaces containing few elements
3274 if(conSurfaces[j].size() < maxNumEl * 1.e-4) {
3275 // find adjacent surface
3276 int f2 = fi;
3277 if(faces.size() == 2)
3278 f2 = (fi + 1) % 2;
3279 else {
3280 for(std::size_t k = 0; k < conSurfaces[j].size(); k++) {
3281 MElement *el = conSurfaces[j][k];
3282 for(int l = 0; l < el->getNumEdges(); l++) {
3283 MEdge me = el->getEdge(l);
3284 auto itl = e2e.lower_bound(me);
3285 for(; itl != e2e.upper_bound(me); itl++) {
3286 if(itl->second != el) break;
3287 }
3288 MElement *el2 = itl->second;
3289 bool sameSurface = false;
3290 for(std::size_t m = 0; m < conSurfaces[j].size(); m++)
3291 if(conSurfaces[j][m] == el2) {
3292 sameSurface = true;
3293 break;
3294 }
3295 if(sameSurface) continue;
3296 for(std::size_t m = 0; m < faces.size(); m++) {
3297 int fm = faces[m];
3298 if(fm == fi) continue;
3299 for(std::size_t n = 0; n < elements[2][fm].size(); n++)
3300 if(elements[2][fm][n] == el2) {
3301 f2 = fm;
3302 break;
3303 }
3304 if(f2 != fi) break;
3305 }
3306 if(f2 != fi) break;
3307 }
3308 if(f2 != fi) break;
3309 }
3310 if(f2 == fi)
3311 Msg::Warning("Element not found for simply connected surfaces");
3312 }
3313 for(std::size_t k = 0; k < conSurfaces[j].size(); k++) {
3314 MElement *el = conSurfaces[j][k];
3315 std::size_t l = 0;
3316 for(; l < elements[2][fi].size(); l++)
3317 if(elements[2][fi][l] == el) break;
3318 elements[2][fi].erase(elements[2][fi].begin() + l);
3319 elements[2][f2].push_back(el);
3320 }
3321 }
3322 }
3323 }
3324 }
3325
buildCutGModel(gLevelset * ls,bool cutElem,bool saveTri)3326 GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
3327 {
3328 if(saveTri)
3329 CTX::instance()->mesh.saveTri = 1;
3330 else
3331 CTX::instance()->mesh.saveTri = 0;
3332
3333 std::map<int, std::vector<MElement *> > elements[11];
3334 std::map<int, std::map<int, std::string> > physicals[4];
3335 std::map<int, MVertex *> vertexMap;
3336
3337 if(cutElem)
3338 Msg::Info("Cutting mesh...");
3339 else
3340 Msg::Info("Splitting mesh...");
3341 double t1 = Cpu(), w1 = TimeOfDay();
3342
3343 GModel *cutGM =
3344 buildCutMesh(this, ls, elements, vertexMap, physicals, cutElem);
3345
3346 if(!cutElem) makeSimplyConnected(elements);
3347
3348 for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
3349 cutGM->_storeElementsInEntities(elements[i]);
3350 cutGM->_associateEntityWithMeshVertices();
3351 cutGM->_storeVerticesInEntities(vertexMap);
3352
3353 for(int i = 0; i < 4; i++) {
3354 cutGM->_storePhysicalTagsInEntities(i, physicals[i]);
3355 auto it = physicals[i].begin();
3356 for(; it != physicals[i].end(); it++) {
3357 auto it2 = it->second.begin();
3358 for(; it2 != it->second.end(); it2++)
3359 if(it2->second != "")
3360 cutGM->setPhysicalName(it2->second, i, it2->first);
3361 }
3362 }
3363
3364 if(cutElem)
3365 Msg::Info("Mesh cutting completed (Wall %gs, CPU %gs)", TimeOfDay() - w1,
3366 Cpu() - t1);
3367 else
3368 Msg::Info("Mesh splitting completed (Wall %gs, CPU %gs)", TimeOfDay() - w1,
3369 Cpu() - t1);
3370
3371 return cutGM;
3372 }
3373
load(const std::string & fileName)3374 void GModel::load(const std::string &fileName)
3375 {
3376 GModel *temp = GModel::current();
3377 GModel::setCurrent(this);
3378 MergeFile(fileName, true);
3379 GModel::setCurrent(temp);
3380 }
3381
save(const std::string & fileName)3382 void GModel::save(const std::string &fileName)
3383 {
3384 GModel *temp = GModel::current();
3385 GModel::setCurrent(this);
3386 CreateOutputFile(fileName, FORMAT_AUTO);
3387 GModel::setCurrent(temp);
3388 }
3389
readGEO(const std::string & name)3390 int GModel::readGEO(const std::string &name)
3391 {
3392 // readGEO is static, because it can create several models
3393 ParseFile(name, true);
3394 // sync OCC first, as GEO_Internals currently contains attributes (physicals)
3395 // that should also be applied to entities from OCC_Internals
3396 if(GModel::current()->getOCCInternals())
3397 GModel::current()->getOCCInternals()->synchronize(GModel::current());
3398 GModel::current()->getGEOInternals()->synchronize(GModel::current());
3399 return 1;
3400 }
3401
setPhysicalNumToEntitiesInBox(int EntityDimension,int PhysicalNumber,const SBoundingBox3d & box)3402 void GModel::setPhysicalNumToEntitiesInBox(int EntityDimension,
3403 int PhysicalNumber,
3404 const SBoundingBox3d &box)
3405 {
3406 std::vector<GEntity *> entities;
3407 getEntitiesInBox(entities, box, EntityDimension);
3408 for(std::size_t i = 0; i < entities.size(); i++)
3409 entities[i]->addPhysicalEntity(PhysicalNumber);
3410 }
3411
setPhysicalNumToEntitiesInBox(int EntityDimension,int PhysicalNumber,std::vector<double> p1,std::vector<double> p2)3412 void GModel::setPhysicalNumToEntitiesInBox(int EntityDimension,
3413 int PhysicalNumber,
3414 std::vector<double> p1,
3415 std::vector<double> p2)
3416 {
3417 if(p1.size() != 3 || p2.size() != 3) return;
3418 SBoundingBox3d box(p1[0], p1[2], p1[2], p2[0], p2[1], p2[3]);
3419 setPhysicalNumToEntitiesInBox(EntityDimension, PhysicalNumber, box);
3420 }
3421
classifySurfaces(double angleThreshold,bool includeBoundary,bool forReparametrization,double curveAngleThreshold)3422 void GModel::classifySurfaces(double angleThreshold, bool includeBoundary,
3423 bool forReparametrization,
3424 double curveAngleThreshold)
3425 {
3426 classifyFaces(this, angleThreshold, includeBoundary, forReparametrization,
3427 curveAngleThreshold);
3428 }
3429
addHomologyRequest(const std::string & type,const std::vector<int> & domain,const std::vector<int> & subdomain,const std::vector<int> & dim)3430 void GModel::addHomologyRequest(const std::string &type,
3431 const std::vector<int> &domain,
3432 const std::vector<int> &subdomain,
3433 const std::vector<int> &dim)
3434 {
3435 typedef std::pair<const std::vector<int>, const std::vector<int> > dpair;
3436 typedef std::pair<const std::string, const std::vector<int> > tpair;
3437 dpair p(domain, subdomain);
3438 tpair p2(type, dim);
3439 _homologyRequests.insert(std::make_pair(p, p2));
3440 }
3441
computeHomology()3442 void GModel::computeHomology()
3443 {
3444 if(_homologyRequests.empty()) return;
3445
3446 #if defined(HAVE_KBIPACK)
3447 Msg::StatusBar(true, "Homology and cohomology computation...");
3448 double t1 = Cpu(), w1 = TimeOfDay();
3449
3450 // find unique domain/subdomain requests
3451 typedef std::pair<const std::vector<int>, const std::vector<int> > dpair;
3452 typedef std::pair<const std::string, const std::vector<int> > tpair;
3453 std::set<dpair> domains;
3454 for(auto it = _homologyRequests.begin(); it != _homologyRequests.end(); it++)
3455 domains.insert(it->first);
3456 Msg::Info("Number of cell complexes to construct: %d", domains.size());
3457
3458 for(auto it = domains.begin(); it != domains.end(); it++) {
3459 std::pair<std::multimap<dpair, tpair>::iterator,
3460 std::multimap<dpair, tpair>::iterator>
3461 itp = _homologyRequests.equal_range(*it);
3462 bool prepareToRestore = (itp.first != --itp.second);
3463 itp.second++;
3464 std::vector<int> imdomain;
3465 Homology *homology =
3466 new Homology(this, itp.first->first.first, itp.first->first.second,
3467 imdomain, prepareToRestore);
3468
3469 for(auto itt = itp.first; itt != itp.second; itt++) {
3470 std::string type = itt->second.first;
3471 std::vector<int> dim0 = itt->second.second;
3472 if(dim0.empty())
3473 for(int i = 0; i < getDim(); i++) dim0.push_back(i);
3474 std::vector<int> dim;
3475
3476 std::stringstream ss;
3477 for(std::size_t i = 0; i < dim0.size(); i++) {
3478 int d = dim0.at(i);
3479 if(d >= 0 && d <= getDim()) {
3480 dim.push_back(d);
3481 ss << "H";
3482 if(type == "Homology") ss << "_";
3483 if(type == "Cohomology") ss << "^";
3484 ss << d;
3485 if(i < dim0.size() - 1 && dim0.at(i + 1) >= 0 &&
3486 dim0.at(i + 1) <= getDim())
3487 ss << ", ";
3488 }
3489 }
3490 std::string dims = ss.str();
3491
3492 if(type != "Homology" && type != "Cohomology" && type != "Betti") {
3493 Msg::Error("Unknown type of homology computation: %s", type.c_str());
3494 }
3495 else if(dim.empty()) {
3496 Msg::Error("Invalid homology computation dimensions given");
3497 }
3498 else if(type == "Betti") {
3499 homology->findBettiNumbers();
3500 }
3501 else if(type == "Homology" && !homology->isHomologyComputed(dim)) {
3502 homology->findHomologyBasis(dim);
3503 Msg::Info("Homology space basis chains to save: %s", dims.c_str());
3504 for(std::size_t i = 0; i < dim.size(); i++) {
3505 homology->addChainsToModel(dim.at(i));
3506 }
3507 }
3508 else if(type == "Cohomology" && !homology->isCohomologyComputed(dim)) {
3509 homology->findCohomologyBasis(dim);
3510 Msg::Info("Cohomology space basis cochains to save: %s", dims.c_str());
3511 for(std::size_t i = 0; i < dim.size(); i++) {
3512 homology->addCochainsToModel(dim.at(i));
3513 }
3514 }
3515 }
3516 pruneMeshVertexAssociations();
3517 delete homology;
3518 }
3519 Msg::Info("");
3520
3521 double t2 = Cpu(), w2 = TimeOfDay();
3522 Msg::StatusBar(true,
3523 "Done homology and cohomology computation "
3524 "(Wall %gs, CPU %gs)",
3525 w2 - w1, t2 - t1);
3526
3527 #else
3528 Msg::Error("Homology computation requires KBIPACK");
3529 #endif
3530 }
3531
computeSizeField()3532 void GModel::computeSizeField()
3533 {
3534 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
3535 FieldManager *fields = getFields();
3536 int myId = fields->newId();
3537 fields->newField(myId, std::string("AutomaticMeshSizeField"));
3538 fields->get(myId)->update();
3539 #else
3540 Msg::Error("Size field computation requires both HXT and P4EST");
3541 #endif
3542 }
3543
3544 #if !defined(HAVE_ACIS)
3545
createACISInternals()3546 void GModel::createACISInternals() {}
3547
deleteACISInternals()3548 void GModel::deleteACISInternals() {}
3549
readACISSAT(const std::string & fn)3550 int GModel::readACISSAT(const std::string &fn)
3551 {
3552 Msg::Error("Gmsh must be compiled with ACIS support to load '%s'",
3553 fn.c_str());
3554 return 0;
3555 }
3556
3557 #endif
3558
3559 #if !defined(HAVE_PARASOLID)
3560
createParasolidInternals()3561 void GModel::createParasolidInternals() {}
3562
deleteParasolidInternals()3563 void GModel::deleteParasolidInternals() {}
3564
readParasolidXMT(const std::string & fn)3565 int GModel::readParasolidXMT(const std::string &fn)
3566 {
3567 Msg::Error("Gmsh must be compiled with Parasolid support to read '%s'",
3568 fn.c_str());
3569 return 0;
3570 }
3571
writeParasolidXMT(const std::string & fn)3572 int GModel::writeParasolidXMT(const std::string &fn)
3573 {
3574 Msg::Error("Gmsh must be compiled with Parasolid support to write '%s'",
3575 fn.c_str());
3576 return 0;
3577 }
3578
readParasolidSTEP(const std::string & fn)3579 int GModel::readParasolidSTEP(const std::string &fn)
3580 {
3581 Msg::Error("Gmsh must be compiled with Parasolid support to read '%s'",
3582 fn.c_str());
3583 return 0;
3584 }
3585
writeParasolidSTEP(const std::string & fn)3586 int GModel::writeParasolidSTEP(const std::string &fn)
3587 {
3588 Msg::Error("Gmsh must be compiled with Parasolid support to write '%s'",
3589 fn.c_str());
3590 return 0;
3591 }
3592
3593 #endif
3594