Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "EFAElement3D.h"
11 :
12 : #include <iomanip>
13 :
14 : #include "EFAFaceNode.h"
15 : #include "EFAVolumeNode.h"
16 : #include "EFANode.h"
17 : #include "EFAEdge.h"
18 : #include "EFAFace.h"
19 : #include "EFAFragment3D.h"
20 : #include "EFAFuncs.h"
21 : #include "EFAError.h"
22 : #include "XFEMFuncs.h"
23 :
24 : #include "libmesh/int_range.h"
25 :
26 47798 : EFAElement3D::EFAElement3D(unsigned int eid, unsigned int n_nodes, unsigned int n_faces)
27 : : EFAElement(eid, n_nodes),
28 47798 : _num_faces(n_faces),
29 47798 : _faces(_num_faces, nullptr),
30 47798 : _face_neighbors(_num_faces),
31 95596 : _face_edge_neighbors(_num_faces)
32 : {
33 47798 : if (_num_faces == 4)
34 : {
35 21368 : _num_vertices = 4;
36 21368 : if (_num_nodes == 14)
37 10684 : _num_interior_face_nodes = 4;
38 10684 : else if (_num_nodes == 10)
39 10684 : _num_interior_face_nodes = 3;
40 0 : else if (_num_nodes == 4)
41 0 : _num_interior_face_nodes = 0;
42 : else
43 0 : EFAError("In EFAelement3D the supported TET element types are TET4, TET10 and TET14");
44 : }
45 26430 : else if (_num_faces == 6)
46 : {
47 26430 : _num_vertices = 8;
48 26430 : if (_num_nodes == 27)
49 388 : _num_interior_face_nodes = 5;
50 26042 : else if (_num_nodes == 20)
51 388 : _num_interior_face_nodes = 4;
52 25654 : else if (_num_nodes == 8)
53 25654 : _num_interior_face_nodes = 0;
54 : else
55 0 : EFAError("In EFAelement3D the supported HEX element types are HEX8, HEX20 and HEX27");
56 : }
57 : else
58 0 : EFAError("In EFAelement3D the supported element types are TET4, TET10, TET14, HEX8, HEX20 and "
59 : "HEX27");
60 47798 : setLocalCoordinates();
61 47798 : }
62 :
63 2680 : EFAElement3D::EFAElement3D(const EFAElement3D * from_elem, bool convert_to_local)
64 2680 : : EFAElement(from_elem->_id, from_elem->_num_nodes),
65 2680 : _num_faces(from_elem->_num_faces),
66 2680 : _faces(_num_faces, nullptr),
67 2680 : _face_neighbors(_num_faces),
68 5360 : _face_edge_neighbors(_num_faces)
69 : {
70 2680 : if (convert_to_local)
71 : {
72 : // build local nodes from global nodes
73 30480 : for (unsigned int i = 0; i < _num_nodes; ++i)
74 : {
75 27800 : if (from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT ||
76 0 : from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_TEMP)
77 : {
78 27800 : _nodes[i] = from_elem->createLocalNodeFromGlobalNode(from_elem->_nodes[i]);
79 27800 : _local_nodes.push_back(_nodes[i]); // convenient to delete local nodes
80 : }
81 : else
82 0 : EFAError("In EFAelement3D ",
83 : from_elem->id(),
84 : " the copy constructor must have from_elem w/ global nodes. node: ",
85 : i,
86 : " category: ",
87 : from_elem->_nodes[i]->category());
88 : }
89 :
90 : // copy faces, fragments and interior nodes from from_elem
91 16200 : for (unsigned int i = 0; i < _num_faces; ++i)
92 13520 : _faces[i] = new EFAFace(*from_elem->_faces[i]);
93 5360 : for (unsigned int i = 0; i < from_elem->_fragments.size(); ++i)
94 2680 : _fragments.push_back(new EFAFragment3D(this, true, from_elem, i));
95 2680 : for (unsigned int i = 0; i < from_elem->_interior_nodes.size(); ++i)
96 0 : _interior_nodes.push_back(new EFAVolumeNode(*from_elem->_interior_nodes[i]));
97 :
98 : // replace all global nodes with local nodes
99 30480 : for (unsigned int i = 0; i < _num_nodes; ++i)
100 : {
101 27800 : if (_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
102 27800 : switchNode(
103 : _nodes[i],
104 : from_elem->_nodes[i],
105 : false); // when save to _cut_elem_map, the EFAelement is not a child of any parent
106 : else
107 0 : EFAError("In EFAelement3D copy constructor this elem's nodes must be local");
108 : }
109 :
110 : // create element face connectivity array (IMPORTANT)
111 2680 : findFacesAdjacentToFaces();
112 :
113 2680 : _local_node_coor = from_elem->_local_node_coor;
114 2680 : _num_interior_face_nodes = from_elem->_num_interior_face_nodes;
115 2680 : _num_vertices = from_elem->_num_vertices;
116 : }
117 : else
118 0 : EFAError("this EFAelement3D constructor only converts global nodes to local nodes");
119 2680 : }
120 :
121 98276 : EFAElement3D::~EFAElement3D()
122 : {
123 65428 : for (unsigned int i = 0; i < _fragments.size(); ++i)
124 14950 : if (_fragments[i])
125 14950 : delete _fragments[i];
126 :
127 308050 : for (unsigned int i = 0; i < _faces.size(); ++i)
128 257572 : if (_faces[i])
129 257572 : delete _faces[i];
130 :
131 50478 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
132 0 : if (_interior_nodes[i])
133 0 : delete _interior_nodes[i];
134 :
135 78278 : for (unsigned int i = 0; i < _local_nodes.size(); ++i)
136 27800 : if (_local_nodes[i])
137 27800 : delete _local_nodes[i];
138 98276 : }
139 :
140 : void
141 47798 : EFAElement3D::setLocalCoordinates()
142 : {
143 47798 : if (_num_faces == 6)
144 : {
145 : /*
146 : HEX27(HEX20): 7 18 6
147 : o--------------o--------------o
148 : /: / /|
149 : / : / / |
150 : / : / / |
151 : 19/ : 25/ 17/ |
152 : o--------------o--------------o |
153 : / : / /| |
154 : / 15o / 23o / | 14o
155 : / : / / | /|
156 : 4/ : 16/ 5/ | / |
157 : o--------------o--------------o | / |
158 : | : | 26 | |/ |
159 : | 24o : | o | 22o |
160 : | : | 10 | /| |
161 : | 3o....|.........o....|../.|....o
162 : | . | | / | / 2
163 : | . 21| 13|/ | /
164 : 12 o--------------o--------------o | /
165 : | . | | |/
166 : | 11o | 20o | o
167 : | . | | / 9
168 : | . | | /
169 : | . | | /
170 : |. | |/
171 : o--------------o--------------o
172 : 0 8 1
173 :
174 : */
175 26430 : _local_node_coor.resize(_num_nodes);
176 26430 : _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
177 26430 : _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
178 26430 : _local_node_coor[2] = EFAPoint(1.0, 1.0, 0.0);
179 26430 : _local_node_coor[3] = EFAPoint(0.0, 1.0, 0.0);
180 26430 : _local_node_coor[4] = EFAPoint(0.0, 0.0, 1.0);
181 26430 : _local_node_coor[5] = EFAPoint(1.0, 0.0, 1.0);
182 26430 : _local_node_coor[6] = EFAPoint(1.0, 1.0, 1.0);
183 26430 : _local_node_coor[7] = EFAPoint(0.0, 1.0, 1.0);
184 :
185 26430 : if (_num_nodes > 8)
186 : {
187 776 : _local_node_coor[8] = EFAPoint(0.5, 0.0, 0.0);
188 776 : _local_node_coor[9] = EFAPoint(1.0, 0.5, 0.0);
189 776 : _local_node_coor[10] = EFAPoint(0.5, 1.0, 0.0);
190 776 : _local_node_coor[11] = EFAPoint(0.0, 0.5, 0.0);
191 776 : _local_node_coor[12] = EFAPoint(0.0, 0.0, 0.5);
192 776 : _local_node_coor[13] = EFAPoint(1.0, 0.0, 0.5);
193 776 : _local_node_coor[14] = EFAPoint(1.0, 1.0, 0.5);
194 776 : _local_node_coor[15] = EFAPoint(0.0, 1.0, 0.5);
195 776 : _local_node_coor[16] = EFAPoint(0.5, 0.0, 1.0);
196 776 : _local_node_coor[17] = EFAPoint(1.0, 0.5, 1.0);
197 776 : _local_node_coor[18] = EFAPoint(0.5, 1.0, 1.0);
198 776 : _local_node_coor[19] = EFAPoint(0.0, 0.5, 1.0);
199 : }
200 :
201 26430 : if (_num_nodes > 20)
202 : {
203 388 : _local_node_coor[20] = EFAPoint(0.5, 0.5, 0.0);
204 388 : _local_node_coor[21] = EFAPoint(0.5, 0.0, 0.5);
205 388 : _local_node_coor[22] = EFAPoint(1.0, 0.5, 0.5);
206 388 : _local_node_coor[23] = EFAPoint(0.5, 1.0, 0.5);
207 388 : _local_node_coor[24] = EFAPoint(0.0, 0.5, 0.5);
208 388 : _local_node_coor[25] = EFAPoint(0.5, 0.5, 1.0);
209 388 : _local_node_coor[26] = EFAPoint(0.5, 0.5, 0.5);
210 : }
211 : }
212 21368 : else if (_num_faces == 4)
213 : {
214 : /*
215 : 3
216 : TET10: o
217 : /|\
218 : / | \
219 : 7 / | \9
220 : o | o
221 : / |8 \
222 : / o \
223 : / 6 | \
224 : 0 o.....o.|.......o 2
225 : \ | /
226 : \ | /
227 : \ | /
228 : 4 o | o 5
229 : \ | /
230 : \ | /
231 : \|/
232 : o
233 : 1
234 :
235 : */
236 : /*
237 : TET14 elements also include four face nodes:
238 : Node 10, centroid on side 0, arithmetic mean of 0/1/2 or 4/5/6
239 : Node 11, centroid on side 1, arithmetic mean of 0/1/3 or 4/7/8
240 : Node 12, centroid on side 2, arithmetic mean of 1/2/3 or 5/8/9
241 : Node 13, centroid on side 3, arithmetic mean of 0/2/3 or 6/7/9
242 : */
243 21368 : _local_node_coor.resize(_num_nodes);
244 21368 : _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
245 21368 : _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
246 21368 : _local_node_coor[2] = EFAPoint(0.0, 1.0, 0.0);
247 21368 : _local_node_coor[3] = EFAPoint(0.0, 0.0, 1.0);
248 :
249 21368 : if (_num_nodes > 4)
250 : {
251 21368 : _local_node_coor[4] = EFAPoint(0.5, 0.0, 0.0);
252 21368 : _local_node_coor[5] = EFAPoint(0.5, 0.5, 0.0);
253 21368 : _local_node_coor[6] = EFAPoint(0.0, 0.5, 0.0);
254 21368 : _local_node_coor[7] = EFAPoint(0.0, 0.0, 0.5);
255 21368 : _local_node_coor[8] = EFAPoint(0.5, 0.0, 0.5);
256 21368 : _local_node_coor[9] = EFAPoint(0.0, 0.5, 0.5);
257 : }
258 :
259 21368 : if (_num_nodes > 10)
260 : {
261 10684 : _local_node_coor[10] = EFAPoint(1 / 3., 1 / 3., 0.0);
262 10684 : _local_node_coor[11] = EFAPoint(1 / 3., 0.0, 1 / 3.);
263 10684 : _local_node_coor[12] = EFAPoint(1 / 3., 1 / 3., 1 / 3.);
264 10684 : _local_node_coor[13] = EFAPoint(0.0, 1 / 3., 1 / 3.);
265 : }
266 : }
267 : else
268 0 : EFAError("EFAElement3D: number of faces should be either 4(TET) or 6(HEX).");
269 47798 : }
270 :
271 : unsigned int
272 682182 : EFAElement3D::numFragments() const
273 : {
274 682182 : return _fragments.size();
275 : }
276 :
277 : bool
278 952090 : EFAElement3D::isPartial() const
279 : {
280 952090 : if (_fragments.size() > 0)
281 : {
282 2364923 : for (unsigned int i = 0; i < _num_vertices; ++i)
283 : {
284 : bool node_in_frag = false;
285 3108429 : for (unsigned int j = 0; j < _fragments.size(); ++j)
286 : {
287 2260631 : if (_fragments[j]->containsNode(_nodes[i]))
288 : {
289 : node_in_frag = true;
290 : break;
291 : }
292 : } // j
293 :
294 2260631 : if (!node_in_frag)
295 : return true;
296 : } // i
297 : }
298 : return false;
299 : }
300 :
301 : void
302 3828 : EFAElement3D::getNonPhysicalNodes(std::set<EFANode *> & non_physical_nodes) const
303 : {
304 : // Any nodes that don't belong to any fragment are non-physical
305 : // First add all nodes in the element to the set
306 42402 : for (unsigned int i = 0; i < _nodes.size(); ++i)
307 38574 : non_physical_nodes.insert(_nodes[i]);
308 :
309 : // Now delete any nodes that are contained in fragments
310 : std::set<EFANode *>::iterator sit;
311 42402 : for (sit = non_physical_nodes.begin(); sit != non_physical_nodes.end();)
312 : {
313 : bool erased = false;
314 63876 : for (unsigned int i = 0; i < _fragments.size(); ++i)
315 : {
316 38574 : if (_fragments[i]->containsNode(*sit))
317 : {
318 13272 : non_physical_nodes.erase(sit++);
319 : erased = true;
320 13272 : break;
321 : }
322 : }
323 38574 : if (!erased)
324 : ++sit;
325 : }
326 3828 : }
327 :
328 : void
329 887706 : EFAElement3D::switchNode(EFANode * new_node, EFANode * old_node, bool descend_to_parent)
330 : {
331 : // We are not switching any embedded nodes here; This is an enhanced version
332 11408636 : for (unsigned int i = 0; i < _num_nodes; ++i)
333 : {
334 10520930 : if (_nodes[i] == old_node)
335 13741 : _nodes[i] = new_node;
336 : }
337 1441777 : for (unsigned int i = 0; i < _fragments.size(); ++i)
338 554071 : _fragments[i]->switchNode(new_node, old_node);
339 :
340 4728470 : for (unsigned int i = 0; i < _faces.size(); ++i)
341 3840764 : _faces[i]->switchNode(new_node, old_node);
342 :
343 887706 : if (_parent && descend_to_parent)
344 : {
345 11793 : _parent->switchNode(new_node, old_node, false);
346 576876 : for (unsigned int i = 0; i < _parent->numGeneralNeighbors(); ++i)
347 : {
348 565083 : EFAElement * neigh_elem = _parent->getGeneralNeighbor(i); // generalized neighbor element
349 1331649 : for (unsigned int k = 0; k < neigh_elem->numChildren(); ++k)
350 766566 : neigh_elem->getChild(k)->switchNode(new_node, old_node, false);
351 : }
352 : }
353 887706 : }
354 :
355 : void
356 0 : EFAElement3D::switchEmbeddedNode(EFANode * new_emb_node, EFANode * old_emb_node)
357 : {
358 0 : for (unsigned int i = 0; i < _num_faces; ++i)
359 0 : _faces[i]->switchNode(new_emb_node, old_emb_node);
360 0 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
361 0 : _interior_nodes[i]->switchNode(new_emb_node, old_emb_node);
362 0 : for (unsigned int i = 0; i < _fragments.size(); ++i)
363 0 : _fragments[i]->switchNode(new_emb_node, old_emb_node);
364 0 : }
365 :
366 : void
367 2831 : EFAElement3D::updateFragmentNode()
368 : {
369 : // In EFAElement3D, updateFragmentNode needs to be implemented
370 2831 : }
371 :
372 : void
373 3181742 : EFAElement3D::getMasterInfo(EFANode * node,
374 : std::vector<EFANode *> & master_nodes,
375 : std::vector<double> & master_weights) const
376 : {
377 : // Given a EFAnode, return its master nodes and weights
378 3181742 : master_nodes.clear();
379 3181742 : master_weights.clear();
380 : bool masters_found = false;
381 6998975 : for (unsigned int i = 0; i < _num_faces; ++i) // check element exterior faces
382 : {
383 6998975 : if (_faces[i]->containsNode(node))
384 : {
385 3181742 : masters_found = _faces[i]->getMasterInfo(node, master_nodes, master_weights);
386 3181742 : if (masters_found)
387 : break;
388 : else
389 0 : EFAError("In getMasterInfo: cannot find master nodes in element faces");
390 : }
391 : }
392 :
393 : if (!masters_found) // check element interior embedded nodes
394 : {
395 0 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
396 : {
397 0 : if (_interior_nodes[i]->getNode() == node)
398 : {
399 0 : std::vector<double> xi_3d(3, -100.0);
400 0 : for (unsigned int j = 0; j < 3; ++j)
401 0 : xi_3d[j] = _interior_nodes[i]->getParametricCoordinates(j);
402 0 : for (unsigned int j = 0; j < _num_nodes; ++j)
403 : {
404 0 : master_nodes.push_back(_nodes[j]);
405 0 : double weight = 0.0;
406 0 : if (_num_nodes == 8)
407 0 : weight = Efa::linearHexShape3D(j, xi_3d);
408 0 : else if (_num_nodes == 4)
409 0 : weight = Efa::linearTetShape3D(j, xi_3d);
410 : else
411 0 : EFAError("unknown 3D element");
412 0 : master_weights.push_back(weight);
413 : }
414 : masters_found = true;
415 : break;
416 0 : }
417 : }
418 : }
419 :
420 3181742 : if (!masters_found)
421 0 : EFAError("In EFAelement3D::getMasterInfo, cannot find the given EFANode");
422 3181742 : }
423 :
424 : unsigned int
425 1829 : EFAElement3D::numInteriorNodes() const
426 : {
427 1829 : return _interior_nodes.size();
428 : }
429 :
430 : bool
431 538058 : EFAElement3D::overlaysElement(const EFAElement3D * other_elem) const
432 : {
433 : bool overlays = false;
434 : const EFAElement3D * other3d = dynamic_cast<const EFAElement3D *>(other_elem);
435 538058 : if (!other3d)
436 0 : EFAError("failed to dynamic cast to other3d");
437 :
438 : // Find indices of common nodes
439 538058 : const auto & common_face_curr = getCommonFaceID(other3d);
440 538058 : if (common_face_curr.size() == 1)
441 : {
442 176792 : unsigned int curr_face_id = common_face_curr[0];
443 176792 : EFAFace * curr_face = _faces[curr_face_id];
444 176792 : unsigned int other_face_id = other3d->getFaceID(curr_face);
445 176792 : EFAFace * other_face = other3d->_faces[other_face_id];
446 176792 : if (curr_face->hasSameOrientation(other_face))
447 : overlays = true;
448 : }
449 361266 : else if (common_face_curr.size() > 1)
450 : {
451 : // TODO: We probably need more error checking here.
452 : overlays = true;
453 : }
454 538058 : return overlays;
455 538058 : }
456 :
457 : unsigned int
458 189763 : EFAElement3D::getNeighborIndex(const EFAElement * neighbor_elem) const
459 : {
460 595566 : for (unsigned int i = 0; i < _num_faces; ++i)
461 890153 : for (unsigned int j = 0; j < _face_neighbors[i].size(); ++j)
462 484350 : if (_face_neighbors[i][j] == neighbor_elem)
463 189763 : return i;
464 0 : EFAError("in getNeighborIndex() element ", _id, " does not have neighbor ", neighbor_elem->id());
465 : }
466 :
467 : void
468 40204 : EFAElement3D::getNeighborEdgeIndex(const EFAElement3D * neighbor_elem,
469 : unsigned int face_id,
470 : unsigned int edge_id,
471 : unsigned int & neigh_face_id,
472 : unsigned int & neigh_edge_id) const
473 : {
474 40204 : EFAEdge * edge = this->getFace(face_id)->getEdge(edge_id);
475 84162 : for (unsigned int i = 0; i < neighbor_elem->numFaces(); ++i)
476 : {
477 281234 : for (unsigned int j = 0; j < neighbor_elem->getFace(i)->numEdges(); ++j)
478 : {
479 237276 : EFAEdge * neigh_edge = neighbor_elem->getFace(i)->getEdge(j);
480 237276 : if (neigh_edge->equivalent(*edge))
481 : {
482 40204 : neigh_face_id = i;
483 40204 : neigh_edge_id = j;
484 40204 : return;
485 : }
486 : }
487 : }
488 0 : EFAError("in getNeighborEdgeIndex() element ",
489 : _id,
490 : " does not share a common edge with element",
491 : neighbor_elem->id());
492 : }
493 :
494 : void
495 45225 : EFAElement3D::clearNeighbors()
496 : {
497 45225 : _general_neighbors.clear();
498 276399 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
499 : {
500 231174 : _face_neighbors[face_iter].clear();
501 1075518 : for (unsigned int edge_iter = 0; edge_iter < _faces[face_iter]->numEdges(); ++edge_iter)
502 844344 : _face_edge_neighbors[face_iter][edge_iter].clear();
503 : }
504 45225 : }
505 :
506 : void
507 45225 : EFAElement3D::setupNeighbors(std::map<EFANode *, std::set<EFAElement *>> & inverse_connectivity_map)
508 : {
509 45225 : findGeneralNeighbors(inverse_connectivity_map);
510 : std::vector<std::pair<unsigned int, unsigned int>> common_ids;
511 1603515 : for (unsigned int eit2 = 0; eit2 < _general_neighbors.size(); ++eit2)
512 : {
513 1558290 : EFAElement3D * neigh_elem = dynamic_cast<EFAElement3D *>(_general_neighbors[eit2]);
514 1558290 : if (!neigh_elem)
515 0 : EFAError("neighbor_elem is not of EFAelement3D type");
516 :
517 1558290 : const auto & common_face_id = getCommonFaceID(neigh_elem);
518 1919556 : if (common_face_id.empty() && getCommonEdgeID(neigh_elem, common_ids) &&
519 361266 : !overlaysElement(neigh_elem))
520 : {
521 : bool is_edge_neighbor = false;
522 :
523 : // Fragments must match up.
524 361266 : if ((_fragments.size() > 1) || (neigh_elem->numFragments() > 1))
525 0 : EFAError("in updateFaceNeighbors: Cannot have more than 1 fragment");
526 :
527 361266 : if ((_fragments.size() == 1) && (neigh_elem->numFragments() == 1))
528 : {
529 32172 : if (_fragments[0]->isEdgeConnected(neigh_elem->getFragment(0)))
530 : is_edge_neighbor = true;
531 : }
532 : else // If there are no fragments to match up, consider them edge neighbors
533 : is_edge_neighbor = true;
534 :
535 : if (is_edge_neighbor)
536 1075256 : for (const auto & [face_id, edge_id] : common_ids)
537 716924 : _face_edge_neighbors[face_id][edge_id].push_back(neigh_elem);
538 : }
539 :
540 1558290 : if (common_face_id.size() == 1 && !overlaysElement(neigh_elem))
541 : {
542 174488 : unsigned int face_id = common_face_id[0];
543 : bool is_face_neighbor = false;
544 :
545 : // Fragments must match up.
546 174488 : if ((_fragments.size() > 1) || (neigh_elem->numFragments() > 1))
547 0 : EFAError("in updateFaceNeighbors: Cannot have more than 1 fragment");
548 :
549 174488 : if ((_fragments.size() == 1) && (neigh_elem->numFragments() == 1))
550 : {
551 18924 : if (_fragments[0]->isConnected(neigh_elem->getFragment(0)))
552 : is_face_neighbor = true;
553 : }
554 : else // If there are no fragments to match up, consider them neighbors
555 : is_face_neighbor = true;
556 :
557 : if (is_face_neighbor)
558 : {
559 174344 : if (_face_neighbors[face_id].size() > 1)
560 0 : EFAError("Element ",
561 : _id,
562 : " already has 2 face neighbors: ",
563 : _face_neighbors[face_id][0]->id(),
564 : " ",
565 : _face_neighbors[face_id][1]->id());
566 174344 : _face_neighbors[face_id].push_back(neigh_elem);
567 : }
568 : }
569 1558290 : }
570 45225 : }
571 :
572 : void
573 45225 : EFAElement3D::neighborSanityCheck() const
574 : {
575 276399 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
576 : {
577 405518 : for (unsigned int en_iter = 0; en_iter < _face_neighbors[face_iter].size(); ++en_iter)
578 : {
579 174344 : EFAElement3D * neigh_elem = _face_neighbors[face_iter][en_iter];
580 174344 : if (neigh_elem != nullptr)
581 : {
582 : bool found_neighbor = false;
583 1074976 : for (unsigned int face_iter2 = 0; face_iter2 < neigh_elem->numFaces(); ++face_iter2)
584 : {
585 1426492 : for (unsigned int en_iter2 = 0; en_iter2 < neigh_elem->numFaceNeighbors(face_iter2);
586 : ++en_iter2)
587 : {
588 700204 : if (neigh_elem->getFaceNeighbor(face_iter2, en_iter2) == this)
589 : {
590 174344 : if ((en_iter2 > 1) && (en_iter > 1))
591 0 : EFAError(
592 : "Element and neighbor element cannot both have >1 neighbors on a common face");
593 : found_neighbor = true;
594 : break;
595 : }
596 : }
597 : }
598 174344 : if (!found_neighbor)
599 0 : EFAError("Neighbor element doesn't recognize current element as neighbor");
600 : }
601 : }
602 : }
603 45225 : }
604 :
605 : void
606 45048 : EFAElement3D::initCrackTip(std::set<EFAElement *> & CrackTipElements)
607 : {
608 45048 : if (isCrackTipElement())
609 : {
610 617 : CrackTipElements.insert(this);
611 4031 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
612 : {
613 3414 : if ((_face_neighbors[face_iter].size() == 2) && (_faces[face_iter]->hasIntersection()))
614 : {
615 : // Neither neighbor overlays current element. We are on the uncut element ahead of the tip.
616 : // Flag neighbors as crack tip split elements and add this element as their crack tip
617 : // neighbor.
618 1256 : if (_face_neighbors[face_iter][0]->overlaysElement(this) ||
619 628 : _face_neighbors[face_iter][1]->overlaysElement(this))
620 0 : EFAError("Element has a neighbor that overlays itself");
621 :
622 : // Make sure the current elment hasn't been flagged as a tip element
623 628 : if (_crack_tip_split_element)
624 0 : EFAError("crack_tip_split_element already flagged. In elem: ",
625 : _id,
626 : " flags: ",
627 : _crack_tip_split_element,
628 : " ",
629 : _face_neighbors[face_iter][0]->isCrackTipSplit(),
630 : " ",
631 : _face_neighbors[face_iter][1]->isCrackTipSplit());
632 :
633 628 : _face_neighbors[face_iter][0]->setCrackTipSplit();
634 628 : _face_neighbors[face_iter][1]->setCrackTipSplit();
635 :
636 628 : _face_neighbors[face_iter][0]->addCrackTipNeighbor(this);
637 628 : _face_neighbors[face_iter][1]->addCrackTipNeighbor(this);
638 : }
639 : } // face_iter
640 : }
641 45048 : }
642 :
643 : bool
644 32081 : EFAElement3D::shouldDuplicateForCrackTip(const std::set<EFAElement *> & CrackTipElements)
645 : {
646 : // This method is called in createChildElements()
647 : // Only duplicate when
648 : // 1) currElem will be a NEW crack tip element
649 : // 2) currElem is a crack tip split element at last time step and the tip will extend
650 : // 3) currElem is the neighbor of a to-be-second-split element which has another neighbor
651 : // sharing a phantom node with currElem
652 : bool should_duplicate = false;
653 32081 : if (_fragments.size() == 1)
654 : {
655 : std::set<EFAElement *>::iterator sit;
656 7275 : sit = CrackTipElements.find(this);
657 7275 : if (sit == CrackTipElements.end() && isCrackTipElement())
658 : should_duplicate = true;
659 4548 : else if (shouldDuplicateCrackTipSplitElement(CrackTipElements))
660 : should_duplicate = true;
661 3828 : else if (shouldDuplicateForPhantomCorner())
662 : should_duplicate = true;
663 : }
664 32081 : return should_duplicate;
665 : }
666 :
667 : bool
668 4548 : EFAElement3D::shouldDuplicateCrackTipSplitElement(const std::set<EFAElement *> & CrackTipElements)
669 : {
670 : // Determine whether element at crack tip should be duplicated. It should be duplicated
671 : // if the crack will extend into the next element, or if it has a non-physical node
672 : // connected to a face where a crack terminates, but will extend.
673 :
674 : bool should_duplicate = false;
675 4548 : if (_fragments.size() == 1)
676 : {
677 : std::vector<unsigned int> split_neighbors;
678 4548 : if (willCrackTipExtend(split_neighbors))
679 : should_duplicate = true;
680 : else
681 : {
682 : // The element may not be at the crack tip, but could have a non-physical node
683 : // connected to a crack tip face (on a neighbor element) that will be split. We need to
684 : // duplicate in that case as well.
685 : std::set<EFANode *> non_physical_nodes;
686 3828 : getNonPhysicalNodes(non_physical_nodes);
687 :
688 133984 : for (unsigned int eit = 0; eit < _general_neighbors.size(); ++eit)
689 : {
690 130156 : EFAElement3D * neigh_elem = dynamic_cast<EFAElement3D *>(_general_neighbors[eit]);
691 130156 : if (!neigh_elem)
692 0 : EFAError("general elem is not of type EFAelement3D");
693 :
694 : // check if a general neighbor is an old crack tip element and will be split
695 : std::set<EFAElement *>::iterator sit;
696 130156 : sit = CrackTipElements.find(neigh_elem);
697 130156 : if (sit != CrackTipElements.end() && neigh_elem->numFragments() > 1)
698 : {
699 0 : for (unsigned int i = 0; i < neigh_elem->numFaces(); ++i)
700 : {
701 0 : std::set<EFANode *> neigh_face_nodes = neigh_elem->getFaceNodes(i);
702 0 : if (neigh_elem->numFaceNeighbors(i) == 2 &&
703 0 : Efa::numCommonElems(neigh_face_nodes, non_physical_nodes) > 0)
704 : {
705 : should_duplicate = true;
706 : break;
707 : }
708 : } // i
709 : }
710 : if (should_duplicate)
711 : break;
712 : } // eit
713 : }
714 4548 : } // IF only one fragment
715 4548 : return should_duplicate;
716 : }
717 :
718 : bool
719 3828 : EFAElement3D::shouldDuplicateForPhantomCorner()
720 : {
721 : // if a partial element will be split for a second time and it has two neighbor elements
722 : // sharing one phantom node with the aforementioned partial element, then the two neighbor
723 : // elements should be duplicated
724 : bool should_duplicate = false;
725 3828 : if (_fragments.size() == 1 && (!_crack_tip_split_element))
726 : {
727 19650 : for (unsigned int i = 0; i < _num_faces; ++i)
728 : {
729 16420 : std::set<EFANode *> phantom_nodes = getPhantomNodeOnFace(i);
730 16420 : if (phantom_nodes.size() > 0 && numFaceNeighbors(i) == 1)
731 : {
732 7752 : EFAElement3D * neighbor_elem = _face_neighbors[i][0];
733 7752 : if (neighbor_elem->numFragments() > 1) // neighbor will be split
734 : {
735 0 : for (unsigned int j = 0; j < neighbor_elem->numFaces(); ++j)
736 : {
737 0 : if (!neighbor_elem->getFace(j)->equivalent(_faces[i]) &&
738 0 : neighbor_elem->numFaceNeighbors(j) > 0)
739 : {
740 0 : std::set<EFANode *> neigh_phantom_nodes = neighbor_elem->getPhantomNodeOnFace(j);
741 0 : if (Efa::numCommonElems(phantom_nodes, neigh_phantom_nodes) > 0)
742 : {
743 : should_duplicate = true;
744 : break;
745 : }
746 : }
747 : } // j
748 : }
749 : }
750 : if (should_duplicate)
751 : break;
752 : } // i
753 : }
754 3828 : return should_duplicate;
755 : }
756 :
757 : bool
758 4548 : EFAElement3D::willCrackTipExtend(std::vector<unsigned int> & split_neighbors) const
759 : {
760 : // Determine whether the current element is a crack tip element for which the crack will
761 : // extend into the next element.
762 : // N.B. this is called at the beginning of createChildElements
763 : bool will_extend = false;
764 4548 : if (_fragments.size() == 1 && _crack_tip_split_element)
765 : {
766 2708 : for (unsigned int i = 0; i < _crack_tip_neighbors.size(); ++i)
767 : {
768 1390 : unsigned int neigh_idx = _crack_tip_neighbors[i]; // essentially a face_id
769 1390 : if (numFaceNeighbors(neigh_idx) != 1)
770 0 : EFAError("in will_crack_tip_extend() element ",
771 : _id,
772 : " has ",
773 : _face_neighbors[neigh_idx].size(),
774 : " neighbors on face ",
775 : neigh_idx);
776 :
777 1390 : EFAElement3D * neighbor_elem = _face_neighbors[neigh_idx][0];
778 1390 : if (neighbor_elem->numFragments() > 2)
779 : {
780 0 : EFAError("in will_crack_tip_extend() element ",
781 : neighbor_elem->id(),
782 : " has ",
783 : neighbor_elem->numFragments(),
784 : " fragments");
785 : }
786 1390 : else if (neighbor_elem->numFragments() == 2)
787 : {
788 720 : EFAFragment3D * neigh_frag1 = neighbor_elem->getFragment(0);
789 720 : EFAFragment3D * neigh_frag2 = neighbor_elem->getFragment(1);
790 720 : std::vector<EFANode *> neigh_cut_nodes = neigh_frag1->getCommonNodes(neigh_frag2);
791 : unsigned int counter = 0; // counter how many common nodes are contained by current face
792 3600 : for (unsigned int j = 0; j < neigh_cut_nodes.size(); ++j)
793 : {
794 2880 : if (_faces[neigh_idx]->containsNode(neigh_cut_nodes[j]))
795 1440 : counter += 1;
796 : }
797 720 : if (counter == 2)
798 : {
799 720 : split_neighbors.push_back(neigh_idx);
800 : will_extend = true;
801 : }
802 720 : }
803 : } // i
804 : }
805 4548 : return will_extend;
806 : }
807 :
808 : bool
809 54824 : EFAElement3D::isCrackTipElement() const
810 : {
811 54824 : return fragmentHasTipFaces();
812 : }
813 :
814 : unsigned int
815 6728 : EFAElement3D::getNumCuts() const
816 : {
817 : unsigned int num_cut_faces = 0;
818 41336 : for (unsigned int i = 0; i < _num_faces; ++i)
819 34608 : if (_faces[i]->hasIntersection())
820 0 : num_cut_faces += 1;
821 6728 : return num_cut_faces;
822 : }
823 :
824 : bool
825 23388 : EFAElement3D::isFinalCut() const
826 : {
827 : // if an element has been cut third times its fragment must have 3 interior faces
828 : // and at this point, we do not want it to be further cut
829 : bool cut_third = false;
830 23388 : if (_fragments.size() > 0)
831 : {
832 : unsigned int num_interior_faces = 0;
833 20074 : for (unsigned int i = 0; i < _fragments[0]->numFaces(); ++i)
834 : {
835 16990 : if (_fragments[0]->isFaceInterior(i))
836 2798 : num_interior_faces += 1;
837 : }
838 3084 : if (num_interior_faces == 3)
839 : cut_third = true;
840 : }
841 23388 : return cut_third;
842 : }
843 :
844 : void
845 30213 : EFAElement3D::updateFragments(const std::set<EFAElement *> & CrackTipElements,
846 : std::map<unsigned int, EFANode *> & EmbeddedNodes)
847 : {
848 : // combine the crack-tip faces in a fragment to a single intersected face
849 : std::set<EFAElement *>::iterator sit;
850 30213 : sit = CrackTipElements.find(this);
851 30213 : if (sit != CrackTipElements.end()) // curr_elem is a crack tip element
852 : {
853 370 : if (_fragments.size() == 1)
854 370 : _fragments[0]->combine_tip_faces();
855 : else
856 0 : EFAError("crack tip elem ", _id, " must have 1 fragment");
857 : }
858 :
859 : // remove the inappropriate embedded nodes on interior faces
860 : // (MUST DO THIS AFTER combine_tip_faces())
861 30213 : if (_fragments.size() == 1)
862 3948 : _fragments[0]->removeInvalidEmbeddedNodes(EmbeddedNodes);
863 :
864 : // for an element with no fragment, create one fragment identical to the element
865 30213 : if (_fragments.size() == 0)
866 26265 : _fragments.push_back(new EFAFragment3D(this, true, this));
867 30213 : if (_fragments.size() != 1)
868 0 : EFAError("Element ", _id, " must have 1 fragment at this point");
869 :
870 : // count fragment's cut faces
871 30213 : unsigned int num_cut_frag_faces = _fragments[0]->getNumCuts();
872 30213 : unsigned int num_frag_faces = _fragments[0]->numFaces();
873 30213 : if (num_cut_frag_faces > _fragments[0]->numFaces())
874 0 : EFAError("In element ", _id, " there are too many cut fragment faces");
875 :
876 : // leave the uncut frag as it is
877 30213 : if (num_cut_frag_faces == 0)
878 : {
879 28384 : if (!isPartial()) // delete the temp frag for an uncut elem
880 : {
881 24806 : delete _fragments[0];
882 24806 : _fragments.clear();
883 : }
884 28384 : return;
885 : }
886 :
887 : // split one fragment into one or two new fragments
888 1829 : std::vector<EFAFragment3D *> new_frags = _fragments[0]->split();
889 1829 : if (new_frags.size() == 1 || new_frags.size() == 2)
890 : {
891 1829 : delete _fragments[0]; // delete the old fragment
892 1829 : _fragments.clear();
893 4910 : for (unsigned int i = 0; i < new_frags.size(); ++i)
894 3081 : _fragments.push_back(new_frags[i]);
895 : }
896 : else
897 0 : EFAError("Number of fragments must be 1 or 2 at this point");
898 :
899 1829 : fragmentSanityCheck(num_frag_faces, num_cut_frag_faces);
900 1829 : }
901 :
902 : void
903 1829 : EFAElement3D::fragmentSanityCheck(unsigned int n_old_frag_faces, unsigned int n_old_frag_cuts) const
904 : {
905 1829 : unsigned int n_interior_nodes = numInteriorNodes();
906 1829 : if (n_interior_nodes > 0 && n_interior_nodes != 1)
907 0 : EFAError("After update_fragments this element has ", n_interior_nodes, " interior nodes");
908 :
909 1829 : if (n_old_frag_cuts == 0)
910 : {
911 0 : if (_fragments.size() != 1 || _fragments[0]->numFaces() != n_old_frag_faces)
912 0 : EFAError("Incorrect link size for element with 0 cuts");
913 : }
914 1829 : else if (fragmentHasTipFaces()) // crack tip case
915 : {
916 577 : if (_fragments.size() != 1 || _fragments[0]->numFaces() != n_old_frag_faces + n_old_frag_cuts)
917 0 : EFAError("Incorrect link size for element with crack-tip faces");
918 : }
919 : else // frag is thoroughly cut
920 : {
921 1252 : if (_fragments.size() != 2 || (_fragments[0]->numFaces() + _fragments[1]->numFaces()) !=
922 1252 : n_old_frag_faces + n_old_frag_cuts + 2)
923 0 : EFAError("Incorrect link size for element that has been completely cut");
924 : }
925 1829 : }
926 :
927 : void
928 6728 : EFAElement3D::restoreFragment(const EFAElement * const from_elem)
929 : {
930 6728 : const EFAElement3D * from_elem3d = dynamic_cast<const EFAElement3D *>(from_elem);
931 6728 : if (!from_elem3d)
932 0 : EFAError("from_elem is not of EFAelement3D type");
933 :
934 : // restore fragments
935 6728 : if (_fragments.size() != 0)
936 0 : EFAError("in restoreFragmentInfo elements must not have any pre-existing fragments");
937 13456 : for (unsigned int i = 0; i < from_elem3d->numFragments(); ++i)
938 6728 : _fragments.push_back(new EFAFragment3D(this, true, from_elem3d, i));
939 :
940 : // restore interior nodes
941 6728 : if (_interior_nodes.size() != 0)
942 0 : EFAError("in restoreFragmentInfo elements must not have any pre-exsiting interior nodes");
943 6728 : for (unsigned int i = 0; i < from_elem3d->_interior_nodes.size(); ++i)
944 0 : _interior_nodes.push_back(new EFAVolumeNode(*from_elem3d->_interior_nodes[i]));
945 :
946 : // restore face intersections
947 6728 : if (getNumCuts() != 0)
948 0 : EFAError("In restoreEdgeIntersection: edge cuts already exist in element ", _id);
949 41336 : for (unsigned int i = 0; i < _num_faces; ++i)
950 34608 : _faces[i]->copyIntersection(*from_elem3d->_faces[i]);
951 :
952 : // replace all local nodes with global nodes
953 74862 : for (unsigned int i = 0; i < from_elem3d->numNodes(); ++i)
954 : {
955 68134 : if (from_elem3d->_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
956 68134 : switchNode(
957 : _nodes[i], from_elem3d->_nodes[i], false); // EFAelement is not a child of any parent
958 : else
959 0 : EFAError("In restoreFragmentInfo all of from_elem's nodes must be local");
960 : }
961 6728 : }
962 :
963 : void
964 30213 : EFAElement3D::createChild(const std::set<EFAElement *> & CrackTipElements,
965 : std::map<unsigned int, EFAElement *> & Elements,
966 : std::map<unsigned int, EFAElement *> & newChildElements,
967 : std::vector<EFAElement *> & ChildElements,
968 : std::vector<EFAElement *> & ParentElements,
969 : std::map<unsigned int, EFANode *> & TempNodes)
970 : {
971 30213 : if (_children.size() != 0)
972 0 : EFAError("Element cannot have existing children in createChildElements");
973 :
974 30213 : if (_fragments.size() > 1 || shouldDuplicateForCrackTip(CrackTipElements))
975 : {
976 1579 : if (_fragments.size() > 2)
977 0 : EFAError("More than 2 fragments not yet supported");
978 :
979 : // set up the children
980 1579 : ParentElements.push_back(this);
981 4410 : for (unsigned int ichild = 0; ichild < _fragments.size(); ++ichild)
982 : {
983 : unsigned int new_elem_id;
984 2831 : if (newChildElements.size() == 0)
985 150 : new_elem_id = Efa::getNewID(Elements);
986 : else
987 2681 : new_elem_id = Efa::getNewID(newChildElements);
988 :
989 2831 : EFAElement3D * childElem = new EFAElement3D(new_elem_id, this->numNodes(), this->numFaces());
990 2831 : newChildElements.insert(std::make_pair(new_elem_id, childElem));
991 :
992 2831 : ChildElements.push_back(childElem);
993 2831 : childElem->setParent(this);
994 2831 : _children.push_back(childElem);
995 :
996 : std::vector<std::vector<EFANode *>> cut_plane_nodes;
997 18260 : for (unsigned int i = 0; i < this->getFragment(ichild)->numFaces(); ++i)
998 : {
999 15429 : if (this->getFragment(ichild)->isFaceInterior(i))
1000 : {
1001 2584 : EFAFace * face = this->getFragment(ichild)->getFace(i);
1002 : std::vector<EFANode *> node_line;
1003 11976 : for (unsigned int j = 0; j < face->numNodes(); ++j)
1004 9392 : node_line.push_back(face->getNode(j));
1005 2584 : cut_plane_nodes.push_back(node_line);
1006 2584 : }
1007 : }
1008 :
1009 : std::vector<EFAPoint> cut_plane_points;
1010 :
1011 2831 : EFAPoint orig(0.0, 0.0, 0.0);
1012 :
1013 : std::vector<EFAPoint> normal_v;
1014 :
1015 2831 : if (cut_plane_nodes.size())
1016 : {
1017 11976 : for (unsigned int i = 0; i < cut_plane_nodes[0].size(); ++i)
1018 : {
1019 : std::vector<EFANode *> master_nodes;
1020 : std::vector<double> master_weights;
1021 :
1022 9392 : this->getMasterInfo(cut_plane_nodes[0][i], master_nodes, master_weights);
1023 9392 : EFAPoint coor(0.0, 0.0, 0.0);
1024 28176 : for (unsigned int i = 0; i < master_nodes.size(); ++i)
1025 : {
1026 18784 : EFANode * local = this->createLocalNodeFromGlobalNode(master_nodes[i]);
1027 18784 : coor += _local_node_coor[local->id()] * master_weights[i];
1028 18784 : delete local;
1029 : }
1030 9392 : cut_plane_points.push_back(coor);
1031 9392 : }
1032 :
1033 11976 : for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
1034 9392 : orig += cut_plane_points[i];
1035 2584 : orig /= cut_plane_points.size();
1036 :
1037 11976 : for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
1038 : {
1039 9392 : unsigned int iplus1 = i < cut_plane_points.size() - 1 ? i + 1 : 0;
1040 9392 : unsigned int iminus1 = i == 0 ? cut_plane_points.size() - 1 : i - 1;
1041 9392 : EFAPoint ray1 = cut_plane_points[iminus1] - cut_plane_points[i];
1042 9392 : EFAPoint ray2 = cut_plane_points[i] - cut_plane_points[iplus1];
1043 18784 : normal_v.push_back(ray1.cross(ray2));
1044 : }
1045 :
1046 11976 : for (unsigned int i = 0; i < normal_v.size(); ++i)
1047 9392 : Xfem::normalizePoint(normal_v[i]);
1048 : }
1049 :
1050 : // get child element's nodes
1051 31839 : for (unsigned int j = 0; j < _num_nodes; ++j)
1052 : {
1053 29008 : EFAPoint p(0.0, 0.0, 0.0);
1054 29008 : p = _local_node_coor[j];
1055 29008 : EFAPoint origin_to_point = p - orig;
1056 :
1057 : bool node_outside_fragment = false;
1058 76432 : for (unsigned int k = 0; k < normal_v.size(); ++k)
1059 60688 : if (origin_to_point * normal_v[k] > Xfem::tol)
1060 : {
1061 : node_outside_fragment = true;
1062 : break;
1063 : }
1064 :
1065 29008 : if (_fragments.size() == 1 && !shouldDuplicateForCrackTip(CrackTipElements))
1066 0 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1067 29008 : else if (!node_outside_fragment)
1068 15744 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1069 : else // parent element's node is not in fragment
1070 : {
1071 13264 : unsigned int new_node_id = Efa::getNewID(TempNodes);
1072 13264 : EFANode * newNode = new EFANode(new_node_id, EFANode::N_CATEGORY_TEMP, _nodes[j]);
1073 13264 : TempNodes.insert(std::make_pair(new_node_id, newNode));
1074 13264 : childElem->setNode(j, newNode); // be a temp node
1075 : }
1076 : }
1077 :
1078 : // get child element's fragments
1079 2831 : EFAFragment3D * new_frag = new EFAFragment3D(childElem, true, this, ichild);
1080 2831 : childElem->_fragments.push_back(new_frag);
1081 :
1082 : // get child element's faces and set up adjacent faces
1083 2831 : childElem->createFaces();
1084 17257 : for (unsigned int j = 0; j < _num_faces; ++j)
1085 14426 : childElem->_faces[j]->copyIntersection(*_faces[j]);
1086 2831 : childElem->removePhantomEmbeddedNode(); // IMPORTANT
1087 :
1088 : // inherit old interior nodes
1089 2831 : for (unsigned int j = 0; j < _interior_nodes.size(); ++j)
1090 0 : childElem->_interior_nodes.push_back(new EFAVolumeNode(*_interior_nodes[j]));
1091 2831 : }
1092 : }
1093 : else // num_links == 1 || num_links == 0
1094 : {
1095 : // child is itself - but don't insert into the list of ChildElements!!!
1096 28634 : _children.push_back(this);
1097 : }
1098 30213 : }
1099 :
1100 : void
1101 2831 : EFAElement3D::removePhantomEmbeddedNode()
1102 : {
1103 : // remove the embedded nodes on faces that are outside the real domain
1104 2831 : if (_fragments.size() > 0)
1105 : {
1106 17257 : for (unsigned int i = 0; i < _num_faces; ++i)
1107 : {
1108 : // get emb nodes to be removed on edges
1109 : std::vector<EFANode *> nodes_to_delete;
1110 67010 : for (unsigned int j = 0; j < _faces[i]->numEdges(); ++j)
1111 : {
1112 : EFAEdge * edge = _faces[i]->getEdge(j);
1113 72380 : for (unsigned int k = 0; k < edge->numEmbeddedNodes(); ++k)
1114 : {
1115 19796 : if (!_fragments[0]->containsNode(edge->getEmbeddedNode(k)))
1116 0 : nodes_to_delete.push_back(edge->getEmbeddedNode(k));
1117 : } // k
1118 : } // j
1119 :
1120 : // get emb nodes to be removed in the face interior
1121 14426 : for (unsigned int j = 0; j < _faces[i]->numInteriorNodes(); ++j)
1122 : {
1123 0 : EFANode * face_node = _faces[i]->getInteriorNode(j)->getNode();
1124 0 : if (!_fragments[0]->containsNode(face_node))
1125 0 : nodes_to_delete.push_back(face_node);
1126 : } // j
1127 :
1128 : // remove all invalid embedded nodes
1129 14426 : for (unsigned int j = 0; j < nodes_to_delete.size(); ++j)
1130 0 : _faces[i]->removeEmbeddedNode(nodes_to_delete[j]);
1131 14426 : } // i
1132 : }
1133 2831 : }
1134 :
1135 : void
1136 2831 : EFAElement3D::connectNeighbors(std::map<unsigned int, EFANode *> & PermanentNodes,
1137 : std::map<unsigned int, EFANode *> & TempNodes,
1138 : std::map<EFANode *, std::set<EFAElement *>> & InverseConnectivityMap,
1139 : bool merge_phantom_faces)
1140 : {
1141 : // N.B. "this" must point to a child element that was just created
1142 2831 : if (!_parent)
1143 0 : EFAError("no parent element for child element ", _id, " in connect_neighbors");
1144 2831 : EFAElement3D * parent3d = dynamic_cast<EFAElement3D *>(_parent);
1145 2831 : if (!parent3d)
1146 0 : EFAError("cannot dynamic cast to parent3d in connect_neighbors");
1147 :
1148 : // First loop through edges and merge nodes with neighbors as appropriate
1149 17257 : for (unsigned int j = 0; j < _num_faces; ++j)
1150 : {
1151 25953 : for (unsigned int k = 0; k < parent3d->numFaceNeighbors(j); ++k)
1152 : {
1153 11527 : EFAElement3D * NeighborElem = parent3d->getFaceNeighbor(j, k);
1154 11527 : unsigned int neighbor_face_id = NeighborElem->getNeighborIndex(parent3d);
1155 :
1156 11527 : if (_faces[j]->hasIntersection())
1157 : {
1158 21011 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1159 : {
1160 : EFAElement3D * childOfNeighborElem =
1161 13732 : dynamic_cast<EFAElement3D *>(NeighborElem->getChild(l));
1162 13732 : if (!childOfNeighborElem)
1163 0 : EFAError("dynamic cast childOfNeighborElem fails");
1164 :
1165 : // Check to see if the nodes are already merged. There's nothing else to do in that case.
1166 13732 : EFAFace * neighborChildFace = childOfNeighborElem->getFace(neighbor_face_id);
1167 13732 : if (_faces[j]->equivalent(neighborChildFace))
1168 4079 : continue;
1169 :
1170 9653 : if (_fragments[0]->isConnected(childOfNeighborElem->getFragment(0)))
1171 : {
1172 15481 : for (unsigned int i = 0; i < _faces[j]->numNodes(); ++i)
1173 : {
1174 : unsigned int childNodeIndex = i;
1175 : unsigned int neighborChildNodeIndex =
1176 12076 : parent3d->getNeighborFaceNodeID(j, childNodeIndex, NeighborElem);
1177 :
1178 12076 : EFANode * childNode = _faces[j]->getNode(childNodeIndex);
1179 12076 : EFANode * childOfNeighborNode = neighborChildFace->getNode(neighborChildNodeIndex);
1180 12076 : mergeNodes(
1181 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1182 : } // i
1183 :
1184 9241 : for (unsigned int m = 0; m < _num_interior_face_nodes; ++m)
1185 : {
1186 : unsigned int childNodeIndex = m;
1187 : unsigned int neighborChildNodeIndex =
1188 5836 : parent3d->getNeighborFaceInteriorNodeID(j, childNodeIndex, NeighborElem);
1189 :
1190 5836 : EFANode * childNode = _faces[j]->getInteriorFaceNode(childNodeIndex);
1191 : EFANode * childOfNeighborNode =
1192 5836 : neighborChildFace->getInteriorFaceNode(neighborChildNodeIndex);
1193 5836 : mergeNodes(
1194 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1195 : } // m
1196 : }
1197 : } // l, loop over NeighborElem's children
1198 : }
1199 : else // No edge intersection -- optionally merge non-material nodes if they share a common
1200 : // parent
1201 : {
1202 4248 : if (merge_phantom_faces)
1203 : {
1204 8496 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1205 : {
1206 : EFAElement3D * childOfNeighborElem =
1207 4248 : dynamic_cast<EFAElement3D *>(NeighborElem->getChild(l));
1208 4248 : if (!childOfNeighborElem)
1209 0 : EFAError("dynamic cast childOfNeighborElem fails");
1210 :
1211 4248 : EFAFace * neighborChildFace = childOfNeighborElem->getFace(neighbor_face_id);
1212 4248 : if (!neighborChildFace
1213 4248 : ->hasIntersection()) // neighbor face must NOT have intersection either
1214 : {
1215 : // Check to see if the nodes are already merged. There's nothing else to do in that
1216 : // case.
1217 4248 : if (_faces[j]->equivalent(neighborChildFace))
1218 2588 : continue;
1219 :
1220 7852 : for (unsigned int i = 0; i < _faces[j]->numNodes(); ++i)
1221 : {
1222 : unsigned int childNodeIndex = i;
1223 : unsigned int neighborChildNodeIndex =
1224 6192 : parent3d->getNeighborFaceNodeID(j, childNodeIndex, NeighborElem);
1225 :
1226 6192 : EFANode * childNode = _faces[j]->getNode(childNodeIndex);
1227 6192 : EFANode * childOfNeighborNode = neighborChildFace->getNode(neighborChildNodeIndex);
1228 :
1229 11761 : if (childNode->parent() != nullptr &&
1230 5569 : childNode->parent() ==
1231 : childOfNeighborNode
1232 5569 : ->parent()) // non-material node and both come from same parent
1233 0 : mergeNodes(childNode,
1234 : childOfNeighborNode,
1235 : childOfNeighborElem,
1236 : PermanentNodes,
1237 : TempNodes);
1238 : } // i
1239 : }
1240 : } // loop over NeighborElem's children
1241 : } // if (merge_phantom_edges)
1242 : } // IF edge-j has_intersection()
1243 : } // k, loop over neighbors on edge j
1244 : } // j, loop over all faces
1245 :
1246 : // Now do a second loop through faces and convert remaining nodes to permanent nodes.
1247 : // If there is no neighbor on that face, also duplicate the embedded node if it exists
1248 31839 : for (unsigned int j = 0; j < _num_nodes; ++j)
1249 : {
1250 29008 : EFANode * childNode = _nodes[j];
1251 29008 : if (childNode->category() == EFANode::N_CATEGORY_TEMP)
1252 : {
1253 : // if current child element does not have siblings, and if current temp node is a lone one
1254 : // this temp node should be merged back to its parent permanent node. Otherwise we would have
1255 : // permanent nodes that are not connected to any element
1256 1620 : std::set<EFAElement *> patch_elems = InverseConnectivityMap[childNode->parent()];
1257 1620 : if (parent3d->numFragments() == 1 && patch_elems.size() == 1)
1258 0 : switchNode(childNode->parent(), childNode, false);
1259 : else
1260 : {
1261 1620 : unsigned int new_node_id = Efa::getNewID(PermanentNodes);
1262 : EFANode * newNode =
1263 1620 : new EFANode(new_node_id, EFANode::N_CATEGORY_PERMANENT, childNode->parent());
1264 1620 : PermanentNodes.insert(std::make_pair(new_node_id, newNode));
1265 1620 : switchNode(newNode, childNode, false);
1266 : }
1267 1620 : if (!Efa::deleteFromMap(TempNodes, childNode))
1268 0 : EFAError(
1269 : "Attempted to delete node: ", childNode->id(), " from TempNodes, but couldn't find it");
1270 : }
1271 : }
1272 2831 : }
1273 :
1274 : void
1275 0 : EFAElement3D::printElement(std::ostream & ostream) const
1276 : {
1277 : // first line: all elem faces
1278 : ostream << std::setw(5);
1279 0 : ostream << _id << "| ";
1280 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1281 : {
1282 0 : for (unsigned int k = 0; k < _faces[j]->numNodes(); ++k)
1283 0 : ostream << std::setw(5) << _faces[j]->getNode(k)->idCatString();
1284 0 : ostream << " | ";
1285 : }
1286 : ostream << std::endl;
1287 :
1288 : // second line: emb nodes in all faces + neighbor of each face
1289 : ostream << std::setw(5);
1290 : ostream << "embed"
1291 0 : << "| ";
1292 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1293 : {
1294 0 : for (unsigned int k = 0; k < _faces[j]->numEdges(); ++k)
1295 : {
1296 : ostream << std::setw(4);
1297 0 : if (_faces[j]->getEdge(k)->hasIntersection())
1298 : {
1299 0 : if (_faces[j]->getEdge(k)->numEmbeddedNodes() > 1)
1300 : {
1301 0 : ostream << "[";
1302 0 : for (unsigned int l = 0; l < _faces[j]->getEdge(k)->numEmbeddedNodes(); ++l)
1303 : {
1304 0 : ostream << _faces[j]->getEdge(k)->getEmbeddedNode(l)->id() << " ";
1305 0 : if (l == _faces[j]->getEdge(k)->numEmbeddedNodes() - 1)
1306 0 : ostream << "]";
1307 : else
1308 0 : ostream << " ";
1309 : } // l
1310 : }
1311 : else
1312 0 : ostream << _faces[j]->getEdge(k)->getEmbeddedNode(0)->id() << " ";
1313 : }
1314 : else
1315 0 : ostream << " -- ";
1316 : } // k
1317 0 : ostream << " | ";
1318 : } // j
1319 : ostream << std::endl;
1320 :
1321 : // third line: neighbors
1322 : ostream << std::setw(5);
1323 : ostream << "neigh"
1324 0 : << "| ";
1325 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1326 : {
1327 : ostream << std::setw(4);
1328 0 : if (numFaceNeighbors(j) > 1)
1329 : {
1330 0 : ostream << "[";
1331 0 : for (unsigned int k = 0; k < numFaceNeighbors(j); ++k)
1332 : {
1333 0 : ostream << getFaceNeighbor(j, k)->id();
1334 0 : if (k == numFaceNeighbors(j) - 1)
1335 0 : ostream << "]";
1336 : else
1337 0 : ostream << " ";
1338 : }
1339 : }
1340 : else
1341 : {
1342 0 : if (numFaceNeighbors(j) == 1)
1343 0 : ostream << getFaceNeighbor(j, 0)->id() << " ";
1344 : else
1345 0 : ostream << " -- ";
1346 : }
1347 : }
1348 : ostream << std::endl;
1349 :
1350 : // fourth line: fragments
1351 0 : for (unsigned int j = 0; j < _fragments.size(); ++j)
1352 : {
1353 : ostream << std::setw(4);
1354 0 : ostream << "frag" << j << "| ";
1355 0 : for (unsigned int k = 0; k < _fragments[j]->numFaces(); ++k)
1356 : {
1357 0 : for (unsigned int l = 0; l < _fragments[j]->getFace(k)->numNodes(); ++l)
1358 0 : ostream << std::setw(5) << _fragments[j]->getFace(k)->getNode(l)->idCatString();
1359 0 : ostream << " | ";
1360 : }
1361 : ostream << std::endl;
1362 : }
1363 : ostream << std::endl;
1364 0 : }
1365 :
1366 : EFAFragment3D *
1367 13669912 : EFAElement3D::getFragment(unsigned int frag_id) const
1368 : {
1369 13669912 : if (frag_id < _fragments.size())
1370 13669912 : return _fragments[frag_id];
1371 : else
1372 0 : EFAError("frag_id out of bounds");
1373 : }
1374 :
1375 : std::set<EFANode *>
1376 0 : EFAElement3D::getFaceNodes(unsigned int face_id) const
1377 : {
1378 : std::set<EFANode *> face_nodes;
1379 0 : for (unsigned int i = 0; i < _faces[face_id]->numNodes(); ++i)
1380 0 : face_nodes.insert(_faces[face_id]->getNode(i));
1381 0 : return face_nodes;
1382 : }
1383 :
1384 : bool
1385 0 : EFAElement3D::getFaceNodeParametricCoordinates(EFANode * node, std::vector<double> & xi_3d) const
1386 : {
1387 : // get the parametric coords of a node in an element face
1388 : unsigned int face_id = std::numeric_limits<unsigned int>::max();
1389 : bool face_found = false;
1390 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1391 : {
1392 0 : if (_faces[i]->containsNode(node))
1393 : {
1394 : face_id = i;
1395 : face_found = true;
1396 : break;
1397 : }
1398 : }
1399 0 : if (face_found)
1400 : {
1401 0 : std::vector<double> xi_2d(2, 0.0);
1402 0 : if (_faces[face_id]->getFaceNodeParametricCoords(node, xi_2d))
1403 0 : mapParametricCoordinateFrom2DTo3D(face_id, xi_2d, xi_3d);
1404 : else
1405 0 : EFAError("failed to get the 2D para coords on the face");
1406 0 : }
1407 0 : return face_found;
1408 : }
1409 :
1410 : EFAVolumeNode *
1411 0 : EFAElement3D::getInteriorNode(unsigned int interior_node_id) const
1412 : {
1413 0 : if (interior_node_id < _interior_nodes.size())
1414 0 : return _interior_nodes[interior_node_id];
1415 : else
1416 0 : EFAError("interior_node_id out of bounds");
1417 : }
1418 :
1419 : void
1420 0 : EFAElement3D::removeEmbeddedNode(EFANode * emb_node, bool remove_for_neighbor)
1421 : {
1422 0 : for (unsigned int i = 0; i < _fragments.size(); ++i)
1423 0 : _fragments[i]->removeEmbeddedNode(emb_node);
1424 :
1425 0 : for (unsigned int i = 0; i < _faces.size(); ++i)
1426 0 : _faces[i]->removeEmbeddedNode(emb_node);
1427 :
1428 0 : if (remove_for_neighbor)
1429 : {
1430 0 : for (unsigned int i = 0; i < numFaces(); ++i)
1431 0 : for (unsigned int j = 0; j < numFaceNeighbors(i); ++j)
1432 0 : getFaceNeighbor(i, j)->removeEmbeddedNode(emb_node, false);
1433 : }
1434 0 : }
1435 :
1436 : unsigned int
1437 22767201 : EFAElement3D::numFaces() const
1438 : {
1439 22767201 : return _faces.size();
1440 : }
1441 :
1442 : void
1443 0 : EFAElement3D::setFace(unsigned int face_id, EFAFace * face)
1444 : {
1445 0 : _faces[face_id] = face;
1446 0 : }
1447 :
1448 : void
1449 47798 : EFAElement3D::createFaces()
1450 : {
1451 : // create element faces based on existing element nodes
1452 47798 : int hex_local_node_indices[6][4] = {
1453 : {0, 3, 2, 1}, {0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}};
1454 47798 : int tet_local_node_indices[4][3] = {{0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}};
1455 :
1456 47798 : int hex_interior_face_node_indices[6][5] = {{8, 9, 10, 11, 20},
1457 : {8, 13, 16, 12, 21},
1458 : {9, 14, 17, 13, 22},
1459 : {10, 14, 18, 15, 23},
1460 : {11, 15, 19, 12, 24},
1461 : {16, 17, 18, 19, 25}};
1462 47798 : int tet_interior_face_node_indices[4][4] = {
1463 : {4, 5, 6, 10}, {4, 7, 8, 11}, {5, 8, 9, 12}, {6, 7, 9, 13}};
1464 :
1465 47798 : _faces = std::vector<EFAFace *>(_num_faces, nullptr);
1466 47798 : if (_num_nodes == 8 || _num_nodes == 20 || _num_nodes == 27)
1467 : {
1468 26430 : if (_num_faces != 6)
1469 0 : EFAError("num_faces of hexes must be 6");
1470 185010 : for (unsigned int i = 0; i < _num_faces; ++i)
1471 : {
1472 158580 : _faces[i] = new EFAFace(4, _num_interior_face_nodes);
1473 792900 : for (unsigned int j = 0; j < 4; ++j)
1474 634320 : _faces[i]->setNode(j, _nodes[hex_local_node_indices[i][j]]);
1475 158580 : _faces[i]->createEdges();
1476 179532 : for (unsigned int k = 0; k < _num_interior_face_nodes; ++k)
1477 20952 : _faces[i]->setInteriorFaceNode(k, _nodes[hex_interior_face_node_indices[i][k]]);
1478 : }
1479 : }
1480 : else if (_num_nodes == 4 || _num_nodes == 10 || _num_nodes == 14)
1481 : {
1482 21368 : if (_num_faces != 4)
1483 0 : EFAError("num_faces of tets must be 4");
1484 106840 : for (unsigned int i = 0; i < _num_faces; ++i)
1485 : {
1486 85472 : _faces[i] = new EFAFace(3, _num_interior_face_nodes);
1487 341888 : for (unsigned int j = 0; j < 3; ++j)
1488 256416 : _faces[i]->setNode(j, _nodes[tet_local_node_indices[i][j]]);
1489 85472 : _faces[i]->createEdges();
1490 384624 : for (unsigned int k = 0; k < _num_interior_face_nodes; ++k)
1491 299152 : _faces[i]->setInteriorFaceNode(k, _nodes[tet_interior_face_node_indices[i][k]]);
1492 : }
1493 : }
1494 : else
1495 0 : EFAError("unknown 3D element type in createFaces()");
1496 :
1497 291850 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
1498 : {
1499 244052 : _face_edge_neighbors[face_iter].resize(_faces[face_iter]->numEdges());
1500 1134788 : for (unsigned int edge_iter = 0; edge_iter < _faces[face_iter]->numEdges(); ++edge_iter)
1501 890736 : _face_edge_neighbors[face_iter][edge_iter].clear();
1502 : }
1503 :
1504 : // create element face connectivity array
1505 47798 : findFacesAdjacentToFaces(); // IMPORTANT
1506 47798 : }
1507 :
1508 : EFAFace *
1509 21258453 : EFAElement3D::getFace(unsigned int face_id) const
1510 : {
1511 21258453 : return _faces[face_id];
1512 : }
1513 :
1514 : unsigned int
1515 286225 : EFAElement3D::getFaceID(EFAFace * face) const
1516 : {
1517 : bool found_face_id = false;
1518 : unsigned int face_id;
1519 931687 : for (unsigned int iface = 0; iface < _num_faces; ++iface)
1520 : {
1521 931687 : if (_faces[iface]->equivalent(face))
1522 : {
1523 : face_id = iface;
1524 : found_face_id = true;
1525 : break;
1526 : }
1527 : }
1528 286225 : if (!found_face_id)
1529 0 : EFAError("input face not found in get_face_id()");
1530 286225 : return face_id;
1531 : }
1532 :
1533 : std::vector<unsigned int>
1534 2096348 : EFAElement3D::getCommonFaceID(const EFAElement3D * other_elem) const
1535 : {
1536 : std::vector<unsigned int> face_id;
1537 11565572 : for (unsigned int i = 0; i < _num_faces; ++i)
1538 52782190 : for (unsigned int j = 0; j < other_elem->_num_faces; ++j)
1539 43665718 : if (_faces[i]->equivalent(other_elem->_faces[j]))
1540 : {
1541 352752 : face_id.push_back(i);
1542 : break;
1543 : }
1544 :
1545 2096348 : return face_id;
1546 0 : }
1547 :
1548 : bool
1549 1382614 : EFAElement3D::getCommonEdgeID(const EFAElement3D * other_elem,
1550 : std::vector<std::pair<unsigned int, unsigned int>> & common_ids) const
1551 : {
1552 : // all edges of the other element
1553 : std::set<std::pair<EFANode *, EFANode *>> other_edges;
1554 7309042 : for (const auto k : index_range(other_elem->_faces))
1555 : {
1556 5926428 : const auto & face = *other_elem->_faces[k];
1557 24893628 : for (const auto l : make_range(other_elem->_faces[k]->numEdges()))
1558 : {
1559 : const auto & edge = *face.getEdge(l);
1560 18967200 : other_edges.insert(edge.getSortedNodes());
1561 : }
1562 : }
1563 :
1564 : // loop over all edges of this element
1565 1382614 : common_ids.clear();
1566 7309042 : for (const auto i : index_range(_faces))
1567 24893628 : for (const auto j : make_range(_faces[i]->numEdges()))
1568 : {
1569 18967200 : const auto & edge = *_faces[i]->getEdge(j);
1570 18967200 : const auto edge_nodes = edge.getSortedNodes();
1571 :
1572 : // is this edge contained in the other element?
1573 18967200 : if (edge.isEmbeddedPermanent() || other_edges.count(edge_nodes) == 0)
1574 18244408 : continue;
1575 :
1576 722792 : common_ids.emplace_back(i, j);
1577 : }
1578 :
1579 1382614 : return common_ids.size() > 0;
1580 : }
1581 :
1582 : unsigned int
1583 18268 : EFAElement3D::getNeighborFaceNodeID(unsigned int face_id,
1584 : unsigned int node_id,
1585 : EFAElement3D * neighbor_elem) const
1586 : {
1587 : // get the corresponding node_id on the corresponding face of neighbor_elem
1588 : bool found_id = false;
1589 : unsigned int neigh_face_node_id;
1590 18268 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1591 18268 : if (common_face_id == face_id)
1592 : {
1593 18268 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1594 18268 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1595 42682 : for (unsigned int i = 0; i < neigh_face->numNodes(); ++i)
1596 : {
1597 42682 : if (_faces[face_id]->getNode(node_id) == neigh_face->getNode(i))
1598 : {
1599 : neigh_face_node_id = i;
1600 : found_id = true;
1601 : break;
1602 : }
1603 : }
1604 : }
1605 : else
1606 0 : EFAError("getNeighborFaceNodeID: neighbor_elem is not a neighbor on face_id");
1607 18268 : if (!found_id)
1608 0 : EFAError("getNeighborFaceNodeID: could not find neighbor face node id");
1609 18268 : return neigh_face_node_id;
1610 : }
1611 :
1612 : unsigned int
1613 5836 : EFAElement3D::getNeighborFaceInteriorNodeID(unsigned int face_id,
1614 : unsigned int node_id,
1615 : EFAElement3D * neighbor_elem) const
1616 : {
1617 : // get the corresponding node_id on the corresponding face of neighbor_elem
1618 : bool found_id = false;
1619 : unsigned int neigh_face_node_id;
1620 5836 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1621 5836 : if (common_face_id == face_id)
1622 : {
1623 5836 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1624 5836 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1625 :
1626 13552 : for (unsigned int i = 0; i < _num_interior_face_nodes; ++i)
1627 : {
1628 13552 : if (_faces[face_id]->getInteriorFaceNode(node_id) == neigh_face->getInteriorFaceNode(i))
1629 : {
1630 : neigh_face_node_id = i;
1631 : found_id = true;
1632 : break;
1633 : }
1634 : }
1635 : }
1636 : else
1637 0 : EFAError("getNeighborFaceNodeID: neighbor_elem is not a neighbor on face_id");
1638 5836 : if (!found_id)
1639 0 : EFAError("getNeighborFaceNodeID: could not find neighbor face node id");
1640 5836 : return neigh_face_node_id;
1641 : }
1642 :
1643 : unsigned int
1644 42924 : EFAElement3D::getNeighborFaceEdgeID(unsigned int face_id,
1645 : unsigned int edge_id,
1646 : EFAElement3D * neighbor_elem) const
1647 : {
1648 : // get the corresponding edge_id on the corresponding face of neighbor_elem
1649 : bool found_id = false;
1650 : unsigned int neigh_face_edge_id;
1651 42924 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1652 42924 : if (common_face_id == face_id)
1653 : {
1654 42924 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1655 42924 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1656 90988 : for (unsigned int i = 0; i < neigh_face->numEdges(); ++i)
1657 : {
1658 90988 : if (_faces[face_id]->getEdge(edge_id)->equivalent(*neigh_face->getEdge(i)))
1659 : {
1660 : neigh_face_edge_id = i;
1661 : found_id = true;
1662 : break;
1663 : }
1664 : }
1665 : }
1666 : else
1667 0 : EFAError("getNeighborFaceEdgeID: neighbor_elem is not a neighbor on face_id");
1668 42924 : if (!found_id)
1669 0 : EFAError("getNeighborFaceEdgeID: could not find neighbor face edge id");
1670 42924 : return neigh_face_edge_id;
1671 : }
1672 :
1673 : void
1674 50478 : EFAElement3D::findFacesAdjacentToFaces()
1675 : {
1676 50478 : _faces_adjacent_to_faces.clear();
1677 308050 : for (unsigned int i = 0; i < _faces.size(); ++i)
1678 : {
1679 257572 : std::vector<EFAFace *> face_adjacents(_faces[i]->numEdges(), nullptr);
1680 1621820 : for (unsigned int j = 0; j < _faces.size(); ++j)
1681 : {
1682 1364248 : if (_faces[j] != _faces[i] && _faces[i]->isAdjacent(_faces[j]))
1683 : {
1684 939696 : unsigned int adj_edge = _faces[i]->adjacentCommonEdge(_faces[j]);
1685 939696 : face_adjacents[adj_edge] = _faces[j];
1686 : }
1687 : }
1688 257572 : _faces_adjacent_to_faces.push_back(face_adjacents);
1689 257572 : }
1690 50478 : }
1691 :
1692 : EFAFace *
1693 109433 : EFAElement3D::getAdjacentFace(unsigned int face_id, unsigned int edge_id) const
1694 : {
1695 109433 : return _faces_adjacent_to_faces[face_id][edge_id];
1696 : }
1697 :
1698 : EFAFace *
1699 149525 : EFAElement3D::getFragmentFace(unsigned int frag_id, unsigned int face_id) const
1700 : {
1701 149525 : if (frag_id < _fragments.size())
1702 149525 : return _fragments[frag_id]->getFace(face_id);
1703 : else
1704 0 : EFAError("frag_id out of bounds in getFragmentFace");
1705 : }
1706 :
1707 : std::set<EFANode *>
1708 16420 : EFAElement3D::getPhantomNodeOnFace(unsigned int face_id) const
1709 : {
1710 : std::set<EFANode *> phantom_nodes;
1711 16420 : if (_fragments.size() > 0)
1712 : {
1713 76180 : for (unsigned int j = 0; j < _faces[face_id]->numNodes(); ++j) // loop ove 2 edge nodes
1714 : {
1715 : bool node_in_frag = false;
1716 86160 : for (unsigned int k = 0; k < _fragments.size(); ++k)
1717 : {
1718 59760 : if (_fragments[k]->containsNode(_faces[face_id]->getNode(j)))
1719 : {
1720 : node_in_frag = true;
1721 : break;
1722 : }
1723 : }
1724 59760 : if (!node_in_frag)
1725 26400 : phantom_nodes.insert(_faces[face_id]->getNode(j));
1726 : }
1727 : }
1728 16420 : return phantom_nodes;
1729 : }
1730 :
1731 : bool
1732 20284 : EFAElement3D::getFragmentFaceID(unsigned int elem_face_id, unsigned int & frag_face_id) const
1733 : {
1734 : // find the fragment face that is contained by given element edge
1735 : // N.B. if the elem edge contains two frag edges, this method will only return
1736 : // the first frag edge ID
1737 : bool frag_face_found = false;
1738 20284 : frag_face_id = std::numeric_limits<unsigned int>::max();
1739 20284 : if (_fragments.size() == 1)
1740 : {
1741 1680 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1742 : {
1743 1680 : if (_faces[elem_face_id]->containsFace(_fragments[0]->getFace(j)))
1744 : {
1745 480 : frag_face_id = j;
1746 : frag_face_found = true;
1747 480 : break;
1748 : }
1749 : }
1750 : }
1751 20284 : return frag_face_found;
1752 : }
1753 :
1754 : bool
1755 10222 : EFAElement3D::getFragmentFaceEdgeID(unsigned int ElemFaceID,
1756 : unsigned int ElemFaceEdgeID,
1757 : unsigned int & FragFaceID,
1758 : unsigned int & FragFaceEdgeID) const
1759 : {
1760 : // Purpose: given an edge of an elem face, find which frag face's edge it contains
1761 : bool frag_edge_found = false;
1762 10222 : FragFaceID = FragFaceEdgeID = std::numeric_limits<unsigned int>::max();
1763 10222 : if (getFragmentFaceID(ElemFaceID, FragFaceID))
1764 : {
1765 320 : EFAEdge * elem_edge = _faces[ElemFaceID]->getEdge(ElemFaceEdgeID);
1766 320 : EFAFace * frag_face = getFragmentFace(0, FragFaceID);
1767 800 : for (unsigned int i = 0; i < frag_face->numEdges(); ++i)
1768 : {
1769 800 : if (elem_edge->containsEdge(*frag_face->getEdge(i)))
1770 : {
1771 320 : FragFaceEdgeID = i;
1772 : frag_edge_found = true;
1773 320 : break;
1774 : }
1775 : }
1776 : }
1777 10222 : return frag_edge_found;
1778 : }
1779 :
1780 : bool
1781 10062 : EFAElement3D::isPhysicalEdgeCut(unsigned int ElemFaceID,
1782 : unsigned int ElemFaceEdgeID,
1783 : double position) const
1784 : {
1785 10062 : unsigned int FragFaceID, FragFaceEdgeID = std::numeric_limits<unsigned int>::max();
1786 : bool is_in_real = false;
1787 10062 : if (_fragments.size() == 0)
1788 : {
1789 : is_in_real = true;
1790 : }
1791 160 : else if (getFragmentFaceEdgeID(ElemFaceID, ElemFaceEdgeID, FragFaceID, FragFaceEdgeID))
1792 : {
1793 160 : EFAEdge * elem_edge = _faces[ElemFaceID]->getEdge(ElemFaceEdgeID);
1794 160 : EFAEdge * frag_edge = getFragmentFace(0, FragFaceID)->getEdge(FragFaceEdgeID);
1795 : double xi[2] = {-1.0, -1.0}; // relative coords of two frag edge nodes
1796 160 : xi[0] = elem_edge->distanceFromNode1(frag_edge->getNode(0));
1797 160 : xi[1] = elem_edge->distanceFromNode1(frag_edge->getNode(1));
1798 160 : if ((position - xi[0]) * (position - xi[1]) <
1799 : 0.0) // the cut to be added is within the real part of the edge
1800 : is_in_real = true;
1801 : }
1802 10062 : return is_in_real;
1803 : }
1804 :
1805 : bool
1806 17958 : EFAElement3D::isFacePhantom(unsigned int face_id) const
1807 : {
1808 : bool is_phantom = false;
1809 17958 : if (_fragments.size() > 0)
1810 : {
1811 : bool contains_frag_face = false;
1812 13610 : for (unsigned int i = 0; i < _fragments.size(); ++i)
1813 : {
1814 38142 : for (unsigned int j = 0; j < _fragments[i]->numFaces(); ++j)
1815 : {
1816 38142 : if (_faces[face_id]->containsFace(_fragments[i]->getFace(j)))
1817 : {
1818 : contains_frag_face = true;
1819 : break;
1820 : }
1821 : }
1822 13610 : if (contains_frag_face)
1823 : break;
1824 : }
1825 13610 : if (!contains_frag_face)
1826 : is_phantom = true;
1827 : }
1828 17958 : return is_phantom;
1829 : }
1830 :
1831 : unsigned int
1832 1567541 : EFAElement3D::numFaceNeighbors(unsigned int face_id) const
1833 : {
1834 1567541 : return _face_neighbors[face_id].size();
1835 : }
1836 :
1837 : unsigned int
1838 76522 : EFAElement3D::numEdgeNeighbors(unsigned int face_id, unsigned int edge_id) const
1839 : {
1840 76522 : return _face_edge_neighbors[face_id][edge_id].size();
1841 : }
1842 :
1843 : EFAElement3D *
1844 755359 : EFAElement3D::getFaceNeighbor(unsigned int face_id, unsigned int neighbor_id) const
1845 : {
1846 755359 : if (_face_neighbors[face_id][0] != nullptr && neighbor_id < _face_neighbors[face_id].size())
1847 755359 : return _face_neighbors[face_id][neighbor_id];
1848 : else
1849 0 : EFAError("edge neighbor does not exist");
1850 : }
1851 :
1852 : EFAElement3D *
1853 40204 : EFAElement3D::getEdgeNeighbor(unsigned int face_id,
1854 : unsigned int edge_id,
1855 : unsigned int neighbor_id) const
1856 : {
1857 40204 : if (neighbor_id < _face_edge_neighbors[face_id][edge_id].size())
1858 40204 : return _face_edge_neighbors[face_id][edge_id][neighbor_id];
1859 : else
1860 0 : EFAError("edge neighbor does not exist");
1861 : }
1862 :
1863 : bool
1864 56653 : EFAElement3D::fragmentHasTipFaces() const
1865 : {
1866 : bool has_tip_faces = false;
1867 56653 : if (_fragments.size() == 1)
1868 : {
1869 94714 : for (unsigned int i = 0; i < _num_faces; ++i)
1870 : {
1871 : unsigned int num_frag_faces = 0; // count how many fragment edges this element edge contains
1872 81638 : if (_faces[i]->hasIntersection())
1873 : {
1874 343872 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1875 : {
1876 291608 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1877 56432 : num_frag_faces += 1;
1878 : } // j
1879 52264 : if (num_frag_faces == 2)
1880 : {
1881 : has_tip_faces = true;
1882 : break;
1883 : }
1884 : }
1885 : } // i
1886 : }
1887 56653 : return has_tip_faces;
1888 : }
1889 :
1890 : std::vector<unsigned int>
1891 0 : EFAElement3D::getTipFaceIDs() const
1892 : {
1893 : // if this element is a crack tip element, returns the crack tip faces' ID
1894 : std::vector<unsigned int> tip_face_id;
1895 0 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1896 : {
1897 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1898 : {
1899 : unsigned int num_frag_faces = 0; // count how many fragment faces this element edge contains
1900 0 : if (_faces[i]->hasIntersection())
1901 : {
1902 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1903 : {
1904 0 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1905 0 : num_frag_faces += 1;
1906 : } // j
1907 0 : if (num_frag_faces == 2) // element face contains two fragment edges
1908 0 : tip_face_id.push_back(i);
1909 : }
1910 : } // i
1911 : }
1912 0 : return tip_face_id;
1913 0 : }
1914 :
1915 : std::set<EFANode *>
1916 0 : EFAElement3D::getTipEmbeddedNodes() const
1917 : {
1918 : // if this element is a crack tip element, returns the crack tip edge's ID
1919 : std::set<EFANode *> tip_emb;
1920 0 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1921 : {
1922 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1923 : {
1924 : std::vector<EFAFace *> frag_faces; // count how many fragment edges this element edge contains
1925 0 : if (_faces[i]->hasIntersection())
1926 : {
1927 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1928 : {
1929 0 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1930 0 : frag_faces.push_back(_fragments[0]->getFace(j));
1931 : } // j
1932 0 : if (frag_faces.size() == 2) // element edge contains two fragment edges
1933 : {
1934 0 : unsigned int edge_id = frag_faces[0]->adjacentCommonEdge(frag_faces[1]);
1935 0 : tip_emb.insert(frag_faces[0]->getEdge(edge_id)->getNode(0));
1936 0 : tip_emb.insert(frag_faces[0]->getEdge(edge_id)->getNode(1));
1937 : }
1938 : }
1939 0 : } // i
1940 : }
1941 0 : return tip_emb;
1942 : }
1943 :
1944 : bool
1945 10062 : EFAElement3D::faceContainsTip(unsigned int face_id) const
1946 : {
1947 : bool contain_tip = false;
1948 10062 : if (_fragments.size() == 1)
1949 : {
1950 : unsigned int num_frag_faces = 0; // count how many fragment faces this element face contains
1951 160 : if (_faces[face_id]->hasIntersection())
1952 : {
1953 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1954 : {
1955 0 : if (_faces[face_id]->containsFace(_fragments[0]->getFace(j)))
1956 0 : num_frag_faces += 1;
1957 : } // j
1958 0 : if (num_frag_faces == 2)
1959 : contain_tip = true;
1960 : }
1961 : }
1962 10062 : return contain_tip;
1963 : }
1964 :
1965 : bool
1966 10062 : EFAElement3D::fragmentFaceAlreadyCut(unsigned int ElemFaceID) const
1967 : {
1968 : // when marking cuts, check if the corresponding frag face already has been cut
1969 : bool has_cut = false;
1970 10062 : if (faceContainsTip(ElemFaceID))
1971 : has_cut = true;
1972 : else
1973 : {
1974 10062 : unsigned int FragFaceID = std::numeric_limits<unsigned int>::max();
1975 10062 : if (getFragmentFaceID(ElemFaceID, FragFaceID))
1976 : {
1977 160 : EFAFace * frag_face = getFragmentFace(0, FragFaceID);
1978 160 : if (frag_face->hasIntersection())
1979 : has_cut = true;
1980 : }
1981 : }
1982 10062 : return has_cut;
1983 : }
1984 :
1985 : void
1986 109433 : EFAElement3D::addFaceEdgeCut(unsigned int face_id,
1987 : unsigned int edge_id,
1988 : double position,
1989 : EFANode * embedded_node,
1990 : std::map<unsigned int, EFANode *> & EmbeddedNodes,
1991 : bool add_to_neighbor,
1992 : bool add_to_adjacent)
1993 : {
1994 : // Purpose: add intersection on Edge edge_id of Face face_id
1995 109433 : EFANode * local_embedded = nullptr;
1996 109433 : EFAEdge * cut_edge = _faces[face_id]->getEdge(edge_id);
1997 109433 : EFANode * edge_node1 = cut_edge->getNode(0);
1998 109433 : if (embedded_node) // use the existing embedded node if it was passed in
1999 73115 : local_embedded = embedded_node;
2000 :
2001 : // get adjacent face info
2002 109433 : EFAFace * adj_face = getAdjacentFace(face_id, edge_id);
2003 109433 : unsigned int adj_face_id = getFaceID(adj_face);
2004 109433 : unsigned int adj_edge_id = adj_face->adjacentCommonEdge(_faces[face_id]);
2005 :
2006 : // check if cut has already been added to this face edge
2007 : bool cut_exist = false;
2008 :
2009 109433 : if (cut_edge->hasIntersectionAtPosition(position, edge_node1))
2010 : {
2011 99371 : unsigned int emb_id = cut_edge->getEmbeddedNodeIndex(position, edge_node1);
2012 99371 : EFANode * old_emb = cut_edge->getEmbeddedNode(emb_id);
2013 99371 : if (embedded_node && embedded_node != old_emb)
2014 0 : EFAError("Attempting to add edge intersection when one already exists with different node.",
2015 : " elem: ",
2016 : _id,
2017 : " edge: ",
2018 : edge_id,
2019 : " position: ",
2020 : position);
2021 99371 : local_embedded = old_emb;
2022 : cut_exist = true;
2023 : }
2024 :
2025 20124 : if (!cut_exist && !fragmentFaceAlreadyCut(face_id) &&
2026 10062 : isPhysicalEdgeCut(face_id, edge_id, position))
2027 : {
2028 : // check if cut has already been added to the neighbor edges
2029 10062 : checkNeighborFaceCut(face_id, edge_id, position, edge_node1, embedded_node, local_embedded);
2030 10062 : checkNeighborFaceCut(
2031 : adj_face_id, adj_edge_id, position, edge_node1, embedded_node, local_embedded);
2032 :
2033 10062 : if (!local_embedded) // need to create new embedded node
2034 : {
2035 1850 : unsigned int new_node_id = Efa::getNewID(EmbeddedNodes);
2036 1850 : local_embedded = new EFANode(new_node_id, EFANode::N_CATEGORY_EMBEDDED);
2037 1850 : EmbeddedNodes.insert(std::make_pair(new_node_id, local_embedded));
2038 : }
2039 :
2040 : // add to elem face edge
2041 10062 : cut_edge->addIntersection(position, local_embedded, edge_node1);
2042 10062 : if (cut_edge->numEmbeddedNodes() > 2)
2043 0 : EFAError("element edge can't have >2 embedded nodes");
2044 :
2045 : // cut fragment faces, which is an essential addition to other code in this method that cuts
2046 : // element faces
2047 : unsigned int FragFaceID;
2048 : unsigned int FragFaceEdgeID;
2049 10062 : if (getFragmentFaceEdgeID(face_id, edge_id, FragFaceID, FragFaceEdgeID))
2050 : {
2051 160 : EFAEdge * elem_edge = _faces[face_id]->getEdge(edge_id);
2052 160 : EFAEdge * frag_edge = getFragmentFace(0, FragFaceID)->getEdge(FragFaceEdgeID);
2053 : double xi[2] = {-1.0, -1.0}; // relative coords of two frag edge nodes
2054 160 : xi[0] = elem_edge->distanceFromNode1(frag_edge->getNode(0));
2055 160 : xi[1] = elem_edge->distanceFromNode1(frag_edge->getNode(1));
2056 160 : double frag_pos = (position - xi[0]) / (xi[1] - xi[0]);
2057 160 : EFANode * frag_edge_node1 = frag_edge->getNode(0);
2058 :
2059 160 : if (!frag_edge->hasIntersection())
2060 160 : frag_edge->addIntersection(frag_pos, local_embedded, frag_edge_node1);
2061 : }
2062 :
2063 : // add to adjacent face edge
2064 10062 : if (add_to_adjacent)
2065 : {
2066 5031 : double adj_pos = 1.0 - position;
2067 5031 : addFaceEdgeCut(
2068 : adj_face_id, adj_edge_id, adj_pos, local_embedded, EmbeddedNodes, false, false);
2069 : }
2070 : }
2071 :
2072 : // add cut to neighbor face edge
2073 109433 : if (add_to_neighbor)
2074 : {
2075 64198 : for (unsigned int en_iter = 0; en_iter < numFaceNeighbors(face_id); ++en_iter)
2076 : {
2077 27880 : EFAElement3D * face_neighbor = getFaceNeighbor(face_id, en_iter);
2078 27880 : unsigned int neigh_face_id = face_neighbor->getNeighborIndex(this);
2079 27880 : unsigned neigh_edge_id = getNeighborFaceEdgeID(face_id, edge_id, face_neighbor);
2080 27880 : double neigh_pos = 1.0 - position; // get emb node's postion on neighbor edge
2081 27880 : face_neighbor->addFaceEdgeCut(
2082 : neigh_face_id, neigh_edge_id, neigh_pos, local_embedded, EmbeddedNodes, false, true);
2083 : }
2084 :
2085 76522 : for (unsigned int en_iter = 0; en_iter < numEdgeNeighbors(face_id, edge_id); ++en_iter)
2086 : {
2087 40204 : EFAElement3D * edge_neighbor = getEdgeNeighbor(face_id, edge_id, en_iter);
2088 : unsigned int neigh_face_id, neigh_edge_id;
2089 40204 : getNeighborEdgeIndex(edge_neighbor, face_id, edge_id, neigh_face_id, neigh_edge_id);
2090 :
2091 : // Check the ordering of the node and the assign the position
2092 : double neigh_pos;
2093 80408 : if (_faces[face_id]->getEdge(edge_id)->getNode(0) ==
2094 40204 : edge_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id)->getNode(0))
2095 : neigh_pos = position;
2096 40122 : else if (_faces[face_id]->getEdge(edge_id)->getNode(1) ==
2097 20061 : edge_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id)->getNode(0))
2098 20061 : neigh_pos = 1.0 - position;
2099 : else
2100 0 : EFAError("The EFANodes on commaon edge are not matched.");
2101 :
2102 40204 : edge_neighbor->addFaceEdgeCut(
2103 : neigh_face_id, neigh_edge_id, neigh_pos, local_embedded, EmbeddedNodes, false, true);
2104 : }
2105 : } // If add_to_neighbor required
2106 109433 : }
2107 :
2108 : void
2109 0 : EFAElement3D::addFragFaceEdgeCut(unsigned int /*frag_face_id*/,
2110 : unsigned int /*frag_edge_id*/,
2111 : double /*position*/,
2112 : std::map<unsigned int, EFANode *> & /*EmbeddedNodes*/,
2113 : bool /*add_to_neighbor*/,
2114 : bool /*add_to_adjacent*/)
2115 : {
2116 : // TODO: mark frag face edges
2117 : // also need to check if cut has been added to this frag face edge or neighbor edge of adjacent
2118 : // face
2119 0 : }
2120 :
2121 : void
2122 20124 : EFAElement3D::checkNeighborFaceCut(unsigned int face_id,
2123 : unsigned int edge_id,
2124 : double position,
2125 : EFANode * from_node,
2126 : EFANode * embedded_node,
2127 : EFANode *& local_embedded)
2128 : {
2129 : // N.B. this is important. We are checking if the corresponding edge of the neighbor face or of
2130 : // the adjacent
2131 : // face's neighbor face has a cut at the same position. If so, use the existing embedded node as
2132 : // local_embedded
2133 35168 : for (unsigned int en_iter = 0; en_iter < numFaceNeighbors(face_id); ++en_iter)
2134 : {
2135 15044 : EFAElement3D * face_neighbor = getFaceNeighbor(face_id, en_iter);
2136 15044 : unsigned int neigh_face_id = face_neighbor->getNeighborIndex(this);
2137 15044 : unsigned neigh_edge_id = getNeighborFaceEdgeID(face_id, edge_id, face_neighbor);
2138 15044 : EFAEdge * neigh_edge = face_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id);
2139 :
2140 15044 : if (neigh_edge->hasIntersectionAtPosition(position, from_node))
2141 : {
2142 7588 : unsigned int emb_id = neigh_edge->getEmbeddedNodeIndex(position, from_node);
2143 7588 : EFANode * old_emb = neigh_edge->getEmbeddedNode(emb_id);
2144 :
2145 7588 : if (embedded_node && embedded_node != old_emb)
2146 0 : EFAError(
2147 : "attempting to add edge intersection when one already exists with different node.");
2148 7588 : if (local_embedded && local_embedded != old_emb)
2149 0 : EFAError("attempting to assign contradictory pointer to local_embedded.");
2150 :
2151 7588 : local_embedded = old_emb;
2152 : }
2153 : } // en_iter
2154 20124 : }
2155 :
2156 : void
2157 0 : EFAElement3D::mapParametricCoordinateFrom2DTo3D(unsigned int face_id,
2158 : std::vector<double> & xi_2d,
2159 : std::vector<double> & xi_3d) const
2160 : {
2161 : // given the coords of a point in a 2D face, translate it to 3D parametric coords
2162 0 : xi_3d.resize(3, 0.0);
2163 0 : if (_num_faces == 6)
2164 : {
2165 : if (face_id == 0)
2166 : {
2167 0 : xi_3d[0] = xi_2d[1];
2168 0 : xi_3d[1] = xi_2d[0];
2169 0 : xi_3d[2] = -1.0;
2170 : }
2171 : else if (face_id == 1)
2172 : {
2173 0 : xi_3d[0] = xi_2d[0];
2174 0 : xi_3d[1] = -1.0;
2175 0 : xi_3d[2] = xi_2d[1];
2176 : }
2177 : else if (face_id == 2)
2178 : {
2179 0 : xi_3d[0] = 1.0;
2180 0 : xi_3d[1] = xi_2d[0];
2181 0 : xi_3d[2] = xi_2d[1];
2182 : }
2183 : else if (face_id == 3)
2184 : {
2185 0 : xi_3d[0] = -xi_2d[0];
2186 0 : xi_3d[1] = 1.0;
2187 0 : xi_3d[2] = xi_2d[1];
2188 : }
2189 : else if (face_id == 4)
2190 : {
2191 0 : xi_3d[0] = -1.0;
2192 0 : xi_3d[1] = -xi_2d[0];
2193 0 : xi_3d[2] = xi_2d[1];
2194 : }
2195 : else if (face_id == 5)
2196 : {
2197 0 : xi_3d[0] = xi_2d[0];
2198 0 : xi_3d[1] = xi_2d[1];
2199 0 : xi_3d[2] = 1.0;
2200 : }
2201 : else
2202 0 : EFAError("face_id out of bounds");
2203 : }
2204 0 : else if (_num_faces == 4)
2205 : {
2206 0 : if (face_id == 0)
2207 : {
2208 0 : xi_3d[0] = xi_2d[0];
2209 0 : xi_3d[1] = xi_2d[1];
2210 0 : xi_3d[2] = 0.0;
2211 : }
2212 0 : else if (face_id == 1)
2213 : {
2214 0 : xi_3d[0] = 0.0;
2215 0 : xi_3d[1] = xi_2d[0];
2216 0 : xi_3d[2] = xi_2d[1];
2217 : }
2218 0 : else if (face_id == 2)
2219 : {
2220 0 : xi_3d[0] = xi_2d[1];
2221 0 : xi_3d[1] = 0.0;
2222 0 : xi_3d[2] = xi_2d[0];
2223 : }
2224 0 : else if (face_id == 3)
2225 : {
2226 0 : xi_3d[0] = xi_2d[0];
2227 0 : xi_3d[1] = xi_2d[2];
2228 0 : xi_3d[2] = xi_2d[1];
2229 : }
2230 : else
2231 0 : EFAError("face_id out of bounds");
2232 : }
2233 : else
2234 0 : EFAError("unknown element for 3D");
2235 0 : }
2236 :
2237 : std::vector<EFANode *>
2238 0 : EFAElement3D::getCommonNodes(const EFAElement3D * other_elem) const
2239 : {
2240 : std::set<EFANode *> e1nodes(_nodes.begin(),
2241 0 : _nodes.begin() + _num_vertices); // only account for corner nodes
2242 : std::set<EFANode *> e2nodes(other_elem->_nodes.begin(),
2243 0 : other_elem->_nodes.begin() + _num_vertices);
2244 0 : std::vector<EFANode *> common_nodes = Efa::getCommonElems(e1nodes, e2nodes);
2245 0 : return common_nodes;
2246 : }
|
