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