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> &regions,
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