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 43474 : EFAElement3D::EFAElement3D(unsigned int eid, unsigned int n_nodes, unsigned int n_faces)
27 : : EFAElement(eid, n_nodes),
28 43474 : _num_faces(n_faces),
29 43474 : _faces(_num_faces, nullptr),
30 43474 : _face_neighbors(_num_faces),
31 86948 : _face_edge_neighbors(_num_faces)
32 : {
33 43474 : 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 22106 : else if (_num_faces == 6)
46 : {
47 22106 : _num_vertices = 8;
48 22106 : if (_num_nodes == 27)
49 388 : _num_interior_face_nodes = 5;
50 21718 : else if (_num_nodes == 20)
51 388 : _num_interior_face_nodes = 4;
52 21330 : else if (_num_nodes == 8)
53 21330 : _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 43474 : setLocalCoordinates();
61 43474 : }
62 :
63 2560 : EFAElement3D::EFAElement3D(const EFAElement3D * from_elem, bool convert_to_local)
64 2560 : : EFAElement(from_elem->_id, from_elem->_num_nodes),
65 2560 : _num_faces(from_elem->_num_faces),
66 2560 : _faces(_num_faces, nullptr),
67 2560 : _face_neighbors(_num_faces),
68 5120 : _face_edge_neighbors(_num_faces)
69 : {
70 2560 : if (convert_to_local)
71 : {
72 : // build local nodes from global nodes
73 29400 : for (unsigned int i = 0; i < _num_nodes; ++i)
74 : {
75 26840 : if (from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT ||
76 0 : from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_TEMP)
77 : {
78 26840 : _nodes[i] = from_elem->createLocalNodeFromGlobalNode(from_elem->_nodes[i]);
79 26840 : _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 15360 : for (unsigned int i = 0; i < _num_faces; ++i)
92 12800 : _faces[i] = new EFAFace(*from_elem->_faces[i]);
93 5120 : for (unsigned int i = 0; i < from_elem->_fragments.size(); ++i)
94 2560 : _fragments.push_back(new EFAFragment3D(this, true, from_elem, i));
95 2560 : 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 29400 : for (unsigned int i = 0; i < _num_nodes; ++i)
100 : {
101 26840 : if (_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
102 26840 : 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 2560 : findFacesAdjacentToFaces();
112 :
113 2560 : _local_node_coor = from_elem->_local_node_coor;
114 2560 : _num_interior_face_nodes = from_elem->_num_interior_face_nodes;
115 2560 : _num_vertices = from_elem->_num_vertices;
116 : }
117 : else
118 0 : EFAError("this EFAelement3D constructor only converts global nodes to local nodes");
119 2560 : }
120 :
121 89508 : EFAElement3D::~EFAElement3D()
122 : {
123 59900 : for (unsigned int i = 0; i < _fragments.size(); ++i)
124 13866 : if (_fragments[i])
125 13866 : delete _fragments[i];
126 :
127 276942 : for (unsigned int i = 0; i < _faces.size(); ++i)
128 230908 : if (_faces[i])
129 230908 : delete _faces[i];
130 :
131 46034 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
132 0 : if (_interior_nodes[i])
133 0 : delete _interior_nodes[i];
134 :
135 72874 : for (unsigned int i = 0; i < _local_nodes.size(); ++i)
136 26840 : if (_local_nodes[i])
137 26840 : delete _local_nodes[i];
138 89508 : }
139 :
140 : void
141 43474 : EFAElement3D::setLocalCoordinates()
142 : {
143 43474 : 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 22106 : _local_node_coor.resize(_num_nodes);
176 22106 : _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
177 22106 : _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
178 22106 : _local_node_coor[2] = EFAPoint(1.0, 1.0, 0.0);
179 22106 : _local_node_coor[3] = EFAPoint(0.0, 1.0, 0.0);
180 22106 : _local_node_coor[4] = EFAPoint(0.0, 0.0, 1.0);
181 22106 : _local_node_coor[5] = EFAPoint(1.0, 0.0, 1.0);
182 22106 : _local_node_coor[6] = EFAPoint(1.0, 1.0, 1.0);
183 22106 : _local_node_coor[7] = EFAPoint(0.0, 1.0, 1.0);
184 :
185 22106 : 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 22106 : 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 43474 : }
270 :
271 : unsigned int
272 625456 : EFAElement3D::numFragments() const
273 : {
274 625456 : return _fragments.size();
275 : }
276 :
277 : bool
278 867958 : EFAElement3D::isPartial() const
279 : {
280 867958 : if (_fragments.size() > 0)
281 : {
282 2068604 : for (unsigned int i = 0; i < _num_vertices; ++i)
283 : {
284 : bool node_in_frag = false;
285 2764230 : for (unsigned int j = 0; j < _fragments.size(); ++j)
286 : {
287 1982438 : if (_fragments[j]->containsNode(_nodes[i]))
288 : {
289 : node_in_frag = true;
290 : break;
291 : }
292 : } // j
293 :
294 1982438 : if (!node_in_frag)
295 : return true;
296 : } // i
297 : }
298 : return false;
299 : }
300 :
301 : void
302 3256 : 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 37254 : for (unsigned int i = 0; i < _nodes.size(); ++i)
307 33998 : non_physical_nodes.insert(_nodes[i]);
308 :
309 : // Now delete any nodes that are contained in fragments
310 : std::set<EFANode *>::iterator sit;
311 37254 : for (sit = non_physical_nodes.begin(); sit != non_physical_nodes.end();)
312 : {
313 : bool erased = false;
314 57444 : for (unsigned int i = 0; i < _fragments.size(); ++i)
315 : {
316 33998 : if (_fragments[i]->containsNode(*sit))
317 : {
318 10552 : non_physical_nodes.erase(sit++);
319 : erased = true;
320 10552 : break;
321 : }
322 : }
323 33998 : if (!erased)
324 : ++sit;
325 : }
326 3256 : }
327 :
328 : void
329 873626 : 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 11281916 : for (unsigned int i = 0; i < _num_nodes; ++i)
333 : {
334 10408290 : if (_nodes[i] == old_node)
335 13357 : _nodes[i] = new_node;
336 : }
337 1417041 : for (unsigned int i = 0; i < _fragments.size(); ++i)
338 543415 : _fragments[i]->switchNode(new_node, old_node);
339 :
340 4629910 : for (unsigned int i = 0; i < _faces.size(); ++i)
341 3756284 : _faces[i]->switchNode(new_node, old_node);
342 :
343 873626 : if (_parent && descend_to_parent)
344 : {
345 11441 : _parent->switchNode(new_node, old_node, false);
346 570988 : for (unsigned int i = 0; i < _parent->numGeneralNeighbors(); ++i)
347 : {
348 559547 : EFAElement * neigh_elem = _parent->getGeneralNeighbor(i); // generalized neighbor element
349 1319713 : for (unsigned int k = 0; k < neigh_elem->numChildren(); ++k)
350 760166 : neigh_elem->getChild(k)->switchNode(new_node, old_node, false);
351 : }
352 : }
353 873626 : }
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 2711 : EFAElement3D::updateFragmentNode()
368 : {
369 : // In EFAElement3D, updateFragmentNode needs to be implemented
370 2711 : }
371 :
372 : void
373 2611614 : 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 2611614 : master_nodes.clear();
379 2611614 : master_weights.clear();
380 : bool masters_found = false;
381 5571533 : for (unsigned int i = 0; i < _num_faces; ++i) // check element exterior faces
382 : {
383 5571533 : if (_faces[i]->containsNode(node))
384 : {
385 2611614 : masters_found = _faces[i]->getMasterInfo(node, master_nodes, master_weights);
386 2611614 : 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 2611614 : if (!masters_found)
421 0 : EFAError("In EFAelement3D::getMasterInfo, cannot find the given EFANode");
422 2611614 : }
423 :
424 : unsigned int
425 1657 : EFAElement3D::numInteriorNodes() const
426 : {
427 1657 : return _interior_nodes.size();
428 : }
429 :
430 : bool
431 495734 : EFAElement3D::overlaysElement(const EFAElement3D * other_elem) const
432 : {
433 : bool overlays = false;
434 : const EFAElement3D * other3d = dynamic_cast<const EFAElement3D *>(other_elem);
435 495734 : if (!other3d)
436 0 : EFAError("failed to dynamic cast to other3d");
437 :
438 : // Find indices of common nodes
439 495734 : const auto & common_face_curr = getCommonFaceID(other3d);
440 495734 : if (common_face_curr.size() == 1)
441 : {
442 158940 : unsigned int curr_face_id = common_face_curr[0];
443 158940 : EFAFace * curr_face = _faces[curr_face_id];
444 158940 : unsigned int other_face_id = other3d->getFaceID(curr_face);
445 158940 : EFAFace * other_face = other3d->_faces[other_face_id];
446 158940 : if (curr_face->hasSameOrientation(other_face))
447 : overlays = true;
448 : }
449 336794 : else if (common_face_curr.size() > 1)
450 : {
451 : // TODO: We probably need more error checking here.
452 : overlays = true;
453 : }
454 495734 : return overlays;
455 495734 : }
456 :
457 : unsigned int
458 176107 : EFAElement3D::getNeighborIndex(const EFAElement * neighbor_elem) const
459 : {
460 544178 : for (unsigned int i = 0; i < _num_faces; ++i)
461 812917 : for (unsigned int j = 0; j < _face_neighbors[i].size(); ++j)
462 444846 : if (_face_neighbors[i][j] == neighbor_elem)
463 176107 : return i;
464 0 : EFAError("in getNeighborIndex() element ", _id, " does not have neighbor ", neighbor_elem->id());
465 : }
466 :
467 : void
468 37900 : 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 37900 : EFAEdge * edge = this->getFace(face_id)->getEdge(edge_id);
475 78514 : for (unsigned int i = 0; i < neighbor_elem->numFaces(); ++i)
476 : {
477 258642 : for (unsigned int j = 0; j < neighbor_elem->getFace(i)->numEdges(); ++j)
478 : {
479 218028 : EFAEdge * neigh_edge = neighbor_elem->getFace(i)->getEdge(j);
480 218028 : if (neigh_edge->equivalent(*edge))
481 : {
482 37900 : neigh_face_id = i;
483 37900 : neigh_edge_id = j;
484 37900 : 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 41021 : EFAElement3D::clearNeighbors()
496 : {
497 41021 : _general_neighbors.clear();
498 246971 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
499 : {
500 205950 : _face_neighbors[face_iter].clear();
501 949398 : for (unsigned int edge_iter = 0; edge_iter < _faces[face_iter]->numEdges(); ++edge_iter)
502 743448 : _face_edge_neighbors[face_iter][edge_iter].clear();
503 : }
504 41021 : }
505 :
506 : void
507 41021 : EFAElement3D::setupNeighbors(std::map<EFANode *, std::set<EFAElement *>> & inverse_connectivity_map)
508 : {
509 41021 : findGeneralNeighbors(inverse_connectivity_map);
510 : std::vector<std::pair<unsigned int, unsigned int>> common_ids;
511 1546163 : for (unsigned int eit2 = 0; eit2 < _general_neighbors.size(); ++eit2)
512 : {
513 1505142 : EFAElement3D * neigh_elem = dynamic_cast<EFAElement3D *>(_general_neighbors[eit2]);
514 1505142 : if (!neigh_elem)
515 0 : EFAError("neighbor_elem is not of EFAelement3D type");
516 :
517 1505142 : const auto & common_face_id = getCommonFaceID(neigh_elem);
518 1841936 : if (common_face_id.empty() && getCommonEdgeID(neigh_elem, common_ids) &&
519 336794 : !overlaysElement(neigh_elem))
520 : {
521 : bool is_edge_neighbor = false;
522 :
523 : // Fragments must match up.
524 336794 : if ((_fragments.size() > 1) || (neigh_elem->numFragments() > 1))
525 0 : EFAError("in updateFaceNeighbors: Cannot have more than 1 fragment");
526 :
527 336794 : if ((_fragments.size() == 1) && (neigh_elem->numFragments() == 1))
528 : {
529 30676 : 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 1001840 : for (const auto & [face_id, edge_id] : common_ids)
537 667980 : _face_edge_neighbors[face_id][edge_id].push_back(neigh_elem);
538 : }
539 :
540 1505142 : if (common_face_id.size() == 1 && !overlaysElement(neigh_elem))
541 : {
542 157196 : unsigned int face_id = common_face_id[0];
543 : bool is_face_neighbor = false;
544 :
545 : // Fragments must match up.
546 157196 : if ((_fragments.size() > 1) || (neigh_elem->numFragments() > 1))
547 0 : EFAError("in updateFaceNeighbors: Cannot have more than 1 fragment");
548 :
549 157196 : if ((_fragments.size() == 1) && (neigh_elem->numFragments() == 1))
550 : {
551 16960 : 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 157052 : 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 157052 : _face_neighbors[face_id].push_back(neigh_elem);
567 : }
568 : }
569 1505142 : }
570 41021 : }
571 :
572 : void
573 41021 : EFAElement3D::neighborSanityCheck() const
574 : {
575 246971 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
576 : {
577 363002 : for (unsigned int en_iter = 0; en_iter < _face_neighbors[face_iter].size(); ++en_iter)
578 : {
579 157052 : EFAElement3D * neigh_elem = _face_neighbors[face_iter][en_iter];
580 157052 : if (neigh_elem != nullptr)
581 : {
582 : bool found_neighbor = false;
583 953932 : for (unsigned int face_iter2 = 0; face_iter2 < neigh_elem->numFaces(); ++face_iter2)
584 : {
585 1266552 : for (unsigned int en_iter2 = 0; en_iter2 < neigh_elem->numFaceNeighbors(face_iter2);
586 : ++en_iter2)
587 : {
588 626724 : if (neigh_elem->getFaceNeighbor(face_iter2, en_iter2) == this)
589 : {
590 157052 : 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 157052 : if (!found_neighbor)
599 0 : EFAError("Neighbor element doesn't recognize current element as neighbor");
600 : }
601 : }
602 : }
603 41021 : }
604 :
605 : void
606 40844 : EFAElement3D::initCrackTip(std::set<EFAElement *> & CrackTipElements)
607 : {
608 40844 : if (isCrackTipElement())
609 : {
610 477 : CrackTipElements.insert(this);
611 3051 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
612 : {
613 2574 : 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 976 : if (_face_neighbors[face_iter][0]->overlaysElement(this) ||
619 488 : _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 488 : 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 488 : _face_neighbors[face_iter][0]->setCrackTipSplit();
634 488 : _face_neighbors[face_iter][1]->setCrackTipSplit();
635 :
636 488 : _face_neighbors[face_iter][0]->addCrackTipNeighbor(this);
637 488 : _face_neighbors[face_iter][1]->addCrackTipNeighbor(this);
638 : }
639 : } // face_iter
640 : }
641 40844 : }
642 :
643 : bool
644 28261 : 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 28261 : if (_fragments.size() == 1)
654 : {
655 : std::set<EFAElement *>::iterator sit;
656 6343 : sit = CrackTipElements.find(this);
657 6343 : if (sit == CrackTipElements.end() && isCrackTipElement())
658 : should_duplicate = true;
659 3832 : else if (shouldDuplicateCrackTipSplitElement(CrackTipElements))
660 : should_duplicate = true;
661 3256 : else if (shouldDuplicateForPhantomCorner())
662 : should_duplicate = true;
663 : }
664 28261 : return should_duplicate;
665 : }
666 :
667 : bool
668 3832 : 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 3832 : if (_fragments.size() == 1)
676 : {
677 : std::vector<unsigned int> split_neighbors;
678 3832 : 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 3256 : getNonPhysicalNodes(non_physical_nodes);
687 :
688 125888 : for (unsigned int eit = 0; eit < _general_neighbors.size(); ++eit)
689 : {
690 122632 : EFAElement3D * neigh_elem = dynamic_cast<EFAElement3D *>(_general_neighbors[eit]);
691 122632 : 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 122632 : sit = CrackTipElements.find(neigh_elem);
697 122632 : 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 3832 : } // IF only one fragment
715 3832 : return should_duplicate;
716 : }
717 :
718 : bool
719 3256 : 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 3256 : if (_fragments.size() == 1 && (!_crack_tip_split_element))
726 : {
727 17158 : for (unsigned int i = 0; i < _num_faces; ++i)
728 : {
729 14284 : std::set<EFANode *> phantom_nodes = getPhantomNodeOnFace(i);
730 14284 : if (phantom_nodes.size() > 0 && numFaceNeighbors(i) == 1)
731 : {
732 7240 : EFAElement3D * neighbor_elem = _face_neighbors[i][0];
733 7240 : 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 3256 : return should_duplicate;
755 : }
756 :
757 : bool
758 3832 : 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 3832 : if (_fragments.size() == 1 && _crack_tip_split_element)
765 : {
766 1988 : for (unsigned int i = 0; i < _crack_tip_neighbors.size(); ++i)
767 : {
768 1030 : unsigned int neigh_idx = _crack_tip_neighbors[i]; // essentially a face_id
769 1030 : 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 1030 : EFAElement3D * neighbor_elem = _face_neighbors[neigh_idx][0];
778 1030 : 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 1030 : else if (neighbor_elem->numFragments() == 2)
787 : {
788 576 : EFAFragment3D * neigh_frag1 = neighbor_elem->getFragment(0);
789 576 : EFAFragment3D * neigh_frag2 = neighbor_elem->getFragment(1);
790 576 : 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 2880 : for (unsigned int j = 0; j < neigh_cut_nodes.size(); ++j)
793 : {
794 2304 : if (_faces[neigh_idx]->containsNode(neigh_cut_nodes[j]))
795 1152 : counter += 1;
796 : }
797 576 : if (counter == 2)
798 : {
799 576 : split_neighbors.push_back(neigh_idx);
800 : will_extend = true;
801 : }
802 576 : }
803 : } // i
804 : }
805 3832 : return will_extend;
806 : }
807 :
808 : bool
809 49676 : EFAElement3D::isCrackTipElement() const
810 : {
811 49676 : return fragmentHasTipFaces();
812 : }
813 :
814 : unsigned int
815 5980 : EFAElement3D::getNumCuts() const
816 : {
817 : unsigned int num_cut_faces = 0;
818 36100 : for (unsigned int i = 0; i < _num_faces; ++i)
819 30120 : if (_faces[i]->hasIntersection())
820 0 : num_cut_faces += 1;
821 5980 : return num_cut_faces;
822 : }
823 :
824 : bool
825 20706 : 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 20706 : if (_fragments.size() > 0)
831 : {
832 : unsigned int num_interior_faces = 0;
833 16822 : for (unsigned int i = 0; i < _fragments[0]->numFaces(); ++i)
834 : {
835 14190 : if (_fragments[0]->isFaceInterior(i))
836 2434 : num_interior_faces += 1;
837 : }
838 2632 : if (num_interior_faces == 3)
839 : cut_third = true;
840 : }
841 20706 : return cut_third;
842 : }
843 :
844 : void
845 26673 : 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 26673 : sit = CrackTipElements.find(this);
851 26673 : if (sit != CrackTipElements.end()) // curr_elem is a crack tip element
852 : {
853 254 : if (_fragments.size() == 1)
854 254 : _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 26673 : if (_fragments.size() == 1)
862 3352 : _fragments[0]->removeInvalidEmbeddedNodes(EmbeddedNodes);
863 :
864 : // for an element with no fragment, create one fragment identical to the element
865 26673 : if (_fragments.size() == 0)
866 23321 : _fragments.push_back(new EFAFragment3D(this, true, this));
867 26673 : if (_fragments.size() != 1)
868 0 : EFAError("Element ", _id, " must have 1 fragment at this point");
869 :
870 : // count fragment's cut faces
871 26673 : unsigned int num_cut_frag_faces = _fragments[0]->getNumCuts();
872 26673 : unsigned int num_frag_faces = _fragments[0]->numFaces();
873 26673 : 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 26673 : if (num_cut_frag_faces == 0)
878 : {
879 25016 : if (!isPartial()) // delete the temp frag for an uncut elem
880 : {
881 21918 : delete _fragments[0];
882 21918 : _fragments.clear();
883 : }
884 25016 : return;
885 : }
886 :
887 : // split one fragment into one or two new fragments
888 1657 : std::vector<EFAFragment3D *> new_frags = _fragments[0]->split();
889 1657 : if (new_frags.size() == 1 || new_frags.size() == 2)
890 : {
891 1657 : delete _fragments[0]; // delete the old fragment
892 1657 : _fragments.clear();
893 4526 : for (unsigned int i = 0; i < new_frags.size(); ++i)
894 2869 : _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 1657 : fragmentSanityCheck(num_frag_faces, num_cut_frag_faces);
900 1657 : }
901 :
902 : void
903 1657 : EFAElement3D::fragmentSanityCheck(unsigned int n_old_frag_faces, unsigned int n_old_frag_cuts) const
904 : {
905 1657 : unsigned int n_interior_nodes = numInteriorNodes();
906 1657 : 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 1657 : 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 1657 : else if (fragmentHasTipFaces()) // crack tip case
915 : {
916 445 : 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 1212 : if (_fragments.size() != 2 || (_fragments[0]->numFaces() + _fragments[1]->numFaces()) !=
922 1212 : n_old_frag_faces + n_old_frag_cuts + 2)
923 0 : EFAError("Incorrect link size for element that has been completely cut");
924 : }
925 1657 : }
926 :
927 : void
928 5980 : EFAElement3D::restoreFragment(const EFAElement * const from_elem)
929 : {
930 5980 : const EFAElement3D * from_elem3d = dynamic_cast<const EFAElement3D *>(from_elem);
931 5980 : if (!from_elem3d)
932 0 : EFAError("from_elem is not of EFAelement3D type");
933 :
934 : // restore fragments
935 5980 : if (_fragments.size() != 0)
936 0 : EFAError("in restoreFragmentInfo elements must not have any pre-existing fragments");
937 11960 : for (unsigned int i = 0; i < from_elem3d->numFragments(); ++i)
938 5980 : _fragments.push_back(new EFAFragment3D(this, true, from_elem3d, i));
939 :
940 : // restore interior nodes
941 5980 : if (_interior_nodes.size() != 0)
942 0 : EFAError("in restoreFragmentInfo elements must not have any pre-exsiting interior nodes");
943 5980 : 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 5980 : if (getNumCuts() != 0)
948 0 : EFAError("In restoreEdgeIntersection: edge cuts already exist in element ", _id);
949 36100 : for (unsigned int i = 0; i < _num_faces; ++i)
950 30120 : _faces[i]->copyIntersection(*from_elem3d->_faces[i]);
951 :
952 : // replace all local nodes with global nodes
953 68130 : for (unsigned int i = 0; i < from_elem3d->numNodes(); ++i)
954 : {
955 62150 : if (from_elem3d->_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
956 62150 : 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 5980 : }
962 :
963 : void
964 26673 : 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 26673 : if (_children.size() != 0)
972 0 : EFAError("Element cannot have existing children in createChildElements");
973 :
974 26673 : if (_fragments.size() > 1 || shouldDuplicateForCrackTip(CrackTipElements))
975 : {
976 1499 : if (_fragments.size() > 2)
977 0 : EFAError("More than 2 fragments not yet supported");
978 :
979 : // set up the children
980 1499 : ParentElements.push_back(this);
981 4210 : for (unsigned int ichild = 0; ichild < _fragments.size(); ++ichild)
982 : {
983 : unsigned int new_elem_id;
984 2711 : if (newChildElements.size() == 0)
985 138 : new_elem_id = Efa::getNewID(Elements);
986 : else
987 2573 : new_elem_id = Efa::getNewID(newChildElements);
988 :
989 2711 : EFAElement3D * childElem = new EFAElement3D(new_elem_id, this->numNodes(), this->numFaces());
990 2711 : newChildElements.insert(std::make_pair(new_elem_id, childElem));
991 :
992 2711 : ChildElements.push_back(childElem);
993 2711 : childElem->setParent(this);
994 2711 : _children.push_back(childElem);
995 :
996 : std::vector<std::vector<EFANode *>> cut_plane_nodes;
997 17396 : for (unsigned int i = 0; i < this->getFragment(ichild)->numFaces(); ++i)
998 : {
999 14685 : if (this->getFragment(ichild)->isFaceInterior(i))
1000 : {
1001 2488 : EFAFace * face = this->getFragment(ichild)->getFace(i);
1002 : std::vector<EFANode *> node_line;
1003 11496 : for (unsigned int j = 0; j < face->numNodes(); ++j)
1004 9008 : node_line.push_back(face->getNode(j));
1005 2488 : cut_plane_nodes.push_back(node_line);
1006 2488 : }
1007 : }
1008 :
1009 : std::vector<EFAPoint> cut_plane_points;
1010 :
1011 2711 : EFAPoint normal(0.0, 0.0, 0.0);
1012 2711 : EFAPoint orig(0.0, 0.0, 0.0);
1013 :
1014 2711 : if (cut_plane_nodes.size())
1015 : {
1016 11496 : for (unsigned int i = 0; i < cut_plane_nodes[0].size(); ++i)
1017 : {
1018 : std::vector<EFANode *> master_nodes;
1019 : std::vector<double> master_weights;
1020 :
1021 9008 : this->getMasterInfo(cut_plane_nodes[0][i], master_nodes, master_weights);
1022 9008 : EFAPoint coor(0.0, 0.0, 0.0);
1023 27024 : for (unsigned int i = 0; i < master_nodes.size(); ++i)
1024 : {
1025 18016 : EFANode * local = this->createLocalNodeFromGlobalNode(master_nodes[i]);
1026 18016 : coor += _local_node_coor[local->id()] * master_weights[i];
1027 18016 : delete local;
1028 : }
1029 9008 : cut_plane_points.push_back(coor);
1030 9008 : }
1031 11496 : for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
1032 9008 : orig += cut_plane_points[i];
1033 2488 : orig /= cut_plane_points.size();
1034 :
1035 2488 : EFAPoint center(0.0, 0.0, 0.0);
1036 11496 : for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
1037 9008 : center += cut_plane_points[i];
1038 2488 : center /= cut_plane_points.size();
1039 :
1040 11496 : for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
1041 : {
1042 9008 : unsigned int iplus1 = i < cut_plane_points.size() - 1 ? i + 1 : 0;
1043 9008 : EFAPoint ray1 = cut_plane_points[i] - center;
1044 9008 : EFAPoint ray2 = cut_plane_points[iplus1] - center;
1045 9008 : normal += ray1.cross(ray2);
1046 : }
1047 2488 : normal /= cut_plane_points.size();
1048 2488 : Xfem::normalizePoint(normal);
1049 : }
1050 :
1051 : // get child element's nodes
1052 30759 : for (unsigned int j = 0; j < _num_nodes; ++j)
1053 : {
1054 28048 : EFAPoint p(0.0, 0.0, 0.0);
1055 28048 : p = _local_node_coor[j];
1056 28048 : EFAPoint origin_to_point = p - orig;
1057 28048 : if (_fragments.size() == 1 && !shouldDuplicateForCrackTip(CrackTipElements))
1058 0 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1059 28048 : else if (origin_to_point * normal < Xfem::tol)
1060 15168 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1061 : else // parent element's node is not in fragment
1062 : {
1063 12880 : unsigned int new_node_id = Efa::getNewID(TempNodes);
1064 12880 : EFANode * newNode = new EFANode(new_node_id, EFANode::N_CATEGORY_TEMP, _nodes[j]);
1065 12880 : TempNodes.insert(std::make_pair(new_node_id, newNode));
1066 12880 : childElem->setNode(j, newNode); // be a temp node
1067 : }
1068 : }
1069 :
1070 : // get child element's fragments
1071 2711 : EFAFragment3D * new_frag = new EFAFragment3D(childElem, true, this, ichild);
1072 2711 : childElem->_fragments.push_back(new_frag);
1073 :
1074 : // get child element's faces and set up adjacent faces
1075 2711 : childElem->createFaces();
1076 16417 : for (unsigned int j = 0; j < _num_faces; ++j)
1077 13706 : childElem->_faces[j]->copyIntersection(*_faces[j]);
1078 2711 : childElem->removePhantomEmbeddedNode(); // IMPORTANT
1079 :
1080 : // inherit old interior nodes
1081 2711 : for (unsigned int j = 0; j < _interior_nodes.size(); ++j)
1082 0 : childElem->_interior_nodes.push_back(new EFAVolumeNode(*_interior_nodes[j]));
1083 2711 : }
1084 : }
1085 : else // num_links == 1 || num_links == 0
1086 : {
1087 : // child is itself - but don't insert into the list of ChildElements!!!
1088 25174 : _children.push_back(this);
1089 : }
1090 26673 : }
1091 :
1092 : void
1093 2711 : EFAElement3D::removePhantomEmbeddedNode()
1094 : {
1095 : // remove the embedded nodes on faces that are outside the real domain
1096 2711 : if (_fragments.size() > 0)
1097 : {
1098 16417 : for (unsigned int i = 0; i < _num_faces; ++i)
1099 : {
1100 : // get emb nodes to be removed on edges
1101 : std::vector<EFANode *> nodes_to_delete;
1102 63410 : for (unsigned int j = 0; j < _faces[i]->numEdges(); ++j)
1103 : {
1104 : EFAEdge * edge = _faces[i]->getEdge(j);
1105 68636 : for (unsigned int k = 0; k < edge->numEmbeddedNodes(); ++k)
1106 : {
1107 18932 : if (!_fragments[0]->containsNode(edge->getEmbeddedNode(k)))
1108 0 : nodes_to_delete.push_back(edge->getEmbeddedNode(k));
1109 : } // k
1110 : } // j
1111 :
1112 : // get emb nodes to be removed in the face interior
1113 13706 : for (unsigned int j = 0; j < _faces[i]->numInteriorNodes(); ++j)
1114 : {
1115 0 : EFANode * face_node = _faces[i]->getInteriorNode(j)->getNode();
1116 0 : if (!_fragments[0]->containsNode(face_node))
1117 0 : nodes_to_delete.push_back(face_node);
1118 : } // j
1119 :
1120 : // remove all invalid embedded nodes
1121 13706 : for (unsigned int j = 0; j < nodes_to_delete.size(); ++j)
1122 0 : _faces[i]->removeEmbeddedNode(nodes_to_delete[j]);
1123 13706 : } // i
1124 : }
1125 2711 : }
1126 :
1127 : void
1128 2711 : EFAElement3D::connectNeighbors(std::map<unsigned int, EFANode *> & PermanentNodes,
1129 : std::map<unsigned int, EFANode *> & TempNodes,
1130 : std::map<EFANode *, std::set<EFAElement *>> & InverseConnectivityMap,
1131 : bool merge_phantom_faces)
1132 : {
1133 : // N.B. "this" must point to a child element that was just created
1134 2711 : if (!_parent)
1135 0 : EFAError("no parent element for child element ", _id, " in connect_neighbors");
1136 2711 : EFAElement3D * parent3d = dynamic_cast<EFAElement3D *>(_parent);
1137 2711 : if (!parent3d)
1138 0 : EFAError("cannot dynamic cast to parent3d in connect_neighbors");
1139 :
1140 : // First loop through edges and merge nodes with neighbors as appropriate
1141 16417 : for (unsigned int j = 0; j < _num_faces; ++j)
1142 : {
1143 24673 : for (unsigned int k = 0; k < parent3d->numFaceNeighbors(j); ++k)
1144 : {
1145 10967 : EFAElement3D * NeighborElem = parent3d->getFaceNeighbor(j, k);
1146 10967 : unsigned int neighbor_face_id = NeighborElem->getNeighborIndex(parent3d);
1147 :
1148 10967 : if (_faces[j]->hasIntersection())
1149 : {
1150 20235 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1151 : {
1152 : EFAElement3D * childOfNeighborElem =
1153 13252 : dynamic_cast<EFAElement3D *>(NeighborElem->getChild(l));
1154 13252 : if (!childOfNeighborElem)
1155 0 : EFAError("dynamic cast childOfNeighborElem fails");
1156 :
1157 : // Check to see if the nodes are already merged. There's nothing else to do in that case.
1158 13252 : EFAFace * neighborChildFace = childOfNeighborElem->getFace(neighbor_face_id);
1159 13252 : if (_faces[j]->equivalent(neighborChildFace))
1160 3935 : continue;
1161 :
1162 9317 : if (_fragments[0]->isConnected(childOfNeighborElem->getFragment(0)))
1163 : {
1164 14681 : for (unsigned int i = 0; i < _faces[j]->numNodes(); ++i)
1165 : {
1166 : unsigned int childNodeIndex = i;
1167 : unsigned int neighborChildNodeIndex =
1168 11436 : parent3d->getNeighborFaceNodeID(j, childNodeIndex, NeighborElem);
1169 :
1170 11436 : EFANode * childNode = _faces[j]->getNode(childNodeIndex);
1171 11436 : EFANode * childOfNeighborNode = neighborChildFace->getNode(neighborChildNodeIndex);
1172 11436 : mergeNodes(
1173 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1174 : } // i
1175 :
1176 9081 : for (unsigned int m = 0; m < _num_interior_face_nodes; ++m)
1177 : {
1178 : unsigned int childNodeIndex = m;
1179 : unsigned int neighborChildNodeIndex =
1180 5836 : parent3d->getNeighborFaceInteriorNodeID(j, childNodeIndex, NeighborElem);
1181 :
1182 5836 : EFANode * childNode = _faces[j]->getInteriorFaceNode(childNodeIndex);
1183 : EFANode * childOfNeighborNode =
1184 5836 : neighborChildFace->getInteriorFaceNode(neighborChildNodeIndex);
1185 5836 : mergeNodes(
1186 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1187 : } // m
1188 : }
1189 : } // l, loop over NeighborElem's children
1190 : }
1191 : else // No edge intersection -- optionally merge non-material nodes if they share a common
1192 : // parent
1193 : {
1194 3984 : if (merge_phantom_faces)
1195 : {
1196 7968 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1197 : {
1198 : EFAElement3D * childOfNeighborElem =
1199 3984 : dynamic_cast<EFAElement3D *>(NeighborElem->getChild(l));
1200 3984 : if (!childOfNeighborElem)
1201 0 : EFAError("dynamic cast childOfNeighborElem fails");
1202 :
1203 3984 : EFAFace * neighborChildFace = childOfNeighborElem->getFace(neighbor_face_id);
1204 3984 : if (!neighborChildFace
1205 3984 : ->hasIntersection()) // neighbor face must NOT have intersection either
1206 : {
1207 : // Check to see if the nodes are already merged. There's nothing else to do in that
1208 : // case.
1209 3984 : if (_faces[j]->equivalent(neighborChildFace))
1210 2404 : continue;
1211 :
1212 7452 : for (unsigned int i = 0; i < _faces[j]->numNodes(); ++i)
1213 : {
1214 : unsigned int childNodeIndex = i;
1215 : unsigned int neighborChildNodeIndex =
1216 5872 : parent3d->getNeighborFaceNodeID(j, childNodeIndex, NeighborElem);
1217 :
1218 5872 : EFANode * childNode = _faces[j]->getNode(childNodeIndex);
1219 5872 : EFANode * childOfNeighborNode = neighborChildFace->getNode(neighborChildNodeIndex);
1220 :
1221 11193 : if (childNode->parent() != nullptr &&
1222 5321 : childNode->parent() ==
1223 : childOfNeighborNode
1224 5321 : ->parent()) // non-material node and both come from same parent
1225 0 : mergeNodes(childNode,
1226 : childOfNeighborNode,
1227 : childOfNeighborElem,
1228 : PermanentNodes,
1229 : TempNodes);
1230 : } // i
1231 : }
1232 : } // loop over NeighborElem's children
1233 : } // if (merge_phantom_edges)
1234 : } // IF edge-j has_intersection()
1235 : } // k, loop over neighbors on edge j
1236 : } // j, loop over all faces
1237 :
1238 : // Now do a second loop through faces and convert remaining nodes to permanent nodes.
1239 : // If there is no neighbor on that face, also duplicate the embedded node if it exists
1240 30759 : for (unsigned int j = 0; j < _num_nodes; ++j)
1241 : {
1242 28048 : EFANode * childNode = _nodes[j];
1243 28048 : if (childNode->category() == EFANode::N_CATEGORY_TEMP)
1244 : {
1245 : // if current child element does not have siblings, and if current temp node is a lone one
1246 : // this temp node should be merged back to its parent permanent node. Otherwise we would have
1247 : // permanent nodes that are not connected to any element
1248 1588 : std::set<EFAElement *> patch_elems = InverseConnectivityMap[childNode->parent()];
1249 1588 : if (parent3d->numFragments() == 1 && patch_elems.size() == 1)
1250 0 : switchNode(childNode->parent(), childNode, false);
1251 : else
1252 : {
1253 1588 : unsigned int new_node_id = Efa::getNewID(PermanentNodes);
1254 : EFANode * newNode =
1255 1588 : new EFANode(new_node_id, EFANode::N_CATEGORY_PERMANENT, childNode->parent());
1256 1588 : PermanentNodes.insert(std::make_pair(new_node_id, newNode));
1257 1588 : switchNode(newNode, childNode, false);
1258 : }
1259 1588 : if (!Efa::deleteFromMap(TempNodes, childNode))
1260 0 : EFAError(
1261 : "Attempted to delete node: ", childNode->id(), " from TempNodes, but couldn't find it");
1262 : }
1263 : }
1264 2711 : }
1265 :
1266 : void
1267 0 : EFAElement3D::printElement(std::ostream & ostream) const
1268 : {
1269 : // first line: all elem faces
1270 : ostream << std::setw(5);
1271 0 : ostream << _id << "| ";
1272 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1273 : {
1274 0 : for (unsigned int k = 0; k < _faces[j]->numNodes(); ++k)
1275 0 : ostream << std::setw(5) << _faces[j]->getNode(k)->idCatString();
1276 0 : ostream << " | ";
1277 : }
1278 : ostream << std::endl;
1279 :
1280 : // second line: emb nodes in all faces + neighbor of each face
1281 : ostream << std::setw(5);
1282 : ostream << "embed"
1283 0 : << "| ";
1284 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1285 : {
1286 0 : for (unsigned int k = 0; k < _faces[j]->numEdges(); ++k)
1287 : {
1288 : ostream << std::setw(4);
1289 0 : if (_faces[j]->getEdge(k)->hasIntersection())
1290 : {
1291 0 : if (_faces[j]->getEdge(k)->numEmbeddedNodes() > 1)
1292 : {
1293 0 : ostream << "[";
1294 0 : for (unsigned int l = 0; l < _faces[j]->getEdge(k)->numEmbeddedNodes(); ++l)
1295 : {
1296 0 : ostream << _faces[j]->getEdge(k)->getEmbeddedNode(l)->id() << " ";
1297 0 : if (l == _faces[j]->getEdge(k)->numEmbeddedNodes() - 1)
1298 0 : ostream << "]";
1299 : else
1300 0 : ostream << " ";
1301 : } // l
1302 : }
1303 : else
1304 0 : ostream << _faces[j]->getEdge(k)->getEmbeddedNode(0)->id() << " ";
1305 : }
1306 : else
1307 0 : ostream << " -- ";
1308 : } // k
1309 0 : ostream << " | ";
1310 : } // j
1311 : ostream << std::endl;
1312 :
1313 : // third line: neighbors
1314 : ostream << std::setw(5);
1315 : ostream << "neigh"
1316 0 : << "| ";
1317 0 : for (unsigned int j = 0; j < _num_faces; ++j)
1318 : {
1319 : ostream << std::setw(4);
1320 0 : if (numFaceNeighbors(j) > 1)
1321 : {
1322 0 : ostream << "[";
1323 0 : for (unsigned int k = 0; k < numFaceNeighbors(j); ++k)
1324 : {
1325 0 : ostream << getFaceNeighbor(j, k)->id();
1326 0 : if (k == numFaceNeighbors(j) - 1)
1327 0 : ostream << "]";
1328 : else
1329 0 : ostream << " ";
1330 : }
1331 : }
1332 : else
1333 : {
1334 0 : if (numFaceNeighbors(j) == 1)
1335 0 : ostream << getFaceNeighbor(j, 0)->id() << " ";
1336 : else
1337 0 : ostream << " -- ";
1338 : }
1339 : }
1340 : ostream << std::endl;
1341 :
1342 : // fourth line: fragments
1343 0 : for (unsigned int j = 0; j < _fragments.size(); ++j)
1344 : {
1345 : ostream << std::setw(4);
1346 0 : ostream << "frag" << j << "| ";
1347 0 : for (unsigned int k = 0; k < _fragments[j]->numFaces(); ++k)
1348 : {
1349 0 : for (unsigned int l = 0; l < _fragments[j]->getFace(k)->numNodes(); ++l)
1350 0 : ostream << std::setw(5) << _fragments[j]->getFace(k)->getNode(l)->idCatString();
1351 0 : ostream << " | ";
1352 : }
1353 : ostream << std::endl;
1354 : }
1355 : ostream << std::endl;
1356 0 : }
1357 :
1358 : EFAFragment3D *
1359 11540970 : EFAElement3D::getFragment(unsigned int frag_id) const
1360 : {
1361 11540970 : if (frag_id < _fragments.size())
1362 11540970 : return _fragments[frag_id];
1363 : else
1364 0 : EFAError("frag_id out of bounds");
1365 : }
1366 :
1367 : std::set<EFANode *>
1368 0 : EFAElement3D::getFaceNodes(unsigned int face_id) const
1369 : {
1370 : std::set<EFANode *> face_nodes;
1371 0 : for (unsigned int i = 0; i < _faces[face_id]->numNodes(); ++i)
1372 0 : face_nodes.insert(_faces[face_id]->getNode(i));
1373 0 : return face_nodes;
1374 : }
1375 :
1376 : bool
1377 0 : EFAElement3D::getFaceNodeParametricCoordinates(EFANode * node, std::vector<double> & xi_3d) const
1378 : {
1379 : // get the parametric coords of a node in an element face
1380 : unsigned int face_id = std::numeric_limits<unsigned int>::max();
1381 : bool face_found = false;
1382 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1383 : {
1384 0 : if (_faces[i]->containsNode(node))
1385 : {
1386 : face_id = i;
1387 : face_found = true;
1388 : break;
1389 : }
1390 : }
1391 0 : if (face_found)
1392 : {
1393 0 : std::vector<double> xi_2d(2, 0.0);
1394 0 : if (_faces[face_id]->getFaceNodeParametricCoords(node, xi_2d))
1395 0 : mapParametricCoordinateFrom2DTo3D(face_id, xi_2d, xi_3d);
1396 : else
1397 0 : EFAError("failed to get the 2D para coords on the face");
1398 0 : }
1399 0 : return face_found;
1400 : }
1401 :
1402 : EFAVolumeNode *
1403 0 : EFAElement3D::getInteriorNode(unsigned int interior_node_id) const
1404 : {
1405 0 : if (interior_node_id < _interior_nodes.size())
1406 0 : return _interior_nodes[interior_node_id];
1407 : else
1408 0 : EFAError("interior_node_id out of bounds");
1409 : }
1410 :
1411 : void
1412 0 : EFAElement3D::removeEmbeddedNode(EFANode * emb_node, bool remove_for_neighbor)
1413 : {
1414 0 : for (unsigned int i = 0; i < _fragments.size(); ++i)
1415 0 : _fragments[i]->removeEmbeddedNode(emb_node);
1416 :
1417 0 : for (unsigned int i = 0; i < _faces.size(); ++i)
1418 0 : _faces[i]->removeEmbeddedNode(emb_node);
1419 :
1420 0 : if (remove_for_neighbor)
1421 : {
1422 0 : for (unsigned int i = 0; i < numFaces(); ++i)
1423 0 : for (unsigned int j = 0; j < numFaceNeighbors(i); ++j)
1424 0 : getFaceNeighbor(i, j)->removeEmbeddedNode(emb_node, false);
1425 : }
1426 0 : }
1427 :
1428 : unsigned int
1429 18829009 : EFAElement3D::numFaces() const
1430 : {
1431 18829009 : return _faces.size();
1432 : }
1433 :
1434 : void
1435 0 : EFAElement3D::setFace(unsigned int face_id, EFAFace * face)
1436 : {
1437 0 : _faces[face_id] = face;
1438 0 : }
1439 :
1440 : void
1441 43474 : EFAElement3D::createFaces()
1442 : {
1443 : // create element faces based on existing element nodes
1444 43474 : int hex_local_node_indices[6][4] = {
1445 : {0, 3, 2, 1}, {0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}};
1446 43474 : int tet_local_node_indices[4][3] = {{0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}};
1447 :
1448 43474 : int hex_interior_face_node_indices[6][5] = {{8, 9, 10, 11, 20},
1449 : {8, 13, 16, 12, 21},
1450 : {9, 14, 17, 13, 22},
1451 : {10, 14, 18, 15, 23},
1452 : {11, 15, 19, 12, 24},
1453 : {16, 17, 18, 19, 25}};
1454 43474 : int tet_interior_face_node_indices[4][4] = {
1455 : {4, 5, 6, 10}, {4, 7, 8, 11}, {5, 8, 9, 12}, {6, 7, 9, 13}};
1456 :
1457 43474 : _faces = std::vector<EFAFace *>(_num_faces, nullptr);
1458 43474 : if (_num_nodes == 8 || _num_nodes == 20 || _num_nodes == 27)
1459 : {
1460 22106 : if (_num_faces != 6)
1461 0 : EFAError("num_faces of hexes must be 6");
1462 154742 : for (unsigned int i = 0; i < _num_faces; ++i)
1463 : {
1464 132636 : _faces[i] = new EFAFace(4, _num_interior_face_nodes);
1465 663180 : for (unsigned int j = 0; j < 4; ++j)
1466 530544 : _faces[i]->setNode(j, _nodes[hex_local_node_indices[i][j]]);
1467 132636 : _faces[i]->createEdges();
1468 153588 : for (unsigned int k = 0; k < _num_interior_face_nodes; ++k)
1469 20952 : _faces[i]->setInteriorFaceNode(k, _nodes[hex_interior_face_node_indices[i][k]]);
1470 : }
1471 : }
1472 : else if (_num_nodes == 4 || _num_nodes == 10 || _num_nodes == 14)
1473 : {
1474 21368 : if (_num_faces != 4)
1475 0 : EFAError("num_faces of tets must be 4");
1476 106840 : for (unsigned int i = 0; i < _num_faces; ++i)
1477 : {
1478 85472 : _faces[i] = new EFAFace(3, _num_interior_face_nodes);
1479 341888 : for (unsigned int j = 0; j < 3; ++j)
1480 256416 : _faces[i]->setNode(j, _nodes[tet_local_node_indices[i][j]]);
1481 85472 : _faces[i]->createEdges();
1482 384624 : for (unsigned int k = 0; k < _num_interior_face_nodes; ++k)
1483 299152 : _faces[i]->setInteriorFaceNode(k, _nodes[tet_interior_face_node_indices[i][k]]);
1484 : }
1485 : }
1486 : else
1487 0 : EFAError("unknown 3D element type in createFaces()");
1488 :
1489 261582 : for (unsigned int face_iter = 0; face_iter < _num_faces; ++face_iter)
1490 : {
1491 218108 : _face_edge_neighbors[face_iter].resize(_faces[face_iter]->numEdges());
1492 1005068 : for (unsigned int edge_iter = 0; edge_iter < _faces[face_iter]->numEdges(); ++edge_iter)
1493 786960 : _face_edge_neighbors[face_iter][edge_iter].clear();
1494 : }
1495 :
1496 : // create element face connectivity array
1497 43474 : findFacesAdjacentToFaces(); // IMPORTANT
1498 43474 : }
1499 :
1500 : EFAFace *
1501 17543635 : EFAElement3D::getFace(unsigned int face_id) const
1502 : {
1503 17543635 : return _faces[face_id];
1504 : }
1505 :
1506 : unsigned int
1507 258333 : EFAElement3D::getFaceID(EFAFace * face) const
1508 : {
1509 : bool found_face_id = false;
1510 : unsigned int face_id;
1511 828467 : for (unsigned int iface = 0; iface < _num_faces; ++iface)
1512 : {
1513 828467 : if (_faces[iface]->equivalent(face))
1514 : {
1515 : face_id = iface;
1516 : found_face_id = true;
1517 : break;
1518 : }
1519 : }
1520 258333 : if (!found_face_id)
1521 0 : EFAError("input face not found in get_face_id()");
1522 258333 : return face_id;
1523 : }
1524 :
1525 : std::vector<unsigned int>
1526 2000876 : EFAElement3D::getCommonFaceID(const EFAElement3D * other_elem) const
1527 : {
1528 : std::vector<unsigned int> face_id;
1529 10897268 : for (unsigned int i = 0; i < _num_faces; ++i)
1530 48895602 : for (unsigned int j = 0; j < other_elem->_num_faces; ++j)
1531 40316538 : if (_faces[i]->equivalent(other_elem->_faces[j]))
1532 : {
1533 317328 : face_id.push_back(i);
1534 : break;
1535 : }
1536 :
1537 2000876 : return face_id;
1538 0 : }
1539 :
1540 : bool
1541 1347038 : EFAElement3D::getCommonEdgeID(const EFAElement3D * other_elem,
1542 : std::vector<std::pair<unsigned int, unsigned int>> & common_ids) const
1543 : {
1544 : // all edges of the other element
1545 : std::set<std::pair<EFANode *, EFANode *>> other_edges;
1546 7060010 : for (const auto k : index_range(other_elem->_faces))
1547 : {
1548 5712972 : const auto & face = *other_elem->_faces[k];
1549 23826348 : for (const auto l : make_range(other_elem->_faces[k]->numEdges()))
1550 : {
1551 : const auto & edge = *face.getEdge(l);
1552 18113376 : other_edges.insert(edge.getSortedNodes());
1553 : }
1554 : }
1555 :
1556 : // loop over all edges of this element
1557 1347038 : common_ids.clear();
1558 7060010 : for (const auto i : index_range(_faces))
1559 23826348 : for (const auto j : make_range(_faces[i]->numEdges()))
1560 : {
1561 18113376 : const auto & edge = *_faces[i]->getEdge(j);
1562 18113376 : const auto edge_nodes = edge.getSortedNodes();
1563 :
1564 : // is this edge contained in the other element?
1565 18113376 : if (edge.isEmbeddedPermanent() || other_edges.count(edge_nodes) == 0)
1566 17439528 : continue;
1567 :
1568 673848 : common_ids.emplace_back(i, j);
1569 : }
1570 :
1571 1347038 : return common_ids.size() > 0;
1572 : }
1573 :
1574 : unsigned int
1575 17308 : EFAElement3D::getNeighborFaceNodeID(unsigned int face_id,
1576 : unsigned int node_id,
1577 : EFAElement3D * neighbor_elem) const
1578 : {
1579 : // get the corresponding node_id on the corresponding face of neighbor_elem
1580 : bool found_id = false;
1581 : unsigned int neigh_face_node_id;
1582 17308 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1583 17308 : if (common_face_id == face_id)
1584 : {
1585 17308 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1586 17308 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1587 40282 : for (unsigned int i = 0; i < neigh_face->numNodes(); ++i)
1588 : {
1589 40282 : if (_faces[face_id]->getNode(node_id) == neigh_face->getNode(i))
1590 : {
1591 : neigh_face_node_id = i;
1592 : found_id = true;
1593 : break;
1594 : }
1595 : }
1596 : }
1597 : else
1598 0 : EFAError("getNeighborFaceNodeID: neighbor_elem is not a neighbor on face_id");
1599 17308 : if (!found_id)
1600 0 : EFAError("getNeighborFaceNodeID: could not find neighbor face node id");
1601 17308 : return neigh_face_node_id;
1602 : }
1603 :
1604 : unsigned int
1605 5836 : EFAElement3D::getNeighborFaceInteriorNodeID(unsigned int face_id,
1606 : unsigned int node_id,
1607 : EFAElement3D * neighbor_elem) const
1608 : {
1609 : // get the corresponding node_id on the corresponding face of neighbor_elem
1610 : bool found_id = false;
1611 : unsigned int neigh_face_node_id;
1612 5836 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1613 5836 : if (common_face_id == face_id)
1614 : {
1615 5836 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1616 5836 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1617 :
1618 13552 : for (unsigned int i = 0; i < _num_interior_face_nodes; ++i)
1619 : {
1620 13552 : if (_faces[face_id]->getInteriorFaceNode(node_id) == neigh_face->getInteriorFaceNode(i))
1621 : {
1622 : neigh_face_node_id = i;
1623 : found_id = true;
1624 : break;
1625 : }
1626 : }
1627 : }
1628 : else
1629 0 : EFAError("getNeighborFaceNodeID: neighbor_elem is not a neighbor on face_id");
1630 5836 : if (!found_id)
1631 0 : EFAError("getNeighborFaceNodeID: could not find neighbor face node id");
1632 5836 : return neigh_face_node_id;
1633 : }
1634 :
1635 : unsigned int
1636 39292 : EFAElement3D::getNeighborFaceEdgeID(unsigned int face_id,
1637 : unsigned int edge_id,
1638 : EFAElement3D * neighbor_elem) const
1639 : {
1640 : // get the corresponding edge_id on the corresponding face of neighbor_elem
1641 : bool found_id = false;
1642 : unsigned int neigh_face_edge_id;
1643 39292 : unsigned int common_face_id = getNeighborIndex(neighbor_elem);
1644 39292 : if (common_face_id == face_id)
1645 : {
1646 39292 : unsigned int neigh_face_id = neighbor_elem->getNeighborIndex(this);
1647 39292 : EFAFace * neigh_face = neighbor_elem->getFace(neigh_face_id);
1648 83108 : for (unsigned int i = 0; i < neigh_face->numEdges(); ++i)
1649 : {
1650 83108 : if (_faces[face_id]->getEdge(edge_id)->equivalent(*neigh_face->getEdge(i)))
1651 : {
1652 : neigh_face_edge_id = i;
1653 : found_id = true;
1654 : break;
1655 : }
1656 : }
1657 : }
1658 : else
1659 0 : EFAError("getNeighborFaceEdgeID: neighbor_elem is not a neighbor on face_id");
1660 39292 : if (!found_id)
1661 0 : EFAError("getNeighborFaceEdgeID: could not find neighbor face edge id");
1662 39292 : return neigh_face_edge_id;
1663 : }
1664 :
1665 : void
1666 46034 : EFAElement3D::findFacesAdjacentToFaces()
1667 : {
1668 46034 : _faces_adjacent_to_faces.clear();
1669 276942 : for (unsigned int i = 0; i < _faces.size(); ++i)
1670 : {
1671 230908 : std::vector<EFAFace *> face_adjacents(_faces[i]->numEdges(), nullptr);
1672 1435172 : for (unsigned int j = 0; j < _faces.size(); ++j)
1673 : {
1674 1204264 : if (_faces[j] != _faces[i] && _faces[i]->isAdjacent(_faces[j]))
1675 : {
1676 833040 : unsigned int adj_edge = _faces[i]->adjacentCommonEdge(_faces[j]);
1677 833040 : face_adjacents[adj_edge] = _faces[j];
1678 : }
1679 : }
1680 230908 : _faces_adjacent_to_faces.push_back(face_adjacents);
1681 230908 : }
1682 46034 : }
1683 :
1684 : EFAFace *
1685 99393 : EFAElement3D::getAdjacentFace(unsigned int face_id, unsigned int edge_id) const
1686 : {
1687 99393 : return _faces_adjacent_to_faces[face_id][edge_id];
1688 : }
1689 :
1690 : EFAFace *
1691 129073 : EFAElement3D::getFragmentFace(unsigned int frag_id, unsigned int face_id) const
1692 : {
1693 129073 : if (frag_id < _fragments.size())
1694 129073 : return _fragments[frag_id]->getFace(face_id);
1695 : else
1696 0 : EFAError("frag_id out of bounds in getFragmentFace");
1697 : }
1698 :
1699 : std::set<EFANode *>
1700 14284 : EFAElement3D::getPhantomNodeOnFace(unsigned int face_id) const
1701 : {
1702 : std::set<EFANode *> phantom_nodes;
1703 14284 : if (_fragments.size() > 0)
1704 : {
1705 65500 : for (unsigned int j = 0; j < _faces[face_id]->numNodes(); ++j) // loop ove 2 edge nodes
1706 : {
1707 : bool node_in_frag = false;
1708 74640 : for (unsigned int k = 0; k < _fragments.size(); ++k)
1709 : {
1710 51216 : if (_fragments[k]->containsNode(_faces[face_id]->getNode(j)))
1711 : {
1712 : node_in_frag = true;
1713 : break;
1714 : }
1715 : }
1716 51216 : if (!node_in_frag)
1717 23424 : phantom_nodes.insert(_faces[face_id]->getNode(j));
1718 : }
1719 : }
1720 14284 : return phantom_nodes;
1721 : }
1722 :
1723 : bool
1724 19484 : EFAElement3D::getFragmentFaceID(unsigned int elem_face_id, unsigned int & frag_face_id) const
1725 : {
1726 : // find the fragment face that is contained by given element edge
1727 : // N.B. if the elem edge contains two frag edges, this method will only return
1728 : // the first frag edge ID
1729 : bool frag_face_found = false;
1730 19484 : frag_face_id = std::numeric_limits<unsigned int>::max();
1731 19484 : if (_fragments.size() == 1)
1732 : {
1733 1344 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1734 : {
1735 1344 : if (_faces[elem_face_id]->containsFace(_fragments[0]->getFace(j)))
1736 : {
1737 384 : frag_face_id = j;
1738 : frag_face_found = true;
1739 384 : break;
1740 : }
1741 : }
1742 : }
1743 19484 : return frag_face_found;
1744 : }
1745 :
1746 : bool
1747 9806 : EFAElement3D::getFragmentFaceEdgeID(unsigned int ElemFaceID,
1748 : unsigned int ElemFaceEdgeID,
1749 : unsigned int & FragFaceID,
1750 : unsigned int & FragFaceEdgeID) const
1751 : {
1752 : // Purpose: given an edge of an elem face, find which frag face's edge it contains
1753 : bool frag_edge_found = false;
1754 9806 : FragFaceID = FragFaceEdgeID = std::numeric_limits<unsigned int>::max();
1755 9806 : if (getFragmentFaceID(ElemFaceID, FragFaceID))
1756 : {
1757 256 : EFAEdge * elem_edge = _faces[ElemFaceID]->getEdge(ElemFaceEdgeID);
1758 256 : EFAFace * frag_face = getFragmentFace(0, FragFaceID);
1759 640 : for (unsigned int i = 0; i < frag_face->numEdges(); ++i)
1760 : {
1761 640 : if (elem_edge->containsEdge(*frag_face->getEdge(i)))
1762 : {
1763 256 : FragFaceEdgeID = i;
1764 : frag_edge_found = true;
1765 256 : break;
1766 : }
1767 : }
1768 : }
1769 9806 : return frag_edge_found;
1770 : }
1771 :
1772 : bool
1773 9678 : EFAElement3D::isPhysicalEdgeCut(unsigned int ElemFaceID,
1774 : unsigned int ElemFaceEdgeID,
1775 : double position) const
1776 : {
1777 9678 : unsigned int FragFaceID, FragFaceEdgeID = std::numeric_limits<unsigned int>::max();
1778 : bool is_in_real = false;
1779 9678 : if (_fragments.size() == 0)
1780 : {
1781 : is_in_real = true;
1782 : }
1783 128 : else if (getFragmentFaceEdgeID(ElemFaceID, ElemFaceEdgeID, FragFaceID, FragFaceEdgeID))
1784 : {
1785 128 : EFAEdge * elem_edge = _faces[ElemFaceID]->getEdge(ElemFaceEdgeID);
1786 128 : EFAEdge * frag_edge = getFragmentFace(0, FragFaceID)->getEdge(FragFaceEdgeID);
1787 : double xi[2] = {-1.0, -1.0}; // relative coords of two frag edge nodes
1788 128 : xi[0] = elem_edge->distanceFromNode1(frag_edge->getNode(0));
1789 128 : xi[1] = elem_edge->distanceFromNode1(frag_edge->getNode(1));
1790 128 : if ((position - xi[0]) * (position - xi[1]) <
1791 : 0.0) // the cut to be added is within the real part of the edge
1792 : is_in_real = true;
1793 : }
1794 9678 : return is_in_real;
1795 : }
1796 :
1797 : bool
1798 15746 : EFAElement3D::isFacePhantom(unsigned int face_id) const
1799 : {
1800 : bool is_phantom = false;
1801 15746 : if (_fragments.size() > 0)
1802 : {
1803 : bool contains_frag_face = false;
1804 11550 : for (unsigned int i = 0; i < _fragments.size(); ++i)
1805 : {
1806 31474 : for (unsigned int j = 0; j < _fragments[i]->numFaces(); ++j)
1807 : {
1808 31474 : if (_faces[face_id]->containsFace(_fragments[i]->getFace(j)))
1809 : {
1810 : contains_frag_face = true;
1811 : break;
1812 : }
1813 : }
1814 11550 : if (contains_frag_face)
1815 : break;
1816 : }
1817 11550 : if (!contains_frag_face)
1818 : is_phantom = true;
1819 : }
1820 15746 : return is_phantom;
1821 : }
1822 :
1823 : unsigned int
1824 1395897 : EFAElement3D::numFaceNeighbors(unsigned int face_id) const
1825 : {
1826 1395897 : return _face_neighbors[face_id].size();
1827 : }
1828 :
1829 : unsigned int
1830 69794 : EFAElement3D::numEdgeNeighbors(unsigned int face_id, unsigned int edge_id) const
1831 : {
1832 69794 : return _face_edge_neighbors[face_id][edge_id].size();
1833 : }
1834 :
1835 : EFAElement3D *
1836 677687 : EFAElement3D::getFaceNeighbor(unsigned int face_id, unsigned int neighbor_id) const
1837 : {
1838 677687 : if (_face_neighbors[face_id][0] != nullptr && neighbor_id < _face_neighbors[face_id].size())
1839 677687 : return _face_neighbors[face_id][neighbor_id];
1840 : else
1841 0 : EFAError("edge neighbor does not exist");
1842 : }
1843 :
1844 : EFAElement3D *
1845 37900 : EFAElement3D::getEdgeNeighbor(unsigned int face_id,
1846 : unsigned int edge_id,
1847 : unsigned int neighbor_id) const
1848 : {
1849 37900 : if (neighbor_id < _face_edge_neighbors[face_id][edge_id].size())
1850 37900 : return _face_edge_neighbors[face_id][edge_id][neighbor_id];
1851 : else
1852 0 : EFAError("edge neighbor does not exist");
1853 : }
1854 :
1855 : bool
1856 51333 : EFAElement3D::fragmentHasTipFaces() const
1857 : {
1858 : bool has_tip_faces = false;
1859 51333 : if (_fragments.size() == 1)
1860 : {
1861 82970 : for (unsigned int i = 0; i < _num_faces; ++i)
1862 : {
1863 : unsigned int num_frag_faces = 0; // count how many fragment edges this element edge contains
1864 71206 : if (_faces[i]->hasIntersection())
1865 : {
1866 303040 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1867 : {
1868 256536 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1869 50160 : num_frag_faces += 1;
1870 : } // j
1871 46504 : if (num_frag_faces == 2)
1872 : {
1873 : has_tip_faces = true;
1874 : break;
1875 : }
1876 : }
1877 : } // i
1878 : }
1879 51333 : return has_tip_faces;
1880 : }
1881 :
1882 : std::vector<unsigned int>
1883 0 : EFAElement3D::getTipFaceIDs() const
1884 : {
1885 : // if this element is a crack tip element, returns the crack tip faces' ID
1886 : std::vector<unsigned int> tip_face_id;
1887 0 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1888 : {
1889 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1890 : {
1891 : unsigned int num_frag_faces = 0; // count how many fragment faces this element edge contains
1892 0 : if (_faces[i]->hasIntersection())
1893 : {
1894 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1895 : {
1896 0 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1897 0 : num_frag_faces += 1;
1898 : } // j
1899 0 : if (num_frag_faces == 2) // element face contains two fragment edges
1900 0 : tip_face_id.push_back(i);
1901 : }
1902 : } // i
1903 : }
1904 0 : return tip_face_id;
1905 0 : }
1906 :
1907 : std::set<EFANode *>
1908 0 : EFAElement3D::getTipEmbeddedNodes() const
1909 : {
1910 : // if this element is a crack tip element, returns the crack tip edge's ID
1911 : std::set<EFANode *> tip_emb;
1912 0 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1913 : {
1914 0 : for (unsigned int i = 0; i < _num_faces; ++i)
1915 : {
1916 : std::vector<EFAFace *> frag_faces; // count how many fragment edges this element edge contains
1917 0 : if (_faces[i]->hasIntersection())
1918 : {
1919 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1920 : {
1921 0 : if (_faces[i]->containsFace(_fragments[0]->getFace(j)))
1922 0 : frag_faces.push_back(_fragments[0]->getFace(j));
1923 : } // j
1924 0 : if (frag_faces.size() == 2) // element edge contains two fragment edges
1925 : {
1926 0 : unsigned int edge_id = frag_faces[0]->adjacentCommonEdge(frag_faces[1]);
1927 0 : tip_emb.insert(frag_faces[0]->getEdge(edge_id)->getNode(0));
1928 0 : tip_emb.insert(frag_faces[0]->getEdge(edge_id)->getNode(1));
1929 : }
1930 : }
1931 0 : } // i
1932 : }
1933 0 : return tip_emb;
1934 : }
1935 :
1936 : bool
1937 9678 : EFAElement3D::faceContainsTip(unsigned int face_id) const
1938 : {
1939 : bool contain_tip = false;
1940 9678 : if (_fragments.size() == 1)
1941 : {
1942 : unsigned int num_frag_faces = 0; // count how many fragment faces this element face contains
1943 128 : if (_faces[face_id]->hasIntersection())
1944 : {
1945 0 : for (unsigned int j = 0; j < _fragments[0]->numFaces(); ++j)
1946 : {
1947 0 : if (_faces[face_id]->containsFace(_fragments[0]->getFace(j)))
1948 0 : num_frag_faces += 1;
1949 : } // j
1950 0 : if (num_frag_faces == 2)
1951 : contain_tip = true;
1952 : }
1953 : }
1954 9678 : return contain_tip;
1955 : }
1956 :
1957 : bool
1958 9678 : EFAElement3D::fragmentFaceAlreadyCut(unsigned int ElemFaceID) const
1959 : {
1960 : // when marking cuts, check if the corresponding frag face already has been cut
1961 : bool has_cut = false;
1962 9678 : if (faceContainsTip(ElemFaceID))
1963 : has_cut = true;
1964 : else
1965 : {
1966 9678 : unsigned int FragFaceID = std::numeric_limits<unsigned int>::max();
1967 9678 : if (getFragmentFaceID(ElemFaceID, FragFaceID))
1968 : {
1969 128 : EFAFace * frag_face = getFragmentFace(0, FragFaceID);
1970 128 : if (frag_face->hasIntersection())
1971 : has_cut = true;
1972 : }
1973 : }
1974 9678 : return has_cut;
1975 : }
1976 :
1977 : void
1978 99393 : EFAElement3D::addFaceEdgeCut(unsigned int face_id,
1979 : unsigned int edge_id,
1980 : double position,
1981 : EFANode * embedded_node,
1982 : std::map<unsigned int, EFANode *> & EmbeddedNodes,
1983 : bool add_to_neighbor,
1984 : bool add_to_adjacent)
1985 : {
1986 : // Purpose: add intersection on Edge edge_id of Face face_id
1987 99393 : EFANode * local_embedded = nullptr;
1988 99393 : EFAEdge * cut_edge = _faces[face_id]->getEdge(edge_id);
1989 99393 : EFANode * edge_node1 = cut_edge->getNode(0);
1990 99393 : if (embedded_node) // use the existing embedded node if it was passed in
1991 67499 : local_embedded = embedded_node;
1992 :
1993 : // get adjacent face info
1994 99393 : EFAFace * adj_face = getAdjacentFace(face_id, edge_id);
1995 99393 : unsigned int adj_face_id = getFaceID(adj_face);
1996 99393 : unsigned int adj_edge_id = adj_face->adjacentCommonEdge(_faces[face_id]);
1997 :
1998 : // check if cut has already been added to this face edge
1999 : bool cut_exist = false;
2000 :
2001 99393 : if (cut_edge->hasIntersectionAtPosition(position, edge_node1))
2002 : {
2003 89715 : unsigned int emb_id = cut_edge->getEmbeddedNodeIndex(position, edge_node1);
2004 89715 : EFANode * old_emb = cut_edge->getEmbeddedNode(emb_id);
2005 89715 : if (embedded_node && embedded_node != old_emb)
2006 0 : EFAError("Attempting to add edge intersection when one already exists with different node.",
2007 : " elem: ",
2008 : _id,
2009 : " edge: ",
2010 : edge_id,
2011 : " position: ",
2012 : position);
2013 89715 : local_embedded = old_emb;
2014 : cut_exist = true;
2015 : }
2016 :
2017 19356 : if (!cut_exist && !fragmentFaceAlreadyCut(face_id) &&
2018 9678 : isPhysicalEdgeCut(face_id, edge_id, position))
2019 : {
2020 : // check if cut has already been added to the neighbor edges
2021 9678 : checkNeighborFaceCut(face_id, edge_id, position, edge_node1, embedded_node, local_embedded);
2022 9678 : checkNeighborFaceCut(
2023 : adj_face_id, adj_edge_id, position, edge_node1, embedded_node, local_embedded);
2024 :
2025 9678 : if (!local_embedded) // need to create new embedded node
2026 : {
2027 1766 : unsigned int new_node_id = Efa::getNewID(EmbeddedNodes);
2028 1766 : local_embedded = new EFANode(new_node_id, EFANode::N_CATEGORY_EMBEDDED);
2029 1766 : EmbeddedNodes.insert(std::make_pair(new_node_id, local_embedded));
2030 : }
2031 :
2032 : // add to elem face edge
2033 9678 : cut_edge->addIntersection(position, local_embedded, edge_node1);
2034 9678 : if (cut_edge->numEmbeddedNodes() > 2)
2035 0 : EFAError("element edge can't have >2 embedded nodes");
2036 :
2037 : // cut fragment faces, which is an essential addition to other code in this method that cuts
2038 : // element faces
2039 : unsigned int FragFaceID;
2040 : unsigned int FragFaceEdgeID;
2041 9678 : if (getFragmentFaceEdgeID(face_id, edge_id, FragFaceID, FragFaceEdgeID))
2042 : {
2043 128 : EFAEdge * elem_edge = _faces[face_id]->getEdge(edge_id);
2044 128 : EFAEdge * frag_edge = getFragmentFace(0, FragFaceID)->getEdge(FragFaceEdgeID);
2045 : double xi[2] = {-1.0, -1.0}; // relative coords of two frag edge nodes
2046 128 : xi[0] = elem_edge->distanceFromNode1(frag_edge->getNode(0));
2047 128 : xi[1] = elem_edge->distanceFromNode1(frag_edge->getNode(1));
2048 128 : double frag_pos = (position - xi[0]) / (xi[1] - xi[0]);
2049 128 : EFANode * frag_edge_node1 = frag_edge->getNode(0);
2050 :
2051 128 : if (!frag_edge->hasIntersection())
2052 128 : frag_edge->addIntersection(frag_pos, local_embedded, frag_edge_node1);
2053 : }
2054 :
2055 : // add to adjacent face edge
2056 9678 : if (add_to_adjacent)
2057 : {
2058 4839 : double adj_pos = 1.0 - position;
2059 4839 : addFaceEdgeCut(
2060 : adj_face_id, adj_edge_id, adj_pos, local_embedded, EmbeddedNodes, false, false);
2061 : }
2062 : }
2063 :
2064 : // add cut to neighbor face edge
2065 99393 : if (add_to_neighbor)
2066 : {
2067 56654 : for (unsigned int en_iter = 0; en_iter < numFaceNeighbors(face_id); ++en_iter)
2068 : {
2069 24760 : EFAElement3D * face_neighbor = getFaceNeighbor(face_id, en_iter);
2070 24760 : unsigned int neigh_face_id = face_neighbor->getNeighborIndex(this);
2071 24760 : unsigned neigh_edge_id = getNeighborFaceEdgeID(face_id, edge_id, face_neighbor);
2072 24760 : double neigh_pos = 1.0 - position; // get emb node's postion on neighbor edge
2073 24760 : face_neighbor->addFaceEdgeCut(
2074 : neigh_face_id, neigh_edge_id, neigh_pos, local_embedded, EmbeddedNodes, false, true);
2075 : }
2076 :
2077 69794 : for (unsigned int en_iter = 0; en_iter < numEdgeNeighbors(face_id, edge_id); ++en_iter)
2078 : {
2079 37900 : EFAElement3D * edge_neighbor = getEdgeNeighbor(face_id, edge_id, en_iter);
2080 : unsigned int neigh_face_id, neigh_edge_id;
2081 37900 : getNeighborEdgeIndex(edge_neighbor, face_id, edge_id, neigh_face_id, neigh_edge_id);
2082 :
2083 : // Check the ordering of the node and the assign the position
2084 : double neigh_pos;
2085 75800 : if (_faces[face_id]->getEdge(edge_id)->getNode(0) ==
2086 37900 : edge_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id)->getNode(0))
2087 : neigh_pos = position;
2088 37818 : else if (_faces[face_id]->getEdge(edge_id)->getNode(1) ==
2089 18909 : edge_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id)->getNode(0))
2090 18909 : neigh_pos = 1.0 - position;
2091 : else
2092 0 : EFAError("The EFANodes on commaon edge are not matched.");
2093 :
2094 37900 : edge_neighbor->addFaceEdgeCut(
2095 : neigh_face_id, neigh_edge_id, neigh_pos, local_embedded, EmbeddedNodes, false, true);
2096 : }
2097 : } // If add_to_neighbor required
2098 99393 : }
2099 :
2100 : void
2101 0 : EFAElement3D::addFragFaceEdgeCut(unsigned int /*frag_face_id*/,
2102 : unsigned int /*frag_edge_id*/,
2103 : double /*position*/,
2104 : std::map<unsigned int, EFANode *> & /*EmbeddedNodes*/,
2105 : bool /*add_to_neighbor*/,
2106 : bool /*add_to_adjacent*/)
2107 : {
2108 : // TODO: mark frag face edges
2109 : // also need to check if cut has been added to this frag face edge or neighbor edge of adjacent
2110 : // face
2111 0 : }
2112 :
2113 : void
2114 19356 : EFAElement3D::checkNeighborFaceCut(unsigned int face_id,
2115 : unsigned int edge_id,
2116 : double position,
2117 : EFANode * from_node,
2118 : EFANode * embedded_node,
2119 : EFANode *& local_embedded)
2120 : {
2121 : // N.B. this is important. We are checking if the corresponding edge of the neighbor face or of
2122 : // the adjacent
2123 : // face's neighbor face has a cut at the same position. If so, use the existing embedded node as
2124 : // local_embedded
2125 33888 : for (unsigned int en_iter = 0; en_iter < numFaceNeighbors(face_id); ++en_iter)
2126 : {
2127 14532 : EFAElement3D * face_neighbor = getFaceNeighbor(face_id, en_iter);
2128 14532 : unsigned int neigh_face_id = face_neighbor->getNeighborIndex(this);
2129 14532 : unsigned neigh_edge_id = getNeighborFaceEdgeID(face_id, edge_id, face_neighbor);
2130 14532 : EFAEdge * neigh_edge = face_neighbor->getFace(neigh_face_id)->getEdge(neigh_edge_id);
2131 :
2132 14532 : if (neigh_edge->hasIntersectionAtPosition(position, from_node))
2133 : {
2134 7332 : unsigned int emb_id = neigh_edge->getEmbeddedNodeIndex(position, from_node);
2135 7332 : EFANode * old_emb = neigh_edge->getEmbeddedNode(emb_id);
2136 :
2137 7332 : if (embedded_node && embedded_node != old_emb)
2138 0 : EFAError(
2139 : "attempting to add edge intersection when one already exists with different node.");
2140 7332 : if (local_embedded && local_embedded != old_emb)
2141 0 : EFAError("attempting to assign contradictory pointer to local_embedded.");
2142 :
2143 7332 : local_embedded = old_emb;
2144 : }
2145 : } // en_iter
2146 19356 : }
2147 :
2148 : void
2149 0 : EFAElement3D::mapParametricCoordinateFrom2DTo3D(unsigned int face_id,
2150 : std::vector<double> & xi_2d,
2151 : std::vector<double> & xi_3d) const
2152 : {
2153 : // given the coords of a point in a 2D face, translate it to 3D parametric coords
2154 0 : xi_3d.resize(3, 0.0);
2155 0 : if (_num_faces == 6)
2156 : {
2157 : if (face_id == 0)
2158 : {
2159 0 : xi_3d[0] = xi_2d[1];
2160 0 : xi_3d[1] = xi_2d[0];
2161 0 : xi_3d[2] = -1.0;
2162 : }
2163 : else if (face_id == 1)
2164 : {
2165 0 : xi_3d[0] = xi_2d[0];
2166 0 : xi_3d[1] = -1.0;
2167 0 : xi_3d[2] = xi_2d[1];
2168 : }
2169 : else if (face_id == 2)
2170 : {
2171 0 : xi_3d[0] = 1.0;
2172 0 : xi_3d[1] = xi_2d[0];
2173 0 : xi_3d[2] = xi_2d[1];
2174 : }
2175 : else if (face_id == 3)
2176 : {
2177 0 : xi_3d[0] = -xi_2d[0];
2178 0 : xi_3d[1] = 1.0;
2179 0 : xi_3d[2] = xi_2d[1];
2180 : }
2181 : else if (face_id == 4)
2182 : {
2183 0 : xi_3d[0] = -1.0;
2184 0 : xi_3d[1] = -xi_2d[0];
2185 0 : xi_3d[2] = xi_2d[1];
2186 : }
2187 : else if (face_id == 5)
2188 : {
2189 0 : xi_3d[0] = xi_2d[0];
2190 0 : xi_3d[1] = xi_2d[1];
2191 0 : xi_3d[2] = 1.0;
2192 : }
2193 : else
2194 0 : EFAError("face_id out of bounds");
2195 : }
2196 0 : else if (_num_faces == 4)
2197 : {
2198 0 : if (face_id == 0)
2199 : {
2200 0 : xi_3d[0] = xi_2d[0];
2201 0 : xi_3d[1] = xi_2d[1];
2202 0 : xi_3d[2] = 0.0;
2203 : }
2204 0 : else if (face_id == 1)
2205 : {
2206 0 : xi_3d[0] = 0.0;
2207 0 : xi_3d[1] = xi_2d[0];
2208 0 : xi_3d[2] = xi_2d[1];
2209 : }
2210 0 : else if (face_id == 2)
2211 : {
2212 0 : xi_3d[0] = xi_2d[1];
2213 0 : xi_3d[1] = 0.0;
2214 0 : xi_3d[2] = xi_2d[0];
2215 : }
2216 0 : else if (face_id == 3)
2217 : {
2218 0 : xi_3d[0] = xi_2d[0];
2219 0 : xi_3d[1] = xi_2d[2];
2220 0 : xi_3d[2] = xi_2d[1];
2221 : }
2222 : else
2223 0 : EFAError("face_id out of bounds");
2224 : }
2225 : else
2226 0 : EFAError("unknown element for 3D");
2227 0 : }
2228 :
2229 : std::vector<EFANode *>
2230 0 : EFAElement3D::getCommonNodes(const EFAElement3D * other_elem) const
2231 : {
2232 : std::set<EFANode *> e1nodes(_nodes.begin(),
2233 0 : _nodes.begin() + _num_vertices); // only account for corner nodes
2234 : std::set<EFANode *> e2nodes(other_elem->_nodes.begin(),
2235 0 : other_elem->_nodes.begin() + _num_vertices);
2236 0 : std::vector<EFANode *> common_nodes = Efa::getCommonElems(e1nodes, e2nodes);
2237 0 : return common_nodes;
2238 : }
|
