33 #include "libmesh/mesh_communication.h" 34 #include "libmesh/partitioner.h" 39 #ifndef LIBMESH_ENABLE_UNIQUE_ID 40 mooseError(
"MOOSE requires unique ids to be enabled in libmesh (configure with " 41 "--enable-unique-id) to use XFEM!");
48 for (std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
_cut_elem_map.begin();
66 std::vector<Point> & crack_front_points)
68 elem_id_crack_tip.clear();
69 crack_front_points.clear();
72 unsigned int crack_tip_index = 0;
75 std::map<unsigned int, const Elem *> elem_id_map;
78 for (std::map<
const Elem *, std::vector<Point>>::iterator mit1 =
83 unsigned int elem_id = mit1->first->id();
84 if (elem_id == std::numeric_limits<unsigned int>::max())
86 elem_id_map[m] = mit1->first;
90 elem_id_map[elem_id] = mit1->first;
93 for (std::map<unsigned int, const Elem *>::iterator mit1 = elem_id_map.begin();
94 mit1 != elem_id_map.end();
97 const Elem * elem = mit1->second;
98 std::map<const Elem *, std::vector<Point>>::iterator mit2 =
102 elem_id_crack_tip[crack_tip_index] = mit2->first;
103 crack_front_points[crack_tip_index] =
113 Elem * elem =
_mesh->elem_ptr(elem_id);
114 std::map<const Elem *, RealVectorValue>::iterator mit;
117 mooseError(
" ERROR: element ", elem->id(),
" already marked for crack growth.");
125 Elem * elem =
_mesh->elem_ptr(elem_id);
126 std::map<const Elem *, unsigned int>::iterator mit;
129 mooseError(
" ERROR: side of element ", elem->id(),
" already marked for crack initiation.");
138 Elem * elem =
_mesh->elem_ptr(elem_id);
139 std::set<const Elem *>::iterator mit;
143 " ERROR: element ", elem->id(),
" already marked for fragment-secondary crack initiation.");
159 const unsigned int interface_id)
161 Elem * elem =
_mesh->elem_ptr(elem_id);
169 const unsigned int interface_id)
171 Elem * elem =
_mesh->elem_ptr(elem_id);
188 std::set<EFAElement *>::iterator sit;
189 for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
191 if (
_mesh->mesh_dimension() == 2)
197 Point origin(0, 0, 0);
198 Point direction(0, 0, 0);
200 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
212 std::vector<Point> tip_data;
213 tip_data.push_back(origin);
214 tip_data.push_back(direction);
215 const Elem * elem =
_mesh->elem_ptr((*sit)->id());
217 std::pair<
const Elem *, std::vector<Point>>(elem, tip_data));
225 bool mesh_changed =
false;
234 _mesh->update_parallel_id_counts();
235 MeshCommunication().make_elems_parallel_consistent(*
_mesh);
236 MeshCommunication().make_nodes_parallel_consistent(*
_mesh);
239 _mesh->allow_renumbering(
false);
240 _mesh->skip_partitioning(
true);
241 _mesh->prepare_for_use();
261 const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
265 mooseError(
"Use of XFEM with distributed mesh is not yet supported");
267 bool mesh_changed =
false;
286 _mesh->allow_renumbering(
false);
287 _mesh->skip_partitioning(
true);
288 _mesh->prepare_for_use();
309 mooseError(
"XFEM does not currently support multiple nonlinear systems");
311 nls[0]->serializeSolution();
324 current_solution.
close();
325 old_solution.
close();
326 older_solution.
close();
327 current_aux_solution.
close();
328 old_aux_solution.
close();
329 older_aux_solution.
close();
341 for (
auto & elem :
_mesh->element_ptr_range())
343 std::vector<unsigned int> quad;
344 for (
unsigned int i = 0; i < elem->n_nodes(); ++i)
345 quad.push_back(elem->node_id(i));
347 if (
_mesh->mesh_dimension() == 2)
349 else if (
_mesh->mesh_dimension() == 3)
356 for (
auto & elem :
_mesh->element_ptr_range())
358 std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
_cut_elem_map.find(elem->unique_id());
377 bool marked_sides =
false;
378 if (
_mesh->mesh_dimension() == 2)
383 else if (
_mesh->mesh_dimension() == 3)
394 bool marked_edges =
false;
395 bool marked_nodes =
false;
399 for (
const auto & gmei : gme.second)
403 for (
unsigned int i = 0; i < gmei._elem_cut_edges.size(); ++i)
406 gmei._elem_cut_edges[i]._host_side_id))
409 gmei._elem_cut_edges[i]._host_side_id,
410 gmei._elem_cut_edges[i]._distance);
415 for (
unsigned int i = 0; i < gmei._elem_cut_nodes.size(); ++i)
421 for (
unsigned int i = 0; i < gmei._frag_cut_edges.size();
425 gmei._frag_cut_edges[i]._host_side_id))
428 gmei._frag_cut_edges[i]._host_side_id,
429 gmei._frag_cut_edges[i]._distance))
440 return marked_edges || marked_nodes;
448 Point crack_tip_origin,
449 Point crack_tip_direction,
450 Real & distance_keep,
451 unsigned int & edge_id_keep,
454 std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
455 Point edge1(0.0, 0.0, 0.0);
456 Point edge2(0.0, 0.0, 0.0);
457 Point left_angle(0.0, 0.0, 0.0);
458 Point right_angle(0.0, 0.0, 0.0);
459 Point left_angle_normal(0.0, 0.0, 0.0);
460 Point right_angle_normal(0.0, 0.0, 0.0);
461 Point crack_direction_normal(0.0, 0.0, 0.0);
462 Point edge1_to_tip(0.0, 0.0, 0.0);
463 Point edge2_to_tip(0.0, 0.0, 0.0);
464 Point edge1_to_tip_normal(0.0, 0.0, 0.0);
465 Point edge2_to_tip_normal(0.0, 0.0, 0.0);
467 Real cos_45 = std::cos(45.0 / 180.0 * 3.14159);
468 Real sin_45 = std::sin(45.0 / 180.0 * 3.14159);
470 left_angle(0) = cos_45 * crack_tip_direction(0) - sin_45 * crack_tip_direction(1);
471 left_angle(1) = sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
473 right_angle(0) = cos_45 * crack_tip_direction(0) + sin_45 * crack_tip_direction(1);
474 right_angle(1) = -sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
476 left_angle_normal(0) = -left_angle(1);
477 left_angle_normal(1) = left_angle(0);
479 right_angle_normal(0) = -right_angle(1);
480 right_angle_normal(1) = right_angle(0);
482 crack_direction_normal(0) = -crack_tip_direction(1);
483 crack_direction_normal(1) = crack_tip_direction(0);
485 Real angle_min = 0.0;
487 unsigned int nsides = CEMElem->
numEdges();
489 for (
unsigned int i = 0; i < nsides; ++i)
496 edge1_to_tip = (edge_ends[0] * 0.95 + edge_ends[1] * 0.05) - crack_tip_origin;
497 edge2_to_tip = (edge_ends[0] * 0.05 + edge_ends[1] * 0.95) - crack_tip_origin;
499 edge1_to_tip /=
pow(edge1_to_tip.norm_sq(), 0.5);
500 edge2_to_tip /=
pow(edge2_to_tip.norm_sq(), 0.5);
502 edge1_to_tip_normal(0) = -edge1_to_tip(1);
503 edge1_to_tip_normal(1) = edge1_to_tip(0);
505 edge2_to_tip_normal(0) = -edge2_to_tip(1);
506 edge2_to_tip_normal(1) = edge2_to_tip(0);
508 Real angle_edge1_normal = edge1_to_tip_normal * normal;
509 Real angle_edge2_normal = edge2_to_tip_normal * normal;
512 (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
515 distance_keep = 0.05;
516 normal_keep = edge1_to_tip_normal;
517 angle_min = angle_edge1_normal;
520 (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
523 distance_keep = 0.95;
524 normal_keep = edge2_to_tip_normal;
525 angle_min = angle_edge2_normal;
529 crack_tip_origin, left_angle_normal, edge_ends[0], edge_ends[1],
distance) &&
533 (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
537 normal_keep = left_angle_normal;
538 angle_min = left_angle_normal * normal;
542 crack_tip_origin, right_angle_normal, edge_ends[0], edge_ends[1],
distance) &&
546 (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
550 normal_keep = right_angle_normal;
551 angle_min = right_angle_normal * normal;
555 crack_direction_normal,
562 (crack_tip_direction * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
566 normal_keep = crack_direction_normal;
567 angle_min = crack_direction_normal * normal;
574 if ((distance_keep - 0.05) < 0.0)
576 distance_keep = 0.05;
578 else if ((distance_keep - 0.95) > 0.0)
580 distance_keep = 0.95;
587 bool marked_edges =
false;
588 for (std::map<const Elem *, RealVectorValue>::iterator pmeit =
_state_marked_elems.begin();
592 const Elem * elem = pmeit->first;
597 if (volfrac_elem < 0.25)
605 unsigned int nsides = CEMElem->
numEdges();
606 unsigned int orig_cut_side_id = std::numeric_limits<unsigned int>::max();
607 Real orig_cut_distance = -1.0;
612 Point crack_tip_origin(0, 0, 0);
613 Point crack_tip_direction(0, 0, 0);
618 if (orig_cut_side_id < nsides)
620 orig_edge = CEMElem->
getEdge(orig_cut_side_id);
624 mooseError(
"element ", elem->id(),
" has no valid crack-tip edge");
627 std::map<const Elem *, std::vector<Point>>::iterator ecodm =
631 crack_tip_origin = (ecodm->second)[0];
632 crack_tip_direction = (ecodm->second)[1];
635 mooseError(
"element ", elem->id(),
" cannot find its crack tip origin and direction.");
639 std::map<const Elem *, unsigned int>::iterator mit1;
641 std::set<const Elem *>::iterator mit2;
646 orig_cut_side_id = mit1->second;
650 orig_cut_distance = 0.5;
652 orig_edge = CEMElem->
getEdge(orig_cut_side_id);
655 Point elem_center(0.0, 0.0, 0.0);
657 for (
unsigned int i = 0; i < nsides; ++i)
662 elem_center /= nsides * 2.0;
666 crack_tip_origin = edge_center;
667 crack_tip_direction = elem_center - edge_center;
668 crack_tip_direction /=
pow(crack_tip_direction.norm_sq(), 0.5);
678 " flagged for a secondary crack, but has ",
682 if (interior_edge_id.size() == 1)
683 orig_cut_side_id = interior_edge_id[0];
687 orig_cut_distance = 0.5;
691 Point elem_center(0.0, 0.0, 0.0);
694 for (
unsigned int i = 0; i < nsides_frag; ++i)
701 elem_center /= nsides_frag * 2.0;
705 crack_tip_origin = edge_center;
706 crack_tip_direction = elem_center - edge_center;
707 crack_tip_direction /=
pow(crack_tip_direction.norm_sq(), 0.5);
712 " flagged for state-based growth, but has no edge intersections");
715 Point cut_origin(0.0, 0.0, 0.0);
719 mooseError(
"element ", elem->id(),
" does not have valid orig_node");
722 std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
723 Point edge1(0.0, 0.0, 0.0);
724 Point edge2(0.0, 0.0, 0.0);
725 Point cut_edge_point(0.0, 0.0, 0.0);
726 bool find_compatible_direction =
false;
727 unsigned int edge_id_keep = 0;
728 Real distance_keep = 0.0;
729 Point normal_keep(0.0, 0.0, 0.0);
731 bool edge_cut =
false;
733 for (
unsigned int i = 0; i < nsides; ++i)
740 crack_tip_origin, normal, edge_ends[0], edge_ends[1],
distance) &&
746 normal_keep = normal;
753 Point between_two_cuts = (cut_edge_point - crack_tip_origin);
754 between_two_cuts /=
pow(between_two_cuts.norm_sq(), 0.5);
755 Real angle_between_two_cuts = between_two_cuts * crack_tip_direction;
757 if (angle_between_two_cuts > std::cos(45.0 / 180.0 * 3.14159))
758 find_compatible_direction =
true;
760 if (!find_compatible_direction && edge_cut)
777 Point growth_direction(0.0, 0.0, 0.0);
779 growth_direction(0) = -normal_keep(1);
780 growth_direction(1) = normal_keep(0);
782 if (growth_direction * crack_tip_direction < 1.0e-10)
783 growth_direction *= (-1.0);
785 Real x0 = crack_tip_origin(0);
786 Real y0 = crack_tip_origin(1);
792 for (
const auto & elem :
_mesh->element_ptr_range())
794 std::vector<CutEdgeForCrackGrowthIncr> elem_cut_edges;
804 for (
unsigned int i = 0; i < elem_cut_edges.size(); ++i)
807 elem_cut_edges[i]._host_side_id))
810 elem->id(), elem_cut_edges[i]._host_side_id, elem_cut_edges[i]._distance);
829 crack_tip_origin, normal, edge_ends[0], edge_ends[1],
distance) &&
851 bool marked_faces =
false;
855 for (
const auto & gmei : gme.second)
859 for (
unsigned int i = 0; i < gmei._elem_cut_faces.size(); ++i)
861 if (!EFAElem->
isFacePhantom(gmei._elem_cut_faces[i]._face_id))
864 gmei._elem_cut_faces[i]._face_id,
865 gmei._elem_cut_faces[i]._face_edge,
866 gmei._elem_cut_faces[i]._position);
871 for (
unsigned int i = 0; i < gmei._frag_cut_faces.size();
877 gmei._frag_cut_faces[i]._face_id,
878 gmei._frag_cut_faces[i]._face_edge,
879 gmei._frag_cut_faces[i]._position);
892 bool marked_faces =
false;
899 Point cut_origin,
RealVectorValue cut_normal, Point & edge_p1, Point & edge_p2, Real & dist)
902 bool does_intersect =
false;
903 Point origin2p1 = edge_p1 - cut_origin;
904 Real plane2p1 = cut_normal(0) * origin2p1(0) + cut_normal(1) * origin2p1(1);
905 Point origin2p2 = edge_p2 - cut_origin;
906 Real plane2p2 = cut_normal(0) * origin2p2(0) + cut_normal(1) * origin2p2(1);
908 if (plane2p1 * plane2p2 < 0.0)
910 dist = -plane2p1 / (plane2p2 - plane2p1);
911 does_intersect =
true;
913 return does_intersect;
919 bool mesh_changed =
false;
921 std::set<Node *> nodes_to_delete;
922 std::set<Node *> nodes_to_delete_displaced;
923 std::set<unsigned int> cutelems_to_delete;
924 unsigned int deleted_elem_count = 0;
925 std::vector<std::string> healed_geometric_cuts;
934 Elem * elem1 =
const_cast<Elem *
>(it.first);
935 Elem * elem2 =
const_cast<Elem *
>(it.second);
937 std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
943 cutelems_to_delete.insert(elem1->unique_id());
945 for (
unsigned int in = 0; in < elem1->n_nodes(); ++in)
947 Node * e1node = elem1->node_ptr(in);
948 Node * e2node = elem2->node_ptr(in);
952 elem1->set_node(in) = e2node;
953 nodes_to_delete.insert(e1node);
955 else if (e1node != e2node)
956 nodes_to_delete.insert(e2node);
960 mooseError(
"Could not find XFEMCutElem for element to be kept in healing");
965 std::vector<const Elem *> healed_elems = {elem1, elem2};
971 for (
auto e : healed_elems)
972 if (elem1->processor_id() ==
_mesh->processor_id() &&
973 e->processor_id() ==
_mesh->processor_id())
981 if (parent_gcsid == gcsid)
990 std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
996 for (
unsigned int in = 0; in < elem1_displaced->n_nodes(); ++in)
998 Node * e1node_displaced = elem1_displaced->node_ptr(in);
999 Node * e2node_displaced = elem2_displaced->node_ptr(in);
1001 e1node_displaced != e2node_displaced)
1003 elem1_displaced->set_node(in) = e2node_displaced;
1004 nodes_to_delete_displaced.insert(e1node_displaced);
1006 else if (e1node_displaced != e2node_displaced)
1007 nodes_to_delete_displaced.insert(e2node_displaced);
1011 mooseError(
"Could not find XFEMCutElem for element to be kept in healing");
1013 elem2_displaced->nullify_neighbors();
1022 cutelems_to_delete.insert(elem2->unique_id());
1023 elem2->nullify_neighbors();
1024 _mesh->get_boundary_info().remove(elem2);
1025 unsigned int deleted_elem_id = elem2->id();
1026 _mesh->delete_elem(elem2);
1029 if (deleted_elem_count == 0)
1031 _console <<
"XFEM healing deleted element: " << deleted_elem_id << std::endl;
1033 ++deleted_elem_count;
1034 mesh_changed =
true;
1039 for (
auto & sit : nodes_to_delete)
1041 Node * node_to_delete = sit;
1042 dof_id_type deleted_node_id = node_to_delete->id();
1043 _mesh->get_boundary_info().remove(node_to_delete);
1044 _mesh->delete_node(node_to_delete);
1046 _console <<
"XFEM healing deleted node: " << deleted_node_id << std::endl;
1051 for (
auto & sit : nodes_to_delete_displaced)
1053 Node * node_to_delete_displaced = sit;
1054 _displaced_mesh->get_boundary_info().remove(node_to_delete_displaced);
1059 for (
auto & ced : cutelems_to_delete)
1087 _console <<
"\nXFEM mesh healing complete\n";
1088 _console <<
"Names of healed geometric cut objects: ";
1089 for (
auto geomcut : healed_geometric_cuts)
1092 _console <<
"# deleted nodes: " << nodes_to_delete.size() <<
"\n";
1093 _console <<
"# deleted elements: " << deleted_elem_count <<
"\n";
1097 return mesh_changed;
1104 if (nls.size() != 1)
1105 mooseError(
"XFEM does not currently support multiple nonlinear systems");
1107 std::map<unsigned int, Node *> efa_id_to_new_node;
1108 std::map<unsigned int, Node *> efa_id_to_new_node2;
1109 std::map<unsigned int, Elem *> efa_id_to_new_elem;
1122 _console <<
"\nXFEM Element fragment algorithm mesh prior to cutting:\n";
1131 _console <<
"\nXFEM Element fragment algorithm mesh after cutting:\n";
1140 bool mesh_changed = (new_nodes.size() + new_elements.size() + delete_elements.size() > 0);
1145 nls[0]->serializeSolution();
1157 std::map<Node *, Node *> new_nodes_to_parents;
1160 for (
unsigned int i = 0; i < new_nodes.size(); ++i)
1162 unsigned int new_node_id = new_nodes[i]->id();
1163 unsigned int parent_id = new_nodes[i]->parent()->id();
1165 Node * parent_node =
_mesh->node_ptr(parent_id);
1166 Node * new_node = Node::build(*parent_node,
_mesh->n_nodes()).release();
1167 _mesh->add_node(new_node);
1169 new_nodes_to_parents[new_node] = parent_node;
1171 new_node->set_n_systems(parent_node->n_systems());
1172 efa_id_to_new_node.insert(std::make_pair(new_node_id, new_node));
1174 _console <<
"XFEM added new node: " << new_node->id() << std::endl;
1178 Node * new_node2 = Node::build(*parent_node2,
_displaced_mesh->n_nodes()).release();
1181 new_node2->set_n_systems(parent_node2->n_systems());
1182 efa_id_to_new_node2.insert(std::make_pair(new_node_id, new_node2));
1187 std::map<unsigned int, std::vector<const Elem *>> temporary_parent_children_map;
1189 std::vector<boundary_id_type> parent_boundary_ids;
1191 for (
unsigned int i = 0; i < new_elements.size(); ++i)
1193 unsigned int parent_id = new_elements[i]->getParent()->id();
1194 unsigned int efa_child_id = new_elements[i]->id();
1196 Elem * parent_elem =
_mesh->elem_ptr(parent_id);
1197 Elem * libmesh_elem = Elem::build(parent_elem->type()).release();
1203 if (parent_elem == it.first)
1204 it.first = libmesh_elem;
1205 else if (parent_elem == it.second)
1206 it.second = libmesh_elem;
1211 if (new_elements[i]->getParent()->numChildren() > 1)
1212 temporary_parent_children_map[parent_elem->id()].push_back(libmesh_elem);
1214 Elem * parent_elem2 =
nullptr;
1215 Elem * libmesh_elem2 =
nullptr;
1219 libmesh_elem2 = Elem::build(parent_elem2->type()).release();
1225 if (parent_elem2 == it.first)
1226 it.first = libmesh_elem2;
1227 else if (parent_elem2 == it.second)
1228 it.second = libmesh_elem2;
1233 for (
unsigned int j = 0;
j < new_elements[i]->numNodes(); ++
j)
1235 unsigned int node_id = new_elements[i]->getNode(
j)->id();
1236 Node * libmesh_node;
1238 std::map<unsigned int, Node *>::iterator nit = efa_id_to_new_node.find(node_id);
1239 if (nit != efa_id_to_new_node.end())
1240 libmesh_node = nit->second;
1242 libmesh_node =
_mesh->node_ptr(node_id);
1244 if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
1245 libmesh_node->processor_id() = parent_elem->processor_id();
1247 libmesh_elem->set_node(
j) = libmesh_node;
1250 if (parent_elem->is_semilocal(
_mesh->processor_id()))
1252 Node * solution_node = libmesh_node;
1253 if (new_nodes_to_parents.find(libmesh_node) != new_nodes_to_parents.end())
1254 solution_node = new_nodes_to_parents[libmesh_node];
1257 (libmesh_node->processor_id() ==
_mesh->processor_id()))
1270 current_aux_solution,
1272 older_aux_solution);
1276 Node * parent_node = parent_elem->node_ptr(
j);
1277 _mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
1278 _mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
1282 std::map<unsigned int, Node *>::iterator nit2 = efa_id_to_new_node2.find(node_id);
1283 if (nit2 != efa_id_to_new_node2.end())
1284 libmesh_node = nit2->second;
1288 if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
1289 libmesh_node->processor_id() = parent_elem2->processor_id();
1291 libmesh_elem2->set_node(
j) = libmesh_node;
1293 parent_node = parent_elem2->node_ptr(
j);
1294 _displaced_mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
1295 _displaced_mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
1299 libmesh_elem->set_p_level(parent_elem->p_level());
1300 libmesh_elem->set_p_refinement_flag(parent_elem->p_refinement_flag());
1301 _mesh->add_elem(libmesh_elem);
1302 libmesh_elem->set_n_systems(parent_elem->n_systems());
1303 libmesh_elem->subdomain_id() = parent_elem->subdomain_id();
1304 libmesh_elem->processor_id() = parent_elem->processor_id();
1308 std::map<const Elem *, std::vector<Point>>::iterator mit =
1318 _console <<
"XFEM added new element: " << libmesh_elem->id() << std::endl;
1321 if (
_mesh->mesh_dimension() == 2)
1324 if (!new_efa_elem2d)
1325 mooseError(
"EFAelem is not of EFAelement2D type");
1329 libmesh_elem->n_sides());
1331 else if (
_mesh->mesh_dimension() == 3)
1334 if (!new_efa_elem3d)
1335 mooseError(
"EFAelem is not of EFAelement3D type");
1339 libmesh_elem->n_sides());
1341 _cut_elem_map.insert(std::pair<unique_id_type, XFEMCutElem *>(libmesh_elem->unique_id(), xfce));
1342 efa_id_to_new_elem.insert(std::make_pair(efa_child_id, libmesh_elem));
1346 libmesh_elem2->set_p_level(parent_elem2->p_level());
1347 libmesh_elem2->set_p_refinement_flag(parent_elem2->p_refinement_flag());
1349 libmesh_elem2->set_n_systems(parent_elem2->n_systems());
1350 libmesh_elem2->subdomain_id() = parent_elem2->subdomain_id();
1351 libmesh_elem2->processor_id() = parent_elem2->processor_id();
1354 unsigned int n_sides = parent_elem->n_sides();
1355 for (
unsigned int side = 0; side < n_sides; ++side)
1357 _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
1358 _mesh->get_boundary_info().add_side(libmesh_elem, side, parent_boundary_ids);
1362 n_sides = parent_elem2->n_sides();
1363 for (
unsigned int side = 0; side < n_sides; ++side)
1365 _displaced_mesh->get_boundary_info().boundary_ids(parent_elem2, side, parent_boundary_ids);
1366 _displaced_mesh->get_boundary_info().add_side(libmesh_elem2, side, parent_boundary_ids);
1370 unsigned int n_edges = parent_elem->n_edges();
1373 _mesh->get_boundary_info().edge_boundary_ids(parent_elem,
edge, parent_boundary_ids);
1374 _mesh->get_boundary_info().add_edge(libmesh_elem,
edge, parent_boundary_ids);
1378 n_edges = parent_elem2->n_edges();
1382 parent_elem2,
edge, parent_boundary_ids);
1388 if (parent_elem->processor_id() ==
_mesh->processor_id())
1390 if (
_material_data[0]->getMaterialPropertyStorage().hasStatefulProperties())
1394 for (
unsigned int side = 0; side < parent_elem->n_sides(); ++side)
1396 _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
1397 std::vector<boundary_id_type>::iterator it_bd = parent_boundary_ids.begin();
1398 for (; it_bd != parent_boundary_ids.end(); ++it_bd)
1418 if (cei.
match(old_cei.second))
1422 _console <<
"XFEM set material properties for element: " << libmesh_elem->id()
1440 current_aux_solution,
1442 older_aux_solution);
1447 for (std::size_t i = 0; i < delete_elements.size(); ++i)
1449 Elem * elem_to_delete =
_mesh->elem_ptr(delete_elements[i]->
id());
1452 std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
1456 delete cemit->second;
1464 elem_to_delete->nullify_neighbors();
1465 _mesh->get_boundary_info().remove(elem_to_delete);
1466 unsigned int deleted_elem_id = elem_to_delete->id();
1467 _mesh->delete_elem(elem_to_delete);
1469 _console <<
"XFEM deleted element: " << deleted_elem_id << std::endl;
1473 Elem * elem_to_delete2 =
_displaced_mesh->elem_ptr(delete_elements[i]->
id());
1474 elem_to_delete2->nullify_neighbors();
1480 for (std::map<
unsigned int, std::vector<const Elem *>>::iterator it =
1481 temporary_parent_children_map.begin();
1482 it != temporary_parent_children_map.end();
1485 std::vector<const Elem *> & sibling_elem_vec = it->second;
1492 if (it->first == elem_id)
1494 std::make_pair(sibling_elem_vec[0], sibling_elem_vec[1]));
1507 std::make_pair(elem, elem_pair));
1513 temporary_parent_children_map.clear();
1521 std::set<EFAElement *>::const_iterator sit;
1522 for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
1524 unsigned int eid = (*sit)->id();
1525 Elem * crack_tip_elem;
1526 std::map<unsigned int, Elem *>::iterator eit = efa_id_to_new_elem.find(eid);
1527 if (eit != efa_id_to_new_elem.end())
1528 crack_tip_elem = eit->second;
1530 crack_tip_elem =
_mesh->elem_ptr(eid);
1539 if ((*sit)->getParent() !=
nullptr)
1541 if (
_mesh->mesh_dimension() == 2)
1545 mooseError(
"EFAelem is not of EFAelement2D type");
1547 for (
unsigned int edge_id = 0; edge_id < efa_elem2d->
numEdges(); ++edge_id)
1549 for (
unsigned int en_iter = 0; en_iter < efa_elem2d->
numEdgeNeighbors(edge_id);
1553 if (edge_neighbor !=
nullptr && edge_neighbor->
id() == mie)
1558 else if (
_mesh->mesh_dimension() == 3)
1562 mooseError(
"EFAelem is not of EFAelement3D type");
1564 for (
unsigned int face_id = 0; face_id < efa_elem3d->
numFaces(); ++face_id)
1566 for (
unsigned int fn_iter = 0; fn_iter < efa_elem3d->
numFaceNeighbors(face_id);
1570 if (face_neighbor !=
nullptr && face_neighbor->
id() == mie)
1583 _console <<
"\nXFEM mesh cutting with element fragment algorithm complete\n";
1584 _console <<
"# new nodes: " << new_nodes.size() <<
"\n";
1585 _console <<
"# new elements: " << new_elements.size() <<
"\n";
1586 _console <<
"# deleted elements: " << delete_elements.size() <<
"\n";
1592 return mesh_changed;
1599 MeshBase * displaced_mesh)
const 1601 Point node_coor(0.0, 0.0, 0.0);
1602 std::vector<EFANode *> master_nodes;
1603 std::vector<Point> master_points;
1604 std::vector<double> master_weights;
1606 CEMElem->
getMasterInfo(CEMnode, master_nodes, master_weights);
1607 for (std::size_t i = 0; i < master_nodes.size(); ++i)
1612 const Node * node = elem->node_ptr(local_node_id);
1614 node = displaced_mesh->node_ptr(node->id());
1615 Point node_p((*node)(0), (*node)(1), (*node)(2));
1616 master_points.push_back(node_p);
1619 mooseError(
"master nodes must be permanent");
1621 for (std::size_t i = 0; i < master_nodes.size(); ++i)
1622 node_coor += master_weights[i] * master_points[i];
1630 Real phys_volfrac = 1.0;
1631 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1644 return phys_volfrac;
1650 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1668 unsigned int plane_id)
const 1671 Point planedata(0.0, 0.0, 0.0);
1672 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1680 if ((
unsigned int)quantity < 3)
1682 unsigned int index = (
unsigned int)quantity;
1684 comp = planedata(index);
1686 else if ((
unsigned int)quantity < 6)
1688 unsigned int index = (
unsigned int)quantity - 3;
1690 comp = planedata(index);
1693 mooseError(
"In get_cut_plane index out of range");
1730 std::vector<std::vector<Point>> & frag_faces,
1731 bool displaced_mesh)
const 1733 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1752 mooseError(
"EFAelem is not of EFAelement2D type");
1764 mooseError(
"EFAelem is not of EFAelement3D type");
1772 std::vector<std::vector<Point>> & frag_edges)
const 1779 mooseError(
"element ", elem->id(),
" has more than one fragment at this point");
1782 std::vector<Point> p_line(2, Point(0.0, 0.0, 0.0));
1785 frag_edges.push_back(p_line);
1793 std::vector<std::vector<Point>> & frag_faces)
const 1800 mooseError(
"element ", elem->id(),
" has more than one fragment at this point");
1804 std::vector<Point> p_line(num_face_nodes, Point(0.0, 0.0, 0.0));
1805 for (
unsigned int j = 0;
j < num_face_nodes; ++
j)
1807 frag_faces.push_back(p_line);
1821 if (xfem_qrule ==
"volfrac")
1823 else if (xfem_qrule ==
"moment_fitting")
1825 else if (xfem_qrule ==
"direct")
1848 bool have_weights =
false;
1852 mooseAssert(xfce !=
nullptr,
"Must have valid XFEMCutElem object here");
1854 have_weights =
true;
1856 return have_weights;
1866 bool have_weights =
false;
1870 mooseAssert(xfce !=
nullptr,
"Must have valid XFEMCutElem object here");
1872 have_weights =
true;
1874 return have_weights;
1879 unsigned int plane_id,
1881 std::vector<Point> & intersectionPoints,
1882 bool displaced_mesh)
const 1884 std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1898 std::vector<Point> & quad_pts,
1899 std::vector<Real> & quad_wts)
const 1901 Point p1 = intersection_points[0];
1902 Point p2 = intersection_points[1];
1905 std::size_t num_qpoints = 2;
1911 quad_wts.resize(num_qpoints);
1912 quad_pts.resize(num_qpoints);
1916 quad_wts[0] = 1.0 * integ_jacobian;
1917 quad_wts[1] = 1.0 * integ_jacobian;
1919 quad_pts[0] = (1.0 - xi0) / 2.0 * p1 + (1.0 + xi0) / 2.0 * p2;
1920 quad_pts[1] = (1.0 - xi1) / 2.0 * p1 + (1.0 + xi1) / 2.0 * p2;
1925 std::vector<Point> & quad_pts,
1926 std::vector<Real> & quad_wts)
const 1928 std::size_t nnd_pe = intersection_points.size();
1929 Point xcrd(0.0, 0.0, 0.0);
1930 for (std::size_t i = 0; i < nnd_pe; ++i)
1931 xcrd += intersection_points[i];
1934 quad_pts.resize(nnd_pe);
1935 quad_wts.resize(nnd_pe);
1939 for (std::size_t
j = 0;
j < nnd_pe; ++
j)
1941 std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
1942 std::vector<Point> subtrig_points(3, Point(0.0, 0.0, 0.0));
1944 int jplus1 =
j < nnd_pe - 1 ?
j + 1 : 0;
1945 subtrig_points[0] = xcrd;
1946 subtrig_points[1] = intersection_points[
j];
1947 subtrig_points[2] = intersection_points[jplus1];
1949 std::vector<std::vector<Real>> sg2;
1951 for (std::size_t l = 0; l < sg2.size(); ++l)
1954 std::vector<Real> tsg_line(3, 0.0);
1955 for (std::size_t
k = 0;
k < 3; ++
k)
1957 tsg_line[0] += shape[
k][2] * subtrig_points[
k](0);
1958 tsg_line[1] += shape[
k][2] * subtrig_points[
k](1);
1959 tsg_line[2] += shape[
k][2] * subtrig_points[
k](2);
1961 quad_pts[
j + l] = Point(tsg_line[0], tsg_line[1], tsg_line[2]);
1962 quad_wts[
j + l] = sg2[l][3] * jac;
1969 const Node * node_to_store_from,
1976 std::vector<dof_id_type> stored_solution_dofs =
getNodeSolutionDofs(node_to_store_from, sys);
1977 std::vector<Real> stored_solution_scratch;
1979 std::size_t stored_solution_size =
1981 stored_solution_scratch.reserve(stored_solution_size);
1985 for (
auto dof : stored_solution_dofs)
1986 stored_solution_scratch.push_back(current_solution(dof));
1990 for (
auto dof : stored_solution_dofs)
1991 stored_solution_scratch.push_back(old_solution(dof));
1993 for (
auto dof : stored_solution_dofs)
1994 stored_solution_scratch.push_back(older_solution(dof));
1997 if (stored_solution_scratch.size() > 0)
1998 stored_solution[node_to_store_to->unique_id()] = stored_solution_scratch;
2003 const Elem * elem_to_store_from,
2011 std::vector<Real> stored_solution_scratch;
2013 std::size_t stored_solution_size =
2015 stored_solution_scratch.reserve(stored_solution_size);
2019 for (
auto dof : stored_solution_dofs)
2020 stored_solution_scratch.push_back(current_solution(dof));
2024 for (
auto dof : stored_solution_dofs)
2025 stored_solution_scratch.push_back(old_solution(dof));
2027 for (
auto dof : stored_solution_dofs)
2028 stored_solution_scratch.push_back(older_solution(dof));
2031 if (stored_solution_scratch.size() > 0)
2032 stored_solution[elem_to_store_to->unique_id()] = stored_solution_scratch;
2037 const std::map<
unique_id_type, std::vector<Real>> & stored_solution,
2042 for (
auto & node :
_mesh->local_node_ptr_range())
2044 auto mit = stored_solution.find(node->unique_id());
2045 if (mit != stored_solution.end())
2047 const std::vector<Real> & stored_node_solution = mit->second;
2050 stored_solution_dofs,
2057 for (
auto & elem :
as_range(
_mesh->local_elements_begin(),
_mesh->local_elements_end()))
2059 auto mit = stored_solution.find(elem->unique_id());
2060 if (mit != stored_solution.end())
2062 const std::vector<Real> & stored_elem_solution = mit->second;
2065 stored_solution_dofs,
2075 const std::vector<dof_id_type> & stored_solution_dofs,
2082 const auto old_solution_offset = stored_solution_dofs.size();
2083 const auto older_solution_offset = old_solution_offset * 2;
2085 for (std::size_t i = 0; i < stored_solution_dofs.size(); ++i)
2087 current_solution.
set(stored_solution_dofs[i], stored_solution[i]);
2090 old_solution.
set(stored_solution_dofs[i], stored_solution[old_solution_offset + i]);
2091 older_solution.
set(stored_solution_dofs[i], stored_solution[older_solution_offset + i]);
2096 std::vector<dof_id_type>
2100 const std::vector<MooseVariableFEBase *> & vars = sys.
getVariables(0);
2101 std::vector<dof_id_type> solution_dofs;
2102 solution_dofs.reserve(vars.size());
2103 for (
auto var : vars)
2105 if (!var->isNodal())
2108 if (var_subdomains.empty() || var_subdomains.find(sid) != var_subdomains.end())
2110 unsigned int n_comp = elem->n_comp(sys.
number(), var->number());
2111 for (
unsigned int icomp = 0; icomp < n_comp; ++icomp)
2114 solution_dofs.push_back(elem_dof);
2119 return solution_dofs;
2122 std::vector<dof_id_type>
2126 const std::vector<MooseVariableFEBase *> & vars = sys.
getVariables(0);
2127 std::vector<dof_id_type> solution_dofs;
2128 solution_dofs.reserve(vars.size());
2129 for (
auto var : vars)
2134 std::set<SubdomainID> intersect;
2135 set_intersection(var_subdomains.begin(),
2136 var_subdomains.end(),
2139 std::inserter(intersect, intersect.begin()));
2140 if (var_subdomains.empty() || !intersect.empty())
2142 unsigned int n_comp = node->n_comp(sys.
number(), var->number());
2143 for (
unsigned int icomp = 0; icomp < n_comp; ++icomp)
2146 solution_dofs.push_back(node_dof);
2151 return solution_dofs;
2159 std::set<unsigned int> elems = gcmit.second;
2160 if (elems.find(elem->id()) != elems.end())
2171 const auto & elem_props = storage.
props(state).at(elem);
2172 auto & serialized_props =
_geom_cut_elems[elem]._elem_material_properties[state - 1];
2173 serialized_props.clear();
2174 for (
const auto & side_props_pair : elem_props)
2176 const auto side = side_props_pair.first;
2177 std::ostringstream oss;
2179 serialized_props[side].assign(oss.str());
2196 bool need_boundary_materials =
false;
2197 for (
unsigned int side = 0; side < child_elem->n_sides(); ++side)
2199 std::vector<boundary_id_type> elem_boundary_ids;
2200 _mesh->get_boundary_info().boundary_ids(child_elem, side, elem_boundary_ids);
2201 for (
auto bdid : elem_boundary_ids)
2203 need_boundary_materials =
true;
2207 if (need_boundary_materials)
2222 const auto & serialized_props = cached_props[state - 1];
2223 for (
const auto & [side, serialized_side_props] : serialized_props)
2225 std::istringstream iss;
2226 iss.str(serialized_side_props);
2239 const Elem * elem_from,
2240 std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei)
const 2243 mooseAssert(cached_cei.count(elem_from) > 0,
"XFEM: Unable to find cached material properties.");
2248 cei._elem_material_properties,
2252 bool need_boundary_materials =
false;
2253 for (
unsigned int side = 0; side < elem->n_sides(); ++side)
2255 std::vector<boundary_id_type> elem_boundary_ids;
2256 _mesh->get_boundary_info().boundary_ids(elem, side, elem_boundary_ids);
2257 for (
auto bdid : elem_boundary_ids)
2259 need_boundary_materials =
true;
2263 if (need_boundary_materials)
2266 cei._bnd_material_properties,
2272 const Elem * cut_elem,
2273 const Elem * parent_elem)
const 2276 parent_elem = cut_elem;
2286 for (
auto i : e0->node_index_range())
2288 return e0->node_ptr(i);
2289 mooseError(
"cannot find a physical node in the current element");
void getCrackTipOrigin(std::map< unsigned int, const Elem *> &elem_id_crack_tip, std::vector< Point > &crack_front_points)
EFAFragment3D * getFragment(unsigned int frag_id) const
void getFragmentFaces(const Elem *elem, std::vector< std::vector< Point >> &frag_faces, bool displaced_mesh=false) const
void correctCrackExtensionDirection(const Elem *elem, EFAElement2D *CEMElem, EFAEdge *orig_edge, Point normal, Point crack_tip_origin, Point crack_tip_direction, Real &distance_keep, unsigned int &edge_id_keep, Point &normal_keep)
bool _use_crack_growth_increment
unsigned int _debug_output_level
Controls amount of debugging output information 0: None 1: Summary 2: Details on modifications to mes...
void setSolutionForDOFs(const std::vector< Real > &stored_solution, const std::vector< dof_id_type > &stored_solution_dofs, NumericVector< Number > ¤t_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
Set the solution for a set of DOFs.
MaterialProperties & setProps(const Elem *elem, unsigned int side, const unsigned int state=0)
const GeometricCutUserObject * getGeometricCutForElem(const Elem *elem) const
Get the GeometricCutUserObject associated with an element.
const std::vector< MooseVariableFieldBase *> & getVariables(THREAD_ID tid)
unsigned int numEdgeNeighbors(unsigned int edge_id) const
bool markCutEdgesByState(Real time)
void addElemNodeIntersection(unsigned int elemid, unsigned int nodeid)
void storeMaterialPropertiesForElement(const Elem *parent_elem, const Elem *child_elem)
Helper function to store the material properties of a healed element.
std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo3D > > _geom_marked_elems_3d
Data structure for storing information about all 3D elements to be cut by geometry.
EFAElement3D * getFaceNeighbor(unsigned int face_id, unsigned int neighbor_id) const
EFANode * getEmbeddedNode(unsigned int index) const
CutSubdomainID getCutSubdomainID(const GeometricCutUserObject *gcuo, const Elem *cut_elem, const Elem *parent_elem=nullptr) const
Determine which cut subdomain the element belongs to relative to the cut.
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Data structure describing geometrically described cut through 3D element.
bool isFacePhantom(unsigned int face_id) const
void addElemEdgeIntersection(unsigned int elemid, unsigned int edgeid, double position)
std::vector< const GeometricCutUserObject * > _geometric_cuts
void setSolution(SystemBase &sys, const std::map< unique_id_type, std::vector< Real >> &stored_solution, NumericVector< Number > ¤t_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
Set the solution for all locally-owned nodes/elements that have stored values.
unsigned int getCrackTipSplitElementID() const
virtual bool update(Real time, const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override
std::map< unsigned int, ElementPairLocator::ElementPairList > _sibling_elems
Real _crack_growth_increment
bool isSecondaryInteriorEdge(unsigned int edge_id) const
bool shouldHealMesh() const
Should the elements cut by this cutting object be healed in the current time step?
virtual bool getXFEMFaceWeights(MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points, unsigned int side) override
void mooseError(Args &&... args)
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=nullptr) const =0
bool isPartialOverlap(const EFAEdge &other) const
bool isDistributedMesh() const
std::vector< dof_id_type > getNodeSolutionDofs(const Node *node, SystemBase &sys) const
Get a vector of the dof indices for all components of all variables associated with a node...
bool isPointInsidePhysicalDomain(const Elem *elem, const Point &point) const
Return true if the point is inside the element physical domain Note: if this element is not cut...
bool isSemiLocal(Node *const node) const
void updateEdgeNeighbors()
void getFragmentEdges(const Elem *elem, EFAElement2D *CEMElem, std::vector< std::vector< Point >> &frag_edges) const
bool match(const CutElemInfo &rhs)
bool cutMeshWithEFA(const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux)
unsigned int getTipEdgeID() const
unsigned int numEdges() const
EFAEdge * getEdge(unsigned int edge_id) const
std::map< unsigned int, ElementPairLocator::ElementPairList > _sibling_displaced_elems
Real getPhysicalVolumeFraction() const
Returns the volume fraction of the element fragment.
void clearGeomMarkedElems()
Clear out the list of elements to be marked for cutting.
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
EFAElement * add3DElement(const std::vector< unsigned int > &quad, unsigned int id)
bool needBoundaryMaterialOnSide(BoundaryID bnd_id, const THREAD_ID tid)
const std::vector< EFAElement * > & getChildElements()
std::set< const Elem * > _crack_tip_elems
bool healMesh()
Potentially heal the mesh by merging some of the pairs of partial elements cut by XFEM back into sing...
XFEM(const InputParameters ¶ms)
virtual bool isPartial() const =0
ElementFragmentAlgorithm _efa_mesh
const std::set< SubdomainID > & getNodeBlockIds(const Node &node) const
bool hasIntersection() const
virtual unsigned int numFragments() const
NumericVector< Number > & solutionOlder()
std::array< std::unordered_map< unsigned int, std::string >, 2 > CachedMaterialProperties
Convenient typedef for local storage of stateful material properties.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
IntRange< unsigned int > statefulIndexRange() const
std::vector< unsigned int > getInteriorEdgeID() const
unsigned int numEdges() const
void getFaceWeightMultipliers(MooseArray< Real > &face_weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points, unsigned int side)
bool isEdgePhantom(unsigned int edge_id) const
virtual unsigned int numFragments() const
Real distance(const Point &p)
std::set< const Elem * > _crack_tip_elems_to_be_healed
void addGeometricCut(GeometricCutUserObject *geometric_cut)
std::map< unique_id_type, std::vector< Real > > _cached_aux_solution
Data structure to store the auxiliary solution for nodes/elements affected by XFEM For each node/elem...
unsigned int numFaces() const
bool markCutEdgesByGeometry()
unsigned int CutSubdomainID
FEProblemBase * _fe_problem
std::vector< MaterialData *> _bnd_material_data
EFANode * getNode(unsigned int node_id) const
void updateTopology(bool mergeUncutVirtualEdges=true)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const std::vector< EFAElement * > & getParentElements()
virtual const EFAElement * getEFAElement() const =0
void updatePhysicalLinksAndFragments()
bool initCutIntersectionEdge(Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
void addGeomMarkedElem3D(const unsigned int elem_id, const Xfem::GeomMarkedElemInfo3D geom_info, const unsigned int interface_id)
Add information about a new cut to be performed on a specific 3d element.
virtual void execute(const ExecFlagType &exec_type)
void addElemFaceIntersection(unsigned int elemid, unsigned int faceid, const std::vector< unsigned int > &edgeid, const std::vector< double > &position)
void storeCrackTipOriginAndDirection()
const Node * pickFirstPhysicalNode(const Elem *e, const Elem *e0) const
Return the first node in the provided element that is found to be in the physical domain...
void loadMaterialPropertiesForElement(const Elem *elem, const Elem *elem_from, std::unordered_map< const Elem *, Xfem::CutElemInfo > &cached_cei) const
Helper function to store the material properties of a healed element.
const std::vector< EFANode * > & getNewNodes()
EFAFragment2D * getFragment(unsigned int frag_id) const
Real getCutPlane(const Elem *elem, const Xfem::XFEM_CUTPLANE_QUANTITY quantity, unsigned int plane_id) const
Get specified component of normal or origin for cut plane for a given element.
unsigned int numNodes() const
std::map< const GeometricCutUserObject *, unsigned int > _geom_marker_id_map
Data structure for storing the GeommetricCutUserObjects and their corresponding id.
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=nullptr) const =0
std::vector< MaterialData *> _material_data
void addGeomMarkedElem2D(const unsigned int elem_id, const Xfem::GeomMarkedElemInfo2D geom_info, const unsigned int interface_id)
Add information about a new cut to be performed on a specific 2d element.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
void storeSolutionForElement(const Elem *elem_to_store_to, const Elem *elem_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > ¤t_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
Store the solution in stored_solution for a given element.
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const =0
bool addFragEdgeIntersection(unsigned int elemid, unsigned int frag_edge_id, double position)
virtual bool isFinalCut() const
Xfem::XFEM_QRULE _XFEM_qrule
const QBase *const & qRule() const
void storeSolutionForNode(const Node *node_to_store_to, const Node *node_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > ¤t_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
Store the solution in stored_solution for a given node.
bool isThirdInteriorFace(unsigned int face_id) const
unsigned int number() const
virtual void getXFEMqRuleOnLine(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
Data structure describing geometrically described cut through 2D element.
void loadMaterialPropertiesForElementHelper(const Elem *elem, const Xfem::CachedMaterialProperties &cached_props, MaterialPropertyStorage &storage) const
Load the material properties.
void setDebugOutputLevel(unsigned int debug_output_level)
Controls amount of debugging information output.
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
virtual void getXFEMIntersectionInfo(const Elem *elem, unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, bool displaced_mesh=false) const
std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo2D > > _geom_marked_elems_2d
Data structure for storing information about all 2D elements to be cut by geometry.
std::map< unsigned int, std::set< unsigned int > > _geom_marker_id_elems
Data structure for storing the elements cut by specific geometric cutters.
void addFragFaceIntersection(unsigned int ElemID, unsigned int FragFaceID, const std::vector< unsigned int > &FragFaceEdgeID, const std::vector< double > &position)
virtual CutSubdomainID getCutSubdomainID(const Node *) const
Get CutSubdomainID telling which side the node belongs to relative to the cut.
bool hasStatefulProperties() const
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual void initSolution(const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override
std::set< const Elem * > _state_marked_frags
EFANode * getTipEmbeddedNode() const
void restoreFragmentInfo(EFAElement *const elem, const EFAElement *const from_elem)
void storeMaterialPropertiesForElementHelper(const Elem *elem, MaterialPropertyStorage &storage)
const PropsType & props(const unsigned int state=0) const
virtual void getMasterInfo(EFANode *node, std::vector< EFANode *> &master_nodes, std::vector< double > &master_weights) const =0
const ExecFlagType EXEC_XFEM_MARK
Exec flag used to execute MooseObjects while elements are being marked for cutting by XFEM...
EFAElement3D * getEFAElem3D(const Elem *elem)
Get the EFAElement3D object for a specified libMesh element.
EFAElement2D * getEdgeNeighbor(unsigned int edge_id, unsigned int neighbor_id) const
virtual void computePhysicalVolumeFraction()=0
Computes the volume fraction of the element fragment.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool isElemAtCrackTip(const Elem *elem) const
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
std::vector< dof_id_type > getElementSolutionDofs(const Elem *elem, SystemBase &sys) const
Get a vector of the dof indices for all components of all variables associated with an element...
virtual System & system() override
Xfem::XFEM_QRULE & getXFEMQRule()
Real getPhysicalVolumeFraction(const Elem *elem) const
Get the volume fraction of an element that is physical.
virtual bool getXFEMWeights(MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points) override
virtual void getXFEMqRuleOnSurface(std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const
EFAElement * add2DElement(const std::vector< unsigned int > &quad, unsigned int id)
bool markCutFacesByGeometry()
unsigned int getLocalNodeIndex(EFANode *node) const
void dataStore(std::ostream &stream, FeatureFloodCount::FeatureData &feature, void *context)
void setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment)
void addStateMarkedElem(unsigned int elem_id, RealVectorValue &normal)
EFAElement * getElemByID(unsigned int id)
MeshBase * _displaced_mesh
const std::set< SubdomainID > & getSubdomainsForVar(unsigned int var_number) const
std::map< unique_id_type, std::vector< Real > > _cached_solution
Data structure to store the nonlinear solution for nodes/elements affected by XFEM For each node/elem...
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
bool isElemCut(const Elem *elem, XFEMCutElem *&xfce) const
virtual void set(const numeric_index_type i, const Number value)=0
virtual bool updateHeal() override
const ConsoleStream _console
virtual void serializeSolution()
void clearStateMarkedElems()
void getWeightMultipliers(MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
virtual bool isTransient() const override
virtual bool cutElementByCrackGrowthIncrement(const Elem *elem, std::vector< CutEdgeForCrackGrowthIncr > &cut_edges, Real time)
void setXFEMQRule(std::string &xfem_qrule)
EFAFace * getFragmentFace(unsigned int frag_id, unsigned int face_id) const
NumericVector< Number > & solutionOld()
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Assembly & assembly(const THREAD_ID tid, const unsigned int nl_sys_num) override
void dataLoad(std::istream &stream, FeatureFloodCount::FeatureData &feature, void *context)
Information about a cut element.
EFANode * getNode(unsigned int index) const
const std::set< EFAElement * > & getCrackTipElements()
static const std::string k
unsigned int numFaceNeighbors(unsigned int face_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=nullptr) const
std::map< const Elem *, RealVectorValue > _state_marked_elems
EFAElement2D * getEFAElem2D(const Elem *elem)
Get the EFAElement2D object for a specified libMesh element.
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const =0
void ErrorVector unsigned int
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
void initCrackTipTopology()
bool markCutFacesByState()
void setInterfaceID(unsigned int interface_id)
Set the interface ID for this cutting object.
void addStateMarkedFrag(unsigned int elem_id, RealVectorValue &normal)
unsigned int numFaces() const
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const =0
std::unordered_map< const Elem *, Xfem::CutElemInfo > _geom_cut_elems
All geometrically cut elements and their CutElemInfo during the current execution of XFEM_MARK...
bool isPointPhysical(const Point &p) const
std::unordered_map< const Elem *, Xfem::CutElemInfo > _old_geom_cut_elems
All geometrically cut elements and their CutElemInfo before the current execution of XFEM_MARK...