1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Matti Pellikka <matti.pellikka@gmail.com>.
7 
8 #include "Cell.h"
9 #include "MTriangle.h"
10 #include "MQuadrangle.h"
11 #include "MTetrahedron.h"
12 #include "MPyramid.h"
13 #include "MHexahedron.h"
14 #include "MPrism.h"
15 
operator ()(const Cell * c1,const Cell * c2) const16 bool CellPtrLessThan::operator()(const Cell *c1, const Cell *c2) const
17 {
18   // If cell complex is done use enumeration (same as vertex order)
19   if(c1->getNum() != 0) return (c1->getNum() < c2->getNum());
20 
21   // Otherwise order by vertex numbering (good heuristic for reduction)
22   if(c1->getNumSortedVertices() != c2->getNumSortedVertices()) {
23     return (c1->getNumSortedVertices() < c2->getNumSortedVertices());
24   }
25   for(int i = 0; i < c1->getNumSortedVertices(); i++) {
26     if(c1->getSortedVertex(i) < c2->getSortedVertex(i))
27       return true;
28     else if(c1->getSortedVertex(i) > c2->getSortedVertex(i))
29       return false;
30   }
31   return false;
32 }
33 
equalVertices(const std::vector<MVertex * > & v1,const std::vector<MVertex * > & v2)34 bool equalVertices(const std::vector<MVertex *> &v1,
35                    const std::vector<MVertex *> &v2)
36 {
37   if(v1.size() != v2.size()) return false;
38   for(std::size_t i = 0; i < v1.size(); i++)
39     if(v1[i]->getNum() != v2[i]->getNum()) return false;
40   return true;
41 }
42 
43 int Cell::_globalNum = 0;
44 
createCell(MElement * element,int domain)45 std::pair<Cell *, bool> Cell::createCell(MElement *element, int domain)
46 {
47   Cell *cell = new Cell();
48   cell->_dim = element->getDim();
49   cell->_domain = domain;
50   cell->_combined = false;
51   cell->_immune = false;
52   cell->_num = 0;
53 
54   for(std::size_t i = 0; i < element->getNumPrimaryVertices(); i++)
55     cell->_v.push_back(element->getVertex(i));
56 
57   return std::make_pair(cell, cell->_sortVertexIndices());
58 }
59 
createCell(Cell * parent,int i)60 std::pair<Cell *, bool> Cell::createCell(Cell *parent, int i)
61 {
62   Cell *cell = new Cell();
63   cell->_dim = parent->getDim() - 1;
64   cell->_domain = parent->getDomain();
65   cell->_combined = false;
66   cell->_immune = false;
67   cell->_num = 0;
68 
69   parent->findBdElement(i, cell->_v);
70   return std::make_pair(cell, cell->_sortVertexIndices());
71 }
72 
Cell(MElement * element,int domain)73 Cell::Cell(MElement *element, int domain)
74 {
75   _dim = element->getDim();
76   _domain = domain;
77   _combined = false;
78   _immune = false;
79   _num = 0;
80 
81   for(std::size_t i = 0; i < element->getNumPrimaryVertices(); i++)
82     _v.push_back(element->getVertex(i));
83 
84   _sortVertexIndices();
85 }
86 
Cell(Cell * parent,int i)87 Cell::Cell(Cell *parent, int i)
88 {
89   _dim = parent->getDim() - 1;
90   _domain = parent->getDomain();
91   _combined = false;
92   _immune = false;
93   _num = 0;
94 
95   parent->findBdElement(i, _v);
96   _sortVertexIndices();
97 }
98 
_sortVertexIndices()99 bool Cell::_sortVertexIndices()
100 {
101   std::map<MVertex *, int, MVertexPtrLessThan> si;
102 
103   bool noinsert = false;
104   for(std::size_t i = 0; i < _v.size(); i++)
105     noinsert = (!si.insert(std::make_pair(_v[i], i)).second || noinsert);
106 
107   if(noinsert == true) {
108     Msg::Warning("The input mesh has degenerate elements, ignored");
109     return false;
110   }
111 
112   for(auto it = si.begin(); it != si.end(); it++) _si.push_back(it->second);
113 
114   return true;
115 }
116 
getSortedVertex(int vertex) const117 inline int Cell::getSortedVertex(int vertex) const
118 {
119   return _v[(int)_si[vertex]]->getNum();
120 }
121 
findBdElement(int i,std::vector<MVertex * > & vertices) const122 void Cell::findBdElement(int i, std::vector<MVertex *> &vertices) const
123 {
124   vertices.clear();
125   switch(_dim) {
126   case 1: vertices.push_back(_v[i]); return;
127   case 2:
128     switch(getNumVertices()) {
129     case 3:
130       for(int j = 0; j < 2; j++)
131         vertices.push_back(_v[MTriangle::edges_tri(i, j)]);
132       return;
133     case 4:
134       for(int j = 0; j < 2; j++)
135         vertices.push_back(_v[MQuadrangle::edges_quad(i, j)]);
136       return;
137     default: return;
138     }
139   case 3:
140     switch(getNumVertices()) {
141     case 4:
142       for(int j = 0; j < 3; j++)
143         vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
144       return;
145     case 5:
146       if(i < 4)
147         for(int j = 0; j < 3; j++)
148           vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
149       else
150         for(int j = 0; j < 4; j++)
151           vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
152       return;
153     case 6:
154       if(i < 2)
155         for(int j = 0; j < 3; j++)
156           vertices.push_back(_v[MPrism::faces_prism(i, j)]);
157       else
158         for(int j = 0; j < 4; j++)
159           vertices.push_back(_v[MPrism::faces_prism(i, j)]);
160       return;
161     case 8:
162       for(int j = 0; j < 4; j++)
163         vertices.push_back(_v[MHexahedron::faces_hexa(i, j)]);
164       return;
165     default: return;
166     }
167   default: return;
168   }
169 }
170 
getNumBdElements() const171 int Cell::getNumBdElements() const
172 {
173   switch(_dim) {
174   case 0: return 0;
175   case 1: return 2;
176   case 2:
177     switch(getNumVertices()) {
178     case 3: return 3;
179     case 4: return 4;
180     default: return 0;
181     }
182   case 3:
183     switch(getNumVertices()) {
184     case 4: return 4;
185     case 5: return 5;
186     case 6: return 5;
187     case 8: return 6;
188     default: return 0;
189     }
190   default: return 0;
191   }
192 }
193 
findBdCellOrientation(Cell * cell,int i) const194 int Cell::findBdCellOrientation(Cell *cell, int i) const
195 {
196   /*std::vector<MVertex*> v1;
197   std::vector<MVertex*> v2;
198   this->findBdElement(i, v1);
199   cell->getMeshVertices(v2);
200 
201   int perm = 1;
202   if(equalVertices(v1, v2)) return perm;
203   while(std::next_permutation(v2.begin(), v2.end(), MVertexPtrLessThan())) {
204     perm *= -1;
205     if(equalVertices(v1, v2)) return perm;
206   }
207   cell->getMeshVertices(v2);
208   perm = 1;
209   while(std::prev_permutation(v2.begin(), v2.end(), MVertexPtrLessThan())) {
210     perm *= -1;
211     if(equalVertices(v1, v2)) return perm;
212   }
213   return 0;*/
214 
215   std::vector<MVertex *> v;
216   cell->getMeshVertices(v);
217   switch(_dim) {
218   case 0: return 0;
219   case 1:
220     if(v[0]->getNum() == _v[0]->getNum())
221       return -1;
222     else if(v[0]->getNum() == _v[1]->getNum())
223       return 1;
224     return 0;
225   case 2:
226     switch(getNumVertices()) {
227     case 3:
228       if(_v[MTriangle::edges_tri(i, 0)]->getNum() == v[0]->getNum() &&
229          _v[MTriangle::edges_tri(i, 1)]->getNum() == v[1]->getNum())
230         return 1;
231       if(_v[MTriangle::edges_tri(i, 1)]->getNum() == v[0]->getNum() &&
232          _v[MTriangle::edges_tri(i, 0)]->getNum() == v[1]->getNum())
233         return -1;
234       break;
235     case 4:
236       if(_v[MQuadrangle::edges_quad(i, 0)]->getNum() == v[0]->getNum() &&
237          _v[MQuadrangle::edges_quad(i, 1)]->getNum() == v[1]->getNum())
238         return 1;
239       if(_v[MQuadrangle::edges_quad(i, 0)]->getNum() == v[1]->getNum() &&
240          _v[MQuadrangle::edges_quad(i, 1)]->getNum() == v[0]->getNum())
241         return -1;
242       break;
243     default: return 0;
244     }
245   case 3:
246     switch(getNumVertices()) {
247     case 4:
248       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[0]->getNum() &&
249          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
250          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[2]->getNum())
251         return 1;
252       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[1]->getNum() &&
253          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[2]->getNum() &&
254          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
255         return 1;
256       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[2]->getNum() &&
257          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[0]->getNum() &&
258          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[1]->getNum())
259         return 1;
260       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[0]->getNum() &&
261          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[2]->getNum() &&
262          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[1]->getNum())
263         return -1;
264       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[1]->getNum() &&
265          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[0]->getNum() &&
266          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[2]->getNum())
267         return -1;
268       if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[2]->getNum() &&
269          _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
270          _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
271         return -1;
272       break;
273     case 5:
274       if(i < 4) {
275         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
276            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
277            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
278           return 1;
279         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
280            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
281            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
282           return 1;
283         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
284            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
285            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
286           return 1;
287         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
288            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
289            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
290           return -1;
291         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
292            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
293            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
294           return -1;
295         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
296            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
297            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
298           return -1;
299       }
300       else {
301         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
302            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
303            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
304            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
305           return 1;
306         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
307            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
308            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
309            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
310           return 1;
311         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
312            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
313            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
314            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
315           return 1;
316         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
317            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
318            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
319            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
320           return 1;
321         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
322            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
323            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
324            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
325           return -1;
326         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
327            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
328            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
329            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
330           return -1;
331         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
332            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
333            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
334            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
335           return -1;
336         if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
337            _v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
338            _v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
339            _v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
340           return -1;
341       }
342       return 0;
343       break;
344     case 6:
345       if(i < 2) {
346         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
347            _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
348            _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
349           return 1;
350         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
351            _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
352            _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
353           return 1;
354         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
355            _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
356            _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
357           return 1;
358         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
359            _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
360            _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
361           return -1;
362         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
363            _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
364            _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
365           return -1;
366         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
367            _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
368            _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
369           return -1;
370       }
371       else {
372         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
373            _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
374            _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
375            _v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
376           return 1;
377         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
378            _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
379            _v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
380            _v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
381           return 1;
382         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
383            _v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
384            _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
385            _v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
386           return 1;
387         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
388            _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
389            _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
390            _v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
391           return 1;
392         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
393            _v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
394            _v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
395            _v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
396           return -1;
397         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
398            _v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
399            _v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
400            _v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
401           return -1;
402         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
403            _v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
404            _v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
405            _v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
406           return -1;
407         if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
408            _v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
409            _v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
410            _v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
411           return -1;
412       }
413       break;
414     case 8:
415       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
416          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
417          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[2]->getNum() &&
418          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[3]->getNum())
419         return 1;
420       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[1]->getNum() &&
421          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[2]->getNum() &&
422          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[3]->getNum() &&
423          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[0]->getNum())
424         return 1;
425       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[2]->getNum() &&
426          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[3]->getNum() &&
427          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[0]->getNum() &&
428          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[1]->getNum())
429         return 1;
430       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[3]->getNum() &&
431          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[0]->getNum() &&
432          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[1]->getNum() &&
433          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[2]->getNum())
434         return 1;
435       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
436          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[3]->getNum() &&
437          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[2]->getNum() &&
438          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[1]->getNum())
439         return -1;
440       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[3]->getNum() &&
441          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[2]->getNum() &&
442          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[1]->getNum() &&
443          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[0]->getNum())
444         return -1;
445       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[2]->getNum() &&
446          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
447          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[0]->getNum() &&
448          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[3]->getNum())
449         return -1;
450       if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[1]->getNum() &&
451          _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[0]->getNum() &&
452          _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[3]->getNum() &&
453          _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[2]->getNum())
454         return -1;
455       break;
456     default: return 0;
457     }
458   default: return 0;
459   }
460 }
461 
getTypeMSH() const462 int Cell::getTypeMSH() const
463 {
464   switch(_dim) {
465   case 0: return MSH_PNT;
466   case 1: return MSH_LIN_2;
467   case 2:
468     switch(getNumVertices()) {
469     case 3: return MSH_TRI_3;
470     case 4: return MSH_QUA_4;
471     default: return 0;
472     }
473   case 3:
474     switch(getNumVertices()) {
475     case 4: return MSH_TET_4;
476     case 5: return MSH_PYR_5;
477     case 6: return MSH_PRI_6;
478     case 8: return MSH_HEX_8;
479     default: return 0;
480     }
481   default: return 0;
482   }
483 }
484 
hasVertex(int vertex) const485 bool Cell::hasVertex(int vertex) const
486 {
487   std::vector<int> v;
488   for(std::size_t i = 0; i < _v.size(); i++) {
489     v.push_back(_v[(int)_si[i]]->getNum());
490   }
491   auto it = std::find(v.begin(), v.end(), vertex);
492   if(it != v.end())
493     return true;
494   else
495     return false;
496 }
497 
hasVertex(int vertex) const498 bool CombinedCell::hasVertex(int vertex) const
499 {
500   for(auto cit = _cells.begin(); cit != _cells.end(); cit++) {
501     if(cit->first->hasVertex(vertex)) return true;
502   }
503   return false;
504 }
505 
printCell()506 void Cell::printCell()
507 {
508   printf("%d-cell %d: \n", getDim(), getNum());
509   printf("  Vertices:");
510   for(int i = 0; i < this->getNumVertices(); i++) {
511     printf(" %lu", this->getMeshVertex(i)->getNum());
512   }
513   printf(", in subdomain: %d, ", inSubdomain());
514   printf("combined: %d. \n", isCombined());
515 };
516 
saveCellBoundary()517 void Cell::saveCellBoundary()
518 {
519   for(auto it = firstCoboundary(); it != lastCoboundary(); it++) {
520     it->second.init();
521   }
522   for(auto it = firstBoundary(); it != lastBoundary(); it++) {
523     it->second.init();
524   }
525 }
526 
restoreCellBoundary()527 void Cell::restoreCellBoundary()
528 {
529   std::vector<Cell *> toRemove;
530   for(auto it = firstCoboundary(true); it != lastCoboundary(); it++) {
531     it->second.reset();
532     if(it->second.get() == 0) toRemove.push_back(it->first);
533   }
534   for(std::size_t i = 0; i < toRemove.size(); i++) _cbd.erase(toRemove[i]);
535   toRemove.clear();
536   for(auto it = firstBoundary(true); it != lastBoundary(); it++) {
537     it->second.reset();
538     if(it->second.get() == 0) toRemove.push_back(it->first);
539   }
540   for(std::size_t i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]);
541 }
542 
addBoundaryCell(int orientation,Cell * cell,bool other)543 void Cell::addBoundaryCell(int orientation, Cell *cell, bool other)
544 {
545   auto it = _bd.find(cell);
546   if(it != _bd.end()) {
547     int newOrientation = it->second.get() + orientation;
548     it->second.set(newOrientation);
549     if(newOrientation == 0) {
550       it->first->removeCoboundaryCell(this, false);
551       if(it->second.geto() == 0) { _bd.erase(it); }
552       return;
553     }
554   }
555   else
556     _bd.insert(std::make_pair(cell, BdInfo(orientation)));
557   if(other) cell->addCoboundaryCell(orientation, this, false);
558 }
559 
addCoboundaryCell(int orientation,Cell * cell,bool other)560 void Cell::addCoboundaryCell(int orientation, Cell *cell, bool other)
561 {
562   auto it = _cbd.find(cell);
563   if(it != _cbd.end()) {
564     int newOrientation = it->second.get() + orientation;
565     it->second.set(newOrientation);
566     if(newOrientation == 0) {
567       it->first->removeBoundaryCell(this, false);
568       if(it->second.geto() == 0) { _cbd.erase(it); }
569       return;
570     }
571   }
572   else
573     _cbd.insert(std::make_pair(cell, BdInfo(orientation)));
574   if(other) cell->addBoundaryCell(orientation, this, false);
575 }
576 
removeBoundaryCell(Cell * cell,bool other)577 void Cell::removeBoundaryCell(Cell *cell, bool other)
578 {
579   auto it = _bd.find(cell);
580   if(it != _bd.end()) {
581     it->second.set(0);
582     if(other) it->first->removeCoboundaryCell(this, false);
583     if(it->second.geto() == 0) _bd.erase(it);
584   }
585 }
586 
removeCoboundaryCell(Cell * cell,bool other)587 void Cell::removeCoboundaryCell(Cell *cell, bool other)
588 {
589   auto it = _cbd.find(cell);
590   if(it != _cbd.end()) {
591     it->second.set(0);
592     if(other) it->first->removeBoundaryCell(this, false);
593     if(it->second.geto() == 0) _cbd.erase(it);
594   }
595 }
596 
hasBoundary(Cell * cell,bool orig)597 bool Cell::hasBoundary(Cell *cell, bool orig)
598 {
599   if(!orig) {
600     auto it = _bd.find(cell);
601     if(it != _bd.end() && it->second.get() != 0) return true;
602     return false;
603   }
604   else {
605     auto it = _bd.find(cell);
606     if(it != _bd.end() && it->second.geto() != 0) return true;
607     return false;
608   }
609 }
610 
hasCoboundary(Cell * cell,bool orig)611 bool Cell::hasCoboundary(Cell *cell, bool orig)
612 {
613   if(!orig) {
614     auto it = _cbd.find(cell);
615     if(it != _cbd.end() && it->second.get() != 0) return true;
616     return false;
617   }
618   else {
619     auto it = _cbd.find(cell);
620     if(it != _cbd.end() && it->second.geto() != 0) return true;
621     return false;
622   }
623 }
624 
firstBoundary(bool orig)625 Cell::biter Cell::firstBoundary(bool orig)
626 {
627   auto it = _bd.begin();
628   if(!orig)
629     while(it->second.get() == 0 && it != _bd.end()) it++;
630   else
631     while(it->second.geto() == 0 && it != _bd.end()) it++;
632   return it;
633 }
634 
lastBoundary()635 Cell::biter Cell::lastBoundary() { return _bd.end(); }
636 
firstCoboundary(bool orig)637 Cell::biter Cell::firstCoboundary(bool orig)
638 {
639   auto it = _cbd.begin();
640   if(!orig)
641     while(it->second.get() == 0 && it != _cbd.end()) it++;
642   else
643     while(it->second.geto() == 0 && it != _cbd.end()) it++;
644   return it;
645 }
646 
lastCoboundary()647 Cell::biter Cell::lastCoboundary() { return _cbd.end(); }
648 
getBoundarySize(bool orig)649 int Cell::getBoundarySize(bool orig)
650 {
651   int size = 0;
652   for(auto bit = _bd.begin(); bit != _bd.end(); bit++) {
653     if(!orig && bit->second.get() != 0)
654       size++;
655     else if(orig && bit->second.geto() != 0)
656       size++;
657   }
658   return size;
659 }
660 
getCoboundarySize(bool orig)661 int Cell::getCoboundarySize(bool orig)
662 {
663   int size = 0;
664   for(auto bit = _cbd.begin(); bit != _cbd.end(); bit++) {
665     if(!orig && bit->second.get() != 0)
666       size++;
667     else if(orig && bit->second.geto() != 0)
668       size++;
669   }
670   return size;
671 }
672 
getBoundary(std::map<Cell *,short int,CellPtrLessThan> & boundary,bool orig)673 void Cell::getBoundary(std::map<Cell *, short int, CellPtrLessThan> &boundary,
674                        bool orig)
675 {
676   boundary.clear();
677   for(auto it = firstBoundary(); it != lastBoundary(); it++) {
678     Cell *cell = it->first;
679     if(!orig && it->second.get() != 0) boundary[cell] = it->second.get();
680     if(orig && it->second.geto() != 0) boundary[cell] = it->second.geto();
681   }
682 }
683 
getCoboundary(std::map<Cell *,short int,CellPtrLessThan> & coboundary,bool orig)684 void Cell::getCoboundary(
685   std::map<Cell *, short int, CellPtrLessThan> &coboundary, bool orig)
686 {
687   coboundary.clear();
688   for(auto it = firstCoboundary(); it != lastCoboundary(); it++) {
689     Cell *cell = it->first;
690     if(!orig && it->second.get() != 0) coboundary[cell] = it->second.get();
691     if(orig && it->second.geto() != 0) coboundary[cell] = it->second.geto();
692   }
693 }
694 
printBoundary()695 void Cell::printBoundary()
696 {
697   for(auto it = firstBoundary(); it != lastBoundary(); it++) {
698     printf("Boundary cell orientation: %d ", it->second.get());
699     Cell *cell2 = it->first;
700     cell2->printCell();
701   }
702   if(firstBoundary() == lastBoundary()) {
703     printf("Cell boundary is empty. \n");
704   }
705 }
706 
printCoboundary()707 void Cell::printCoboundary()
708 {
709   for(auto it = firstCoboundary(); it != lastCoboundary(); it++) {
710     printf("Coboundary cell orientation: %d, ", it->second.get());
711     Cell *cell2 = it->first;
712     cell2->printCell();
713     if(firstCoboundary() == lastCoboundary()) {
714       printf("Cell coboundary is empty. \n");
715     }
716   }
717 }
718 
CombinedCell(Cell * c1,Cell * c2,bool orMatch,bool co)719 CombinedCell::CombinedCell(Cell *c1, Cell *c2, bool orMatch, bool co)
720 {
721   // use "smaller" cell as c2
722   if(c1->getNumCells() < c2->getNumCells()) {
723     Cell *temp = c1;
724     c1 = c2;
725     c2 = temp;
726   }
727 
728   _num = ++_globalNum;
729   _domain = c1->getDomain();
730   _combined = true;
731   _immune = (c1->getImmune() || c2->getImmune());
732 
733   // cells
734   c1->getCells(_cells);
735   std::map<Cell *, int, CellPtrLessThan> c2Cells;
736   c2->getCells(c2Cells);
737   for(auto cit = c2Cells.begin(); cit != c2Cells.end(); cit++) {
738     if(!orMatch) (*cit).second = -1 * (*cit).second;
739     _cells.insert(*cit);
740   }
741 
742   // boundary cells
743   for(auto it = c1->firstBoundary(); it != c1->lastBoundary(); it++) {
744     Cell *cell = it->first;
745     int ori = it->second.get();
746     if(ori == 0) continue;
747     cell->removeCoboundaryCell(c1, false);
748     this->addBoundaryCell(ori, cell, true);
749   }
750   for(auto it = c2->firstBoundary(); it != c2->lastBoundary(); it++) {
751     Cell *cell = it->first;
752     if(!orMatch) it->second.set(-1 * it->second.get());
753     int ori = it->second.get();
754     if(ori == 0) continue;
755     cell->removeCoboundaryCell(c2, false);
756     if(co && !c1->hasBoundary(cell)) { this->addBoundaryCell(ori, cell, true); }
757     else if(!co)
758       this->addBoundaryCell(ori, cell, true);
759   }
760 
761   // coboundary cells
762   for(auto it = c1->firstCoboundary(); it != c1->lastCoboundary(); it++) {
763     Cell *cell = it->first;
764     int ori = it->second.get();
765     if(ori == 0) continue;
766     cell->removeBoundaryCell(c1, false);
767     this->addCoboundaryCell(ori, cell, true);
768   }
769   for(auto it = c2->firstCoboundary(); it != c2->lastCoboundary(); it++) {
770     Cell *cell = it->first;
771     if(!orMatch) it->second.set(-1 * it->second.get());
772     int ori = it->second.get();
773     if(ori == 0) continue;
774     cell->removeBoundaryCell(c2, false);
775     if(!co && !c1->hasCoboundary(cell)) {
776       this->addCoboundaryCell(ori, cell, true);
777     }
778     else if(co)
779       this->addCoboundaryCell(ori, cell, true);
780   }
781 }
782 
CombinedCell(std::vector<Cell * > & cells)783 CombinedCell::CombinedCell(std::vector<Cell *> &cells)
784 {
785   _num = ++_globalNum;
786   _domain = cells.at(0)->getDomain();
787   _combined = true;
788   _immune = false;
789 
790   // cells
791   for(std::size_t i = 0; i < cells.size(); i++) {
792     Cell *c = cells.at(i);
793     if(c->getImmune()) _immune = true;
794     _cells[c] = 1;
795   }
796 }
797