www.mooseframework.org
XFEM.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "XFEM.h"
11 
12 // XFEM includes
13 #include "XFEMAppTypes.h"
14 #include "XFEMCutElem2D.h"
15 #include "XFEMCutElem3D.h"
16 #include "XFEMFuncs.h"
17 #include "EFANode.h"
18 #include "EFAEdge.h"
19 #include "EFAFace.h"
20 #include "EFAFragment2D.h"
21 #include "EFAFragment3D.h"
22 #include "EFAFuncs.h"
23 
24 // MOOSE includes
25 #include "AuxiliarySystem.h"
26 #include "MooseVariable.h"
27 #include "NonlinearSystem.h"
28 #include "FEProblem.h"
29 #include "Assembly.h"
30 #include "MooseUtils.h"
32 
33 #include "libmesh/mesh_communication.h"
34 #include "libmesh/partitioner.h"
35 
36 XFEM::XFEM(const InputParameters & params)
37  : XFEMInterface(params), _efa_mesh(Moose::out), _debug_output_level(1)
38 {
39 #ifndef LIBMESH_ENABLE_UNIQUE_ID
40  mooseError("MOOSE requires unique ids to be enabled in libmesh (configure with "
41  "--enable-unique-id) to use XFEM!");
42 #endif
43  _has_secondary_cut = false;
44 }
45 
47 {
48  for (std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.begin();
49  cemit != _cut_elem_map.end();
50  ++cemit)
51  delete cemit->second;
52 }
53 
54 void
56 {
57  _geometric_cuts.push_back(geometric_cut);
58 
59  geometric_cut->setInterfaceID(_geometric_cuts.size() - 1);
60 
61  _geom_marker_id_map[geometric_cut] = _geometric_cuts.size() - 1;
62 }
63 
64 void
65 XFEM::getCrackTipOrigin(std::map<unsigned int, const Elem *> & elem_id_crack_tip,
66  std::vector<Point> & crack_front_points)
67 {
68  elem_id_crack_tip.clear();
69  crack_front_points.clear();
70  crack_front_points.resize(_elem_crack_origin_direction_map.size());
71 
72  unsigned int crack_tip_index = 0;
73  // This map is used to sort the order in _elem_crack_origin_direction_map such that every process
74  // has same order
75  std::map<unsigned int, const Elem *> elem_id_map;
76 
77  int m = -1;
78  for (std::map<const Elem *, std::vector<Point>>::iterator mit1 =
81  ++mit1)
82  {
83  unsigned int elem_id = mit1->first->id();
84  if (elem_id == std::numeric_limits<unsigned int>::max())
85  {
86  elem_id_map[m] = mit1->first;
87  m--;
88  }
89  else
90  elem_id_map[elem_id] = mit1->first;
91  }
92 
93  for (std::map<unsigned int, const Elem *>::iterator mit1 = elem_id_map.begin();
94  mit1 != elem_id_map.end();
95  mit1++)
96  {
97  const Elem * elem = mit1->second;
98  std::map<const Elem *, std::vector<Point>>::iterator mit2 =
100  if (mit2 != _elem_crack_origin_direction_map.end())
101  {
102  elem_id_crack_tip[crack_tip_index] = mit2->first;
103  crack_front_points[crack_tip_index] =
104  (mit2->second)[0]; // [0] stores origin coordinates and [1] stores direction
105  crack_tip_index++;
106  }
107  }
108 }
109 
110 void
111 XFEM::addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal)
112 {
113  Elem * elem = _mesh->elem_ptr(elem_id);
114  std::map<const Elem *, RealVectorValue>::iterator mit;
115  mit = _state_marked_elems.find(elem);
116  if (mit != _state_marked_elems.end())
117  mooseError(" ERROR: element ", elem->id(), " already marked for crack growth.");
118  _state_marked_elems[elem] = normal;
119 }
120 
121 void
122 XFEM::addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal, unsigned int marked_side)
123 {
124  addStateMarkedElem(elem_id, normal);
125  Elem * elem = _mesh->elem_ptr(elem_id);
126  std::map<const Elem *, unsigned int>::iterator mit;
127  mit = _state_marked_elem_sides.find(elem);
128  if (mit != _state_marked_elem_sides.end())
129  mooseError(" ERROR: side of element ", elem->id(), " already marked for crack initiation.");
130 
131  _state_marked_elem_sides[elem] = marked_side;
132 }
133 
134 void
135 XFEM::addStateMarkedFrag(unsigned int elem_id, RealVectorValue & normal)
136 {
137  addStateMarkedElem(elem_id, normal);
138  Elem * elem = _mesh->elem_ptr(elem_id);
139  std::set<const Elem *>::iterator mit;
140  mit = _state_marked_frags.find(elem);
141  if (mit != _state_marked_frags.end())
142  mooseError(
143  " ERROR: element ", elem->id(), " already marked for fragment-secondary crack initiation.");
144 
145  _state_marked_frags.insert(elem);
146 }
147 
148 void
150 {
151  _state_marked_elems.clear();
152  _state_marked_frags.clear();
153  _state_marked_elem_sides.clear();
154 }
155 
156 void
157 XFEM::addGeomMarkedElem2D(const unsigned int elem_id,
158  const Xfem::GeomMarkedElemInfo2D geom_info,
159  const unsigned int interface_id)
160 {
161  Elem * elem = _mesh->elem_ptr(elem_id);
162  _geom_marked_elems_2d[elem].push_back(geom_info);
163  _geom_marker_id_elems[interface_id].insert(elem_id);
164 }
165 
166 void
167 XFEM::addGeomMarkedElem3D(const unsigned int elem_id,
168  const Xfem::GeomMarkedElemInfo3D geom_info,
169  const unsigned int interface_id)
170 {
171  Elem * elem = _mesh->elem_ptr(elem_id);
172  _geom_marked_elems_3d[elem].push_back(geom_info);
173  _geom_marker_id_elems[interface_id].insert(elem_id);
174 }
175 
176 void
178 {
179  _geom_marked_elems_2d.clear();
180  _geom_marked_elems_3d.clear();
181 }
182 
183 void
185 {
187  std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
188  std::set<EFAElement *>::iterator sit;
189  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
190  {
191  if (_mesh->mesh_dimension() == 2)
192  {
193  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(*sit);
194  EFANode * tip_node = CEMElem->getTipEmbeddedNode();
195  unsigned int cts_id = CEMElem->getCrackTipSplitElementID();
196 
197  Point origin(0, 0, 0);
198  Point direction(0, 0, 0);
199 
200  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
201  it = _cut_elem_map.find(_mesh->elem_ptr(cts_id)->unique_id());
202  if (it != _cut_elem_map.end())
203  {
204  const XFEMCutElem * xfce = it->second;
205  const EFAElement * EFAelem = xfce->getEFAElement();
206  if (EFAelem->isPartial()) // exclude the full crack tip elements
207  {
208  xfce->getCrackTipOriginAndDirection(tip_node->id(), origin, direction);
209  }
210  }
211 
212  std::vector<Point> tip_data;
213  tip_data.push_back(origin);
214  tip_data.push_back(direction);
215  const Elem * elem = _mesh->elem_ptr((*sit)->id());
217  std::pair<const Elem *, std::vector<Point>>(elem, tip_data));
218  }
219  }
220 }
221 
222 bool
224 {
225  bool mesh_changed = false;
226 
227  mesh_changed = healMesh();
228 
229  if (mesh_changed)
230  buildEFAMesh();
231 
232  if (mesh_changed)
233  {
234  _mesh->update_parallel_id_counts();
235  MeshCommunication().make_elems_parallel_consistent(*_mesh);
236  MeshCommunication().make_nodes_parallel_consistent(*_mesh);
237  // _mesh->find_neighbors();
238  // _mesh->contract();
239  _mesh->allow_renumbering(false);
240  _mesh->skip_partitioning(true);
241  _mesh->prepare_for_use();
242 
243  if (_displaced_mesh)
244  {
245  _displaced_mesh->update_parallel_id_counts();
246  MeshCommunication().make_elems_parallel_consistent(*_displaced_mesh);
247  MeshCommunication().make_nodes_parallel_consistent(*_displaced_mesh);
248  _displaced_mesh->allow_renumbering(false);
249  _displaced_mesh->skip_partitioning(true);
250  _displaced_mesh->prepare_for_use();
251  }
252  }
253 
254  _geom_marker_id_elems.clear();
255 
256  return mesh_changed;
257 }
258 
259 bool
260 XFEM::update(Real time,
261  const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
262  AuxiliarySystem & aux)
263 {
265  mooseError("Use of XFEM with distributed mesh is not yet supported");
266 
267  bool mesh_changed = false;
268 
269  buildEFAMesh();
270 
272 
274 
275  if (markCuts(time))
276  mesh_changed = cutMeshWithEFA(nl, aux);
277 
278  if (mesh_changed)
279  {
280  buildEFAMesh();
282  }
283 
284  if (mesh_changed)
285  {
286  _mesh->allow_renumbering(false);
287  _mesh->skip_partitioning(true);
288  _mesh->prepare_for_use();
289 
290  if (_displaced_mesh)
291  {
292  _displaced_mesh->allow_renumbering(false);
293  _displaced_mesh->skip_partitioning(true);
294  _displaced_mesh->prepare_for_use();
295  }
296  }
297 
300 
301  return mesh_changed;
302 }
303 
304 void
305 XFEM::initSolution(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nls,
306  AuxiliarySystem & aux)
307 {
308  if (nls.size() != 1)
309  mooseError("XFEM does not currently support multiple nonlinear systems");
310 
311  nls[0]->serializeSolution();
312  aux.serializeSolution();
313  NumericVector<Number> & current_solution = *nls[0]->system().current_local_solution;
314  NumericVector<Number> & old_solution = nls[0]->solutionOld();
315  NumericVector<Number> & older_solution = nls[0]->solutionOlder();
316  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
317  NumericVector<Number> & old_aux_solution = aux.solutionOld();
318  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
319 
320  setSolution(*nls[0], _cached_solution, current_solution, old_solution, older_solution);
321  setSolution(
322  aux, _cached_aux_solution, current_aux_solution, old_aux_solution, older_aux_solution);
323 
324  current_solution.close();
325  old_solution.close();
326  older_solution.close();
327  current_aux_solution.close();
328  old_aux_solution.close();
329  older_aux_solution.close();
330 
331  _cached_solution.clear();
332  _cached_aux_solution.clear();
333 }
334 
335 void
337 {
338  _efa_mesh.reset();
339 
340  // Load all existing elements in to EFA mesh
341  for (auto & elem : _mesh->element_ptr_range())
342  {
343  std::vector<unsigned int> quad;
344  for (unsigned int i = 0; i < elem->n_nodes(); ++i)
345  quad.push_back(elem->node_id(i));
346 
347  if (_mesh->mesh_dimension() == 2)
348  _efa_mesh.add2DElement(quad, elem->id());
349  else if (_mesh->mesh_dimension() == 3)
350  _efa_mesh.add3DElement(quad, elem->id());
351  else
352  mooseError("XFEM only works for 2D and 3D");
353  }
354 
355  // Restore fragment information for elements that have been previously cut
356  for (auto & elem : _mesh->element_ptr_range())
357  {
358  std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.find(elem->unique_id());
359  if (cemit != _cut_elem_map.end())
360  {
361  XFEMCutElem * xfce = cemit->second;
362  EFAElement * CEMElem = _efa_mesh.getElemByID(elem->id());
363  _efa_mesh.restoreFragmentInfo(CEMElem, xfce->getEFAElement());
364  }
365  }
366 
367  // Must update edge neighbors before restore edge intersections. Otherwise, when we
368  // add edge intersections, we do not have neighbor information to use.
369  // Correction: no need to use neighbor info now
372 }
373 
374 bool
375 XFEM::markCuts(Real time)
376 {
377  bool marked_sides = false;
378  if (_mesh->mesh_dimension() == 2)
379  {
380  marked_sides = markCutEdgesByGeometry();
381  marked_sides |= markCutEdgesByState(time);
382  }
383  else if (_mesh->mesh_dimension() == 3)
384  {
385  marked_sides = markCutFacesByGeometry();
386  marked_sides |= markCutFacesByState();
387  }
388  return marked_sides;
389 }
390 
391 bool
393 {
394  bool marked_edges = false;
395  bool marked_nodes = false;
396 
397  for (const auto & gme : _geom_marked_elems_2d)
398  {
399  for (const auto & gmei : gme.second)
400  {
401  EFAElement2D * EFAElem = getEFAElem2D(gme.first);
402 
403  for (unsigned int i = 0; i < gmei._elem_cut_edges.size(); ++i) // mark element edges
404  {
405  if (!EFAElem->isEdgePhantom(
406  gmei._elem_cut_edges[i]._host_side_id)) // must not be phantom edge
407  {
408  _efa_mesh.addElemEdgeIntersection(gme.first->id(),
409  gmei._elem_cut_edges[i]._host_side_id,
410  gmei._elem_cut_edges[i]._distance);
411  marked_edges = true;
412  }
413  }
414 
415  for (unsigned int i = 0; i < gmei._elem_cut_nodes.size(); ++i) // mark element edges
416  {
417  _efa_mesh.addElemNodeIntersection(gme.first->id(), gmei._elem_cut_nodes[i]._host_id);
418  marked_nodes = true;
419  }
420 
421  for (unsigned int i = 0; i < gmei._frag_cut_edges.size();
422  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
423  {
424  if (!EFAElem->getFragment(0)->isSecondaryInteriorEdge(
425  gmei._frag_cut_edges[i]._host_side_id))
426  {
427  if (_efa_mesh.addFragEdgeIntersection(gme.first->id(),
428  gmei._frag_cut_edges[i]._host_side_id,
429  gmei._frag_cut_edges[i]._distance))
430  {
431  marked_edges = true;
432  if (!isElemAtCrackTip(gme.first))
433  _has_secondary_cut = true;
434  }
435  }
436  }
437  }
438  }
439 
440  return marked_edges || marked_nodes;
441 }
442 
443 void
445  EFAElement2D * CEMElem,
446  EFAEdge * orig_edge,
447  Point normal,
448  Point crack_tip_origin,
449  Point crack_tip_direction,
450  Real & distance_keep,
451  unsigned int & edge_id_keep,
452  Point & normal_keep)
453 {
454  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
455  Point edge1(0.0, 0.0, 0.0);
456  Point edge2(0.0, 0.0, 0.0);
457  Point left_angle(0.0, 0.0, 0.0);
458  Point right_angle(0.0, 0.0, 0.0);
459  Point left_angle_normal(0.0, 0.0, 0.0);
460  Point right_angle_normal(0.0, 0.0, 0.0);
461  Point crack_direction_normal(0.0, 0.0, 0.0);
462  Point edge1_to_tip(0.0, 0.0, 0.0);
463  Point edge2_to_tip(0.0, 0.0, 0.0);
464  Point edge1_to_tip_normal(0.0, 0.0, 0.0);
465  Point edge2_to_tip_normal(0.0, 0.0, 0.0);
466 
467  Real cos_45 = std::cos(45.0 / 180.0 * 3.14159);
468  Real sin_45 = std::sin(45.0 / 180.0 * 3.14159);
469 
470  left_angle(0) = cos_45 * crack_tip_direction(0) - sin_45 * crack_tip_direction(1);
471  left_angle(1) = sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
472 
473  right_angle(0) = cos_45 * crack_tip_direction(0) + sin_45 * crack_tip_direction(1);
474  right_angle(1) = -sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
475 
476  left_angle_normal(0) = -left_angle(1);
477  left_angle_normal(1) = left_angle(0);
478 
479  right_angle_normal(0) = -right_angle(1);
480  right_angle_normal(1) = right_angle(0);
481 
482  crack_direction_normal(0) = -crack_tip_direction(1);
483  crack_direction_normal(1) = crack_tip_direction(0);
484 
485  Real angle_min = 0.0;
486  Real distance = 0.0;
487  unsigned int nsides = CEMElem->numEdges();
488 
489  for (unsigned int i = 0; i < nsides; ++i)
490  {
491  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
492  {
493  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
494  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
495 
496  edge1_to_tip = (edge_ends[0] * 0.95 + edge_ends[1] * 0.05) - crack_tip_origin;
497  edge2_to_tip = (edge_ends[0] * 0.05 + edge_ends[1] * 0.95) - crack_tip_origin;
498 
499  edge1_to_tip /= pow(edge1_to_tip.norm_sq(), 0.5);
500  edge2_to_tip /= pow(edge2_to_tip.norm_sq(), 0.5);
501 
502  edge1_to_tip_normal(0) = -edge1_to_tip(1);
503  edge1_to_tip_normal(1) = edge1_to_tip(0);
504 
505  edge2_to_tip_normal(0) = -edge2_to_tip(1);
506  edge2_to_tip_normal(1) = edge2_to_tip(0);
507 
508  Real angle_edge1_normal = edge1_to_tip_normal * normal;
509  Real angle_edge2_normal = edge2_to_tip_normal * normal;
510 
511  if (std::abs(angle_edge1_normal) > std::abs(angle_min) &&
512  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
513  {
514  edge_id_keep = i;
515  distance_keep = 0.05;
516  normal_keep = edge1_to_tip_normal;
517  angle_min = angle_edge1_normal;
518  }
519  else if (std::abs(angle_edge2_normal) > std::abs(angle_min) &&
520  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
521  {
522  edge_id_keep = i;
523  distance_keep = 0.95;
524  normal_keep = edge2_to_tip_normal;
525  angle_min = angle_edge2_normal;
526  }
527 
529  crack_tip_origin, left_angle_normal, edge_ends[0], edge_ends[1], distance) &&
530  (!CEMElem->isEdgePhantom(i)))
531  {
532  if (std::abs(left_angle_normal * normal) > std::abs(angle_min) &&
533  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
534  {
535  edge_id_keep = i;
536  distance_keep = distance;
537  normal_keep = left_angle_normal;
538  angle_min = left_angle_normal * normal;
539  }
540  }
541  else if (initCutIntersectionEdge(
542  crack_tip_origin, right_angle_normal, edge_ends[0], edge_ends[1], distance) &&
543  (!CEMElem->isEdgePhantom(i)))
544  {
545  if (std::abs(right_angle_normal * normal) > std::abs(angle_min) &&
546  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
547  {
548  edge_id_keep = i;
549  distance_keep = distance;
550  normal_keep = right_angle_normal;
551  angle_min = right_angle_normal * normal;
552  }
553  }
554  else if (initCutIntersectionEdge(crack_tip_origin,
555  crack_direction_normal,
556  edge_ends[0],
557  edge_ends[1],
558  distance) &&
559  (!CEMElem->isEdgePhantom(i)))
560  {
561  if (std::abs(crack_direction_normal * normal) > std::abs(angle_min) &&
562  (crack_tip_direction * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
563  {
564  edge_id_keep = i;
565  distance_keep = distance;
566  normal_keep = crack_direction_normal;
567  angle_min = crack_direction_normal * normal;
568  }
569  }
570  }
571  }
572 
573  // avoid small volume fraction cut
574  if ((distance_keep - 0.05) < 0.0)
575  {
576  distance_keep = 0.05;
577  }
578  else if ((distance_keep - 0.95) > 0.0)
579  {
580  distance_keep = 0.95;
581  }
582 }
583 
584 bool
586 {
587  bool marked_edges = false;
588  for (std::map<const Elem *, RealVectorValue>::iterator pmeit = _state_marked_elems.begin();
589  pmeit != _state_marked_elems.end();
590  ++pmeit)
591  {
592  const Elem * elem = pmeit->first;
593  RealVectorValue normal = pmeit->second;
594  EFAElement2D * CEMElem = getEFAElem2D(elem);
595 
596  Real volfrac_elem = getPhysicalVolumeFraction(elem);
597  if (volfrac_elem < 0.25)
598  continue;
599 
600  // continue if elem is already cut twice - IMPORTANT
601  if (CEMElem->isFinalCut())
602  continue;
603 
604  // find the first cut edge
605  unsigned int nsides = CEMElem->numEdges();
606  unsigned int orig_cut_side_id = std::numeric_limits<unsigned int>::max();
607  Real orig_cut_distance = -1.0;
608  EFANode * orig_node = nullptr;
609  EFAEdge * orig_edge = nullptr;
610 
611  // crack tip origin coordinates and direction
612  Point crack_tip_origin(0, 0, 0);
613  Point crack_tip_direction(0, 0, 0);
614 
615  if (isElemAtCrackTip(elem)) // crack tip element's crack intiation
616  {
617  orig_cut_side_id = CEMElem->getTipEdgeID();
618  if (orig_cut_side_id < nsides) // valid crack-tip edge found
619  {
620  orig_edge = CEMElem->getEdge(orig_cut_side_id);
621  orig_node = CEMElem->getTipEmbeddedNode();
622  }
623  else
624  mooseError("element ", elem->id(), " has no valid crack-tip edge");
625 
626  // obtain the crack tip origin coordinates and direction.
627  std::map<const Elem *, std::vector<Point>>::iterator ecodm =
629  if (ecodm != _elem_crack_origin_direction_map.end())
630  {
631  crack_tip_origin = (ecodm->second)[0];
632  crack_tip_direction = (ecodm->second)[1];
633  }
634  else
635  mooseError("element ", elem->id(), " cannot find its crack tip origin and direction.");
636  }
637  else
638  {
639  std::map<const Elem *, unsigned int>::iterator mit1;
640  mit1 = _state_marked_elem_sides.find(elem);
641  std::set<const Elem *>::iterator mit2;
642  mit2 = _state_marked_frags.find(elem);
643 
644  if (mit1 != _state_marked_elem_sides.end()) // specified boundary crack initiation
645  {
646  orig_cut_side_id = mit1->second;
647  if (!CEMElem->isEdgePhantom(orig_cut_side_id) &&
648  !CEMElem->getEdge(orig_cut_side_id)->hasIntersection())
649  {
650  orig_cut_distance = 0.5;
651  _efa_mesh.addElemEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
652  orig_edge = CEMElem->getEdge(orig_cut_side_id);
653  orig_node = orig_edge->getEmbeddedNode(0);
654  // get a virtual crack tip direction
655  Point elem_center(0.0, 0.0, 0.0);
656  Point edge_center;
657  for (unsigned int i = 0; i < nsides; ++i)
658  {
659  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
660  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
661  }
662  elem_center /= nsides * 2.0;
663  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
664  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
665  edge_center /= 2.0;
666  crack_tip_origin = edge_center;
667  crack_tip_direction = elem_center - edge_center;
668  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
669  }
670  else
671  continue; // skip this elem if specified boundary edge is phantom
672  }
673  else if (mit2 != _state_marked_frags.end()) // cut-surface secondary crack initiation
674  {
675  if (CEMElem->numFragments() != 1)
676  mooseError("element ",
677  elem->id(),
678  " flagged for a secondary crack, but has ",
679  CEMElem->numFragments(),
680  " fragments");
681  std::vector<unsigned int> interior_edge_id = CEMElem->getFragment(0)->getInteriorEdgeID();
682  if (interior_edge_id.size() == 1)
683  orig_cut_side_id = interior_edge_id[0];
684  else
685  continue; // skip this elem if more than one interior edges found (i.e. elem's been cut
686  // twice)
687  orig_cut_distance = 0.5;
688  _efa_mesh.addFragEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
689  orig_edge = CEMElem->getFragmentEdge(0, orig_cut_side_id);
690  orig_node = orig_edge->getEmbeddedNode(0); // must be an interior embedded node
691  Point elem_center(0.0, 0.0, 0.0);
692  Point edge_center;
693  unsigned int nsides_frag = CEMElem->getFragment(0)->numEdges();
694  for (unsigned int i = 0; i < nsides_frag; ++i)
695  {
696  elem_center +=
697  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
698  elem_center +=
699  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
700  }
701  elem_center /= nsides_frag * 2.0;
702  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
703  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
704  edge_center /= 2.0;
705  crack_tip_origin = edge_center;
706  crack_tip_direction = elem_center - edge_center;
707  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
708  }
709  else
710  mooseError("element ",
711  elem->id(),
712  " flagged for state-based growth, but has no edge intersections");
713  }
714 
715  Point cut_origin(0.0, 0.0, 0.0);
716  if (orig_node)
717  cut_origin = getEFANodeCoords(orig_node, CEMElem, elem); // cutting plane origin's coords
718  else
719  mooseError("element ", elem->id(), " does not have valid orig_node");
720 
721  // loop through element edges to add possible second cut points
722  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
723  Point edge1(0.0, 0.0, 0.0);
724  Point edge2(0.0, 0.0, 0.0);
725  Point cut_edge_point(0.0, 0.0, 0.0);
726  bool find_compatible_direction = false;
727  unsigned int edge_id_keep = 0;
728  Real distance_keep = 0.0;
729  Point normal_keep(0.0, 0.0, 0.0);
730  Real distance = 0.0;
731  bool edge_cut = false;
732 
733  for (unsigned int i = 0; i < nsides; ++i)
734  {
735  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
736  {
737  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
738  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
740  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
741  (!CEMElem->isEdgePhantom(i))))
742  {
743  cut_edge_point = distance * edge_ends[1] + (1.0 - distance) * edge_ends[0];
744  distance_keep = distance;
745  edge_id_keep = i;
746  normal_keep = normal;
747  edge_cut = true;
748  break;
749  }
750  }
751  }
752 
753  Point between_two_cuts = (cut_edge_point - crack_tip_origin);
754  between_two_cuts /= pow(between_two_cuts.norm_sq(), 0.5);
755  Real angle_between_two_cuts = between_two_cuts * crack_tip_direction;
756 
757  if (angle_between_two_cuts > std::cos(45.0 / 180.0 * 3.14159)) // original cut direction is good
758  find_compatible_direction = true;
759 
760  if (!find_compatible_direction && edge_cut)
762  CEMElem,
763  orig_edge,
764  normal,
765  crack_tip_origin,
766  crack_tip_direction,
767  distance_keep,
768  edge_id_keep,
769  normal_keep);
770 
771  if (edge_cut)
772  {
774  _efa_mesh.addElemEdgeIntersection(elem->id(), edge_id_keep, distance_keep);
775  else
776  {
777  Point growth_direction(0.0, 0.0, 0.0);
778 
779  growth_direction(0) = -normal_keep(1);
780  growth_direction(1) = normal_keep(0);
781 
782  if (growth_direction * crack_tip_direction < 1.0e-10)
783  growth_direction *= (-1.0);
784 
785  Real x0 = crack_tip_origin(0);
786  Real y0 = crack_tip_origin(1);
787  Real x1 = x0 + _crack_growth_increment * growth_direction(0);
788  Real y1 = y0 + _crack_growth_increment * growth_direction(1);
789 
790  XFEMCrackGrowthIncrement2DCut geometric_cut(x0, y0, x1, y1, time * 0.9, time * 0.9);
791 
792  for (const auto & elem : _mesh->element_ptr_range())
793  {
794  std::vector<CutEdgeForCrackGrowthIncr> elem_cut_edges;
795  EFAElement2D * CEMElem = getEFAElem2D(elem);
796 
797  // continue if elem has been already cut twice - IMPORTANT
798  if (CEMElem->isFinalCut())
799  continue;
800 
801  // mark cut edges for the element and its fragment
802  geometric_cut.cutElementByCrackGrowthIncrement(elem, elem_cut_edges, time);
803 
804  for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
805  {
806  if (!CEMElem->isEdgePhantom(
807  elem_cut_edges[i]._host_side_id)) // must not be phantom edge
808  {
810  elem->id(), elem_cut_edges[i]._host_side_id, elem_cut_edges[i]._distance);
811  }
812  }
813  }
814  }
815  }
816  // loop though framgent boundary edges to add possible second cut points
817  // N.B. must do this after marking element edges
818  if (CEMElem->numFragments() > 0 && !edge_cut)
819  {
820  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
821  {
822  if (!orig_edge->isPartialOverlap(*CEMElem->getFragmentEdge(0, i)))
823  {
824  edge_ends[0] =
825  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
826  edge_ends[1] =
827  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
829  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
830  (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(i)))
831  {
832  if (_efa_mesh.addFragEdgeIntersection(elem->id(), edge_id_keep, distance_keep))
833  if (!isElemAtCrackTip(elem))
834  _has_secondary_cut = true;
835  break;
836  }
837  }
838  }
839  }
840 
841  marked_edges = true;
842 
843  } // loop over all state_marked_elems
844 
845  return marked_edges;
846 }
847 
848 bool
850 {
851  bool marked_faces = false;
852 
853  for (const auto & gme : _geom_marked_elems_3d)
854  {
855  for (const auto & gmei : gme.second)
856  {
857  EFAElement3D * EFAElem = getEFAElem3D(gme.first);
858 
859  for (unsigned int i = 0; i < gmei._elem_cut_faces.size(); ++i) // mark element faces
860  {
861  if (!EFAElem->isFacePhantom(gmei._elem_cut_faces[i]._face_id)) // must not be phantom face
862  {
863  _efa_mesh.addElemFaceIntersection(gme.first->id(),
864  gmei._elem_cut_faces[i]._face_id,
865  gmei._elem_cut_faces[i]._face_edge,
866  gmei._elem_cut_faces[i]._position);
867  marked_faces = true;
868  }
869  }
870 
871  for (unsigned int i = 0; i < gmei._frag_cut_faces.size();
872  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
873  {
874  if (!EFAElem->getFragment(0)->isThirdInteriorFace(gmei._frag_cut_faces[i]._face_id))
875  {
876  _efa_mesh.addFragFaceIntersection(gme.first->id(),
877  gmei._frag_cut_faces[i]._face_id,
878  gmei._frag_cut_faces[i]._face_edge,
879  gmei._frag_cut_faces[i]._position);
880  marked_faces = true;
881  }
882  }
883  }
884  }
885 
886  return marked_faces;
887 }
888 
889 bool
891 {
892  bool marked_faces = false;
893  // TODO: need to finish this for 3D problems
894  return marked_faces;
895 }
896 
897 bool
899  Point cut_origin, RealVectorValue cut_normal, Point & edge_p1, Point & edge_p2, Real & dist)
900 {
901  dist = 0.0;
902  bool does_intersect = false;
903  Point origin2p1 = edge_p1 - cut_origin;
904  Real plane2p1 = cut_normal(0) * origin2p1(0) + cut_normal(1) * origin2p1(1);
905  Point origin2p2 = edge_p2 - cut_origin;
906  Real plane2p2 = cut_normal(0) * origin2p2(0) + cut_normal(1) * origin2p2(1);
907 
908  if (plane2p1 * plane2p2 < 0.0)
909  {
910  dist = -plane2p1 / (plane2p2 - plane2p1);
911  does_intersect = true;
912  }
913  return does_intersect;
914 }
915 
916 bool
918 {
919  bool mesh_changed = false;
920 
921  std::set<Node *> nodes_to_delete;
922  std::set<Node *> nodes_to_delete_displaced;
923  std::set<unsigned int> cutelems_to_delete;
924  unsigned int deleted_elem_count = 0;
925  std::vector<std::string> healed_geometric_cuts;
926 
927  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
928  {
929  if (_geometric_cuts[i]->shouldHealMesh())
930  {
931  healed_geometric_cuts.push_back(_geometric_cuts[i]->name());
932  for (auto & it : _sibling_elems[_geometric_cuts[i]->getInterfaceID()])
933  {
934  Elem * elem1 = const_cast<Elem *>(it.first);
935  Elem * elem2 = const_cast<Elem *>(it.second);
936 
937  std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
938  _cut_elem_map.find(elem1->unique_id());
939  if (cemit != _cut_elem_map.end())
940  {
941  const XFEMCutElem * xfce = cemit->second;
942 
943  cutelems_to_delete.insert(elem1->unique_id());
944 
945  for (unsigned int in = 0; in < elem1->n_nodes(); ++in)
946  {
947  Node * e1node = elem1->node_ptr(in);
948  Node * e2node = elem2->node_ptr(in);
949  if (!xfce->isPointPhysical(*e1node) &&
950  e1node != e2node) // This would happen at the crack tip
951  {
952  elem1->set_node(in) = e2node;
953  nodes_to_delete.insert(e1node);
954  }
955  else if (e1node != e2node)
956  nodes_to_delete.insert(e2node);
957  }
958  }
959  else
960  mooseError("Could not find XFEMCutElem for element to be kept in healing");
961 
962  // Store the material properties of the elements to be healed. So that if the element is
963  // immediately re-cut, we can restore the material properties (especially those stateful
964  // ones).
965  std::vector<const Elem *> healed_elems = {elem1, elem2};
966 
967  if (_geometric_cuts[i]->shouldHealMesh())
968  // If the parent element will not be re-cut, then all of its nodes must have the same
969  // CutSubdomainID. Therefore, just query the first node in this parent element to
970  // get its CutSubdomainID.
971  for (auto e : healed_elems)
972  if (elem1->processor_id() == _mesh->processor_id() &&
973  e->processor_id() == _mesh->processor_id())
974  {
975  storeMaterialPropertiesForElement(/*parent_elem = */ elem1, /*child_elem = */ e);
976  // In case the healed element is not re-cut, copy the corresponding material
977  // properties to the parent element now than later.
978  CutSubdomainID parent_gcsid =
979  _geometric_cuts[i]->getCutSubdomainID(elem1->node_ptr(0));
980  CutSubdomainID gcsid = _geom_cut_elems[e]._cut_subdomain_id;
981  if (parent_gcsid == gcsid)
983  }
984 
985  if (_displaced_mesh)
986  {
987  Elem * elem1_displaced = _displaced_mesh->elem_ptr(it.first->id());
988  Elem * elem2_displaced = _displaced_mesh->elem_ptr(it.second->id());
989 
990  std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
991  _cut_elem_map.find(elem1_displaced->unique_id());
992  if (cemit != _cut_elem_map.end())
993  {
994  const XFEMCutElem * xfce = cemit->second;
995 
996  for (unsigned int in = 0; in < elem1_displaced->n_nodes(); ++in)
997  {
998  Node * e1node_displaced = elem1_displaced->node_ptr(in);
999  Node * e2node_displaced = elem2_displaced->node_ptr(in);
1000  if (!xfce->isPointPhysical(*elem1->node_ptr(in)) &&
1001  e1node_displaced != e2node_displaced)
1002  {
1003  elem1_displaced->set_node(in) = e2node_displaced;
1004  nodes_to_delete_displaced.insert(e1node_displaced);
1005  }
1006  else if (e1node_displaced != e2node_displaced)
1007  nodes_to_delete_displaced.insert(e2node_displaced);
1008  }
1009  }
1010  else
1011  mooseError("Could not find XFEMCutElem for element to be kept in healing");
1012 
1013  elem2_displaced->nullify_neighbors();
1014  _displaced_mesh->get_boundary_info().remove(elem2_displaced);
1015  _displaced_mesh->delete_elem(elem2_displaced);
1016  }
1017 
1018  // remove the property storage of deleted element/side
1019  _material_data[0]->eraseProperty(elem2);
1020  _bnd_material_data[0]->eraseProperty(elem2);
1021 
1022  cutelems_to_delete.insert(elem2->unique_id());
1023  elem2->nullify_neighbors();
1024  _mesh->get_boundary_info().remove(elem2);
1025  unsigned int deleted_elem_id = elem2->id();
1026  _mesh->delete_elem(elem2);
1027  if (_debug_output_level > 1)
1028  {
1029  if (deleted_elem_count == 0)
1030  _console << "\n";
1031  _console << "XFEM healing deleted element: " << deleted_elem_id << std::endl;
1032  }
1033  ++deleted_elem_count;
1034  mesh_changed = true;
1035  }
1036  }
1037  }
1038 
1039  for (auto & sit : nodes_to_delete)
1040  {
1041  Node * node_to_delete = sit;
1042  dof_id_type deleted_node_id = node_to_delete->id();
1043  _mesh->get_boundary_info().remove(node_to_delete);
1044  _mesh->delete_node(node_to_delete);
1045  if (_debug_output_level > 1)
1046  _console << "XFEM healing deleted node: " << deleted_node_id << std::endl;
1047  }
1048 
1049  if (_displaced_mesh)
1050  {
1051  for (auto & sit : nodes_to_delete_displaced)
1052  {
1053  Node * node_to_delete_displaced = sit;
1054  _displaced_mesh->get_boundary_info().remove(node_to_delete_displaced);
1055  _displaced_mesh->delete_node(node_to_delete_displaced);
1056  }
1057  }
1058 
1059  for (auto & ced : cutelems_to_delete)
1060  if (_cut_elem_map.find(ced) != _cut_elem_map.end())
1061  {
1062  delete _cut_elem_map.find(ced)->second;
1063  _cut_elem_map.erase(ced);
1064  }
1065 
1066  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
1067  if (_geometric_cuts[i]->shouldHealMesh())
1068  _sibling_elems[_geometric_cuts[i]->getInterfaceID()].clear();
1069 
1070  if (_displaced_mesh)
1071  {
1072  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
1073  if (_geometric_cuts[i]->shouldHealMesh())
1074  _sibling_displaced_elems[_geometric_cuts[i]->getInterfaceID()].clear();
1075  }
1076 
1077  for (auto & ceh : _crack_tip_elems_to_be_healed)
1078  {
1079  _crack_tip_elems.erase(ceh);
1081  delete _cut_elem_map.find(ceh->unique_id())->second;
1082  _cut_elem_map.erase(ceh->unique_id());
1083  }
1084 
1085  if (!healed_geometric_cuts.empty() && _debug_output_level > 0)
1086  {
1087  _console << "\nXFEM mesh healing complete\n";
1088  _console << "Names of healed geometric cut objects: ";
1089  for (auto geomcut : healed_geometric_cuts)
1090  _console << geomcut << " ";
1091  _console << "\n";
1092  _console << "# deleted nodes: " << nodes_to_delete.size() << "\n";
1093  _console << "# deleted elements: " << deleted_elem_count << "\n";
1094  _console << std::flush;
1095  }
1096 
1097  return mesh_changed;
1098 }
1099 
1100 bool
1101 XFEM::cutMeshWithEFA(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nls,
1102  AuxiliarySystem & aux)
1103 {
1104  if (nls.size() != 1)
1105  mooseError("XFEM does not currently support multiple nonlinear systems");
1106 
1107  std::map<unsigned int, Node *> efa_id_to_new_node;
1108  std::map<unsigned int, Node *> efa_id_to_new_node2;
1109  std::map<unsigned int, Elem *> efa_id_to_new_elem;
1110  _cached_solution.clear();
1111  _cached_aux_solution.clear();
1112 
1113  // Copy the current geometric cut element info (from last time) into the
1114  // _old_geom_cut_elems.
1116  _geom_cut_elems.clear();
1117 
1119 
1120  if (_debug_output_level > 2)
1121  {
1122  _console << "\nXFEM Element fragment algorithm mesh prior to cutting:\n";
1123  _console << std::flush;
1124  _efa_mesh.printMesh();
1125  }
1126 
1128 
1129  if (_debug_output_level > 2)
1130  {
1131  _console << "\nXFEM Element fragment algorithm mesh after cutting:\n";
1132  _console << std::flush;
1133  _efa_mesh.printMesh();
1134  }
1135 
1136  const std::vector<EFANode *> new_nodes = _efa_mesh.getNewNodes();
1137  const std::vector<EFAElement *> new_elements = _efa_mesh.getChildElements();
1138  const std::vector<EFAElement *> delete_elements = _efa_mesh.getParentElements();
1139 
1140  bool mesh_changed = (new_nodes.size() + new_elements.size() + delete_elements.size() > 0);
1141 
1142  // Prepare to cache solution on DOFs modified by XFEM
1143  if (mesh_changed)
1144  {
1145  nls[0]->serializeSolution();
1146  aux.serializeSolution();
1147  if (_debug_output_level > 1)
1148  _console << "\n";
1149  }
1150  NumericVector<Number> & current_solution = *nls[0]->system().current_local_solution;
1151  NumericVector<Number> & old_solution = nls[0]->solutionOld();
1152  NumericVector<Number> & older_solution = nls[0]->solutionOlder();
1153  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
1154  NumericVector<Number> & old_aux_solution = aux.solutionOld();
1155  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
1156 
1157  std::map<Node *, Node *> new_nodes_to_parents;
1158 
1159  // Add new nodes
1160  for (unsigned int i = 0; i < new_nodes.size(); ++i)
1161  {
1162  unsigned int new_node_id = new_nodes[i]->id();
1163  unsigned int parent_id = new_nodes[i]->parent()->id();
1164 
1165  Node * parent_node = _mesh->node_ptr(parent_id);
1166  Node * new_node = Node::build(*parent_node, _mesh->n_nodes()).release();
1167  _mesh->add_node(new_node);
1168 
1169  new_nodes_to_parents[new_node] = parent_node;
1170 
1171  new_node->set_n_systems(parent_node->n_systems());
1172  efa_id_to_new_node.insert(std::make_pair(new_node_id, new_node));
1173  if (_debug_output_level > 1)
1174  _console << "XFEM added new node: " << new_node->id() << std::endl;
1175  if (_displaced_mesh)
1176  {
1177  const Node * parent_node2 = _displaced_mesh->node_ptr(parent_id);
1178  Node * new_node2 = Node::build(*parent_node2, _displaced_mesh->n_nodes()).release();
1179  _displaced_mesh->add_node(new_node2);
1180 
1181  new_node2->set_n_systems(parent_node2->n_systems());
1182  efa_id_to_new_node2.insert(std::make_pair(new_node_id, new_node2));
1183  }
1184  }
1185 
1186  // Add new elements
1187  std::map<unsigned int, std::vector<const Elem *>> temporary_parent_children_map;
1188 
1189  std::vector<boundary_id_type> parent_boundary_ids;
1190 
1191  for (unsigned int i = 0; i < new_elements.size(); ++i)
1192  {
1193  unsigned int parent_id = new_elements[i]->getParent()->id();
1194  unsigned int efa_child_id = new_elements[i]->id();
1195 
1196  Elem * parent_elem = _mesh->elem_ptr(parent_id);
1197  Elem * libmesh_elem = Elem::build(parent_elem->type()).release();
1198 
1199  for (unsigned int m = 0; m < _geometric_cuts.size(); ++m)
1200  {
1201  for (auto & it : _sibling_elems[_geometric_cuts[m]->getInterfaceID()])
1202  {
1203  if (parent_elem == it.first)
1204  it.first = libmesh_elem;
1205  else if (parent_elem == it.second)
1206  it.second = libmesh_elem;
1207  }
1208  }
1209 
1210  // parent has at least two children
1211  if (new_elements[i]->getParent()->numChildren() > 1)
1212  temporary_parent_children_map[parent_elem->id()].push_back(libmesh_elem);
1213 
1214  Elem * parent_elem2 = nullptr;
1215  Elem * libmesh_elem2 = nullptr;
1216  if (_displaced_mesh)
1217  {
1218  parent_elem2 = _displaced_mesh->elem_ptr(parent_id);
1219  libmesh_elem2 = Elem::build(parent_elem2->type()).release();
1220 
1221  for (unsigned int m = 0; m < _geometric_cuts.size(); ++m)
1222  {
1223  for (auto & it : _sibling_displaced_elems[_geometric_cuts[m]->getInterfaceID()])
1224  {
1225  if (parent_elem2 == it.first)
1226  it.first = libmesh_elem2;
1227  else if (parent_elem2 == it.second)
1228  it.second = libmesh_elem2;
1229  }
1230  }
1231  }
1232 
1233  for (unsigned int j = 0; j < new_elements[i]->numNodes(); ++j)
1234  {
1235  unsigned int node_id = new_elements[i]->getNode(j)->id();
1236  Node * libmesh_node;
1237 
1238  std::map<unsigned int, Node *>::iterator nit = efa_id_to_new_node.find(node_id);
1239  if (nit != efa_id_to_new_node.end())
1240  libmesh_node = nit->second;
1241  else
1242  libmesh_node = _mesh->node_ptr(node_id);
1243 
1244  if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
1245  libmesh_node->processor_id() = parent_elem->processor_id();
1246 
1247  libmesh_elem->set_node(j) = libmesh_node;
1248 
1249  // Store solution for all nodes affected by XFEM (even existing nodes)
1250  if (parent_elem->is_semilocal(_mesh->processor_id()))
1251  {
1252  Node * solution_node = libmesh_node; // Node from which to store solution
1253  if (new_nodes_to_parents.find(libmesh_node) != new_nodes_to_parents.end())
1254  solution_node = new_nodes_to_parents[libmesh_node];
1255 
1256  if ((_moose_mesh->isSemiLocal(solution_node)) ||
1257  (libmesh_node->processor_id() == _mesh->processor_id()))
1258  {
1259  storeSolutionForNode(libmesh_node,
1260  solution_node,
1261  *nls[0],
1263  current_solution,
1264  old_solution,
1265  older_solution);
1266  storeSolutionForNode(libmesh_node,
1267  solution_node,
1268  aux,
1270  current_aux_solution,
1271  old_aux_solution,
1272  older_aux_solution);
1273  }
1274  }
1275 
1276  Node * parent_node = parent_elem->node_ptr(j);
1277  _mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
1278  _mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
1279 
1280  if (_displaced_mesh)
1281  {
1282  std::map<unsigned int, Node *>::iterator nit2 = efa_id_to_new_node2.find(node_id);
1283  if (nit2 != efa_id_to_new_node2.end())
1284  libmesh_node = nit2->second;
1285  else
1286  libmesh_node = _displaced_mesh->node_ptr(node_id);
1287 
1288  if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
1289  libmesh_node->processor_id() = parent_elem2->processor_id();
1290 
1291  libmesh_elem2->set_node(j) = libmesh_node;
1292 
1293  parent_node = parent_elem2->node_ptr(j);
1294  _displaced_mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
1295  _displaced_mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
1296  }
1297  }
1298 
1299  libmesh_elem->set_p_level(parent_elem->p_level());
1300  libmesh_elem->set_p_refinement_flag(parent_elem->p_refinement_flag());
1301  _mesh->add_elem(libmesh_elem);
1302  libmesh_elem->set_n_systems(parent_elem->n_systems());
1303  libmesh_elem->subdomain_id() = parent_elem->subdomain_id();
1304  libmesh_elem->processor_id() = parent_elem->processor_id();
1305 
1306  // The crack tip origin map is stored before cut, thus the elem should be updated with new
1307  // element.
1308  std::map<const Elem *, std::vector<Point>>::iterator mit =
1309  _elem_crack_origin_direction_map.find(parent_elem);
1310  if (mit != _elem_crack_origin_direction_map.end())
1311  {
1312  std::vector<Point> crack_data = _elem_crack_origin_direction_map[parent_elem];
1314  _elem_crack_origin_direction_map[libmesh_elem] = crack_data;
1315  }
1316 
1317  if (_debug_output_level > 1)
1318  _console << "XFEM added new element: " << libmesh_elem->id() << std::endl;
1319 
1320  XFEMCutElem * xfce = nullptr;
1321  if (_mesh->mesh_dimension() == 2)
1322  {
1323  EFAElement2D * new_efa_elem2d = dynamic_cast<EFAElement2D *>(new_elements[i]);
1324  if (!new_efa_elem2d)
1325  mooseError("EFAelem is not of EFAelement2D type");
1326  xfce = new XFEMCutElem2D(libmesh_elem,
1327  new_efa_elem2d,
1328  _fe_problem->assembly(0, /*nl_sys_num=*/0).qRule()->n_points(),
1329  libmesh_elem->n_sides());
1330  }
1331  else if (_mesh->mesh_dimension() == 3)
1332  {
1333  EFAElement3D * new_efa_elem3d = dynamic_cast<EFAElement3D *>(new_elements[i]);
1334  if (!new_efa_elem3d)
1335  mooseError("EFAelem is not of EFAelement3D type");
1336  xfce = new XFEMCutElem3D(libmesh_elem,
1337  new_efa_elem3d,
1338  _fe_problem->assembly(0, /*nl_sys_num=*/0).qRule()->n_points(),
1339  libmesh_elem->n_sides());
1340  }
1341  _cut_elem_map.insert(std::pair<unique_id_type, XFEMCutElem *>(libmesh_elem->unique_id(), xfce));
1342  efa_id_to_new_elem.insert(std::make_pair(efa_child_id, libmesh_elem));
1343 
1344  if (_displaced_mesh)
1345  {
1346  libmesh_elem2->set_p_level(parent_elem2->p_level());
1347  libmesh_elem2->set_p_refinement_flag(parent_elem2->p_refinement_flag());
1348  _displaced_mesh->add_elem(libmesh_elem2);
1349  libmesh_elem2->set_n_systems(parent_elem2->n_systems());
1350  libmesh_elem2->subdomain_id() = parent_elem2->subdomain_id();
1351  libmesh_elem2->processor_id() = parent_elem2->processor_id();
1352  }
1353 
1354  unsigned int n_sides = parent_elem->n_sides();
1355  for (unsigned int side = 0; side < n_sides; ++side)
1356  {
1357  _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
1358  _mesh->get_boundary_info().add_side(libmesh_elem, side, parent_boundary_ids);
1359  }
1360  if (_displaced_mesh)
1361  {
1362  n_sides = parent_elem2->n_sides();
1363  for (unsigned int side = 0; side < n_sides; ++side)
1364  {
1365  _displaced_mesh->get_boundary_info().boundary_ids(parent_elem2, side, parent_boundary_ids);
1366  _displaced_mesh->get_boundary_info().add_side(libmesh_elem2, side, parent_boundary_ids);
1367  }
1368  }
1369 
1370  unsigned int n_edges = parent_elem->n_edges();
1371  for (unsigned int edge = 0; edge < n_edges; ++edge)
1372  {
1373  _mesh->get_boundary_info().edge_boundary_ids(parent_elem, edge, parent_boundary_ids);
1374  _mesh->get_boundary_info().add_edge(libmesh_elem, edge, parent_boundary_ids);
1375  }
1376  if (_displaced_mesh)
1377  {
1378  n_edges = parent_elem2->n_edges();
1379  for (unsigned int edge = 0; edge < n_edges; ++edge)
1380  {
1381  _displaced_mesh->get_boundary_info().edge_boundary_ids(
1382  parent_elem2, edge, parent_boundary_ids);
1383  _displaced_mesh->get_boundary_info().add_edge(libmesh_elem2, edge, parent_boundary_ids);
1384  }
1385  }
1386 
1387  // TODO: Also need to copy neighbor material data
1388  if (parent_elem->processor_id() == _mesh->processor_id())
1389  {
1390  if (_material_data[0]->getMaterialPropertyStorage().hasStatefulProperties())
1391  _material_data[0]->copy(*libmesh_elem, *parent_elem, 0);
1392 
1393  if (_bnd_material_data[0]->getMaterialPropertyStorage().hasStatefulProperties())
1394  for (unsigned int side = 0; side < parent_elem->n_sides(); ++side)
1395  {
1396  _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
1397  std::vector<boundary_id_type>::iterator it_bd = parent_boundary_ids.begin();
1398  for (; it_bd != parent_boundary_ids.end(); ++it_bd)
1399  {
1400  if (_fe_problem->needBoundaryMaterialOnSide(*it_bd, 0))
1401  _bnd_material_data[0]->copy(*libmesh_elem, *parent_elem, side);
1402  }
1403  }
1404 
1405  // Store the current information about the geometrically cut element, and load cached material
1406  // properties into the new child element, if any.
1407  const GeometricCutUserObject * gcuo = getGeometricCutForElem(parent_elem);
1408  if (gcuo && gcuo->shouldHealMesh())
1409  {
1410  CutSubdomainID gcsid = getCutSubdomainID(gcuo, libmesh_elem, parent_elem);
1411  Xfem::CutElemInfo cei(parent_elem, gcuo, gcsid);
1412  _geom_cut_elems.emplace(libmesh_elem, cei);
1413  // Find the element to copy data from.
1414  // Iterate through the old geometrically cut elements, if its parent element AND the
1415  // geometric cut user object AND the cut subdomain ID are the same as the
1416  // current element, then that must be it.
1417  for (auto old_cei : _old_geom_cut_elems)
1418  if (cei.match(old_cei.second))
1419  {
1420  loadMaterialPropertiesForElement(libmesh_elem, old_cei.first, _old_geom_cut_elems);
1421  if (_debug_output_level > 1)
1422  _console << "XFEM set material properties for element: " << libmesh_elem->id()
1423  << "\n";
1424  break;
1425  }
1426  }
1427 
1428  // Store solution for all elements affected by XFEM
1429  storeSolutionForElement(libmesh_elem,
1430  parent_elem,
1431  *nls[0],
1433  current_solution,
1434  old_solution,
1435  older_solution);
1436  storeSolutionForElement(libmesh_elem,
1437  parent_elem,
1438  aux,
1440  current_aux_solution,
1441  old_aux_solution,
1442  older_aux_solution);
1443  }
1444  }
1445 
1446  // delete elements
1447  for (std::size_t i = 0; i < delete_elements.size(); ++i)
1448  {
1449  Elem * elem_to_delete = _mesh->elem_ptr(delete_elements[i]->id());
1450 
1451  // delete the XFEMCutElem object for any elements that are to be deleted
1452  std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
1453  _cut_elem_map.find(elem_to_delete->unique_id());
1454  if (cemit != _cut_elem_map.end())
1455  {
1456  delete cemit->second;
1457  _cut_elem_map.erase(cemit);
1458  }
1459 
1460  // remove the property storage of deleted element/side
1461  _material_data[0]->eraseProperty(elem_to_delete);
1462  _bnd_material_data[0]->eraseProperty(elem_to_delete);
1463 
1464  elem_to_delete->nullify_neighbors();
1465  _mesh->get_boundary_info().remove(elem_to_delete);
1466  unsigned int deleted_elem_id = elem_to_delete->id();
1467  _mesh->delete_elem(elem_to_delete);
1468  if (_debug_output_level > 1)
1469  _console << "XFEM deleted element: " << deleted_elem_id << std::endl;
1470 
1471  if (_displaced_mesh)
1472  {
1473  Elem * elem_to_delete2 = _displaced_mesh->elem_ptr(delete_elements[i]->id());
1474  elem_to_delete2->nullify_neighbors();
1475  _displaced_mesh->get_boundary_info().remove(elem_to_delete2);
1476  _displaced_mesh->delete_elem(elem_to_delete2);
1477  }
1478  }
1479 
1480  for (std::map<unsigned int, std::vector<const Elem *>>::iterator it =
1481  temporary_parent_children_map.begin();
1482  it != temporary_parent_children_map.end();
1483  ++it)
1484  {
1485  std::vector<const Elem *> & sibling_elem_vec = it->second;
1486  // TODO: for cut-node case, how to find the sibling elements?
1487  // if (sibling_elem_vec.size() != 2)
1488  // mooseError("Must have exactly 2 sibling elements");
1489 
1490  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
1491  for (auto const & elem_id : _geom_marker_id_elems[_geometric_cuts[i]->getInterfaceID()])
1492  if (it->first == elem_id)
1493  _sibling_elems[_geometric_cuts[i]->getInterfaceID()].push_back(
1494  std::make_pair(sibling_elem_vec[0], sibling_elem_vec[1]));
1495  }
1496 
1497  // add sibling elems on displaced mesh
1498  if (_displaced_mesh)
1499  {
1500  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
1501  {
1502  for (auto & se : _sibling_elems[_geometric_cuts[i]->getInterfaceID()])
1503  {
1504  Elem * elem = _displaced_mesh->elem_ptr(se.first->id());
1505  Elem * elem_pair = _displaced_mesh->elem_ptr(se.second->id());
1506  _sibling_displaced_elems[_geometric_cuts[i]->getInterfaceID()].push_back(
1507  std::make_pair(elem, elem_pair));
1508  }
1509  }
1510  }
1511 
1512  // clear the temporary map
1513  temporary_parent_children_map.clear();
1514 
1515  // Store information about crack tip elements
1516  if (mesh_changed)
1517  {
1518  _crack_tip_elems.clear();
1520  const std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
1521  std::set<EFAElement *>::const_iterator sit;
1522  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
1523  {
1524  unsigned int eid = (*sit)->id();
1525  Elem * crack_tip_elem;
1526  std::map<unsigned int, Elem *>::iterator eit = efa_id_to_new_elem.find(eid);
1527  if (eit != efa_id_to_new_elem.end())
1528  crack_tip_elem = eit->second;
1529  else
1530  crack_tip_elem = _mesh->elem_ptr(eid);
1531  _crack_tip_elems.insert(crack_tip_elem);
1532 
1533  // Store the crack tip elements which are going to be healed
1534  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
1535  {
1536  if (_geometric_cuts[i]->shouldHealMesh())
1537  {
1538  for (auto const & mie : _geom_marker_id_elems[_geometric_cuts[i]->getInterfaceID()])
1539  if ((*sit)->getParent() != nullptr)
1540  {
1541  if (_mesh->mesh_dimension() == 2)
1542  {
1543  EFAElement2D * efa_elem2d = dynamic_cast<EFAElement2D *>((*sit)->getParent());
1544  if (!efa_elem2d)
1545  mooseError("EFAelem is not of EFAelement2D type");
1546 
1547  for (unsigned int edge_id = 0; edge_id < efa_elem2d->numEdges(); ++edge_id)
1548  {
1549  for (unsigned int en_iter = 0; en_iter < efa_elem2d->numEdgeNeighbors(edge_id);
1550  ++en_iter)
1551  {
1552  EFAElement2D * edge_neighbor = efa_elem2d->getEdgeNeighbor(edge_id, en_iter);
1553  if (edge_neighbor != nullptr && edge_neighbor->id() == mie)
1554  _crack_tip_elems_to_be_healed.insert(crack_tip_elem);
1555  }
1556  }
1557  }
1558  else if (_mesh->mesh_dimension() == 3)
1559  {
1560  EFAElement3D * efa_elem3d = dynamic_cast<EFAElement3D *>((*sit)->getParent());
1561  if (!efa_elem3d)
1562  mooseError("EFAelem is not of EFAelement3D type");
1563 
1564  for (unsigned int face_id = 0; face_id < efa_elem3d->numFaces(); ++face_id)
1565  {
1566  for (unsigned int fn_iter = 0; fn_iter < efa_elem3d->numFaceNeighbors(face_id);
1567  ++fn_iter)
1568  {
1569  EFAElement3D * face_neighbor = efa_elem3d->getFaceNeighbor(face_id, fn_iter);
1570  if (face_neighbor != nullptr && face_neighbor->id() == mie)
1571  _crack_tip_elems_to_be_healed.insert(crack_tip_elem);
1572  }
1573  }
1574  }
1575  }
1576  }
1577  }
1578  }
1579  }
1580 
1581  if (_debug_output_level > 0)
1582  {
1583  _console << "\nXFEM mesh cutting with element fragment algorithm complete\n";
1584  _console << "# new nodes: " << new_nodes.size() << "\n";
1585  _console << "# new elements: " << new_elements.size() << "\n";
1586  _console << "# deleted elements: " << delete_elements.size() << "\n";
1587  _console << std::flush;
1588  }
1589 
1590  // store virtual nodes
1591  // store cut edge info
1592  return mesh_changed;
1593 }
1594 
1595 Point
1597  EFAElement * CEMElem,
1598  const Elem * elem,
1599  MeshBase * displaced_mesh) const
1600 {
1601  Point node_coor(0.0, 0.0, 0.0);
1602  std::vector<EFANode *> master_nodes;
1603  std::vector<Point> master_points;
1604  std::vector<double> master_weights;
1605 
1606  CEMElem->getMasterInfo(CEMnode, master_nodes, master_weights);
1607  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1608  {
1609  if (master_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT)
1610  {
1611  unsigned int local_node_id = CEMElem->getLocalNodeIndex(master_nodes[i]);
1612  const Node * node = elem->node_ptr(local_node_id);
1613  if (displaced_mesh)
1614  node = displaced_mesh->node_ptr(node->id());
1615  Point node_p((*node)(0), (*node)(1), (*node)(2));
1616  master_points.push_back(node_p);
1617  }
1618  else
1619  mooseError("master nodes must be permanent");
1620  }
1621  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1622  node_coor += master_weights[i] * master_points[i];
1623 
1624  return node_coor;
1625 }
1626 
1627 Real
1628 XFEM::getPhysicalVolumeFraction(const Elem * elem) const
1629 {
1630  Real phys_volfrac = 1.0;
1631  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1632  it = _cut_elem_map.find(elem->unique_id());
1633  if (it != _cut_elem_map.end())
1634  {
1635  XFEMCutElem * xfce = it->second;
1636  const EFAElement * EFAelem = xfce->getEFAElement();
1637  if (EFAelem->isPartial())
1638  { // exclude the full crack tip elements
1640  phys_volfrac = xfce->getPhysicalVolumeFraction();
1641  }
1642  }
1643 
1644  return phys_volfrac;
1645 }
1646 
1647 bool
1648 XFEM::isPointInsidePhysicalDomain(const Elem * elem, const Point & point) const
1649 {
1650  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1651  it = _cut_elem_map.find(elem->unique_id());
1652  if (it != _cut_elem_map.end())
1653  {
1654  XFEMCutElem * xfce = it->second;
1655 
1656  if (xfce->isPointPhysical(point))
1657  return true;
1658  }
1659  else
1660  return true;
1661 
1662  return false;
1663 }
1664 
1665 Real
1666 XFEM::getCutPlane(const Elem * elem,
1667  const Xfem::XFEM_CUTPLANE_QUANTITY quantity,
1668  unsigned int plane_id) const
1669 {
1670  Real comp = 0.0;
1671  Point planedata(0.0, 0.0, 0.0);
1672  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1673  it = _cut_elem_map.find(elem->unique_id());
1674  if (it != _cut_elem_map.end())
1675  {
1676  const XFEMCutElem * xfce = it->second;
1677  const EFAElement * EFAelem = xfce->getEFAElement();
1678  if (EFAelem->isPartial()) // exclude the full crack tip elements
1679  {
1680  if ((unsigned int)quantity < 3)
1681  {
1682  unsigned int index = (unsigned int)quantity;
1683  planedata = xfce->getCutPlaneOrigin(plane_id, _displaced_mesh);
1684  comp = planedata(index);
1685  }
1686  else if ((unsigned int)quantity < 6)
1687  {
1688  unsigned int index = (unsigned int)quantity - 3;
1689  planedata = xfce->getCutPlaneNormal(plane_id, _displaced_mesh);
1690  comp = planedata(index);
1691  }
1692  else
1693  mooseError("In get_cut_plane index out of range");
1694  }
1695  }
1696  return comp;
1697 }
1698 
1699 bool
1700 XFEM::isElemAtCrackTip(const Elem * elem) const
1701 {
1702  return (_crack_tip_elems.find(elem) != _crack_tip_elems.end());
1703 }
1704 
1705 bool
1706 XFEM::isElemCut(const Elem * elem, XFEMCutElem *& xfce) const
1707 {
1708  const auto it = _cut_elem_map.find(elem->unique_id());
1709  if (it != _cut_elem_map.end())
1710  {
1711  xfce = it->second;
1712  const EFAElement * EFAelem = xfce->getEFAElement();
1713  if (EFAelem->isPartial()) // exclude the full crack tip elements
1714  return true;
1715  }
1716 
1717  xfce = nullptr;
1718  return false;
1719 }
1720 
1721 bool
1722 XFEM::isElemCut(const Elem * elem) const
1723 {
1724  XFEMCutElem * xfce;
1725  return isElemCut(elem, xfce);
1726 }
1727 
1728 void
1729 XFEM::getFragmentFaces(const Elem * elem,
1730  std::vector<std::vector<Point>> & frag_faces,
1731  bool displaced_mesh) const
1732 {
1733  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1734  it = _cut_elem_map.find(elem->unique_id());
1735  if (it != _cut_elem_map.end())
1736  {
1737  const XFEMCutElem * xfce = it->second;
1738  if (displaced_mesh)
1739  xfce->getFragmentFaces(frag_faces, _displaced_mesh);
1740  else
1741  xfce->getFragmentFaces(frag_faces);
1742  }
1743 }
1744 
1745 EFAElement2D *
1746 XFEM::getEFAElem2D(const Elem * elem)
1747 {
1748  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
1749  EFAElement2D * EFAelem2D = dynamic_cast<EFAElement2D *>(EFAelem);
1750 
1751  if (!EFAelem2D)
1752  mooseError("EFAelem is not of EFAelement2D type");
1753 
1754  return EFAelem2D;
1755 }
1756 
1757 EFAElement3D *
1758 XFEM::getEFAElem3D(const Elem * elem)
1759 {
1760  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
1761  EFAElement3D * EFAelem3D = dynamic_cast<EFAElement3D *>(EFAelem);
1762 
1763  if (!EFAelem3D)
1764  mooseError("EFAelem is not of EFAelement3D type");
1765 
1766  return EFAelem3D;
1767 }
1768 
1769 void
1770 XFEM::getFragmentEdges(const Elem * elem,
1771  EFAElement2D * CEMElem,
1772  std::vector<std::vector<Point>> & frag_edges) const
1773 {
1774  // N.B. CEMElem here has global EFAnode
1775  frag_edges.clear();
1776  if (CEMElem->numFragments() > 0)
1777  {
1778  if (CEMElem->numFragments() > 1)
1779  mooseError("element ", elem->id(), " has more than one fragment at this point");
1780  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
1781  {
1782  std::vector<Point> p_line(2, Point(0.0, 0.0, 0.0));
1783  p_line[0] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
1784  p_line[1] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
1785  frag_edges.push_back(p_line);
1786  }
1787  }
1788 }
1789 
1790 void
1791 XFEM::getFragmentFaces(const Elem * elem,
1792  EFAElement3D * CEMElem,
1793  std::vector<std::vector<Point>> & frag_faces) const
1794 {
1795  // N.B. CEMElem here has global EFAnode
1796  frag_faces.clear();
1797  if (CEMElem->numFragments() > 0)
1798  {
1799  if (CEMElem->numFragments() > 1)
1800  mooseError("element ", elem->id(), " has more than one fragment at this point");
1801  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numFaces(); ++i)
1802  {
1803  unsigned int num_face_nodes = CEMElem->getFragmentFace(0, i)->numNodes();
1804  std::vector<Point> p_line(num_face_nodes, Point(0.0, 0.0, 0.0));
1805  for (unsigned int j = 0; j < num_face_nodes; ++j)
1806  p_line[j] = getEFANodeCoords(CEMElem->getFragmentFace(0, i)->getNode(j), CEMElem, elem);
1807  frag_faces.push_back(p_line);
1808  }
1809  }
1810 }
1811 
1814 {
1815  return _XFEM_qrule;
1816 }
1817 
1818 void
1819 XFEM::setXFEMQRule(std::string & xfem_qrule)
1820 {
1821  if (xfem_qrule == "volfrac")
1823  else if (xfem_qrule == "moment_fitting")
1825  else if (xfem_qrule == "direct")
1827 }
1828 
1829 void
1830 XFEM::setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment)
1831 {
1832  _use_crack_growth_increment = use_crack_growth_increment;
1833  _crack_growth_increment = crack_growth_increment;
1834 }
1835 
1836 void
1837 XFEM::setDebugOutputLevel(unsigned int debug_output_level)
1838 {
1839  _debug_output_level = debug_output_level;
1840 }
1841 
1842 bool
1844  const Elem * elem,
1845  QBase * qrule,
1846  const MooseArray<Point> & q_points)
1847 {
1848  bool have_weights = false;
1849  XFEMCutElem * xfce = nullptr;
1850  if (isElemCut(elem, xfce))
1851  {
1852  mooseAssert(xfce != nullptr, "Must have valid XFEMCutElem object here");
1853  xfce->getWeightMultipliers(weights, qrule, getXFEMQRule(), q_points);
1854  have_weights = true;
1855  }
1856  return have_weights;
1857 }
1858 
1859 bool
1861  const Elem * elem,
1862  QBase * qrule,
1863  const MooseArray<Point> & q_points,
1864  unsigned int side)
1865 {
1866  bool have_weights = false;
1867  XFEMCutElem * xfce = nullptr;
1868  if (isElemCut(elem, xfce))
1869  {
1870  mooseAssert(xfce != nullptr, "Must have valid XFEMCutElem object here");
1871  xfce->getFaceWeightMultipliers(weights, qrule, getXFEMQRule(), q_points, side);
1872  have_weights = true;
1873  }
1874  return have_weights;
1875 }
1876 
1877 void
1879  unsigned int plane_id,
1880  Point & normal,
1881  std::vector<Point> & intersectionPoints,
1882  bool displaced_mesh) const
1883 {
1884  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1885  it = _cut_elem_map.find(elem->unique_id());
1886  if (it != _cut_elem_map.end())
1887  {
1888  const XFEMCutElem * xfce = it->second;
1889  if (displaced_mesh)
1890  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints, _displaced_mesh);
1891  else
1892  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints);
1893  }
1894 }
1895 
1896 void
1897 XFEM::getXFEMqRuleOnLine(std::vector<Point> & intersection_points,
1898  std::vector<Point> & quad_pts,
1899  std::vector<Real> & quad_wts) const
1900 {
1901  Point p1 = intersection_points[0];
1902  Point p2 = intersection_points[1];
1903 
1904  // number of quadrature points
1905  std::size_t num_qpoints = 2;
1906 
1907  // quadrature coordinates
1908  Real xi0 = -std::sqrt(1.0 / 3.0);
1909  Real xi1 = std::sqrt(1.0 / 3.0);
1910 
1911  quad_wts.resize(num_qpoints);
1912  quad_pts.resize(num_qpoints);
1913 
1914  Real integ_jacobian = pow((p1 - p2).norm_sq(), 0.5) * 0.5;
1915 
1916  quad_wts[0] = 1.0 * integ_jacobian;
1917  quad_wts[1] = 1.0 * integ_jacobian;
1918 
1919  quad_pts[0] = (1.0 - xi0) / 2.0 * p1 + (1.0 + xi0) / 2.0 * p2;
1920  quad_pts[1] = (1.0 - xi1) / 2.0 * p1 + (1.0 + xi1) / 2.0 * p2;
1921 }
1922 
1923 void
1924 XFEM::getXFEMqRuleOnSurface(std::vector<Point> & intersection_points,
1925  std::vector<Point> & quad_pts,
1926  std::vector<Real> & quad_wts) const
1927 {
1928  std::size_t nnd_pe = intersection_points.size();
1929  Point xcrd(0.0, 0.0, 0.0);
1930  for (std::size_t i = 0; i < nnd_pe; ++i)
1931  xcrd += intersection_points[i];
1932  xcrd /= nnd_pe;
1933 
1934  quad_pts.resize(nnd_pe);
1935  quad_wts.resize(nnd_pe);
1936 
1937  Real jac = 0.0;
1938 
1939  for (std::size_t j = 0; j < nnd_pe; ++j) // loop all sub-tris
1940  {
1941  std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
1942  std::vector<Point> subtrig_points(3, Point(0.0, 0.0, 0.0)); // sub-trig nodal coords
1943 
1944  int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
1945  subtrig_points[0] = xcrd;
1946  subtrig_points[1] = intersection_points[j];
1947  subtrig_points[2] = intersection_points[jplus1];
1948 
1949  std::vector<std::vector<Real>> sg2;
1950  Xfem::stdQuadr2D(3, 1, sg2); // get sg2
1951  for (std::size_t l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-trig
1952  {
1953  Xfem::shapeFunc2D(3, sg2[l], subtrig_points, shape, jac, true); // Get shape
1954  std::vector<Real> tsg_line(3, 0.0);
1955  for (std::size_t k = 0; k < 3; ++k) // loop sub-trig nodes
1956  {
1957  tsg_line[0] += shape[k][2] * subtrig_points[k](0);
1958  tsg_line[1] += shape[k][2] * subtrig_points[k](1);
1959  tsg_line[2] += shape[k][2] * subtrig_points[k](2);
1960  }
1961  quad_pts[j + l] = Point(tsg_line[0], tsg_line[1], tsg_line[2]);
1962  quad_wts[j + l] = sg2[l][3] * jac;
1963  }
1964  }
1965 }
1966 
1967 void
1968 XFEM::storeSolutionForNode(const Node * node_to_store_to,
1969  const Node * node_to_store_from,
1970  SystemBase & sys,
1971  std::map<unique_id_type, std::vector<Real>> & stored_solution,
1972  const NumericVector<Number> & current_solution,
1973  const NumericVector<Number> & old_solution,
1974  const NumericVector<Number> & older_solution)
1975 {
1976  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node_to_store_from, sys);
1977  std::vector<Real> stored_solution_scratch;
1978  // Size for current solution, as well as for old, and older solution only for transient case
1979  std::size_t stored_solution_size =
1980  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
1981  stored_solution_scratch.reserve(stored_solution_size);
1982 
1983  // Store in the order defined in stored_solution_dofs first for the current, then for old and
1984  // older if applicable
1985  for (auto dof : stored_solution_dofs)
1986  stored_solution_scratch.push_back(current_solution(dof));
1987 
1988  if (_fe_problem->isTransient())
1989  {
1990  for (auto dof : stored_solution_dofs)
1991  stored_solution_scratch.push_back(old_solution(dof));
1992 
1993  for (auto dof : stored_solution_dofs)
1994  stored_solution_scratch.push_back(older_solution(dof));
1995  }
1996 
1997  if (stored_solution_scratch.size() > 0)
1998  stored_solution[node_to_store_to->unique_id()] = stored_solution_scratch;
1999 }
2000 
2001 void
2002 XFEM::storeSolutionForElement(const Elem * elem_to_store_to,
2003  const Elem * elem_to_store_from,
2004  SystemBase & sys,
2005  std::map<unique_id_type, std::vector<Real>> & stored_solution,
2006  const NumericVector<Number> & current_solution,
2007  const NumericVector<Number> & old_solution,
2008  const NumericVector<Number> & older_solution)
2009 {
2010  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem_to_store_from, sys);
2011  std::vector<Real> stored_solution_scratch;
2012  // Size for current solution, as well as for old, and older solution only for transient case
2013  std::size_t stored_solution_size =
2014  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
2015  stored_solution_scratch.reserve(stored_solution_size);
2016 
2017  // Store in the order defined in stored_solution_dofs first for the current, then for old and
2018  // older if applicable
2019  for (auto dof : stored_solution_dofs)
2020  stored_solution_scratch.push_back(current_solution(dof));
2021 
2022  if (_fe_problem->isTransient())
2023  {
2024  for (auto dof : stored_solution_dofs)
2025  stored_solution_scratch.push_back(old_solution(dof));
2026 
2027  for (auto dof : stored_solution_dofs)
2028  stored_solution_scratch.push_back(older_solution(dof));
2029  }
2030 
2031  if (stored_solution_scratch.size() > 0)
2032  stored_solution[elem_to_store_to->unique_id()] = stored_solution_scratch;
2033 }
2034 
2035 void
2037  const std::map<unique_id_type, std::vector<Real>> & stored_solution,
2038  NumericVector<Number> & current_solution,
2039  NumericVector<Number> & old_solution,
2040  NumericVector<Number> & older_solution)
2041 {
2042  for (auto & node : _mesh->local_node_ptr_range())
2043  {
2044  auto mit = stored_solution.find(node->unique_id());
2045  if (mit != stored_solution.end())
2046  {
2047  const std::vector<Real> & stored_node_solution = mit->second;
2048  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node, sys);
2049  setSolutionForDOFs(stored_node_solution,
2050  stored_solution_dofs,
2051  current_solution,
2052  old_solution,
2053  older_solution);
2054  }
2055  }
2056 
2057  for (auto & elem : as_range(_mesh->local_elements_begin(), _mesh->local_elements_end()))
2058  {
2059  auto mit = stored_solution.find(elem->unique_id());
2060  if (mit != stored_solution.end())
2061  {
2062  const std::vector<Real> & stored_elem_solution = mit->second;
2063  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem, sys);
2064  setSolutionForDOFs(stored_elem_solution,
2065  stored_solution_dofs,
2066  current_solution,
2067  old_solution,
2068  older_solution);
2069  }
2070  }
2071 }
2072 
2073 void
2074 XFEM::setSolutionForDOFs(const std::vector<Real> & stored_solution,
2075  const std::vector<dof_id_type> & stored_solution_dofs,
2076  NumericVector<Number> & current_solution,
2077  NumericVector<Number> & old_solution,
2078  NumericVector<Number> & older_solution)
2079 {
2080  // Solution vector is stored first for current, then old and older solutions.
2081  // These are the offsets to the beginning of the old and older solutions in the vector.
2082  const auto old_solution_offset = stored_solution_dofs.size();
2083  const auto older_solution_offset = old_solution_offset * 2;
2084 
2085  for (std::size_t i = 0; i < stored_solution_dofs.size(); ++i)
2086  {
2087  current_solution.set(stored_solution_dofs[i], stored_solution[i]);
2088  if (_fe_problem->isTransient())
2089  {
2090  old_solution.set(stored_solution_dofs[i], stored_solution[old_solution_offset + i]);
2091  older_solution.set(stored_solution_dofs[i], stored_solution[older_solution_offset + i]);
2092  }
2093  }
2094 }
2095 
2096 std::vector<dof_id_type>
2097 XFEM::getElementSolutionDofs(const Elem * elem, SystemBase & sys) const
2098 {
2099  SubdomainID sid = elem->subdomain_id();
2100  const std::vector<MooseVariableFEBase *> & vars = sys.getVariables(0);
2101  std::vector<dof_id_type> solution_dofs;
2102  solution_dofs.reserve(vars.size()); // just an approximation
2103  for (auto var : vars)
2104  {
2105  if (!var->isNodal())
2106  {
2107  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
2108  if (var_subdomains.empty() || var_subdomains.find(sid) != var_subdomains.end())
2109  {
2110  unsigned int n_comp = elem->n_comp(sys.number(), var->number());
2111  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
2112  {
2113  dof_id_type elem_dof = elem->dof_number(sys.number(), var->number(), icomp);
2114  solution_dofs.push_back(elem_dof);
2115  }
2116  }
2117  }
2118  }
2119  return solution_dofs;
2120 }
2121 
2122 std::vector<dof_id_type>
2123 XFEM::getNodeSolutionDofs(const Node * node, SystemBase & sys) const
2124 {
2125  const std::set<SubdomainID> & sids = _moose_mesh->getNodeBlockIds(*node);
2126  const std::vector<MooseVariableFEBase *> & vars = sys.getVariables(0);
2127  std::vector<dof_id_type> solution_dofs;
2128  solution_dofs.reserve(vars.size()); // just an approximation
2129  for (auto var : vars)
2130  {
2131  if (var->isNodal())
2132  {
2133  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
2134  std::set<SubdomainID> intersect;
2135  set_intersection(var_subdomains.begin(),
2136  var_subdomains.end(),
2137  sids.begin(),
2138  sids.end(),
2139  std::inserter(intersect, intersect.begin()));
2140  if (var_subdomains.empty() || !intersect.empty())
2141  {
2142  unsigned int n_comp = node->n_comp(sys.number(), var->number());
2143  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
2144  {
2145  dof_id_type node_dof = node->dof_number(sys.number(), var->number(), icomp);
2146  solution_dofs.push_back(node_dof);
2147  }
2148  }
2149  }
2150  }
2151  return solution_dofs;
2152 }
2153 
2154 const GeometricCutUserObject *
2155 XFEM::getGeometricCutForElem(const Elem * elem) const
2156 {
2157  for (auto gcmit : _geom_marker_id_elems)
2158  {
2159  std::set<unsigned int> elems = gcmit.second;
2160  if (elems.find(elem->id()) != elems.end())
2161  return _geometric_cuts[gcmit.first];
2162  }
2163  return nullptr;
2164 }
2165 
2166 void
2168 {
2169  for (const auto state : storage.statefulIndexRange())
2170  {
2171  const auto & elem_props = storage.props(state).at(elem);
2172  auto & serialized_props = _geom_cut_elems[elem]._elem_material_properties[state - 1];
2173  serialized_props.clear();
2174  for (const auto & side_props_pair : elem_props)
2175  {
2176  const auto side = side_props_pair.first;
2177  std::ostringstream oss;
2178  dataStore(oss, storage.setProps(elem, side, state), nullptr);
2179  serialized_props[side].assign(oss.str());
2180  }
2181  }
2182 }
2183 
2184 void
2185 XFEM::storeMaterialPropertiesForElement(const Elem * parent_elem, const Elem * child_elem)
2186 {
2187  // Set the parent element so that it is consistent post-healing
2188  _geom_cut_elems[child_elem]._parent_elem = parent_elem;
2189 
2190  // Locally store the element material properties
2192  _material_data[0]->getMaterialPropertyStorageForXFEM({}));
2193 
2194  // Locally store the boundary material properties
2195  // First check if any of the side need material properties
2196  bool need_boundary_materials = false;
2197  for (unsigned int side = 0; side < child_elem->n_sides(); ++side)
2198  {
2199  std::vector<boundary_id_type> elem_boundary_ids;
2200  _mesh->get_boundary_info().boundary_ids(child_elem, side, elem_boundary_ids);
2201  for (auto bdid : elem_boundary_ids)
2203  need_boundary_materials = true;
2204  }
2205 
2206  // If boundary material properties are needed for this element, then store them.
2207  if (need_boundary_materials)
2209  child_elem, _bnd_material_data[0]->getMaterialPropertyStorageForXFEM({}));
2210 }
2211 
2212 void
2214  const Xfem::CachedMaterialProperties & cached_props,
2215  MaterialPropertyStorage & storage) const
2216 {
2217  if (!storage.hasStatefulProperties())
2218  return;
2219 
2220  for (const auto state : storage.statefulIndexRange())
2221  {
2222  const auto & serialized_props = cached_props[state - 1];
2223  for (const auto & [side, serialized_side_props] : serialized_props)
2224  {
2225  std::istringstream iss;
2226  iss.str(serialized_side_props);
2227  iss.clear();
2228 
2229  // This is very dirty. We should not write to MOOSE's stateful properties.
2230  // Please remove me :(
2231  dataLoad(iss, storage.setProps(elem, side, state), nullptr);
2232  }
2233  }
2234 }
2235 
2236 void
2238  const Elem * elem,
2239  const Elem * elem_from,
2240  std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei) const
2241 {
2242  // Restore the element material properties
2243  mooseAssert(cached_cei.count(elem_from) > 0, "XFEM: Unable to find cached material properties.");
2244  Xfem::CutElemInfo & cei = cached_cei[elem_from];
2245 
2246  // Load element material properties from cached properties
2248  cei._elem_material_properties,
2249  _material_data[0]->getMaterialPropertyStorageForXFEM({}));
2250 
2251  // Check if any of the element side need material properties
2252  bool need_boundary_materials = false;
2253  for (unsigned int side = 0; side < elem->n_sides(); ++side)
2254  {
2255  std::vector<boundary_id_type> elem_boundary_ids;
2256  _mesh->get_boundary_info().boundary_ids(elem, side, elem_boundary_ids);
2257  for (auto bdid : elem_boundary_ids)
2259  need_boundary_materials = true;
2260  }
2261 
2262  // Load boundary material properties from cached properties
2263  if (need_boundary_materials)
2265  elem,
2266  cei._bnd_material_properties,
2267  _bnd_material_data[0]->getMaterialPropertyStorageForXFEM({}));
2268 }
2269 
2272  const Elem * cut_elem,
2273  const Elem * parent_elem) const
2274 {
2275  if (!parent_elem)
2276  parent_elem = cut_elem;
2277  // Pick any node from the parent element that is inside the physical domain and return its
2278  // CutSubdomainID.
2279  const Node * node = pickFirstPhysicalNode(cut_elem, parent_elem);
2280  return gcuo->getCutSubdomainID(node);
2281 }
2282 
2283 const Node *
2284 XFEM::pickFirstPhysicalNode(const Elem * e, const Elem * e0) const
2285 {
2286  for (auto i : e0->node_index_range())
2287  if (isPointInsidePhysicalDomain(e, e0->node_ref(i)))
2288  return e0->node_ptr(i);
2289  mooseError("cannot find a physical node in the current element");
2290  return nullptr;
2291 }
~XFEM()
Definition: XFEM.C:46
void getCrackTipOrigin(std::map< unsigned int, const Elem *> &elem_id_crack_tip, std::vector< Point > &crack_front_points)
Definition: XFEM.C:65
EFAFragment3D * getFragment(unsigned int frag_id) const
void getFragmentFaces(const Elem *elem, std::vector< std::vector< Point >> &frag_faces, bool displaced_mesh=false) const
Definition: XFEM.C:1729
void correctCrackExtensionDirection(const Elem *elem, EFAElement2D *CEMElem, EFAEdge *orig_edge, Point normal, Point crack_tip_origin, Point crack_tip_direction, Real &distance_keep, unsigned int &edge_id_keep, Point &normal_keep)
Definition: XFEM.C:444
bool _use_crack_growth_increment
Definition: XFEM.h:324
unsigned int _debug_output_level
Controls amount of debugging output information 0: None 1: Summary 2: Details on modifications to mes...
Definition: XFEM.h:362
void setSolutionForDOFs(const std::vector< Real > &stored_solution, const std::vector< dof_id_type > &stored_solution_dofs, NumericVector< Number > &current_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
Set the solution for a set of DOFs.
Definition: XFEM.C:2074
MaterialProperties & setProps(const Elem *elem, unsigned int side, const unsigned int state=0)
const GeometricCutUserObject * getGeometricCutForElem(const Elem *elem) const
Get the GeometricCutUserObject associated with an element.
Definition: XFEM.C:2155
const std::vector< MooseVariableFieldBase *> & getVariables(THREAD_ID tid)
unsigned int numEdgeNeighbors(unsigned int edge_id) const
bool markCutEdgesByState(Real time)
Definition: XFEM.C:585
void addElemNodeIntersection(unsigned int elemid, unsigned int nodeid)
void storeMaterialPropertiesForElement(const Elem *parent_elem, const Elem *child_elem)
Helper function to store the material properties of a healed element.
Definition: XFEM.C:2185
std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo3D > > _geom_marked_elems_3d
Data structure for storing information about all 3D elements to be cut by geometry.
Definition: XFEM.h:347
EFAElement3D * getFaceNeighbor(unsigned int face_id, unsigned int neighbor_id) const
EFANode * getEmbeddedNode(unsigned int index) const
Definition: EFAEdge.C:332
CutSubdomainID getCutSubdomainID(const GeometricCutUserObject *gcuo, const Elem *cut_elem, const Elem *parent_elem=nullptr) const
Determine which cut subdomain the element belongs to relative to the cut.
Definition: XFEM.C:2271
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Definition: XFEMFuncs.C:269
Data structure describing geometrically described cut through 3D element.
XFEM_QRULE
Definition: XFEM.h:37
bool isFacePhantom(unsigned int face_id) const
void addElemEdgeIntersection(unsigned int elemid, unsigned int edgeid, double position)
std::vector< const GeometricCutUserObject * > _geometric_cuts
Definition: XFEM.h:327
void setSolution(SystemBase &sys, const std::map< unique_id_type, std::vector< Real >> &stored_solution, NumericVector< Number > &current_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
Set the solution for all locally-owned nodes/elements that have stored values.
Definition: XFEM.C:2036
unsigned int getCrackTipSplitElementID() const
Definition: EFAElement2D.C:599
unsigned int id() const
Definition: EFAElement.C:28
virtual bool update(Real time, const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override
Definition: XFEM.C:260
std::map< unsigned int, ElementPairLocator::ElementPairList > _sibling_elems
Definition: XFEM.h:332
Real _crack_growth_increment
Definition: XFEM.h:325
bool isSecondaryInteriorEdge(unsigned int edge_id) const
bool shouldHealMesh() const
Should the elements cut by this cutting object be healed in the current time step?
virtual bool getXFEMFaceWeights(MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points, unsigned int side) override
Definition: XFEM.C:1860
void mooseError(Args &&... args)
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=nullptr) const =0
bool isPartialOverlap(const EFAEdge &other) const
Definition: EFAEdge.C:59
bool isDistributedMesh() const
std::vector< dof_id_type > getNodeSolutionDofs(const Node *node, SystemBase &sys) const
Get a vector of the dof indices for all components of all variables associated with a node...
Definition: XFEM.C:2123
bool isPointInsidePhysicalDomain(const Elem *elem, const Point &point) const
Return true if the point is inside the element physical domain Note: if this element is not cut...
Definition: XFEM.C:1648
bool isSemiLocal(Node *const node) const
void getFragmentEdges(const Elem *elem, EFAElement2D *CEMElem, std::vector< std::vector< Point >> &frag_edges) const
Definition: XFEM.C:1770
bool match(const CutElemInfo &rhs)
Definition: XFEM.h:83
bool cutMeshWithEFA(const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux)
Definition: XFEM.C:1101
unsigned int getTipEdgeID() const
unsigned int numEdges() const
EFAEdge * getEdge(unsigned int edge_id) const
std::map< unsigned int, ElementPairLocator::ElementPairList > _sibling_displaced_elems
Definition: XFEM.h:333
Real getPhysicalVolumeFraction() const
Returns the volume fraction of the element fragment.
Definition: XFEMCutElem.C:51
void clearGeomMarkedElems()
Clear out the list of elements to be marked for cutting.
Definition: XFEM.C:177
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:329
EFAElement * add3DElement(const std::vector< unsigned int > &quad, unsigned int id)
bool needBoundaryMaterialOnSide(BoundaryID bnd_id, const THREAD_ID tid)
const std::vector< EFAElement * > & getChildElements()
std::set< const Elem * > _crack_tip_elems
Definition: XFEM.h:330
bool healMesh()
Potentially heal the mesh by merging some of the pairs of partial elements cut by XFEM back into sing...
Definition: XFEM.C:917
XFEM(const InputParameters &params)
Definition: XFEM.C:36
virtual bool isPartial() const =0
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:355
const std::set< SubdomainID > & getNodeBlockIds(const Node &node) const
bool hasIntersection() const
Definition: EFAEdge.C:200
virtual unsigned int numFragments() const
Definition: EFAElement3D.C:272
NumericVector< Number > & solutionOlder()
std::array< std::unordered_map< unsigned int, std::string >, 2 > CachedMaterialProperties
Convenient typedef for local storage of stateful material properties.
Definition: XFEM.h:50
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
IntRange< unsigned int > statefulIndexRange() const
std::vector< unsigned int > getInteriorEdgeID() const
unsigned int numEdges() const
void getFaceWeightMultipliers(MooseArray< Real > &face_weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points, unsigned int side)
Definition: XFEMCutElem.C:77
bool isEdgePhantom(unsigned int edge_id) const
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:207
Real distance(const Point &p)
std::set< const Elem * > _crack_tip_elems_to_be_healed
Definition: XFEM.h:331
bool _has_secondary_cut
Definition: XFEM.h:320
void addGeometricCut(GeometricCutUserObject *geometric_cut)
Definition: XFEM.C:55
std::map< unique_id_type, std::vector< Real > > _cached_aux_solution
Data structure to store the auxiliary solution for nodes/elements affected by XFEM For each node/elem...
Definition: XFEM.h:380
unsigned int numFaces() const
bool markCutEdgesByGeometry()
Definition: XFEM.C:392
unsigned int CutSubdomainID
Definition: XFEMAppTypes.h:18
FEProblemBase * _fe_problem
std::vector< MaterialData *> _bnd_material_data
EFANode * getNode(unsigned int node_id) const
Definition: EFAFace.C:99
void updateTopology(bool mergeUncutVirtualEdges=true)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const std::vector< EFAElement * > & getParentElements()
virtual const EFAElement * getEFAElement() const =0
bool initCutIntersectionEdge(Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
Definition: XFEM.C:898
void addGeomMarkedElem3D(const unsigned int elem_id, const Xfem::GeomMarkedElemInfo3D geom_info, const unsigned int interface_id)
Add information about a new cut to be performed on a specific 3d element.
Definition: XFEM.C:167
virtual void execute(const ExecFlagType &exec_type)
void addElemFaceIntersection(unsigned int elemid, unsigned int faceid, const std::vector< unsigned int > &edgeid, const std::vector< double > &position)
void storeCrackTipOriginAndDirection()
Definition: XFEM.C:184
const Node * pickFirstPhysicalNode(const Elem *e, const Elem *e0) const
Return the first node in the provided element that is found to be in the physical domain...
Definition: XFEM.C:2284
void loadMaterialPropertiesForElement(const Elem *elem, const Elem *elem_from, std::unordered_map< const Elem *, Xfem::CutElemInfo > &cached_cei) const
Helper function to store the material properties of a healed element.
Definition: XFEM.C:2237
const std::vector< EFANode * > & getNewNodes()
EFAFragment2D * getFragment(unsigned int frag_id) const
Real getCutPlane(const Elem *elem, const Xfem::XFEM_CUTPLANE_QUANTITY quantity, unsigned int plane_id) const
Get specified component of normal or origin for cut plane for a given element.
Definition: XFEM.C:1666
bool markCuts(Real time)
Definition: XFEM.C:375
unsigned int numNodes() const
Definition: EFAFace.C:87
std::map< const GeometricCutUserObject *, unsigned int > _geom_marker_id_map
Data structure for storing the GeommetricCutUserObjects and their corresponding id.
Definition: XFEM.h:353
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=nullptr) const =0
std::vector< MaterialData *> _material_data
const std::string name
Definition: Setup.h:20
void addGeomMarkedElem2D(const unsigned int elem_id, const Xfem::GeomMarkedElemInfo2D geom_info, const unsigned int interface_id)
Add information about a new cut to be performed on a specific 2d element.
Definition: XFEM.C:157
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
void storeSolutionForElement(const Elem *elem_to_store_to, const Elem *elem_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > &current_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
Store the solution in stored_solution for a given element.
Definition: XFEM.C:2002
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const =0
bool addFragEdgeIntersection(unsigned int elemid, unsigned int frag_edge_id, double position)
virtual bool isFinalCut() const
Definition: EFAElement2D.C:796
Xfem::XFEM_QRULE _XFEM_qrule
Definition: XFEM.h:322
const QBase *const & qRule() const
void storeSolutionForNode(const Node *node_to_store_to, const Node *node_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > &current_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
Store the solution in stored_solution for a given node.
Definition: XFEM.C:1968
bool isThirdInteriorFace(unsigned int face_id) const
unsigned int number() const
virtual void getXFEMqRuleOnLine(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
Definition: XFEM.C:1897
Data structure describing geometrically described cut through 2D element.
void loadMaterialPropertiesForElementHelper(const Elem *elem, const Xfem::CachedMaterialProperties &cached_props, MaterialPropertyStorage &storage) const
Load the material properties.
Definition: XFEM.C:2213
void setDebugOutputLevel(unsigned int debug_output_level)
Controls amount of debugging information output.
Definition: XFEM.C:1837
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:335
virtual void close()=0
virtual void getXFEMIntersectionInfo(const Elem *elem, unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, bool displaced_mesh=false) const
Definition: XFEM.C:1878
std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo2D > > _geom_marked_elems_2d
Data structure for storing information about all 2D elements to be cut by geometry.
Definition: XFEM.h:344
std::map< unsigned int, std::set< unsigned int > > _geom_marker_id_elems
Data structure for storing the elements cut by specific geometric cutters.
Definition: XFEM.h:350
void addFragFaceIntersection(unsigned int ElemID, unsigned int FragFaceID, const std::vector< unsigned int > &FragFaceEdgeID, const std::vector< double > &position)
XFEM_CUTPLANE_QUANTITY
Definition: XFEM.h:27
virtual CutSubdomainID getCutSubdomainID(const Node *) const
Get CutSubdomainID telling which side the node belongs to relative to the cut.
bool hasStatefulProperties() const
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual void initSolution(const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override
Definition: XFEM.C:305
std::set< const Elem * > _state_marked_frags
Definition: XFEM.h:340
EFANode * getTipEmbeddedNode() const
void restoreFragmentInfo(EFAElement *const elem, const EFAElement *const from_elem)
void storeMaterialPropertiesForElementHelper(const Elem *elem, MaterialPropertyStorage &storage)
Definition: XFEM.C:2167
MeshBase * _mesh
const PropsType & props(const unsigned int state=0) const
virtual void getMasterInfo(EFANode *node, std::vector< EFANode *> &master_nodes, std::vector< double > &master_weights) const =0
const ExecFlagType EXEC_XFEM_MARK
Exec flag used to execute MooseObjects while elements are being marked for cutting by XFEM...
Definition: XFEMAppTypes.C:13
EFAElement3D * getEFAElem3D(const Elem *elem)
Get the EFAElement3D object for a specified libMesh element.
Definition: XFEM.C:1758
EFAElement2D * getEdgeNeighbor(unsigned int edge_id, unsigned int neighbor_id) const
virtual void computePhysicalVolumeFraction()=0
Computes the volume fraction of the element fragment.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool isElemAtCrackTip(const Elem *elem) const
Definition: XFEM.C:1700
MooseMesh * _moose_mesh
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:96
std::vector< dof_id_type > getElementSolutionDofs(const Elem *elem, SystemBase &sys) const
Get a vector of the dof indices for all components of all variables associated with an element...
Definition: XFEM.C:2097
virtual System & system() override
auto norm_sq(const T &a) -> decltype(std::norm(a))
Xfem::XFEM_QRULE & getXFEMQRule()
Definition: XFEM.C:1813
OStreamProxy out
Real getPhysicalVolumeFraction(const Elem *elem) const
Get the volume fraction of an element that is physical.
Definition: XFEM.C:1628
virtual bool getXFEMWeights(MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points) override
Definition: XFEM.C:1843
virtual void getXFEMqRuleOnSurface(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
Definition: XFEM.C:1924
EFAElement * add2DElement(const std::vector< unsigned int > &quad, unsigned int id)
bool markCutFacesByGeometry()
Definition: XFEM.C:849
unsigned int getLocalNodeIndex(EFANode *node) const
Definition: EFAElement.C:111
void dataStore(std::ostream &stream, FeatureFloodCount::FeatureData &feature, void *context)
void setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment)
Definition: XFEM.C:1830
void addStateMarkedElem(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:111
EFAElement * getElemByID(unsigned int id)
MeshBase * _displaced_mesh
const std::set< SubdomainID > & getSubdomainsForVar(unsigned int var_number) const
std::map< unique_id_type, std::vector< Real > > _cached_solution
Data structure to store the nonlinear solution for nodes/elements affected by XFEM For each node/elem...
Definition: XFEM.h:371
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
bool isElemCut(const Elem *elem, XFEMCutElem *&xfce) const
Definition: XFEM.C:1706
virtual void set(const numeric_index_type i, const Number value)=0
virtual bool updateHeal() override
Definition: XFEM.C:223
const ConsoleStream _console
virtual void serializeSolution()
void clearStateMarkedElems()
Definition: XFEM.C:149
void getWeightMultipliers(MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:63
virtual bool isTransient() const override
virtual bool cutElementByCrackGrowthIncrement(const Elem *elem, std::vector< CutEdgeForCrackGrowthIncr > &cut_edges, Real time)
void setXFEMQRule(std::string &xfem_qrule)
Definition: XFEM.C:1819
void buildEFAMesh()
Definition: XFEM.C:336
EFAFace * getFragmentFace(unsigned int frag_id, unsigned int face_id) const
NumericVector< Number > & solutionOld()
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Definition: XFEM.h:341
Assembly & assembly(const THREAD_ID tid, const unsigned int nl_sys_num) override
void dataLoad(std::istream &stream, FeatureFloodCount::FeatureData &feature, void *context)
Information about a cut element.
Definition: XFEM.h:57
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:181
uint8_t unique_id_type
const std::set< EFAElement * > & getCrackTipElements()
static const std::string k
Definition: NS.h:124
unsigned int numFaceNeighbors(unsigned int face_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=nullptr) const
Definition: XFEM.C:1596
unsigned int id() const
Definition: EFANode.C:36
std::map< const Elem *, RealVectorValue > _state_marked_elems
Definition: XFEM.h:339
EFAElement2D * getEFAElem2D(const Elem *elem)
Get the EFAElement2D object for a specified libMesh element.
Definition: XFEM.C:1746
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const =0
void ErrorVector unsigned int
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
bool markCutFacesByState()
Definition: XFEM.C:890
void setInterfaceID(unsigned int interface_id)
Set the interface ID for this cutting object.
void addStateMarkedFrag(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:135
unsigned int numFaces() const
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const =0
uint8_t dof_id_type
std::unordered_map< const Elem *, Xfem::CutElemInfo > _geom_cut_elems
All geometrically cut elements and their CutElemInfo during the current execution of XFEM_MARK...
Definition: XFEM.h:386
bool isPointPhysical(const Point &p) const
Definition: XFEMCutElem.C:193
std::unordered_map< const Elem *, Xfem::CutElemInfo > _old_geom_cut_elems
All geometrically cut elements and their CutElemInfo before the current execution of XFEM_MARK...
Definition: XFEM.h:392