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