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