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