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);
276 if (candidate_elements.empty())
277 mooseError(
"No proximate elements found at node ",
280 static_cast<const Point &
>(*closest_node),
283 ". This should never happen.");
285 for (
const Elem * elem : candidate_elements)
287 for (
auto s : elem->side_index_range())
290 located_elem_ids.push_back(elem->id());
295 if (located_elem_ids.empty())
296 mooseError(
"No proximate elements found at node ",
299 static_cast<const Point &
>(*closest_node),
302 " share that boundary. This may happen if the mesh uses the same boundary id " 303 "for a nodeset and an unrelated sideset.");
305 closest_elems = &located_elem_ids;
311 "Missing entry in node to elem map");
312 closest_elems = &(node_to_elem_pair->second);
315 for (
const auto & elem_id : *closest_elems)
319 std::vector<PenetrationInfo *> thisElemInfo;
321 std::vector<const Node *> nodesThatMustBeOnSide;
332 nodesThatMustBeOnSide.push_back(closest_node);
344 for (
unsigned int i = 0; i < p_info.size(); ++i)
346 const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
347 const Real distance_sq = (closest_point - node).
norm_sq();
348 if (distance_sq < min_distance_sq)
350 min_distance_sq = distance_sq;
351 best_point = closest_point;
356 p_info[best_i]->_closest_point = best_point;
357 p_info[best_i]->_distance =
358 (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) *
std::sqrt(min_distance_sq);
360 mooseError(
"Normal smoothing not implemented with point locator code");
361 Point normal = (best_point - node).unit();
362 const Real dot = normal * p_info[best_i]->_normal;
365 p_info[best_i]->_normal = normal;
372 if (p_info.size() == 1)
380 else if (p_info.size() > 1)
383 std::vector<RidgeData> ridgeDataVec;
384 for (
unsigned int i = 0; i + 1 < p_info.size(); ++i)
385 for (
unsigned int j = i + 1; j < p_info.size(); ++j)
388 Real tangential_distance(0.0);
389 const Node * closest_node_on_ridge = NULL;
390 unsigned int index = 0;
391 Point closest_coor_ref;
394 closest_node_on_ridge,
400 if (found_ridge_contact_point)
408 ridgeDataVec.push_back(rpd);
412 if (ridgeDataVec.size() > 0)
416 std::vector<RidgeSetData> ridgeSetDataVec;
417 for (
unsigned int i = 0; i < ridgeDataVec.size(); ++i)
419 bool foundSetWithMatchingNode =
false;
420 for (
unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
422 if (ridgeDataVec[i]._closest_node != NULL &&
423 ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
425 foundSetWithMatchingNode =
true;
426 ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
430 if (!foundSetWithMatchingNode)
436 ridgeSetDataVec.push_back(rsd);
440 for (
unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
442 if (ridgeSetDataVec[i]._closest_node !=
445 if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1)
447 if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
450 ridgeSetDataVec[i]._closest_coor =
451 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
452 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
453 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
459 ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
460 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
461 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
466 ridgeSetDataVec[i]._closest_coor =
467 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
468 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
469 ridgeSetDataVec[i]._distance = contact_point_vec.
norm();
473 unsigned int closest_ridge_set_index(0);
474 Real closest_distance(ridgeSetDataVec[0]._distance);
475 Point closest_point(ridgeSetDataVec[0]._closest_coor);
476 for (
unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
478 if (ridgeSetDataVec[i]._distance < closest_distance)
480 closest_ridge_set_index = i;
481 closest_distance = ridgeSetDataVec[i]._distance;
482 closest_point = ridgeSetDataVec[i]._closest_coor;
486 if (closest_distance <
492 for (
unsigned int i = 0;
493 i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
496 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
497 face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
501 "face_index invalid");
503 p_info[face_index]->_closest_point = closest_point;
504 p_info[face_index]->_distance =
505 (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
511 Point normal(closest_point - node);
517 const Real dot(normal * p_info[face_index]->_normal);
522 p_info[face_index]->_normal = normal;
524 p_info[face_index]->_tangential_distance = 0.0;
526 Point closest_point_ref;
527 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
530 p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
532 ._tangential_distance;
533 p_info[face_index]->_closest_point_ref =
534 ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
538 const Node * closest_node_on_face;
540 closest_node_on_face,
541 p_info[face_index]->_side);
544 if (closest_node_on_face !=
545 ridgeSetDataVec[closest_ridge_set_index]._closest_node)
547 mooseError(
"Closest node when restricting point to face != closest node from " 554 std::vector<Point> points(1);
555 points[0] = p_info[face_index]->_closest_point_ref;
556 fe->
reinit(p_info[face_index]->_side, &points);
557 p_info[face_index]->_side_phi = fe->
get_phi();
558 p_info[face_index]->_side_grad_phi = fe->
get_dphi();
574 unsigned int best(0), i(1);
592 }
while (i < p_info.size() && best < p_info.size());
593 if (best < p_info.size())
626 for (
unsigned int j = 0; j < p_info.size(); ++j)
648 mooseAssert(infoNew != NULL,
"infoNew object is null");
775 const std::vector<const Node *> & off_edge_nodes1 = pi1->
_off_edge_nodes;
776 const std::vector<const Node *> & off_edge_nodes2 = pi2->
_off_edge_nodes;
777 const unsigned dim1 = pi1->
_side->
dim();
781 mooseAssert(pi2->
_side->
dim() == 1,
"Incompatible dimensions.");
782 mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
783 "off_edge_nodes size should be <2 for 2D contact");
784 if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
785 off_edge_nodes1[0] == off_edge_nodes2[0])
790 mooseAssert(dim1 == 2 && pi2->
_side->
dim() == 2,
"Incompatible dimensions.");
791 mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
792 "off_edge_nodes size should be <3 for 3D contact");
793 if (off_edge_nodes1.size() == 1)
795 if (off_edge_nodes2.size() == 1)
797 if (off_edge_nodes1[0] == off_edge_nodes2[0])
800 else if (off_edge_nodes2.size() == 2)
802 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
806 else if (off_edge_nodes1.size() == 2)
808 if (off_edge_nodes2.size() == 1)
810 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
813 else if (off_edge_nodes2.size() == 2)
815 if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
816 off_edge_nodes1[1] == off_edge_nodes2[1]) ||
817 (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
827 Real & tangential_distance,
828 const Node *& closest_node,
829 unsigned int & index,
830 Point & contact_point_ref,
831 std::vector<PenetrationInfo *> & p_info,
832 const unsigned int index1,
833 const unsigned int index2)
835 tangential_distance = 0.0;
839 const unsigned sidedim(pi1->
_side->
dim());
840 mooseAssert(sidedim == pi2->
_side->
dim(),
"Incompatible dimensionalities");
843 std::vector<const Node *> side1_nodes;
845 std::vector<const Node *> side2_nodes;
848 std::sort(side1_nodes.begin(), side1_nodes.end());
849 std::sort(side2_nodes.begin(), side2_nodes.end());
852 std::vector<const Node *> common_nodes;
853 std::set_intersection(side1_nodes.begin(),
857 std::inserter(common_nodes, common_nodes.end()));
859 if (common_nodes.size() != sidedim)
862 bool found_point1, found_point2;
864 const Node * closest_node1;
866 closest_coor_ref1, closest_node1, pi1->
_side, common_nodes);
869 const Node * closest_node2;
871 closest_coor_ref2, closest_node2, pi2->
_side, common_nodes);
873 if (!found_point1 || !found_point2)
886 std::vector<Point> points(1);
893 contact_point_ref = closest_coor_ref1;
894 points[0] = closest_coor_ref1;
901 contact_point_ref = closest_coor_ref2;
902 points[0] = closest_coor_ref2;
907 contact_point = fe->
get_xyz()[0];
913 mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
914 "If off edge of ridge, closest node must be the same on both elements");
915 closest_node = closest_node1;
918 tangential_distance = off_face.
norm();
929 corner_nodes.clear();
931 corner_nodes.push_back(side->
node_ptr(0));
932 corner_nodes.push_back(side->
node_ptr(1));
946 corner_nodes.push_back(side->
node_ptr(2));
954 corner_nodes.push_back(side->
node_ptr(2));
955 corner_nodes.push_back(side->
node_ptr(3));
969 const Node *& closest_node,
971 const std::vector<const Node *> & edge_nodes)
978 std::vector<unsigned int> local_node_indices;
979 for (
const auto & edge_node : edge_nodes)
984 local_node_indices.push_back(local_index);
986 mooseAssert(local_node_indices.size() == side->
dim(),
987 "Number of edge nodes must match side dimensionality");
988 std::sort(local_node_indices.begin(), local_node_indices.end());
990 bool off_of_this_edge =
false;
998 if (local_node_indices[0] == 0)
1003 off_of_this_edge =
true;
1007 else if (local_node_indices[0] == 1)
1012 off_of_this_edge =
true;
1027 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1032 off_of_this_edge =
true;
1039 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1041 if ((xi + eta) > 1.0)
1043 Real delta = (xi + eta - 1.0) / 2.0;
1046 off_of_this_edge =
true;
1053 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1058 off_of_this_edge =
true;
1077 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1082 off_of_this_edge =
true;
1089 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1094 off_of_this_edge =
true;
1101 else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1106 off_of_this_edge =
true;
1113 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1118 off_of_this_edge =
true;
1138 return off_of_this_edge;
1147 closest_node = NULL;
1149 bool off_of_this_face(
false);
1160 off_of_this_face =
true;
1166 off_of_this_face =
true;
1179 off_of_this_face =
true;
1193 else if ((xi + eta) > 1.0)
1195 Real delta = (xi + eta - 1.0) / 2.0;
1198 off_of_this_face =
true;
1221 off_of_this_face =
true;
1245 off_of_this_face =
true;
1262 off_of_this_face =
true;
1279 off_of_this_face =
true;
1296 off_of_this_face =
true;
1319 return off_of_this_face;
1326 const Point * secondary_point,
1327 const Real tangential_tolerance)
1329 unsigned int dim = primary_elem->
dim();
1331 const std::vector<Point> & phys_point = fe->
get_xyz();
1333 const std::vector<RealGradient> & dxyz_dxi = fe->
get_dxyzdxi();
1334 const std::vector<RealGradient> & dxyz_deta = fe->
get_dxyzdeta();
1338 std::vector<Point> points(1);
1340 fe->
reinit(side, &points);
1344 const Real twosqrt2 = 2.8284;
1345 Real max_face_length = side->
hmax() + twosqrt2 * tangential_tolerance;
1350 normal = dxyz_dxi[0].
cross(dxyz_deta[0]);
1352 else if (
dim - 1 == 1)
1354 const Node *
const * elem_nodes = primary_elem->
get_nodes();
1355 const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1356 const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1358 Point out_of_plane_normal = in_plane_vector1.
cross(in_plane_vector2);
1359 out_of_plane_normal /= out_of_plane_normal.
norm();
1361 normal = dxyz_dxi[0].
cross(out_of_plane_normal);
1367 normal /= normal.
norm();
1369 const Real dot(d * normal);
1374 const Real tangdist = tangcomp.
norm();
1378 const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.
norm());
1380 bool isReasonableCandidate =
true;
1381 if (tangdist > faceExpansionFactor * max_face_length)
1383 isReasonableCandidate =
false;
1385 return isReasonableCandidate;
1393 std::vector<Point> points(1);
1394 points[0] =
info._starting_closest_point_ref;
1396 fe.
reinit(&side, &points);
1397 const std::vector<Point> & starting_point = fe.
get_xyz();
1398 info._incremental_slip =
info._closest_point - starting_point[0];
1399 if (
info.isCaptured())
1401 info._frictional_energy =
1402 info._frictional_energy_old +
info._contact_force *
info._incremental_slip;
1403 info._accumulated_slip =
info._accumulated_slip_old +
info._incremental_slip.norm();
1409 std::vector<PenetrationInfo *> & p_info,
1418 std::vector<Real> edge_face_weights;
1419 std::vector<PenetrationInfo *> edge_face_info;
1423 mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1424 "edge_face_info.size() != edge_face_weights.size()");
1426 if (edge_face_info.size() > 0)
1430 Real this_face_weight = 1.0;
1432 for (
unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1436 new_normal += npi->
_normal * edge_face_weights[efwi];
1438 this_face_weight -= edge_face_weights[efwi];
1440 mooseAssert(this_face_weight >= (0.25 - 1e-8),
1441 "Sum of weights of other faces shouldn't exceed 0.75");
1442 new_normal +=
info->_normal * this_face_weight;
1444 const Real len = new_normal.
norm();
1448 info->_normal = new_normal;
1458 const Real len(
info->_normal.norm());
1460 info->_normal /= len;
1467 std::vector<PenetrationInfo *> & edge_face_info,
1468 std::vector<Real> & edge_face_weights,
1469 std::vector<PenetrationInfo *> & p_info,
1470 const Node & secondary_node)
1473 const Point & p =
info->_closest_point_ref;
1474 std::set<dof_id_type> elems_to_exclude;
1475 elems_to_exclude.insert(
info->_elem->id());
1477 std::vector<std::vector<const Node *>> edge_nodes;
1481 std::vector<Elem *> edge_neighbor_elems;
1482 edge_face_info.resize(edge_nodes.size(), NULL);
1484 std::vector<unsigned int> edges_without_neighbors;
1486 for (
unsigned int i = 0; i < edge_nodes.size(); ++i)
1489 std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1491 std::vector<PenetrationInfo *> face_info_comm_edge;
1493 &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1495 if (face_info_comm_edge.size() == 0)
1496 edges_without_neighbors.push_back(i);
1497 else if (face_info_comm_edge.size() > 1)
1498 mooseError(
"Only one neighbor allowed per edge");
1500 edge_face_info[i] = face_info_comm_edge[0];
1504 std::vector<unsigned int>::reverse_iterator rit;
1505 for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1507 unsigned int index = *rit;
1508 edge_nodes.erase(edge_nodes.begin() + index);
1509 edge_face_weights.erase(edge_face_weights.begin() + index);
1510 edge_face_info.erase(edge_face_info.begin() + index);
1514 if (edge_nodes.size() > 1)
1516 if (edge_nodes.size() != 2)
1517 mooseError(
"Invalid number of smoothing edges");
1520 std::vector<const Node *> common_nodes;
1521 std::set_intersection(edge_nodes[0].begin(),
1522 edge_nodes[0].end(),
1523 edge_nodes[1].begin(),
1524 edge_nodes[1].end(),
1525 std::inserter(common_nodes, common_nodes.end()));
1527 if (common_nodes.size() != 1)
1528 mooseError(
"Invalid number of common nodes");
1530 for (
const auto & pinfo : edge_face_info)
1531 elems_to_exclude.insert(pinfo->_elem->id());
1533 std::vector<PenetrationInfo *> face_info_comm_edge;
1535 &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1537 unsigned int num_corner_neighbors = face_info_comm_edge.size();
1539 if (num_corner_neighbors > 0)
1541 Real fw0 = edge_face_weights[0];
1542 Real fw1 = edge_face_weights[1];
1545 Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1548 edge_face_weights[0] *= (1.0 - fw1);
1549 edge_face_weights[1] *= (1.0 - fw0);
1551 for (
unsigned int i = 0; i < num_corner_neighbors; ++i)
1553 edge_face_weights.push_back(fw_corner);
1554 edge_face_info.push_back(face_info_comm_edge[i]);
1564 std::vector<std::vector<const Node *>> & edge_nodes,
1565 std::vector<Real> & edge_face_weights)
1568 const Real & xi = p(0);
1569 const Real & eta = p(1);
1579 if (xi < -smooth_limit)
1581 std::vector<const Node *> en;
1583 edge_nodes.push_back(en);
1587 edge_face_weights.push_back(fw);
1589 else if (xi > smooth_limit)
1591 std::vector<const Node *> en;
1593 edge_nodes.push_back(en);
1597 edge_face_weights.push_back(fw);
1606 if (eta < -smooth_limit)
1608 std::vector<const Node *> en;
1611 edge_nodes.push_back(en);
1615 edge_face_weights.push_back(fw);
1617 if ((xi + eta) > smooth_limit)
1619 std::vector<const Node *> en;
1622 edge_nodes.push_back(en);
1626 edge_face_weights.push_back(fw);
1628 if (xi < -smooth_limit)
1630 std::vector<const Node *> en;
1633 edge_nodes.push_back(en);
1637 edge_face_weights.push_back(fw);
1646 if (eta < -smooth_limit)
1648 std::vector<const Node *> en;
1651 edge_nodes.push_back(en);
1655 edge_face_weights.push_back(fw);
1657 if (xi > smooth_limit)
1659 std::vector<const Node *> en;
1662 edge_nodes.push_back(en);
1666 edge_face_weights.push_back(fw);
1668 if (eta > smooth_limit)
1670 std::vector<const Node *> en;
1673 edge_nodes.push_back(en);
1677 edge_face_weights.push_back(fw);
1679 if (xi < -smooth_limit)
1681 std::vector<const Node *> en;
1684 edge_nodes.push_back(en);
1688 edge_face_weights.push_back(fw);
1703 const Node * secondary_node,
1704 const std::set<dof_id_type> & elems_to_exclude,
1705 const std::vector<const Node *> edge_nodes,
1706 std::vector<PenetrationInfo *> & face_info_comm_edge,
1707 std::vector<PenetrationInfo *> & p_info)
1713 mooseAssert(node_to_elem_pair !=
_node_to_elem_map.end(),
"Missing entry in node to elem map");
1714 const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1716 std::vector<const Elem *> elems_connected_to_edge;
1718 for (
unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1720 if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1724 std::vector<const Node *> nodevec;
1725 for (
unsigned int ni = 0; ni < elem->
n_nodes(); ++ni)
1727 nodevec.push_back(elem->
node_ptr(ni));
1729 std::vector<const Node *> common_nodes;
1730 std::sort(nodevec.begin(), nodevec.end());
1731 std::set_intersection(edge_nodes.begin(),
1735 std::inserter(common_nodes, common_nodes.end()));
1737 if (common_nodes.size() == edge_nodes.size())
1738 elems_connected_to_edge.push_back(elem);
1741 if (elems_connected_to_edge.size() > 0)
1751 bool allowMultipleNeighbors =
false;
1753 if (elems_connected_to_edge[0]->
dim() == 3)
1755 if (edge_nodes.size() == 1)
1757 allowMultipleNeighbors =
true;
1761 for (
unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1763 std::vector<PenetrationInfo *> thisElemInfo;
1764 getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1765 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1767 if (thisElemInfo.size() > 1)
1769 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1770 face_info_comm_edge.push_back(thisElemInfo[0]);
1775 thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1776 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1778 if (thisElemInfo.size() > 1)
1780 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1781 face_info_comm_edge.push_back(thisElemInfo[0]);
1785 for (
unsigned int j = 0; j < thisElemInfo.size(); ++j)
1786 face_info_comm_edge.push_back(thisElemInfo[j]);
1793 std::vector<PenetrationInfo *> & p_info,
1796 for (
const auto &
pi : p_info)
1801 if (
pi->_elem == elem)
1802 thisElemInfo.push_back(
pi);
1808 std::vector<PenetrationInfo *> & p_info,
1809 const Node * secondary_node,
1811 const std::vector<const Node *> & nodes_that_must_be_on_side,
1812 const bool check_whether_reasonable)
1822 bool already_have_info_this_side =
false;
1823 for (
const auto &
pi : thisElemInfo)
1824 if (
pi->_side_num == s)
1826 already_have_info_this_side =
true;
1830 if (already_have_info_this_side)
1837 std::vector<const Node *> nodevec;
1838 for (
unsigned int ni = 0; ni < side->
n_nodes(); ++ni)
1839 nodevec.push_back(side->
node_ptr(ni));
1841 std::sort(nodevec.begin(), nodevec.end());
1842 std::vector<const Node *> common_nodes;
1843 std::set_intersection(nodes_that_must_be_on_side.begin(),
1844 nodes_that_must_be_on_side.end(),
1847 std::inserter(common_nodes, common_nodes.end()));
1848 if (common_nodes.size() != nodes_that_must_be_on_side.size())
1859 if (check_whether_reasonable)
1868 Point contact_on_face_ref;
1870 Real tangential_distance = 0.;
1872 bool contact_point_on_side;
1873 std::vector<const Node *> off_edge_nodes;
1874 std::vector<std::vector<Real>> side_phi;
1875 std::vector<std::vector<RealGradient>> side_grad_phi;
1876 std::vector<RealGradient> dxyzdxi;
1877 std::vector<RealGradient> dxyzdeta;
1878 std::vector<RealGradient> d2xyzdxideta;
1880 std::unique_ptr<PenetrationInfo> pen_info =
1881 std::make_unique<PenetrationInfo>(elem,
1886 tangential_distance,
1889 contact_on_face_ref,
1897 bool search_succeeded =
false;
1905 contact_point_on_side,
1909 if (search_succeeded)
1911 thisElemInfo.push_back(pen_info.get());
1912 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