https://mooseframework.inl.gov
PenetrationThread.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 // Moose
11 #include "PenetrationThread.h"
12 #include "ParallelUniqueId.h"
13 #include "FindContactPoint.h"
14 #include "NearestNodeLocator.h"
15 #include "SubProblem.h"
16 #include "MooseVariableFE.h"
17 #include "MooseMesh.h"
18 #include "MooseUtils.h"
19 
20 #include "libmesh/threads.h"
21 
22 #include <algorithm>
23 
24 using namespace libMesh;
25 
26 // Anonymous namespace for helper functions that ought to be moved
27 // into libMesh
28 namespace
29 {
30 Point
31 closest_point_to_edge(const Point & src, const Point & p0, const Point & p1)
32 {
33  const Point line01 = p1 - p0;
34  const Real line0c_xi = ((src - p0) * line01) / line01.norm_sq();
35  // The projection would be behind p0; p0 is closest
36  if (line0c_xi <= 0)
37  return p0;
38  // The projection would be past p1; p1 is closest
39  if (line0c_xi >= 1)
40  return p1;
41  // The projection is on the segment between p0 to p1.
42  return p0 + line0c_xi * line01;
43 }
44 
45 Point
46 closest_point_to_side(const Point & src, const Elem & side)
47 {
48  switch (side.type())
49  {
50  case EDGE2:
51  case EDGE3:
52  case EDGE4:
53  mooseAssert(side.has_affine_map(),
54  "Penetration of elements with curved sides not implemented");
55  return closest_point_to_edge(src, side.point(0), side.point(1));
56  case TRI3:
57  case TRI6:
58  {
59  mooseAssert(side.has_affine_map(),
60  "Penetration of elements with curved sides not implemented");
61  const Point p0 = side.point(0), p1 = side.point(1), p2 = side.point(2);
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;
67  // If we're outside the triangle past line 01, our closest point
68  // is on that line.
69  if (planar_offset.cross(l01) * tri_normal > 0)
70  return closest_point_to_edge(src, p0, p1);
71  // If we're outside the triangle past line 02, our closest point
72  // is on that line.
73  if (planar_offset.cross(l02) * tri_normal < 0)
74  return closest_point_to_edge(src, p0, p2);
75  // If we're outside the triangle past line 12, our closest point
76  // is on that line.
77  if ((in_plane - p1).cross(p2 - p1) * tri_normal > 0)
78  return closest_point_to_edge(src, p1, p2);
79  // We must be inside the triangle!
80  return in_plane;
81  }
82  case QUAD4:
83  case QUAD8:
84  case QUAD9:
85  case C0POLYGON:
86  mooseError("Not implemented");
87  default:
88  mooseError("Side type not recognized");
89  break;
90  }
91 }
92 
93 } // anonymous namespace
94 
95 // Mutex to use when accessing _penetration_info;
97 
99  SubProblem & subproblem,
100  const MooseMesh & mesh,
101  BoundaryID primary_boundary,
102  BoundaryID secondary_boundary,
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,
109  PenetrationLocator::NORMAL_SMOOTHING_METHOD normal_smoothing_method,
110  bool use_point_locator,
111  std::vector<std::vector<FEBase *>> & fes,
112  FEType & fe_type,
113  NearestNodeLocator & nearest_node,
114  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map)
115  : _subproblem(subproblem),
116  _mesh(mesh),
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),
130  _fes(fes),
131  _fe_type(fe_type),
132  _nearest_node(nearest_node),
133  _node_to_elem_map(node_to_elem_map)
134 {
135 }
136 
137 // Splitting Constructor
139  : _subproblem(x._subproblem),
140  _mesh(x._mesh),
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),
151  _fes(x._fes),
152  _fe_type(x._fe_type),
153  _nearest_node(x._nearest_node),
154  _node_to_elem_map(x._node_to_elem_map)
155 {
156 }
157 
158 void
160 {
161  ParallelUniqueId puid;
162  _tid = puid.id;
163 
164  // Must get the variables every time this is run because _tid can change
165  if (_do_normal_smoothing &&
167  {
168  _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
169  _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
170  _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
171  }
172 
173  const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
174  std::unique_ptr<PointLocatorBase> point_locator;
175  if (_use_point_locator)
176  point_locator = _mesh.getPointLocator();
177 
178  for (const auto & node_id : range)
179  {
180  const Node & node = _mesh.nodeRef(node_id);
181 
182  // We're going to get a reference to the pointer for the pinfo for this node
183  // This will allow us to manipulate this pointer without having to go through
184  // the _penetration_info map... meaning this is the only mutex we'll have to do!
185  pinfo_mutex.lock();
188 
189  std::vector<PenetrationInfo *> p_info;
190  bool info_set(false);
191 
192  // See if we already have info about this node
193  if (info)
194  {
195  FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
196  FEBase * fe_side = _fes[_tid][info->_side->dim()];
197 
198  if (!_update_location && (info->_distance >= 0 || info->isCaptured()))
199  {
200  const Point contact_ref = info->_closest_point_ref;
201  bool contact_point_on_side(false);
202 
203  // Secondary position must be the previous contact point
204  // Use the previous reference coordinates
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;
209 
211  fe_elem,
212  fe_side,
213  _fe_type,
214  secondary_pos[0],
215  false,
217  contact_point_on_side,
218  search_succeeded);
219 
220  // Restore the original reference coordinates
221  info->_closest_point_ref = contact_ref;
222  // Just calculated as the distance of the contact point off the surface (0). Set to 0 to
223  // avoid round-off.
224  info->_distance = 0.0;
225  info_set = true;
226  }
227  else
228  {
229  Real old_tangential_distance(info->_tangential_distance);
230  bool contact_point_on_side(false);
231  bool search_succeeded = false;
232 
234  fe_elem,
235  fe_side,
236  _fe_type,
237  node,
238  false,
240  contact_point_on_side,
241  search_succeeded);
242 
243  if (contact_point_on_side)
244  {
245  if (info->_tangential_distance <= 0.0) // on the face
246  {
247  info_set = true;
248  }
249  else if (info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
250  { // off the face but within tolerance, was that way on the last step too
251  if (info->_side->dim() == 2 && info->_off_edge_nodes.size() < 2)
252  { // Closest point on face is on a node rather than an edge. Another
253  // face might be a better candidate.
254  }
255  else
256  {
257  info_set = true;
258  }
259  }
260  }
261  }
262  }
263 
264  if (!info_set)
265  {
266  const Node * closest_node = _nearest_node.nearestNode(node.id());
267 
268  std::vector<dof_id_type> located_elem_ids;
269  const std::vector<dof_id_type> * closest_elems;
270 
271  if (_use_point_locator)
272  {
273  std::set<const Elem *> candidate_elements;
274  (*point_locator)(*closest_node, candidate_elements);
275  for (const Elem * elem : candidate_elements)
276  {
277  for (auto s : elem->side_index_range())
278  if (boundary_info.has_boundary_id(elem, s, _primary_boundary))
279  {
280  located_elem_ids.push_back(elem->id());
281  break;
282  }
283  }
284  closest_elems = &located_elem_ids;
285  }
286  else
287  {
288  auto node_to_elem_pair = _node_to_elem_map.find(closest_node->id());
289  mooseAssert(node_to_elem_pair != _node_to_elem_map.end(),
290  "Missing entry in node to elem map");
291  closest_elems = &(node_to_elem_pair->second);
292  }
293 
294  for (const auto & elem_id : *closest_elems)
295  {
296  const Elem * elem = _mesh.elemPtr(elem_id);
297 
298  std::vector<PenetrationInfo *> thisElemInfo;
299 
300  std::vector<const Node *> nodesThatMustBeOnSide;
301  // If we have a disconnected mesh, we might not have *any*
302  // nodes that must be on a side we check; we'll rely on
303  // boundary info to find valid sides, then rely on comparing
304  // closest points from each to find the best.
305  //
306  // If we don't have a disconnected mesh, then for maximum
307  // backwards compatibility we're still using the older ridge
308  // and peak detection code, which depends on us ruling out
309  // sides that don't touch closest_node.
310  if (!_use_point_locator)
311  nodesThatMustBeOnSide.push_back(closest_node);
313  thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
314  }
315 
316  if (_use_point_locator)
317  {
318  Real min_distance_sq = std::numeric_limits<Real>::max();
319  Point best_point;
320  unsigned int best_i = invalid_uint;
321 
322  // Find closest point in all p_info to the node of interest
323  for (unsigned int i = 0; i < p_info.size(); ++i)
324  {
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)
328  {
329  min_distance_sq = distance_sq;
330  best_point = closest_point;
331  best_i = i;
332  }
333  }
334 
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;
342  if (dot < 0)
343  normal *= -1;
344  p_info[best_i]->_normal = normal;
345 
346  switchInfo(info, p_info[best_i]);
347  info_set = true;
348  }
349  else
350  {
351  if (p_info.size() == 1)
352  {
353  if (p_info[0]->_tangential_distance <= _tangential_tolerance)
354  {
355  switchInfo(info, p_info[0]);
356  info_set = true;
357  }
358  }
359  else if (p_info.size() > 1)
360  {
361  // Loop through all pairs of faces, and check for contact on ridge between each face pair
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)
365  {
366  Point closest_coor;
367  Real tangential_distance(0.0);
368  const Node * closest_node_on_ridge = NULL;
369  unsigned int index = 0;
370  Point closest_coor_ref;
371  bool found_ridge_contact_point = findRidgeContactPoint(closest_coor,
372  tangential_distance,
373  closest_node_on_ridge,
374  index,
375  closest_coor_ref,
376  p_info,
377  i,
378  j);
379  if (found_ridge_contact_point)
380  {
381  RidgeData rpd;
382  rpd._closest_coor = closest_coor;
383  rpd._tangential_distance = tangential_distance;
384  rpd._closest_node = closest_node_on_ridge;
385  rpd._index = index;
386  rpd._closest_coor_ref = closest_coor_ref;
387  ridgeDataVec.push_back(rpd);
388  }
389  }
390 
391  if (ridgeDataVec.size() > 0) // Either find the ridge pair that is the best or find a peak
392  {
393  // Group together ridges for which we are off the edge of a common node.
394  // Those are peaks.
395  std::vector<RidgeSetData> ridgeSetDataVec;
396  for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
397  {
398  bool foundSetWithMatchingNode = false;
399  for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
400  {
401  if (ridgeDataVec[i]._closest_node != NULL &&
402  ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
403  {
404  foundSetWithMatchingNode = true;
405  ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
406  break;
407  }
408  }
409  if (!foundSetWithMatchingNode)
410  {
411  RidgeSetData rsd;
413  rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
414  rsd._closest_node = ridgeDataVec[i]._closest_node;
415  ridgeSetDataVec.push_back(rsd);
416  }
417  }
418  // Compute distance to each set of ridges
419  for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
420  {
421  if (ridgeSetDataVec[i]._closest_node !=
422  NULL) // Either a peak or off the edge of single ridge
423  {
424  if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
425  {
426  if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
427  _tangential_tolerance) // off within tolerance
428  {
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();
433  }
434  }
435  else // several ridges join at common node to make peak. The common node is the
436  // contact point
437  {
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();
441  }
442  }
443  else // on a single ridge
444  {
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();
449  }
450  }
451  // Find the set of ridges closest to us.
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)
456  {
457  if (ridgeSetDataVec[i]._distance < closest_distance)
458  {
459  closest_ridge_set_index = i;
460  closest_distance = ridgeSetDataVec[i]._distance;
461  closest_point = ridgeSetDataVec[i]._closest_coor;
462  }
463  }
464 
465  if (closest_distance <
466  std::numeric_limits<Real>::max()) // contact point is on the closest ridge set
467  {
468  // find the face in the ridge set with the smallest index, assign that one to the
469  // interaction
470  unsigned int face_index(std::numeric_limits<unsigned int>::max());
471  for (unsigned int i = 0;
472  i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
473  ++i)
474  {
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;
477  }
478 
479  mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
480  "face_index invalid");
481 
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;
485  // Calculate the normal as the vector from the ridge to the point only if we're not
486  // doing normal
487  // smoothing. Normal smoothing will average out the normals on its own.
489  {
490  Point normal(closest_point - node);
491  const Real len(normal.norm());
492  if (len > 0)
493  {
494  normal /= len;
495  }
496  const Real dot(normal * p_info[face_index]->_normal);
497  if (dot < 0)
498  {
499  normal *= -1;
500  }
501  p_info[face_index]->_normal = normal;
502  }
503  p_info[face_index]->_tangential_distance = 0.0;
504 
505  Point closest_point_ref;
506  if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
507  1) // contact with a single ridge rather than a peak
508  {
509  p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
510  ._ridge_data_vec[0]
511  ._tangential_distance;
512  p_info[face_index]->_closest_point_ref =
513  ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
514  }
515  else
516  { // peak
517  const Node * closest_node_on_face;
518  bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
519  closest_node_on_face,
520  p_info[face_index]->_side);
521  if (restricted)
522  {
523  if (closest_node_on_face !=
524  ridgeSetDataVec[closest_ridge_set_index]._closest_node)
525  {
526  mooseError("Closest node when restricting point to face != closest node from "
527  "RidgeSetData");
528  }
529  }
530  }
531 
532  FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
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();
538  p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
539  p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
540  p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
541 
542  switchInfo(info, p_info[face_index]);
543  info_set = true;
544  }
545  else
546  { // todo:remove invalid ridge cases so they don't mess up individual face
547  // competition????
548  }
549  }
550 
551  if (!info_set) // contact wasn't on a ridge -- compete individual interactions
552  {
553  unsigned int best(0), i(1);
554  do
555  {
556  CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
557  if (CIResult == FIRST_WINS)
558  {
559  i++;
560  }
561  else if (CIResult == SECOND_WINS)
562  {
563  best = i;
564  i++;
565  }
566  else if (CIResult == NEITHER_WINS)
567  {
568  best = i + 1;
569  i += 2;
570  }
571  } while (i < p_info.size() && best < p_info.size());
572  if (best < p_info.size())
573  {
574  // Ensure final info is within the tangential tolerance
575  if (p_info[best]->_tangential_distance <= _tangential_tolerance)
576  {
577  switchInfo(info, p_info[best]);
578  info_set = true;
579  }
580  }
581  }
582  }
583  }
584  }
585 
586  if (!info_set)
587  {
588  // If penetration is not detected within the saved patch, it is possible
589  // that the secondary node has moved outside the saved patch. So, the patch
590  // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
591  // and the penetration detection has to be re-run on the updated patch.
592 
593  _recheck_secondary_nodes.push_back(node_id);
594 
595  delete info;
596  info = NULL;
597  }
598  else
599  {
600  smoothNormal(info, p_info, node);
601  FEBase * fe = _fes[_tid][info->_side->dim()];
602  computeSlip(*fe, *info);
603  }
604 
605  for (unsigned int j = 0; j < p_info.size(); ++j)
606  {
607  if (p_info[j])
608  {
609  delete p_info[j];
610  p_info[j] = NULL;
611  }
612  }
613  }
614 }
615 
616 void
618 {
620  other._recheck_secondary_nodes.begin(),
621  other._recheck_secondary_nodes.end());
622 }
623 
624 void
626 {
627  mooseAssert(infoNew != NULL, "infoNew object is null");
628  if (info)
629  {
630  infoNew->_starting_elem = info->_starting_elem;
631  infoNew->_starting_side_num = info->_starting_side_num;
632  infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
633  infoNew->_incremental_slip = info->_incremental_slip;
634  infoNew->_accumulated_slip = info->_accumulated_slip;
635  infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
636  infoNew->_frictional_energy = info->_frictional_energy;
637  infoNew->_frictional_energy_old = info->_frictional_energy_old;
638  infoNew->_contact_force = info->_contact_force;
639  infoNew->_contact_force_old = info->_contact_force_old;
640  infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
641  infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
642  infoNew->_locked_this_step = info->_locked_this_step;
643  infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
644  infoNew->_mech_status = info->_mech_status;
645  infoNew->_mech_status_old = info->_mech_status_old;
646  }
647  else
648  {
649  infoNew->_starting_elem = infoNew->_elem;
650  infoNew->_starting_side_num = infoNew->_side_num;
652  }
653  delete info;
654  info = infoNew;
655  infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
656 }
657 
660 {
661 
663 
665  pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
666  result = NEITHER_WINS;
667 
668  else if (pi1->_tangential_distance == 0.0 &&
669  pi2->_tangential_distance > 0.0) // on face 1, off face 2
670  result = FIRST_WINS;
671 
672  else if (pi2->_tangential_distance == 0.0 &&
673  pi1->_tangential_distance > 0.0) // on face 2, off face 1
674  result = SECOND_WINS;
675 
676  else if (pi1->_tangential_distance <= _tangential_tolerance &&
677  pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
678  result = FIRST_WINS;
679 
680  else if (pi2->_tangential_distance <= _tangential_tolerance &&
681  pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
682  result = SECOND_WINS;
683 
684  else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
685  result = competeInteractionsBothOnFace(pi1, pi2);
686 
687  else if (pi1->_tangential_distance <= _tangential_tolerance &&
688  pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
689  {
691  if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
692  {
693  // We already checked for ridges, and it got rejected, so neither face must be valid
694  result = NEITHER_WINS;
695  // mooseError("Erroneously encountered ridge case");
696  }
697  else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
698  // the off-side face
699  {
700  if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
701  mooseError("Invalid off_edge_nodes counts");
702 
703  else if (pi1->_off_edge_nodes.size() == 2)
704  result = FIRST_WINS;
705 
706  else if (pi2->_off_edge_nodes.size() == 2)
707  result = SECOND_WINS;
708 
709  else
710  mooseError("Invalid off_edge_nodes counts");
711  }
712  else // The node projects to both faces within tangential tolerance.
713  result = competeInteractionsBothOnFace(pi1, pi2);
714  }
715 
716  return result;
717 }
718 
721 {
723 
724  if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
725  result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
726 
727  else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
728  result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
729 
730  // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
731  // movement. If there is some way to smooth the transition, we should do it.
733  result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
734 
736  result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
737 
738  else // Equal within tolerance. Favor the one with a smaller element id (for repeatibility)
739  {
740  if (pi1->_elem->id() < pi2->_elem->id())
741  result = FIRST_WINS;
742 
743  else
744  result = SECOND_WINS;
745  }
746 
747  return result;
748 }
749 
752 {
753  CommonEdgeResult common_edge(NO_COMMON);
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();
757 
758  if (dim1 == 1)
759  {
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])
765  common_edge = COMMON_EDGE;
766  }
767  else
768  {
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)
773  {
774  if (off_edge_nodes2.size() == 1)
775  {
776  if (off_edge_nodes1[0] == off_edge_nodes2[0])
777  common_edge = COMMON_NODE;
778  }
779  else if (off_edge_nodes2.size() == 2)
780  {
781  if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
782  common_edge = EDGE_AND_COMMON_NODE;
783  }
784  }
785  else if (off_edge_nodes1.size() == 2)
786  {
787  if (off_edge_nodes2.size() == 1)
788  {
789  if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
790  common_edge = EDGE_AND_COMMON_NODE;
791  }
792  else if (off_edge_nodes2.size() == 2)
793  {
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]))
797  common_edge = COMMON_EDGE;
798  }
799  }
800  }
801  return common_edge;
802 }
803 
804 bool
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)
813 {
814  tangential_distance = 0.0;
815  closest_node = NULL;
816  PenetrationInfo * pi1 = p_info[index1];
817  PenetrationInfo * pi2 = p_info[index2];
818  const unsigned sidedim(pi1->_side->dim());
819  mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
820 
821  // Nodes on faces for the two interactions
822  std::vector<const Node *> side1_nodes;
823  getSideCornerNodes(pi1->_side, side1_nodes);
824  std::vector<const Node *> side2_nodes;
825  getSideCornerNodes(pi2->_side, side2_nodes);
826 
827  std::sort(side1_nodes.begin(), side1_nodes.end());
828  std::sort(side2_nodes.begin(), side2_nodes.end());
829 
830  // Find nodes shared by the two faces
831  std::vector<const Node *> common_nodes;
832  std::set_intersection(side1_nodes.begin(),
833  side1_nodes.end(),
834  side2_nodes.begin(),
835  side2_nodes.end(),
836  std::inserter(common_nodes, common_nodes.end()));
837 
838  if (common_nodes.size() != sidedim)
839  return false;
840 
841  bool found_point1, found_point2;
842  Point closest_coor_ref1(pi1->_closest_point_ref);
843  const Node * closest_node1;
844  found_point1 = restrictPointToSpecifiedEdgeOfFace(
845  closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
846 
847  Point closest_coor_ref2(pi2->_closest_point_ref);
848  const Node * closest_node2;
849  found_point2 = restrictPointToSpecifiedEdgeOfFace(
850  closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
851 
852  if (!found_point1 || !found_point2)
853  return false;
854 
855  // if (sidedim == 2)
856  // {
857  // TODO:
858  // We have the parametric coordinates of the closest intersection point for both faces.
859  // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
860  // Find that point by taking dot products of vector from contact point to secondary node point
861  // with face normal vectors to see which face we're closer to.
862  // }
863 
864  FEBase * fe = NULL;
865  std::vector<Point> points(1);
866 
867  // We have to pick one of the two faces to own the contact point. It doesn't really matter
868  // which one we pick. For repeatibility, pick the face with the lowest index.
869  if (index1 < index2)
870  {
871  fe = _fes[_tid][pi1->_side->dim()];
872  contact_point_ref = closest_coor_ref1;
873  points[0] = closest_coor_ref1;
874  fe->reinit(pi1->_side, &points);
875  index = index1;
876  }
877  else
878  {
879  fe = _fes[_tid][pi2->_side->dim()];
880  contact_point_ref = closest_coor_ref2;
881  points[0] = closest_coor_ref2;
882  fe->reinit(pi2->_side, &points);
883  index = index2;
884  }
885 
886  contact_point = fe->get_xyz()[0];
887 
888  if (sidedim == 2)
889  {
890  if (closest_node1) // point is off the ridge between the two elements
891  {
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;
895 
896  RealGradient off_face = *closest_node1 - contact_point;
897  tangential_distance = off_face.norm();
898  }
899  }
900 
901  return true;
902 }
903 
904 void
905 PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
906 {
907  const ElemType t(side->type());
908  corner_nodes.clear();
909 
910  corner_nodes.push_back(side->node_ptr(0));
911  corner_nodes.push_back(side->node_ptr(1));
912  switch (t)
913  {
914  case EDGE2:
915  case EDGE3:
916  case EDGE4:
917  {
918  break;
919  }
920 
921  case TRI3:
922  case TRI6:
923  case TRI7:
924  {
925  corner_nodes.push_back(side->node_ptr(2));
926  break;
927  }
928 
929  case QUAD4:
930  case QUAD8:
931  case QUAD9:
932  {
933  corner_nodes.push_back(side->node_ptr(2));
934  corner_nodes.push_back(side->node_ptr(3));
935  break;
936  }
937 
938  default:
939  {
940  mooseError("Unsupported face type: ", t);
941  break;
942  }
943  }
944 }
945 
946 bool
948  const Node *& closest_node,
949  const Elem * side,
950  const std::vector<const Node *> & edge_nodes)
951 {
952  const ElemType t = side->type();
953  Real & xi = p(0);
954  Real & eta = p(1);
955  closest_node = NULL;
956 
957  std::vector<unsigned int> local_node_indices;
958  for (const auto & edge_node : edge_nodes)
959  {
960  unsigned int local_index = side->get_node_index(edge_node);
961  if (local_index == libMesh::invalid_uint)
962  mooseError("Side does not contain node");
963  local_node_indices.push_back(local_index);
964  }
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());
968 
969  bool off_of_this_edge = false;
970 
971  switch (t)
972  {
973  case EDGE2:
974  case EDGE3:
975  case EDGE4:
976  {
977  if (local_node_indices[0] == 0)
978  {
979  if (xi <= -1.0)
980  {
981  xi = -1.0;
982  off_of_this_edge = true;
983  closest_node = side->node_ptr(0);
984  }
985  }
986  else if (local_node_indices[0] == 1)
987  {
988  if (xi >= 1.0)
989  {
990  xi = 1.0;
991  off_of_this_edge = true;
992  closest_node = side->node_ptr(1);
993  }
994  }
995  else
996  {
997  mooseError("Invalid local node indices");
998  }
999  break;
1000  }
1001 
1002  case TRI3:
1003  case TRI6:
1004  case TRI7:
1005  {
1006  if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1007  {
1008  if (eta <= 0.0)
1009  {
1010  eta = 0.0;
1011  off_of_this_edge = true;
1012  if (xi < 0.0)
1013  closest_node = side->node_ptr(0);
1014  else if (xi > 1.0)
1015  closest_node = side->node_ptr(1);
1016  }
1017  }
1018  else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1019  {
1020  if ((xi + eta) > 1.0)
1021  {
1022  Real delta = (xi + eta - 1.0) / 2.0;
1023  xi -= delta;
1024  eta -= delta;
1025  off_of_this_edge = true;
1026  if (xi > 1.0)
1027  closest_node = side->node_ptr(1);
1028  else if (xi < 0.0)
1029  closest_node = side->node_ptr(2);
1030  }
1031  }
1032  else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1033  {
1034  if (xi <= 0.0)
1035  {
1036  xi = 0.0;
1037  off_of_this_edge = true;
1038  if (eta > 1.0)
1039  closest_node = side->node_ptr(2);
1040  else if (eta < 0.0)
1041  closest_node = side->node_ptr(0);
1042  }
1043  }
1044  else
1045  {
1046  mooseError("Invalid local node indices");
1047  }
1048 
1049  break;
1050  }
1051 
1052  case QUAD4:
1053  case QUAD8:
1054  case QUAD9:
1055  {
1056  if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1057  {
1058  if (eta <= -1.0)
1059  {
1060  eta = -1.0;
1061  off_of_this_edge = true;
1062  if (xi < -1.0)
1063  closest_node = side->node_ptr(0);
1064  else if (xi > 1.0)
1065  closest_node = side->node_ptr(1);
1066  }
1067  }
1068  else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1069  {
1070  if (xi >= 1.0)
1071  {
1072  xi = 1.0;
1073  off_of_this_edge = true;
1074  if (eta < -1.0)
1075  closest_node = side->node_ptr(1);
1076  else if (eta > 1.0)
1077  closest_node = side->node_ptr(2);
1078  }
1079  }
1080  else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1081  {
1082  if (eta >= 1.0)
1083  {
1084  eta = 1.0;
1085  off_of_this_edge = true;
1086  if (xi < -1.0)
1087  closest_node = side->node_ptr(3);
1088  else if (xi > 1.0)
1089  closest_node = side->node_ptr(2);
1090  }
1091  }
1092  else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1093  {
1094  if (xi <= -1.0)
1095  {
1096  xi = -1.0;
1097  off_of_this_edge = true;
1098  if (eta < -1.0)
1099  closest_node = side->node_ptr(0);
1100  else if (eta > 1.0)
1101  closest_node = side->node_ptr(3);
1102  }
1103  }
1104  else
1105  {
1106  mooseError("Invalid local node indices");
1107  }
1108  break;
1109  }
1110 
1111  default:
1112  {
1113  mooseError("Unsupported face type: ", t);
1114  break;
1115  }
1116  }
1117  return off_of_this_edge;
1118 }
1119 
1120 bool
1121 PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
1122 {
1123  const ElemType t(side->type());
1124  Real & xi = p(0);
1125  Real & eta = p(1);
1126  closest_node = NULL;
1127 
1128  bool off_of_this_face(false);
1129 
1130  switch (t)
1131  {
1132  case EDGE2:
1133  case EDGE3:
1134  case EDGE4:
1135  {
1136  if (xi < -1.0)
1137  {
1138  xi = -1.0;
1139  off_of_this_face = true;
1140  closest_node = side->node_ptr(0);
1141  }
1142  else if (xi > 1.0)
1143  {
1144  xi = 1.0;
1145  off_of_this_face = true;
1146  closest_node = side->node_ptr(1);
1147  }
1148  break;
1149  }
1150 
1151  case TRI3:
1152  case TRI6:
1153  case TRI7:
1154  {
1155  if (eta < 0.0)
1156  {
1157  eta = 0.0;
1158  off_of_this_face = true;
1159  if (xi < 0.5)
1160  {
1161  closest_node = side->node_ptr(0);
1162  if (xi < 0.0)
1163  xi = 0.0;
1164  }
1165  else
1166  {
1167  closest_node = side->node_ptr(1);
1168  if (xi > 1.0)
1169  xi = 1.0;
1170  }
1171  }
1172  else if ((xi + eta) > 1.0)
1173  {
1174  Real delta = (xi + eta - 1.0) / 2.0;
1175  xi -= delta;
1176  eta -= delta;
1177  off_of_this_face = true;
1178  if (xi > 0.5)
1179  {
1180  closest_node = side->node_ptr(1);
1181  if (xi > 1.0)
1182  {
1183  xi = 1.0;
1184  eta = 0.0;
1185  }
1186  }
1187  else
1188  {
1189  closest_node = side->node_ptr(2);
1190  if (xi < 0.0)
1191  {
1192  xi = 0.0;
1193  eta = 1.0;
1194  }
1195  }
1196  }
1197  else if (xi < 0.0)
1198  {
1199  xi = 0.0;
1200  off_of_this_face = true;
1201  if (eta > 0.5)
1202  {
1203  closest_node = side->node_ptr(2);
1204  if (eta > 1.0)
1205  eta = 1.0;
1206  }
1207  else
1208  {
1209  closest_node = side->node_ptr(0);
1210  if (eta < 0.0)
1211  eta = 0.0;
1212  }
1213  }
1214  break;
1215  }
1216 
1217  case QUAD4:
1218  case QUAD8:
1219  case QUAD9:
1220  {
1221  if (eta < -1.0)
1222  {
1223  eta = -1.0;
1224  off_of_this_face = true;
1225  if (xi < 0.0)
1226  {
1227  closest_node = side->node_ptr(0);
1228  if (xi < -1.0)
1229  xi = -1.0;
1230  }
1231  else
1232  {
1233  closest_node = side->node_ptr(1);
1234  if (xi > 1.0)
1235  xi = 1.0;
1236  }
1237  }
1238  else if (xi > 1.0)
1239  {
1240  xi = 1.0;
1241  off_of_this_face = true;
1242  if (eta < 0.0)
1243  {
1244  closest_node = side->node_ptr(1);
1245  if (eta < -1.0)
1246  eta = -1.0;
1247  }
1248  else
1249  {
1250  closest_node = side->node_ptr(2);
1251  if (eta > 1.0)
1252  eta = 1.0;
1253  }
1254  }
1255  else if (eta > 1.0)
1256  {
1257  eta = 1.0;
1258  off_of_this_face = true;
1259  if (xi < 0.0)
1260  {
1261  closest_node = side->node_ptr(3);
1262  if (xi < -1.0)
1263  xi = -1.0;
1264  }
1265  else
1266  {
1267  closest_node = side->node_ptr(2);
1268  if (xi > 1.0)
1269  xi = 1.0;
1270  }
1271  }
1272  else if (xi < -1.0)
1273  {
1274  xi = -1.0;
1275  off_of_this_face = true;
1276  if (eta < 0.0)
1277  {
1278  closest_node = side->node_ptr(0);
1279  if (eta < -1.0)
1280  eta = -1.0;
1281  }
1282  else
1283  {
1284  closest_node = side->node_ptr(3);
1285  if (eta > 1.0)
1286  eta = 1.0;
1287  }
1288  }
1289  break;
1290  }
1291 
1292  default:
1293  {
1294  mooseError("Unsupported face type: ", t);
1295  break;
1296  }
1297  }
1298  return off_of_this_face;
1299 }
1300 
1301 bool
1303  const Elem * side,
1304  FEBase * fe,
1305  const Point * secondary_point,
1306  const Real tangential_tolerance)
1307 {
1308  unsigned int dim = primary_elem->dim();
1309 
1310  const std::vector<Point> & phys_point = fe->get_xyz();
1311 
1312  const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1313  const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1314 
1315  Point ref_point;
1316 
1317  std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1318 
1319  fe->reinit(side, &points);
1320 
1321  RealGradient d = *secondary_point - phys_point[0];
1322 
1323  const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1324  Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1325 
1326  RealVectorValue normal;
1327  if (dim - 1 == 2)
1328  {
1329  normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1330  }
1331  else if (dim - 1 == 1)
1332  {
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];
1336 
1337  Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1338  out_of_plane_normal /= out_of_plane_normal.norm();
1339 
1340  normal = dxyz_dxi[0].cross(out_of_plane_normal);
1341  }
1342  else
1343  {
1344  return true;
1345  }
1346  normal /= normal.norm();
1347 
1348  const Real dot(d * normal);
1349 
1350  const RealGradient normcomp = dot * normal;
1351  const RealGradient tangcomp = d - normcomp;
1352 
1353  const Real tangdist = tangcomp.norm();
1354 
1355  // Increase the size of the zone that we consider if the vector from the face
1356  // to the node has a larger normal component
1357  const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1358 
1359  bool isReasonableCandidate = true;
1360  if (tangdist > faceExpansionFactor * max_face_length)
1361  {
1362  isReasonableCandidate = false;
1363  }
1364  return isReasonableCandidate;
1365 }
1366 
1367 void
1369 {
1370  // Slip is current projected position of secondary node minus
1371  // original projected position of secondary node
1372  std::vector<Point> points(1);
1373  points[0] = info._starting_closest_point_ref;
1374  const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
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())
1379  {
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();
1383  }
1384 }
1385 
1386 void
1388  std::vector<PenetrationInfo *> & p_info,
1389  const Node & node)
1390 {
1392  {
1394  {
1395  // If we are within the smoothing distance of any edges or corners, find the
1396  // corner nodes for those edges/corners, and weights from distance to edge/corner
1397  std::vector<Real> edge_face_weights;
1398  std::vector<PenetrationInfo *> edge_face_info;
1399 
1400  getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
1401 
1402  mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1403  "edge_face_info.size() != edge_face_weights.size()");
1404 
1405  if (edge_face_info.size() > 0)
1406  {
1407  // Smooth the normal using the weighting functions for all participating faces.
1408  RealVectorValue new_normal;
1409  Real this_face_weight = 1.0;
1410 
1411  for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1412  {
1413  PenetrationInfo * npi = edge_face_info[efwi];
1414  if (npi)
1415  new_normal += npi->_normal * edge_face_weights[efwi];
1416 
1417  this_face_weight -= edge_face_weights[efwi];
1418  }
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;
1422 
1423  const Real len = new_normal.norm();
1424  if (len > 0)
1425  new_normal /= len;
1426 
1427  info->_normal = new_normal;
1428  }
1429  }
1431  {
1432  // params.addParam<VariableName>("var_name","description");
1433  // getParam<VariableName>("var_name")
1434  info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1435  info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1436  info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1437  const Real len(info->_normal.norm());
1438  if (len > 0)
1439  info->_normal /= len;
1440  }
1441  }
1442 }
1443 
1444 void
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)
1450 {
1451  const Elem * side = info->_side;
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());
1455 
1456  std::vector<std::vector<const Node *>> edge_nodes;
1457 
1458  // Get the pairs of nodes along every edge that we are close enough to smooth with
1459  getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1460  std::vector<Elem *> edge_neighbor_elems;
1461  edge_face_info.resize(edge_nodes.size(), NULL);
1462 
1463  std::vector<unsigned int> edges_without_neighbors;
1464 
1465  for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1466  {
1467  // Sort all sets of edge nodes (needed for comparing edges)
1468  std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1469 
1470  std::vector<PenetrationInfo *> face_info_comm_edge;
1472  &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1473 
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");
1478  else
1479  edge_face_info[i] = face_info_comm_edge[0];
1480  }
1481 
1482  // Remove edges without neighbors from the vector, starting from end
1483  std::vector<unsigned int>::reverse_iterator rit;
1484  for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1485  {
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);
1490  }
1491 
1492  // Handle corner case
1493  if (edge_nodes.size() > 1)
1494  {
1495  if (edge_nodes.size() != 2)
1496  mooseError("Invalid number of smoothing edges");
1497 
1498  // find common node
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()));
1505 
1506  if (common_nodes.size() != 1)
1507  mooseError("Invalid number of common nodes");
1508 
1509  for (const auto & pinfo : edge_face_info)
1510  elems_to_exclude.insert(pinfo->_elem->id());
1511 
1512  std::vector<PenetrationInfo *> face_info_comm_edge;
1514  &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1515 
1516  unsigned int num_corner_neighbors = face_info_comm_edge.size();
1517 
1518  if (num_corner_neighbors > 0)
1519  {
1520  Real fw0 = edge_face_weights[0];
1521  Real fw1 = edge_face_weights[1];
1522 
1523  // Corner weight is product of edge weights. Spread out over multiple neighbors.
1524  Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1525 
1526  // Adjust original edge weights
1527  edge_face_weights[0] *= (1.0 - fw1);
1528  edge_face_weights[1] *= (1.0 - fw0);
1529 
1530  for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1531  {
1532  edge_face_weights.push_back(fw_corner);
1533  edge_face_info.push_back(face_info_comm_edge[i]);
1534  }
1535  }
1536  }
1537 }
1538 
1539 void
1541  const Point & p,
1542  const Elem * side,
1543  std::vector<std::vector<const Node *>> & edge_nodes,
1544  std::vector<Real> & edge_face_weights)
1545 {
1546  const ElemType t(side->type());
1547  const Real & xi = p(0);
1548  const Real & eta = p(1);
1549 
1550  Real smooth_limit = 1.0 - _normal_smoothing_distance;
1551 
1552  switch (t)
1553  {
1554  case EDGE2:
1555  case EDGE3:
1556  case EDGE4:
1557  {
1558  if (xi < -smooth_limit)
1559  {
1560  std::vector<const Node *> en;
1561  en.push_back(side->node_ptr(0));
1562  edge_nodes.push_back(en);
1563  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1564  if (fw > 0.5)
1565  fw = 0.5;
1566  edge_face_weights.push_back(fw);
1567  }
1568  else if (xi > smooth_limit)
1569  {
1570  std::vector<const Node *> en;
1571  en.push_back(side->node_ptr(1));
1572  edge_nodes.push_back(en);
1573  Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1574  if (fw > 0.5)
1575  fw = 0.5;
1576  edge_face_weights.push_back(fw);
1577  }
1578  break;
1579  }
1580 
1581  case TRI3:
1582  case TRI6:
1583  case TRI7:
1584  {
1585  if (eta < -smooth_limit)
1586  {
1587  std::vector<const Node *> en;
1588  en.push_back(side->node_ptr(0));
1589  en.push_back(side->node_ptr(1));
1590  edge_nodes.push_back(en);
1591  Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1592  if (fw > 0.5)
1593  fw = 0.5;
1594  edge_face_weights.push_back(fw);
1595  }
1596  if ((xi + eta) > smooth_limit)
1597  {
1598  std::vector<const Node *> en;
1599  en.push_back(side->node_ptr(1));
1600  en.push_back(side->node_ptr(2));
1601  edge_nodes.push_back(en);
1602  Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
1603  if (fw > 0.5)
1604  fw = 0.5;
1605  edge_face_weights.push_back(fw);
1606  }
1607  if (xi < -smooth_limit)
1608  {
1609  std::vector<const Node *> en;
1610  en.push_back(side->node_ptr(2));
1611  en.push_back(side->node_ptr(0));
1612  edge_nodes.push_back(en);
1613  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1614  if (fw > 0.5)
1615  fw = 0.5;
1616  edge_face_weights.push_back(fw);
1617  }
1618  break;
1619  }
1620 
1621  case QUAD4:
1622  case QUAD8:
1623  case QUAD9:
1624  {
1625  if (eta < -smooth_limit)
1626  {
1627  std::vector<const Node *> en;
1628  en.push_back(side->node_ptr(0));
1629  en.push_back(side->node_ptr(1));
1630  edge_nodes.push_back(en);
1631  Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1632  if (fw > 0.5)
1633  fw = 0.5;
1634  edge_face_weights.push_back(fw);
1635  }
1636  if (xi > smooth_limit)
1637  {
1638  std::vector<const Node *> en;
1639  en.push_back(side->node_ptr(1));
1640  en.push_back(side->node_ptr(2));
1641  edge_nodes.push_back(en);
1642  Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1643  if (fw > 0.5)
1644  fw = 0.5;
1645  edge_face_weights.push_back(fw);
1646  }
1647  if (eta > smooth_limit)
1648  {
1649  std::vector<const Node *> en;
1650  en.push_back(side->node_ptr(2));
1651  en.push_back(side->node_ptr(3));
1652  edge_nodes.push_back(en);
1653  Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1654  if (fw > 0.5)
1655  fw = 0.5;
1656  edge_face_weights.push_back(fw);
1657  }
1658  if (xi < -smooth_limit)
1659  {
1660  std::vector<const Node *> en;
1661  en.push_back(side->node_ptr(3));
1662  en.push_back(side->node_ptr(0));
1663  edge_nodes.push_back(en);
1664  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1665  if (fw > 0.5)
1666  fw = 0.5;
1667  edge_face_weights.push_back(fw);
1668  }
1669  break;
1670  }
1671 
1672  default:
1673  {
1674  mooseError("Unsupported face type: ", t);
1675  break;
1676  }
1677  }
1678 }
1679 
1680 void
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)
1687 {
1688  // elems connected to a node on this edge, find one that has the same corners as this, and is not
1689  // the current elem
1690  auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
1691  // nodes
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;
1694 
1695  std::vector<const Elem *> elems_connected_to_edge;
1696 
1697  for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1698  {
1699  if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1700  continue;
1701  const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1702 
1703  std::vector<const Node *> nodevec;
1704  for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1705  if (elem->is_vertex(ni))
1706  nodevec.push_back(elem->node_ptr(ni));
1707 
1708  std::vector<const Node *> common_nodes;
1709  std::sort(nodevec.begin(), nodevec.end());
1710  std::set_intersection(edge_nodes.begin(),
1711  edge_nodes.end(),
1712  nodevec.begin(),
1713  nodevec.end(),
1714  std::inserter(common_nodes, common_nodes.end()));
1715 
1716  if (common_nodes.size() == edge_nodes.size())
1717  elems_connected_to_edge.push_back(elem);
1718  }
1719 
1720  if (elems_connected_to_edge.size() > 0)
1721  {
1722 
1723  // There are potentially multiple elements that share a common edge
1724  // 2D:
1725  // There can only be one element on the same surface
1726  // 3D:
1727  // If there are two edge nodes, there can only be one element on the same surface
1728  // If there is only one edge node (a corner), there could be multiple elements on the same
1729  // surface
1730  bool allowMultipleNeighbors = false;
1731 
1732  if (elems_connected_to_edge[0]->dim() == 3)
1733  {
1734  if (edge_nodes.size() == 1)
1735  {
1736  allowMultipleNeighbors = true;
1737  }
1738  }
1739 
1740  for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1741  {
1742  std::vector<PenetrationInfo *> thisElemInfo;
1743  getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1744  if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1745  {
1746  if (thisElemInfo.size() > 1)
1747  mooseError(
1748  "Found multiple neighbors to current edge/face on surface when only one is allowed");
1749  face_info_comm_edge.push_back(thisElemInfo[0]);
1750  break;
1751  }
1752 
1754  thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1755  if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1756  {
1757  if (thisElemInfo.size() > 1)
1758  mooseError(
1759  "Found multiple neighbors to current edge/face on surface when only one is allowed");
1760  face_info_comm_edge.push_back(thisElemInfo[0]);
1761  break;
1762  }
1763 
1764  for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1765  face_info_comm_edge.push_back(thisElemInfo[j]);
1766  }
1767  }
1768 }
1769 
1770 void
1771 PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1772  std::vector<PenetrationInfo *> & p_info,
1773  const Elem * elem)
1774 {
1775  for (const auto & pi : p_info)
1776  {
1777  if (!pi)
1778  continue;
1779 
1780  if (pi->_elem == elem)
1781  thisElemInfo.push_back(pi);
1782  }
1783 }
1784 
1785 void
1786 PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1787  std::vector<PenetrationInfo *> & p_info,
1788  const Node * secondary_node,
1789  const Elem * elem,
1790  const std::vector<const Node *> & nodes_that_must_be_on_side,
1791  const bool check_whether_reasonable)
1792 {
1793  const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
1794 
1795  for (auto s : elem->side_index_range())
1796  {
1797  if (!boundary_info.has_boundary_id(elem, s, _primary_boundary))
1798  continue;
1799 
1800  // Don't create info for this side if one already exists
1801  bool already_have_info_this_side = false;
1802  for (const auto & pi : thisElemInfo)
1803  if (pi->_side_num == s)
1804  {
1805  already_have_info_this_side = true;
1806  break;
1807  }
1808 
1809  if (already_have_info_this_side)
1810  break;
1811 
1812  const Elem * side = elem->build_side_ptr(s).release();
1813 
1814  // Only continue with creating info for this side if the side contains
1815  // all of the nodes in nodes_that_must_be_on_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));
1819 
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(),
1824  nodevec.begin(),
1825  nodevec.end(),
1826  std::inserter(common_nodes, common_nodes.end()));
1827  if (common_nodes.size() != nodes_that_must_be_on_side.size())
1828  {
1829  delete side;
1830  break;
1831  }
1832 
1833  FEBase * fe_elem = _fes[_tid][elem->dim()];
1834  FEBase * fe_side = _fes[_tid][side->dim()];
1835 
1836  // Optionally check to see whether face is reasonable candidate based on an
1837  // estimate of how closely it is likely to project to the face
1838  if (check_whether_reasonable)
1839  if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1840  {
1841  delete side;
1842  break;
1843  }
1844 
1845  Point contact_phys;
1846  Point contact_ref;
1847  Point contact_on_face_ref;
1848  Real distance = 0.;
1849  Real tangential_distance = 0.;
1850  RealGradient normal;
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;
1858 
1859  std::unique_ptr<PenetrationInfo> pen_info =
1860  std::make_unique<PenetrationInfo>(elem,
1861  side,
1862  s,
1863  normal,
1864  distance,
1865  tangential_distance,
1866  contact_phys,
1867  contact_ref,
1868  contact_on_face_ref,
1869  off_edge_nodes,
1870  side_phi,
1871  side_grad_phi,
1872  dxyzdxi,
1873  dxyzdeta,
1874  d2xyzdxideta);
1875 
1876  bool search_succeeded = false;
1877  Moose::findContactPoint(*pen_info,
1878  fe_elem,
1879  fe_side,
1880  _fe_type,
1881  *secondary_node,
1882  true,
1884  contact_point_on_side,
1885  search_succeeded);
1886 
1887  // Do not add contact info from failed searches
1888  if (search_succeeded)
1889  {
1890  thisElemInfo.push_back(pen_info.get());
1891  p_info.push_back(pen_info.release());
1892  }
1893  }
1894 }
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)
Definition: EigenADReal.h:42
ElemType
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
RealVectorValue _normal
virtual_for_inffe const std::vector< RealGradient > & get_dxyzdeta() const
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:3153
MPI_Info info
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...
Definition: MooseError.h:323
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
const Elem * _starting_elem
MeshBase & mesh
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:159
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
Definition: MooseMesh.C:849
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
SubProblem & _subproblem
const std::vector< std::vector< OutputGradient > > & get_dphi() const
const MooseMesh & _mesh
auto norm_sq() const -> decltype(std::norm(Real()))
dof_id_type id() const
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3488
virtual unsigned int n_nodes() const=0
const Elem * _elem
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 *const * get_nodes() const
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:92
void computeSlip(libMesh::FEBase &fe, PenetrationInfo &info)
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...
unsigned int _side_num
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.
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)
const Elem * _side
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
RealVectorValue _contact_force
MECH_STATUS_ENUM _mech_status
void join(const PenetrationThread &other)
std::vector< const Node * > _off_edge_nodes
Point _starting_closest_point_ref
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.
Definition: SubProblem.h:78
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
auto norm_sq(const T &a) -> decltype(std::norm(a))
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)
TRI7
MECH_STATUS_ENUM _mech_status_old
virtual std::unique_ptr< libMesh::PointLocatorBase > getPointLocator() const
Proxy function to get a (sub)PointLocator from either the underlying libMesh mesh (default)...
Definition: MooseMesh.C:3773
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...
Definition: MooseUtils.h:629
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
uint8_t dof_id_type
const Real pi
const std::vector< std::vector< OutputShape > > & get_phi() const
std::vector< RidgeData > _ridge_data_vec