20 #include "libmesh/threads.h" 31 closest_point_to_edge(
const Point & src,
const Point & p0,
const Point & p1)
33 const Point line01 = p1 - p0;
34 const Real line0c_xi = ((src - p0) * line01) / line01.
norm_sq();
42 return p0 + line0c_xi * line01;
46 closest_point_to_side(
const Point & src,
const Elem & side)
54 "Penetration of elements with curved sides not implemented");
55 return closest_point_to_edge(src, side.
point(0), side.
point(1));
60 "Penetration of elements with curved sides not implemented");
62 const Point l01 = p1 - p0, l02 = p2 - p0;
63 const Point tri_normal = (l01.
cross(l02)).unit();
64 const Point linecs = ((src - p0) * tri_normal) / tri_normal.
norm_sq() * tri_normal;
65 const Point in_plane = src - linecs;
66 const Point planar_offset = in_plane - p0;
69 if (planar_offset.
cross(l01) * tri_normal > 0)
70 return closest_point_to_edge(src, p0, p1);
73 if (planar_offset.
cross(l02) * tri_normal < 0)
74 return closest_point_to_edge(src, p0, p2);
77 if ((in_plane - p1).cross(p2 - p1) * tri_normal > 0)
78 return closest_point_to_edge(src, p1, p2);
103 std::map<dof_id_type, PenetrationInfo *> & penetration_info,
104 bool check_whether_reasonable,
105 bool update_location,
106 Real tangential_tolerance,
107 bool do_normal_smoothing,
108 Real normal_smoothing_distance,
110 bool use_point_locator,
111 std::vector<std::vector<FEBase *>> & fes,
114 const std::map<
dof_id_type, std::vector<dof_id_type>> & node_to_elem_map)
115 : _subproblem(subproblem),
117 _primary_boundary(primary_boundary),
118 _secondary_boundary(secondary_boundary),
119 _penetration_info(penetration_info),
120 _check_whether_reasonable(check_whether_reasonable),
121 _update_location(update_location),
122 _tangential_tolerance(tangential_tolerance),
123 _do_normal_smoothing(do_normal_smoothing),
124 _normal_smoothing_distance(normal_smoothing_distance),
125 _normal_smoothing_method(normal_smoothing_method),
126 _use_point_locator(use_point_locator),
127 _nodal_normal_x(NULL),
128 _nodal_normal_y(NULL),
129 _nodal_normal_z(NULL),
132 _nearest_node(nearest_node),
133 _node_to_elem_map(node_to_elem_map)
139 : _subproblem(x._subproblem),
141 _primary_boundary(x._primary_boundary),
142 _secondary_boundary(x._secondary_boundary),
143 _penetration_info(x._penetration_info),
144 _check_whether_reasonable(x._check_whether_reasonable),
145 _update_location(x._update_location),
146 _tangential_tolerance(x._tangential_tolerance),
147 _do_normal_smoothing(x._do_normal_smoothing),
148 _normal_smoothing_distance(x._normal_smoothing_distance),
149 _normal_smoothing_method(x._normal_smoothing_method),
150 _use_point_locator(x._use_point_locator),
152 _fe_type(x._fe_type),
153 _nearest_node(x._nearest_node),
154 _node_to_elem_map(x._node_to_elem_map)
174 std::unique_ptr<PointLocatorBase> point_locator;
178 for (
const auto & node_id : range)
189 std::vector<PenetrationInfo *> p_info;
190 bool info_set(
false);
200 const Point contact_ref =
info->_closest_point_ref;
201 bool contact_point_on_side(
false);
205 std::vector<Point> points(1);
206 points[0] = contact_ref;
207 const std::vector<Point> & secondary_pos = fe_side->
get_xyz();
208 bool search_succeeded =
false;
217 contact_point_on_side,
221 info->_closest_point_ref = contact_ref;
224 info->_distance = 0.0;
229 Real old_tangential_distance(
info->_tangential_distance);
230 bool contact_point_on_side(
false);
231 bool search_succeeded =
false;
240 contact_point_on_side,
243 if (contact_point_on_side)
245 if (
info->_tangential_distance <= 0.0)
249 else if (
info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
251 if (
info->_side->dim() == 2 &&
info->_off_edge_nodes.size() < 2)
268 std::vector<dof_id_type> located_elem_ids;
269 const std::vector<dof_id_type> * closest_elems;
273 std::set<const Elem *> candidate_elements;
274 (*point_locator)(*closest_node, candidate_elements);
275 for (
const Elem * elem : candidate_elements)
277 for (
auto s : elem->side_index_range())
280 located_elem_ids.push_back(elem->id());
284 closest_elems = &located_elem_ids;
290 "Missing entry in node to elem map");
291 closest_elems = &(node_to_elem_pair->second);
294 for (
const auto & elem_id : *closest_elems)
298 std::vector<PenetrationInfo *> thisElemInfo;
300 std::vector<const Node *> nodesThatMustBeOnSide;
311 nodesThatMustBeOnSide.push_back(closest_node);
323 for (
unsigned int i = 0; i < p_info.size(); ++i)
325 const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
326 const Real distance_sq = (closest_point - node).
norm_sq();
327 if (distance_sq < min_distance_sq)
329 min_distance_sq = distance_sq;
330 best_point = closest_point;
335 p_info[best_i]->_closest_point = best_point;
336 p_info[best_i]->_distance =
337 (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) *
std::sqrt(min_distance_sq);
339 mooseError(
"Normal smoothing not implemented with point locator code");
340 Point normal = (best_point - node).unit();
341 const Real dot = normal * p_info[best_i]->_normal;
344 p_info[best_i]->_normal = normal;
351 if (p_info.size() == 1)
359 else if (p_info.size() > 1)
362 std::vector<RidgeData> ridgeDataVec;
363 for (
unsigned int i = 0; i + 1 < p_info.size(); ++i)
364 for (
unsigned int j = i + 1; j < p_info.size(); ++j)
367 Real tangential_distance(0.0);
368 const Node * closest_node_on_ridge = NULL;
369 unsigned int index = 0;
370 Point closest_coor_ref;
373 closest_node_on_ridge,
379 if (found_ridge_contact_point)
387 ridgeDataVec.push_back(rpd);
391 if (ridgeDataVec.size() > 0)
395 std::vector<RidgeSetData> ridgeSetDataVec;
396 for (
unsigned int i = 0; i < ridgeDataVec.size(); ++i)
398 bool foundSetWithMatchingNode =
false;
399 for (
unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
401 if (ridgeDataVec[i]._closest_node != NULL &&
402 ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
404 foundSetWithMatchingNode =
true;
405 ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
409 if (!foundSetWithMatchingNode)
415 ridgeSetDataVec.push_back(rsd);
419 for (
unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
421 if (ridgeSetDataVec[i]._closest_node !=
424 if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1)
426 if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
429 ridgeSetDataVec[i]._closest_coor =
430 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
431 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
432 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
438 ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
439 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
440 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
445 ridgeSetDataVec[i]._closest_coor =
446 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
447 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
448 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
452 unsigned int closest_ridge_set_index(0);
453 Real closest_distance(ridgeSetDataVec[0]._distance);
454 Point closest_point(ridgeSetDataVec[0]._closest_coor);
455 for (
unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
457 if (ridgeSetDataVec[i]._distance < closest_distance)
459 closest_ridge_set_index = i;
460 closest_distance = ridgeSetDataVec[i]._distance;
461 closest_point = ridgeSetDataVec[i]._closest_coor;
465 if (closest_distance <
471 for (
unsigned int i = 0;
472 i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
475 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
476 face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
480 "face_index invalid");
482 p_info[face_index]->_closest_point = closest_point;
483 p_info[face_index]->_distance =
484 (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
490 Point normal(closest_point - node);
496 const Real dot(normal * p_info[face_index]->_normal);
501 p_info[face_index]->_normal = normal;
503 p_info[face_index]->_tangential_distance = 0.0;
505 Point closest_point_ref;
506 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
509 p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
511 ._tangential_distance;
512 p_info[face_index]->_closest_point_ref =
513 ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
517 const Node * closest_node_on_face;
519 closest_node_on_face,
520 p_info[face_index]->_side);
523 if (closest_node_on_face !=
524 ridgeSetDataVec[closest_ridge_set_index]._closest_node)
526 mooseError(
"Closest node when restricting point to face != closest node from " 533 std::vector<Point> points(1);
534 points[0] = p_info[face_index]->_closest_point_ref;
535 fe->
reinit(p_info[face_index]->_side, &points);
536 p_info[face_index]->_side_phi = fe->
get_phi();
537 p_info[face_index]->_side_grad_phi = fe->
get_dphi();
553 unsigned int best(0), i(1);
571 }
while (i < p_info.size() && best < p_info.size());
572 if (best < p_info.size())
605 for (
unsigned int j = 0; j < p_info.size(); ++j)
627 mooseAssert(infoNew != NULL,
"infoNew object is null");
754 const std::vector<const Node *> & off_edge_nodes1 = pi1->
_off_edge_nodes;
755 const std::vector<const Node *> & off_edge_nodes2 = pi2->
_off_edge_nodes;
756 const unsigned dim1 = pi1->
_side->
dim();
760 mooseAssert(pi2->
_side->
dim() == 1,
"Incompatible dimensions.");
761 mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
762 "off_edge_nodes size should be <2 for 2D contact");
763 if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
764 off_edge_nodes1[0] == off_edge_nodes2[0])
769 mooseAssert(dim1 == 2 && pi2->
_side->
dim() == 2,
"Incompatible dimensions.");
770 mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
771 "off_edge_nodes size should be <3 for 3D contact");
772 if (off_edge_nodes1.size() == 1)
774 if (off_edge_nodes2.size() == 1)
776 if (off_edge_nodes1[0] == off_edge_nodes2[0])
779 else if (off_edge_nodes2.size() == 2)
781 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
785 else if (off_edge_nodes1.size() == 2)
787 if (off_edge_nodes2.size() == 1)
789 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
792 else if (off_edge_nodes2.size() == 2)
794 if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
795 off_edge_nodes1[1] == off_edge_nodes2[1]) ||
796 (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
806 Real & tangential_distance,
807 const Node *& closest_node,
808 unsigned int & index,
809 Point & contact_point_ref,
810 std::vector<PenetrationInfo *> & p_info,
811 const unsigned int index1,
812 const unsigned int index2)
814 tangential_distance = 0.0;
818 const unsigned sidedim(pi1->
_side->
dim());
819 mooseAssert(sidedim == pi2->
_side->
dim(),
"Incompatible dimensionalities");
822 std::vector<const Node *> side1_nodes;
824 std::vector<const Node *> side2_nodes;
827 std::sort(side1_nodes.begin(), side1_nodes.end());
828 std::sort(side2_nodes.begin(), side2_nodes.end());
831 std::vector<const Node *> common_nodes;
832 std::set_intersection(side1_nodes.begin(),
836 std::inserter(common_nodes, common_nodes.end()));
838 if (common_nodes.size() != sidedim)
841 bool found_point1, found_point2;
843 const Node * closest_node1;
845 closest_coor_ref1, closest_node1, pi1->
_side, common_nodes);
848 const Node * closest_node2;
850 closest_coor_ref2, closest_node2, pi2->
_side, common_nodes);
852 if (!found_point1 || !found_point2)
865 std::vector<Point> points(1);
872 contact_point_ref = closest_coor_ref1;
873 points[0] = closest_coor_ref1;
880 contact_point_ref = closest_coor_ref2;
881 points[0] = closest_coor_ref2;
886 contact_point = fe->
get_xyz()[0];
892 mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
893 "If off edge of ridge, closest node must be the same on both elements");
894 closest_node = closest_node1;
897 tangential_distance = off_face.
norm();
908 corner_nodes.clear();
910 corner_nodes.push_back(side->
node_ptr(0));
911 corner_nodes.push_back(side->
node_ptr(1));
925 corner_nodes.push_back(side->
node_ptr(2));
933 corner_nodes.push_back(side->
node_ptr(2));
934 corner_nodes.push_back(side->
node_ptr(3));
948 const Node *& closest_node,
950 const std::vector<const Node *> & edge_nodes)
957 std::vector<unsigned int> local_node_indices;
958 for (
const auto & edge_node : edge_nodes)
963 local_node_indices.push_back(local_index);
965 mooseAssert(local_node_indices.size() == side->
dim(),
966 "Number of edge nodes must match side dimensionality");
967 std::sort(local_node_indices.begin(), local_node_indices.end());
969 bool off_of_this_edge =
false;
977 if (local_node_indices[0] == 0)
982 off_of_this_edge =
true;
986 else if (local_node_indices[0] == 1)
991 off_of_this_edge =
true;
1006 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1011 off_of_this_edge =
true;
1018 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1020 if ((xi + eta) > 1.0)
1022 Real delta = (xi + eta - 1.0) / 2.0;
1025 off_of_this_edge =
true;
1032 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1037 off_of_this_edge =
true;
1056 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1061 off_of_this_edge =
true;
1068 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1073 off_of_this_edge =
true;
1080 else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1085 off_of_this_edge =
true;
1092 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1097 off_of_this_edge =
true;
1117 return off_of_this_edge;
1126 closest_node = NULL;
1128 bool off_of_this_face(
false);
1139 off_of_this_face =
true;
1145 off_of_this_face =
true;
1158 off_of_this_face =
true;
1172 else if ((xi + eta) > 1.0)
1174 Real delta = (xi + eta - 1.0) / 2.0;
1177 off_of_this_face =
true;
1200 off_of_this_face =
true;
1224 off_of_this_face =
true;
1241 off_of_this_face =
true;
1258 off_of_this_face =
true;
1275 off_of_this_face =
true;
1298 return off_of_this_face;
1305 const Point * secondary_point,
1306 const Real tangential_tolerance)
1308 unsigned int dim = primary_elem->
dim();
1310 const std::vector<Point> & phys_point = fe->
get_xyz();
1312 const std::vector<RealGradient> & dxyz_dxi = fe->
get_dxyzdxi();
1313 const std::vector<RealGradient> & dxyz_deta = fe->
get_dxyzdeta();
1317 std::vector<Point> points(1);
1319 fe->
reinit(side, &points);
1323 const Real twosqrt2 = 2.8284;
1324 Real max_face_length = side->
hmax() + twosqrt2 * tangential_tolerance;
1329 normal = dxyz_dxi[0].
cross(dxyz_deta[0]);
1331 else if (
dim - 1 == 1)
1333 const Node *
const * elem_nodes = primary_elem->
get_nodes();
1334 const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1335 const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1337 Point out_of_plane_normal = in_plane_vector1.
cross(in_plane_vector2);
1338 out_of_plane_normal /= out_of_plane_normal.
norm();
1340 normal = dxyz_dxi[0].
cross(out_of_plane_normal);
1346 normal /= normal.
norm();
1348 const Real dot(d * normal);
1353 const Real tangdist = tangcomp.
norm();
1357 const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.
norm());
1359 bool isReasonableCandidate =
true;
1360 if (tangdist > faceExpansionFactor * max_face_length)
1362 isReasonableCandidate =
false;
1364 return isReasonableCandidate;
1372 std::vector<Point> points(1);
1373 points[0] =
info._starting_closest_point_ref;
1375 fe.
reinit(&side, &points);
1376 const std::vector<Point> & starting_point = fe.
get_xyz();
1377 info._incremental_slip =
info._closest_point - starting_point[0];
1378 if (
info.isCaptured())
1380 info._frictional_energy =
1381 info._frictional_energy_old +
info._contact_force *
info._incremental_slip;
1382 info._accumulated_slip =
info._accumulated_slip_old +
info._incremental_slip.norm();
1388 std::vector<PenetrationInfo *> & p_info,
1397 std::vector<Real> edge_face_weights;
1398 std::vector<PenetrationInfo *> edge_face_info;
1402 mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1403 "edge_face_info.size() != edge_face_weights.size()");
1405 if (edge_face_info.size() > 0)
1409 Real this_face_weight = 1.0;
1411 for (
unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1415 new_normal += npi->
_normal * edge_face_weights[efwi];
1417 this_face_weight -= edge_face_weights[efwi];
1419 mooseAssert(this_face_weight >= (0.25 - 1e-8),
1420 "Sum of weights of other faces shouldn't exceed 0.75");
1421 new_normal +=
info->_normal * this_face_weight;
1423 const Real len = new_normal.
norm();
1427 info->_normal = new_normal;
1437 const Real len(
info->_normal.norm());
1439 info->_normal /= len;
1446 std::vector<PenetrationInfo *> & edge_face_info,
1447 std::vector<Real> & edge_face_weights,
1448 std::vector<PenetrationInfo *> & p_info,
1449 const Node & secondary_node)
1452 const Point & p =
info->_closest_point_ref;
1453 std::set<dof_id_type> elems_to_exclude;
1454 elems_to_exclude.insert(
info->_elem->id());
1456 std::vector<std::vector<const Node *>> edge_nodes;
1460 std::vector<Elem *> edge_neighbor_elems;
1461 edge_face_info.resize(edge_nodes.size(), NULL);
1463 std::vector<unsigned int> edges_without_neighbors;
1465 for (
unsigned int i = 0; i < edge_nodes.size(); ++i)
1468 std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1470 std::vector<PenetrationInfo *> face_info_comm_edge;
1472 &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1474 if (face_info_comm_edge.size() == 0)
1475 edges_without_neighbors.push_back(i);
1476 else if (face_info_comm_edge.size() > 1)
1477 mooseError(
"Only one neighbor allowed per edge");
1479 edge_face_info[i] = face_info_comm_edge[0];
1483 std::vector<unsigned int>::reverse_iterator rit;
1484 for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1486 unsigned int index = *rit;
1487 edge_nodes.erase(edge_nodes.begin() + index);
1488 edge_face_weights.erase(edge_face_weights.begin() + index);
1489 edge_face_info.erase(edge_face_info.begin() + index);
1493 if (edge_nodes.size() > 1)
1495 if (edge_nodes.size() != 2)
1496 mooseError(
"Invalid number of smoothing edges");
1499 std::vector<const Node *> common_nodes;
1500 std::set_intersection(edge_nodes[0].begin(),
1501 edge_nodes[0].end(),
1502 edge_nodes[1].begin(),
1503 edge_nodes[1].end(),
1504 std::inserter(common_nodes, common_nodes.end()));
1506 if (common_nodes.size() != 1)
1507 mooseError(
"Invalid number of common nodes");
1509 for (
const auto & pinfo : edge_face_info)
1510 elems_to_exclude.insert(pinfo->_elem->id());
1512 std::vector<PenetrationInfo *> face_info_comm_edge;
1514 &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1516 unsigned int num_corner_neighbors = face_info_comm_edge.size();
1518 if (num_corner_neighbors > 0)
1520 Real fw0 = edge_face_weights[0];
1521 Real fw1 = edge_face_weights[1];
1524 Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1527 edge_face_weights[0] *= (1.0 - fw1);
1528 edge_face_weights[1] *= (1.0 - fw0);
1530 for (
unsigned int i = 0; i < num_corner_neighbors; ++i)
1532 edge_face_weights.push_back(fw_corner);
1533 edge_face_info.push_back(face_info_comm_edge[i]);
1543 std::vector<std::vector<const Node *>> & edge_nodes,
1544 std::vector<Real> & edge_face_weights)
1547 const Real & xi = p(0);
1548 const Real & eta = p(1);
1558 if (xi < -smooth_limit)
1560 std::vector<const Node *> en;
1562 edge_nodes.push_back(en);
1566 edge_face_weights.push_back(fw);
1568 else if (xi > smooth_limit)
1570 std::vector<const Node *> en;
1572 edge_nodes.push_back(en);
1576 edge_face_weights.push_back(fw);
1585 if (eta < -smooth_limit)
1587 std::vector<const Node *> en;
1590 edge_nodes.push_back(en);
1594 edge_face_weights.push_back(fw);
1596 if ((xi + eta) > smooth_limit)
1598 std::vector<const Node *> en;
1601 edge_nodes.push_back(en);
1605 edge_face_weights.push_back(fw);
1607 if (xi < -smooth_limit)
1609 std::vector<const Node *> en;
1612 edge_nodes.push_back(en);
1616 edge_face_weights.push_back(fw);
1625 if (eta < -smooth_limit)
1627 std::vector<const Node *> en;
1630 edge_nodes.push_back(en);
1634 edge_face_weights.push_back(fw);
1636 if (xi > smooth_limit)
1638 std::vector<const Node *> en;
1641 edge_nodes.push_back(en);
1645 edge_face_weights.push_back(fw);
1647 if (eta > smooth_limit)
1649 std::vector<const Node *> en;
1652 edge_nodes.push_back(en);
1656 edge_face_weights.push_back(fw);
1658 if (xi < -smooth_limit)
1660 std::vector<const Node *> en;
1663 edge_nodes.push_back(en);
1667 edge_face_weights.push_back(fw);
1682 const Node * secondary_node,
1683 const std::set<dof_id_type> & elems_to_exclude,
1684 const std::vector<const Node *> edge_nodes,
1685 std::vector<PenetrationInfo *> & face_info_comm_edge,
1686 std::vector<PenetrationInfo *> & p_info)
1692 mooseAssert(node_to_elem_pair !=
_node_to_elem_map.end(),
"Missing entry in node to elem map");
1693 const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1695 std::vector<const Elem *> elems_connected_to_edge;
1697 for (
unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1699 if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1703 std::vector<const Node *> nodevec;
1704 for (
unsigned int ni = 0; ni < elem->
n_nodes(); ++ni)
1706 nodevec.push_back(elem->
node_ptr(ni));
1708 std::vector<const Node *> common_nodes;
1709 std::sort(nodevec.begin(), nodevec.end());
1710 std::set_intersection(edge_nodes.begin(),
1714 std::inserter(common_nodes, common_nodes.end()));
1716 if (common_nodes.size() == edge_nodes.size())
1717 elems_connected_to_edge.push_back(elem);
1720 if (elems_connected_to_edge.size() > 0)
1730 bool allowMultipleNeighbors =
false;
1732 if (elems_connected_to_edge[0]->
dim() == 3)
1734 if (edge_nodes.size() == 1)
1736 allowMultipleNeighbors =
true;
1740 for (
unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1742 std::vector<PenetrationInfo *> thisElemInfo;
1743 getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1744 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1746 if (thisElemInfo.size() > 1)
1748 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1749 face_info_comm_edge.push_back(thisElemInfo[0]);
1754 thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1755 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1757 if (thisElemInfo.size() > 1)
1759 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1760 face_info_comm_edge.push_back(thisElemInfo[0]);
1764 for (
unsigned int j = 0; j < thisElemInfo.size(); ++j)
1765 face_info_comm_edge.push_back(thisElemInfo[j]);
1772 std::vector<PenetrationInfo *> & p_info,
1775 for (
const auto &
pi : p_info)
1780 if (
pi->_elem == elem)
1781 thisElemInfo.push_back(
pi);
1787 std::vector<PenetrationInfo *> & p_info,
1788 const Node * secondary_node,
1790 const std::vector<const Node *> & nodes_that_must_be_on_side,
1791 const bool check_whether_reasonable)
1801 bool already_have_info_this_side =
false;
1802 for (
const auto &
pi : thisElemInfo)
1803 if (
pi->_side_num == s)
1805 already_have_info_this_side =
true;
1809 if (already_have_info_this_side)
1816 std::vector<const Node *> nodevec;
1817 for (
unsigned int ni = 0; ni < side->
n_nodes(); ++ni)
1818 nodevec.push_back(side->
node_ptr(ni));
1820 std::sort(nodevec.begin(), nodevec.end());
1821 std::vector<const Node *> common_nodes;
1822 std::set_intersection(nodes_that_must_be_on_side.begin(),
1823 nodes_that_must_be_on_side.end(),
1826 std::inserter(common_nodes, common_nodes.end()));
1827 if (common_nodes.size() != nodes_that_must_be_on_side.size())
1838 if (check_whether_reasonable)
1847 Point contact_on_face_ref;
1849 Real tangential_distance = 0.;
1851 bool contact_point_on_side;
1852 std::vector<const Node *> off_edge_nodes;
1853 std::vector<std::vector<Real>> side_phi;
1854 std::vector<std::vector<RealGradient>> side_grad_phi;
1855 std::vector<RealGradient> dxyzdxi;
1856 std::vector<RealGradient> dxyzdeta;
1857 std::vector<RealGradient> d2xyzdxideta;
1859 std::unique_ptr<PenetrationInfo> pen_info =
1860 std::make_unique<PenetrationInfo>(elem,
1865 tangential_distance,
1868 contact_on_face_ref,
1876 bool search_succeeded =
false;
1884 contact_point_on_side,
1888 if (search_succeeded)
1890 thisElemInfo.push_back(pen_info.get());
1891 p_info.push_back(pen_info.release());
MooseVariable * _nodal_normal_z
PenetrationThread(SubProblem &subproblem, const MooseMesh &mesh, BoundaryID primary_boundary, BoundaryID secondary_boundary, std::map< dof_id_type, PenetrationInfo *> &penetration_info, bool check_whether_reasonable, bool update_location, Real tangential_tolerance, bool do_normal_smoothing, Real normal_smoothing_distance, PenetrationLocator::NORMAL_SMOOTHING_METHOD normal_smoothing_method, bool use_point_locator, std::vector< std::vector< libMesh::FEBase *>> &fes, libMesh::FEType &fe_type, NearestNodeLocator &nearest_node, const std::map< dof_id_type, std::vector< dof_id_type >> &node_to_elem_map)
void getSmoothingEdgeNodesAndWeights(const libMesh::Point &p, const Elem *side, std::vector< std::vector< const Node *>> &edge_nodes, std::vector< Real > &edge_face_weights)
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
auto norm() const -> decltype(std::norm(Real()))
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
const unsigned int invalid_uint
unsigned int get_node_index(const Node *node_ptr) const
Real _accumulated_slip_old
virtual_for_inffe const std::vector< RealGradient > & get_dxyzdeta() const
virtual Elem * elemPtr(const dof_id_type i)
bool findRidgeContactPoint(libMesh::Point &contact_point, Real &tangential_distance, const Node *&closest_node, unsigned int &index, libMesh::Point &contact_point_ref, std::vector< PenetrationInfo *> &p_info, const unsigned int index1, const unsigned int index2)
OutputType getValue(const Elem *elem, const std::vector< std::vector< OutputShape >> &phi) const
Compute the variable value at a point on an element.
IntRange< unsigned short > side_index_range() const
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Data structure used to hold penetration information.
NearestNodeLocator & _nearest_node
Finds the nearest node to each node in boundary1 to each node in boundary2 and the other way around...
virtual bool has_affine_map() const
Real _normal_smoothing_distance
const Elem * _starting_elem
bool _check_whether_reasonable
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
virtual_for_inffe const std::vector< RealGradient > & get_d2xyzdxideta() const
void createInfoForElem(std::vector< PenetrationInfo *> &thisElemInfo, std::vector< PenetrationInfo *> &p_info, const Node *secondary_node, const Elem *elem, const std::vector< const Node *> &nodes_that_must_be_on_side, const bool check_whether_reasonable=false)
unsigned int _stick_locked_this_step
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
void getSideCornerNodes(const Elem *side, std::vector< const Node *> &corner_nodes)
Threads::spin_mutex pinfo_mutex
virtual Real hmax() const
const BoundaryInfo & get_boundary_info() const
RealVectorValue _contact_force_old
Real distance(const Point &p)
virtual const Node & nodeRef(const dof_id_type i) const
void smoothNormal(PenetrationInfo *info, std::vector< PenetrationInfo *> &p_info, const Node &node)
libMesh::FEType & _fe_type
void getInfoForFacesWithCommonNodes(const Node *secondary_node, const std::set< dof_id_type > &elems_to_exclude, const std::vector< const Node *> edge_nodes, std::vector< PenetrationInfo *> &face_info_comm_edge, std::vector< PenetrationInfo *> &p_info)
auto max(const L &left, const R &right)
void getSmoothingFacesAndWeights(PenetrationInfo *info, std::vector< PenetrationInfo *> &edge_face_info, std::vector< Real > &edge_face_weights, std::vector< PenetrationInfo *> &p_info, const Node &secondary_node)
unsigned int _locked_this_step
BoundaryID _primary_boundary
Real _tangential_distance
const std::vector< std::vector< OutputGradient > > & get_dphi() const
auto norm_sq() const -> decltype(std::norm(Real()))
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
virtual unsigned int n_nodes() const=0
unsigned int _starting_side_num
boundary_id_type BoundaryID
CompeteInteractionResult competeInteractions(PenetrationInfo *pi1, PenetrationInfo *pi2)
When interactions are identified between a node and two faces, compete between the faces to determine...
void operator()(const NodeIdRange &range)
const Node * _closest_node
const Node *const * get_nodes() const
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
void computeSlip(libMesh::FEBase &fe, PenetrationInfo &info)
const Node * _closest_node
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
void findContactPoint(PenetrationInfo &p_info, libMesh::FEBase *fe_elem, libMesh::FEBase *fe_side, libMesh::FEType &fe_side_type, const libMesh::Point &secondary_point, bool start_with_centroid, const Real tangential_tolerance, bool &contact_point_on_side, bool &search_succeeded)
Finds the closest point (called the contact point) on the primary_elem on side "side" to the secondar...
std::vector< dof_id_type > _recheck_secondary_nodes
List of secondary nodes for which penetration was not detected in the current patch and for which pat...
virtual MooseVariable & getStandardVariable(const THREAD_ID tid, const std::string &var_name)=0
Returns the variable reference for requested MooseVariable which may be in any system.
libMesh::Point _closest_coor
virtual_for_inffe const std::vector< Point > & get_xyz() const
bool restrictPointToSpecifiedEdgeOfFace(libMesh::Point &p, const Node *&closest_node, const Elem *side, const std::vector< const Node *> &edge_nodes)
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
RealVectorValue _contact_force
Real _frictional_energy_old
MECH_STATUS_ENUM _mech_status
void join(const PenetrationThread &other)
Real _tangential_tolerance
std::vector< const Node * > _off_edge_nodes
Point _starting_closest_point_ref
Real _tangential_distance
bool restrictPointToFace(libMesh::Point &p, const Node *&closest_node, const Elem *side)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Generic class for solving transient nonlinear problems.
Real _lagrange_multiplier
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
const Node * nearestNode(dof_id_type node_id)
Valid to call this after findNodes() has been called to get a pointer to the nearest node...
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
virtual unsigned short dim() const=0
MooseVariable * _nodal_normal_y
const Node * node_ptr(const unsigned int i) const
virtual bool is_vertex(const unsigned int i) const=0
libMesh::ElemSideBuilder _elem_side_builder
Helper for building element sides without extraneous allocation.
void getInfoForElem(std::vector< PenetrationInfo *> &thisElemInfo, std::vector< PenetrationInfo *> &p_info, const Elem *elem)
bool _do_normal_smoothing
MECH_STATUS_ENUM _mech_status_old
libMesh::Point _closest_coor_ref
virtual std::unique_ptr< libMesh::PointLocatorBase > getPointLocator() const
Proxy function to get a (sub)PointLocator from either the underlying libMesh mesh (default)...
const std::map< dof_id_type, std::vector< dof_id_type > > & _node_to_elem_map
void switchInfo(PenetrationInfo *&info, PenetrationInfo *&infoNew)
virtual_for_inffe const std::vector< RealGradient > & get_dxyzdxi() const
RealVectorValue _lagrange_multiplier_slip
MooseVariable * _nodal_normal_x
bool relativeFuzzyLessThan(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether a variable is less than another variable within a relative tolerance...
std::vector< std::vector< libMesh::FEBase * > > & _fes
virtual ElemType type() const=0
CompeteInteractionResult competeInteractionsBothOnFace(PenetrationInfo *pi1, PenetrationInfo *pi2)
Determine whether first (pi1) or second (pi2) interaction is stronger when it is known that the node ...
const Point & point(const unsigned int i) const
CommonEdgeResult interactionsOffCommonEdge(PenetrationInfo *pi1, PenetrationInfo *pi2)
bool isFaceReasonableCandidate(const Elem *primary_elem, const Elem *side, libMesh::FEBase *fe, const libMesh::Point *secondary_point, const Real tangential_tolerance)
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method
const std::vector< std::vector< OutputShape > > & get_phi() const
std::vector< RidgeData > _ridge_data_vec