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