20 #include "libmesh/threads.h" 32 std::map<dof_id_type, PenetrationInfo *> & penetration_info,
33 bool check_whether_reasonable,
35 Real tangential_tolerance,
36 bool do_normal_smoothing,
37 Real normal_smoothing_distance,
39 std::vector<std::vector<FEBase *>> & fes,
42 const std::map<
dof_id_type, std::vector<dof_id_type>> & node_to_elem_map,
43 const std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> & bc_tuples)
44 : _subproblem(subproblem),
46 _primary_boundary(primary_boundary),
47 _secondary_boundary(secondary_boundary),
48 _penetration_info(penetration_info),
49 _check_whether_reasonable(check_whether_reasonable),
50 _update_location(update_location),
51 _tangential_tolerance(tangential_tolerance),
52 _do_normal_smoothing(do_normal_smoothing),
53 _normal_smoothing_distance(normal_smoothing_distance),
54 _normal_smoothing_method(normal_smoothing_method),
55 _nodal_normal_x(NULL),
56 _nodal_normal_y(NULL),
57 _nodal_normal_z(NULL),
60 _nearest_node(nearest_node),
61 _node_to_elem_map(node_to_elem_map),
68 : _subproblem(x._subproblem),
70 _primary_boundary(x._primary_boundary),
71 _secondary_boundary(x._secondary_boundary),
72 _penetration_info(x._penetration_info),
73 _check_whether_reasonable(x._check_whether_reasonable),
74 _update_location(x._update_location),
75 _tangential_tolerance(x._tangential_tolerance),
76 _do_normal_smoothing(x._do_normal_smoothing),
77 _normal_smoothing_distance(x._normal_smoothing_distance),
78 _normal_smoothing_method(x._normal_smoothing_method),
81 _nearest_node(x._nearest_node),
82 _node_to_elem_map(x._node_to_elem_map),
83 _bc_tuples(x._bc_tuples)
102 for (
const auto & node_id : range)
113 std::vector<PenetrationInfo *> p_info;
114 bool info_set(
false);
124 const Point contact_ref =
info->_closest_point_ref;
125 bool contact_point_on_side(
false);
129 std::vector<Point> points(1);
130 points[0] = contact_ref;
131 const std::vector<Point> & secondary_pos = fe_side->get_xyz();
140 contact_point_on_side);
143 info->_closest_point_ref = contact_ref;
146 info->_distance = 0.0;
151 Real old_tangential_distance(
info->_tangential_distance);
152 bool contact_point_on_side(
false);
161 contact_point_on_side);
163 if (contact_point_on_side)
165 if (
info->_tangential_distance <= 0.0)
169 else if (
info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
171 if (
info->_side->dim() == 2 &&
info->_off_edge_nodes.size() < 2)
189 "Missing entry in node to elem map");
190 const std::vector<dof_id_type> & closest_elems = node_to_elem_pair->second;
192 for (
const auto & elem_id : closest_elems)
196 std::vector<PenetrationInfo *> thisElemInfo;
197 std::vector<const Node *> nodesThatMustBeOnSide;
198 nodesThatMustBeOnSide.push_back(closest_node);
203 if (p_info.size() == 1)
211 else if (p_info.size() > 1)
215 std::vector<RidgeData> ridgeDataVec;
216 for (
unsigned int i = 0; i + 1 < p_info.size(); ++i)
217 for (
unsigned int j = i + 1; j < p_info.size(); ++j)
220 Real tangential_distance(0.0);
221 const Node * closest_node_on_ridge = NULL;
222 unsigned int index = 0;
223 Point closest_coor_ref;
226 closest_node_on_ridge,
232 if (found_ridge_contact_point)
240 ridgeDataVec.push_back(rpd);
244 if (ridgeDataVec.size() > 0)
248 std::vector<RidgeSetData> ridgeSetDataVec;
249 for (
unsigned int i = 0; i < ridgeDataVec.size(); ++i)
251 bool foundSetWithMatchingNode =
false;
252 for (
unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
254 if (ridgeDataVec[i]._closest_node != NULL &&
255 ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
257 foundSetWithMatchingNode =
true;
258 ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
262 if (!foundSetWithMatchingNode)
268 ridgeSetDataVec.push_back(rsd);
272 for (
unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
274 if (ridgeSetDataVec[i]._closest_node !=
277 if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1)
279 if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
282 ridgeSetDataVec[i]._closest_coor =
283 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
284 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
285 ridgeSetDataVec[i]._distance = contact_point_vec.norm();
291 ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
292 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
293 ridgeSetDataVec[i]._distance = contact_point_vec.norm();
298 ridgeSetDataVec[i]._closest_coor =
299 ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
300 Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
301 ridgeSetDataVec[i]._distance = contact_point_vec.norm();
305 unsigned int closest_ridge_set_index(0);
306 Real closest_distance(ridgeSetDataVec[0]._distance);
307 Point closest_point(ridgeSetDataVec[0]._closest_coor);
308 for (
unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
310 if (ridgeSetDataVec[i]._distance < closest_distance)
312 closest_ridge_set_index = i;
313 closest_distance = ridgeSetDataVec[i]._distance;
314 closest_point = ridgeSetDataVec[i]._closest_coor;
318 if (closest_distance <
324 for (
unsigned int i = 0;
325 i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
328 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
329 face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
333 "face_index invalid");
335 p_info[face_index]->_closest_point = closest_point;
336 p_info[face_index]->_distance =
337 (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
343 Point normal(closest_point - node);
344 const Real len(normal.norm());
349 const Real dot(normal * p_info[face_index]->_normal);
354 p_info[face_index]->_normal = normal;
356 p_info[face_index]->_tangential_distance = 0.0;
358 Point closest_point_ref;
359 if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
362 p_info[face_index]->_tangential_distance =
363 ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._tangential_distance;
364 p_info[face_index]->_closest_point_ref =
365 ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
369 const Node * closest_node_on_face;
371 closest_node_on_face,
372 p_info[face_index]->_side);
375 if (closest_node_on_face != ridgeSetDataVec[closest_ridge_set_index]._closest_node)
377 mooseError(
"Closest node when restricting point to face != closest node from " 384 std::vector<Point> points(1);
385 points[0] = p_info[face_index]->_closest_point_ref;
386 fe->reinit(p_info[face_index]->_side, &points);
387 p_info[face_index]->_side_phi = fe->get_phi();
388 p_info[face_index]->_side_grad_phi = fe->get_dphi();
389 p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
390 p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
391 p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
403 unsigned int best(0), i(1);
421 }
while (i < p_info.size() && best < p_info.size());
422 if (best < p_info.size())
450 for (
unsigned int j = 0; j < p_info.size(); ++j)
472 mooseAssert(infoNew != NULL,
"infoNew object is null");
599 const std::vector<const Node *> & off_edge_nodes1 = pi1->
_off_edge_nodes;
600 const std::vector<const Node *> & off_edge_nodes2 = pi2->
_off_edge_nodes;
601 const unsigned dim1 = pi1->
_side->dim();
605 mooseAssert(pi2->
_side->dim() == 1,
"Incompatible dimensions.");
606 mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
607 "off_edge_nodes size should be <2 for 2D contact");
608 if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
609 off_edge_nodes1[0] == off_edge_nodes2[0])
614 mooseAssert(dim1 == 2 && pi2->
_side->dim() == 2,
"Incompatible dimensions.");
615 mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
616 "off_edge_nodes size should be <3 for 3D contact");
617 if (off_edge_nodes1.size() == 1)
619 if (off_edge_nodes2.size() == 1)
621 if (off_edge_nodes1[0] == off_edge_nodes2[0])
624 else if (off_edge_nodes2.size() == 2)
626 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
630 else if (off_edge_nodes1.size() == 2)
632 if (off_edge_nodes2.size() == 1)
634 if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
637 else if (off_edge_nodes2.size() == 2)
639 if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
640 off_edge_nodes1[1] == off_edge_nodes2[1]) ||
641 (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
651 Real & tangential_distance,
652 const Node *& closest_node,
653 unsigned int & index,
654 Point & contact_point_ref,
655 std::vector<PenetrationInfo *> & p_info,
656 const unsigned int index1,
657 const unsigned int index2)
659 tangential_distance = 0.0;
663 const unsigned sidedim(pi1->
_side->dim());
664 mooseAssert(sidedim == pi2->
_side->dim(),
"Incompatible dimensionalities");
667 std::vector<const Node *> side1_nodes;
669 std::vector<const Node *> side2_nodes;
672 std::sort(side1_nodes.begin(), side1_nodes.end());
673 std::sort(side2_nodes.begin(), side2_nodes.end());
676 std::vector<const Node *> common_nodes;
677 std::set_intersection(side1_nodes.begin(),
681 std::inserter(common_nodes, common_nodes.end()));
683 if (common_nodes.size() != sidedim)
686 bool found_point1, found_point2;
688 const Node * closest_node1;
690 closest_coor_ref1, closest_node1, pi1->
_side, common_nodes);
693 const Node * closest_node2;
695 closest_coor_ref2, closest_node2, pi2->
_side, common_nodes);
697 if (!found_point1 || !found_point2)
710 std::vector<Point> points(1);
717 contact_point_ref = closest_coor_ref1;
718 points[0] = closest_coor_ref1;
719 fe->reinit(pi1->
_side, &points);
725 contact_point_ref = closest_coor_ref2;
726 points[0] = closest_coor_ref2;
727 fe->reinit(pi2->
_side, &points);
731 contact_point = fe->get_xyz()[0];
737 mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
738 "If off edge of ridge, closest node must be the same on both elements");
739 closest_node = closest_node1;
742 tangential_distance = off_face.
norm();
753 corner_nodes.clear();
755 corner_nodes.push_back(side->node_ptr(0));
756 corner_nodes.push_back(side->node_ptr(1));
770 corner_nodes.push_back(side->node_ptr(2));
778 corner_nodes.push_back(side->node_ptr(2));
779 corner_nodes.push_back(side->node_ptr(3));
793 const Node *& closest_node,
795 const std::vector<const Node *> & edge_nodes)
802 std::vector<unsigned int> local_node_indices;
803 for (
const auto & edge_node : edge_nodes)
805 unsigned int local_index = side->get_node_index(edge_node);
808 local_node_indices.push_back(local_index);
810 mooseAssert(local_node_indices.size() == side->dim(),
811 "Number of edge nodes must match side dimensionality");
812 std::sort(local_node_indices.begin(), local_node_indices.end());
814 bool off_of_this_edge =
false;
822 if (local_node_indices[0] == 0)
827 off_of_this_edge =
true;
828 closest_node = side->node_ptr(0);
831 else if (local_node_indices[0] == 1)
836 off_of_this_edge =
true;
837 closest_node = side->node_ptr(1);
851 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
856 off_of_this_edge =
true;
858 closest_node = side->node_ptr(0);
860 closest_node = side->node_ptr(1);
863 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
865 if ((xi + eta) > 1.0)
867 Real delta = (xi + eta - 1.0) / 2.0;
870 off_of_this_edge =
true;
872 closest_node = side->node_ptr(1);
874 closest_node = side->node_ptr(2);
877 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
882 off_of_this_edge =
true;
884 closest_node = side->node_ptr(2);
886 closest_node = side->node_ptr(0);
901 if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
906 off_of_this_edge =
true;
908 closest_node = side->node_ptr(0);
910 closest_node = side->node_ptr(1);
913 else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
918 off_of_this_edge =
true;
920 closest_node = side->node_ptr(1);
922 closest_node = side->node_ptr(2);
925 else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
930 off_of_this_edge =
true;
932 closest_node = side->node_ptr(3);
934 closest_node = side->node_ptr(2);
937 else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
942 off_of_this_edge =
true;
944 closest_node = side->node_ptr(0);
946 closest_node = side->node_ptr(3);
962 return off_of_this_edge;
973 bool off_of_this_face(
false);
984 off_of_this_face =
true;
985 closest_node = side->node_ptr(0);
990 off_of_this_face =
true;
991 closest_node = side->node_ptr(1);
1003 off_of_this_face =
true;
1006 closest_node = side->node_ptr(0);
1012 closest_node = side->node_ptr(1);
1017 else if ((xi + eta) > 1.0)
1019 Real delta = (xi + eta - 1.0) / 2.0;
1022 off_of_this_face =
true;
1025 closest_node = side->node_ptr(1);
1034 closest_node = side->node_ptr(2);
1045 off_of_this_face =
true;
1048 closest_node = side->node_ptr(2);
1054 closest_node = side->node_ptr(0);
1069 off_of_this_face =
true;
1072 closest_node = side->node_ptr(0);
1078 closest_node = side->node_ptr(1);
1086 off_of_this_face =
true;
1089 closest_node = side->node_ptr(1);
1095 closest_node = side->node_ptr(2);
1103 off_of_this_face =
true;
1106 closest_node = side->node_ptr(3);
1112 closest_node = side->node_ptr(2);
1120 off_of_this_face =
true;
1123 closest_node = side->node_ptr(0);
1129 closest_node = side->node_ptr(3);
1143 return off_of_this_face;
1150 const Point * secondary_point,
1151 const Real tangential_tolerance)
1153 unsigned int dim = primary_elem->dim();
1155 const std::vector<Point> & phys_point = fe->get_xyz();
1157 const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1158 const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1162 std::vector<Point> points(1);
1164 fe->reinit(side, &points);
1168 const Real twosqrt2 = 2.8284;
1169 Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1174 normal = dxyz_dxi[0].
cross(dxyz_deta[0]);
1176 else if (
dim - 1 == 1)
1178 const Node *
const * elem_nodes = primary_elem->get_nodes();
1179 const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1180 const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1182 Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1183 out_of_plane_normal /= out_of_plane_normal.norm();
1185 normal = dxyz_dxi[0].
cross(out_of_plane_normal);
1191 normal /= normal.
norm();
1193 const Real dot(d * normal);
1198 const Real tangdist = tangcomp.
norm();
1202 const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.
norm());
1204 bool isReasonableCandidate =
true;
1205 if (tangdist > faceExpansionFactor * max_face_length)
1207 isReasonableCandidate =
false;
1209 return isReasonableCandidate;
1217 std::vector<Point> points(1);
1218 points[0] =
info._starting_closest_point_ref;
1220 fe.reinit(&side, &points);
1221 const std::vector<Point> & starting_point = fe.get_xyz();
1222 info._incremental_slip =
info._closest_point - starting_point[0];
1223 if (
info.isCaptured())
1225 info._frictional_energy =
1226 info._frictional_energy_old +
info._contact_force *
info._incremental_slip;
1227 info._accumulated_slip =
info._accumulated_slip_old +
info._incremental_slip.norm();
1233 std::vector<PenetrationInfo *> & p_info,
1242 std::vector<Real> edge_face_weights;
1243 std::vector<PenetrationInfo *> edge_face_info;
1247 mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1248 "edge_face_info.size() != edge_face_weights.size()");
1250 if (edge_face_info.size() > 0)
1254 Real this_face_weight = 1.0;
1256 for (
unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1260 new_normal += npi->
_normal * edge_face_weights[efwi];
1262 this_face_weight -= edge_face_weights[efwi];
1264 mooseAssert(this_face_weight >= (0.25 - 1e-8),
1265 "Sum of weights of other faces shouldn't exceed 0.75");
1266 new_normal +=
info->_normal * this_face_weight;
1268 const Real len = new_normal.
norm();
1272 info->_normal = new_normal;
1282 const Real len(
info->_normal.norm());
1284 info->_normal /= len;
1291 std::vector<PenetrationInfo *> & edge_face_info,
1292 std::vector<Real> & edge_face_weights,
1293 std::vector<PenetrationInfo *> & p_info,
1294 const Node & secondary_node)
1296 const Elem * side =
info->_side;
1297 const Point & p =
info->_closest_point_ref;
1298 std::set<dof_id_type> elems_to_exclude;
1299 elems_to_exclude.insert(
info->_elem->id());
1301 std::vector<std::vector<const Node *>> edge_nodes;
1305 std::vector<Elem *> edge_neighbor_elems;
1306 edge_face_info.resize(edge_nodes.size(), NULL);
1308 std::vector<unsigned int> edges_without_neighbors;
1310 for (
unsigned int i = 0; i < edge_nodes.size(); ++i)
1313 std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1315 std::vector<PenetrationInfo *> face_info_comm_edge;
1317 &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1319 if (face_info_comm_edge.size() == 0)
1320 edges_without_neighbors.push_back(i);
1321 else if (face_info_comm_edge.size() > 1)
1322 mooseError(
"Only one neighbor allowed per edge");
1324 edge_face_info[i] = face_info_comm_edge[0];
1328 std::vector<unsigned int>::reverse_iterator rit;
1329 for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1331 unsigned int index = *rit;
1332 edge_nodes.erase(edge_nodes.begin() + index);
1333 edge_face_weights.erase(edge_face_weights.begin() + index);
1334 edge_face_info.erase(edge_face_info.begin() + index);
1338 if (edge_nodes.size() > 1)
1340 if (edge_nodes.size() != 2)
1341 mooseError(
"Invalid number of smoothing edges");
1344 std::vector<const Node *> common_nodes;
1345 std::set_intersection(edge_nodes[0].begin(),
1346 edge_nodes[0].end(),
1347 edge_nodes[1].begin(),
1348 edge_nodes[1].end(),
1349 std::inserter(common_nodes, common_nodes.end()));
1351 if (common_nodes.size() != 1)
1352 mooseError(
"Invalid number of common nodes");
1354 for (
const auto & pinfo : edge_face_info)
1355 elems_to_exclude.insert(pinfo->_elem->id());
1357 std::vector<PenetrationInfo *> face_info_comm_edge;
1359 &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1361 unsigned int num_corner_neighbors = face_info_comm_edge.size();
1363 if (num_corner_neighbors > 0)
1365 Real fw0 = edge_face_weights[0];
1366 Real fw1 = edge_face_weights[1];
1369 Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1372 edge_face_weights[0] *= (1.0 - fw1);
1373 edge_face_weights[1] *= (1.0 - fw0);
1375 for (
unsigned int i = 0; i < num_corner_neighbors; ++i)
1377 edge_face_weights.push_back(fw_corner);
1378 edge_face_info.push_back(face_info_comm_edge[i]);
1388 std::vector<std::vector<const Node *>> & edge_nodes,
1389 std::vector<Real> & edge_face_weights)
1392 const Real & xi = p(0);
1393 const Real & eta = p(1);
1403 if (xi < -smooth_limit)
1405 std::vector<const Node *> en;
1406 en.push_back(side->node_ptr(0));
1407 edge_nodes.push_back(en);
1411 edge_face_weights.push_back(fw);
1413 else if (xi > smooth_limit)
1415 std::vector<const Node *> en;
1416 en.push_back(side->node_ptr(1));
1417 edge_nodes.push_back(en);
1421 edge_face_weights.push_back(fw);
1430 if (eta < -smooth_limit)
1432 std::vector<const Node *> en;
1433 en.push_back(side->node_ptr(0));
1434 en.push_back(side->node_ptr(1));
1435 edge_nodes.push_back(en);
1439 edge_face_weights.push_back(fw);
1441 if ((xi + eta) > smooth_limit)
1443 std::vector<const Node *> en;
1444 en.push_back(side->node_ptr(1));
1445 en.push_back(side->node_ptr(2));
1446 edge_nodes.push_back(en);
1450 edge_face_weights.push_back(fw);
1452 if (xi < -smooth_limit)
1454 std::vector<const Node *> en;
1455 en.push_back(side->node_ptr(2));
1456 en.push_back(side->node_ptr(0));
1457 edge_nodes.push_back(en);
1461 edge_face_weights.push_back(fw);
1470 if (eta < -smooth_limit)
1472 std::vector<const Node *> en;
1473 en.push_back(side->node_ptr(0));
1474 en.push_back(side->node_ptr(1));
1475 edge_nodes.push_back(en);
1479 edge_face_weights.push_back(fw);
1481 if (xi > smooth_limit)
1483 std::vector<const Node *> en;
1484 en.push_back(side->node_ptr(1));
1485 en.push_back(side->node_ptr(2));
1486 edge_nodes.push_back(en);
1490 edge_face_weights.push_back(fw);
1492 if (eta > smooth_limit)
1494 std::vector<const Node *> en;
1495 en.push_back(side->node_ptr(2));
1496 en.push_back(side->node_ptr(3));
1497 edge_nodes.push_back(en);
1501 edge_face_weights.push_back(fw);
1503 if (xi < -smooth_limit)
1505 std::vector<const Node *> en;
1506 en.push_back(side->node_ptr(3));
1507 en.push_back(side->node_ptr(0));
1508 edge_nodes.push_back(en);
1512 edge_face_weights.push_back(fw);
1527 const Node * secondary_node,
1528 const std::set<dof_id_type> & elems_to_exclude,
1529 const std::vector<const Node *> edge_nodes,
1530 std::vector<PenetrationInfo *> & face_info_comm_edge,
1531 std::vector<PenetrationInfo *> & p_info)
1537 mooseAssert(node_to_elem_pair !=
_node_to_elem_map.end(),
"Missing entry in node to elem map");
1538 const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1540 std::vector<const Elem *> elems_connected_to_edge;
1542 for (
unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1544 if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1546 const Elem * elem =
_mesh.
elemPtr(elems_connected_to_node[ecni]);
1548 std::vector<const Node *> nodevec;
1549 for (
unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1550 if (elem->is_vertex(ni))
1551 nodevec.push_back(elem->node_ptr(ni));
1553 std::vector<const Node *> common_nodes;
1554 std::sort(nodevec.begin(), nodevec.end());
1555 std::set_intersection(edge_nodes.begin(),
1559 std::inserter(common_nodes, common_nodes.end()));
1561 if (common_nodes.size() == edge_nodes.size())
1562 elems_connected_to_edge.push_back(elem);
1565 if (elems_connected_to_edge.size() > 0)
1575 bool allowMultipleNeighbors =
false;
1577 if (elems_connected_to_edge[0]->
dim() == 3)
1579 if (edge_nodes.size() == 1)
1581 allowMultipleNeighbors =
true;
1585 for (
unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1587 std::vector<PenetrationInfo *> thisElemInfo;
1588 getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1589 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1591 if (thisElemInfo.size() > 1)
1593 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1594 face_info_comm_edge.push_back(thisElemInfo[0]);
1599 thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1600 if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1602 if (thisElemInfo.size() > 1)
1604 "Found multiple neighbors to current edge/face on surface when only one is allowed");
1605 face_info_comm_edge.push_back(thisElemInfo[0]);
1609 for (
unsigned int j = 0; j < thisElemInfo.size(); ++j)
1610 face_info_comm_edge.push_back(thisElemInfo[j]);
1617 std::vector<PenetrationInfo *> & p_info,
1620 for (
const auto &
pi : p_info)
1625 if (
pi->_elem == elem)
1626 thisElemInfo.push_back(
pi);
1632 std::vector<PenetrationInfo *> & p_info,
1633 const Node * secondary_node,
1635 const std::vector<const Node *> & nodes_that_must_be_on_side,
1636 const bool check_whether_reasonable)
1638 std::vector<unsigned int> sides;
1649 for (
unsigned int i = 0; i < sides.size(); ++i)
1652 bool already_have_info_this_side =
false;
1653 for (
const auto &
pi : thisElemInfo)
1654 if (
pi->_side_num == sides[i])
1656 already_have_info_this_side =
true;
1660 if (already_have_info_this_side)
1663 const Elem * side = (elem->build_side_ptr(sides[i],
false)).release();
1667 std::vector<const Node *> nodevec;
1668 for (
unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1669 nodevec.push_back(side->node_ptr(ni));
1671 std::sort(nodevec.begin(), nodevec.end());
1672 std::vector<const Node *> common_nodes;
1673 std::set_intersection(nodes_that_must_be_on_side.begin(),
1674 nodes_that_must_be_on_side.end(),
1677 std::inserter(common_nodes, common_nodes.end()));
1678 if (common_nodes.size() != nodes_that_must_be_on_side.size())
1689 if (check_whether_reasonable)
1698 Point contact_on_face_ref;
1700 Real tangential_distance = 0.;
1702 bool contact_point_on_side;
1703 std::vector<const Node *> off_edge_nodes;
1704 std::vector<std::vector<Real>> side_phi;
1705 std::vector<std::vector<RealGradient>> side_grad_phi;
1706 std::vector<RealGradient> dxyzdxi;
1707 std::vector<RealGradient> dxyzdeta;
1708 std::vector<RealGradient> d2xyzdxideta;
1715 tangential_distance,
1718 contact_on_face_ref,
1733 contact_point_on_side);
1735 thisElemInfo.push_back(pen_info);
1737 p_info.push_back(pen_info);
1745 const Elem *
const elem)
1753 return std::get<0>(tup) <
id;
1757 return id < std::get<0>(tup);
1763 for (
auto & t = range.first; t != range.second; ++t)
1765 sides.push_back(std::get<1>(*t));
void getSidesOnPrimaryBoundary(std::vector< unsigned int > &sides, const Elem *const elem)
MooseVariable * _nodal_normal_z
std::tuple< dof_id_type, unsigned short int, boundary_id_type > BCTuple
void getSmoothingEdgeNodesAndWeights(const Point &p, const Elem *side, std::vector< std::vector< const Node *>> &edge_nodes, std::vector< Real > &edge_face_weights)
auto norm() const -> decltype(std::norm(Real()))
const unsigned int invalid_uint
Real _accumulated_slip_old
virtual Elem * elemPtr(const dof_id_type i)
OutputType getValue(const Elem *elem, const std::vector< std::vector< OutputShape >> &phi) const
Compute the variable value at a point on an element.
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
StoredRange< std::vector< dof_id_type >::iterator, dof_id_type > NodeIdRange
Finds the nearest node to each node in boundary1 to each node in boundary2 and the other way around...
Real _normal_smoothing_distance
const Elem * _starting_elem
bool _check_whether_reasonable
const std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > & _bc_tuples
bool findRidgeContactPoint(Point &contact_point, Real &tangential_distance, const Node *&closest_node, unsigned int &index, Point &contact_point_ref, std::vector< PenetrationInfo *> &p_info, const unsigned int index1, const unsigned int index2)
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
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)
bool restrictPointToFace(Point &p, const Node *&closest_node, const Elem *side)
unsigned int _stick_locked_this_step
void getSideCornerNodes(const Elem *side, std::vector< const Node *> &corner_nodes)
Threads::spin_mutex pinfo_mutex
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)
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
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Real _tangential_distance
std::vector< std::vector< FEBase * > > & _fes
void findContactPoint(PenetrationInfo &p_info, FEBase *fe_elem, FEBase *fe_side, FEType &fe_side_type, const Point &secondary_point, bool start_with_centroid, const Real tangential_tolerance, bool &contact_point_on_side)
Finds the closest point (called the contact point) on the primary_elem on side "side" to the secondar...
bool restrictPointToSpecifiedEdgeOfFace(Point &p, const Node *&closest_node, const Elem *side, const std::vector< const Node *> &edge_nodes)
bool isFaceReasonableCandidate(const Elem *primary_elem, const Elem *side, FEBase *fe, const Point *secondary_point, const Real tangential_tolerance)
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
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
const Node * _closest_node
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.
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
FEGenericBase< Real > FEBase
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
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
ElemSideBuilder _elem_side_builder
Helper for building element sides without extraneous allocation.
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...
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, std::vector< std::vector< FEBase *>> &fes, FEType &fe_type, NearestNodeLocator &nearest_node, const std::map< dof_id_type, std::vector< dof_id_type >> &node_to_elem_map, const std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type >> &bc_tuples)
MooseVariable * _nodal_normal_y
void getInfoForElem(std::vector< PenetrationInfo *> &thisElemInfo, std::vector< PenetrationInfo *> &p_info, const Elem *elem)
bool _do_normal_smoothing
void computeSlip(FEBase &fe, PenetrationInfo &info)
MECH_STATUS_ENUM _mech_status_old
const std::map< dof_id_type, std::vector< dof_id_type > > & _node_to_elem_map
void switchInfo(PenetrationInfo *&info, PenetrationInfo *&infoNew)
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...
CompeteInteractionResult competeInteractionsBothOnFace(PenetrationInfo *pi1, PenetrationInfo *pi2)
Determine whether first (pi1) or second (pi2) interaction is stronger when it is known that the node ...
CommonEdgeResult interactionsOffCommonEdge(PenetrationInfo *pi1, PenetrationInfo *pi2)
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method
std::vector< RidgeData > _ridge_data_vec