1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Anthony Royer.
7
8 #include <vector>
9 #include <set>
10 #include <sstream>
11 #include <algorithm>
12 #include <ctime>
13 #include <limits>
14 #include <stack>
15 #include <cstdlib>
16 #include <map>
17 #include <unordered_map>
18 #include "GmshConfig.h"
19 #include "GmshMessage.h"
20 #include "GModel.h"
21 #include "ElementType.h"
22
23 struct OriGEntityPtrFullLessThan {
operator ()OriGEntityPtrFullLessThan24 bool operator()(const std::pair<int, GEntity *> &p1,
25 const std::pair<int, GEntity *> &p2) const
26 {
27 if(p1.first != p2.first) return p1.first < p2.first;
28 if(p1.second->dim() != p2.second->dim())
29 return p1.second->dim() < p2.second->dim();
30 return p1.second->tag() < p2.second->tag();
31 }
32 };
33
34 typedef std::set<std::pair<int, GEntity *>, OriGEntityPtrFullLessThan>
35 setorientity;
36
37 #define hashmap std::unordered_map
38 #define hashmapentity \
39 std::unordered_map<GEntity *, setorientity, GEntityPtrFullHash, \
40 GEntityPtrFullEqual>
41 #define hashmapelement \
42 std::unordered_map<MElement *, GEntity *, MElementPtrHash, MElementPtrEqual>
43 #define hashmapelementpart \
44 std::unordered_map<MElement *, int, MElementPtrHash, MElementPtrEqual>
45 #define hashmapface \
46 std::unordered_map<MFace, \
47 std::vector<std::pair<MElement *, std::vector<int> > >, \
48 MFaceHash, MFaceEqual>
49 #define hashmapedge \
50 std::unordered_map<MEdge, \
51 std::vector<std::pair<MElement *, std::vector<int> > >, \
52 MEdgeHash, MEdgeEqual>
53 #define hashmapvertex \
54 std::unordered_map<MVertex *, \
55 std::vector<std::pair<MElement *, std::vector<int> > >, \
56 MVertexPtrHash, MVertexPtrEqual>
57
58 #if defined(HAVE_METIS)
59
60 #include "OS.h"
61 #include "Context.h"
62 #include "partitionRegion.h"
63 #include "partitionFace.h"
64 #include "partitionEdge.h"
65 #include "partitionVertex.h"
66 #include "ghostRegion.h"
67 #include "ghostFace.h"
68 #include "ghostEdge.h"
69 #include "MFaceHash.h"
70 #include "MEdgeHash.h"
71 #include "MTriangle.h"
72 #include "MQuadrangle.h"
73 #include "MTetrahedron.h"
74 #include "MHexahedron.h"
75 #include "MPrism.h"
76 #include "MPyramid.h"
77 #include "MTrihedron.h"
78 #include "MElementCut.h"
79 #include "MPoint.h"
80
81 extern "C" {
82 #include <metis.h>
83 }
84
85 // Graph of the mesh for partitioning purposes.
86 class Graph {
87 private:
88 // The GModel
89 GModel *_model;
90 // The number of partitions
91 std::size_t _nparts;
92 // The number of elements
93 std::size_t _ne;
94 // The number of nodes
95 std::size_t _nn;
96 // The dimension of the mesh
97 int _dim;
98 // The list of nodes belonging to the ith element of the mesh are stored in
99 // consecutive locations of eind starting at position eptr[i] up to (but not
100 // including) position eptr[i+1]. The size of the eind array is of size equal
101 // to the sum of the number of nodes in all the elements of the mesh.
102 std::vector<idx_t> _eind;
103 // The size of the eptr array is n + 1, where n is the number of elements in
104 // the mesh.
105 std::vector<idx_t> _eptr;
106 // The metis graph structure
107 idx_t *_xadj, *_adjncy;
108 // The weight associated to each elements of eptr
109 idx_t *_vwgt;
110 // Element corresponding to each graph element in eptr
111 std::vector<MElement *> _element;
112 // Vertex corresponding to each graph vertex in eptr
113 std::vector<idx_t> _vertex;
114 // The partitions output from the partitioner, in an integer type independent
115 // from METIS
116 std::vector<int> _partition;
117
118 public:
Graph(GModel * model)119 Graph(GModel *model)
120 : _model(model), _nparts(0), _ne(0), _nn(0), _dim(0), _xadj(nullptr),
121 _adjncy(nullptr), _vwgt(nullptr)
122 {
123 }
~Graph()124 ~Graph() { clear(); }
nparts() const125 std::size_t nparts() const { return _nparts; };
ne() const126 std::size_t ne() const { return _ne; };
nn() const127 std::size_t nn() const { return _nn; };
dim() const128 int dim() const { return _dim; };
eind(std::size_t i) const129 std::size_t eind(std::size_t i) const { return _eind[i]; };
eptr(std::size_t i) const130 std::size_t eptr(std::size_t i) const { return _eptr[i]; };
xadj(std::size_t i) const131 idx_t xadj(std::size_t i) const { return _xadj[i]; };
xadj() const132 idx_t *xadj() const { return _xadj; };
adjncy(std::size_t i) const133 idx_t adjncy(std::size_t i) const { return _adjncy[i]; };
adjncy() const134 idx_t *adjncy() const { return _adjncy; };
vwgt() const135 idx_t *vwgt() const { return _vwgt; };
element(std::size_t i) const136 MElement *element(std::size_t i) const { return _element[i]; };
vertex(std::size_t i) const137 idx_t vertex(std::size_t i) const { return _vertex[i]; };
partition(std::size_t i) const138 int partition(std::size_t i) const { return _partition[i]; };
numNodes() const139 std::size_t numNodes() const { return _ne; };
numEdges() const140 std::size_t numEdges() const { return _xadj[_ne] / 2; };
nparts(std::size_t nparts)141 void nparts(std::size_t nparts) { _nparts = nparts; };
ne(std::size_t ne)142 void ne(std::size_t ne) { _ne = ne; };
nn(std::size_t nn)143 void nn(std::size_t nn) { _nn = nn; };
dim(int dim)144 void dim(int dim) { _dim = dim; };
eindResize(std::size_t size)145 void eindResize(std::size_t size)
146 {
147 _eind.clear();
148 _eind.resize(size, 0);
149 }
eind(std::size_t i,idx_t eind)150 void eind(std::size_t i, idx_t eind) { _eind[i] = eind; };
eptrResize(std::size_t size)151 void eptrResize(std::size_t size)
152 {
153 _eptr.clear();
154 _eptr.resize(size, 0);
155 }
eptr(std::size_t i,idx_t eptr)156 void eptr(std::size_t i, idx_t eptr) { _eptr[i] = eptr; };
elementResize(std::size_t size)157 void elementResize(std::size_t size)
158 {
159 _element.clear();
160 _element.resize(size, nullptr);
161 }
element(std::size_t i,MElement * element)162 void element(std::size_t i, MElement *element) { _element[i] = element; };
vertexResize(std::size_t size)163 void vertexResize(std::size_t size)
164 {
165 _vertex.clear();
166 _vertex.resize(size, -1);
167 }
adjncy(std::size_t i,idx_t adjncy)168 void adjncy(std::size_t i, idx_t adjncy) { _adjncy[i] = adjncy; };
vertex(std::size_t i,idx_t vertex)169 void vertex(std::size_t i, idx_t vertex) { _vertex[i] = vertex; };
partition(const std::vector<idx_t> & epart)170 void partition(const std::vector<idx_t> &epart)
171 {
172 // converts into METIS-independent integer type
173 _partition.resize(epart.size());
174 for(std::size_t i = 0; i < epart.size(); i++) _partition[i] = epart[i];
175 };
clear()176 void clear()
177 {
178 if(_xadj) {
179 delete[] _xadj;
180 _xadj = nullptr;
181 }
182 if(_adjncy) {
183 delete[] _adjncy;
184 _adjncy = nullptr;
185 }
186 if(_vwgt) {
187 delete[] _vwgt;
188 _vwgt = nullptr;
189 }
190 }
clearDualGraph()191 void clearDualGraph()
192 {
193 if(_xadj) {
194 delete[] _xadj;
195 _xadj = nullptr;
196 }
197 if(_adjncy) {
198 delete[] _adjncy;
199 _adjncy = nullptr;
200 }
201 }
eraseVertex()202 void eraseVertex()
203 {
204 for(std::size_t i = 0; i < _vertex.size(); i++) _vertex[i] = -1;
205 }
206 std::vector<std::set<MElement *, MElementPtrLessThan> >
getBoundaryElements(idx_t size=0)207 getBoundaryElements(idx_t size = 0)
208 {
209 std::vector<std::set<MElement *, MElementPtrLessThan> > elements(
210 (size ? size : _nparts), std::set<MElement *, MElementPtrLessThan>());
211 for(std::size_t i = 0; i < _ne; i++) {
212 for(idx_t j = _xadj[i]; j < _xadj[i + 1]; j++) {
213 if(_partition[i] != _partition[_adjncy[j]]) {
214 if(_element[i]->getDim() == _dim) {
215 elements[_partition[i]].insert(_element[i]);
216 }
217 }
218 }
219 }
220
221 return elements;
222 }
createGhostEntities()223 std::vector<GEntity *> createGhostEntities() {
224 std::vector<GEntity *> ghostEntities(_nparts, (GEntity *)nullptr);
225 int elementaryNumber = _model->getMaxElementaryNumber(_dim);
226 for(std::size_t i = 1; i <= _nparts; i++) {
227 switch(_dim) {
228 case 1:
229 ghostEntities[i - 1] = new ghostEdge(_model, ++elementaryNumber, i);
230 _model->add(static_cast<ghostEdge *>(ghostEntities[i - 1]));
231 break;
232 case 2:
233 ghostEntities[i - 1] = new ghostFace(_model, ++elementaryNumber, i);
234 _model->add(static_cast<ghostFace *>(ghostEntities[i - 1]));
235 break;
236 case 3:
237 ghostEntities[i - 1] = new ghostRegion(_model, ++elementaryNumber, i);
238 _model->add(static_cast<ghostRegion *>(ghostEntities[i - 1]));
239 break;
240 default: break;
241 }
242 }
243 return ghostEntities;
244 }
assignGhostCells()245 void assignGhostCells()
246 {
247 std::vector<GEntity *> ghostEntities = createGhostEntities();
248 for(std::size_t i = 0; i < _ne; i++) {
249 std::set<int> ghostCellsPartition;
250 for(idx_t j = _xadj[i]; j < _xadj[i + 1]; j++) {
251 if(_partition[i] != _partition[_adjncy[j]] &&
252 ghostCellsPartition.find(_partition[_adjncy[j]]) ==
253 ghostCellsPartition.end()) {
254 if(_element[i]->getDim() == _dim) {
255 switch(_dim) {
256 case 1:
257 static_cast<ghostEdge *>(ghostEntities[_partition[_adjncy[j]]])
258 ->addElement(_element[i]->getType(), _element[i],
259 _partition[i] + 1);
260 break;
261 case 2:
262 static_cast<ghostFace *>(ghostEntities[_partition[_adjncy[j]]])
263 ->addElement(_element[i]->getType(), _element[i],
264 _partition[i] + 1);
265 break;
266 case 3:
267 static_cast<ghostRegion *>(ghostEntities[_partition[_adjncy[j]]])
268 ->addElement(_element[i]->getType(), _element[i],
269 _partition[i] + 1);
270 break;
271 default: break;
272 }
273 ghostCellsPartition.insert(_partition[_adjncy[j]]);
274 }
275 }
276 }
277 }
278 }
createDualGraph(bool connectedAll)279 void createDualGraph(bool connectedAll)
280 {
281 std::vector<idx_t> nptr(_nn + 1, 0);
282 std::vector<idx_t> nind(_eptr[_ne], 0);
283
284 for(std::size_t i = 0; i < _ne; i++) {
285 for(idx_t j = _eptr[i]; j < _eptr[i + 1]; j++) { nptr[_eind[j]]++; }
286 }
287
288 for(std::size_t i = 1; i < _nn; i++) nptr[i] += nptr[i - 1];
289 for(std::size_t i = _nn; i > 0; i--) nptr[i] = nptr[i - 1];
290 nptr[0] = 0;
291
292 for(std::size_t i = 0; i < _ne; i++) {
293 for(idx_t j = _eptr[i]; j < _eptr[i + 1]; j++) {
294 nind[nptr[_eind[j]]++] = i;
295 }
296 }
297
298 for(std::size_t i = _nn; i > 0; i--) nptr[i] = nptr[i - 1];
299 nptr[0] = 0;
300
301 _xadj = new idx_t[_ne + 1];
302 for(std::size_t i = 0; i < _ne + 1; i++) _xadj[i] = 0;
303
304 std::vector<idx_t> nbrs(_ne, 0);
305 std::vector<idx_t> marker(_ne, 0);
306 for(std::size_t i = 0; i < _ne; i++) {
307 std::size_t l = 0;
308 for(idx_t j = _eptr[i]; j < _eptr[i + 1]; j++) {
309 for(idx_t k = nptr[_eind[j]]; k < nptr[_eind[j] + 1]; k++) {
310 if(nind[k] != (idx_t)i) {
311 if(marker[nind[k]] == 0) nbrs[l++] = nind[k];
312 marker[nind[k]]++;
313 }
314 }
315 }
316
317 std::size_t nbrsNeighbors = 0;
318 for(std::size_t j = 0; j < l; j++) {
319 if(marker[nbrs[j]] >=
320 (connectedAll ?
321 1 :
322 _element[i]->numCommonNodesInDualGraph(_element[nbrs[j]])))
323 nbrsNeighbors++;
324 marker[nbrs[j]] = 0;
325 nbrs[j] = 0;
326 }
327
328 _xadj[i] = nbrsNeighbors;
329 }
330
331 for(std::size_t i = 1; i < _ne; i++) _xadj[i] = _xadj[i] + _xadj[i - 1];
332 for(std::size_t i = _ne; i > 0; i--) _xadj[i] = _xadj[i - 1];
333 _xadj[0] = 0;
334
335 _adjncy = new idx_t[_xadj[_ne]];
336 for(idx_t i = 0; i < _xadj[_ne]; i++) _adjncy[i] = 0;
337
338 for(std::size_t i = 0; i < _ne; i++) {
339 std::size_t l = 0;
340 for(idx_t j = _eptr[i]; j < _eptr[i + 1]; j++) {
341 for(idx_t k = nptr[_eind[j]]; k < nptr[_eind[j] + 1]; k++) {
342 if(nind[k] != (idx_t)i) {
343 if(marker[nind[k]] == 0) nbrs[l++] = nind[k];
344 marker[nind[k]]++;
345 }
346 }
347 }
348
349 for(std::size_t j = 0; j < l; j++) {
350 if(marker[nbrs[j]] >=
351 (connectedAll ?
352 1 :
353 _element[i]->numCommonNodesInDualGraph(_element[nbrs[j]]))) {
354 _adjncy[_xadj[i]] = nbrs[j];
355 _xadj[i] = _xadj[i] + 1;
356 }
357 marker[nbrs[j]] = 0;
358 nbrs[j] = 0;
359 }
360 }
361 for(std::size_t i = _ne; i > 0; i--) _xadj[i] = _xadj[i - 1];
362 _xadj[0] = 0;
363 }
fillDefaultWeights()364 void fillDefaultWeights()
365 {
366 if(CTX::instance()->mesh.partitionLinWeight == 1 &&
367 CTX::instance()->mesh.partitionTriWeight == 1 &&
368 CTX::instance()->mesh.partitionQuaWeight == 1 &&
369 CTX::instance()->mesh.partitionTetWeight == 1 &&
370 CTX::instance()->mesh.partitionPyrWeight == 1 &&
371 CTX::instance()->mesh.partitionPriWeight == 1 &&
372 CTX::instance()->mesh.partitionHexWeight == 1)
373 return;
374
375 _vwgt = new idx_t[_ne];
376 if(CTX::instance()->mesh.partitionLinWeight == -1 ||
377 CTX::instance()->mesh.partitionTriWeight == -1 ||
378 CTX::instance()->mesh.partitionQuaWeight == -1 ||
379 CTX::instance()->mesh.partitionTetWeight == -1 ||
380 CTX::instance()->mesh.partitionPyrWeight == -1 ||
381 CTX::instance()->mesh.partitionPriWeight == -1 ||
382 CTX::instance()->mesh.partitionHexWeight == -1) {
383 for(std::size_t i = 0; i < _ne; i++) {
384 if(!_element[i]) { _vwgt[i] = 1; }
385 else {
386 _vwgt[i] = (_element[i]->getDim() == _dim ? 1 : 0);
387 }
388 }
389 }
390 else {
391 for(std::size_t i = 0; i < _ne; i++) {
392 if(!_element[i]) { _vwgt[i] = 1; }
393 else {
394 switch(_element[i]->getType()) {
395 case TYPE_LIN:
396 _vwgt[i] = CTX::instance()->mesh.partitionLinWeight;
397 break;
398 case TYPE_TRI:
399 _vwgt[i] = CTX::instance()->mesh.partitionTriWeight;
400 break;
401 case TYPE_QUA:
402 _vwgt[i] = CTX::instance()->mesh.partitionQuaWeight;
403 break;
404 case TYPE_TET:
405 _vwgt[i] = CTX::instance()->mesh.partitionTetWeight;
406 break;
407 case TYPE_PYR:
408 _vwgt[i] = CTX::instance()->mesh.partitionPyrWeight;
409 break;
410 case TYPE_PRI:
411 _vwgt[i] = CTX::instance()->mesh.partitionPriWeight;
412 break;
413 case TYPE_HEX:
414 _vwgt[i] = CTX::instance()->mesh.partitionHexWeight;
415 break;
416 default: _vwgt[i] = 1; break;
417 }
418 }
419 }
420 }
421 }
422 };
423
424 template <class ITERATOR>
fillElementsToNodesMap(Graph & graph,GEntity * entity,idx_t & eptrIndex,idx_t & eindIndex,idx_t & numVertex,ITERATOR it_beg,ITERATOR it_end)425 static void fillElementsToNodesMap(Graph &graph, GEntity *entity,
426 idx_t &eptrIndex, idx_t &eindIndex,
427 idx_t &numVertex, ITERATOR it_beg,
428 ITERATOR it_end)
429 {
430 for(ITERATOR it = it_beg; it != it_end; ++it) {
431 const std::size_t numVertices = (*it)->getNumPrimaryVertices();
432 graph.element(eptrIndex++, *it);
433 graph.eptr(eptrIndex, graph.eptr(eptrIndex - 1) + numVertices);
434 for(std::size_t i = 0; i < numVertices; i++) {
435 if(graph.vertex((*it)->getVertex(i)->getNum() - 1) == -1) {
436 graph.vertex((*it)->getVertex(i)->getNum() - 1, numVertex++);
437 }
438 graph.eind(eindIndex, graph.vertex((*it)->getVertex(i)->getNum() - 1));
439 eindIndex++;
440 }
441 }
442 }
443
getSizeOfEind(GModel * model)444 static std::size_t getSizeOfEind(GModel *model)
445 {
446 std::size_t size = 0;
447 // Loop over volumes
448 for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
449 size += 4 * (*it)->tetrahedra.size();
450 size += 8 * (*it)->hexahedra.size();
451 size += 6 * (*it)->prisms.size();
452 size += 5 * (*it)->pyramids.size();
453 size += 4 * (*it)->trihedra.size();
454 }
455
456 // Loop over surfaces
457 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
458 size += 3 * (*it)->triangles.size();
459 size += 4 * (*it)->quadrangles.size();
460 }
461
462 // Loop over curves
463 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
464 size += 2 * (*it)->lines.size();
465 }
466
467 // Loop over points
468 for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
469 size += 1 * (*it)->points.size();
470 }
471
472 return size;
473 }
474
475 // Creates a mesh data structure used by Metis routines. Returns: 0 = success, 1
476 // = no elements found, 2 = error.
makeGraph(GModel * model,Graph & graph,int selectDim)477 static int makeGraph(GModel *model, Graph &graph, int selectDim)
478 {
479 std::size_t eindSize = 0;
480 if(selectDim < 0) {
481 graph.ne(model->getNumMeshElements());
482 graph.nn(model->getNumMeshVertices());
483 graph.dim(model->getMeshDim());
484 graph.elementResize(graph.ne());
485 graph.vertexResize(model->getMaxVertexNumber());
486 graph.eptrResize(graph.ne() + 1);
487 graph.eptr(0, 0);
488 eindSize = getSizeOfEind(model);
489 graph.eindResize(eindSize);
490 }
491 else {
492 GModel *tmp = new GModel();
493 std::vector<GEntity *> entities;
494 model->getEntities(entities);
495
496 std::set<MVertex *> vertices;
497 for(std::size_t i = 0; i < entities.size(); i++) {
498 if(entities[i]->dim() == selectDim) {
499 switch(entities[i]->dim()) {
500 case 3: tmp->add(static_cast<GRegion *>(entities[i])); break;
501 case 2: tmp->add(static_cast<GFace *>(entities[i])); break;
502 case 1: tmp->add(static_cast<GEdge *>(entities[i])); break;
503 case 0: tmp->add(static_cast<GVertex *>(entities[i])); break;
504 default: break;
505 }
506 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
507 for(std::size_t k = 0;
508 k < entities[i]->getMeshElement(j)->getNumVertices(); k++) {
509 vertices.insert(entities[i]->getMeshElement(j)->getVertex(k));
510 }
511 }
512 }
513 }
514
515 graph.ne(tmp->getNumMeshElements());
516 graph.nn(vertices.size());
517 graph.dim(tmp->getMeshDim());
518 graph.elementResize(graph.ne());
519 graph.vertexResize(model->getMaxVertexNumber());
520 graph.eptrResize(graph.ne() + 1);
521 graph.eptr(0, 0);
522 eindSize = getSizeOfEind(tmp);
523 graph.eindResize(eindSize);
524
525 tmp->remove();
526 delete tmp;
527 }
528
529 idx_t eptrIndex = 0;
530 idx_t eindIndex = 0;
531 idx_t numVertex = 0;
532
533 if(graph.ne() == 0) {
534 Msg::Error("No mesh elements were found");
535 return 1;
536 }
537 if(graph.dim() == 0) {
538 Msg::Error("Cannot partition a point");
539 return 1;
540 }
541
542 // Loop over volumes
543 if(selectDim < 0 || selectDim == 3) {
544 for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
545 GRegion *r = *it;
546 fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
547 r->tetrahedra.begin(), r->tetrahedra.end());
548 fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
549 r->hexahedra.begin(), r->hexahedra.end());
550 fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
551 r->prisms.begin(), r->prisms.end());
552 fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
553 r->pyramids.begin(), r->pyramids.end());
554 fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
555 r->trihedra.begin(), r->trihedra.end());
556 }
557 }
558
559 // Loop over surfaces
560 if(selectDim < 0 || selectDim == 2) {
561 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
562 GFace *f = *it;
563 fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
564 f->triangles.begin(), f->triangles.end());
565 fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
566 f->quadrangles.begin(), f->quadrangles.end());
567 }
568 }
569
570 // Loop over curves
571 if(selectDim < 0 || selectDim == 1) {
572 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
573 GEdge *e = *it;
574 fillElementsToNodesMap(graph, e, eptrIndex, eindIndex, numVertex,
575 e->lines.begin(), e->lines.end());
576 }
577 }
578
579 // Loop over points
580 if(selectDim < 0 || selectDim == 0) {
581 for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
582 GVertex *v = *it;
583 fillElementsToNodesMap(graph, v, eptrIndex, eindIndex, numVertex,
584 v->points.begin(), v->points.end());
585 }
586 }
587
588 return 0;
589 }
590
591 // Partition a graph created by makeGraph using Metis library. Returns: 0 =
592 // success, 1 = error, 2 = exception thrown.
partitionGraph(Graph & graph,bool verbose)593 static int partitionGraph(Graph &graph, bool verbose)
594 {
595 #ifdef HAVE_METIS
596 std::stringstream opt;
597 try {
598 idx_t metisOptions[METIS_NOPTIONS];
599 METIS_SetDefaultOptions(metisOptions);
600
601 opt << 8 * sizeof(idx_t) << " bit indices";
602
603 opt << ", ptype:";
604 switch(CTX::instance()->mesh.metisAlgorithm) {
605 case 1: // Recursive
606 metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_RB;
607 opt << "rb";
608 break;
609 case 2: // K-way
610 metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
611 opt << "kway";
612 break;
613 default: opt << "default"; break;
614 }
615
616 opt << ", ufactor:";
617 if(CTX::instance()->mesh.metisMaxLoadImbalance >= 0) {
618 metisOptions[METIS_OPTION_UFACTOR] =
619 CTX::instance()->mesh.metisMaxLoadImbalance;
620 opt << CTX::instance()->mesh.metisMaxLoadImbalance;
621 }
622 else {
623 opt << "default";
624 }
625
626 opt << ", ctype:";
627 switch(CTX::instance()->mesh.metisEdgeMatching) {
628 case 1: // Random matching
629 metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_RM;
630 opt << "rm";
631 break;
632 case 2: // Sorted heavy-edge matching
633 metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM;
634 opt << "shem";
635 break;
636 default: opt << "default"; break;
637 }
638
639 opt << ", rtype:";
640 switch(CTX::instance()->mesh.metisRefinementAlgorithm) {
641 case 1: // FM-based cut refinement
642 metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_FM;
643 opt << "fm";
644 break;
645 case 2: // Greedy boundary refinement
646 metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY;
647 opt << "greedy";
648 break;
649 case 3: // Two-sided node FM refinement
650 metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP2SIDED;
651 opt << "sep2sided";
652 break;
653 case 4: // One-sided node FM refinement
654 metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP1SIDED;
655 opt << "sep1sided";
656 break;
657 default: opt << "default"; break;
658 }
659
660 opt << ", objtype:";
661 switch(CTX::instance()->mesh.metisObjective) {
662 case 1: // Min. cut
663 metisOptions[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
664 opt << "cut";
665 break;
666 case 2: // Min. communication volume (slower)
667 metisOptions[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
668 opt << "vol";
669 break;
670 default: opt << "default"; break;
671 }
672
673 opt << ", minconn:";
674 switch(CTX::instance()->mesh.metisMinConn) {
675 case 0:
676 metisOptions[METIS_OPTION_MINCONN] = 0;
677 opt << 0;
678 break;
679 case 1:
680 metisOptions[METIS_OPTION_MINCONN] = 1;
681 opt << 1;
682 break;
683 default: opt << "default"; break;
684 }
685
686 if(verbose) Msg::Info("Running METIS with %s", opt.str().c_str());
687
688 // C numbering
689 metisOptions[METIS_OPTION_NUMBERING] = 0;
690
691 idx_t objval;
692 std::vector<idx_t> epart(graph.ne());
693 idx_t ne = graph.ne();
694 idx_t numPart = graph.nparts();
695 idx_t ncon = 1;
696 graph.fillDefaultWeights();
697
698 int metisError = 0;
699 graph.createDualGraph(false);
700
701 if(metisOptions[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY) {
702 metisError = METIS_PartGraphKway(
703 &ne, &ncon, graph.xadj(), graph.adjncy(), graph.vwgt(), nullptr,
704 nullptr, &numPart, nullptr, nullptr, metisOptions, &objval, &epart[0]);
705 }
706 else {
707 metisError = METIS_PartGraphRecursive(
708 &ne, &ncon, graph.xadj(), graph.adjncy(), graph.vwgt(), nullptr,
709 nullptr, &numPart, nullptr, nullptr, metisOptions, &objval, &epart[0]);
710 }
711
712 switch(metisError) {
713 case METIS_OK: break;
714 case METIS_ERROR_INPUT: Msg::Error("METIS input error"); return 1;
715 case METIS_ERROR_MEMORY: Msg::Error("METIS memory error"); return 1;
716 case METIS_ERROR:
717 default: Msg::Error("METIS error"); return 1;
718 }
719
720 // Check and correct the topology
721 for(int i = 1; i < 4; i++) {
722 for(std::size_t j = 0; j < graph.ne(); j++) {
723 if(graph.element(j)->getDim() == graph.dim()) continue;
724
725 for(idx_t k = graph.xadj(j); k < graph.xadj(j + 1); k++) {
726 if(graph.element(j)->getDim() ==
727 graph.element(graph.adjncy(k))->getDim() - i) {
728 if(epart[j] != epart[graph.adjncy(k)]) {
729 epart[j] = epart[graph.adjncy(k)];
730 break;
731 }
732 }
733 }
734 }
735 }
736 graph.partition(epart);
737 if(verbose) Msg::Info("%d partitions, %d total edge-cuts", numPart, objval);
738 } catch(...) {
739 Msg::Error("METIS exception");
740 return 2;
741 }
742 #endif
743
744 return 0;
745 }
746
747 template <class ENTITY, class ITERATOR>
748 static void
assignElementsToEntities(GModel * model,hashmapelementpart & elmToPartition,std::vector<ENTITY * > & newEntities,ITERATOR it_beg,ITERATOR it_end,int & elementaryNumber)749 assignElementsToEntities(GModel *model, hashmapelementpart &elmToPartition,
750 std::vector<ENTITY *> &newEntities, ITERATOR it_beg,
751 ITERATOR it_end, int &elementaryNumber)
752 {
753 for(ITERATOR it = it_beg; it != it_end; ++it) {
754 int partition = elmToPartition[(*it)] - 1;
755
756 if(!newEntities[partition]) {
757 std::vector<int> partitions;
758 partitions.push_back(partition + 1);
759 ENTITY *de = new ENTITY(model, ++elementaryNumber, partitions);
760 model->add(de);
761 newEntities[partition] = de;
762 }
763
764 newEntities[partition]->addElement((*it)->getType(), *it);
765 }
766 }
767
768 template <class ITERATOR>
setVerticesToEntity(GEntity * entity,ITERATOR it_beg,ITERATOR it_end)769 static void setVerticesToEntity(GEntity *entity, ITERATOR it_beg,
770 ITERATOR it_end)
771 {
772 for(ITERATOR it = it_beg; it != it_end; ++it) {
773 for(std::size_t i = 0; i < (*it)->getNumVertices(); i++) {
774 if(!(*it)->getVertex(i)->onWhat()) {
775 (*it)->getVertex(i)->setEntity(entity);
776 entity->addMeshVertex((*it)->getVertex(i));
777 }
778 }
779 }
780 }
781
782 template <class ITERATOR>
removeVerticesEntity(ITERATOR it_beg,ITERATOR it_end)783 static void removeVerticesEntity(ITERATOR it_beg, ITERATOR it_end)
784 {
785 for(ITERATOR it = it_beg; it != it_end; ++it) {
786 for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) {
787 for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices();
788 j++) {
789 (*it)->getMeshElement(i)->getVertex(j)->setEntity(nullptr);
790 }
791 }
792 (*it)->mesh_vertices.clear();
793 }
794 }
795
796 // Assign the vertices to its corresponding entity
assignMeshVertices(GModel * model)797 static void assignMeshVertices(GModel *model)
798 {
799 removeVerticesEntity(model->firstVertex(), model->lastVertex());
800 removeVerticesEntity(model->firstEdge(), model->lastEdge());
801 removeVerticesEntity(model->firstFace(), model->lastFace());
802 removeVerticesEntity(model->firstRegion(), model->lastRegion());
803
804 // Loop over points
805 for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
806 setVerticesToEntity(*it, (*it)->points.begin(), (*it)->points.end());
807 }
808
809 // Loop over curves
810 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
811 setVerticesToEntity(*it, (*it)->lines.begin(), (*it)->lines.end());
812 }
813
814 // Loop over surfaces
815 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
816 setVerticesToEntity(*it, (*it)->triangles.begin(), (*it)->triangles.end());
817 setVerticesToEntity(*it, (*it)->quadrangles.begin(),
818 (*it)->quadrangles.end());
819 }
820
821 // Loop over volumes
822 for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
823 setVerticesToEntity(*it, (*it)->tetrahedra.begin(),
824 (*it)->tetrahedra.end());
825 setVerticesToEntity(*it, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
826 setVerticesToEntity(*it, (*it)->prisms.begin(), (*it)->prisms.end());
827 setVerticesToEntity(*it, (*it)->pyramids.begin(), (*it)->pyramids.end());
828 setVerticesToEntity(*it, (*it)->trihedra.begin(), (*it)->trihedra.end());
829 }
830 }
831
fillConnectedElements(std::vector<std::set<MElement *,MElementPtrLessThan>> & connectedElements,Graph & graph)832 static void fillConnectedElements(
833 std::vector<std::set<MElement *, MElementPtrLessThan> > &connectedElements,
834 Graph &graph)
835 {
836 std::stack<idx_t> elementStack;
837 std::set<MElement *, MElementPtrLessThan> elements;
838 idx_t startElement = 0;
839 bool stop = true;
840 std::size_t size = 0;
841 idx_t isolatedElements = 0;
842
843 do {
844 // Inititalization
845 elementStack.push(startElement);
846 elements.insert(graph.element(startElement));
847
848 while(elementStack.size() != 0) {
849 idx_t top = elementStack.top();
850 elementStack.pop();
851 elements.insert(graph.element(top));
852
853 for(idx_t i = graph.xadj(top); i < graph.xadj(top + 1); i++) {
854 if(graph.adjncy(i) != 0) {
855 elementStack.push(graph.adjncy(i));
856 graph.adjncy(i, 0);
857 }
858 }
859 }
860 connectedElements.push_back(elements);
861 size += elements.size();
862 elements.clear();
863
864 stop = (size == graph.ne() ? true : false);
865
866 startElement = 0;
867 if(!stop) {
868 for(std::size_t i = 0; i < graph.ne(); i++) {
869 for(idx_t j = graph.xadj(i); j < graph.xadj(i + 1); j++) {
870 if(graph.adjncy(j) != 0) {
871 startElement = i;
872 i = graph.ne();
873 break;
874 }
875 }
876 }
877 if(startElement == 0) {
878 idx_t skipIsolatedElements = 0;
879 for(std::size_t i = 1; i < graph.ne(); i++) {
880 if(graph.xadj(i) == graph.xadj(i + 1)) {
881 if(skipIsolatedElements == isolatedElements) {
882 startElement = i;
883 isolatedElements++;
884 break;
885 }
886 skipIsolatedElements++;
887 }
888 }
889 }
890 }
891 } while(!stop);
892 }
893
894 static bool
divideNonConnectedEntities(GModel * model,int dim,std::set<GRegion *,GEntityPtrLessThan> & regions,std::set<GFace *,GEntityPtrLessThan> & faces,std::set<GEdge *,GEntityPtrLessThan> & edges,std::set<GVertex *,GEntityPtrLessThan> & vertices)895 divideNonConnectedEntities(GModel *model, int dim,
896 std::set<GRegion *, GEntityPtrLessThan> ®ions,
897 std::set<GFace *, GEntityPtrLessThan> &faces,
898 std::set<GEdge *, GEntityPtrLessThan> &edges,
899 std::set<GVertex *, GEntityPtrLessThan> &vertices)
900 {
901 bool ret = false;
902 // Loop over points
903 if(dim < 0 || dim == 0) {
904 int elementaryNumber = model->getMaxElementaryNumber(0);
905
906 for(auto it = vertices.begin(); it != vertices.end(); ++it) {
907 if((*it)->geomType() == GEntity::PartitionPoint) {
908 partitionVertex *vertex = static_cast<partitionVertex *>(*it);
909
910 if(vertex->getNumMeshElements() > 1) {
911 ret = true;
912 for(std::size_t i = 0; i < vertex->getNumMeshElements(); i++) {
913 // Create the new partitionVertex
914 partitionVertex *pvertex = new partitionVertex(
915 model, ++elementaryNumber, vertex->getPartitions());
916 // Assign parent entity
917 pvertex->setParentEntity(vertex->getParentEntity());
918 // Add to model
919 model->add(pvertex);
920 // Add elements
921 pvertex->addElement(vertex->getMeshElement(i)->getType(),
922 vertex->getMeshElement(i));
923 // Move B-Rep
924 std::vector<GEdge *> BRepEdges = vertex->edges();
925 if(!BRepEdges.empty()) {
926 for(auto itBRep = BRepEdges.begin(); itBRep != BRepEdges.end();
927 ++itBRep) {
928 if(vertex == (*itBRep)->getBeginVertex()) {
929 (*itBRep)->setVertex(pvertex, 1);
930 pvertex->addEdge(*itBRep);
931 }
932 if(vertex == (*itBRep)->getEndVertex()) {
933 (*itBRep)->setVertex(pvertex, -1);
934 pvertex->addEdge(*itBRep);
935 }
936 }
937 }
938 }
939
940 model->remove(vertex);
941 vertex->points.clear();
942 vertex->mesh_vertices.clear();
943 delete vertex;
944 }
945 }
946 }
947 }
948
949 // Loop over curves
950 if(dim < 0 || dim == 1) {
951 // We build a graph
952 Graph graph(model);
953 graph.ne(model->getNumMeshElements(1));
954 graph.nn(model->getNumMeshVertices(1));
955 graph.dim(model->getMeshDim());
956 graph.elementResize(graph.ne());
957 graph.vertexResize(model->getMaxVertexNumber());
958 graph.eptrResize(graph.ne() + 1);
959 graph.eptr(0, 0);
960 std::size_t eindSize = getSizeOfEind(model);
961 graph.eindResize(eindSize);
962
963 int elementaryNumber = model->getMaxElementaryNumber(1);
964
965 for(auto it = edges.begin(); it != edges.end(); ++it) {
966 if((*it)->geomType() == GEntity::PartitionCurve) {
967 partitionEdge *edge = static_cast<partitionEdge *>(*it);
968
969 graph.ne(edge->getNumMeshElements());
970 graph.dim(1);
971 graph.eptr(0, 0);
972 graph.clearDualGraph();
973 graph.eraseVertex();
974
975 idx_t eptrIndex = 0;
976 idx_t eindIndex = 0;
977 idx_t numVertex = 0;
978
979 fillElementsToNodesMap(graph, edge, eptrIndex, eindIndex, numVertex,
980 edge->lines.begin(), edge->lines.end());
981 graph.nn(numVertex);
982 graph.createDualGraph(false);
983
984 // if a graph contains more than ((n-1)*(n-2))/2 edges (where n is the
985 // number of nodes), then it is connected.
986 if(((graph.numNodes() - 1) * (graph.numNodes() - 2)) / 2 <
987 graph.numEdges()) {
988 continue;
989 }
990
991 std::vector<std::set<MElement *, MElementPtrLessThan> >
992 connectedElements;
993 fillConnectedElements(connectedElements, graph);
994
995 if(connectedElements.size() > 1) {
996 ret = true;
997 std::vector<GFace *> BRepFaces = edge->faces();
998
999 std::vector<int> oldOrientations;
1000 oldOrientations.reserve(BRepFaces.size());
1001
1002 if(!BRepFaces.empty()) {
1003 for(auto itBRep = BRepFaces.begin(); itBRep != BRepFaces.end();
1004 ++itBRep) {
1005 oldOrientations.push_back((*itBRep)->delEdge(edge));
1006 }
1007 }
1008
1009 for(std::size_t i = 0; i < connectedElements.size(); i++) {
1010 // Create the new partitionEdge
1011 partitionEdge *pedge =
1012 new partitionEdge(model, ++elementaryNumber, nullptr, nullptr,
1013 edge->getPartitions());
1014 // Assign parent entity
1015 pedge->setParentEntity(edge->getParentEntity());
1016 // Add to model
1017 model->add(pedge);
1018 for(auto itSet = connectedElements[i].begin();
1019 itSet != connectedElements[i].end(); ++itSet) {
1020 // Add elements
1021 pedge->addElement((*itSet)->getType(), (*itSet));
1022 }
1023 // Move B-Rep
1024 if(BRepFaces.size() > 0) {
1025 std::size_t i = 0;
1026 for(auto itBRep = BRepFaces.begin(); itBRep != BRepFaces.end();
1027 ++itBRep) {
1028 (*itBRep)->setEdge(pedge, oldOrientations[i]);
1029 pedge->addFace(*itBRep);
1030 i++;
1031 }
1032 }
1033 }
1034
1035 model->remove(edge);
1036 edge->lines.clear();
1037 edge->mesh_vertices.clear();
1038 delete edge;
1039 }
1040
1041 connectedElements.clear();
1042 }
1043 }
1044 }
1045
1046 // Loop over surfaces
1047 if(dim < 0 || dim == 2) {
1048 // We build a graph
1049 Graph graph(model);
1050 graph.ne(model->getNumMeshElements(2));
1051 graph.nn(model->getNumMeshVertices(2));
1052 graph.dim(model->getMeshDim());
1053 graph.elementResize(graph.ne());
1054 graph.vertexResize(model->getMaxVertexNumber());
1055 graph.eptrResize(graph.ne() + 1);
1056 graph.eptr(0, 0);
1057 std::size_t eindSize = getSizeOfEind(model);
1058 graph.eindResize(eindSize);
1059
1060 int elementaryNumber = model->getMaxElementaryNumber(2);
1061
1062 for(auto it = faces.begin(); it != faces.end(); ++it) {
1063 if((*it)->geomType() == GEntity::PartitionSurface) {
1064 partitionFace *face = static_cast<partitionFace *>(*it);
1065
1066 graph.ne(face->getNumMeshElements());
1067 graph.dim(2);
1068 graph.eptr(0, 0);
1069 graph.clearDualGraph();
1070 graph.eraseVertex();
1071
1072 idx_t eptrIndex = 0;
1073 idx_t eindIndex = 0;
1074 idx_t numVertex = 0;
1075
1076 fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
1077 face->triangles.begin(), face->triangles.end());
1078 fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
1079 face->quadrangles.begin(),
1080 face->quadrangles.end());
1081 graph.nn(numVertex);
1082 graph.createDualGraph(false);
1083
1084 // if a graph contains more than ((n-1)*(n-2))/2 edges
1085 // (where n is the number of nodes), then it is connected.
1086 if(((graph.numNodes() - 1) * (graph.numNodes() - 2)) / 2 <
1087 graph.numEdges()) {
1088 continue;
1089 }
1090
1091 std::vector<std::set<MElement *, MElementPtrLessThan> >
1092 connectedElements;
1093 fillConnectedElements(connectedElements, graph);
1094
1095 if(connectedElements.size() > 1) {
1096 ret = true;
1097 std::list<GRegion *> BRepRegions = face->regions();
1098 std::vector<int> oldOrientations;
1099 if(BRepRegions.size() > 0) {
1100 for(auto itBRep = BRepRegions.begin(); itBRep != BRepRegions.end();
1101 ++itBRep) {
1102 oldOrientations.push_back((*itBRep)->delFace(face));
1103 }
1104 }
1105
1106 for(std::size_t i = 0; i < connectedElements.size(); i++) {
1107 // Create the new partitionFace
1108 partitionFace *pface = new partitionFace(model, ++elementaryNumber,
1109 face->getPartitions());
1110 // Assign parent entity
1111 pface->setParentEntity(face->getParentEntity());
1112 // Add to model
1113 model->add(pface);
1114 for(auto itSet = connectedElements[i].begin();
1115 itSet != connectedElements[i].end(); ++itSet) {
1116 // Add elements
1117 pface->addElement((*itSet)->getType(), (*itSet));
1118 }
1119 // Move B-Rep
1120 if(BRepRegions.size() > 0) {
1121 std::size_t i = 0;
1122 for(auto itBRep = BRepRegions.begin();
1123 itBRep != BRepRegions.end(); ++itBRep) {
1124 (*itBRep)->setFace(pface, oldOrientations[i]);
1125 pface->addRegion(*itBRep);
1126 i++;
1127 }
1128 }
1129 }
1130
1131 model->remove(face);
1132 face->triangles.clear();
1133 face->quadrangles.clear();
1134 face->mesh_vertices.clear();
1135 delete face;
1136 }
1137
1138 connectedElements.clear();
1139 }
1140 }
1141 }
1142
1143 // Loop over volumes
1144 if(dim < 0 || dim == 3) {
1145 // We build a graph
1146 Graph graph(model);
1147 graph.ne(model->getNumMeshElements(3));
1148 graph.nn(model->getNumMeshVertices(3));
1149 graph.dim(model->getMeshDim());
1150 graph.elementResize(graph.ne());
1151 graph.vertexResize(model->getMaxVertexNumber());
1152 graph.eptrResize(graph.ne() + 1);
1153 graph.eptr(0, 0);
1154 std::size_t eindSize = getSizeOfEind(model);
1155 graph.eindResize(eindSize);
1156
1157 int elementaryNumber = model->getMaxElementaryNumber(3);
1158
1159 for(auto it = regions.begin(); it != regions.end(); ++it) {
1160 if((*it)->geomType() == GEntity::PartitionVolume) {
1161 partitionRegion *region = static_cast<partitionRegion *>(*it);
1162
1163 graph.ne(region->getNumMeshElements());
1164 graph.dim(3);
1165 graph.eptr(0, 0);
1166 graph.clearDualGraph();
1167 graph.eraseVertex();
1168
1169 idx_t eptrIndex = 0;
1170 idx_t eindIndex = 0;
1171 idx_t numVertex = 0;
1172
1173 fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
1174 region->tetrahedra.begin(),
1175 region->tetrahedra.end());
1176 fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
1177 region->hexahedra.begin(),
1178 region->hexahedra.end());
1179 fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
1180 region->prisms.begin(), region->prisms.end());
1181 fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
1182 region->pyramids.begin(),
1183 region->pyramids.end());
1184 fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
1185 region->trihedra.begin(),
1186 region->trihedra.end());
1187 graph.nn(numVertex);
1188 graph.createDualGraph(false);
1189
1190 // if a graph contains more than ((n-1)*(n-2))/2 edges (where n is the
1191 // number of nodes), then it is connected.
1192 if(((graph.numNodes() - 1) * (graph.numNodes() - 2)) / 2 <
1193 graph.numEdges()) {
1194 continue;
1195 }
1196
1197 std::vector<std::set<MElement *, MElementPtrLessThan> >
1198 connectedElements;
1199 fillConnectedElements(connectedElements, graph);
1200
1201 if(connectedElements.size() > 1) {
1202 ret = true;
1203 for(std::size_t i = 0; i < connectedElements.size(); i++) {
1204 // Create the new partitionRegion
1205 partitionRegion *pregion = new partitionRegion(
1206 model, ++elementaryNumber, region->getPartitions());
1207 // Assign d parent entity
1208 pregion->setParentEntity(region->getParentEntity());
1209 // Add to model
1210 model->add(pregion);
1211 for(auto itSet = connectedElements[i].begin();
1212 itSet != connectedElements[i].end(); ++itSet) {
1213 // Add elements
1214 pregion->addElement((*itSet)->getType(), (*itSet));
1215 }
1216 }
1217
1218 model->remove(region);
1219 region->tetrahedra.clear();
1220 region->hexahedra.clear();
1221 region->prisms.clear();
1222 region->pyramids.clear();
1223 region->trihedra.clear();
1224 region->mesh_vertices.clear();
1225 delete region;
1226 }
1227
1228 connectedElements.clear();
1229 }
1230 }
1231 }
1232
1233 return ret;
1234 }
1235
1236 // Create the new volume entities (omega)
createNewEntities(GModel * model,hashmapelementpart & elmToPartition)1237 static void createNewEntities(GModel *model, hashmapelementpart &elmToPartition)
1238 {
1239 std::set<GRegion *, GEntityPtrLessThan> regions = model->getRegions();
1240 std::set<GFace *, GEntityPtrLessThan> faces = model->getFaces();
1241 std::set<GEdge *, GEntityPtrLessThan> edges = model->getEdges();
1242 std::set<GVertex *, GEntityPtrLessThan> vertices = model->getVertices();
1243
1244 int elementaryNumber = model->getMaxElementaryNumber(0);
1245 for(auto it = vertices.begin(); it != vertices.end(); ++it) {
1246 std::vector<partitionVertex *> newVertices(model->getNumPartitions(),
1247 nullptr);
1248
1249 assignElementsToEntities(model, elmToPartition, newVertices,
1250 (*it)->points.begin(), (*it)->points.end(),
1251 elementaryNumber);
1252
1253 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1254 if(newVertices[i]) {
1255 static_cast<partitionVertex *>(newVertices[i])->setParentEntity((*it));
1256 }
1257 }
1258
1259 (*it)->mesh_vertices.clear();
1260
1261 (*it)->points.clear();
1262 }
1263
1264 elementaryNumber = model->getMaxElementaryNumber(1);
1265 for(auto it = edges.begin(); it != edges.end(); ++it) {
1266 std::vector<partitionEdge *> newEdges(model->getNumPartitions(), nullptr);
1267
1268 assignElementsToEntities(model, elmToPartition, newEdges,
1269 (*it)->lines.begin(), (*it)->lines.end(),
1270 elementaryNumber);
1271
1272 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1273 if(newEdges[i]) {
1274 static_cast<partitionEdge *>(newEdges[i])->setParentEntity(*it);
1275 }
1276 }
1277
1278 (*it)->mesh_vertices.clear();
1279
1280 (*it)->lines.clear();
1281 }
1282
1283 elementaryNumber = model->getMaxElementaryNumber(2);
1284 for(auto it = faces.begin(); it != faces.end(); ++it) {
1285 std::vector<partitionFace *> newFaces(model->getNumPartitions(), nullptr);
1286
1287 assignElementsToEntities(model, elmToPartition, newFaces,
1288 (*it)->triangles.begin(), (*it)->triangles.end(),
1289 elementaryNumber);
1290 assignElementsToEntities(model, elmToPartition, newFaces,
1291 (*it)->quadrangles.begin(),
1292 (*it)->quadrangles.end(), elementaryNumber);
1293
1294 std::list<GRegion *> BRepRegions = (*it)->regions();
1295 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1296 if(newFaces[i]) {
1297 static_cast<partitionFace *>(newFaces[i])->setParentEntity(*it);
1298 }
1299 }
1300
1301 (*it)->mesh_vertices.clear();
1302
1303 (*it)->triangles.clear();
1304 (*it)->quadrangles.clear();
1305 }
1306
1307 elementaryNumber = model->getMaxElementaryNumber(3);
1308 for(auto it = regions.begin(); it != regions.end(); ++it) {
1309 std::vector<partitionRegion *> newRegions(model->getNumPartitions(),
1310 nullptr);
1311
1312 assignElementsToEntities(model, elmToPartition, newRegions,
1313 (*it)->tetrahedra.begin(), (*it)->tetrahedra.end(),
1314 elementaryNumber);
1315 assignElementsToEntities(model, elmToPartition, newRegions,
1316 (*it)->hexahedra.begin(), (*it)->hexahedra.end(),
1317 elementaryNumber);
1318 assignElementsToEntities(model, elmToPartition, newRegions,
1319 (*it)->prisms.begin(), (*it)->prisms.end(),
1320 elementaryNumber);
1321 assignElementsToEntities(model, elmToPartition, newRegions,
1322 (*it)->pyramids.begin(), (*it)->pyramids.end(),
1323 elementaryNumber);
1324 assignElementsToEntities(model, elmToPartition, newRegions,
1325 (*it)->trihedra.begin(), (*it)->trihedra.end(),
1326 elementaryNumber);
1327
1328 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1329 if(newRegions[i]) {
1330 static_cast<partitionRegion *>(newRegions[i])->setParentEntity(*it);
1331 }
1332 }
1333
1334 (*it)->mesh_vertices.clear();
1335
1336 (*it)->tetrahedra.clear();
1337 (*it)->hexahedra.clear();
1338 (*it)->prisms.clear();
1339 (*it)->pyramids.clear();
1340 (*it)->trihedra.clear();
1341 }
1342
1343 // If we don't create the partition topology let's just assume that the user
1344 // does not care about multiply connected partitions or partition boundaries.
1345 if(!CTX::instance()->mesh.partitionCreateTopology) return;
1346 regions = model->getRegions();
1347 faces = model->getFaces();
1348 edges = model->getEdges();
1349 vertices = model->getVertices();
1350 divideNonConnectedEntities(model, -1, regions, faces, edges, vertices);
1351 }
1352
fillElementToEntity(GModel * model,hashmapelement & elmToEntity,int dim)1353 static void fillElementToEntity(GModel *model, hashmapelement &elmToEntity,
1354 int dim)
1355 {
1356 // Loop over volumes
1357 if(dim < 0 || dim == 3) {
1358 for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
1359 for(auto itElm = (*it)->tetrahedra.begin();
1360 itElm != (*it)->tetrahedra.end(); ++itElm)
1361 elmToEntity.insert(std::make_pair(*itElm, *it));
1362 for(auto itElm = (*it)->hexahedra.begin();
1363 itElm != (*it)->hexahedra.end(); ++itElm)
1364 elmToEntity.insert(std::make_pair(*itElm, *it));
1365 for(auto itElm = (*it)->prisms.begin(); itElm != (*it)->prisms.end();
1366 ++itElm)
1367 elmToEntity.insert(std::make_pair(*itElm, *it));
1368 for(auto itElm = (*it)->pyramids.begin(); itElm != (*it)->pyramids.end();
1369 ++itElm)
1370 elmToEntity.insert(std::make_pair(*itElm, *it));
1371 for(auto itElm = (*it)->trihedra.begin(); itElm != (*it)->trihedra.end();
1372 ++itElm)
1373 elmToEntity.insert(std::make_pair(*itElm, *it));
1374 }
1375 }
1376
1377 // Loop over surfaces
1378 if(dim < 0 || dim == 2) {
1379 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
1380 for(auto itElm = (*it)->triangles.begin();
1381 itElm != (*it)->triangles.end(); ++itElm)
1382 elmToEntity.insert(std::make_pair(*itElm, *it));
1383 for(auto itElm = (*it)->quadrangles.begin();
1384 itElm != (*it)->quadrangles.end(); ++itElm)
1385 elmToEntity.insert(std::make_pair(*itElm, *it));
1386 }
1387 }
1388
1389 // Loop over curves
1390 if(dim < 0 || dim == 1) {
1391 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
1392 for(auto itElm = (*it)->lines.begin(); itElm != (*it)->lines.end();
1393 ++itElm)
1394 elmToEntity.insert(std::make_pair(*itElm, *it));
1395 }
1396 }
1397
1398 // Loop over points
1399 if(dim < 0 || dim == 0) {
1400 for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
1401 for(auto itElm = (*it)->points.begin(); itElm != (*it)->points.end();
1402 ++itElm)
1403 elmToEntity.insert(std::make_pair(*itElm, *it));
1404 }
1405 }
1406 }
1407
getReferenceElement(const std::vector<std::pair<MElement *,std::vector<int>>> & elementPairs)1408 static MElement *getReferenceElement(
1409 const std::vector<std::pair<MElement *, std::vector<int> > > &elementPairs)
1410 {
1411 int min = std::numeric_limits<int>::max();
1412 std::vector<std::pair<MElement *, std::vector<int> > > minSizeElementPairs;
1413 std::vector<std::pair<MElement *, std::vector<int> > > minSizeElementPairsTmp;
1414
1415 // Take only the elements having the less partition in commun. For exemple we
1416 // take (1,2) and (3,8) but not (2,5,9) or (1,4,5,7)
1417 for(std::size_t i = 0; i < elementPairs.size(); i++) {
1418 if(min > (int)elementPairs[i].second.size()) {
1419 minSizeElementPairs.clear();
1420 min = elementPairs[i].second.size();
1421 minSizeElementPairs.push_back(elementPairs[i]);
1422 }
1423 else if(min == (int)elementPairs[i].second.size()) {
1424 minSizeElementPairs.push_back(elementPairs[i]);
1425 }
1426 }
1427
1428 // Check if the element separated different partitions
1429 if(minSizeElementPairs.size() == elementPairs.size()) {
1430 bool isEqual = true;
1431 for(std::size_t i = 1; i < minSizeElementPairs.size(); i++) {
1432 if(minSizeElementPairs[i].second != minSizeElementPairs[0].second) {
1433 isEqual = false;
1434 break;
1435 }
1436 }
1437 if(isEqual) return nullptr;
1438 }
1439
1440 while(minSizeElementPairs.size() > 1) {
1441 min = std::numeric_limits<int>::max();
1442 for(std::size_t i = 0; i < minSizeElementPairs.size(); i++) {
1443 // The partition vector is sorted thus we can check only the first element
1444 if(minSizeElementPairs[i].second.size() == 0)
1445 return minSizeElementPairs[0].first;
1446 if(min > minSizeElementPairs[i].second[0]) {
1447 min = minSizeElementPairs[i].second[0];
1448 }
1449 }
1450
1451 for(std::size_t i = 0; i < minSizeElementPairs.size(); i++) {
1452 if(min == minSizeElementPairs[i].second[0]) {
1453 minSizeElementPairs[i].second.erase(
1454 minSizeElementPairs[i].second.begin());
1455 minSizeElementPairsTmp.push_back(minSizeElementPairs[i]);
1456 }
1457 }
1458 minSizeElementPairs.clear();
1459 minSizeElementPairs = std::move(minSizeElementPairsTmp);
1460 minSizeElementPairsTmp.clear();
1461 }
1462
1463 return minSizeElementPairs[0].first;
1464 }
1465
getPartitionInVector(std::vector<int> & partitions,const std::vector<std::pair<MElement *,std::vector<int>>> & boundaryPair)1466 static void getPartitionInVector(
1467 std::vector<int> &partitions,
1468 const std::vector<std::pair<MElement *, std::vector<int> > > &boundaryPair)
1469 {
1470 for(std::size_t i = 0; i < boundaryPair.size(); i++) {
1471 for(std::size_t j = 0; j < boundaryPair[i].second.size(); j++) {
1472 if(std::find(partitions.begin(), partitions.end(),
1473 boundaryPair[i].second[j]) == partitions.end()) {
1474 partitions.push_back(boundaryPair[i].second[j]);
1475 }
1476 }
1477 }
1478
1479 std::sort(partitions.begin(), partitions.end());
1480 }
1481
1482 template <class PART_ENTITY, class LESS_PART_ENTITY>
createPartitionEntity(std::pair<typename std::multimap<PART_ENTITY *,GEntity *,LESS_PART_ENTITY>::iterator,typename std::multimap<PART_ENTITY *,GEntity *,LESS_PART_ENTITY>::iterator> & ret,GModel * model,int & numEntity,const std::vector<int> & partitions,GEntity * referenceEntity,PART_ENTITY ** newEntity,typename std::multimap<PART_ENTITY *,GEntity *,LESS_PART_ENTITY> & pentities)1483 static PART_ENTITY *createPartitionEntity(
1484 std::pair<typename std::multimap<PART_ENTITY *, GEntity *,
1485 LESS_PART_ENTITY>::iterator,
1486 typename std::multimap<PART_ENTITY *, GEntity *,
1487 LESS_PART_ENTITY>::iterator> &ret,
1488 GModel *model, int &numEntity, const std::vector<int> &partitions,
1489 GEntity *referenceEntity, PART_ENTITY **newEntity,
1490 typename std::multimap<PART_ENTITY *, GEntity *, LESS_PART_ENTITY> &pentities)
1491 {
1492 PART_ENTITY *ppe = nullptr;
1493 // Create the new partition entity for the mesh
1494 if(ret.first == ret.second) {
1495 // Create new entity and add them to the model
1496 ppe = new PART_ENTITY(model, ++numEntity, partitions);
1497 ppe->setParentEntity(referenceEntity->getParentEntity());
1498 pentities.insert(std::make_pair(ppe, referenceEntity));
1499 model->add(ppe);
1500 *newEntity = ppe;
1501 }
1502 else {
1503 for(auto it = ret.first; it != ret.second; ++it) {
1504 if(referenceEntity == it->second) { ppe = it->first; }
1505 }
1506 if(!ppe) {
1507 // Create new entity and add them to the model
1508 ppe = new PART_ENTITY(model, ++numEntity, partitions);
1509 ppe->setParentEntity(referenceEntity->getParentEntity());
1510 pentities.insert(std::make_pair(ppe, referenceEntity));
1511 model->add(ppe);
1512 *newEntity = ppe;
1513 }
1514 }
1515
1516 return ppe;
1517 }
1518
assignPartitionBoundary(GModel * model,MFace & me,MElement * reference,const std::vector<int> & partitions,std::multimap<partitionFace *,GEntity *,partitionFacePtrLessThan> & pfaces,hashmapelement & elementToEntity,int & numEntity)1519 static partitionFace *assignPartitionBoundary(
1520 GModel *model, MFace &me, MElement *reference,
1521 const std::vector<int> &partitions,
1522 std::multimap<partitionFace *, GEntity *, partitionFacePtrLessThan> &pfaces,
1523 hashmapelement &elementToEntity, int &numEntity)
1524 {
1525 partitionFace *newEntity = nullptr;
1526 partitionFace pf(model, partitions);
1527 std::pair<std::multimap<partitionFace *, GEntity *,
1528 partitionFacePtrLessThan>::iterator,
1529 std::multimap<partitionFace *, GEntity *,
1530 partitionFacePtrLessThan>::iterator>
1531 ret = pfaces.equal_range(&pf);
1532
1533 partitionFace *ppf =
1534 createPartitionEntity(ret, model, numEntity, partitions,
1535 elementToEntity[reference], &newEntity, pfaces);
1536 int numFace = 0;
1537 for(int i = 0; i < reference->getNumFaces(); i++) {
1538 if(reference->getFace(i) == me) {
1539 numFace = i;
1540 break;
1541 }
1542 }
1543
1544 if(me.getNumVertices() == 3) {
1545 std::vector<MVertex *> verts;
1546 reference->getFaceVertices(numFace, verts);
1547
1548 if(verts.size() == 3) {
1549 MTriangle *element = new MTriangle(verts);
1550 ppf->addTriangle(element);
1551 }
1552 else if(verts.size() == 6) {
1553 MTriangle6 *element = new MTriangle6(verts);
1554 ppf->addTriangle(element);
1555 }
1556 else {
1557 MTriangleN *element =
1558 new MTriangleN(verts, verts.back()->getPolynomialOrder());
1559 ppf->addTriangle(element);
1560 }
1561 }
1562 else if(me.getNumVertices() == 4) {
1563 std::vector<MVertex *> verts;
1564 reference->getFaceVertices(numFace, verts);
1565
1566 if(verts.size() == 4) {
1567 MQuadrangle *element = new MQuadrangle(verts);
1568 ppf->addQuadrangle(element);
1569 }
1570 else if(verts.size() == 8) {
1571 MQuadrangle8 *element = new MQuadrangle8(verts);
1572 ppf->addQuadrangle(element);
1573 }
1574 else if(verts.size() == 9) {
1575 MQuadrangle9 *element = new MQuadrangle9(verts);
1576 ppf->addQuadrangle(element);
1577 }
1578 else {
1579 MQuadrangleN *element =
1580 new MQuadrangleN(verts, verts.back()->getPolynomialOrder());
1581 ppf->addQuadrangle(element);
1582 }
1583 }
1584
1585 return newEntity;
1586 }
1587
assignPartitionBoundary(GModel * model,MEdge & me,MElement * reference,const std::vector<int> & partitions,std::multimap<partitionEdge *,GEntity *,partitionEdgePtrLessThan> & pedges,hashmapelement & elementToEntity,int & numEntity)1588 static partitionEdge *assignPartitionBoundary(
1589 GModel *model, MEdge &me, MElement *reference,
1590 const std::vector<int> &partitions,
1591 std::multimap<partitionEdge *, GEntity *, partitionEdgePtrLessThan> &pedges,
1592 hashmapelement &elementToEntity, int &numEntity)
1593 {
1594 partitionEdge *newEntity = nullptr;
1595 partitionEdge pe(model, partitions);
1596 std::pair<std::multimap<partitionEdge *, GEntity *,
1597 partitionEdgePtrLessThan>::iterator,
1598 std::multimap<partitionEdge *, GEntity *,
1599 partitionEdgePtrLessThan>::iterator>
1600 ret = pedges.equal_range(&pe);
1601
1602 partitionEdge *ppe =
1603 createPartitionEntity(ret, model, numEntity, partitions,
1604 elementToEntity[reference], &newEntity, pedges);
1605
1606 int numEdge = 0;
1607 for(int i = 0; i < reference->getNumEdges(); i++) {
1608 if(reference->getEdge(i) == me) {
1609 numEdge = i;
1610 break;
1611 }
1612 }
1613
1614 if(me.getNumVertices() == 2) {
1615 std::vector<MVertex *> verts;
1616 reference->getEdgeVertices(numEdge, verts);
1617
1618 if(verts.size() == 2) {
1619 MLine *element = new MLine(verts);
1620 ppe->addLine(element);
1621 }
1622 else if(verts.size() == 3) {
1623 MLine3 *element = new MLine3(verts);
1624 ppe->addLine(element);
1625 }
1626 else {
1627 MLineN *element = new MLineN(verts);
1628 ppe->addLine(element);
1629 }
1630 }
1631
1632 return newEntity;
1633 }
1634
assignPartitionBoundary(GModel * model,MVertex * ve,MElement * reference,const std::vector<int> & partitions,std::multimap<partitionVertex *,GEntity *,partitionVertexPtrLessThan> & pvertices,hashmapelement & elementToEntity,int & numEntity)1635 static partitionVertex *assignPartitionBoundary(
1636 GModel *model, MVertex *ve, MElement *reference,
1637 const std::vector<int> &partitions,
1638 std::multimap<partitionVertex *, GEntity *, partitionVertexPtrLessThan>
1639 &pvertices,
1640 hashmapelement &elementToEntity, int &numEntity)
1641 {
1642 partitionVertex *newEntity = nullptr;
1643 partitionVertex pv(model, partitions);
1644 std::pair<std::multimap<partitionVertex *, GEntity *,
1645 partitionVertexPtrLessThan>::iterator,
1646 std::multimap<partitionVertex *, GEntity *,
1647 partitionVertexPtrLessThan>::iterator>
1648 ret = pvertices.equal_range(&pv);
1649
1650 partitionVertex *ppv =
1651 createPartitionEntity(ret, model, numEntity, partitions,
1652 elementToEntity[reference], &newEntity, pvertices);
1653
1654 ppv->addPoint(new MPoint(ve));
1655
1656 return newEntity;
1657 }
1658
computeOrientation(MElement * reference,MElement * element)1659 static int computeOrientation(MElement *reference, MElement *element)
1660 {
1661 if(element->getDim() == 2) {
1662 std::vector<MVertex *> vertices;
1663 element->getVertices(vertices);
1664 MFace face = element->getFace(0);
1665 for(int i = 0; i < reference->getNumFaces(); i++) {
1666 if(reference->getFace(i) == face) {
1667 std::vector<MVertex *> referenceVertices;
1668 reference->getFaceVertices(i, referenceVertices);
1669
1670 if(referenceVertices == vertices)
1671 return 1;
1672 else
1673 return -1;
1674 }
1675 }
1676 }
1677 else if(element->getDim() == 1) {
1678 std::vector<MVertex *> vertices;
1679 element->getVertices(vertices);
1680 MEdge face = element->getEdge(0);
1681 for(int i = 0; i < reference->getNumEdges(); i++) {
1682 if(reference->getEdge(i) == face) {
1683 std::vector<MVertex *> referenceVertices;
1684 reference->getEdgeVertices(i, referenceVertices);
1685
1686 if(referenceVertices == vertices)
1687 return 1;
1688 else
1689 return -1;
1690 }
1691 }
1692 }
1693 else if(element->getDim() == 0) {
1694 std::vector<MVertex *> vertices;
1695 element->getVertices(vertices);
1696
1697 std::vector<MVertex *> referenceVertices;
1698 reference->getVertices(referenceVertices);
1699
1700 if(referenceVertices[0] == vertices[0]) return 1;
1701 if(referenceVertices[1] == vertices[0]) return -1;
1702 }
1703
1704 return 0;
1705 }
1706
assignBrep(GModel * model,std::map<GEntity *,MElement *,GEntityPtrFullLessThan> & boundaryEntityAndRefElement,GEntity * e)1707 static void assignBrep(GModel *model,
1708 std::map<GEntity *, MElement *, GEntityPtrFullLessThan>
1709 &boundaryEntityAndRefElement,
1710 GEntity *e)
1711 {
1712 if(e->dim() == 2) {
1713 partitionFace *entity = static_cast<partitionFace *>(e);
1714
1715 for(auto it = boundaryEntityAndRefElement.begin();
1716 it != boundaryEntityAndRefElement.end(); ++it) {
1717 static_cast<GRegion *>(it->first)->setFace(
1718 entity, computeOrientation(it->second, entity->getMeshElement(0)));
1719 entity->addRegion(static_cast<GRegion *>(it->first));
1720 }
1721 }
1722 else if(e->dim() == 1) {
1723 partitionEdge *entity = static_cast<partitionEdge *>(e);
1724
1725 for(auto it = boundaryEntityAndRefElement.begin();
1726 it != boundaryEntityAndRefElement.end(); ++it) {
1727 static_cast<GFace *>(it->first)->setEdge(
1728 entity, computeOrientation(it->second, entity->getMeshElement(0)));
1729 entity->addFace(static_cast<GFace *>(it->first));
1730 }
1731 }
1732 else if(e->dim() == 0) {
1733 partitionVertex *entity = static_cast<partitionVertex *>(e);
1734
1735 for(auto it = boundaryEntityAndRefElement.begin();
1736 it != boundaryEntityAndRefElement.end(); ++it) {
1737 static_cast<GEdge *>(it->first)->setVertex(
1738 entity, computeOrientation(it->second, entity->getMeshElement(0)));
1739 entity->addEdge(static_cast<GEdge *>(it->first));
1740 }
1741 }
1742 }
1743
assignNewEntityBRep(Graph & graph,hashmapelement & elementToEntity)1744 static void assignNewEntityBRep(Graph &graph, hashmapelement &elementToEntity)
1745 {
1746 std::set<std::pair<GEntity *, GEntity *> > brepWithoutOri;
1747 hashmapentity brep;
1748 for(std::size_t i = 0; i < graph.ne(); i++) {
1749 MElement *current = graph.element(i);
1750 for(idx_t j = graph.xadj(i); j < graph.xadj(i + 1); j++) {
1751 if(current->getDim() == graph.element(graph.adjncy(j))->getDim() + 1) {
1752 GEntity *g1 = elementToEntity[current];
1753 GEntity *g2 = elementToEntity[graph.element(graph.adjncy(j))];
1754 if(brepWithoutOri.find(std::pair<GEntity *, GEntity *>(g1, g2)) ==
1755 brepWithoutOri.end()) {
1756 const int ori =
1757 computeOrientation(current, graph.element(graph.adjncy(j)));
1758 brepWithoutOri.insert(std::make_pair(g1, g2));
1759 brep[g1].insert(std::make_pair(ori, g2));
1760 }
1761 }
1762 }
1763 }
1764
1765 for(auto it = brep.begin(); it != brep.end(); ++it) {
1766 switch(it->first->dim()) {
1767 case 3:
1768 for(auto itSet = it->second.begin(); itSet != it->second.end(); ++itSet) {
1769 static_cast<GRegion *>(it->first)->setFace(
1770 static_cast<GFace *>(itSet->second), itSet->first);
1771 static_cast<GFace *>(itSet->second)
1772 ->addRegion(static_cast<GRegion *>(it->first));
1773 }
1774 break;
1775 case 2:
1776 for(auto itSet = it->second.begin(); itSet != it->second.end(); ++itSet) {
1777 static_cast<GFace *>(it->first)->setEdge(
1778 static_cast<GEdge *>(itSet->second), itSet->first);
1779 static_cast<GEdge *>(itSet->second)
1780 ->addFace(static_cast<GFace *>(it->first));
1781 }
1782 break;
1783 case 1:
1784 for(auto itSet = it->second.begin(); itSet != it->second.end(); ++itSet) {
1785 static_cast<GEdge *>(it->first)->setVertex(
1786 static_cast<GVertex *>(itSet->second), itSet->first);
1787 static_cast<GVertex *>(itSet->second)
1788 ->addEdge(static_cast<GEdge *>(it->first));
1789 }
1790 break;
1791 default: break;
1792 }
1793 }
1794 }
1795
1796 // Create the new entities between each partitions (sigma and bndSigma).
createPartitionTopology(GModel * model,const std::vector<std::set<MElement *,MElementPtrLessThan>> & boundaryElements,Graph & meshGraph)1797 static void createPartitionTopology(
1798 GModel *model,
1799 const std::vector<std::set<MElement *, MElementPtrLessThan> >
1800 &boundaryElements,
1801 Graph &meshGraph)
1802 {
1803 int meshDim = model->getMeshDim();
1804 hashmapelement elementToEntity;
1805 fillElementToEntity(model, elementToEntity, -1);
1806 assignNewEntityBRep(meshGraph, elementToEntity);
1807
1808 std::multimap<partitionFace *, GEntity *, partitionFacePtrLessThan> pfaces;
1809 std::multimap<partitionEdge *, GEntity *, partitionEdgePtrLessThan> pedges;
1810 std::multimap<partitionVertex *, GEntity *, partitionVertexPtrLessThan>
1811 pvertices;
1812
1813 hashmapface faceToElement;
1814 hashmapedge edgeToElement;
1815 hashmapvertex vertexToElement;
1816
1817 std::set<GRegion *, GEntityPtrLessThan> regions = model->getRegions();
1818 std::set<GFace *, GEntityPtrLessThan> faces = model->getFaces();
1819 std::set<GEdge *, GEntityPtrLessThan> edges = model->getEdges();
1820 std::set<GVertex *, GEntityPtrLessThan> vertices = model->getVertices();
1821
1822 if(meshDim >= 3) {
1823 Msg::Info(" - Creating partition surfaces");
1824
1825 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1826 for(auto it = boundaryElements[i].begin();
1827 it != boundaryElements[i].end(); ++it) {
1828 for(int j = 0; j < (*it)->getNumFaces(); j++) {
1829 faceToElement[(*it)->getFace(j)].push_back(
1830 std::make_pair(*it, std::vector<int>(1, i + 1)));
1831 }
1832 }
1833 }
1834 int numFaceEntity = model->getMaxElementaryNumber(2);
1835 for(auto it = faceToElement.begin(); it != faceToElement.end(); ++it) {
1836 MFace f = it->first;
1837
1838 std::vector<int> partitions;
1839 getPartitionInVector(partitions, it->second);
1840 if(partitions.size() < 2) continue;
1841
1842 MElement *reference = getReferenceElement(it->second);
1843 if(!reference) continue;
1844
1845 partitionFace *pf =
1846 assignPartitionBoundary(model, f, reference, partitions, pfaces,
1847 elementToEntity, numFaceEntity);
1848 if(pf) {
1849 std::map<GEntity *, MElement *, GEntityPtrFullLessThan>
1850 boundaryEntityAndRefElement;
1851 for(std::size_t i = 0; i < it->second.size(); i++)
1852 boundaryEntityAndRefElement.insert(std::make_pair(
1853 elementToEntity[it->second[i].first], it->second[i].first));
1854
1855 assignBrep(model, boundaryEntityAndRefElement, pf);
1856 }
1857 }
1858 faceToElement.clear();
1859
1860 faces = model->getFaces();
1861 divideNonConnectedEntities(model, 2, regions, faces, edges, vertices);
1862 elementToEntity.clear();
1863 fillElementToEntity(model, elementToEntity, 2);
1864 }
1865
1866 if(meshDim >= 2) {
1867 Msg::Info(" - Creating partition curves");
1868
1869 if(meshDim == 2) {
1870 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1871 for(auto it = boundaryElements[i].begin();
1872 it != boundaryElements[i].end(); ++it) {
1873 for(int j = 0; j < (*it)->getNumEdges(); j++) {
1874 edgeToElement[(*it)->getEdge(j)].push_back(
1875 std::make_pair(*it, std::vector<int>(1, i + 1)));
1876 }
1877 }
1878 }
1879 }
1880 else {
1881 Graph subGraph(model);
1882 makeGraph(model, subGraph, 2);
1883 subGraph.createDualGraph(false);
1884 std::vector<idx_t> part(subGraph.ne());
1885 int partIndex = 0;
1886
1887 std::map<idx_t, std::vector<int> > mapOfPartitions;
1888 idx_t mapOfPartitionsTag = 0;
1889 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
1890 if((*it)->geomType() == GEntity::PartitionSurface) {
1891 std::vector<int> partitions =
1892 static_cast<partitionFace *>(*it)->getPartitions();
1893 mapOfPartitions.insert(std::make_pair(mapOfPartitionsTag, partitions));
1894 // Must absolutely be in the same order as in the makeGraph function
1895 for(auto itElm = (*it)->triangles.begin();
1896 itElm != (*it)->triangles.end(); ++itElm)
1897 part[partIndex++] = mapOfPartitionsTag;
1898 for(auto itElm = (*it)->quadrangles.begin();
1899 itElm != (*it)->quadrangles.end(); ++itElm)
1900 part[partIndex++] = mapOfPartitionsTag;
1901 mapOfPartitionsTag++;
1902 }
1903 }
1904 subGraph.partition(part);
1905
1906 std::vector<std::set<MElement *, MElementPtrLessThan> >
1907 subBoundaryElements = subGraph.getBoundaryElements(mapOfPartitionsTag);
1908
1909 for(idx_t i = 0; i < mapOfPartitionsTag; i++) {
1910 for(auto it = subBoundaryElements[i].begin();
1911 it != subBoundaryElements[i].end(); ++it) {
1912 for(int j = 0; j < (*it)->getNumEdges(); j++) {
1913 edgeToElement[(*it)->getEdge(j)].push_back(
1914 std::make_pair(*it, mapOfPartitions[i]));
1915 }
1916 }
1917 }
1918 }
1919
1920 int numEdgeEntity = model->getMaxElementaryNumber(1);
1921 for(auto it = edgeToElement.begin(); it != edgeToElement.end(); ++it) {
1922 MEdge e = it->first;
1923
1924 std::vector<int> partitions;
1925 getPartitionInVector(partitions, it->second);
1926 if(partitions.size() < 2) continue;
1927
1928 MElement *reference = getReferenceElement(it->second);
1929 if(!reference) continue;
1930
1931 partitionEdge *pe =
1932 assignPartitionBoundary(model, e, reference, partitions, pedges,
1933 elementToEntity, numEdgeEntity);
1934 if(pe) {
1935 std::map<GEntity *, MElement *, GEntityPtrFullLessThan>
1936 boundaryEntityAndRefElement;
1937 for(std::size_t i = 0; i < it->second.size(); i++) {
1938 boundaryEntityAndRefElement.insert(std::make_pair(
1939 elementToEntity[it->second[i].first], it->second[i].first));
1940 }
1941
1942 assignBrep(model, boundaryEntityAndRefElement, pe);
1943 }
1944 }
1945 edgeToElement.clear();
1946
1947 edges = model->getEdges();
1948 divideNonConnectedEntities(model, 1, regions, faces, edges, vertices);
1949 elementToEntity.clear();
1950 fillElementToEntity(model, elementToEntity, 1);
1951 }
1952
1953 if(meshDim >= 1) {
1954 Msg::Info(" - Creating partition points");
1955 if(meshDim == 1) {
1956 for(std::size_t i = 0; i < model->getNumPartitions(); i++) {
1957 for(auto it = boundaryElements[i].begin();
1958 it != boundaryElements[i].end(); ++it) {
1959 for(std::size_t j = 0; j < (*it)->getNumPrimaryVertices(); j++) {
1960 vertexToElement[(*it)->getVertex(j)].push_back(
1961 std::make_pair(*it, std::vector<int>(1, i + 1)));
1962 }
1963 }
1964 }
1965 }
1966 else {
1967 Graph subGraph(model);
1968 makeGraph(model, subGraph, 1);
1969 subGraph.createDualGraph(false);
1970 std::vector<idx_t> part(subGraph.ne());
1971 int partIndex = 0;
1972
1973 std::map<idx_t, std::vector<int> > mapOfPartitions;
1974 idx_t mapOfPartitionsTag = 0;
1975 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
1976 if((*it)->geomType() == GEntity::PartitionCurve) {
1977 std::vector<int> partitions =
1978 static_cast<partitionEdge *>(*it)->getPartitions();
1979 mapOfPartitions.insert(std::make_pair(mapOfPartitionsTag, partitions));
1980 // Must absolutely be in the same order as in the makeGraph function
1981 for(auto itElm = (*it)->lines.begin(); itElm != (*it)->lines.end();
1982 ++itElm)
1983 part[partIndex++] = mapOfPartitionsTag;
1984 mapOfPartitionsTag++;
1985 }
1986 }
1987 subGraph.partition(part);
1988
1989 std::vector<std::set<MElement *, MElementPtrLessThan> >
1990 subBoundaryElements = subGraph.getBoundaryElements(mapOfPartitionsTag);
1991
1992 for(idx_t i = 0; i < mapOfPartitionsTag; i++) {
1993 for(auto it = subBoundaryElements[i].begin();
1994 it != subBoundaryElements[i].end(); ++it) {
1995 for(std::size_t j = 0; j < (*it)->getNumPrimaryVertices(); j++) {
1996 vertexToElement[(*it)->getVertex(j)].push_back(
1997 std::make_pair(*it, mapOfPartitions[i]));
1998 }
1999 }
2000 }
2001 }
2002 int numVertexEntity = model->getMaxElementaryNumber(0);
2003 for(auto it = vertexToElement.begin(); it != vertexToElement.end(); ++it) {
2004 MVertex *v = it->first;
2005
2006 std::vector<int> partitions;
2007 getPartitionInVector(partitions, it->second);
2008 if(partitions.size() < 2) continue;
2009
2010 MElement *reference = getReferenceElement(it->second);
2011 if(!reference) continue;
2012
2013 partitionVertex *pv =
2014 assignPartitionBoundary(model, v, reference, partitions, pvertices,
2015 elementToEntity, numVertexEntity);
2016 if(pv) {
2017 std::map<GEntity *, MElement *, GEntityPtrFullLessThan>
2018 boundaryEntityAndRefElement;
2019 for(std::size_t i = 0; i < it->second.size(); i++)
2020 boundaryEntityAndRefElement.insert(std::make_pair(
2021 elementToEntity[it->second[i].first], it->second[i].first));
2022
2023 assignBrep(model, boundaryEntityAndRefElement, pv);
2024 }
2025 }
2026 vertexToElement.clear();
2027
2028 vertices = model->getVertices();
2029 divideNonConnectedEntities(model, 0, regions, faces, edges, vertices);
2030 }
2031 }
2032
addPhysical(GModel * model,GEntity * entity,hashmap<std::string,int> & nameToNumber,std::vector<GModel::piter> & iterators,int & numPhysical)2033 static void addPhysical(GModel *model, GEntity *entity,
2034 hashmap<std::string, int> &nameToNumber,
2035 std::vector<GModel::piter> &iterators, int &numPhysical)
2036 {
2037 GEntity *parent = entity->getParentEntity();
2038 if(parent == nullptr) return;
2039
2040 if(!CTX::instance()->mesh.partitionCreatePhysicals ||
2041 CTX::instance()->mesh.partitionOldStyleMsh2) {
2042 // reuse physicals from parent entity
2043 entity->physicals = parent->physicals;
2044 return;
2045 }
2046
2047 std::size_t numPartitions = 0;
2048 if(entity->dim() == 3) {
2049 numPartitions = static_cast<partitionRegion *>(entity)->numPartitions();
2050 }
2051 else if(entity->dim() == 2) {
2052 numPartitions = static_cast<partitionFace *>(entity)->numPartitions();
2053 }
2054 else if(entity->dim() == 1) {
2055 numPartitions = static_cast<partitionEdge *>(entity)->numPartitions();
2056 }
2057 else if(entity->dim() == 0) {
2058 numPartitions = static_cast<partitionVertex *>(entity)->numPartitions();
2059 }
2060
2061 std::vector<int> physical = parent->getPhysicalEntities();
2062 int dim = entity->dim();
2063 for(size_t phys = 0; phys < physical.size(); ++phys) {
2064 std::string name = "_part{";
2065
2066 for(std::size_t i = 0; i < numPartitions; i++) {
2067 if(i > 0) name += ",";
2068 int partition = 0;
2069 if(entity->dim() == 3) {
2070 partition = static_cast<partitionRegion *>(entity)->getPartition(i);
2071 }
2072 else if(entity->dim() == 2) {
2073 partition = static_cast<partitionFace *>(entity)->getPartition(i);
2074 }
2075 else if(entity->dim() == 1) {
2076 partition = static_cast<partitionEdge *>(entity)->getPartition(i);
2077 }
2078 else if(entity->dim() == 0) {
2079 partition = static_cast<partitionVertex *>(entity)->getPartition(i);
2080 }
2081 name += std::to_string(partition);
2082 }
2083 name += "}_physical{";
2084 name +=
2085 std::to_string(physical[phys]) + "}_dim{" + std::to_string(dim) + "}";
2086
2087 int number = 0;
2088 auto it = nameToNumber.find(name);
2089 if(it == nameToNumber.end()) {
2090 number = ++numPhysical;
2091 iterators[entity->dim()] = model->setPhysicalName(
2092 iterators[entity->dim()], name, entity->dim(), number);
2093 nameToNumber.insert(std::make_pair(name, number));
2094 }
2095 else {
2096 number = it->second;
2097 }
2098 entity->addPhysicalEntity(number);
2099 }
2100
2101 if(physical.size() == 0) {
2102 std::string name = "_part{";
2103
2104 for(std::size_t i = 0; i < numPartitions; i++) {
2105 if(i > 0) name += ",";
2106 int partition = 0;
2107 if(entity->dim() == 3) {
2108 partition = static_cast<partitionRegion *>(entity)->getPartition(i);
2109 }
2110 else if(entity->dim() == 2) {
2111 partition = static_cast<partitionFace *>(entity)->getPartition(i);
2112 }
2113 else if(entity->dim() == 1) {
2114 partition = static_cast<partitionEdge *>(entity)->getPartition(i);
2115 }
2116 else if(entity->dim() == 0) {
2117 partition = static_cast<partitionVertex *>(entity)->getPartition(i);
2118 }
2119 name += std::to_string(partition);
2120 }
2121 name += "}_";
2122 name += "physical{0}_dim{" + std::to_string(dim) + "}";
2123
2124 int number = 0;
2125 auto it = nameToNumber.find(name);
2126 if(it == nameToNumber.end()) {
2127 number = ++numPhysical;
2128 iterators[entity->dim()] = model->setPhysicalName(
2129 iterators[entity->dim()], name, entity->dim(), number);
2130 nameToNumber.insert(std::make_pair(name, number));
2131 }
2132 else {
2133 number = it->second;
2134 }
2135 entity->addPhysicalEntity(number);
2136 }
2137 }
2138
2139 // Assign physical group information
assignPhysicals(GModel * model)2140 static void assignPhysicals(GModel *model)
2141 {
2142 hashmap<std::string, int> nameToNumber;
2143 std::vector<GModel::piter> iterators;
2144 model->getInnerPhysicalNamesIterators(iterators);
2145 int numPhysical = model->getMaxPhysicalNumber(-1);
2146 // Loop over volumes
2147 for(auto it = model->firstRegion(); it != model->lastRegion(); ++it) {
2148 if((*it)->geomType() == GEntity::PartitionVolume) {
2149 addPhysical(model, *it, nameToNumber, iterators, numPhysical);
2150 }
2151 }
2152
2153 // Loop over surfaces
2154 for(auto it = model->firstFace(); it != model->lastFace(); ++it) {
2155 if((*it)->geomType() == GEntity::PartitionSurface) {
2156 addPhysical(model, *it, nameToNumber, iterators, numPhysical);
2157 }
2158 }
2159
2160 // Loop over curves
2161 for(auto it = model->firstEdge(); it != model->lastEdge(); ++it) {
2162 if((*it)->geomType() == GEntity::PartitionCurve) {
2163 addPhysical(model, *it, nameToNumber, iterators, numPhysical);
2164 }
2165 }
2166
2167 // Loop over points
2168 for(auto it = model->firstVertex(); it != model->lastVertex(); ++it) {
2169 if((*it)->geomType() == GEntity::PartitionPoint) {
2170 addPhysical(model, *it, nameToNumber, iterators, numPhysical);
2171 }
2172 }
2173 }
2174
cmp_hedges(const std::pair<MEdge,size_t> & he0,const std::pair<MEdge,size_t> & he1)2175 bool cmp_hedges(const std::pair<MEdge, size_t> &he0,
2176 const std::pair<MEdge, size_t> &he1)
2177 {
2178 MEdgeLessThan cmp;
2179 return cmp(he0.first, he1.first);
2180 }
2181
PartitionFaceMinEdgeLength(GFace * gf,int np,double tol)2182 int PartitionFaceMinEdgeLength(GFace *gf, int np, double tol)
2183 {
2184 std::vector<std::pair<MEdge, size_t> > halfEdges;
2185 halfEdges.reserve(gf->triangles.size() * 3);
2186 for(size_t i = 0; i < gf->triangles.size(); ++i) {
2187 for(size_t j = 0; j < 3; ++j) {
2188 halfEdges.push_back(std::make_pair(gf->triangles[i]->getEdge(j), i));
2189 }
2190 }
2191 std::sort(halfEdges.begin(), halfEdges.end(), cmp_hedges);
2192 std::vector<idx_t> neighbors(gf->triangles.size() * 3, -1);
2193 std::vector<double> neighborsWeight(gf->triangles.size() * 3, -1);
2194 MEdgeEqual eq;
2195 double minEdgeLength = std::numeric_limits<double>::max();
2196 for(size_t i = 0; i + 1 < halfEdges.size(); ++i) {
2197 if(eq(halfEdges[i].first, halfEdges[i + 1].first)) {
2198 size_t t0 = halfEdges[i].second;
2199 size_t t1 = halfEdges[i + 1].second;
2200 double l = halfEdges[i].first.length();
2201 minEdgeLength = std::min(l, minEdgeLength);
2202 for(int j = 0; j < 3; ++j) {
2203 if(neighbors[t0 * 3 + j] == -1) {
2204 neighbors[t0 * 3 + j] = t1;
2205 neighborsWeight[t0 * 3 + j] = l;
2206 break;
2207 }
2208 }
2209 for(int j = 0; j < 3; ++j) {
2210 if(neighbors[t1 * 3 + j] == -1) {
2211 neighbors[t1 * 3 + j] = t0;
2212 neighborsWeight[t1 * 3 + j] = l;
2213 break;
2214 }
2215 }
2216 i++;
2217 }
2218 }
2219 std::vector<idx_t> adjncy;
2220 std::vector<idx_t> xadjncy;
2221 std::vector<idx_t> adjncyw;
2222 xadjncy.push_back(0);
2223 for(size_t i = 0; i < gf->triangles.size(); ++i) {
2224 for(size_t j = 0; j < 3; ++j) {
2225 if(neighbors[i * 3 + j] == -1) break;
2226 adjncy.push_back(neighbors[i * 3 + j]);
2227 adjncyw.push_back(
2228 (idx_t)(neighborsWeight[i * 3 + j] / minEdgeLength * 10));
2229 }
2230 xadjncy.push_back(adjncy.size());
2231 }
2232 idx_t nvtxs = gf->triangles.size(), ncon = 1, nparts = np, objval;
2233 std::vector<idx_t> epart(gf->triangles.size());
2234 real_t ubvec = tol;
2235 METIS_PartGraphKway(&nvtxs, &ncon, &xadjncy[0], &adjncy[0], nullptr, nullptr,
2236 &adjncyw[0], &nparts, nullptr, &ubvec, nullptr, &objval,
2237 &epart[0]);
2238 for(size_t i = 0; i < gf->triangles.size(); ++i) {
2239 gf->triangles[i]->setPartition(epart[i]);
2240 }
2241 return 0;
2242 }
2243
2244 // Partition a mesh into n parts. Returns: 0 = success, 1 = error
2245
PartitionMesh(GModel * model,int numPart)2246 int PartitionMesh(GModel *model, int numPart)
2247 {
2248 if(numPart <= 0) return 0;
2249
2250 Msg::StatusBar(true, "Partitioning mesh...");
2251 double t1 = Cpu(), w1 = TimeOfDay();
2252
2253 Graph graph(model);
2254 if(makeGraph(model, graph, -1)) return 1;
2255 graph.nparts(numPart);
2256 if(partitionGraph(graph, true)) return 1;
2257
2258 std::vector<std::size_t> elmCount[TYPE_MAX_NUM + 1];
2259 for(int i = 0; i < TYPE_MAX_NUM + 1; i++) { elmCount[i].resize(numPart, 0); }
2260
2261 // Assign partitions to elements
2262 hashmapelementpart elmToPartition;
2263 for(std::size_t i = 0; i < graph.ne(); i++) {
2264 if(graph.element(i)) {
2265 if(graph.nparts() > 1) {
2266 elmToPartition.insert(std::make_pair(
2267 graph.element(i), graph.partition(i) + 1));
2268 elmCount[graph.element(i)->getType()][graph.partition(i)]++;
2269 // Should be removed
2270 graph.element(i)->setPartition(graph.partition(i) + 1);
2271 }
2272 else {
2273 elmToPartition.insert(std::make_pair(graph.element(i), 1));
2274 // Should be removed
2275 graph.element(i)->setPartition(1);
2276 }
2277 }
2278 }
2279 model->setNumPartitions(graph.nparts());
2280
2281 createNewEntities(model, elmToPartition);
2282 elmToPartition.clear();
2283
2284 double t2 = Cpu(), w2 = TimeOfDay();
2285 Msg::StatusBar(true, "Done partitioning mesh (Wall %gs, CPU %gs)", w2 - w1,
2286 t2 - t1);
2287
2288 for(std::size_t i = 0; i < TYPE_MAX_NUM + 1; i++) {
2289 std::vector<std::size_t> &count = elmCount[i];
2290 std::size_t minCount = std::numeric_limits<std::size_t>::max();
2291 std::size_t maxCount = 0;
2292 std::size_t totCount = 0;
2293 for(std::size_t j = 0; j < count.size(); j++) {
2294 minCount = std::min(count[j], minCount);
2295 maxCount = std::max(count[j], maxCount);
2296 totCount += count[j];
2297 }
2298 if(totCount > 0) {
2299 Msg::Info(" - Repartition of %d %s: %lu(min) %lu(max) %g(avg)", totCount,
2300 ElementType::nameOfParentType(i, totCount > 1).c_str(),
2301 minCount, maxCount, totCount / (double)numPart);
2302 }
2303 }
2304
2305 if(CTX::instance()->mesh.partitionCreateTopology) {
2306 Msg::StatusBar(true, "Creating partition topology...");
2307 std::vector<std::set<MElement *, MElementPtrLessThan> > boundaryElements =
2308 graph.getBoundaryElements();
2309 createPartitionTopology(model, boundaryElements, graph);
2310 boundaryElements.clear();
2311 double t3 = Cpu(), w3 = TimeOfDay();
2312 Msg::StatusBar(true, "Done creating partition topology (Wall %gs, CPU %gs)",
2313 w3 - w2, t3 - t2);
2314 }
2315
2316 assignPhysicals(model);
2317 assignMeshVertices(model);
2318
2319 if(CTX::instance()->mesh.partitionCreateGhostCells) {
2320 double t4 = Cpu(), w4 = TimeOfDay();
2321 Msg::StatusBar(true, "Creating ghost cells...");
2322 graph.clearDualGraph();
2323 graph.createDualGraph(true);
2324 graph.assignGhostCells();
2325 double t5 = Cpu(), w5 = TimeOfDay();
2326 Msg::StatusBar(true, "Done creating ghost cells (Wall %gs, CPU %gs)",
2327 w5 - w4, t5 - t4);
2328 }
2329
2330 return 0;
2331 }
2332
2333 template <class ITERATOR, class PART_ENTITY>
assignToParent(std::set<MVertex * > & verts,PART_ENTITY * entity,ITERATOR it_beg,ITERATOR it_end)2334 static void assignToParent(std::set<MVertex *> &verts, PART_ENTITY *entity,
2335 ITERATOR it_beg, ITERATOR it_end)
2336 {
2337 for(ITERATOR it = it_beg; it != it_end; ++it) {
2338 entity->getParentEntity()->addElement((*it)->getType(), *it);
2339 (*it)->setPartition(0);
2340
2341 for(std::size_t i = 0; i < (*it)->getNumVertices(); i++) {
2342 if(verts.find((*it)->getVertex(i)) == verts.end()) {
2343 (*it)->getVertex(i)->setEntity(entity->getParentEntity());
2344 entity->getParentEntity()->addMeshVertex((*it)->getVertex(i));
2345 verts.insert((*it)->getVertex(i));
2346 }
2347 }
2348 }
2349 }
2350
2351 // Un-partition a mesh and return to the initial mesh geomerty. Returns: 0 =
2352 // success, 1 = error.
UnpartitionMesh(GModel * model)2353 int UnpartitionMesh(GModel *model)
2354 {
2355 // make a copy so we can iterate safely (we will remove some entities)
2356 std::set<GRegion *, GEntityPtrLessThan> regions = model->getRegions();
2357 std::set<GFace *, GEntityPtrLessThan> faces = model->getFaces();
2358 std::set<GEdge *, GEntityPtrLessThan> edges = model->getEdges();
2359 std::set<GVertex *, GEntityPtrLessThan> vertices = model->getVertices();
2360
2361 std::set<MVertex *> verts;
2362
2363 // Loop over points
2364 for(auto it = vertices.begin(); it != vertices.end(); ++it) {
2365 GVertex *vertex = *it;
2366
2367 if(vertex->geomType() == GEntity::PartitionPoint) {
2368 partitionVertex *pvertex = static_cast<partitionVertex *>(vertex);
2369 if(pvertex->getParentEntity() && pvertex->getParentEntity()->dim() == 0) {
2370 assignToParent(verts, pvertex, pvertex->points.begin(),
2371 pvertex->points.end());
2372 }
2373 else {
2374 for(std::size_t j = 0; j < pvertex->points.size(); j++)
2375 delete pvertex->points[j];
2376 }
2377 pvertex->points.clear();
2378 pvertex->mesh_vertices.clear();
2379
2380 model->remove(pvertex);
2381 delete pvertex;
2382 }
2383 }
2384
2385 // Loop over curves
2386 for(auto it = edges.begin(); it != edges.end(); ++it) {
2387 GEdge *edge = *it;
2388 if(edge->geomType() == GEntity::PartitionCurve) {
2389 partitionEdge *pedge = static_cast<partitionEdge *>(edge);
2390 if(pedge->getParentEntity() && pedge->getParentEntity()->dim() == 1) {
2391 assignToParent(verts, pedge, pedge->lines.begin(), pedge->lines.end());
2392 }
2393 else {
2394 for(std::size_t j = 0; j < pedge->lines.size(); j++)
2395 delete pedge->lines[j];
2396 }
2397 pedge->lines.clear();
2398 pedge->mesh_vertices.clear();
2399 pedge->setBeginVertex(nullptr);
2400 pedge->setEndVertex(nullptr);
2401
2402 model->remove(pedge);
2403 delete pedge;
2404 }
2405 else if(edge->geomType() == GEntity::GhostCurve) {
2406 model->remove(edge);
2407 delete edge;
2408 }
2409 }
2410
2411 // Loop over surfaces
2412 for(auto it = faces.begin(); it != faces.end(); ++it) {
2413 GFace *face = *it;
2414
2415 if(face->geomType() == GEntity::PartitionSurface) {
2416 partitionFace *pface = static_cast<partitionFace *>(face);
2417 if(pface->getParentEntity() && pface->getParentEntity()->dim() == 2) {
2418 assignToParent(verts, pface, pface->triangles.begin(),
2419 pface->triangles.end());
2420 assignToParent(verts, pface, pface->quadrangles.begin(),
2421 pface->quadrangles.end());
2422 }
2423 else {
2424 for(std::size_t j = 0; j < pface->triangles.size(); j++)
2425 delete pface->triangles[j];
2426 for(std::size_t j = 0; j < pface->quadrangles.size(); j++)
2427 delete pface->quadrangles[j];
2428 }
2429 pface->triangles.clear();
2430 pface->quadrangles.clear();
2431 pface->mesh_vertices.clear();
2432 pface->set(std::vector<GEdge *>());
2433 pface->setOrientations(std::vector<int>());
2434
2435 model->remove(pface);
2436 delete pface;
2437 }
2438 else if(face->geomType() == GEntity::GhostSurface) {
2439 model->remove(face);
2440 delete face;
2441 }
2442 }
2443
2444 // Loop over volumes
2445 for(auto it = regions.begin(); it != regions.end(); ++it) {
2446 GRegion *region = *it;
2447
2448 if(region->geomType() == GEntity::PartitionVolume) {
2449 partitionRegion *pregion = static_cast<partitionRegion *>(region);
2450 if(pregion->getParentEntity() && pregion->getParentEntity()->dim() == 3) {
2451 assignToParent(verts, pregion, pregion->tetrahedra.begin(),
2452 pregion->tetrahedra.end());
2453 assignToParent(verts, pregion, pregion->hexahedra.begin(),
2454 pregion->hexahedra.end());
2455 assignToParent(verts, pregion, pregion->prisms.begin(),
2456 pregion->prisms.end());
2457 assignToParent(verts, pregion, pregion->pyramids.begin(),
2458 pregion->pyramids.end());
2459 assignToParent(verts, pregion, pregion->trihedra.begin(),
2460 pregion->trihedra.end());
2461 }
2462 else {
2463 for(std::size_t j = 0; j < pregion->tetrahedra.size(); j++)
2464 delete pregion->tetrahedra[j];
2465 for(std::size_t j = 0; j < pregion->hexahedra.size(); j++)
2466 delete pregion->hexahedra[j];
2467 for(std::size_t j = 0; j < pregion->prisms.size(); j++)
2468 delete pregion->prisms[j];
2469 for(std::size_t j = 0; j < pregion->pyramids.size(); j++)
2470 delete pregion->pyramids[j];
2471 for(std::size_t j = 0; j < pregion->trihedra.size(); j++)
2472 delete pregion->trihedra[j];
2473 }
2474 pregion->tetrahedra.clear();
2475 pregion->hexahedra.clear();
2476 pregion->prisms.clear();
2477 pregion->pyramids.clear();
2478 pregion->trihedra.clear();
2479 pregion->mesh_vertices.clear();
2480 pregion->set(std::vector<GFace *>());
2481 pregion->setOrientations(std::vector<int>());
2482
2483 model->remove(pregion);
2484 delete pregion;
2485 }
2486 else if(region->geomType() == GEntity::GhostVolume) {
2487 model->remove(region);
2488 delete region;
2489 }
2490 }
2491
2492 model->setNumPartitions(0);
2493
2494 std::map<std::pair<int, int>, std::string> physicalNames =
2495 model->getPhysicalNames();
2496 for(auto it = physicalNames.begin(); it != physicalNames.end(); ++it) {
2497 std::size_t found = it->second.find("_");
2498 if(found != std::string::npos) {
2499 model->removePhysicalGroup(it->first.first, it->first.second);
2500 }
2501 }
2502
2503 return 0;
2504 }
2505
2506 // Create the partition according to the element split given by elmToPartition
2507 // Returns: 0 = success, 1 = no elements found/invalid data.
PartitionUsingThisSplit(GModel * model,std::vector<std::pair<MElement *,int>> & elmToPart)2508 int PartitionUsingThisSplit(GModel *model,
2509 std::vector<std::pair<MElement *, int> > &elmToPart)
2510 {
2511 Graph graph(model);
2512 if(makeGraph(model, graph, -1)) return 1;
2513
2514 int numPart = 0;
2515 hashmapelementpart elmToPartition;
2516 std::vector<std::pair<MElement*, int> > elmGhosts;
2517 for(auto item : elmToPart) {
2518 MElement *el = item.first;
2519 int part = item.second;
2520 if(part == 0) {
2521 Msg::Error("Partition tag cannot be 0");
2522 return 1;
2523 }
2524 if(part > 0) elmToPartition[el] = part;
2525 else elmGhosts.push_back(std::make_pair(el, -part));
2526 numPart = std::max(std::abs(part), numPart);
2527 }
2528
2529 graph.createDualGraph(false);
2530 graph.nparts(numPart);
2531
2532 if(elmToPartition.size() != graph.ne()) {
2533 Msg::Error("All elements are not partitioned");
2534 return 1;
2535 }
2536
2537 std::vector<idx_t> part(graph.ne());
2538 for(std::size_t i = 0; i < graph.ne(); i++) {
2539 if(graph.element(i)) { part[i] = elmToPartition[graph.element(i)] - 1; }
2540 }
2541
2542 // Check and correct the topology
2543 for(int i = 1; i < 4; i++) {
2544 for(std::size_t j = 0; j < graph.ne(); j++) {
2545 if(graph.element(j)->getDim() == graph.dim()) continue;
2546
2547 for(idx_t k = graph.xadj(j); k < graph.xadj(j + 1); k++) {
2548 if(graph.element(j)->getDim() ==
2549 graph.element(graph.adjncy(k))->getDim() - i) {
2550 if(part[j] != part[graph.adjncy(k)]) {
2551 part[j] = part[graph.adjncy(k)];
2552 break;
2553 }
2554 }
2555 }
2556 }
2557 }
2558 graph.partition(part);
2559
2560 model->setNumPartitions(graph.nparts());
2561
2562 createNewEntities(model, elmToPartition);
2563
2564 if(CTX::instance()->mesh.partitionCreateTopology) {
2565 Msg::StatusBar(true, "Creating partition topology...");
2566 std::vector<std::set<MElement *, MElementPtrLessThan> > boundaryElements =
2567 graph.getBoundaryElements();
2568 createPartitionTopology(model, boundaryElements, graph);
2569 boundaryElements.clear();
2570 Msg::StatusBar(true, "Done creating partition topology");
2571 }
2572
2573 assignPhysicals(model);
2574 assignMeshVertices(model);
2575
2576 if (!elmGhosts.empty()) {
2577 std::sort(elmGhosts.begin(), elmGhosts.end());
2578 auto last = std::unique(elmGhosts.begin(), elmGhosts.end());
2579 elmGhosts.erase(last, elmGhosts.end());
2580 std::vector<GEntity *> ghostEntities = graph.createGhostEntities();
2581 for (auto elmGhost : elmGhosts) {
2582 MElement *el = elmGhost.first;
2583 int part = elmGhost.second;
2584 if(el->getDim() == graph.dim()) {
2585 switch(graph.dim()) {
2586 case 1:
2587 static_cast<ghostEdge *>(ghostEntities[part - 1])
2588 ->addElement(el->getType(), el, elmToPartition[el]);
2589 break;
2590 case 2:
2591 static_cast<ghostFace *>(ghostEntities[part - 1])
2592 ->addElement(el->getType(), el, elmToPartition[el]);
2593 break;
2594 case 3:
2595 static_cast<ghostRegion *>(ghostEntities[part - 1])
2596 ->addElement(el->getType(), el, elmToPartition[el]);
2597 break;
2598 default: break;
2599 }
2600 }
2601 }
2602 }
2603 else if(CTX::instance()->mesh.partitionCreateGhostCells) {
2604 graph.clearDualGraph();
2605 graph.createDualGraph(false);
2606 graph.assignGhostCells();
2607 }
2608 elmToPartition.clear();
2609
2610 return 0;
2611 }
2612
2613 // Import a mesh partitionned by a tag given to the element and create the
2614 // topology (omega, sigma, bndSigma, ...). Returns: 0 = success, 1 = no elements
2615 // found.
ConvertOldPartitioningToNewOne(GModel * model)2616 int ConvertOldPartitioningToNewOne(GModel *model)
2617 {
2618 Msg::StatusBar(true, "Converting old partitioning...");
2619
2620 std::vector<std::pair<MElement *, int> > elmToPartition;
2621 std::set<int> partitions;
2622 std::vector<GEntity *> entities;
2623 model->getEntities(entities);
2624 for(std::size_t i = 0; i < entities.size(); i++) {
2625 for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
2626 MElement *e = entities[i]->getMeshElement(j);
2627 idx_t part = e->getPartition();
2628 if(part < 0) part = -part;
2629 if(!part) part = 1;
2630 elmToPartition.push_back(std::make_pair(e, part));
2631 partitions.insert(part);
2632 }
2633 }
2634
2635 return PartitionUsingThisSplit(model, elmToPartition);
2636 }
2637
2638 #else
2639
PartitionMesh(GModel * model,int numPart)2640 int PartitionMesh(GModel *model, int numPart)
2641 {
2642 Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
2643 return 0;
2644 }
2645
UnpartitionMesh(GModel * model)2646 int UnpartitionMesh(GModel *model) { return 0; }
2647
ConvertOldPartitioningToNewOne(GModel * model)2648 int ConvertOldPartitioningToNewOne(GModel *model) { return 0; }
2649
PartitionUsingThisSplit(GModel * model,std::vector<std::pair<MElement *,int>> & elmToPartition)2650 int PartitionUsingThisSplit(
2651 GModel *model, std::vector<std::pair<MElement *, int> > &elmToPartition)
2652 {
2653 Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
2654 return 0;
2655 }
2656
PartitionFaceMinEdgeLength(GFace * gf,int np,double tol)2657 int PartitionFaceMinEdgeLength(GFace *gf, int np, double tol)
2658 {
2659 Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
2660 return 0;
2661 }
2662
2663 #endif
2664