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 "EFAElement2D.h"
11 :
12 : #include <iomanip>
13 :
14 : #include "EFAFaceNode.h"
15 : #include "EFANode.h"
16 : #include "EFAEdge.h"
17 : #include "EFAFace.h"
18 : #include "EFAFragment2D.h"
19 : #include "EFAFuncs.h"
20 : #include "EFAError.h"
21 : #include "XFEMFuncs.h"
22 :
23 394707 : EFAElement2D::EFAElement2D(unsigned int eid, unsigned int n_nodes) : EFAElement(eid, n_nodes)
24 : {
25 394707 : if (n_nodes == 4 || n_nodes == 8 || n_nodes == 9)
26 370523 : _num_edges = 4;
27 24184 : else if (n_nodes == 3 || n_nodes == 6 || n_nodes == 7)
28 24184 : _num_edges = 3;
29 : else
30 0 : EFAError(
31 : "In EFAelement2D the supported element types are QUAD4, QUAD8, QUAD9, TRI3, TRI6 and TRI7");
32 394707 : setLocalCoordinates();
33 394707 : _edges = std::vector<EFAEdge *>(_num_edges, nullptr);
34 394707 : _edge_neighbors = std::vector<std::vector<EFAElement2D *>>(_num_edges, {nullptr});
35 394707 : }
36 :
37 9298 : EFAElement2D::EFAElement2D(const EFAElement2D * from_elem, bool convert_to_local)
38 9298 : : EFAElement(from_elem->_id, from_elem->_num_nodes),
39 9298 : _num_edges(from_elem->_num_edges),
40 9298 : _edges(_num_edges, nullptr),
41 9298 : _edge_neighbors(_num_edges, {nullptr})
42 : {
43 9298 : if (convert_to_local)
44 : {
45 : // build local nodes from global nodes
46 50134 : for (unsigned int i = 0; i < _num_nodes; ++i)
47 : {
48 40968 : if (from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT ||
49 40968 : from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_TEMP ||
50 132 : from_elem->_nodes[i]->category() == EFANode::N_CATEGORY_EMBEDDED_PERMANENT)
51 : {
52 40836 : _nodes[i] = from_elem->createLocalNodeFromGlobalNode(from_elem->_nodes[i]);
53 40836 : _local_nodes.push_back(_nodes[i]); // convenient to delete local nodes
54 : }
55 : else
56 0 : EFAError("In EFAelement2D ",
57 : from_elem->id(),
58 : " the copy constructor must have from_elem w/ global nodes. node: ",
59 : i,
60 : " category: ",
61 : from_elem->_nodes[i]->category());
62 : }
63 :
64 : // copy edges, fragments and interior nodes from from_elem
65 45370 : for (unsigned int i = 0; i < _num_edges; ++i)
66 36072 : _edges[i] = new EFAEdge(*from_elem->_edges[i]);
67 18596 : for (unsigned int i = 0; i < from_elem->_fragments.size(); ++i)
68 9298 : _fragments.push_back(new EFAFragment2D(this, true, from_elem, i));
69 9346 : for (unsigned int i = 0; i < from_elem->_interior_nodes.size(); ++i)
70 48 : _interior_nodes.push_back(new EFAFaceNode(*from_elem->_interior_nodes[i]));
71 :
72 : // replace all global nodes with local nodes
73 50134 : for (unsigned int i = 0; i < _num_nodes; ++i)
74 : {
75 40836 : if (_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
76 40836 : switchNode(
77 : _nodes[i],
78 : from_elem->_nodes[i],
79 : false); // when save to _cut_elem_map, the EFAelement is not a child of any parent
80 : else
81 0 : EFAError("In EFAelement2D copy constructor this elem's nodes must be local");
82 : }
83 :
84 9298 : _local_node_coor = from_elem->_local_node_coor;
85 : }
86 : else
87 0 : EFAError("this EFAelement2D constructor only converts global nodes to local nodes");
88 9298 : }
89 :
90 0 : EFAElement2D::EFAElement2D(const EFAFace * from_face)
91 : : EFAElement(0, from_face->numNodes()),
92 0 : _num_edges(from_face->numEdges()),
93 0 : _edges(_num_edges, nullptr),
94 0 : _edge_neighbors(_num_edges, {nullptr})
95 : {
96 0 : for (unsigned int i = 0; i < _num_nodes; ++i)
97 0 : _nodes[i] = from_face->getNode(i);
98 0 : for (unsigned int i = 0; i < _num_edges; ++i)
99 0 : _edges[i] = new EFAEdge(*from_face->getEdge(i));
100 0 : for (unsigned int i = 0; i < from_face->numInteriorNodes(); ++i)
101 0 : _interior_nodes.push_back(new EFAFaceNode(*from_face->getInteriorNode(i)));
102 0 : }
103 :
104 798312 : EFAElement2D::~EFAElement2D()
105 : {
106 466748 : for (unsigned int i = 0; i < _fragments.size(); ++i)
107 : {
108 62943 : if (_fragments[i])
109 : {
110 62943 : delete _fragments[i];
111 62943 : _fragments[i] = nullptr;
112 : }
113 : }
114 1993721 : for (unsigned int i = 0; i < _edges.size(); ++i)
115 : {
116 1589916 : if (_edges[i])
117 : {
118 1589916 : delete _edges[i];
119 1589916 : _edges[i] = nullptr;
120 : }
121 : }
122 404081 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
123 : {
124 276 : if (_interior_nodes[i])
125 : {
126 276 : delete _interior_nodes[i];
127 276 : _interior_nodes[i] = nullptr;
128 : }
129 : }
130 444641 : for (unsigned int i = 0; i < _local_nodes.size(); ++i)
131 : {
132 40836 : if (_local_nodes[i])
133 : {
134 40836 : delete _local_nodes[i];
135 40836 : _local_nodes[i] = nullptr;
136 : }
137 : }
138 798312 : }
139 :
140 : void
141 394707 : EFAElement2D::setLocalCoordinates()
142 : {
143 394707 : if (_num_edges == 4)
144 : {
145 : /*
146 : 3 6 2
147 : QUAD9(QUAD8): o-----o-----o
148 : | |
149 : | 8 |
150 : 7 o o o 5
151 : | |
152 : | |
153 : o-----o-----o
154 : 0 4 1
155 :
156 : */
157 370523 : _local_node_coor.resize(_num_nodes);
158 370523 : _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
159 370523 : _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
160 370523 : _local_node_coor[2] = EFAPoint(1.0, 1.0, 0.0);
161 370523 : _local_node_coor[3] = EFAPoint(0.0, 1.0, 0.0);
162 :
163 370523 : if (_num_nodes > 4)
164 : {
165 8496 : _local_node_coor[4] = EFAPoint(0.5, 0.0, 0.0);
166 8496 : _local_node_coor[5] = EFAPoint(1.0, 0.5, 0.0);
167 8496 : _local_node_coor[6] = EFAPoint(0.5, 1.0, 0.0);
168 8496 : _local_node_coor[7] = EFAPoint(0.0, 0.5, 0.0);
169 : }
170 :
171 370523 : if (_num_nodes > 8)
172 4320 : _local_node_coor[8] = EFAPoint(0.5, 0.5, 0.0);
173 : }
174 : else
175 : {
176 : /*
177 : TRI7: 2
178 : o
179 : | \
180 : | \
181 : 5 o o 4
182 : | o \
183 : | 6 \
184 : o-----o-----o
185 : 0 3 1
186 : */
187 24184 : _local_node_coor.resize(_num_nodes);
188 24184 : _local_node_coor[0] = EFAPoint(0.0, 0.0, 0.0);
189 24184 : _local_node_coor[1] = EFAPoint(1.0, 0.0, 0.0);
190 24184 : _local_node_coor[2] = EFAPoint(0.0, 1.0, 0.0);
191 :
192 24184 : if (_num_nodes > 3)
193 : {
194 16664 : _local_node_coor[3] = EFAPoint(0.5, 0.0, 0.0);
195 16664 : _local_node_coor[4] = EFAPoint(0.5, 0.5, 0.0);
196 16664 : _local_node_coor[5] = EFAPoint(0.0, 0.5, 0.0);
197 : }
198 :
199 24184 : if (_num_nodes > 6)
200 : {
201 8332 : _local_node_coor[6] = EFAPoint(1 / 3., 1 / 3., 0.0);
202 : }
203 : }
204 394707 : }
205 :
206 : unsigned int
207 2067148 : EFAElement2D::numFragments() const
208 : {
209 2067148 : return _fragments.size();
210 : }
211 :
212 : bool
213 1483694 : EFAElement2D::isPartial() const
214 : {
215 : bool partial = false;
216 1483694 : if (_fragments.size() > 0)
217 : {
218 3525904 : for (unsigned int i = 0; i < _num_edges; ++i)
219 : {
220 : bool node_in_frag = false;
221 4471106 : for (unsigned int j = 0; j < _fragments.size(); ++j)
222 : {
223 3256658 : if (_fragments[j]->containsNode(_nodes[i]))
224 : {
225 : node_in_frag = true;
226 : break;
227 : }
228 : } // j
229 3256658 : if (!node_in_frag)
230 : {
231 : partial = true;
232 : break;
233 : }
234 : } // i
235 : }
236 1483694 : return partial;
237 : }
238 :
239 : void
240 20577 : EFAElement2D::getNonPhysicalNodes(std::set<EFANode *> & non_physical_nodes) const
241 : {
242 : // Any nodes that don't belong to any fragment are non-physical
243 : // First add all nodes in the element to the set
244 114393 : for (unsigned int i = 0; i < _nodes.size(); ++i)
245 93816 : non_physical_nodes.insert(_nodes[i]);
246 :
247 : // Now delete any nodes that are contained in fragments
248 : std::set<EFANode *>::iterator sit;
249 114393 : for (sit = non_physical_nodes.begin(); sit != non_physical_nodes.end();)
250 : {
251 : bool erased = false;
252 146846 : for (unsigned int i = 0; i < _fragments.size(); ++i)
253 : {
254 93816 : if (_fragments[i]->containsNode(*sit))
255 : {
256 40786 : non_physical_nodes.erase(sit++);
257 : erased = true;
258 40786 : break;
259 : }
260 : }
261 93816 : if (!erased)
262 : ++sit;
263 : }
264 20577 : }
265 :
266 : void
267 393757 : EFAElement2D::switchNode(EFANode * new_node, EFANode * old_node, bool descend_to_parent)
268 : {
269 : // We are not switching any embedded nodes here
270 2255089 : for (unsigned int i = 0; i < _num_nodes; ++i)
271 : {
272 1861332 : if (_nodes[i] == old_node)
273 19921 : _nodes[i] = new_node;
274 : }
275 719730 : for (unsigned int i = 0; i < _fragments.size(); ++i)
276 325973 : _fragments[i]->switchNode(new_node, old_node);
277 :
278 1892933 : for (unsigned int i = 0; i < _edges.size(); ++i)
279 1499176 : _edges[i]->switchNode(new_node, old_node);
280 :
281 393757 : if (_parent && descend_to_parent)
282 : {
283 15468 : _parent->switchNode(new_node, old_node, false);
284 141337 : for (unsigned int i = 0; i < _parent->numGeneralNeighbors(); ++i)
285 : {
286 125869 : EFAElement * neigh_elem = _parent->getGeneralNeighbor(i); // generalized neighbor element
287 283202 : for (unsigned int k = 0; k < neigh_elem->numChildren(); ++k)
288 157333 : neigh_elem->getChild(k)->switchNode(new_node, old_node, false);
289 : }
290 : }
291 393757 : }
292 :
293 : void
294 0 : EFAElement2D::switchEmbeddedNode(EFANode * new_emb_node, EFANode * old_emb_node)
295 : {
296 0 : for (unsigned int i = 0; i < _num_edges; ++i)
297 0 : _edges[i]->switchNode(new_emb_node, old_emb_node);
298 0 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
299 0 : _interior_nodes[i]->switchNode(new_emb_node, old_emb_node);
300 0 : for (unsigned int i = 0; i < _fragments.size(); ++i)
301 0 : _fragments[i]->switchNode(new_emb_node, old_emb_node);
302 0 : }
303 :
304 : void
305 4375567 : EFAElement2D::getMasterInfo(EFANode * node,
306 : std::vector<EFANode *> & master_nodes,
307 : std::vector<double> & master_weights) const
308 : {
309 : // Given a EFANode, find the element edge or fragment edge that contains it
310 : // Return its master nodes and weights
311 4375567 : master_nodes.clear();
312 4375567 : master_weights.clear();
313 : bool masters_found = false;
314 11655400 : for (unsigned int i = 0; i < _num_edges; ++i) // check element exterior edges
315 : {
316 11649178 : if (_edges[i]->containsNode(node))
317 : {
318 4369345 : masters_found = _edges[i]->getNodeMasters(node, master_nodes, master_weights);
319 4369345 : if (masters_found)
320 : break;
321 : else
322 0 : EFAError("In getMasterInfo: cannot find master nodes in element edges");
323 : }
324 : }
325 :
326 : if (!masters_found) // check element interior embedded nodes
327 : {
328 6222 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
329 : {
330 6222 : if (_interior_nodes[i]->getNode() == node)
331 : {
332 6222 : std::vector<double> emb_xi(2, 0.0);
333 6222 : emb_xi[0] = _interior_nodes[i]->getParametricCoordinates(0);
334 6222 : emb_xi[1] = _interior_nodes[i]->getParametricCoordinates(1);
335 27942 : for (unsigned int j = 0; j < _num_edges; ++j)
336 : {
337 21720 : master_nodes.push_back(_nodes[j]);
338 21720 : double weight = 0.0;
339 21720 : if (_num_edges == 4)
340 12216 : weight = Efa::linearQuadShape2D(j, emb_xi);
341 9504 : else if (_num_edges == 3)
342 9504 : weight = Efa::linearTriShape2D(j, emb_xi);
343 : else
344 0 : EFAError("unknown 2D element");
345 21720 : master_weights.push_back(weight);
346 : }
347 : masters_found = true;
348 : break;
349 6222 : }
350 : }
351 : }
352 :
353 4369345 : if (!masters_found)
354 0 : EFAError("In EFAelement2D::getMaterInfo, cannot find the given EFAnode");
355 4375567 : }
356 :
357 : unsigned int
358 40761 : EFAElement2D::numInteriorNodes() const
359 : {
360 40761 : return _interior_nodes.size();
361 : }
362 :
363 : bool
364 1384624 : EFAElement2D::overlaysElement(const EFAElement2D * other_elem) const
365 : {
366 : bool overlays = false;
367 :
368 1384624 : if (other_elem->numEdges() != _num_edges)
369 : return overlays;
370 :
371 1384624 : std::vector<EFANode *> common_nodes = getCommonNodes(other_elem);
372 :
373 : // Find indices of common nodes
374 1384624 : if (common_nodes.size() == 2)
375 : {
376 1383520 : std::vector<EFANode *> common_nodes_vec(common_nodes.begin(), common_nodes.end());
377 :
378 1383520 : unsigned int e1n1idx = _num_edges + 1;
379 : unsigned int e1n2idx = _num_edges + 1;
380 6855608 : for (unsigned int i = 0; i < _num_edges; ++i)
381 : {
382 5472088 : if (_nodes[i] == common_nodes_vec[0])
383 : {
384 : e1n1idx = i;
385 : }
386 4088568 : else if (_nodes[i] == common_nodes_vec[1])
387 : {
388 : e1n2idx = i;
389 : }
390 : }
391 :
392 1383520 : if (e1n1idx > _num_edges || e1n2idx > _num_edges)
393 0 : EFAError("in overlays_elem() couldn't find common node");
394 :
395 : bool e1ascend = false;
396 1383520 : unsigned int e1n1idx_plus1(e1n1idx < (_num_edges - 1) ? e1n1idx + 1 : 0);
397 1383520 : unsigned int e1n1idx_minus1(e1n1idx > 0 ? e1n1idx - 1 : _num_edges - 1);
398 1383520 : if (e1n2idx == e1n1idx_plus1)
399 : {
400 : e1ascend = true;
401 : }
402 : else
403 : {
404 691620 : if (e1n2idx != e1n1idx_minus1)
405 0 : EFAError("in overlays_elem() common nodes must be adjacent to each other");
406 : }
407 :
408 : unsigned int e2n1idx = _num_edges + 1;
409 : unsigned int e2n2idx = _num_edges + 1;
410 6855608 : for (unsigned int i = 0; i < _num_edges; ++i)
411 : {
412 5472088 : if (other_elem->getNode(i) == common_nodes_vec[0])
413 : {
414 : e2n1idx = i;
415 : }
416 4088568 : else if (other_elem->getNode(i) == common_nodes_vec[1])
417 : {
418 : e2n2idx = i;
419 : }
420 : }
421 1383520 : if (e2n1idx > other_elem->numNodes() || e2n2idx > other_elem->numNodes())
422 0 : EFAError("in overlays_elem() couldn't find common node");
423 :
424 : bool e2ascend = false;
425 1383520 : unsigned int e2n1idx_plus1(e2n1idx < (_num_edges - 1) ? e2n1idx + 1 : 0);
426 1383520 : unsigned int e2n1idx_minus1(e2n1idx > 0 ? e2n1idx - 1 : _num_edges - 1);
427 1383520 : if (e2n2idx == e2n1idx_plus1)
428 : {
429 : e2ascend = true;
430 : }
431 : else
432 : {
433 691772 : if (e2n2idx != e2n1idx_minus1)
434 0 : EFAError("in overlays_elem() common nodes must be adjacent to each other");
435 : }
436 :
437 : // if indices both ascend or descend, they overlay
438 1383520 : if ((e1ascend && e2ascend) || (!e1ascend && !e2ascend))
439 : {
440 : overlays = true;
441 : }
442 1383520 : }
443 1104 : else if (common_nodes.size() > 2)
444 : {
445 : // TODO: We probably need more error checking here.
446 : overlays = true;
447 : }
448 : return overlays;
449 1384624 : }
450 :
451 : unsigned int
452 82343 : EFAElement2D::getNeighborIndex(const EFAElement * neighbor_elem) const
453 : {
454 211970 : for (unsigned int i = 0; i < _num_edges; ++i)
455 343549 : for (unsigned int j = 0; j < _edge_neighbors[i].size(); ++j)
456 213922 : if (_edge_neighbors[i][j] == neighbor_elem)
457 82343 : return i;
458 :
459 0 : EFAError(
460 : "in get_neighbor_index() element: ", _id, " does not have neighbor: ", neighbor_elem->id());
461 : }
462 :
463 : void
464 385423 : EFAElement2D::clearNeighbors()
465 : {
466 385423 : _general_neighbors.clear();
467 1904051 : for (unsigned int edge_iter = 0; edge_iter < _num_edges; ++edge_iter)
468 1518628 : _edge_neighbors[edge_iter] = {nullptr};
469 385423 : }
470 :
471 : void
472 385423 : EFAElement2D::setupNeighbors(std::map<EFANode *, std::set<EFAElement *>> & inverse_connectivity_map)
473 : {
474 385423 : findGeneralNeighbors(inverse_connectivity_map);
475 3135301 : for (unsigned int eit2 = 0; eit2 < _general_neighbors.size(); ++eit2)
476 : {
477 2749878 : EFAElement2D * neigh_elem = dynamic_cast<EFAElement2D *>(_general_neighbors[eit2]);
478 2749878 : if (!neigh_elem)
479 0 : EFAError("neighbor_elem is not of EFAelement2D type");
480 :
481 2749878 : std::vector<EFANode *> common_nodes = getCommonNodes(neigh_elem);
482 2749878 : if (common_nodes.size() >= 2)
483 : {
484 6836588 : for (unsigned int edge_iter = 0; edge_iter < _num_edges; ++edge_iter)
485 : {
486 5456888 : std::set<EFANode *> edge_nodes = getEdgeNodes(edge_iter);
487 : bool is_edge_neighbor = false;
488 :
489 : // Must share nodes on this edge
490 5456888 : if (Efa::numCommonElems(edge_nodes, common_nodes) == 2 && (!overlaysElement(neigh_elem)))
491 : {
492 : // Fragments must match up.
493 1375880 : if ((_fragments.size() > 1) || (neigh_elem->numFragments() > 1))
494 0 : EFAError("in updateEdgeNeighbors: Cannot have more than 1 fragment");
495 1375880 : else if ((_fragments.size() == 1) && (neigh_elem->numFragments() == 1))
496 : {
497 67698 : if (_fragments[0]->isConnected(neigh_elem->getFragment(0)))
498 : is_edge_neighbor = true;
499 : }
500 : else // If there are no fragments to match up, consider them neighbors
501 : is_edge_neighbor = true;
502 : }
503 :
504 : if (is_edge_neighbor)
505 : {
506 1375880 : if (_edge_neighbors[edge_iter][0])
507 : {
508 2324 : if (_edge_neighbors[edge_iter].size() > 1)
509 : {
510 0 : EFAError("Element ",
511 : _id,
512 : " already has 2 edge neighbors: ",
513 : _edge_neighbors[edge_iter][0]->id(),
514 : " ",
515 : _edge_neighbors[edge_iter][1]->id());
516 : }
517 2324 : _edge_neighbors[edge_iter].push_back(neigh_elem);
518 : }
519 : else
520 1373556 : _edge_neighbors[edge_iter][0] = neigh_elem;
521 : }
522 : }
523 : }
524 2749878 : }
525 385423 : }
526 :
527 : void
528 385423 : EFAElement2D::neighborSanityCheck() const
529 : {
530 1904051 : for (unsigned int edge_iter = 0; edge_iter < _num_edges; ++edge_iter)
531 : {
532 3039580 : for (unsigned int en_iter = 0; en_iter < _edge_neighbors[edge_iter].size(); ++en_iter)
533 : {
534 1520952 : EFAElement2D * neigh_elem = _edge_neighbors[edge_iter][en_iter];
535 1520952 : if (neigh_elem != nullptr)
536 : {
537 : bool found_neighbor = false;
538 6817568 : for (unsigned int edge_iter2 = 0; edge_iter2 < neigh_elem->numEdges(); ++edge_iter2)
539 : {
540 9133600 : for (unsigned int en_iter2 = 0; en_iter2 < neigh_elem->numEdgeNeighbors(edge_iter2);
541 : ++en_iter2)
542 : {
543 5067792 : if (neigh_elem->getEdgeNeighbor(edge_iter2, en_iter2) == this)
544 : {
545 1375880 : if ((en_iter2 > 1) && (en_iter > 1))
546 0 : EFAError(
547 : "Element and neighbor element cannot both have >1 neighbors on a common edge");
548 : found_neighbor = true;
549 : break;
550 : }
551 : }
552 : }
553 1375880 : if (!found_neighbor)
554 0 : EFAError("Neighbor element doesn't recognize current element as neighbor");
555 : }
556 : }
557 : }
558 385423 : }
559 :
560 : void
561 385375 : EFAElement2D::initCrackTip(std::set<EFAElement *> & CrackTipElements)
562 : {
563 385375 : if (isCrackTipElement())
564 : {
565 2048 : CrackTipElements.insert(this);
566 10200 : for (unsigned int edge_iter = 0; edge_iter < _num_edges; ++edge_iter)
567 : {
568 8152 : if ((_edge_neighbors[edge_iter].size() == 2) && (_edges[edge_iter]->hasIntersection()))
569 : {
570 : // Neither neighbor overlays current element. We are on the uncut element ahead of the tip.
571 : // Flag neighbors as crack tip split elements and add this element as their crack tip
572 : // neighbor.
573 4096 : if (_edge_neighbors[edge_iter][0]->overlaysElement(this) ||
574 2048 : _edge_neighbors[edge_iter][1]->overlaysElement(this))
575 0 : EFAError("Element has a neighbor that overlays itself");
576 :
577 : // Make sure the current elment hasn't been flagged as a tip element
578 2048 : if (_crack_tip_split_element)
579 0 : EFAError("crack_tip_split_element already flagged. In elem: ",
580 : _id,
581 : " flags: ",
582 : _crack_tip_split_element,
583 : " ",
584 : _edge_neighbors[edge_iter][0]->isCrackTipSplit(),
585 : " ",
586 : _edge_neighbors[edge_iter][1]->isCrackTipSplit());
587 :
588 2048 : _edge_neighbors[edge_iter][0]->setCrackTipSplit();
589 2048 : _edge_neighbors[edge_iter][1]->setCrackTipSplit();
590 :
591 2048 : _edge_neighbors[edge_iter][0]->addCrackTipNeighbor(this);
592 2048 : _edge_neighbors[edge_iter][1]->addCrackTipNeighbor(this);
593 : }
594 : } // edge_iter
595 : }
596 385375 : }
597 :
598 : unsigned int
599 2011 : EFAElement2D::getCrackTipSplitElementID() const
600 : {
601 2011 : if (isCrackTipElement())
602 : {
603 6092 : for (unsigned int edge_iter = 0; edge_iter < _num_edges; ++edge_iter)
604 : {
605 6092 : if ((_edge_neighbors[edge_iter].size() == 2) && (_edges[edge_iter]->hasIntersection()))
606 : {
607 4022 : if (_edge_neighbors[edge_iter][0] != nullptr &&
608 2011 : _edge_neighbors[edge_iter][0]->isCrackTipSplit())
609 : {
610 2011 : return _edge_neighbors[edge_iter][0]->id();
611 : }
612 : }
613 : }
614 : }
615 0 : EFAError("In getCrackTipSplitElementID could not find element id");
616 : return 0;
617 : }
618 :
619 : bool
620 235747 : EFAElement2D::shouldDuplicateForCrackTip(const std::set<EFAElement *> & CrackTipElements)
621 : {
622 : // This method is called in createChildElements()
623 : // Only duplicate when
624 : // 1) currElem will be a NEW crack tip element
625 : // 2) currElem is a crack tip split element at last time step and the tip will extend
626 : // 3) currElem is the neighbor of a to-be-second-split element which has another neighbor
627 : // sharing a phantom node with currElem
628 : bool should_duplicate = false;
629 235747 : if (_fragments.size() == 1)
630 : {
631 : std::set<EFAElement *>::iterator sit;
632 33436 : sit = CrackTipElements.find(this);
633 33436 : if (sit == CrackTipElements.end() && isCrackTipElement())
634 : should_duplicate = true;
635 27515 : else if (shouldDuplicateCrackTipSplitElement(CrackTipElements))
636 : should_duplicate = true;
637 20060 : else if (shouldDuplicateForPhantomCorner())
638 : should_duplicate = true;
639 : }
640 235747 : return should_duplicate;
641 : }
642 :
643 : bool
644 27515 : EFAElement2D::shouldDuplicateCrackTipSplitElement(const std::set<EFAElement *> & CrackTipElements)
645 : {
646 : // Determine whether element at crack tip should be duplicated. It should be duplicated
647 : // if the crack will extend into the next element, or if it has a non-physical node
648 : // connected to a face where a crack terminates, but will extend.
649 :
650 : bool should_duplicate = false;
651 27515 : if (_fragments.size() == 1)
652 : {
653 : std::vector<unsigned int> split_neighbors;
654 27515 : if (willCrackTipExtend(split_neighbors))
655 : should_duplicate = true;
656 : else
657 : {
658 : // The element may not be at the crack tip, but could have a non-physical node
659 : // connected to a crack tip face (on a neighbor element) that will be split. We need to
660 : // duplicate in that case as well.
661 : std::set<EFANode *> non_physical_nodes;
662 20577 : getNonPhysicalNodes(non_physical_nodes);
663 :
664 138212 : for (unsigned int eit = 0; eit < _general_neighbors.size(); ++eit)
665 : {
666 118152 : EFAElement2D * neigh_elem = dynamic_cast<EFAElement2D *>(_general_neighbors[eit]);
667 118152 : if (!neigh_elem)
668 0 : EFAError("general elem is not of type EFAelement2D");
669 :
670 : // check if a general neighbor is an old crack tip element and will be split
671 : std::set<EFAElement *>::iterator sit;
672 118152 : sit = CrackTipElements.find(neigh_elem);
673 118152 : if (sit != CrackTipElements.end() && neigh_elem->numFragments() > 1)
674 : {
675 1455 : for (unsigned int i = 0; i < neigh_elem->numEdges(); ++i)
676 : {
677 1398 : std::set<EFANode *> neigh_edge_nodes = neigh_elem->getEdgeNodes(i);
678 1972 : if (neigh_elem->numEdgeNeighbors(i) == 2 &&
679 574 : Efa::numCommonElems(neigh_edge_nodes, non_physical_nodes) > 0)
680 : {
681 : should_duplicate = true;
682 : break;
683 : }
684 : } // i
685 : }
686 : if (should_duplicate)
687 : break;
688 : } // eit
689 : }
690 27515 : } // IF one fragment
691 27515 : return should_duplicate;
692 : }
693 :
694 : bool
695 20060 : EFAElement2D::shouldDuplicateForPhantomCorner()
696 : {
697 : // if a partial element will be split for a second time and it has two neighbor elements
698 : // sharing one phantom node with the aforementioned partial element, then the two neighbor
699 : // elements should be duplicated
700 : bool should_duplicate = false;
701 20060 : if (_fragments.size() == 1 && (!_crack_tip_split_element))
702 : {
703 88098 : for (unsigned int i = 0; i < _num_edges; ++i)
704 : {
705 70132 : std::set<EFANode *> phantom_nodes = getPhantomNodeOnEdge(i);
706 70132 : if (phantom_nodes.size() > 0 && numEdgeNeighbors(i) == 1)
707 : {
708 32511 : EFAElement2D * neighbor_elem = _edge_neighbors[i][0];
709 32511 : if (neighbor_elem->numFragments() > 1) // neighbor will be split
710 : {
711 1480 : for (unsigned int j = 0; j < neighbor_elem->numEdges(); ++j)
712 : {
713 2628 : if (!neighbor_elem->getEdge(j)->equivalent(*_edges[i]) &&
714 1152 : neighbor_elem->numEdgeNeighbors(j) > 0)
715 : {
716 1148 : std::set<EFANode *> neigh_phantom_nodes = neighbor_elem->getPhantomNodeOnEdge(j);
717 1148 : if (Efa::numCommonElems(phantom_nodes, neigh_phantom_nodes) > 0)
718 : {
719 : should_duplicate = true;
720 : break;
721 : }
722 : }
723 : } // j
724 : }
725 : }
726 : if (should_duplicate)
727 : break;
728 : } // i
729 : }
730 20060 : return should_duplicate;
731 : }
732 :
733 : bool
734 27515 : EFAElement2D::willCrackTipExtend(std::vector<unsigned int> & split_neighbors) const
735 : {
736 : // Determine whether the current element is a crack tip element for which the crack will
737 : // extend into the next element.
738 : // N.B. this is called at the beginning of createChildElements
739 : bool will_extend = false;
740 27515 : if (_fragments.size() == 1 && _crack_tip_split_element)
741 : {
742 17076 : for (unsigned int i = 0; i < _crack_tip_neighbors.size(); ++i)
743 : {
744 8684 : unsigned int neigh_idx = _crack_tip_neighbors[i];
745 8684 : if (numEdgeNeighbors(neigh_idx) != 1)
746 0 : EFAError("in will_crack_tip_extend() element: ",
747 : _id,
748 : " has: ",
749 : _edge_neighbors[neigh_idx].size(),
750 : " on edge: ",
751 : neigh_idx);
752 :
753 8684 : EFAElement2D * neighbor_elem = _edge_neighbors[neigh_idx][0];
754 8684 : if (neighbor_elem->numFragments() > 2)
755 0 : EFAError("in will_crack_tip_extend() element: ",
756 : neighbor_elem->id(),
757 : " has: ",
758 : neighbor_elem->numFragments(),
759 : " fragments");
760 8684 : else if (neighbor_elem->numFragments() == 2)
761 : {
762 7082 : EFAFragment2D * frag1 = neighbor_elem->getFragment(0);
763 7082 : EFAFragment2D * frag2 = neighbor_elem->getFragment(1);
764 7082 : std::vector<EFANode *> neigh_cut_nodes = frag1->getCommonNodes(frag2);
765 7082 : if (neigh_cut_nodes.size() != 2)
766 0 : EFAError("2 frags in a elem does not share 2 common nodes");
767 10868 : if (_edges[neigh_idx]->isEmbeddedNode(neigh_cut_nodes[0]) ||
768 3786 : _edges[neigh_idx]->isEmbeddedNode(neigh_cut_nodes[1]))
769 : {
770 7082 : split_neighbors.push_back(neigh_idx);
771 : will_extend = true;
772 : }
773 7082 : }
774 : } // i
775 : }
776 27515 : return will_extend;
777 : }
778 :
779 : bool
780 429394 : EFAElement2D::isCrackTipElement() const
781 : {
782 429394 : return fragmentHasTipEdges();
783 : }
784 :
785 : unsigned int
786 36295 : EFAElement2D::getNumCuts() const
787 : {
788 : unsigned int num_cuts = 0;
789 176467 : for (unsigned int i = 0; i < _num_edges; ++i)
790 140172 : if (_edges[i]->hasIntersection())
791 0 : num_cuts += _edges[i]->numEmbeddedNodes();
792 36295 : return num_cuts;
793 : }
794 :
795 : bool
796 193983 : EFAElement2D::isFinalCut() const
797 : {
798 : // if an element has been cut twice its fragment must have two interior edges
799 : bool cut_twice = false;
800 193983 : if (_fragments.size() > 0)
801 : {
802 : unsigned int num_interior_edges = 0;
803 83012 : for (unsigned int i = 0; i < _fragments[0]->numEdges(); ++i)
804 : {
805 66261 : if (_fragments[0]->isEdgeInterior(i))
806 15800 : num_interior_edges += 1;
807 : }
808 16751 : if (num_interior_edges == 2)
809 : cut_twice = true;
810 : }
811 193983 : return cut_twice;
812 : }
813 :
814 : void
815 226846 : EFAElement2D::updateFragments(const std::set<EFAElement *> & CrackTipElements,
816 : std::map<unsigned int, EFANode *> & EmbeddedNodes)
817 : {
818 : // combine the crack-tip edges in a fragment to a single intersected edge
819 : std::set<EFAElement *>::iterator sit;
820 226846 : sit = CrackTipElements.find(this);
821 226846 : if (sit != CrackTipElements.end()) // curr_elem is a crack tip element
822 : {
823 1181 : if (_fragments.size() == 1)
824 1181 : _fragments[0]->combineTipEdges();
825 : else
826 0 : EFAError("crack tip elem ", _id, " must have 1 fragment");
827 : }
828 :
829 : // if a fragment only has 1 intersection which is in an interior edge
830 : // remove this embedded node (MUST DO THIS AFTER combine_tip_edges())
831 226846 : if (_fragments.size() == 1)
832 20453 : _fragments[0]->removeInvalidEmbeddedNodes(EmbeddedNodes);
833 :
834 : // for an element with no fragment, create one fragment identical to the element
835 226846 : if (_fragments.size() == 0)
836 206393 : _fragments.push_back(new EFAFragment2D(this, true, this));
837 226846 : if (_fragments.size() != 1)
838 0 : EFAError("Element ", _id, " must have 1 fragment at this point");
839 :
840 : // count fragment's cut edges
841 226846 : unsigned int num_cut_frag_edges = _fragments[0]->getNumCuts();
842 226846 : unsigned int num_cut_nodes = _fragments[0]->getNumCutNodes();
843 226846 : unsigned int num_frag_edges = _fragments[0]->numEdges();
844 226846 : if (num_cut_frag_edges > 3)
845 0 : EFAError("In element ", _id, " there are more than 2 cut fragment edges");
846 :
847 226846 : if (num_cut_frag_edges == 0 && num_cut_nodes == 0)
848 : {
849 221415 : if (!isPartial()) // delete the temp frag for an uncut elem
850 : {
851 202311 : delete _fragments[0];
852 202311 : _fragments.clear();
853 : }
854 : // Element has already been cut. Don't recreate fragments because we
855 : // would create multiple fragments to cover the entire element and
856 : // lose the information about what part of this element is physical.
857 221415 : return;
858 : }
859 :
860 : // split one fragment into one, two or three new fragments
861 : std::vector<EFAFragment2D *> new_frags;
862 5431 : if (num_cut_frag_edges == 3)
863 2 : new_frags = branchingSplit(EmbeddedNodes);
864 : else
865 10860 : new_frags = _fragments[0]->split();
866 :
867 5431 : delete _fragments[0]; // delete the old fragment
868 5431 : _fragments.clear();
869 14757 : for (unsigned int i = 0; i < new_frags.size(); ++i)
870 9326 : _fragments.push_back(new_frags[i]);
871 :
872 5431 : fragmentSanityCheck(num_frag_edges, num_cut_frag_edges);
873 5431 : }
874 :
875 : void
876 5431 : EFAElement2D::fragmentSanityCheck(unsigned int n_old_frag_edges, unsigned int n_old_frag_cuts) const
877 : {
878 5431 : if (n_old_frag_cuts > 3)
879 0 : EFAError("Sanity check: in element ", _id, " frag has more than 3 cut edges");
880 :
881 : // count permanent and embedded nodes for new fragments
882 : std::vector<unsigned int> num_emb;
883 : std::vector<unsigned int> num_perm;
884 : std::vector<unsigned int> num_emb_perm;
885 14757 : for (unsigned int i = 0; i < _fragments.size(); ++i)
886 : {
887 9326 : num_emb.push_back(0);
888 9326 : num_perm.push_back(0);
889 9326 : num_emb_perm.push_back(0);
890 : std::set<EFANode *> perm_nodes;
891 : std::set<EFANode *> emb_nodes;
892 : std::set<EFANode *> emb_perm_nodes;
893 47500 : for (unsigned int j = 0; j < _fragments[i]->numEdges(); ++j)
894 : {
895 114522 : for (unsigned int k = 0; k < 2; ++k)
896 : {
897 76348 : EFANode * temp_node = _fragments[i]->getEdge(j)->getNode(k);
898 76348 : if (temp_node->category() == EFANode::N_CATEGORY_PERMANENT)
899 41844 : perm_nodes.insert(temp_node);
900 34504 : else if (temp_node->category() == EFANode::N_CATEGORY_EMBEDDED)
901 34024 : emb_nodes.insert(temp_node);
902 480 : else if (temp_node->category() == EFANode::N_CATEGORY_EMBEDDED_PERMANENT)
903 480 : emb_perm_nodes.insert(temp_node);
904 : else
905 0 : EFAError("Invalid node category");
906 : }
907 : }
908 9326 : num_perm[i] = perm_nodes.size();
909 9326 : num_emb[i] = emb_nodes.size();
910 9326 : num_emb_perm[i] = emb_perm_nodes.size();
911 : }
912 :
913 : // TODO: For cut-node case, how to check fragment sanity
914 14605 : for (unsigned int i = 0; i < _fragments.size(); ++i)
915 9302 : if (num_emb_perm[i] != 0)
916 : return;
917 :
918 5303 : unsigned int n_interior_nodes = numInteriorNodes();
919 5303 : if (n_interior_nodes > 0 && n_interior_nodes != 1)
920 0 : EFAError("After update_fragments this element has ", n_interior_nodes, " interior nodes");
921 :
922 5303 : if (n_old_frag_cuts == 0)
923 : {
924 0 : if (_fragments.size() != 1 || _fragments[0]->numEdges() != n_old_frag_edges)
925 0 : EFAError("Incorrect link size for element with 0 cuts");
926 : }
927 5303 : else if (n_old_frag_cuts == 1) // crack tip case
928 : {
929 1433 : if (_fragments.size() != 1 || _fragments[0]->numEdges() != n_old_frag_edges + 1)
930 0 : EFAError("Incorrect link size for element with 1 cut");
931 : }
932 3870 : else if (n_old_frag_cuts == 2)
933 : {
934 3869 : if (_fragments.size() != 2 ||
935 3869 : (_fragments[0]->numEdges() + _fragments[1]->numEdges()) != n_old_frag_edges + 4)
936 0 : EFAError("Incorrect link size for element with 2 cuts");
937 : }
938 : else if (n_old_frag_cuts == 3)
939 : {
940 1 : if (_fragments.size() != 3 || (_fragments[0]->numEdges() + _fragments[1]->numEdges() +
941 1 : _fragments[2]->numEdges()) != n_old_frag_edges + 9)
942 0 : EFAError("Incorrect link size for element with 3 cuts");
943 : }
944 : else
945 : EFAError("Unexpected number of old fragment cuts");
946 5431 : }
947 :
948 : void
949 36295 : EFAElement2D::restoreFragment(const EFAElement * const from_elem)
950 : {
951 36295 : const EFAElement2D * from_elem2d = dynamic_cast<const EFAElement2D *>(from_elem);
952 36295 : if (!from_elem2d)
953 0 : EFAError("from_elem is not of EFAelement2D type");
954 :
955 : // restore fragments
956 36295 : if (_fragments.size() != 0)
957 0 : EFAError("in restoreFragmentInfo elements must not have any pre-existing fragments");
958 72590 : for (unsigned int i = 0; i < from_elem2d->numFragments(); ++i)
959 36295 : _fragments.push_back(new EFAFragment2D(this, true, from_elem2d, i));
960 :
961 : // restore interior nodes
962 36295 : if (_interior_nodes.size() != 0)
963 0 : EFAError("in restoreFragmentInfo elements must not have any pre-exsiting interior nodes");
964 36439 : for (unsigned int i = 0; i < from_elem2d->_interior_nodes.size(); ++i)
965 144 : _interior_nodes.push_back(new EFAFaceNode(*from_elem2d->_interior_nodes[i]));
966 :
967 : // restore edge intersections
968 36295 : if (getNumCuts() != 0)
969 0 : EFAError("In restoreEdgeIntersection: edge cuts already exist in element ", _id);
970 176467 : for (unsigned int i = 0; i < _num_edges; ++i)
971 : {
972 140172 : if (from_elem2d->_edges[i]->hasIntersection())
973 69899 : _edges[i]->copyIntersection(*from_elem2d->_edges[i], 0);
974 140172 : if (_edges[i]->numEmbeddedNodes() > 2)
975 0 : EFAError("elem ", _id, " has an edge with >2 cuts");
976 : }
977 :
978 : // replace all local nodes with global nodes
979 196503 : for (unsigned int i = 0; i < from_elem2d->numNodes(); ++i)
980 : {
981 160208 : if (from_elem2d->_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
982 160208 : switchNode(
983 : _nodes[i], from_elem2d->_nodes[i], false); // EFAelement is not a child of any parent
984 : else
985 0 : EFAError("In restoreFragmentInfo all of from_elem's nodes must be local");
986 : }
987 36295 : }
988 :
989 : void
990 226849 : EFAElement2D::createChild(const std::set<EFAElement *> & CrackTipElements,
991 : std::map<unsigned int, EFAElement *> & Elements,
992 : std::map<unsigned int, EFAElement *> & newChildElements,
993 : std::vector<EFAElement *> & ChildElements,
994 : std::vector<EFAElement *> & ParentElements,
995 : std::map<unsigned int, EFANode *> & TempNodes)
996 : {
997 226849 : if (_children.size() != 0)
998 0 : EFAError("Element cannot have existing children in createChildElements");
999 :
1000 : bool shouldDuplicateForCutNodeElement = false;
1001 1176941 : for (unsigned int j = 0; j < _num_nodes; ++j)
1002 : {
1003 950092 : if (_nodes[j]->category() == EFANode::N_CATEGORY_EMBEDDED_PERMANENT)
1004 : shouldDuplicateForCutNodeElement = true;
1005 : }
1006 :
1007 226849 : if (_fragments.size() > 1 || shouldDuplicateForCrackTip(CrackTipElements) ||
1008 : shouldDuplicateForCutNodeElement)
1009 : {
1010 5478 : if (_fragments.size() > 3)
1011 0 : EFAError("More than 3 fragments not yet supported");
1012 :
1013 : // set up the children
1014 5478 : ParentElements.push_back(this);
1015 14851 : for (unsigned int ichild = 0; ichild < _fragments.size(); ++ichild)
1016 : {
1017 : unsigned int new_elem_id;
1018 9373 : if (newChildElements.size() == 0)
1019 871 : new_elem_id = Efa::getNewID(Elements);
1020 : else
1021 8502 : new_elem_id = Efa::getNewID(newChildElements);
1022 :
1023 9373 : EFAElement2D * childElem = new EFAElement2D(new_elem_id, this->numNodes());
1024 9373 : newChildElements.insert(std::make_pair(new_elem_id, childElem));
1025 :
1026 9373 : ChildElements.push_back(childElem);
1027 9373 : childElem->setParent(this);
1028 9373 : _children.push_back(childElem);
1029 :
1030 : std::vector<EFAPoint> local_embedded_node_coor;
1031 :
1032 46854 : for (unsigned int i = 0; i < this->getFragment(ichild)->numEdges(); ++i)
1033 : {
1034 37481 : if (this->getFragment(ichild)->isEdgeInterior(i))
1035 : {
1036 : std::vector<EFANode *> master_nodes;
1037 : std::vector<double> master_weights;
1038 :
1039 26118 : for (unsigned int j = 0; j < 2; ++j)
1040 : {
1041 17412 : this->getMasterInfo(
1042 : this->getFragmentEdge(ichild, i)->getNode(j), master_nodes, master_weights);
1043 17412 : EFAPoint coor(0.0, 0.0, 0.0);
1044 52336 : for (unsigned int k = 0; k < master_nodes.size(); ++k)
1045 : {
1046 34924 : EFANode * local = this->createLocalNodeFromGlobalNode(master_nodes[k]);
1047 34924 : coor += _local_node_coor[local->id()] * master_weights[k];
1048 34924 : delete local;
1049 : }
1050 17412 : local_embedded_node_coor.push_back(coor);
1051 : }
1052 8706 : }
1053 : }
1054 :
1055 9373 : EFAPoint normal(0.0, 0.0, 0.0);
1056 9373 : EFAPoint origin(0.0, 0.0, 0.0);
1057 9373 : EFAPoint normal2(0.0, 0.0, 0.0);
1058 9373 : EFAPoint origin2(0.0, 0.0, 0.0);
1059 :
1060 9373 : if (local_embedded_node_coor.size())
1061 : {
1062 8637 : EFAPoint cut_line = local_embedded_node_coor[1] - local_embedded_node_coor[0];
1063 8637 : normal = EFAPoint(cut_line(1), -cut_line(0), 0.0);
1064 8637 : Xfem::normalizePoint(normal);
1065 8637 : origin = (local_embedded_node_coor[0] + local_embedded_node_coor[1]) * 0.5;
1066 : }
1067 :
1068 9373 : if (local_embedded_node_coor.size() == 4)
1069 : {
1070 69 : EFAPoint cut_line = local_embedded_node_coor[3] - local_embedded_node_coor[2];
1071 69 : normal2 = EFAPoint(cut_line(1), -cut_line(0), 0.0);
1072 69 : Xfem::normalizePoint(normal2);
1073 69 : origin2 = (local_embedded_node_coor[2] + local_embedded_node_coor[3]) * 0.5;
1074 : }
1075 :
1076 : // get child element's nodes
1077 50509 : for (unsigned int j = 0; j < _num_nodes; ++j)
1078 : {
1079 41136 : EFAPoint p(0.0, 0.0, 0.0);
1080 41136 : p = _local_node_coor[j];
1081 :
1082 41136 : EFAPoint origin_to_point = p - origin;
1083 41136 : EFAPoint origin2_to_point = p - origin2;
1084 :
1085 47820 : if (_fragments.size() == 1 &&
1086 6684 : _nodes[j]->category() == EFANode::N_CATEGORY_EMBEDDED_PERMANENT) // create temp node for
1087 : // embedded permanent
1088 : // node
1089 : {
1090 160 : unsigned int new_node_id = Efa::getNewID(TempNodes);
1091 160 : EFANode * newNode = new EFANode(new_node_id, EFANode::N_CATEGORY_TEMP, _nodes[j]);
1092 160 : TempNodes.insert(std::make_pair(new_node_id, newNode));
1093 160 : childElem->setNode(j, newNode); // be a temp node
1094 : }
1095 40976 : else if (_fragments.size() == 1 && !shouldDuplicateForCrackTip(CrackTipElements))
1096 : {
1097 256 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1098 : }
1099 40720 : else if (std::abs(origin_to_point * normal) < Xfem::tol &&
1100 : _fragments.size() > 1) // cut through node case
1101 : {
1102 80 : unsigned int new_node_id = Efa::getNewID(TempNodes);
1103 80 : EFANode * newNode = new EFANode(new_node_id, EFANode::N_CATEGORY_TEMP, _nodes[j]);
1104 80 : TempNodes.insert(std::make_pair(new_node_id, newNode));
1105 80 : childElem->setNode(j, newNode); // be a temp node
1106 : }
1107 40640 : else if (origin_to_point * normal < Xfem::tol && origin2_to_point * normal2 < Xfem::tol &&
1108 4321 : (_fragments.size() > 1 || shouldDuplicateForCrackTip(CrackTipElements)))
1109 : {
1110 21433 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1111 : }
1112 19207 : else if (normal.norm() < Xfem::tol && normal2.norm() < Xfem::tol &&
1113 : _fragments.size() == 1) // cut along edge case
1114 : {
1115 0 : childElem->setNode(j, _nodes[j]); // inherit parent's node
1116 : }
1117 21154 : else if ((_fragments.size() > 1 ||
1118 1947 : shouldDuplicateForCrackTip(
1119 : CrackTipElements))) // parent element's node is not in fragment
1120 : {
1121 19207 : unsigned int new_node_id = Efa::getNewID(TempNodes);
1122 19207 : EFANode * newNode = new EFANode(new_node_id, EFANode::N_CATEGORY_TEMP, _nodes[j]);
1123 19207 : TempNodes.insert(std::make_pair(new_node_id, newNode));
1124 19207 : childElem->setNode(j, newNode); // be a temp node
1125 : }
1126 : }
1127 :
1128 : // get child element's fragments
1129 9373 : EFAFragment2D * new_frag = new EFAFragment2D(childElem, true, this, ichild);
1130 9373 : childElem->_fragments.push_back(new_frag);
1131 :
1132 : // get child element's edges
1133 45745 : for (unsigned int j = 0; j < _num_edges; ++j)
1134 : {
1135 36372 : unsigned int jplus1(j < (_num_edges - 1) ? j + 1 : 0);
1136 36372 : EFAEdge * new_edge = new EFAEdge(childElem->getNode(j), childElem->getNode(jplus1));
1137 36372 : if (_edges[j]->hasIntersection())
1138 17882 : new_edge->copyIntersection(*_edges[j], 0);
1139 36372 : if ((_num_edges == 4 && _num_nodes > 4) || (_num_edges == 3 && _num_nodes > 3))
1140 4136 : new_edge->setInteriorNode(childElem->getNode(_num_edges + j));
1141 36372 : childElem->setEdge(j, new_edge);
1142 : }
1143 9373 : childElem->removePhantomEmbeddedNode(); // IMPORTANT
1144 :
1145 : // inherit old interior nodes
1146 9430 : for (unsigned int j = 0; j < _interior_nodes.size(); ++j)
1147 57 : childElem->_interior_nodes.push_back(new EFAFaceNode(*_interior_nodes[j]));
1148 9373 : }
1149 : }
1150 : else // num_links == 1 || num_links == 0
1151 : // child is itself - but don't insert into the list of ChildElements!!!
1152 221371 : _children.push_back(this);
1153 226849 : }
1154 :
1155 : void
1156 9373 : EFAElement2D::removePhantomEmbeddedNode()
1157 : {
1158 : // remove the embedded nodes on edge that are not inside the real part
1159 9373 : if (_fragments.size() > 0)
1160 : {
1161 45745 : for (unsigned int i = 0; i < _num_edges; ++i)
1162 : {
1163 : std::vector<EFANode *> nodes_to_delete;
1164 54302 : for (unsigned int j = 0; j < _edges[i]->numEmbeddedNodes(); ++j)
1165 : {
1166 17930 : if (!_fragments[0]->containsNode(_edges[i]->getEmbeddedNode(j)))
1167 79 : nodes_to_delete.push_back(_edges[i]->getEmbeddedNode(j));
1168 : }
1169 36451 : for (unsigned int j = 0; j < nodes_to_delete.size(); ++j)
1170 79 : _edges[i]->removeEmbeddedNode(nodes_to_delete[j]);
1171 36372 : } // i
1172 : }
1173 9373 : }
1174 :
1175 : void
1176 9373 : EFAElement2D::connectNeighbors(std::map<unsigned int, EFANode *> & PermanentNodes,
1177 : std::map<unsigned int, EFANode *> & TempNodes,
1178 : std::map<EFANode *, std::set<EFAElement *>> & InverseConnectivityMap,
1179 : bool merge_phantom_edges)
1180 : {
1181 : // N.B. "this" must point to a child element that was just created
1182 9373 : if (!_parent)
1183 0 : EFAError("no parent element for child element ", _id, " in connect_neighbors");
1184 9373 : EFAElement2D * parent2d = dynamic_cast<EFAElement2D *>(_parent);
1185 9373 : if (!parent2d)
1186 0 : EFAError("cannot dynamic cast to parent2d in connect_neighbors");
1187 :
1188 : // First loop through edges and merge nodes with neighbors as appropriate
1189 45745 : for (unsigned int j = 0; j < _num_edges; ++j)
1190 : {
1191 69596 : for (unsigned int k = 0; k < parent2d->numEdgeNeighbors(j); ++k)
1192 : {
1193 33224 : EFAElement2D * NeighborElem = parent2d->getEdgeNeighbor(j, k);
1194 33224 : unsigned int neighbor_edge_id = NeighborElem->getNeighborIndex(parent2d);
1195 :
1196 33224 : if (_edges[j]->hasIntersection())
1197 : {
1198 45877 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1199 : {
1200 : EFAElement2D * childOfNeighborElem =
1201 29271 : dynamic_cast<EFAElement2D *>(NeighborElem->getChild(l));
1202 29271 : if (!childOfNeighborElem)
1203 0 : EFAError("dynamic cast childOfNeighborElem fails");
1204 :
1205 : // Check to see if the nodes are already merged. There's nothing else to do in that case.
1206 29271 : EFAEdge * neighborChildEdge = childOfNeighborElem->getEdge(neighbor_edge_id);
1207 29271 : if (_edges[j]->equivalent(*neighborChildEdge))
1208 7828 : continue;
1209 :
1210 21443 : if (_fragments[0]->isConnected(childOfNeighborElem->getFragment(0)))
1211 : {
1212 : unsigned int num_edge_nodes = 2;
1213 25575 : for (unsigned int i = 0; i < num_edge_nodes; ++i)
1214 : {
1215 : unsigned int childNodeIndex = i;
1216 17050 : unsigned int neighborChildNodeIndex = num_edge_nodes - 1 - childNodeIndex;
1217 :
1218 17050 : EFANode * childNode = _edges[j]->getNode(childNodeIndex);
1219 17050 : EFANode * childOfNeighborNode = neighborChildEdge->getNode(neighborChildNodeIndex);
1220 :
1221 17050 : mergeNodes(
1222 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1223 : }
1224 :
1225 8525 : EFANode * childNode = _edges[j]->getInteriorNode();
1226 8525 : EFANode * childOfNeighborNode = neighborChildEdge->getInteriorNode();
1227 :
1228 8525 : mergeNodes(
1229 : childNode, childOfNeighborNode, childOfNeighborElem, PermanentNodes, TempNodes);
1230 : }
1231 : } // l, loop over NeighborElem's children
1232 : }
1233 : else // No edge intersection -- optionally merge non-material nodes if they share a common
1234 : // parent
1235 : {
1236 16618 : if (merge_phantom_edges)
1237 : {
1238 33229 : for (unsigned int l = 0; l < NeighborElem->numChildren(); ++l)
1239 : {
1240 : EFAElement2D * childOfNeighborElem =
1241 16618 : dynamic_cast<EFAElement2D *>(NeighborElem->getChild(l));
1242 16618 : if (!childOfNeighborElem)
1243 0 : EFAError("dynamic cast childOfNeighborElem fails");
1244 :
1245 16618 : EFAEdge * neighborChildEdge = childOfNeighborElem->getEdge(neighbor_edge_id);
1246 16618 : if (!neighborChildEdge->hasIntersection()) // neighbor edge must NOT have
1247 : // intersection either
1248 : {
1249 : // Check to see if the nodes are already merged. There's nothing else to do in that
1250 : // case.
1251 : unsigned int num_edge_nodes = 2;
1252 49680 : for (unsigned int i = 0; i < num_edge_nodes; ++i)
1253 : {
1254 : unsigned int childNodeIndex = i;
1255 33120 : unsigned int neighborChildNodeIndex = num_edge_nodes - 1 - childNodeIndex;
1256 :
1257 33120 : EFANode * childNode = _edges[j]->getNode(childNodeIndex);
1258 33120 : EFANode * childOfNeighborNode = neighborChildEdge->getNode(neighborChildNodeIndex);
1259 :
1260 45921 : if (childNode->parent() != nullptr &&
1261 12801 : childNode->parent() ==
1262 : childOfNeighborNode
1263 12801 : ->parent()) // non-material node and both come from same parent
1264 4 : mergeNodes(childNode,
1265 : childOfNeighborNode,
1266 : childOfNeighborElem,
1267 : PermanentNodes,
1268 : TempNodes);
1269 : } // i
1270 : }
1271 : } // loop over NeighborElem's children
1272 : } // if (merge_phantom_edges)
1273 : } // IF edge-j has_intersection()
1274 : } // k, loop over neighbors on edge j
1275 : } // j, loop over all edges
1276 :
1277 : // Now do a second loop through edges and convert remaining nodes to permanent nodes.
1278 : // Important: temp nodes are not shared by any neighbor, so no need to switch nodes on neighbors
1279 50509 : for (unsigned int j = 0; j < _num_nodes; ++j)
1280 : {
1281 41136 : EFANode * childNode = _nodes[j];
1282 41136 : if (childNode->category() == EFANode::N_CATEGORY_TEMP)
1283 : {
1284 : // if current child element does not have siblings, and if current temp node is a lone one
1285 : // this temp node should be merged back to its parent permanent node. Otherwise we would have
1286 : // permanent nodes that are not connected to any element
1287 4204 : std::set<EFAElement *> patch_elems = InverseConnectivityMap[childNode->parent()];
1288 4204 : if (parent2d->numFragments() == 1 && patch_elems.size() == 1)
1289 301 : switchNode(childNode->parent(), childNode, false);
1290 : else
1291 : {
1292 3903 : unsigned int new_node_id = Efa::getNewID(PermanentNodes);
1293 : EFANode * newNode =
1294 3903 : new EFANode(new_node_id, EFANode::N_CATEGORY_PERMANENT, childNode->parent());
1295 3903 : PermanentNodes.insert(std::make_pair(new_node_id, newNode));
1296 3903 : switchNode(newNode, childNode, false);
1297 : }
1298 4204 : if (!Efa::deleteFromMap(TempNodes, childNode))
1299 0 : EFAError(
1300 : "Attempted to delete node: ", childNode->id(), " from TempNodes, but couldn't find it");
1301 : }
1302 : }
1303 9373 : }
1304 :
1305 : void
1306 9373 : EFAElement2D::updateFragmentNode()
1307 : {
1308 50509 : for (unsigned int j = 0; j < _num_nodes; ++j)
1309 : {
1310 58188 : if (_nodes[j]->parent() != nullptr &&
1311 17052 : _nodes[j]->parent()->category() == EFANode::N_CATEGORY_EMBEDDED_PERMANENT)
1312 240 : switchNode(_nodes[j], _nodes[j]->parent(), false);
1313 : }
1314 9373 : }
1315 :
1316 : void
1317 603 : EFAElement2D::printElement(std::ostream & ostream) const
1318 : {
1319 : ostream << std::setw(4);
1320 603 : ostream << _id << " | ";
1321 3015 : for (unsigned int j = 0; j < _num_nodes; ++j)
1322 : {
1323 4824 : ostream << std::setw(5) << _nodes[j]->idCatString();
1324 : }
1325 :
1326 603 : ostream << " | ";
1327 3015 : for (unsigned int j = 0; j < _num_edges; ++j)
1328 : {
1329 : ostream << std::setw(4);
1330 2412 : if (_edges[j]->hasIntersection())
1331 : {
1332 579 : if (_edges[j]->numEmbeddedNodes() > 1)
1333 : {
1334 0 : ostream << "[";
1335 0 : for (unsigned int k = 0; k < _edges[j]->numEmbeddedNodes(); ++k)
1336 : {
1337 0 : ostream << _edges[j]->getEmbeddedNode(k)->id();
1338 0 : if (k == _edges[j]->numEmbeddedNodes() - 1)
1339 0 : ostream << "]";
1340 : else
1341 0 : ostream << " ";
1342 : }
1343 : }
1344 : else
1345 1158 : ostream << _edges[j]->getEmbeddedNode(0)->id() << " ";
1346 : }
1347 : else
1348 1833 : ostream << " -- ";
1349 : }
1350 603 : ostream << " | ";
1351 3015 : for (unsigned int j = 0; j < _num_edges; ++j)
1352 : {
1353 2412 : if (numEdgeNeighbors(j) > 1)
1354 : {
1355 12 : ostream << "[";
1356 36 : for (unsigned int k = 0; k < numEdgeNeighbors(j); ++k)
1357 : {
1358 24 : ostream << getEdgeNeighbor(j, k)->id();
1359 24 : if (k == numEdgeNeighbors(j) - 1)
1360 12 : ostream << "]";
1361 : else
1362 12 : ostream << " ";
1363 : }
1364 12 : ostream << " ";
1365 : }
1366 : else
1367 : {
1368 : ostream << std::setw(4);
1369 2400 : if (numEdgeNeighbors(j) == 1)
1370 2904 : ostream << getEdgeNeighbor(j, 0)->id() << " ";
1371 : else
1372 948 : ostream << " -- ";
1373 : }
1374 : }
1375 603 : ostream << " | ";
1376 1020 : for (unsigned int j = 0; j < _fragments.size(); ++j)
1377 : {
1378 : ostream << std::setw(4);
1379 834 : ostream << " " << j << " | ";
1380 2112 : for (unsigned int k = 0; k < _fragments[j]->numEdges(); ++k)
1381 : {
1382 1695 : EFANode * prt_node = getFragmentEdge(j, k)->getNode(0);
1383 1695 : unsigned int kprev(k > 0 ? k - 1 : _fragments[j]->numEdges() - 1);
1384 1695 : if (!getFragmentEdge(j, kprev)->containsNode(prt_node))
1385 0 : prt_node = getFragmentEdge(j, k)->getNode(1);
1386 3390 : ostream << std::setw(5) << prt_node->idCatString();
1387 : }
1388 : }
1389 : ostream << std::endl;
1390 603 : }
1391 :
1392 : EFAFragment2D *
1393 22374899 : EFAElement2D::getFragment(unsigned int frag_id) const
1394 : {
1395 22374899 : if (frag_id < _fragments.size())
1396 22374899 : return _fragments[frag_id];
1397 : else
1398 0 : EFAError("frag_id out of bounds");
1399 : }
1400 :
1401 : std::set<EFANode *>
1402 5458286 : EFAElement2D::getEdgeNodes(unsigned int edge_id) const
1403 : {
1404 : std::set<EFANode *> edge_nodes;
1405 5458286 : edge_nodes.insert(_edges[edge_id]->getNode(0));
1406 5458286 : edge_nodes.insert(_edges[edge_id]->getNode(1));
1407 5458286 : return edge_nodes;
1408 : }
1409 :
1410 : bool
1411 199 : EFAElement2D::getEdgeNodeParametricCoordinate(EFANode * node, std::vector<double> & para_coor) const
1412 : {
1413 : // get the parametric coords of a node in an element edge
1414 : unsigned int edge_id = std::numeric_limits<unsigned int>::max();
1415 : bool edge_found = false;
1416 446 : for (unsigned int i = 0; i < _num_edges; ++i)
1417 : {
1418 446 : if (_edges[i]->containsNode(node))
1419 : {
1420 : edge_id = i;
1421 : edge_found = true;
1422 : break;
1423 : }
1424 : }
1425 199 : if (edge_found)
1426 : {
1427 199 : double rel_dist = _edges[edge_id]->distanceFromNode1(node);
1428 199 : double xi_1d = 2.0 * rel_dist - 1.0; // translate to [-1,1] parent coord syst
1429 199 : mapParametricCoordFrom1Dto2D(edge_id, xi_1d, para_coor);
1430 : }
1431 199 : return edge_found;
1432 : }
1433 :
1434 : EFAFaceNode *
1435 0 : EFAElement2D::getInteriorNode(unsigned int interior_node_id) const
1436 : {
1437 0 : if (interior_node_id < _interior_nodes.size())
1438 0 : return _interior_nodes[interior_node_id];
1439 : else
1440 0 : EFAError("interior_node_id out of bounds");
1441 : }
1442 :
1443 : void
1444 72 : EFAElement2D::deleteInteriorNodes()
1445 : {
1446 144 : for (unsigned int i = 0; i < _interior_nodes.size(); ++i)
1447 72 : delete _interior_nodes[i];
1448 72 : _interior_nodes.clear();
1449 72 : }
1450 :
1451 : unsigned int
1452 39379974 : EFAElement2D::numEdges() const
1453 : {
1454 39379974 : return _num_edges;
1455 : }
1456 :
1457 : void
1458 36372 : EFAElement2D::setEdge(unsigned int edge_id, EFAEdge * edge)
1459 : {
1460 36372 : _edges[edge_id] = edge;
1461 36372 : }
1462 :
1463 : void
1464 385334 : EFAElement2D::createEdges()
1465 : {
1466 1903606 : for (unsigned int i = 0; i < _num_edges; ++i)
1467 : {
1468 1518272 : unsigned int i_plus1(i < (_num_edges - 1) ? i + 1 : 0);
1469 1518272 : EFAEdge * new_edge = new EFAEdge(_nodes[i], _nodes[i_plus1]);
1470 :
1471 1518272 : if ((_num_edges == 4 && _num_nodes > 4) || (_num_edges == 3 && _num_nodes > 3))
1472 79840 : new_edge->setInteriorNode(
1473 79840 : _nodes[i + _num_edges]); // '_num_edges' is the offset of interior edge node numbering
1474 :
1475 1518272 : _edges[i] = new_edge;
1476 : }
1477 385334 : }
1478 :
1479 : EFAEdge *
1480 28633850 : EFAElement2D::getEdge(unsigned int edge_id) const
1481 : {
1482 28633850 : return _edges[edge_id];
1483 : }
1484 :
1485 : EFAEdge *
1486 5603222 : EFAElement2D::getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
1487 : {
1488 5603222 : if (frag_id < _fragments.size())
1489 5603222 : return _fragments[frag_id]->getEdge(edge_id);
1490 : else
1491 0 : EFAError("frag_id out of bounds in get_frag_edge()");
1492 : }
1493 :
1494 : std::set<EFANode *>
1495 71280 : EFAElement2D::getPhantomNodeOnEdge(unsigned int edge_id) const
1496 : {
1497 : std::set<EFANode *> phantom_nodes;
1498 71280 : if (_fragments.size() > 0)
1499 : {
1500 213840 : for (unsigned int j = 0; j < 2; ++j) // loop ove 2 edge nodes
1501 : {
1502 : bool node_in_frag = false;
1503 211040 : for (unsigned int k = 0; k < _fragments.size(); ++k)
1504 : {
1505 144350 : if (_fragments[k]->containsNode(_edges[edge_id]->getNode(j)))
1506 : {
1507 : node_in_frag = true;
1508 : break;
1509 : }
1510 : }
1511 142560 : if (!node_in_frag)
1512 66690 : phantom_nodes.insert(_edges[edge_id]->getNode(j));
1513 : }
1514 : }
1515 71280 : return phantom_nodes;
1516 : }
1517 :
1518 : bool
1519 8789 : EFAElement2D::getFragmentEdgeID(unsigned int elem_edge_id, unsigned int & frag_edge_id) const
1520 : {
1521 : // find the fragment edge that is contained by given element edge
1522 : // N.B. if the elem edge contains two frag edges, this method will only return
1523 : // the first frag edge ID
1524 : bool frag_edge_found = false;
1525 8789 : frag_edge_id = std::numeric_limits<unsigned int>::max();
1526 8789 : if (_fragments.size() == 1)
1527 : {
1528 3176 : for (unsigned int j = 0; j < _fragments[0]->numEdges(); ++j)
1529 : {
1530 3176 : if (_edges[elem_edge_id]->containsEdge(*_fragments[0]->getEdge(j)))
1531 : {
1532 1244 : frag_edge_id = j;
1533 : frag_edge_found = true;
1534 1244 : break;
1535 : }
1536 : }
1537 : }
1538 8789 : return frag_edge_found;
1539 : }
1540 :
1541 : bool
1542 47247 : EFAElement2D::isEdgePhantom(unsigned int edge_id) const
1543 : {
1544 : bool is_phantom = false;
1545 47247 : if (_fragments.size() > 0)
1546 : {
1547 : bool contain_frag_edge = false;
1548 39801 : for (unsigned int i = 0; i < _fragments.size(); ++i)
1549 : {
1550 93954 : for (unsigned int j = 0; j < _fragments[i]->numEdges(); ++j)
1551 : {
1552 93858 : if (_edges[edge_id]->containsEdge(*_fragments[i]->getEdge(j)))
1553 : {
1554 : contain_frag_edge = true;
1555 : break;
1556 : }
1557 : } // j
1558 39705 : if (contain_frag_edge)
1559 : break;
1560 : } // i
1561 39705 : if (!contain_frag_edge)
1562 : is_phantom = true;
1563 : }
1564 47247 : return is_phantom;
1565 : }
1566 :
1567 : unsigned int
1568 9365040 : EFAElement2D::numEdgeNeighbors(unsigned int edge_id) const
1569 : {
1570 : unsigned int num_neighbors = 0;
1571 9365040 : if (_edge_neighbors[edge_id][0])
1572 8955027 : num_neighbors = _edge_neighbors[edge_id].size();
1573 9365040 : return num_neighbors;
1574 : }
1575 :
1576 : EFAElement2D *
1577 5149099 : EFAElement2D::getEdgeNeighbor(unsigned int edge_id, unsigned int neighbor_id) const
1578 : {
1579 5149099 : if (_edge_neighbors[edge_id][0] != nullptr && neighbor_id < _edge_neighbors[edge_id].size())
1580 5149099 : return _edge_neighbors[edge_id][neighbor_id];
1581 : else
1582 0 : EFAError("edge neighbor does not exist");
1583 : }
1584 :
1585 : bool
1586 429394 : EFAElement2D::fragmentHasTipEdges() const
1587 : {
1588 : bool has_tip_edges = false;
1589 429394 : if (_fragments.size() == 1)
1590 : {
1591 370661 : for (unsigned int i = 0; i < _num_edges; ++i)
1592 : {
1593 : unsigned int num_frag_edges = 0; // count how many fragment edges this element edge contains
1594 300888 : if (_edges[i]->hasIntersection())
1595 : {
1596 743004 : for (unsigned int j = 0; j < _fragments[0]->numEdges(); ++j)
1597 : {
1598 594166 : if (_edges[i]->containsEdge(*_fragments[0]->getEdge(j)))
1599 159371 : num_frag_edges += 1;
1600 : } // j
1601 148838 : if (num_frag_edges == 2)
1602 : {
1603 : has_tip_edges = true;
1604 : break;
1605 : }
1606 : }
1607 : } // i
1608 : }
1609 429394 : return has_tip_edges;
1610 : }
1611 :
1612 : unsigned int
1613 9 : EFAElement2D::getTipEdgeID() const
1614 : {
1615 : // if this element is a crack tip element, returns the crack tip edge's ID
1616 : unsigned int tip_edge_id = std::numeric_limits<unsigned int>::max();
1617 9 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1618 : {
1619 18 : for (unsigned int i = 0; i < _num_edges; ++i)
1620 : {
1621 : unsigned int num_frag_edges = 0; // count how many fragment edges this element edge contains
1622 18 : if (_edges[i]->hasIntersection())
1623 : {
1624 54 : for (unsigned int j = 0; j < _fragments[0]->numEdges(); ++j)
1625 : {
1626 45 : if (_edges[i]->containsEdge(*_fragments[0]->getEdge(j)))
1627 18 : num_frag_edges += 1;
1628 : }
1629 9 : if (num_frag_edges == 2) // element edge contains two fragment edges
1630 : {
1631 : tip_edge_id = i;
1632 : break;
1633 : }
1634 : }
1635 : }
1636 : }
1637 9 : return tip_edge_id;
1638 : }
1639 :
1640 : EFANode *
1641 2020 : EFAElement2D::getTipEmbeddedNode() const
1642 : {
1643 : // if this element is a crack tip element, returns the crack tip edge's ID
1644 : EFANode * tip_emb = nullptr;
1645 2020 : if (_fragments.size() == 1) // crack tip element with a partial fragment saved
1646 : {
1647 6110 : for (unsigned int i = 0; i < _num_edges; ++i)
1648 : {
1649 : std::vector<EFAEdge *> frag_edges; // count how many fragment edges this element edge contains
1650 6110 : if (_edges[i]->hasIntersection())
1651 : {
1652 12080 : for (unsigned int j = 0; j < _fragments[0]->numEdges(); ++j)
1653 : {
1654 10060 : if (_edges[i]->containsEdge(*_fragments[0]->getEdge(j)))
1655 4040 : frag_edges.push_back(_fragments[0]->getEdge(j));
1656 : } // j
1657 2020 : if (frag_edges.size() == 2) // element edge contains two fragment edges
1658 : {
1659 2020 : if (frag_edges[1]->containsNode(frag_edges[0]->getNode(1)))
1660 2020 : tip_emb = frag_edges[0]->getNode(1);
1661 0 : else if (frag_edges[1]->containsNode(frag_edges[0]->getNode(0)))
1662 0 : tip_emb = frag_edges[0]->getNode(0);
1663 : else
1664 0 : EFAError("Common node can't be found between 2 tip frag edges");
1665 : break;
1666 : }
1667 : }
1668 6110 : }
1669 : }
1670 2020 : return tip_emb;
1671 : }
1672 :
1673 : bool
1674 622 : EFAElement2D::edgeContainsTip(unsigned int edge_id) const
1675 : {
1676 : bool contains_tip = false;
1677 622 : if (_fragments.size() == 1)
1678 : {
1679 : unsigned int num_frag_edges = 0; // count how many fragment edges this element edge contains
1680 622 : if (_edges[edge_id]->hasIntersection())
1681 : {
1682 936 : for (unsigned int j = 0; j < _fragments[0]->numEdges(); ++j)
1683 : {
1684 720 : if (_edges[edge_id]->containsEdge(*_fragments[0]->getEdge(j)))
1685 216 : num_frag_edges += 1;
1686 : } // j
1687 216 : if (num_frag_edges == 2)
1688 : contains_tip = true;
1689 : }
1690 : }
1691 622 : return contains_tip;
1692 : }
1693 :
1694 : bool
1695 622 : EFAElement2D::fragmentEdgeAlreadyCut(unsigned int ElemEdgeID) const
1696 : {
1697 : // when marking cuts, check if the corresponding frag edge already has been cut
1698 : bool has_cut = false;
1699 622 : if (edgeContainsTip(ElemEdgeID))
1700 : has_cut = true;
1701 : else
1702 : {
1703 622 : unsigned int FragEdgeID = std::numeric_limits<unsigned int>::max();
1704 622 : if (getFragmentEdgeID(ElemEdgeID, FragEdgeID))
1705 : {
1706 622 : EFAEdge * frag_edge = getFragmentEdge(0, FragEdgeID);
1707 622 : if (frag_edge->hasIntersection())
1708 : has_cut = true;
1709 : }
1710 : }
1711 622 : return has_cut;
1712 : }
1713 :
1714 : void
1715 92206 : EFAElement2D::addEdgeCut(unsigned int edge_id,
1716 : double position,
1717 : EFANode * embedded_node,
1718 : std::map<unsigned int, EFANode *> & EmbeddedNodes,
1719 : bool add_to_neighbor)
1720 : {
1721 : EFANode * local_embedded = nullptr;
1722 92206 : EFANode * edge_node1 = _edges[edge_id]->getNode(0);
1723 : if (embedded_node) // use the existing embedded node if it was passed in
1724 : local_embedded = embedded_node;
1725 :
1726 92206 : if (_edges[edge_id]->hasIntersectionAtPosition(position, edge_node1) && position > Xfem::tol &&
1727 : position < 1.0 - Xfem::tol)
1728 : {
1729 84039 : unsigned int emb_id = _edges[edge_id]->getEmbeddedNodeIndex(position, edge_node1);
1730 84039 : EFANode * old_emb = _edges[edge_id]->getEmbeddedNode(emb_id);
1731 84039 : if (embedded_node && embedded_node != old_emb)
1732 : {
1733 0 : EFAError("Attempting to add edge intersection when one already exists with different node.",
1734 : " elem: ",
1735 : _id,
1736 : " edge: ",
1737 : edge_id,
1738 : " position: ",
1739 : position);
1740 : }
1741 : local_embedded = old_emb;
1742 : }
1743 : else // if no cut exists at the input position
1744 : {
1745 : bool add2elem = true;
1746 :
1747 : // check if it is necessary to add cuts to fragment
1748 : // id of partially overlapping fragment edge
1749 8167 : unsigned int frag_edge_id = std::numeric_limits<unsigned int>::max();
1750 : EFAEdge * frag_edge = nullptr;
1751 : EFANode * frag_edge_node1 = nullptr;
1752 : double frag_pos = -1.0;
1753 : bool add2frag = false;
1754 :
1755 8167 : if (getFragmentEdgeID(edge_id, frag_edge_id)) // elem edge contains a frag edge
1756 : {
1757 622 : frag_edge = getFragmentEdge(0, frag_edge_id);
1758 622 : if (!fragmentEdgeAlreadyCut(edge_id))
1759 : {
1760 : double xi[2] = {-1.0, -1.0}; // relative coords of two frag edge nodes
1761 622 : xi[0] = _edges[edge_id]->distanceFromNode1(frag_edge->getNode(0));
1762 622 : xi[1] = _edges[edge_id]->distanceFromNode1(frag_edge->getNode(1));
1763 622 : if ((position - xi[0]) * (position - xi[1]) <
1764 : 0.0) // the cut to be added is within the real part of the edge
1765 : {
1766 430 : frag_edge_node1 = frag_edge->getNode(0);
1767 430 : frag_pos = (position - xi[0]) / (xi[1] - xi[0]);
1768 : add2frag = true;
1769 : }
1770 : else // the emb node to be added is in the phantom part of the elem edge
1771 : add2elem = false; // DO NOT ADD INTERSECT IN THIS CASE
1772 : }
1773 : else
1774 : {
1775 : EFAWarning("attempting to add new cut to a cut fragment edge");
1776 : add2elem = false; // DO NOT ADD INTERSECT IN THIS CASE
1777 : }
1778 : }
1779 :
1780 : // If elem edge has 2 cuts but they have not been restored yet, it's OK because
1781 : // getFragmentEdgeID = false so we won't add anything to the restored fragment.
1782 : // add to elem edge (IMPORTANT to do it AFTER the above fragment check)
1783 8167 : if (add2elem)
1784 : {
1785 7975 : if (!local_embedded) // need to create new embedded node
1786 : {
1787 4555 : unsigned int new_node_id = Efa::getNewID(EmbeddedNodes);
1788 4555 : local_embedded = new EFANode(new_node_id, EFANode::N_CATEGORY_EMBEDDED);
1789 4555 : EmbeddedNodes.insert(std::make_pair(new_node_id, local_embedded));
1790 : }
1791 7975 : _edges[edge_id]->addIntersection(position, local_embedded, edge_node1);
1792 7975 : if (_edges[edge_id]->numEmbeddedNodes() > 2)
1793 0 : EFAError("element edge can't have >2 embedded nodes");
1794 : }
1795 :
1796 : // add to frag edge
1797 8167 : if (add2frag)
1798 : {
1799 430 : frag_edge->addIntersection(frag_pos, local_embedded, frag_edge_node1);
1800 430 : if (frag_edge->numEmbeddedNodes() > 1)
1801 0 : EFAError("fragment edge can't have >1 embedded nodes");
1802 : }
1803 : } // IF the input embedded node already exists on this elem edge
1804 :
1805 92206 : if (add_to_neighbor)
1806 : {
1807 92206 : for (unsigned int en_iter = 0; en_iter < numEdgeNeighbors(edge_id); ++en_iter)
1808 : {
1809 45023 : EFAElement2D * edge_neighbor = getEdgeNeighbor(edge_id, en_iter);
1810 45023 : unsigned int neighbor_edge_id = edge_neighbor->getNeighborIndex(this);
1811 45023 : if (edge_neighbor->getEdge(neighbor_edge_id)->getNode(0) == edge_node1) // same direction
1812 0 : EFAError("neighbor edge has the same direction as this edge");
1813 45023 : double neigh_pos = 1.0 - position; // get emb node's postion on neighbor edge
1814 45023 : edge_neighbor->addEdgeCut(neighbor_edge_id, neigh_pos, local_embedded, EmbeddedNodes, false);
1815 : }
1816 : } // If add_to_neighbor required
1817 92206 : }
1818 :
1819 : void
1820 136 : EFAElement2D::addNodeCut(unsigned int node_id,
1821 : EFANode * embedded_permanent_node,
1822 : std::map<unsigned int, EFANode *> & PermanentNodes,
1823 : std::map<unsigned int, EFANode *> & EmbeddedPermanentNodes)
1824 : {
1825 : EFANode * local_embedded_permanent = nullptr;
1826 136 : EFANode * node = _nodes[node_id];
1827 : if (embedded_permanent_node) // use the existing embedded node if it was passed in
1828 : local_embedded_permanent = embedded_permanent_node;
1829 :
1830 136 : if (node->category() == EFANode::N_CATEGORY_PERMANENT)
1831 : {
1832 60 : node->setCategory(EFANode::N_CATEGORY_EMBEDDED_PERMANENT);
1833 : local_embedded_permanent = node;
1834 60 : EmbeddedPermanentNodes.insert(std::make_pair(node->id(), local_embedded_permanent));
1835 60 : if (!Efa::deleteFromMap(PermanentNodes, local_embedded_permanent, false))
1836 0 : EFAError("Attempted to delete node: ",
1837 : local_embedded_permanent->id(),
1838 : " from PermanentNodes, but couldn't find it");
1839 : }
1840 136 : }
1841 :
1842 : bool
1843 35290 : EFAElement2D::addFragmentEdgeCut(unsigned int frag_edge_id,
1844 : double position,
1845 : std::map<unsigned int, EFANode *> & EmbeddedNodes)
1846 : {
1847 35290 : if (_fragments.size() != 1)
1848 0 : EFAError("Element: ", _id, " should have only 1 fragment in addFragEdgeIntersection");
1849 : EFANode * local_embedded = nullptr;
1850 :
1851 : // check if this intersection coincide with any embedded node on this edge
1852 : bool isValidIntersection = true;
1853 35290 : EFAEdge * frag_edge = getFragmentEdge(0, frag_edge_id); // we're considering this edge
1854 35290 : EFANode * edge_node1 = frag_edge->getNode(0);
1855 35290 : EFANode * edge_node2 = frag_edge->getNode(1);
1856 35290 : if ((std::abs(position) < Xfem::tol && edge_node1->category() == EFANode::N_CATEGORY_EMBEDDED) ||
1857 35365 : (std::abs(1.0 - position) < Xfem::tol &&
1858 17338 : edge_node2->category() == EFANode::N_CATEGORY_EMBEDDED))
1859 : isValidIntersection = false;
1860 :
1861 : // TODO: do not allow to cut fragment's node
1862 35290 : if (std::abs(position) < Xfem::tol || std::abs(1.0 - position) < Xfem::tol)
1863 : isValidIntersection = false;
1864 :
1865 : // add valid intersection point to an edge
1866 514 : if (isValidIntersection)
1867 : {
1868 514 : if (frag_edge->hasIntersection())
1869 : {
1870 416 : if (!frag_edge->hasIntersectionAtPosition(position, edge_node1))
1871 0 : EFAError("Attempting to add fragment edge intersection when one already exists with "
1872 : "different position.",
1873 : " elem: ",
1874 : _id,
1875 : " edge: ",
1876 : frag_edge_id,
1877 : " position: ",
1878 : position,
1879 : " old position: ",
1880 : frag_edge->getIntersection(0, edge_node1));
1881 : }
1882 : else // blank edge - in fact, it can only be a blank element interior edge
1883 : {
1884 196 : if (!_fragments[0]->isEdgeInterior(frag_edge_id) ||
1885 98 : _fragments[0]->isSecondaryInteriorEdge(frag_edge_id))
1886 0 : EFAError("Attemping to add intersection to an invalid fragment edge. Element: ",
1887 : _id,
1888 : " fragment_edge: ",
1889 : frag_edge_id);
1890 :
1891 : // create the embedded node and add it to the fragment's boundary edge
1892 98 : unsigned int new_node_id = Efa::getNewID(EmbeddedNodes);
1893 98 : local_embedded = new EFANode(new_node_id, EFANode::N_CATEGORY_EMBEDDED);
1894 98 : EmbeddedNodes.insert(std::make_pair(new_node_id, local_embedded));
1895 98 : frag_edge->addIntersection(position, local_embedded, edge_node1);
1896 :
1897 : // save this interior embedded node to FaceNodes
1898 : // TODO: for unstructured elements, the following calution gives you inaccurate position of
1899 : // face nodes
1900 : // must solve this issue for 3D!
1901 98 : std::vector<double> node1_para_coor(2, 0.0);
1902 98 : std::vector<double> node2_para_coor(2, 0.0);
1903 196 : if (getEdgeNodeParametricCoordinate(edge_node1, node1_para_coor) &&
1904 98 : getEdgeNodeParametricCoordinate(edge_node2, node2_para_coor))
1905 : {
1906 98 : double xi = (1.0 - position) * node1_para_coor[0] + position * node2_para_coor[0];
1907 98 : double eta = (1.0 - position) * node1_para_coor[1] + position * node2_para_coor[1];
1908 98 : _interior_nodes.push_back(new EFAFaceNode(local_embedded, xi, eta));
1909 : }
1910 : else
1911 0 : EFAError("elem: ", _id, " cannot get the parametric coords of two end embedded nodes");
1912 98 : }
1913 : // no need to add intersection for neighbor fragment - if this fragment has a
1914 : // neighbor fragment, the neighbor has already been treated in addEdgeIntersection;
1915 : // for an interior edge, there is no neighbor fragment
1916 : }
1917 :
1918 35290 : return isValidIntersection;
1919 : }
1920 :
1921 : std::vector<EFAFragment2D *>
1922 1 : EFAElement2D::branchingSplit(std::map<unsigned int, EFANode *> & EmbeddedNodes)
1923 : {
1924 1 : if (isPartial())
1925 0 : EFAError("branching is only allowed for an uncut element");
1926 :
1927 : // collect all emb nodes counterclockwise
1928 : std::vector<EFANode *> three_nodes;
1929 5 : for (unsigned int i = 0; i < _edges.size(); ++i)
1930 : {
1931 4 : EFANode * node1 = _edges[i]->getNode(0);
1932 4 : if (_edges[i]->numEmbeddedNodes() == 1)
1933 3 : three_nodes.push_back(_edges[i]->getEmbeddedNode(0));
1934 1 : else if (_edges[i]->numEmbeddedNodes() == 2)
1935 : {
1936 : unsigned int id0(
1937 0 : _edges[i]->getIntersection(0, node1) < _edges[i]->getIntersection(1, node1) ? 0 : 1);
1938 0 : unsigned int id1 = 1 - id0;
1939 0 : three_nodes.push_back(_edges[i]->getEmbeddedNode(id0));
1940 0 : three_nodes.push_back(_edges[i]->getEmbeddedNode(id1));
1941 : }
1942 : }
1943 1 : if (three_nodes.size() != 3)
1944 0 : EFAError("three_nodes.size() != 3");
1945 :
1946 : // get the parent coords of the braycenter of the three nodes
1947 : // TODO: may need a better way to compute this "branching point"
1948 1 : std::vector<double> center_xi(2, 0.0);
1949 4 : for (unsigned int i = 0; i < 3; ++i)
1950 : {
1951 3 : std::vector<double> xi_2d(2, 0.0);
1952 3 : getEdgeNodeParametricCoordinate(three_nodes[i], xi_2d);
1953 3 : center_xi[0] += xi_2d[0];
1954 3 : center_xi[1] += xi_2d[1];
1955 3 : }
1956 1 : center_xi[0] /= 3.0;
1957 1 : center_xi[1] /= 3.0;
1958 :
1959 : // create a new interior node for current element
1960 1 : unsigned int new_node_id = Efa::getNewID(EmbeddedNodes);
1961 1 : EFANode * new_emb = new EFANode(new_node_id, EFANode::N_CATEGORY_EMBEDDED);
1962 1 : EmbeddedNodes.insert(std::make_pair(new_node_id, new_emb));
1963 1 : _interior_nodes.push_back(new EFAFaceNode(new_emb, center_xi[0], center_xi[1]));
1964 :
1965 : // generate the three fragments
1966 : std::vector<EFAFragment2D *> new_fragments;
1967 4 : for (unsigned int i = 0; i < 3; ++i) // loop over 3 sectors
1968 : {
1969 3 : EFAFragment2D * new_frag = new EFAFragment2D(this, false, nullptr);
1970 3 : unsigned int iplus1(i < 2 ? i + 1 : 0);
1971 3 : new_frag->addEdge(new EFAEdge(three_nodes[iplus1], new_emb));
1972 3 : new_frag->addEdge(new EFAEdge(new_emb, three_nodes[i]));
1973 :
1974 : unsigned int iedge = 0;
1975 : bool add_more_edges = true;
1976 6 : for (unsigned int j = 0; j < _edges.size(); ++j)
1977 : {
1978 6 : if (_edges[j]->containsNode(three_nodes[i]))
1979 : {
1980 3 : if (_edges[j]->containsNode(three_nodes[iplus1]))
1981 : {
1982 0 : new_frag->addEdge(new EFAEdge(three_nodes[i], three_nodes[iplus1]));
1983 : add_more_edges = false;
1984 : }
1985 : else
1986 : {
1987 3 : new_frag->addEdge(new EFAEdge(three_nodes[i], _edges[j]->getNode(1)));
1988 : }
1989 : iedge = j;
1990 : break;
1991 : }
1992 : } // j
1993 4 : while (add_more_edges)
1994 : {
1995 4 : iedge += 1;
1996 4 : if (iedge == _edges.size())
1997 : iedge = 0;
1998 4 : if (_edges[iedge]->containsNode(three_nodes[iplus1]))
1999 : {
2000 3 : new_frag->addEdge(new EFAEdge(_edges[iedge]->getNode(0), three_nodes[iplus1]));
2001 : add_more_edges = false;
2002 : }
2003 : else
2004 1 : new_frag->addEdge(new EFAEdge(_edges[iedge]->getNode(0), _edges[iedge]->getNode(1)));
2005 : }
2006 3 : new_fragments.push_back(new_frag);
2007 : } // i
2008 1 : return new_fragments;
2009 1 : }
2010 :
2011 : void
2012 199 : EFAElement2D::mapParametricCoordFrom1Dto2D(unsigned int edge_id,
2013 : double xi_1d,
2014 : std::vector<double> & para_coor) const
2015 : {
2016 199 : para_coor.resize(2, 0.0);
2017 199 : if (_num_edges == 4)
2018 : {
2019 103 : if (edge_id == 0)
2020 : {
2021 51 : para_coor[0] = xi_1d;
2022 51 : para_coor[1] = -1.0;
2023 : }
2024 52 : else if (edge_id == 1)
2025 : {
2026 1 : para_coor[0] = 1.0;
2027 1 : para_coor[1] = xi_1d;
2028 : }
2029 51 : else if (edge_id == 2)
2030 : {
2031 3 : para_coor[0] = -xi_1d;
2032 3 : para_coor[1] = 1.0;
2033 : }
2034 48 : else if (edge_id == 3)
2035 : {
2036 48 : para_coor[0] = -1.0;
2037 48 : para_coor[1] = -xi_1d;
2038 : }
2039 : else
2040 0 : EFAError("edge_id out of bounds");
2041 : }
2042 96 : else if (_num_edges == 3)
2043 : {
2044 96 : if (edge_id == 0)
2045 : {
2046 48 : para_coor[0] = 0.5 * (1.0 - xi_1d);
2047 48 : para_coor[1] = 0.5 * (1.0 + xi_1d);
2048 : }
2049 48 : else if (edge_id == 1)
2050 : {
2051 0 : para_coor[0] = 0.0;
2052 0 : para_coor[1] = 0.5 * (1.0 - xi_1d);
2053 : }
2054 48 : else if (edge_id == 2)
2055 : {
2056 48 : para_coor[0] = 0.5 * (1.0 + xi_1d);
2057 48 : para_coor[1] = 0.0;
2058 : }
2059 : else
2060 0 : EFAError("edge_id out of bounds");
2061 : }
2062 : else
2063 0 : EFAError("unknown element for 2D");
2064 199 : }
2065 :
2066 : std::vector<EFANode *>
2067 4134502 : EFAElement2D::getCommonNodes(const EFAElement2D * other_elem) const
2068 : {
2069 : std::set<EFANode *> e1nodes(_nodes.begin(),
2070 4134502 : _nodes.begin() + _num_edges); // only account for corner nodes
2071 4134502 : std::set<EFANode *> e2nodes(other_elem->_nodes.begin(), other_elem->_nodes.begin() + _num_edges);
2072 4134502 : std::vector<EFANode *> common_nodes = Efa::getCommonElems(e1nodes, e2nodes);
2073 4134502 : return common_nodes;
2074 : }
|
