www.mooseframework.org
XFEM.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "XFEM.h"
9 
10 // XFEM includes
11 #include "XFEMCutElem2D.h"
12 #include "XFEMCutElem3D.h"
13 #include "XFEMFuncs.h"
14 #include "EFANode.h"
15 #include "EFAEdge.h"
16 #include "EFAFace.h"
17 #include "EFAFragment2D.h"
18 #include "EFAFragment3D.h"
19 #include "EFAFuncs.h"
20 
21 // MOOSE includes
22 #include "AuxiliarySystem.h"
23 #include "MooseVariable.h"
24 #include "NonlinearSystem.h"
25 #include "FEProblem.h"
26 
27 #include "libmesh/mesh_communication.h"
28 
29 XFEM::XFEM(const InputParameters & params) : XFEMInterface(params), _efa_mesh(Moose::out)
30 {
31 #ifndef LIBMESH_ENABLE_UNIQUE_ID
32  mooseError("MOOSE requires unique ids to be enabled in libmesh (configure with "
33  "--enable-unique-id) to use XFEM!");
34 #endif
35  _has_secondary_cut = false;
36 }
37 
39 {
40  for (std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.begin();
41  cemit != _cut_elem_map.end();
42  ++cemit)
43  delete cemit->second;
44 }
45 
46 void
48 {
49  _geometric_cuts.push_back(geometric_cut);
50 }
51 
52 void
53 XFEM::getCrackTipOrigin(std::map<unsigned int, const Elem *> & elem_id_crack_tip,
54  std::vector<Point> & crack_front_points)
55 {
56  elem_id_crack_tip.clear();
57  crack_front_points.clear();
58  crack_front_points.resize(_elem_crack_origin_direction_map.size());
59 
60  unsigned int crack_tip_index = 0;
61  // This map is used to sort the order in _elem_crack_origin_direction_map such that every process
62  // has same order
63  std::map<unsigned int, const Elem *> elem_id_map;
64 
65  int m = -1;
66  for (std::map<const Elem *, std::vector<Point>>::iterator mit1 =
69  ++mit1)
70  {
71  unsigned int elem_id = mit1->first->id();
72  if (elem_id > 999999)
73  {
74  elem_id_map[m] = mit1->first;
75  m--;
76  }
77  else
78  {
79  elem_id_map[elem_id] = mit1->first;
80  }
81  }
82 
83  for (std::map<unsigned int, const Elem *>::iterator mit1 = elem_id_map.begin();
84  mit1 != elem_id_map.end();
85  mit1++)
86  {
87  const Elem * elem = mit1->second;
88  std::map<const Elem *, std::vector<Point>>::iterator mit2 =
90  if (mit2 != _elem_crack_origin_direction_map.end())
91  {
92  elem_id_crack_tip[crack_tip_index] = mit2->first;
93  crack_front_points[crack_tip_index] =
94  (mit2->second)[0]; // [0] stores origin coordinates and [1] stores direction
95  crack_tip_index++;
96  }
97  }
98 }
99 
100 void
101 XFEM::addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal)
102 {
103  Elem * elem = _mesh->elem(elem_id);
104  std::map<const Elem *, RealVectorValue>::iterator mit;
105  mit = _state_marked_elems.find(elem);
106  if (mit != _state_marked_elems.end())
107  mooseError(" ERROR: element ", elem->id(), " already marked for crack growth.");
108  _state_marked_elems[elem] = normal;
109 }
110 
111 void
112 XFEM::addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal, unsigned int marked_side)
113 {
114  addStateMarkedElem(elem_id, normal);
115  Elem * elem = _mesh->elem(elem_id);
116  std::map<const Elem *, unsigned int>::iterator mit;
117  mit = _state_marked_elem_sides.find(elem);
118  if (mit != _state_marked_elem_sides.end())
119  {
120  mooseError(" ERROR: side of element ", elem->id(), " already marked for crack initiation.");
121  exit(1);
122  }
123  _state_marked_elem_sides[elem] = marked_side;
124 }
125 
126 void
127 XFEM::addStateMarkedFrag(unsigned int elem_id, RealVectorValue & normal)
128 {
129  addStateMarkedElem(elem_id, normal);
130  Elem * elem = _mesh->elem(elem_id);
131  std::set<const Elem *>::iterator mit;
132  mit = _state_marked_frags.find(elem);
133  if (mit != _state_marked_frags.end())
134  {
135  mooseError(
136  " ERROR: element ", elem->id(), " already marked for fragment-secondary crack initiation.");
137  exit(1);
138  }
139  _state_marked_frags.insert(elem);
140 }
141 
142 void
144 {
145  _state_marked_elems.clear();
146  _state_marked_frags.clear();
147  _state_marked_elem_sides.clear();
148 }
149 
150 void
152 {
154  std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
155  std::set<EFAElement *>::iterator sit;
156  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
157  {
158  if (_mesh->mesh_dimension() == 2)
159  {
160  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(*sit);
161  EFANode * tip_node = CEMElem->getTipEmbeddedNode();
162  unsigned int cts_id = CEMElem->getCrackTipSplitElementID();
163 
164  Point origin(0, 0, 0);
165  Point direction(0, 0, 0);
166 
167  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
168  it = _cut_elem_map.find(_mesh->elem(cts_id)->unique_id());
169  if (it != _cut_elem_map.end())
170  {
171  const XFEMCutElem * xfce = it->second;
172  const EFAElement * EFAelem = xfce->getEFAElement();
173  if (EFAelem->isPartial()) // exclude the full crack tip elements
174  {
175  xfce->getCrackTipOriginAndDirection(tip_node->id(), origin, direction);
176  }
177  }
178 
179  std::vector<Point> tip_data;
180  tip_data.push_back(origin);
181  tip_data.push_back(direction);
182  const Elem * elem = _mesh->elem((*sit)->id());
184  std::pair<const Elem *, std::vector<Point>>(elem, tip_data));
185  }
186  }
187 }
188 
189 bool
190 XFEM::update(Real time, NonlinearSystemBase & nl, AuxiliarySystem & aux)
191 {
192  bool mesh_changed = false;
193 
194  buildEFAMesh();
195 
197 
198  if (markCuts(time))
199  mesh_changed = cutMeshWithEFA(nl, aux);
200 
201  if (mesh_changed)
202  {
203  buildEFAMesh();
205  }
206 
207  if (mesh_changed)
208  {
209  _mesh->update_parallel_id_counts();
210  MeshCommunication().make_elems_parallel_consistent(*_mesh);
211  MeshCommunication().make_nodes_parallel_consistent(*_mesh);
212  // _mesh->find_neighbors();
213  // _mesh->contract();
214  _mesh->allow_renumbering(false);
215  _mesh->skip_partitioning(true);
216  _mesh->prepare_for_use();
217  // _mesh->prepare_for_use(true,true); //doing this preserves the numbering, but generates
218  // warning
219 
220  if (_displaced_mesh)
221  {
222  _displaced_mesh->update_parallel_id_counts();
223  MeshCommunication().make_elems_parallel_consistent(*_displaced_mesh);
224  MeshCommunication().make_nodes_parallel_consistent(*_displaced_mesh);
225  _displaced_mesh->allow_renumbering(false);
226  _displaced_mesh->skip_partitioning(true);
227  _displaced_mesh->prepare_for_use();
228  // _displaced_mesh->prepare_for_use(true,true);
229  }
230  }
231 
233 
234  return mesh_changed;
235 }
236 
237 void
238 XFEM::initSolution(NonlinearSystemBase & nl, AuxiliarySystem & aux)
239 {
240  nl.serializeSolution();
241  aux.serializeSolution();
242  NumericVector<Number> & current_solution = *nl.system().current_local_solution;
243  NumericVector<Number> & old_solution = nl.solutionOld();
244  NumericVector<Number> & older_solution = nl.solutionOlder();
245  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
246  NumericVector<Number> & old_aux_solution = aux.solutionOld();
247  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
248 
249  setSolution(nl, _cached_solution, current_solution, old_solution, older_solution);
250  setSolution(
251  aux, _cached_aux_solution, current_aux_solution, old_aux_solution, older_aux_solution);
252 
253  current_solution.close();
254  old_solution.close();
255  older_solution.close();
256  current_aux_solution.close();
257  old_aux_solution.close();
258  older_aux_solution.close();
259 
260  _cached_solution.clear();
261  _cached_aux_solution.clear();
262 }
263 
264 void
266 {
267  _efa_mesh.reset();
268 
269  MeshBase::element_iterator elem_it = _mesh->elements_begin();
270  const MeshBase::element_iterator elem_end = _mesh->elements_end();
271 
272  // Load all existing elements in to EFA mesh
273  for (elem_it = _mesh->elements_begin(); elem_it != elem_end; ++elem_it)
274  {
275  Elem * elem = *elem_it;
276  std::vector<unsigned int> quad;
277  for (unsigned int i = 0; i < elem->n_nodes(); ++i)
278  quad.push_back(elem->node(i));
279  if (_mesh->mesh_dimension() == 2)
280  _efa_mesh.add2DElement(quad, elem->id());
281  else if (_mesh->mesh_dimension() == 3)
282  _efa_mesh.add3DElement(quad, elem->id());
283  else
284  mooseError("XFEM only works for 2D and 3D");
285  }
286 
287  // Restore fragment information for elements that have been previously cut
288  for (elem_it = _mesh->elements_begin(); elem_it != elem_end; ++elem_it)
289  {
290  Elem * elem = *elem_it;
291  std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.find(elem->unique_id());
292  if (cemit != _cut_elem_map.end())
293  {
294  XFEMCutElem * xfce = cemit->second;
295  EFAElement * CEMElem = _efa_mesh.getElemByID(elem->id());
296  _efa_mesh.restoreFragmentInfo(CEMElem, xfce->getEFAElement());
297  }
298  }
299 
300  // Must update edge neighbors before restore edge intersections. Otherwise, when we
301  // add edge intersections, we do not have neighbor information to use.
302  // Correction: no need to use neighbor info now
305 }
306 
307 bool
308 XFEM::markCuts(Real time)
309 {
310  bool marked_sides = false;
311  if (_mesh->mesh_dimension() == 2)
312  {
313  marked_sides = markCutEdgesByGeometry(time);
314  marked_sides |= markCutEdgesByState(time);
315  }
316  else if (_mesh->mesh_dimension() == 3)
317  {
318  marked_sides = markCutFacesByGeometry(time);
319  marked_sides |= markCutFacesByState();
320  }
321  return marked_sides;
322 }
323 
324 bool
326 {
327  bool marked_edges = false;
328  bool marked_nodes = false;
329 
330  std::vector<const GeometricCutUserObject *> active_geometric_cuts;
331  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
332  if (_geometric_cuts[i]->active(time))
333  {
334  active_geometric_cuts.push_back(_geometric_cuts[i]);
335  }
336 
337  if (active_geometric_cuts.size() > 0)
338  {
339  for (MeshBase::element_iterator elem_it = _mesh->elements_begin();
340  elem_it != _mesh->elements_end();
341  ++elem_it)
342  {
343  const Elem * elem = *elem_it;
344  std::vector<CutEdge> elem_cut_edges;
345  std::vector<CutNode> elem_cut_nodes;
346  std::vector<CutEdge> frag_cut_edges;
347  std::vector<std::vector<Point>> frag_edges;
348  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
349  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
350 
351  if (!CEMElem)
352  mooseError("EFAelem is not of EFAelement2D type");
353 
354  // continue if elem has been already cut twice - IMPORTANT
355  if (CEMElem->isFinalCut())
356  continue;
357 
358  // get fragment edges
359  getFragmentEdges(elem, CEMElem, frag_edges);
360 
361  // mark cut edges for the element and its fragment
362  for (unsigned int i = 0; i < active_geometric_cuts.size(); ++i)
363  {
364  active_geometric_cuts[i]->cutElementByGeometry(elem, elem_cut_edges, elem_cut_nodes, time);
365  if (CEMElem->numFragments() > 0)
366  active_geometric_cuts[i]->cutFragmentByGeometry(frag_edges, frag_cut_edges, time);
367  }
368 
369  for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
370  {
371  if (!CEMElem->isEdgePhantom(elem_cut_edges[i].host_side_id)) // must not be phantom edge
372  {
374  elem->id(), elem_cut_edges[i].host_side_id, elem_cut_edges[i].distance);
375  marked_edges = true;
376  }
377  }
378 
379  for (unsigned int i = 0; i < elem_cut_nodes.size(); ++i) // mark element edges
380  {
381  _efa_mesh.addElemNodeIntersection(elem->id(), elem_cut_nodes[i].host_id);
382  marked_nodes = true;
383  }
384 
385  for (unsigned int i = 0; i < frag_cut_edges.size();
386  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
387  {
388  if (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(frag_cut_edges[i].host_side_id))
389  {
391  elem->id(), frag_cut_edges[i].host_side_id, frag_cut_edges[i].distance))
392  {
393  marked_edges = true;
394 
395  if (!isElemAtCrackTip(elem))
396  _has_secondary_cut = true;
397  }
398  }
399  }
400  }
401  }
402 
403  return marked_edges || marked_nodes;
404 }
405 
406 void
408  EFAElement2D * CEMElem,
409  EFAEdge * orig_edge,
410  Point normal,
411  Point crack_tip_origin,
412  Point crack_tip_direction,
413  Real & distance_keep,
414  unsigned int & edge_id_keep,
415  Point & normal_keep)
416 {
417  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
418  Point edge1(0.0, 0.0, 0.0);
419  Point edge2(0.0, 0.0, 0.0);
420  Point left_angle(0.0, 0.0, 0.0);
421  Point right_angle(0.0, 0.0, 0.0);
422  Point left_angle_normal(0.0, 0.0, 0.0);
423  Point right_angle_normal(0.0, 0.0, 0.0);
424  Point crack_direction_normal(0.0, 0.0, 0.0);
425  Point edge1_to_tip(0.0, 0.0, 0.0);
426  Point edge2_to_tip(0.0, 0.0, 0.0);
427  Point edge1_to_tip_normal(0.0, 0.0, 0.0);
428  Point edge2_to_tip_normal(0.0, 0.0, 0.0);
429 
430  Real cos_45 = std::cos(45.0 / 180.0 * 3.14159);
431  Real sin_45 = std::sin(45.0 / 180.0 * 3.14159);
432 
433  left_angle(0) = cos_45 * crack_tip_direction(0) - sin_45 * crack_tip_direction(1);
434  left_angle(1) = sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
435 
436  right_angle(0) = cos_45 * crack_tip_direction(0) + sin_45 * crack_tip_direction(1);
437  right_angle(1) = -sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
438 
439  left_angle_normal(0) = -left_angle(1);
440  left_angle_normal(1) = left_angle(0);
441 
442  right_angle_normal(0) = -right_angle(1);
443  right_angle_normal(1) = right_angle(0);
444 
445  crack_direction_normal(0) = -crack_tip_direction(1);
446  crack_direction_normal(1) = crack_tip_direction(0);
447 
448  Real angle_min = 0.0;
449  Real distance = 0.0;
450  unsigned int nsides = CEMElem->numEdges();
451 
452  for (unsigned int i = 0; i < nsides; ++i)
453  {
454  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
455  {
456  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
457  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
458 
459  edge1_to_tip = (edge_ends[0] * 0.95 + edge_ends[1] * 0.05) - crack_tip_origin;
460  edge2_to_tip = (edge_ends[0] * 0.05 + edge_ends[1] * 0.95) - crack_tip_origin;
461 
462  edge1_to_tip /= pow(edge1_to_tip.norm_sq(), 0.5);
463  edge2_to_tip /= pow(edge2_to_tip.norm_sq(), 0.5);
464 
465  edge1_to_tip_normal(0) = -edge1_to_tip(1);
466  edge1_to_tip_normal(1) = edge1_to_tip(0);
467 
468  edge2_to_tip_normal(0) = -edge2_to_tip(1);
469  edge2_to_tip_normal(1) = edge2_to_tip(0);
470 
471  Real angle_edge1_normal = edge1_to_tip_normal * normal;
472  Real angle_edge2_normal = edge2_to_tip_normal * normal;
473 
474  if (std::abs(angle_edge1_normal) > std::abs(angle_min) &&
475  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
476  {
477  edge_id_keep = i;
478  distance_keep = 0.05;
479  normal_keep = edge1_to_tip_normal;
480  angle_min = angle_edge1_normal;
481  }
482  else if (std::abs(angle_edge2_normal) > std::abs(angle_min) &&
483  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
484  {
485  edge_id_keep = i;
486  distance_keep = 0.95;
487  normal_keep = edge2_to_tip_normal;
488  angle_min = angle_edge2_normal;
489  }
490 
492  crack_tip_origin, left_angle_normal, edge_ends[0], edge_ends[1], distance) &&
493  (!CEMElem->isEdgePhantom(i)))
494  {
495  if (std::abs(left_angle_normal * normal) > std::abs(angle_min) &&
496  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
497  {
498  edge_id_keep = i;
499  distance_keep = distance;
500  normal_keep = left_angle_normal;
501  angle_min = left_angle_normal * normal;
502  }
503  }
504  else if (initCutIntersectionEdge(
505  crack_tip_origin, right_angle_normal, edge_ends[0], edge_ends[1], distance) &&
506  (!CEMElem->isEdgePhantom(i)))
507  {
508  if (std::abs(right_angle_normal * normal) > std::abs(angle_min) &&
509  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
510  {
511  edge_id_keep = i;
512  distance_keep = distance;
513  normal_keep = right_angle_normal;
514  angle_min = right_angle_normal * normal;
515  }
516  }
517  else if (initCutIntersectionEdge(crack_tip_origin,
518  crack_direction_normal,
519  edge_ends[0],
520  edge_ends[1],
521  distance) &&
522  (!CEMElem->isEdgePhantom(i)))
523  {
524  if (std::abs(crack_direction_normal * normal) > std::abs(angle_min) &&
525  (crack_tip_direction * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
526  {
527  edge_id_keep = i;
528  distance_keep = distance;
529  normal_keep = crack_direction_normal;
530  angle_min = crack_direction_normal * normal;
531  }
532  }
533  }
534  }
535 
536  // avoid small volume fraction cut
537  if ((distance_keep - 0.05) < 0.0)
538  {
539  distance_keep = 0.05;
540  }
541  else if ((distance_keep - 0.95) > 0.0)
542  {
543  distance_keep = 0.95;
544  }
545 }
546 
547 bool
549 {
550  bool marked_edges = false;
551  for (std::map<const Elem *, RealVectorValue>::iterator pmeit = _state_marked_elems.begin();
552  pmeit != _state_marked_elems.end();
553  ++pmeit)
554  {
555  const Elem * elem = pmeit->first;
556  RealVectorValue normal = pmeit->second;
557  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
558  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
559 
560  Real volfrac_elem = getPhysicalVolumeFraction(elem);
561  if (volfrac_elem < 0.25)
562  continue;
563 
564  if (!CEMElem)
565  mooseError("EFAelem is not of EFAelement2D type");
566 
567  // continue if elem is already cut twice - IMPORTANT
568  if (CEMElem->isFinalCut())
569  continue;
570 
571  // find the first cut edge
572  unsigned int nsides = CEMElem->numEdges();
573  unsigned int orig_cut_side_id = 999999;
574  Real orig_cut_distance = -1.0;
575  EFANode * orig_node = NULL;
576  EFAEdge * orig_edge = NULL;
577 
578  // crack tip origin coordinates and direction
579  Point crack_tip_origin(0, 0, 0);
580  Point crack_tip_direction(0, 0, 0);
581 
582  if (isElemAtCrackTip(elem)) // crack tip element's crack intiation
583  {
584  orig_cut_side_id = CEMElem->getTipEdgeID();
585  if (orig_cut_side_id < nsides) // valid crack-tip edge found
586  {
587  orig_edge = CEMElem->getEdge(orig_cut_side_id);
588  orig_node = CEMElem->getTipEmbeddedNode();
589  }
590  else
591  mooseError("element ", elem->id(), " has no valid crack-tip edge");
592 
593  // obtain the crack tip origin coordinates and direction.
594  std::map<const Elem *, std::vector<Point>>::iterator ecodm =
596  if (ecodm != _elem_crack_origin_direction_map.end())
597  {
598  crack_tip_origin = (ecodm->second)[0];
599  crack_tip_direction = (ecodm->second)[1];
600  }
601  else
602  mooseError("element ", elem->id(), " cannot find its crack tip origin and direction.");
603  }
604  else
605  {
606  std::map<const Elem *, unsigned int>::iterator mit1;
607  mit1 = _state_marked_elem_sides.find(elem);
608  std::set<const Elem *>::iterator mit2;
609  mit2 = _state_marked_frags.find(elem);
610 
611  if (mit1 != _state_marked_elem_sides.end()) // specified boundary crack initiation
612  {
613  orig_cut_side_id = mit1->second;
614  if (!CEMElem->isEdgePhantom(orig_cut_side_id) &&
615  !CEMElem->getEdge(orig_cut_side_id)->hasIntersection())
616  {
617  orig_cut_distance = 0.5;
618  _efa_mesh.addElemEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
619  orig_edge = CEMElem->getEdge(orig_cut_side_id);
620  orig_node = orig_edge->getEmbeddedNode(0);
621  // get a virtual crack tip direction
622  Point elem_center(0.0, 0.0, 0.0);
623  Point edge_center;
624  for (unsigned int i = 0; i < nsides; ++i)
625  {
626  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
627  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
628  }
629  elem_center /= nsides * 2.0;
630  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
631  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
632  edge_center /= 2.0;
633  crack_tip_origin = edge_center;
634  crack_tip_direction = elem_center - edge_center;
635  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
636  }
637  else
638  continue; // skip this elem if specified boundary edge is phantom
639  }
640  else if (mit2 != _state_marked_frags.end()) // cut-surface secondary crack initiation
641  {
642  if (CEMElem->numFragments() != 1)
643  mooseError("element ",
644  elem->id(),
645  " flagged for a secondary crack, but has ",
646  CEMElem->numFragments(),
647  " fragments");
648  std::vector<unsigned int> interior_edge_id = CEMElem->getFragment(0)->getInteriorEdgeID();
649  if (interior_edge_id.size() == 1)
650  orig_cut_side_id = interior_edge_id[0];
651  else
652  continue; // skip this elem if more than one interior edges found (i.e. elem's been cut
653  // twice)
654  orig_cut_distance = 0.5;
655  _efa_mesh.addFragEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
656  orig_edge = CEMElem->getFragmentEdge(0, orig_cut_side_id);
657  orig_node = orig_edge->getEmbeddedNode(0); // must be an interior embedded node
658  Point elem_center(0.0, 0.0, 0.0);
659  Point edge_center;
660  unsigned int nsides_frag = CEMElem->getFragment(0)->numEdges();
661  for (unsigned int i = 0; i < nsides_frag; ++i)
662  {
663  elem_center +=
664  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
665  elem_center +=
666  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
667  }
668  elem_center /= nsides_frag * 2.0;
669  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
670  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
671  edge_center /= 2.0;
672  crack_tip_origin = edge_center;
673  crack_tip_direction = elem_center - edge_center;
674  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
675  }
676  else
677  mooseError("element ",
678  elem->id(),
679  " flagged for state-based growth, but has no edge intersections");
680  }
681 
682  Point cut_origin(0.0, 0.0, 0.0);
683  if (orig_node)
684  cut_origin = getEFANodeCoords(orig_node, CEMElem, elem); // cutting plane origin's coords
685  else
686  mooseError("element ", elem->id(), " does not have valid orig_node");
687 
688  // loop through element edges to add possible second cut points
689  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
690  Point edge1(0.0, 0.0, 0.0);
691  Point edge2(0.0, 0.0, 0.0);
692  Point cut_edge_point(0.0, 0.0, 0.0);
693  bool find_compatible_direction = false;
694  unsigned int edge_id_keep = 0;
695  Real distance_keep = 0.0;
696  Point normal_keep(0.0, 0.0, 0.0);
697  Real distance = 0.0;
698  bool edge_cut = false;
699 
700  for (unsigned int i = 0; i < nsides; ++i)
701  {
702  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
703  {
704  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
705  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
707  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
708  (!CEMElem->isEdgePhantom(i))))
709  {
710  cut_edge_point = distance * edge_ends[1] + (1.0 - distance) * edge_ends[0];
711  distance_keep = distance;
712  edge_id_keep = i;
713  normal_keep = normal;
714  edge_cut = true;
715  break;
716  }
717  }
718  }
719 
720  Point between_two_cuts = (cut_edge_point - crack_tip_origin);
721  between_two_cuts /= pow(between_two_cuts.norm_sq(), 0.5);
722  Real angle_between_two_cuts = between_two_cuts * crack_tip_direction;
723 
724  if (angle_between_two_cuts > std::cos(45.0 / 180.0 * 3.14159)) // original cut direction is good
725  find_compatible_direction = true;
726 
727  if (!find_compatible_direction && edge_cut)
729  CEMElem,
730  orig_edge,
731  normal,
732  crack_tip_origin,
733  crack_tip_direction,
734  distance_keep,
735  edge_id_keep,
736  normal_keep);
737 
738  if (edge_cut)
739  {
741  _efa_mesh.addElemEdgeIntersection(elem->id(), edge_id_keep, distance_keep);
742  else
743  {
744  MeshBase::element_iterator elem_it = _mesh->elements_begin();
745  const MeshBase::element_iterator elem_end = _mesh->elements_end();
746 
747  Point growth_direction(0.0, 0.0, 0.0);
748 
749  growth_direction(0) = -normal_keep(1);
750  growth_direction(1) = normal_keep(0);
751 
752  if (growth_direction * crack_tip_direction < 1.0e-10)
753  growth_direction *= (-1.0);
754 
755  Real x0 = crack_tip_origin(0);
756  Real y0 = crack_tip_origin(1);
757  Real x1 = x0 + _crack_growth_increment * growth_direction(0);
758  Real y1 = y0 + _crack_growth_increment * growth_direction(1);
759 
760  XFEMCrackGrowthIncrement2DCut geometric_cut(x0, y0, x1, y1, time * 0.9, time * 0.9);
761 
762  for (; elem_it != elem_end; ++elem_it)
763  {
764  const Elem * elem = *elem_it;
765  std::vector<CutEdgeForCrackGrowthIncr> elem_cut_edges;
766  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
767  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
768 
769  if (!CEMElem)
770  mooseError("EFAelem is not of EFAelement2D type");
771 
772  // continue if elem has been already cut twice - IMPORTANT
773  if (CEMElem->isFinalCut())
774  continue;
775 
776  // mark cut edges for the element and its fragment
777  geometric_cut.cutElementByCrackGrowthIncrement(elem, elem_cut_edges, time);
778 
779  for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
780  {
781  if (!CEMElem->isEdgePhantom(elem_cut_edges[i].host_side_id)) // must not be phantom edge
782  {
784  elem->id(), elem_cut_edges[i].host_side_id, elem_cut_edges[i].distance);
785  }
786  }
787  }
788  }
789  }
790  // loop though framgent boundary edges to add possible second cut points
791  // N.B. must do this after marking element edges
792  if (CEMElem->numFragments() > 0 && !edge_cut)
793  {
794  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
795  {
796  if (!orig_edge->isPartialOverlap(*CEMElem->getFragmentEdge(0, i)))
797  {
798  edge_ends[0] =
799  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
800  edge_ends[1] =
801  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
803  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
804  (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(i)))
805  {
806  if (_efa_mesh.addFragEdgeIntersection(elem->id(), edge_id_keep, distance_keep))
807  if (!isElemAtCrackTip(elem))
808  _has_secondary_cut = true;
809  break;
810  }
811  }
812  }
813  }
814 
815  marked_edges = true;
816 
817  } // loop over all state_marked_elems
818 
819  return marked_edges;
820 }
821 
822 bool
824 {
825  bool marked_faces = false;
826 
827  MeshBase::element_iterator elem_it = _mesh->elements_begin();
828  const MeshBase::element_iterator elem_end = _mesh->elements_end();
829 
830  std::vector<const GeometricCutUserObject *> active_geometric_cuts;
831  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
832  if (_geometric_cuts[i]->active(time))
833  active_geometric_cuts.push_back(_geometric_cuts[i]);
834 
835  if (active_geometric_cuts.size() > 0)
836  {
837  for (MeshBase::element_iterator elem_it = _mesh->elements_begin();
838  elem_it != _mesh->elements_end();
839  ++elem_it)
840  {
841  const Elem * elem = *elem_it;
842  std::vector<CutFace> elem_cut_faces;
843  std::vector<CutFace> frag_cut_faces;
844  std::vector<std::vector<Point>> frag_faces;
845  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
846  EFAElement3D * CEMElem = dynamic_cast<EFAElement3D *>(EFAelem);
847  if (!CEMElem)
848  mooseError("EFAelem is not of EFAelement3D type");
849 
850  // continue if elem has been already cut twice - IMPORTANT
851  if (CEMElem->isFinalCut())
852  continue;
853 
854  // get fragment faces
855  getFragmentFaces(elem, CEMElem, frag_faces);
856 
857  // mark cut faces for the element and its fragment
858  for (unsigned int i = 0; i < active_geometric_cuts.size(); ++i)
859  {
860  active_geometric_cuts[i]->cutElementByGeometry(elem, elem_cut_faces, time);
861  // TODO: This would be done for branching, which is not yet supported in 3D
862  // if (CEMElem->numFragments() > 0)
863  // active_geometric_cuts[i]->cutFragmentByGeometry(frag_faces, frag_cut_faces, time);
864  }
865 
866  for (unsigned int i = 0; i < elem_cut_faces.size(); ++i) // mark element faces
867  {
868  if (!CEMElem->isFacePhantom(elem_cut_faces[i].face_id)) // must not be phantom face
869  {
871  elem_cut_faces[i].face_id,
872  elem_cut_faces[i].face_edge,
873  elem_cut_faces[i].position);
874  marked_faces = true;
875  }
876  }
877 
878  for (unsigned int i = 0; i < frag_cut_faces.size();
879  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
880  {
881  if (!CEMElem->getFragment(0)->isThirdInteriorFace(frag_cut_faces[i].face_id))
882  {
884  frag_cut_faces[i].face_id,
885  frag_cut_faces[i].face_edge,
886  frag_cut_faces[i].position);
887  marked_faces = true;
888  }
889  }
890  }
891  }
892 
893  return marked_faces;
894 }
895 
896 bool
898 {
899  bool marked_faces = false;
900  // TODO: need to finish this for 3D problems
901  return marked_faces;
902 }
903 
904 bool
906  Point cut_origin, RealVectorValue cut_normal, Point & edge_p1, Point & edge_p2, Real & dist)
907 {
908  dist = 0.0;
909  bool does_intersect = false;
910  Point origin2p1 = edge_p1 - cut_origin;
911  Real plane2p1 = cut_normal(0) * origin2p1(0) + cut_normal(1) * origin2p1(1);
912  Point origin2p2 = edge_p2 - cut_origin;
913  Real plane2p2 = cut_normal(0) * origin2p2(0) + cut_normal(1) * origin2p2(1);
914 
915  if (plane2p1 * plane2p2 < 0.0)
916  {
917  dist = -plane2p1 / (plane2p2 - plane2p1);
918  does_intersect = true;
919  }
920  return does_intersect;
921 }
922 
923 bool
924 XFEM::cutMeshWithEFA(NonlinearSystemBase & nl, AuxiliarySystem & aux)
925 {
926  std::map<unsigned int, Node *> efa_id_to_new_node;
927  std::map<unsigned int, Node *> efa_id_to_new_node2;
928  std::map<unsigned int, Elem *> efa_id_to_new_elem;
929  _cached_solution.clear();
930  _cached_aux_solution.clear();
931 
933  // DEBUG
934  //_efa_mesh.printMesh();
935 
937  // DEBUG
938  //_efa_mesh.printMesh();
939 
940  const std::vector<EFANode *> new_nodes = _efa_mesh.getNewNodes();
941  const std::vector<EFAElement *> new_elements = _efa_mesh.getChildElements();
942  const std::vector<EFAElement *> delete_elements = _efa_mesh.getParentElements();
943 
944  bool mesh_changed = (new_nodes.size() + new_elements.size() + delete_elements.size() > 0);
945 
946  // Prepare to cache solution on DOFs modified by XFEM
947  if (mesh_changed)
948  {
949  nl.serializeSolution();
950  aux.serializeSolution();
951  }
952  NumericVector<Number> & current_solution = *nl.system().current_local_solution;
953  NumericVector<Number> & old_solution = nl.solutionOld();
954  NumericVector<Number> & older_solution = nl.solutionOlder();
955  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
956  NumericVector<Number> & old_aux_solution = aux.solutionOld();
957  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
958 
959  std::map<Node *, Node *> new_nodes_to_parents;
960 
961  // Add new nodes
962  for (unsigned int i = 0; i < new_nodes.size(); ++i)
963  {
964  unsigned int new_node_id = new_nodes[i]->id();
965  unsigned int parent_id = new_nodes[i]->parent()->id();
966 
967  Node * parent_node = _mesh->node_ptr(parent_id);
968  Node * new_node = Node::build(*parent_node, _mesh->n_nodes()).release();
969  new_node->processor_id() = parent_node->processor_id();
970  _mesh->add_node(new_node);
971 
972  new_nodes_to_parents[new_node] = parent_node;
973 
974  new_node->set_n_systems(parent_node->n_systems());
975  efa_id_to_new_node.insert(std::make_pair(new_node_id, new_node));
976  _console << "XFEM added new node: " << new_node->id() << "\n";
977  if (_displaced_mesh)
978  {
979  const Node * parent_node2 = _displaced_mesh->node_ptr(parent_id);
980  Node * new_node2 = Node::build(*parent_node2, _displaced_mesh->n_nodes()).release();
981  new_node2->processor_id() = parent_node2->processor_id();
982  _displaced_mesh->add_node(new_node2);
983 
984  new_node2->set_n_systems(parent_node2->n_systems());
985  efa_id_to_new_node2.insert(std::make_pair(new_node_id, new_node2));
986  }
987  }
988 
989  // Add new elements
990  std::map<unsigned int, std::vector<const Elem *>> temporary_parent_children_map;
991 
992  for (unsigned int i = 0; i < new_elements.size(); ++i)
993  {
994  unsigned int parent_id = new_elements[i]->getParent()->id();
995  unsigned int efa_child_id = new_elements[i]->id();
996 
997  Elem * parent_elem = _mesh->elem(parent_id);
998  Elem * libmesh_elem = Elem::build(parent_elem->type()).release();
999 
1000  for (ElementPairLocator::ElementPairList::iterator it = _sibling_elems.begin();
1001  it != _sibling_elems.end();
1002  ++it)
1003  {
1004  if (parent_elem == it->first)
1005  it->first = libmesh_elem;
1006  else if (parent_elem == it->second)
1007  it->second = libmesh_elem;
1008  }
1009 
1010  // parent has at least two children
1011  if (new_elements[i]->getParent()->numChildren() > 1)
1012  temporary_parent_children_map[parent_id].push_back(libmesh_elem);
1013 
1014  Elem * parent_elem2 = NULL;
1015  Elem * libmesh_elem2 = NULL;
1016  if (_displaced_mesh)
1017  {
1018  parent_elem2 = _displaced_mesh->elem(parent_id);
1019  libmesh_elem2 = Elem::build(parent_elem2->type()).release();
1020 
1021  for (ElementPairLocator::ElementPairList::iterator it = _sibling_displaced_elems.begin();
1022  it != _sibling_displaced_elems.end();
1023  ++it)
1024  {
1025  if (parent_elem2 == it->first)
1026  it->first = libmesh_elem2;
1027  else if (parent_elem2 == it->second)
1028  it->second = libmesh_elem2;
1029  }
1030  }
1031 
1032  for (unsigned int j = 0; j < new_elements[i]->numNodes(); ++j)
1033  {
1034  unsigned int node_id = new_elements[i]->getNode(j)->id();
1035  Node * libmesh_node;
1036 
1037  std::map<unsigned int, Node *>::iterator nit = efa_id_to_new_node.find(node_id);
1038  if (nit != efa_id_to_new_node.end())
1039  libmesh_node = nit->second;
1040  else
1041  libmesh_node = _mesh->node_ptr(node_id);
1042 
1043  libmesh_elem->set_node(j) = libmesh_node;
1044 
1045  // Store solution for all nodes affected by XFEM (even existing nodes)
1046  if (parent_elem->is_semilocal(_mesh->processor_id()))
1047  {
1048  Node * solution_node = libmesh_node; // Node from which to store solution
1049  if (new_nodes_to_parents.find(libmesh_node) != new_nodes_to_parents.end())
1050  solution_node = new_nodes_to_parents[libmesh_node];
1051 
1052  if ((_moose_mesh->isSemiLocal(solution_node)) ||
1053  (solution_node->processor_id() == _mesh->processor_id()))
1054  {
1055  storeSolutionForNode(libmesh_node,
1056  solution_node,
1057  nl,
1059  current_solution,
1060  old_solution,
1061  older_solution);
1062  storeSolutionForNode(libmesh_node,
1063  solution_node,
1064  aux,
1066  current_aux_solution,
1067  old_aux_solution,
1068  older_aux_solution);
1069  }
1070  }
1071 
1072  Node * parent_node = parent_elem->get_node(j);
1073  std::vector<boundary_id_type> parent_node_boundary_ids =
1074  _mesh->boundary_info->boundary_ids(parent_node);
1075  _mesh->boundary_info->add_node(libmesh_node, parent_node_boundary_ids);
1076 
1077  if (_displaced_mesh)
1078  {
1079  std::map<unsigned int, Node *>::iterator nit2 = efa_id_to_new_node2.find(node_id);
1080  if (nit2 != efa_id_to_new_node2.end())
1081  libmesh_node = nit2->second;
1082  else
1083  libmesh_node = _displaced_mesh->node_ptr(node_id);
1084 
1085  libmesh_elem2->set_node(j) = libmesh_node;
1086 
1087  parent_node = parent_elem2->get_node(j);
1088  parent_node_boundary_ids.clear();
1089  parent_node_boundary_ids = _displaced_mesh->boundary_info->boundary_ids(parent_node);
1090  _displaced_mesh->boundary_info->add_node(libmesh_node, parent_node_boundary_ids);
1091  }
1092  }
1093 
1094  libmesh_elem->set_p_level(parent_elem->p_level());
1095  libmesh_elem->set_p_refinement_flag(parent_elem->p_refinement_flag());
1096  _mesh->add_elem(libmesh_elem);
1097  libmesh_elem->set_n_systems(parent_elem->n_systems());
1098  libmesh_elem->subdomain_id() = parent_elem->subdomain_id();
1099  libmesh_elem->processor_id() = parent_elem->processor_id();
1100 
1101  // TODO: The 0 here is the thread ID. Need to sort out how to do this correctly
1102  // TODO: Also need to copy neighbor material data
1103  if (parent_elem->processor_id() == _mesh->processor_id())
1104  {
1105  (*_material_data)[0]->copy(*libmesh_elem, *parent_elem, 0);
1106  for (unsigned int side = 0; side < parent_elem->n_sides(); ++side)
1107  {
1108  std::vector<boundary_id_type> parent_elem_boundary_ids =
1109  _mesh->boundary_info->boundary_ids(parent_elem, side);
1110  std::vector<boundary_id_type>::iterator it_bd = parent_elem_boundary_ids.begin();
1111  for (; it_bd != parent_elem_boundary_ids.end(); ++it_bd)
1112  {
1113  if (_fe_problem->needMaterialOnSide(*it_bd, 0))
1114  (*_bnd_material_data)[0]->copy(*libmesh_elem, *parent_elem, side);
1115  }
1116  }
1117 
1118  // Store solution for all elements affected by XFEM
1119  storeSolutionForElement(libmesh_elem,
1120  parent_elem,
1121  nl,
1123  current_solution,
1124  old_solution,
1125  older_solution);
1126  storeSolutionForElement(libmesh_elem,
1127  parent_elem,
1128  aux,
1130  current_aux_solution,
1131  old_aux_solution,
1132  older_aux_solution);
1133  }
1134 
1135  // The crack tip origin map is stored before cut, thus the elem should be updated with new
1136  // element.
1137  std::map<const Elem *, std::vector<Point>>::iterator mit =
1138  _elem_crack_origin_direction_map.find(parent_elem);
1139  if (mit != _elem_crack_origin_direction_map.end())
1140  {
1141  std::vector<Point> crack_data = _elem_crack_origin_direction_map[parent_elem];
1143  _elem_crack_origin_direction_map[libmesh_elem] = crack_data;
1144  }
1145 
1146  _console << "XFEM added new element: " << libmesh_elem->id() << "\n";
1147 
1148  XFEMCutElem * xfce = NULL;
1149  if (_mesh->mesh_dimension() == 2)
1150  {
1151  EFAElement2D * new_efa_elem2d = dynamic_cast<EFAElement2D *>(new_elements[i]);
1152  if (!new_efa_elem2d)
1153  mooseError("EFAelem is not of EFAelement2D type");
1154  xfce = new XFEMCutElem2D(libmesh_elem, new_efa_elem2d, (*_material_data)[0]->nQPoints());
1155  }
1156  else if (_mesh->mesh_dimension() == 3)
1157  {
1158  EFAElement3D * new_efa_elem3d = dynamic_cast<EFAElement3D *>(new_elements[i]);
1159  if (!new_efa_elem3d)
1160  mooseError("EFAelem is not of EFAelement3D type");
1161  xfce = new XFEMCutElem3D(libmesh_elem, new_efa_elem3d, (*_material_data)[0]->nQPoints());
1162  }
1163  _cut_elem_map.insert(std::pair<unique_id_type, XFEMCutElem *>(libmesh_elem->unique_id(), xfce));
1164  efa_id_to_new_elem.insert(std::make_pair(efa_child_id, libmesh_elem));
1165 
1166  if (_displaced_mesh)
1167  {
1168  libmesh_elem2->set_p_level(parent_elem2->p_level());
1169  libmesh_elem2->set_p_refinement_flag(parent_elem2->p_refinement_flag());
1170  _displaced_mesh->add_elem(libmesh_elem2);
1171  libmesh_elem2->set_n_systems(parent_elem2->n_systems());
1172  libmesh_elem2->subdomain_id() = parent_elem2->subdomain_id();
1173  libmesh_elem2->processor_id() = parent_elem2->processor_id();
1174  }
1175 
1176  unsigned int n_sides = parent_elem->n_sides();
1177  for (unsigned int side = 0; side < n_sides; ++side)
1178  {
1179  std::vector<boundary_id_type> parent_elem_boundary_ids =
1180  _mesh->boundary_info->boundary_ids(parent_elem, side);
1181  _mesh->boundary_info->add_side(libmesh_elem, side, parent_elem_boundary_ids);
1182  }
1183  if (_displaced_mesh)
1184  {
1185  n_sides = parent_elem2->n_sides();
1186  for (unsigned int side = 0; side < n_sides; ++side)
1187  {
1188  std::vector<boundary_id_type> parent_elem_boundary_ids =
1189  _displaced_mesh->boundary_info->boundary_ids(parent_elem2, side);
1190  _displaced_mesh->boundary_info->add_side(libmesh_elem2, side, parent_elem_boundary_ids);
1191  }
1192  }
1193 
1194  unsigned int n_edges = parent_elem->n_edges();
1195  for (unsigned int edge = 0; edge < n_edges; ++edge)
1196  {
1197  std::vector<boundary_id_type> parent_elem_boundary_ids =
1198  _mesh->boundary_info->edge_boundary_ids(parent_elem, edge);
1199  _mesh->boundary_info->add_edge(libmesh_elem, edge, parent_elem_boundary_ids);
1200  }
1201  if (_displaced_mesh)
1202  {
1203  n_edges = parent_elem2->n_edges();
1204  for (unsigned int edge = 0; edge < n_edges; ++edge)
1205  {
1206  std::vector<boundary_id_type> parent_elem_boundary_ids =
1207  _displaced_mesh->boundary_info->edge_boundary_ids(parent_elem2, edge);
1208  _displaced_mesh->boundary_info->add_edge(libmesh_elem2, edge, parent_elem_boundary_ids);
1209  }
1210  }
1211  }
1212 
1213  // delete elements
1214  for (std::size_t i = 0; i < delete_elements.size(); ++i)
1215  {
1216  Elem * elem_to_delete = _mesh->elem(delete_elements[i]->id());
1217 
1218  // delete the XFEMCutElem object for any elements that are to be deleted
1219  std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
1220  _cut_elem_map.find(elem_to_delete->unique_id());
1221  if (cemit != _cut_elem_map.end())
1222  {
1223  delete cemit->second;
1224  _cut_elem_map.erase(cemit);
1225  }
1226 
1227  elem_to_delete->nullify_neighbors();
1228  _mesh->boundary_info->remove(elem_to_delete);
1229  unsigned int deleted_elem_id = elem_to_delete->id();
1230  _mesh->delete_elem(elem_to_delete);
1231  _console << "XFEM deleted element: " << deleted_elem_id << "\n";
1232 
1233  if (_displaced_mesh)
1234  {
1235  Elem * elem_to_delete2 = _displaced_mesh->elem(delete_elements[i]->id());
1236  elem_to_delete2->nullify_neighbors();
1237  _displaced_mesh->boundary_info->remove(elem_to_delete2);
1238  _displaced_mesh->delete_elem(elem_to_delete2);
1239  }
1240  }
1241 
1242  for (std::map<unsigned int, std::vector<const Elem *>>::iterator it =
1243  temporary_parent_children_map.begin();
1244  it != temporary_parent_children_map.end();
1245  ++it)
1246  {
1247  std::vector<const Elem *> & sibling_elem_vec = it->second;
1248  // TODO: for cut-node case, how to find the sibling elements?
1249  // if (sibling_elem_vec.size() != 2)
1250  // mooseError("Must have exactly 2 sibling elements");
1251  _sibling_elems.push_back(std::make_pair(sibling_elem_vec[0], sibling_elem_vec[1]));
1252  }
1253 
1254  // add sibling elems on displaced mesh
1255  if (_displaced_mesh)
1256  {
1257  for (ElementPairLocator::ElementPairList::iterator it = _sibling_elems.begin();
1258  it != _sibling_elems.end();
1259  ++it)
1260  {
1261  Elem * elem = _displaced_mesh->elem(it->first->id());
1262  Elem * elem_pair = _displaced_mesh->elem(it->second->id());
1263  _sibling_displaced_elems.push_back(std::make_pair(elem, elem_pair));
1264  }
1265  }
1266 
1267  // clear the temporary map
1268  temporary_parent_children_map.clear();
1269 
1270  // Store information about crack tip elements
1271  if (mesh_changed)
1272  {
1273  _crack_tip_elems.clear();
1274  const std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
1275  std::set<EFAElement *>::const_iterator sit;
1276  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
1277  {
1278  unsigned int eid = (*sit)->id();
1279  Elem * crack_tip_elem;
1280  std::map<unsigned int, Elem *>::iterator eit = efa_id_to_new_elem.find(eid);
1281  if (eit != efa_id_to_new_elem.end())
1282  crack_tip_elem = eit->second;
1283  else
1284  crack_tip_elem = _mesh->elem(eid);
1285  _crack_tip_elems.insert(crack_tip_elem);
1286  }
1287  }
1288  _console << std::flush;
1289 
1290  // store virtual nodes
1291  // store cut edge info
1292  return mesh_changed;
1293 }
1294 
1295 Point
1297  EFAElement * CEMElem,
1298  const Elem * elem,
1299  MeshBase * displaced_mesh) const
1300 {
1301  Point node_coor(0.0, 0.0, 0.0);
1302  std::vector<EFANode *> master_nodes;
1303  std::vector<Point> master_points;
1304  std::vector<double> master_weights;
1305 
1306  CEMElem->getMasterInfo(CEMnode, master_nodes, master_weights);
1307  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1308  {
1309  if (master_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT)
1310  {
1311  unsigned int local_node_id = CEMElem->getLocalNodeIndex(master_nodes[i]);
1312  Node * node = elem->get_node(local_node_id);
1313  if (displaced_mesh)
1314  node = displaced_mesh->node_ptr(node->id());
1315  Point node_p((*node)(0), (*node)(1), (*node)(2));
1316  master_points.push_back(node_p);
1317  }
1318  else
1319  mooseError("master nodes must be permanent");
1320  }
1321  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1322  node_coor += master_weights[i] * master_points[i];
1323 
1324  return node_coor;
1325 }
1326 
1327 Real
1328 XFEM::getPhysicalVolumeFraction(const Elem * elem) const
1329 {
1330  Real phys_volfrac = 1.0;
1331  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1332  it = _cut_elem_map.find(elem->unique_id());
1333  if (it != _cut_elem_map.end())
1334  {
1335  XFEMCutElem * xfce = it->second;
1336  const EFAElement * EFAelem = xfce->getEFAElement();
1337  if (EFAelem->isPartial())
1338  { // exclude the full crack tip elements
1340  phys_volfrac = xfce->getPhysicalVolumeFraction();
1341  }
1342  }
1343 
1344  return phys_volfrac;
1345 }
1346 
1347 Real
1348 XFEM::getCutPlane(const Elem * elem,
1349  const Xfem::XFEM_CUTPLANE_QUANTITY quantity,
1350  unsigned int plane_id) const
1351 {
1352  Real comp = 0.0;
1353  Point planedata(0.0, 0.0, 0.0);
1354  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1355  it = _cut_elem_map.find(elem->unique_id());
1356  if (it != _cut_elem_map.end())
1357  {
1358  const XFEMCutElem * xfce = it->second;
1359  const EFAElement * EFAelem = xfce->getEFAElement();
1360  if (EFAelem->isPartial()) // exclude the full crack tip elements
1361  {
1362  if ((unsigned int)quantity < 3)
1363  {
1364  unsigned int index = (unsigned int)quantity;
1365  planedata = xfce->getCutPlaneOrigin(plane_id, _displaced_mesh);
1366  comp = planedata(index);
1367  }
1368  else if ((unsigned int)quantity < 6)
1369  {
1370  unsigned int index = (unsigned int)quantity - 3;
1371  planedata = xfce->getCutPlaneNormal(plane_id, _displaced_mesh);
1372  comp = planedata(index);
1373  }
1374  else
1375  mooseError("In get_cut_plane index out of range");
1376  }
1377  }
1378  return comp;
1379 }
1380 
1381 bool
1382 XFEM::isElemAtCrackTip(const Elem * elem) const
1383 {
1384  return (_crack_tip_elems.find(elem) != _crack_tip_elems.end());
1385 }
1386 
1387 bool
1388 XFEM::isElemCut(const Elem * elem, XFEMCutElem *& xfce) const
1389 {
1390  xfce = NULL;
1391  bool is_cut = false;
1392  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1393  it = _cut_elem_map.find(elem->unique_id());
1394  if (it != _cut_elem_map.end())
1395  {
1396  xfce = it->second;
1397  const EFAElement * EFAelem = xfce->getEFAElement();
1398  if (EFAelem->isPartial()) // exclude the full crack tip elements
1399  is_cut = true;
1400  }
1401  return is_cut;
1402 }
1403 
1404 bool
1405 XFEM::isElemCut(const Elem * elem) const
1406 {
1407  XFEMCutElem * xfce = NULL;
1408  return isElemCut(elem, xfce);
1409 }
1410 
1411 void
1412 XFEM::getFragmentFaces(const Elem * elem,
1413  std::vector<std::vector<Point>> & frag_faces,
1414  bool displaced_mesh) const
1415 {
1416  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1417  it = _cut_elem_map.find(elem->unique_id());
1418  if (it != _cut_elem_map.end())
1419  {
1420  const XFEMCutElem * xfce = it->second;
1421  if (displaced_mesh)
1422  xfce->getFragmentFaces(frag_faces, _displaced_mesh);
1423  else
1424  xfce->getFragmentFaces(frag_faces);
1425  }
1426 }
1427 
1428 void
1429 XFEM::getFragmentEdges(const Elem * elem,
1430  EFAElement2D * CEMElem,
1431  std::vector<std::vector<Point>> & frag_edges) const
1432 {
1433  // N.B. CEMElem here has global EFAnode
1434  frag_edges.clear();
1435  if (CEMElem->numFragments() > 0)
1436  {
1437  if (CEMElem->numFragments() > 1)
1438  mooseError("element ", elem->id(), " has more than one fragments at this point");
1439  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
1440  {
1441  std::vector<Point> p_line(2, Point(0.0, 0.0, 0.0));
1442  p_line[0] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
1443  p_line[1] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
1444  frag_edges.push_back(p_line);
1445  }
1446  }
1447 }
1448 
1449 void
1450 XFEM::getFragmentFaces(const Elem * elem,
1451  EFAElement3D * CEMElem,
1452  std::vector<std::vector<Point>> & frag_faces) const
1453 {
1454  // N.B. CEMElem here has global EFAnode
1455  frag_faces.clear();
1456  if (CEMElem->numFragments() > 0)
1457  {
1458  if (CEMElem->numFragments() > 1)
1459  mooseError("element ", elem->id(), " has more than one fragments at this point");
1460  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numFaces(); ++i)
1461  {
1462  unsigned int num_face_nodes = CEMElem->getFragmentFace(0, i)->numNodes();
1463  std::vector<Point> p_line(num_face_nodes, Point(0.0, 0.0, 0.0));
1464  for (unsigned int j = 0; j < num_face_nodes; ++j)
1465  p_line[j] = getEFANodeCoords(CEMElem->getFragmentFace(0, i)->getNode(j), CEMElem, elem);
1466  frag_faces.push_back(p_line);
1467  }
1468  }
1469 }
1470 
1473 {
1474  return _XFEM_qrule;
1475 }
1476 
1477 void
1478 XFEM::setXFEMQRule(std::string & xfem_qrule)
1479 {
1480  if (xfem_qrule == "volfrac")
1482  else if (xfem_qrule == "moment_fitting")
1484  else if (xfem_qrule == "direct")
1486 }
1487 
1488 void
1489 XFEM::setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment)
1490 {
1491  _use_crack_growth_increment = use_crack_growth_increment;
1492  _crack_growth_increment = crack_growth_increment;
1493 }
1494 
1495 bool
1496 XFEM::getXFEMWeights(MooseArray<Real> & weights,
1497  const Elem * elem,
1498  QBase * qrule,
1499  const MooseArray<Point> & q_points)
1500 {
1501  bool have_weights = false;
1502  XFEMCutElem * xfce = NULL;
1503  if (isElemCut(elem, xfce))
1504  {
1505  mooseAssert(xfce != NULL, "Must have valid XFEMCutElem object here");
1506  xfce->getWeightMultipliers(weights, qrule, getXFEMQRule(), q_points);
1507  have_weights = true;
1508  }
1509  return have_weights;
1510 }
1511 
1512 void
1514  unsigned int plane_id,
1515  Point & normal,
1516  std::vector<Point> & intersectionPoints,
1517  bool displaced_mesh) const
1518 {
1519  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1520  it = _cut_elem_map.find(elem->unique_id());
1521  if (it != _cut_elem_map.end())
1522  {
1523  const XFEMCutElem * xfce = it->second;
1524  if (displaced_mesh)
1525  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints, _displaced_mesh);
1526  else
1527  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints);
1528  }
1529 }
1530 
1531 void
1532 XFEM::getXFEMqRuleOnLine(std::vector<Point> & intersection_points,
1533  std::vector<Point> & quad_pts,
1534  std::vector<Real> & quad_wts) const
1535 {
1536  Point p1 = intersection_points[0];
1537  Point p2 = intersection_points[1];
1538 
1539  // number of quadrature points
1540  std::size_t num_qpoints = 2;
1541 
1542  // quadrature coordinates
1543  Real xi0 = -std::sqrt(1.0 / 3.0);
1544  Real xi1 = std::sqrt(1.0 / 3.0);
1545 
1546  quad_wts.resize(num_qpoints);
1547  quad_pts.resize(num_qpoints);
1548 
1549  Real integ_jacobian = pow((p1 - p2).size_sq(), 0.5) * 0.5;
1550 
1551  quad_wts[0] = 1.0 * integ_jacobian;
1552  quad_wts[1] = 1.0 * integ_jacobian;
1553 
1554  quad_pts[0] = (1.0 - xi0) / 2.0 * p1 + (1.0 + xi0) / 2.0 * p2;
1555  quad_pts[1] = (1.0 - xi1) / 2.0 * p1 + (1.0 + xi1) / 2.0 * p2;
1556 }
1557 
1558 void
1559 XFEM::getXFEMqRuleOnSurface(std::vector<Point> & intersection_points,
1560  std::vector<Point> & quad_pts,
1561  std::vector<Real> & quad_wts) const
1562 {
1563  std::size_t nnd_pe = intersection_points.size();
1564  Point xcrd(0.0, 0.0, 0.0);
1565  for (std::size_t i = 0; i < nnd_pe; ++i)
1566  xcrd += intersection_points[i];
1567  xcrd /= nnd_pe;
1568 
1569  quad_pts.resize(nnd_pe);
1570  quad_wts.resize(nnd_pe);
1571 
1572  Real jac = 0.0;
1573 
1574  for (std::size_t j = 0; j < nnd_pe; ++j) // loop all sub-tris
1575  {
1576  std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
1577  std::vector<Point> subtrig_points(3, Point(0.0, 0.0, 0.0)); // sub-trig nodal coords
1578 
1579  int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
1580  subtrig_points[0] = xcrd;
1581  subtrig_points[1] = intersection_points[j];
1582  subtrig_points[2] = intersection_points[jplus1];
1583 
1584  std::vector<std::vector<Real>> sg2;
1585  Xfem::stdQuadr2D(3, 1, sg2); // get sg2
1586  for (std::size_t l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-trig
1587  {
1588  Xfem::shapeFunc2D(3, sg2[l], subtrig_points, shape, jac, true); // Get shape
1589  std::vector<Real> tsg_line(3, 0.0);
1590  for (std::size_t k = 0; k < 3; ++k) // loop sub-trig nodes
1591  {
1592  tsg_line[0] += shape[k][2] * subtrig_points[k](0);
1593  tsg_line[1] += shape[k][2] * subtrig_points[k](1);
1594  tsg_line[2] += shape[k][2] * subtrig_points[k](2);
1595  }
1596  quad_pts[j + l] = Point(tsg_line[0], tsg_line[1], tsg_line[2]);
1597  quad_wts[j + l] = sg2[l][3] * jac;
1598  }
1599  }
1600 }
1601 
1602 void
1603 XFEM::storeSolutionForNode(const Node * node_to_store_to,
1604  const Node * node_to_store_from,
1605  SystemBase & sys,
1606  std::map<unique_id_type, std::vector<Real>> & stored_solution,
1607  const NumericVector<Number> & current_solution,
1608  const NumericVector<Number> & old_solution,
1609  const NumericVector<Number> & older_solution)
1610 {
1611  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node_to_store_from, sys);
1612  std::vector<Real> stored_solution_scratch;
1613  // Size for current solution, as well as for old, and older solution only for transient case
1614  std::size_t stored_solution_size =
1615  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
1616  stored_solution_scratch.reserve(stored_solution_size);
1617 
1618  // Store in the order defined in stored_solution_dofs first for the current, then for old and
1619  // older if applicable
1620  for (auto dof : stored_solution_dofs)
1621  stored_solution_scratch.push_back(current_solution(dof));
1622 
1623  if (_fe_problem->isTransient())
1624  {
1625  for (auto dof : stored_solution_dofs)
1626  stored_solution_scratch.push_back(old_solution(dof));
1627 
1628  for (auto dof : stored_solution_dofs)
1629  stored_solution_scratch.push_back(older_solution(dof));
1630  }
1631 
1632  if (stored_solution_scratch.size() > 0)
1633  stored_solution[node_to_store_to->unique_id()] = stored_solution_scratch;
1634 }
1635 
1636 void
1637 XFEM::storeSolutionForElement(const Elem * elem_to_store_to,
1638  const Elem * elem_to_store_from,
1639  SystemBase & sys,
1640  std::map<unique_id_type, std::vector<Real>> & stored_solution,
1641  const NumericVector<Number> & current_solution,
1642  const NumericVector<Number> & old_solution,
1643  const NumericVector<Number> & older_solution)
1644 {
1645  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem_to_store_from, sys);
1646  std::vector<Real> stored_solution_scratch;
1647  // Size for current solution, as well as for old, and older solution only for transient case
1648  std::size_t stored_solution_size =
1649  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
1650  stored_solution_scratch.reserve(stored_solution_size);
1651 
1652  // Store in the order defined in stored_solution_dofs first for the current, then for old and
1653  // older if applicable
1654  for (auto dof : stored_solution_dofs)
1655  stored_solution_scratch.push_back(current_solution(dof));
1656 
1657  if (_fe_problem->isTransient())
1658  {
1659  for (auto dof : stored_solution_dofs)
1660  stored_solution_scratch.push_back(old_solution(dof));
1661 
1662  for (auto dof : stored_solution_dofs)
1663  stored_solution_scratch.push_back(older_solution(dof));
1664  }
1665 
1666  if (stored_solution_scratch.size() > 0)
1667  stored_solution[elem_to_store_to->unique_id()] = stored_solution_scratch;
1668 }
1669 
1670 void
1671 XFEM::setSolution(SystemBase & sys,
1672  const std::map<unique_id_type, std::vector<Real>> & stored_solution,
1673  NumericVector<Number> & current_solution,
1674  NumericVector<Number> & old_solution,
1675  NumericVector<Number> & older_solution)
1676 {
1677  const auto nodes_end = _mesh->local_nodes_end();
1678  for (auto node_it = _mesh->local_nodes_begin(); node_it != nodes_end; ++node_it)
1679  {
1680  Node * node = *node_it;
1681  auto mit = stored_solution.find(node->unique_id());
1682  if (mit != stored_solution.end())
1683  {
1684  const std::vector<Real> & stored_node_solution = mit->second;
1685  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node, sys);
1686  setSolutionForDOFs(stored_node_solution,
1687  stored_solution_dofs,
1688  current_solution,
1689  old_solution,
1690  older_solution);
1691  }
1692  }
1693 
1694  const auto elems_end = _mesh->local_elements_end();
1695  for (auto elem_it = _mesh->local_elements_begin(); elem_it != elems_end; ++elem_it)
1696  {
1697  Elem * elem = *elem_it;
1698  auto mit = stored_solution.find(elem->unique_id());
1699  if (mit != stored_solution.end())
1700  {
1701  const std::vector<Real> & stored_elem_solution = mit->second;
1702  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem, sys);
1703  setSolutionForDOFs(stored_elem_solution,
1704  stored_solution_dofs,
1705  current_solution,
1706  old_solution,
1707  older_solution);
1708  }
1709  }
1710 }
1711 
1712 void
1713 XFEM::setSolutionForDOFs(const std::vector<Real> & stored_solution,
1714  const std::vector<dof_id_type> & stored_solution_dofs,
1715  NumericVector<Number> & current_solution,
1716  NumericVector<Number> & old_solution,
1717  NumericVector<Number> & older_solution)
1718 {
1719  // Solution vector is stored first for current, then old and older solutions.
1720  // These are the offsets to the beginning of the old and older solutions in the vector.
1721  const auto old_solution_offset = stored_solution_dofs.size();
1722  const auto older_solution_offset = old_solution_offset * 2;
1723 
1724  for (std::size_t i = 0; i < stored_solution_dofs.size(); ++i)
1725  {
1726  current_solution.set(stored_solution_dofs[i], stored_solution[i]);
1727  if (_fe_problem->isTransient())
1728  {
1729  old_solution.set(stored_solution_dofs[i], stored_solution[old_solution_offset + i]);
1730  older_solution.set(stored_solution_dofs[i], stored_solution[older_solution_offset + i]);
1731  }
1732  }
1733 }
1734 
1735 std::vector<dof_id_type>
1736 XFEM::getElementSolutionDofs(const Elem * elem, SystemBase & sys) const
1737 {
1738  SubdomainID sid = elem->subdomain_id();
1739  const std::vector<MooseVariable *> & vars = sys.getVariables(0);
1740  std::vector<dof_id_type> solution_dofs;
1741  solution_dofs.reserve(vars.size()); // just an approximation
1742  for (auto var : vars)
1743  {
1744  if (!var->isNodal())
1745  {
1746  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
1747  if (var_subdomains.empty() || var_subdomains.find(sid) != var_subdomains.end())
1748  {
1749  unsigned int n_comp = elem->n_comp(sys.number(), var->number());
1750  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
1751  {
1752  dof_id_type elem_dof = elem->dof_number(sys.number(), var->number(), icomp);
1753  solution_dofs.push_back(elem_dof);
1754  }
1755  }
1756  }
1757  }
1758  return solution_dofs;
1759 }
1760 
1761 std::vector<dof_id_type>
1762 XFEM::getNodeSolutionDofs(const Node * node, SystemBase & sys) const
1763 {
1764  const std::set<SubdomainID> & sids = _moose_mesh->getNodeBlockIds(*node);
1765  const std::vector<MooseVariable *> & vars = sys.getVariables(0);
1766  std::vector<dof_id_type> solution_dofs;
1767  solution_dofs.reserve(vars.size()); // just an approximation
1768  for (auto var : vars)
1769  {
1770  if (var->isNodal())
1771  {
1772  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
1773  std::set<SubdomainID> intersect;
1774  set_intersection(var_subdomains.begin(),
1775  var_subdomains.end(),
1776  sids.begin(),
1777  sids.end(),
1778  std::inserter(intersect, intersect.begin()));
1779  if (var_subdomains.empty() || !intersect.empty())
1780  {
1781  unsigned int n_comp = node->n_comp(sys.number(), var->number());
1782  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
1783  {
1784  dof_id_type node_dof = node->dof_number(sys.number(), var->number(), icomp);
1785  solution_dofs.push_back(node_dof);
1786  }
1787  }
1788  }
1789  }
1790  return solution_dofs;
1791 }
std::vector< unsigned int > getInteriorEdgeID() const
unsigned int numNodes() const
Definition: EFAFace.C:85
~XFEM()
Definition: XFEM.C:38
EFANode * getEmbeddedNode(unsigned int index) const
Definition: EFAEdge.C:328
unsigned int numEdges() const
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:407
bool _use_crack_growth_increment
Definition: XFEM.h:181
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:1713
EFAFragment2D * getFragment(unsigned int frag_id) const
unsigned int getLocalNodeIndex(EFANode *node) const
Definition: EFAElement.C:109
bool markCutEdgesByState(Real time)
Definition: XFEM.C:548
void addElemNodeIntersection(unsigned int elemid, unsigned int nodeid)
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:274
XFEM_QRULE
Definition: XFEM.h:36
ElementPairLocator::ElementPairList _sibling_elems
Definition: XFEM.h:188
void addElemEdgeIntersection(unsigned int elemid, unsigned int edgeid, double position)
std::vector< const GeometricCutUserObject * > _geometric_cuts
Definition: XFEM.h:184
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:1671
Real _crack_growth_increment
Definition: XFEM.h:182
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const =0
void getFragmentFaces(const Elem *elem, std::vector< std::vector< Point >> &frag_faces, bool displaced_mesh=false) const
Definition: XFEM.C:1412
void getCrackTipOrigin(std::map< unsigned int, const Elem * > &elem_id_crack_tip, std::vector< Point > &crack_front_points)
Definition: XFEM.C:53
virtual bool isPartial() const =0
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
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:1348
Real getPhysicalVolumeFraction() const
Definition: XFEMCutElem.C:31
const std::vector< EFAElement * > & getChildElements()
std::set< const Elem * > _crack_tip_elems
Definition: XFEM.h:187
XFEM(const InputParameters &params)
Definition: XFEM.C:29
EFAFragment3D * getFragment(unsigned int frag_id) const
bool isFacePhantom(unsigned int face_id) const
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
bool isThirdInteriorFace(unsigned int face_id) const
virtual bool isFinalCut() const
Definition: EFAElement2D.C:789
virtual bool getXFEMWeights(MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points)
Definition: XFEM.C:1496
virtual void getMasterInfo(EFANode *node, std::vector< EFANode * > &master_nodes, std::vector< double > &master_weights) const =0
void addGeometricCut(const GeometricCutUserObject *geometric_cut)
Definition: XFEM.C:47
bool _has_secondary_cut
Definition: XFEM.h:177
EFAEdge * getEdge(unsigned int edge_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
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:217
EFANode * getTipEmbeddedNode() const
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
void addFragFaceIntersection(unsigned int ElemID, unsigned int FragFaceID, std::vector< unsigned int > FragFaceEdgeID, std::vector< double > position)
void updateTopology(bool mergeUncutVirtualEdges=true)
const std::vector< EFAElement * > & getParentElements()
virtual void getXFEMqRuleOnLine(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
Definition: XFEM.C:1532
bool initCutIntersectionEdge(Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
Definition: XFEM.C:905
unsigned int getCrackTipSplitElementID() const
Definition: EFAElement2D.C:592
void storeCrackTipOriginAndDirection()
Definition: XFEM.C:151
bool isSecondaryInteriorEdge(unsigned int edge_id) const
const std::vector< EFANode * > & getNewNodes()
void addElemFaceIntersection(unsigned int elemid, unsigned int faceid, std::vector< unsigned int > edgeid, std::vector< double > position)
bool markCuts(Real time)
Definition: XFEM.C:308
virtual bool update(Real time, NonlinearSystemBase &nl, AuxiliarySystem &aux)
Method to update the mesh due to modified cut planes.
Definition: XFEM.C:190
virtual bool isFinalCut() const
Definition: EFAElement3D.C:777
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:1637
bool addFragEdgeIntersection(unsigned int elemid, unsigned int frag_edge_id, double position)
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:1762
EFAFace * getFragmentFace(unsigned int frag_id, unsigned int face_id) const
Xfem::XFEM_QRULE _XFEM_qrule
Definition: XFEM.h:179
bool isEdgePhantom(unsigned int edge_id) 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:1603
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=NULL) const =0
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
Definition: XFEM.C:1296
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:191
XFEM_CUTPLANE_QUANTITY
Definition: XFEM.h:26
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::set< const Elem * > _state_marked_frags
Definition: XFEM.h:196
void restoreFragmentInfo(EFAElement *const elem, const EFAElement *const from_elem)
virtual void getXFEMqRuleOnSurface(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
Definition: XFEM.C:1559
bool isElemCut(const Elem *elem, XFEMCutElem *&xfce) const
Definition: XFEM.C:1388
bool hasIntersection() const
Definition: EFAEdge.C:196
virtual void computePhysicalVolumeFraction()=0
bool cutMeshWithEFA(NonlinearSystemBase &nl, AuxiliarySystem &aux)
Definition: XFEM.C:924
bool markCutFacesByGeometry(Real time)
Definition: XFEM.C:823
unsigned int id() const
Definition: EFANode.C:34
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:101
virtual void initSolution(NonlinearSystemBase &nl, AuxiliarySystem &aux)
Initialize the solution on newly created nodes.
Definition: XFEM.C:238
ElementPairLocator::ElementPairList _sibling_displaced_elems
Definition: XFEM.h:189
virtual unsigned int numFragments() const
Definition: EFAElement3D.C:263
Xfem::XFEM_QRULE & getXFEMQRule()
Definition: XFEM.C:1472
bool isPartialOverlap(const EFAEdge &other) const
Definition: EFAEdge.C:55
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:1736
void getFragmentEdges(const Elem *elem, EFAElement2D *CEMElem, std::vector< std::vector< Point >> &frag_edges) const
Definition: XFEM.C:1429
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:200
EFAElement * add2DElement(std::vector< unsigned int > quad, unsigned int id)
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=NULL) const =0
void setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment)
Definition: XFEM.C:1489
void addStateMarkedElem(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:101
EFAElement * getElemByID(unsigned int id)
bool markCutEdgesByGeometry(Real time)
Definition: XFEM.C:325
virtual const EFAElement * getEFAElement() const =0
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:208
bool isElemAtCrackTip(const Elem *elem) const
Definition: XFEM.C:1382
Real getPhysicalVolumeFraction(const Elem *elem) const
Get the volume fraction of an element that is physical.
Definition: XFEM.C:1328
void clearStateMarkedElems()
Definition: XFEM.C:143
void getWeightMultipliers(MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:37
virtual bool cutElementByCrackGrowthIncrement(const Elem *elem, std::vector< CutEdgeForCrackGrowthIncr > &cut_edges, Real time)
unsigned int numFaces() const
unsigned int numEdges() const
void setXFEMQRule(std::string &xfem_qrule)
Definition: XFEM.C:1478
void buildEFAMesh()
Definition: XFEM.C:265
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Definition: XFEM.h:197
EFAElement * add3DElement(std::vector< unsigned int > quad, unsigned int id)
EFANode * getNode(unsigned int node_id) const
Definition: EFAFace.C:97
virtual void getXFEMIntersectionInfo(const Elem *elem, unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, bool displaced_mesh=false) const
Definition: XFEM.C:1513
unsigned int getTipEdgeID() const
const std::set< EFAElement * > & getCrackTipElements()
std::map< const Elem *, RealVectorValue > _state_marked_elems
Definition: XFEM.h:195
bool markCutFacesByState()
Definition: XFEM.C:897
void addStateMarkedFrag(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:127