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