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 
276  if (candidate_elements.empty())
277  mooseError("No proximate elements found at node ",
278  closest_node->id(),
279  " at ",
280  static_cast<const Point &>(*closest_node),
281  " on boundary ",
283  ". This should never happen.");
284 
285  for (const Elem * elem : candidate_elements)
286  {
287  for (auto s : elem->side_index_range())
288  if (boundary_info.has_boundary_id(elem, s, _primary_boundary))
289  {
290  located_elem_ids.push_back(elem->id());
291  break;
292  }
293  }
294 
295  if (located_elem_ids.empty())
296  mooseError("No proximate elements found at node ",
297  closest_node->id(),
298  " at ",
299  static_cast<const Point &>(*closest_node),
300  " on boundary ",
302  " share that boundary. This may happen if the mesh uses the same boundary id "
303  "for a nodeset and an unrelated sideset.");
304 
305  closest_elems = &located_elem_ids;
306  }
307  else
308  {
309  auto node_to_elem_pair = _node_to_elem_map.find(closest_node->id());
310  mooseAssert(node_to_elem_pair != _node_to_elem_map.end(),
311  "Missing entry in node to elem map");
312  closest_elems = &(node_to_elem_pair->second);
313  }
314 
315  for (const auto & elem_id : *closest_elems)
316  {
317  const Elem * elem = _mesh.elemPtr(elem_id);
318 
319  std::vector<PenetrationInfo *> thisElemInfo;
320 
321  std::vector<const Node *> nodesThatMustBeOnSide;
322  // If we have a disconnected mesh, we might not have *any*
323  // nodes that must be on a side we check; we'll rely on
324  // boundary info to find valid sides, then rely on comparing
325  // closest points from each to find the best.
326  //
327  // If we don't have a disconnected mesh, then for maximum
328  // backwards compatibility we're still using the older ridge
329  // and peak detection code, which depends on us ruling out
330  // sides that don't touch closest_node.
331  if (!_use_point_locator)
332  nodesThatMustBeOnSide.push_back(closest_node);
334  thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
335  }
336 
337  if (_use_point_locator)
338  {
339  Real min_distance_sq = std::numeric_limits<Real>::max();
340  Point best_point;
341  unsigned int best_i = invalid_uint;
342 
343  // Find closest point in all p_info to the node of interest
344  for (unsigned int i = 0; i < p_info.size(); ++i)
345  {
346  const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
347  const Real distance_sq = (closest_point - node).norm_sq();
348  if (distance_sq < min_distance_sq)
349  {
350  min_distance_sq = distance_sq;
351  best_point = closest_point;
352  best_i = i;
353  }
354  }
355 
356  p_info[best_i]->_closest_point = best_point;
357  p_info[best_i]->_distance =
358  (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) * std::sqrt(min_distance_sq);
360  mooseError("Normal smoothing not implemented with point locator code");
361  Point normal = (best_point - node).unit();
362  const Real dot = normal * p_info[best_i]->_normal;
363  if (dot < 0)
364  normal *= -1;
365  p_info[best_i]->_normal = normal;
366 
367  switchInfo(info, p_info[best_i]);
368  info_set = true;
369  }
370  else
371  {
372  if (p_info.size() == 1)
373  {
374  if (p_info[0]->_tangential_distance <= _tangential_tolerance)
375  {
376  switchInfo(info, p_info[0]);
377  info_set = true;
378  }
379  }
380  else if (p_info.size() > 1)
381  {
382  // Loop through all pairs of faces, and check for contact on ridge between each face pair
383  std::vector<RidgeData> ridgeDataVec;
384  for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
385  for (unsigned int j = i + 1; j < p_info.size(); ++j)
386  {
387  Point closest_coor;
388  Real tangential_distance(0.0);
389  const Node * closest_node_on_ridge = NULL;
390  unsigned int index = 0;
391  Point closest_coor_ref;
392  bool found_ridge_contact_point = findRidgeContactPoint(closest_coor,
393  tangential_distance,
394  closest_node_on_ridge,
395  index,
396  closest_coor_ref,
397  p_info,
398  i,
399  j);
400  if (found_ridge_contact_point)
401  {
402  RidgeData rpd;
403  rpd._closest_coor = closest_coor;
404  rpd._tangential_distance = tangential_distance;
405  rpd._closest_node = closest_node_on_ridge;
406  rpd._index = index;
407  rpd._closest_coor_ref = closest_coor_ref;
408  ridgeDataVec.push_back(rpd);
409  }
410  }
411 
412  if (ridgeDataVec.size() > 0) // Either find the ridge pair that is the best or find a peak
413  {
414  // Group together ridges for which we are off the edge of a common node.
415  // Those are peaks.
416  std::vector<RidgeSetData> ridgeSetDataVec;
417  for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
418  {
419  bool foundSetWithMatchingNode = false;
420  for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
421  {
422  if (ridgeDataVec[i]._closest_node != NULL &&
423  ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
424  {
425  foundSetWithMatchingNode = true;
426  ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
427  break;
428  }
429  }
430  if (!foundSetWithMatchingNode)
431  {
432  RidgeSetData rsd;
434  rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
435  rsd._closest_node = ridgeDataVec[i]._closest_node;
436  ridgeSetDataVec.push_back(rsd);
437  }
438  }
439  // Compute distance to each set of ridges
440  for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
441  {
442  if (ridgeSetDataVec[i]._closest_node !=
443  NULL) // Either a peak or off the edge of single ridge
444  {
445  if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
446  {
447  if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
448  _tangential_tolerance) // off within tolerance
449  {
450  ridgeSetDataVec[i]._closest_coor =
451  ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
452  Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
453  ridgeSetDataVec[i]._distance = contact_point_vec.norm();
454  }
455  }
456  else // several ridges join at common node to make peak. The common node is the
457  // contact point
458  {
459  ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
460  Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
461  ridgeSetDataVec[i]._distance = contact_point_vec.norm();
462  }
463  }
464  else // on a single ridge
465  {
466  ridgeSetDataVec[i]._closest_coor =
467  ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
468  Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
469  ridgeSetDataVec[i]._distance = contact_point_vec.norm();
470  }
471  }
472  // Find the set of ridges closest to us.
473  unsigned int closest_ridge_set_index(0);
474  Real closest_distance(ridgeSetDataVec[0]._distance);
475  Point closest_point(ridgeSetDataVec[0]._closest_coor);
476  for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
477  {
478  if (ridgeSetDataVec[i]._distance < closest_distance)
479  {
480  closest_ridge_set_index = i;
481  closest_distance = ridgeSetDataVec[i]._distance;
482  closest_point = ridgeSetDataVec[i]._closest_coor;
483  }
484  }
485 
486  if (closest_distance <
487  std::numeric_limits<Real>::max()) // contact point is on the closest ridge set
488  {
489  // find the face in the ridge set with the smallest index, assign that one to the
490  // interaction
491  unsigned int face_index(std::numeric_limits<unsigned int>::max());
492  for (unsigned int i = 0;
493  i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
494  ++i)
495  {
496  if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
497  face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
498  }
499 
500  mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
501  "face_index invalid");
502 
503  p_info[face_index]->_closest_point = closest_point;
504  p_info[face_index]->_distance =
505  (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
506  // Calculate the normal as the vector from the ridge to the point only if we're not
507  // doing normal
508  // smoothing. Normal smoothing will average out the normals on its own.
510  {
511  Point normal(closest_point - node);
512  const Real len(normal.norm());
513  if (len > 0)
514  {
515  normal /= len;
516  }
517  const Real dot(normal * p_info[face_index]->_normal);
518  if (dot < 0)
519  {
520  normal *= -1;
521  }
522  p_info[face_index]->_normal = normal;
523  }
524  p_info[face_index]->_tangential_distance = 0.0;
525 
526  Point closest_point_ref;
527  if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
528  1) // contact with a single ridge rather than a peak
529  {
530  p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
531  ._ridge_data_vec[0]
532  ._tangential_distance;
533  p_info[face_index]->_closest_point_ref =
534  ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
535  }
536  else
537  { // peak
538  const Node * closest_node_on_face;
539  bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
540  closest_node_on_face,
541  p_info[face_index]->_side);
542  if (restricted)
543  {
544  if (closest_node_on_face !=
545  ridgeSetDataVec[closest_ridge_set_index]._closest_node)
546  {
547  mooseError("Closest node when restricting point to face != closest node from "
548  "RidgeSetData");
549  }
550  }
551  }
552 
553  FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
554  std::vector<Point> points(1);
555  points[0] = p_info[face_index]->_closest_point_ref;
556  fe->reinit(p_info[face_index]->_side, &points);
557  p_info[face_index]->_side_phi = fe->get_phi();
558  p_info[face_index]->_side_grad_phi = fe->get_dphi();
559  p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
560  p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
561  p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
562 
563  switchInfo(info, p_info[face_index]);
564  info_set = true;
565  }
566  else
567  { // todo:remove invalid ridge cases so they don't mess up individual face
568  // competition????
569  }
570  }
571 
572  if (!info_set) // contact wasn't on a ridge -- compete individual interactions
573  {
574  unsigned int best(0), i(1);
575  do
576  {
577  CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
578  if (CIResult == FIRST_WINS)
579  {
580  i++;
581  }
582  else if (CIResult == SECOND_WINS)
583  {
584  best = i;
585  i++;
586  }
587  else if (CIResult == NEITHER_WINS)
588  {
589  best = i + 1;
590  i += 2;
591  }
592  } while (i < p_info.size() && best < p_info.size());
593  if (best < p_info.size())
594  {
595  // Ensure final info is within the tangential tolerance
596  if (p_info[best]->_tangential_distance <= _tangential_tolerance)
597  {
598  switchInfo(info, p_info[best]);
599  info_set = true;
600  }
601  }
602  }
603  }
604  }
605  }
606 
607  if (!info_set)
608  {
609  // If penetration is not detected within the saved patch, it is possible
610  // that the secondary node has moved outside the saved patch. So, the patch
611  // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
612  // and the penetration detection has to be re-run on the updated patch.
613 
614  _recheck_secondary_nodes.push_back(node_id);
615 
616  delete info;
617  info = NULL;
618  }
619  else
620  {
621  smoothNormal(info, p_info, node);
622  FEBase * fe = _fes[_tid][info->_side->dim()];
623  computeSlip(*fe, *info);
624  }
625 
626  for (unsigned int j = 0; j < p_info.size(); ++j)
627  {
628  if (p_info[j])
629  {
630  delete p_info[j];
631  p_info[j] = NULL;
632  }
633  }
634  }
635 }
636 
637 void
639 {
641  other._recheck_secondary_nodes.begin(),
642  other._recheck_secondary_nodes.end());
643 }
644 
645 void
647 {
648  mooseAssert(infoNew != NULL, "infoNew object is null");
649  if (info)
650  {
651  infoNew->_starting_elem = info->_starting_elem;
652  infoNew->_starting_side_num = info->_starting_side_num;
653  infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
654  infoNew->_incremental_slip = info->_incremental_slip;
655  infoNew->_accumulated_slip = info->_accumulated_slip;
656  infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
657  infoNew->_frictional_energy = info->_frictional_energy;
658  infoNew->_frictional_energy_old = info->_frictional_energy_old;
659  infoNew->_contact_force = info->_contact_force;
660  infoNew->_contact_force_old = info->_contact_force_old;
661  infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
662  infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
663  infoNew->_locked_this_step = info->_locked_this_step;
664  infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
665  infoNew->_mech_status = info->_mech_status;
666  infoNew->_mech_status_old = info->_mech_status_old;
667  }
668  else
669  {
670  infoNew->_starting_elem = infoNew->_elem;
671  infoNew->_starting_side_num = infoNew->_side_num;
673  }
674  delete info;
675  info = infoNew;
676  infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
677 }
678 
681 {
682 
684 
686  pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
687  result = NEITHER_WINS;
688 
689  else if (pi1->_tangential_distance == 0.0 &&
690  pi2->_tangential_distance > 0.0) // on face 1, off face 2
691  result = FIRST_WINS;
692 
693  else if (pi2->_tangential_distance == 0.0 &&
694  pi1->_tangential_distance > 0.0) // on face 2, off face 1
695  result = SECOND_WINS;
696 
697  else if (pi1->_tangential_distance <= _tangential_tolerance &&
698  pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
699  result = FIRST_WINS;
700 
701  else if (pi2->_tangential_distance <= _tangential_tolerance &&
702  pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
703  result = SECOND_WINS;
704 
705  else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
706  result = competeInteractionsBothOnFace(pi1, pi2);
707 
708  else if (pi1->_tangential_distance <= _tangential_tolerance &&
709  pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
710  {
712  if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
713  {
714  // We already checked for ridges, and it got rejected, so neither face must be valid
715  result = NEITHER_WINS;
716  // mooseError("Erroneously encountered ridge case");
717  }
718  else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
719  // the off-side face
720  {
721  if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
722  mooseError("Invalid off_edge_nodes counts");
723 
724  else if (pi1->_off_edge_nodes.size() == 2)
725  result = FIRST_WINS;
726 
727  else if (pi2->_off_edge_nodes.size() == 2)
728  result = SECOND_WINS;
729 
730  else
731  mooseError("Invalid off_edge_nodes counts");
732  }
733  else // The node projects to both faces within tangential tolerance.
734  result = competeInteractionsBothOnFace(pi1, pi2);
735  }
736 
737  return result;
738 }
739 
742 {
744 
745  if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
746  result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
747 
748  else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
749  result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
750 
751  // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
752  // movement. If there is some way to smooth the transition, we should do it.
754  result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
755 
757  result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
758 
759  else // Equal within tolerance. Favor the one with a smaller element id (for repeatibility)
760  {
761  if (pi1->_elem->id() < pi2->_elem->id())
762  result = FIRST_WINS;
763 
764  else
765  result = SECOND_WINS;
766  }
767 
768  return result;
769 }
770 
773 {
774  CommonEdgeResult common_edge(NO_COMMON);
775  const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
776  const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
777  const unsigned dim1 = pi1->_side->dim();
778 
779  if (dim1 == 1)
780  {
781  mooseAssert(pi2->_side->dim() == 1, "Incompatible dimensions.");
782  mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
783  "off_edge_nodes size should be <2 for 2D contact");
784  if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
785  off_edge_nodes1[0] == off_edge_nodes2[0])
786  common_edge = COMMON_EDGE;
787  }
788  else
789  {
790  mooseAssert(dim1 == 2 && pi2->_side->dim() == 2, "Incompatible dimensions.");
791  mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
792  "off_edge_nodes size should be <3 for 3D contact");
793  if (off_edge_nodes1.size() == 1)
794  {
795  if (off_edge_nodes2.size() == 1)
796  {
797  if (off_edge_nodes1[0] == off_edge_nodes2[0])
798  common_edge = COMMON_NODE;
799  }
800  else if (off_edge_nodes2.size() == 2)
801  {
802  if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
803  common_edge = EDGE_AND_COMMON_NODE;
804  }
805  }
806  else if (off_edge_nodes1.size() == 2)
807  {
808  if (off_edge_nodes2.size() == 1)
809  {
810  if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
811  common_edge = EDGE_AND_COMMON_NODE;
812  }
813  else if (off_edge_nodes2.size() == 2)
814  {
815  if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
816  off_edge_nodes1[1] == off_edge_nodes2[1]) ||
817  (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
818  common_edge = COMMON_EDGE;
819  }
820  }
821  }
822  return common_edge;
823 }
824 
825 bool
827  Real & tangential_distance,
828  const Node *& closest_node,
829  unsigned int & index,
830  Point & contact_point_ref,
831  std::vector<PenetrationInfo *> & p_info,
832  const unsigned int index1,
833  const unsigned int index2)
834 {
835  tangential_distance = 0.0;
836  closest_node = NULL;
837  PenetrationInfo * pi1 = p_info[index1];
838  PenetrationInfo * pi2 = p_info[index2];
839  const unsigned sidedim(pi1->_side->dim());
840  mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
841 
842  // Nodes on faces for the two interactions
843  std::vector<const Node *> side1_nodes;
844  getSideCornerNodes(pi1->_side, side1_nodes);
845  std::vector<const Node *> side2_nodes;
846  getSideCornerNodes(pi2->_side, side2_nodes);
847 
848  std::sort(side1_nodes.begin(), side1_nodes.end());
849  std::sort(side2_nodes.begin(), side2_nodes.end());
850 
851  // Find nodes shared by the two faces
852  std::vector<const Node *> common_nodes;
853  std::set_intersection(side1_nodes.begin(),
854  side1_nodes.end(),
855  side2_nodes.begin(),
856  side2_nodes.end(),
857  std::inserter(common_nodes, common_nodes.end()));
858 
859  if (common_nodes.size() != sidedim)
860  return false;
861 
862  bool found_point1, found_point2;
863  Point closest_coor_ref1(pi1->_closest_point_ref);
864  const Node * closest_node1;
865  found_point1 = restrictPointToSpecifiedEdgeOfFace(
866  closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
867 
868  Point closest_coor_ref2(pi2->_closest_point_ref);
869  const Node * closest_node2;
870  found_point2 = restrictPointToSpecifiedEdgeOfFace(
871  closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
872 
873  if (!found_point1 || !found_point2)
874  return false;
875 
876  // if (sidedim == 2)
877  // {
878  // TODO:
879  // We have the parametric coordinates of the closest intersection point for both faces.
880  // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
881  // Find that point by taking dot products of vector from contact point to secondary node point
882  // with face normal vectors to see which face we're closer to.
883  // }
884 
885  FEBase * fe = NULL;
886  std::vector<Point> points(1);
887 
888  // We have to pick one of the two faces to own the contact point. It doesn't really matter
889  // which one we pick. For repeatibility, pick the face with the lowest index.
890  if (index1 < index2)
891  {
892  fe = _fes[_tid][pi1->_side->dim()];
893  contact_point_ref = closest_coor_ref1;
894  points[0] = closest_coor_ref1;
895  fe->reinit(pi1->_side, &points);
896  index = index1;
897  }
898  else
899  {
900  fe = _fes[_tid][pi2->_side->dim()];
901  contact_point_ref = closest_coor_ref2;
902  points[0] = closest_coor_ref2;
903  fe->reinit(pi2->_side, &points);
904  index = index2;
905  }
906 
907  contact_point = fe->get_xyz()[0];
908 
909  if (sidedim == 2)
910  {
911  if (closest_node1) // point is off the ridge between the two elements
912  {
913  mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
914  "If off edge of ridge, closest node must be the same on both elements");
915  closest_node = closest_node1;
916 
917  RealGradient off_face = *closest_node1 - contact_point;
918  tangential_distance = off_face.norm();
919  }
920  }
921 
922  return true;
923 }
924 
925 void
926 PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
927 {
928  const ElemType t(side->type());
929  corner_nodes.clear();
930 
931  corner_nodes.push_back(side->node_ptr(0));
932  corner_nodes.push_back(side->node_ptr(1));
933  switch (t)
934  {
935  case EDGE2:
936  case EDGE3:
937  case EDGE4:
938  {
939  break;
940  }
941 
942  case TRI3:
943  case TRI6:
944  case TRI7:
945  {
946  corner_nodes.push_back(side->node_ptr(2));
947  break;
948  }
949 
950  case QUAD4:
951  case QUAD8:
952  case QUAD9:
953  {
954  corner_nodes.push_back(side->node_ptr(2));
955  corner_nodes.push_back(side->node_ptr(3));
956  break;
957  }
958 
959  default:
960  {
961  mooseError("Unsupported face type: ", t);
962  break;
963  }
964  }
965 }
966 
967 bool
969  const Node *& closest_node,
970  const Elem * side,
971  const std::vector<const Node *> & edge_nodes)
972 {
973  const ElemType t = side->type();
974  Real & xi = p(0);
975  Real & eta = p(1);
976  closest_node = NULL;
977 
978  std::vector<unsigned int> local_node_indices;
979  for (const auto & edge_node : edge_nodes)
980  {
981  unsigned int local_index = side->get_node_index(edge_node);
982  if (local_index == libMesh::invalid_uint)
983  mooseError("Side does not contain node");
984  local_node_indices.push_back(local_index);
985  }
986  mooseAssert(local_node_indices.size() == side->dim(),
987  "Number of edge nodes must match side dimensionality");
988  std::sort(local_node_indices.begin(), local_node_indices.end());
989 
990  bool off_of_this_edge = false;
991 
992  switch (t)
993  {
994  case EDGE2:
995  case EDGE3:
996  case EDGE4:
997  {
998  if (local_node_indices[0] == 0)
999  {
1000  if (xi <= -1.0)
1001  {
1002  xi = -1.0;
1003  off_of_this_edge = true;
1004  closest_node = side->node_ptr(0);
1005  }
1006  }
1007  else if (local_node_indices[0] == 1)
1008  {
1009  if (xi >= 1.0)
1010  {
1011  xi = 1.0;
1012  off_of_this_edge = true;
1013  closest_node = side->node_ptr(1);
1014  }
1015  }
1016  else
1017  {
1018  mooseError("Invalid local node indices");
1019  }
1020  break;
1021  }
1022 
1023  case TRI3:
1024  case TRI6:
1025  case TRI7:
1026  {
1027  if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1028  {
1029  if (eta <= 0.0)
1030  {
1031  eta = 0.0;
1032  off_of_this_edge = true;
1033  if (xi < 0.0)
1034  closest_node = side->node_ptr(0);
1035  else if (xi > 1.0)
1036  closest_node = side->node_ptr(1);
1037  }
1038  }
1039  else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1040  {
1041  if ((xi + eta) > 1.0)
1042  {
1043  Real delta = (xi + eta - 1.0) / 2.0;
1044  xi -= delta;
1045  eta -= delta;
1046  off_of_this_edge = true;
1047  if (xi > 1.0)
1048  closest_node = side->node_ptr(1);
1049  else if (xi < 0.0)
1050  closest_node = side->node_ptr(2);
1051  }
1052  }
1053  else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1054  {
1055  if (xi <= 0.0)
1056  {
1057  xi = 0.0;
1058  off_of_this_edge = true;
1059  if (eta > 1.0)
1060  closest_node = side->node_ptr(2);
1061  else if (eta < 0.0)
1062  closest_node = side->node_ptr(0);
1063  }
1064  }
1065  else
1066  {
1067  mooseError("Invalid local node indices");
1068  }
1069 
1070  break;
1071  }
1072 
1073  case QUAD4:
1074  case QUAD8:
1075  case QUAD9:
1076  {
1077  if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1078  {
1079  if (eta <= -1.0)
1080  {
1081  eta = -1.0;
1082  off_of_this_edge = true;
1083  if (xi < -1.0)
1084  closest_node = side->node_ptr(0);
1085  else if (xi > 1.0)
1086  closest_node = side->node_ptr(1);
1087  }
1088  }
1089  else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1090  {
1091  if (xi >= 1.0)
1092  {
1093  xi = 1.0;
1094  off_of_this_edge = true;
1095  if (eta < -1.0)
1096  closest_node = side->node_ptr(1);
1097  else if (eta > 1.0)
1098  closest_node = side->node_ptr(2);
1099  }
1100  }
1101  else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1102  {
1103  if (eta >= 1.0)
1104  {
1105  eta = 1.0;
1106  off_of_this_edge = true;
1107  if (xi < -1.0)
1108  closest_node = side->node_ptr(3);
1109  else if (xi > 1.0)
1110  closest_node = side->node_ptr(2);
1111  }
1112  }
1113  else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1114  {
1115  if (xi <= -1.0)
1116  {
1117  xi = -1.0;
1118  off_of_this_edge = true;
1119  if (eta < -1.0)
1120  closest_node = side->node_ptr(0);
1121  else if (eta > 1.0)
1122  closest_node = side->node_ptr(3);
1123  }
1124  }
1125  else
1126  {
1127  mooseError("Invalid local node indices");
1128  }
1129  break;
1130  }
1131 
1132  default:
1133  {
1134  mooseError("Unsupported face type: ", t);
1135  break;
1136  }
1137  }
1138  return off_of_this_edge;
1139 }
1140 
1141 bool
1142 PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
1143 {
1144  const ElemType t(side->type());
1145  Real & xi = p(0);
1146  Real & eta = p(1);
1147  closest_node = NULL;
1148 
1149  bool off_of_this_face(false);
1150 
1151  switch (t)
1152  {
1153  case EDGE2:
1154  case EDGE3:
1155  case EDGE4:
1156  {
1157  if (xi < -1.0)
1158  {
1159  xi = -1.0;
1160  off_of_this_face = true;
1161  closest_node = side->node_ptr(0);
1162  }
1163  else if (xi > 1.0)
1164  {
1165  xi = 1.0;
1166  off_of_this_face = true;
1167  closest_node = side->node_ptr(1);
1168  }
1169  break;
1170  }
1171 
1172  case TRI3:
1173  case TRI6:
1174  case TRI7:
1175  {
1176  if (eta < 0.0)
1177  {
1178  eta = 0.0;
1179  off_of_this_face = true;
1180  if (xi < 0.5)
1181  {
1182  closest_node = side->node_ptr(0);
1183  if (xi < 0.0)
1184  xi = 0.0;
1185  }
1186  else
1187  {
1188  closest_node = side->node_ptr(1);
1189  if (xi > 1.0)
1190  xi = 1.0;
1191  }
1192  }
1193  else if ((xi + eta) > 1.0)
1194  {
1195  Real delta = (xi + eta - 1.0) / 2.0;
1196  xi -= delta;
1197  eta -= delta;
1198  off_of_this_face = true;
1199  if (xi > 0.5)
1200  {
1201  closest_node = side->node_ptr(1);
1202  if (xi > 1.0)
1203  {
1204  xi = 1.0;
1205  eta = 0.0;
1206  }
1207  }
1208  else
1209  {
1210  closest_node = side->node_ptr(2);
1211  if (xi < 0.0)
1212  {
1213  xi = 0.0;
1214  eta = 1.0;
1215  }
1216  }
1217  }
1218  else if (xi < 0.0)
1219  {
1220  xi = 0.0;
1221  off_of_this_face = true;
1222  if (eta > 0.5)
1223  {
1224  closest_node = side->node_ptr(2);
1225  if (eta > 1.0)
1226  eta = 1.0;
1227  }
1228  else
1229  {
1230  closest_node = side->node_ptr(0);
1231  if (eta < 0.0)
1232  eta = 0.0;
1233  }
1234  }
1235  break;
1236  }
1237 
1238  case QUAD4:
1239  case QUAD8:
1240  case QUAD9:
1241  {
1242  if (eta < -1.0)
1243  {
1244  eta = -1.0;
1245  off_of_this_face = true;
1246  if (xi < 0.0)
1247  {
1248  closest_node = side->node_ptr(0);
1249  if (xi < -1.0)
1250  xi = -1.0;
1251  }
1252  else
1253  {
1254  closest_node = side->node_ptr(1);
1255  if (xi > 1.0)
1256  xi = 1.0;
1257  }
1258  }
1259  else if (xi > 1.0)
1260  {
1261  xi = 1.0;
1262  off_of_this_face = true;
1263  if (eta < 0.0)
1264  {
1265  closest_node = side->node_ptr(1);
1266  if (eta < -1.0)
1267  eta = -1.0;
1268  }
1269  else
1270  {
1271  closest_node = side->node_ptr(2);
1272  if (eta > 1.0)
1273  eta = 1.0;
1274  }
1275  }
1276  else if (eta > 1.0)
1277  {
1278  eta = 1.0;
1279  off_of_this_face = true;
1280  if (xi < 0.0)
1281  {
1282  closest_node = side->node_ptr(3);
1283  if (xi < -1.0)
1284  xi = -1.0;
1285  }
1286  else
1287  {
1288  closest_node = side->node_ptr(2);
1289  if (xi > 1.0)
1290  xi = 1.0;
1291  }
1292  }
1293  else if (xi < -1.0)
1294  {
1295  xi = -1.0;
1296  off_of_this_face = true;
1297  if (eta < 0.0)
1298  {
1299  closest_node = side->node_ptr(0);
1300  if (eta < -1.0)
1301  eta = -1.0;
1302  }
1303  else
1304  {
1305  closest_node = side->node_ptr(3);
1306  if (eta > 1.0)
1307  eta = 1.0;
1308  }
1309  }
1310  break;
1311  }
1312 
1313  default:
1314  {
1315  mooseError("Unsupported face type: ", t);
1316  break;
1317  }
1318  }
1319  return off_of_this_face;
1320 }
1321 
1322 bool
1324  const Elem * side,
1325  FEBase * fe,
1326  const Point * secondary_point,
1327  const Real tangential_tolerance)
1328 {
1329  unsigned int dim = primary_elem->dim();
1330 
1331  const std::vector<Point> & phys_point = fe->get_xyz();
1332 
1333  const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1334  const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1335 
1336  Point ref_point;
1337 
1338  std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1339 
1340  fe->reinit(side, &points);
1341 
1342  RealGradient d = *secondary_point - phys_point[0];
1343 
1344  const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1345  Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1346 
1347  RealVectorValue normal;
1348  if (dim - 1 == 2)
1349  {
1350  normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1351  }
1352  else if (dim - 1 == 1)
1353  {
1354  const Node * const * elem_nodes = primary_elem->get_nodes();
1355  const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1356  const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1357 
1358  Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1359  out_of_plane_normal /= out_of_plane_normal.norm();
1360 
1361  normal = dxyz_dxi[0].cross(out_of_plane_normal);
1362  }
1363  else
1364  {
1365  return true;
1366  }
1367  normal /= normal.norm();
1368 
1369  const Real dot(d * normal);
1370 
1371  const RealGradient normcomp = dot * normal;
1372  const RealGradient tangcomp = d - normcomp;
1373 
1374  const Real tangdist = tangcomp.norm();
1375 
1376  // Increase the size of the zone that we consider if the vector from the face
1377  // to the node has a larger normal component
1378  const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1379 
1380  bool isReasonableCandidate = true;
1381  if (tangdist > faceExpansionFactor * max_face_length)
1382  {
1383  isReasonableCandidate = false;
1384  }
1385  return isReasonableCandidate;
1386 }
1387 
1388 void
1390 {
1391  // Slip is current projected position of secondary node minus
1392  // original projected position of secondary node
1393  std::vector<Point> points(1);
1394  points[0] = info._starting_closest_point_ref;
1395  const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
1396  fe.reinit(&side, &points);
1397  const std::vector<Point> & starting_point = fe.get_xyz();
1398  info._incremental_slip = info._closest_point - starting_point[0];
1399  if (info.isCaptured())
1400  {
1401  info._frictional_energy =
1402  info._frictional_energy_old + info._contact_force * info._incremental_slip;
1403  info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
1404  }
1405 }
1406 
1407 void
1409  std::vector<PenetrationInfo *> & p_info,
1410  const Node & node)
1411 {
1413  {
1415  {
1416  // If we are within the smoothing distance of any edges or corners, find the
1417  // corner nodes for those edges/corners, and weights from distance to edge/corner
1418  std::vector<Real> edge_face_weights;
1419  std::vector<PenetrationInfo *> edge_face_info;
1420 
1421  getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
1422 
1423  mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1424  "edge_face_info.size() != edge_face_weights.size()");
1425 
1426  if (edge_face_info.size() > 0)
1427  {
1428  // Smooth the normal using the weighting functions for all participating faces.
1429  RealVectorValue new_normal;
1430  Real this_face_weight = 1.0;
1431 
1432  for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1433  {
1434  PenetrationInfo * npi = edge_face_info[efwi];
1435  if (npi)
1436  new_normal += npi->_normal * edge_face_weights[efwi];
1437 
1438  this_face_weight -= edge_face_weights[efwi];
1439  }
1440  mooseAssert(this_face_weight >= (0.25 - 1e-8),
1441  "Sum of weights of other faces shouldn't exceed 0.75");
1442  new_normal += info->_normal * this_face_weight;
1443 
1444  const Real len = new_normal.norm();
1445  if (len > 0)
1446  new_normal /= len;
1447 
1448  info->_normal = new_normal;
1449  }
1450  }
1452  {
1453  // params.addParam<VariableName>("var_name","description");
1454  // getParam<VariableName>("var_name")
1455  info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1456  info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1457  info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1458  const Real len(info->_normal.norm());
1459  if (len > 0)
1460  info->_normal /= len;
1461  }
1462  }
1463 }
1464 
1465 void
1467  std::vector<PenetrationInfo *> & edge_face_info,
1468  std::vector<Real> & edge_face_weights,
1469  std::vector<PenetrationInfo *> & p_info,
1470  const Node & secondary_node)
1471 {
1472  const Elem * side = info->_side;
1473  const Point & p = info->_closest_point_ref;
1474  std::set<dof_id_type> elems_to_exclude;
1475  elems_to_exclude.insert(info->_elem->id());
1476 
1477  std::vector<std::vector<const Node *>> edge_nodes;
1478 
1479  // Get the pairs of nodes along every edge that we are close enough to smooth with
1480  getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1481  std::vector<Elem *> edge_neighbor_elems;
1482  edge_face_info.resize(edge_nodes.size(), NULL);
1483 
1484  std::vector<unsigned int> edges_without_neighbors;
1485 
1486  for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1487  {
1488  // Sort all sets of edge nodes (needed for comparing edges)
1489  std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1490 
1491  std::vector<PenetrationInfo *> face_info_comm_edge;
1493  &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1494 
1495  if (face_info_comm_edge.size() == 0)
1496  edges_without_neighbors.push_back(i);
1497  else if (face_info_comm_edge.size() > 1)
1498  mooseError("Only one neighbor allowed per edge");
1499  else
1500  edge_face_info[i] = face_info_comm_edge[0];
1501  }
1502 
1503  // Remove edges without neighbors from the vector, starting from end
1504  std::vector<unsigned int>::reverse_iterator rit;
1505  for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1506  {
1507  unsigned int index = *rit;
1508  edge_nodes.erase(edge_nodes.begin() + index);
1509  edge_face_weights.erase(edge_face_weights.begin() + index);
1510  edge_face_info.erase(edge_face_info.begin() + index);
1511  }
1512 
1513  // Handle corner case
1514  if (edge_nodes.size() > 1)
1515  {
1516  if (edge_nodes.size() != 2)
1517  mooseError("Invalid number of smoothing edges");
1518 
1519  // find common node
1520  std::vector<const Node *> common_nodes;
1521  std::set_intersection(edge_nodes[0].begin(),
1522  edge_nodes[0].end(),
1523  edge_nodes[1].begin(),
1524  edge_nodes[1].end(),
1525  std::inserter(common_nodes, common_nodes.end()));
1526 
1527  if (common_nodes.size() != 1)
1528  mooseError("Invalid number of common nodes");
1529 
1530  for (const auto & pinfo : edge_face_info)
1531  elems_to_exclude.insert(pinfo->_elem->id());
1532 
1533  std::vector<PenetrationInfo *> face_info_comm_edge;
1535  &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1536 
1537  unsigned int num_corner_neighbors = face_info_comm_edge.size();
1538 
1539  if (num_corner_neighbors > 0)
1540  {
1541  Real fw0 = edge_face_weights[0];
1542  Real fw1 = edge_face_weights[1];
1543 
1544  // Corner weight is product of edge weights. Spread out over multiple neighbors.
1545  Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1546 
1547  // Adjust original edge weights
1548  edge_face_weights[0] *= (1.0 - fw1);
1549  edge_face_weights[1] *= (1.0 - fw0);
1550 
1551  for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1552  {
1553  edge_face_weights.push_back(fw_corner);
1554  edge_face_info.push_back(face_info_comm_edge[i]);
1555  }
1556  }
1557  }
1558 }
1559 
1560 void
1562  const Point & p,
1563  const Elem * side,
1564  std::vector<std::vector<const Node *>> & edge_nodes,
1565  std::vector<Real> & edge_face_weights)
1566 {
1567  const ElemType t(side->type());
1568  const Real & xi = p(0);
1569  const Real & eta = p(1);
1570 
1571  Real smooth_limit = 1.0 - _normal_smoothing_distance;
1572 
1573  switch (t)
1574  {
1575  case EDGE2:
1576  case EDGE3:
1577  case EDGE4:
1578  {
1579  if (xi < -smooth_limit)
1580  {
1581  std::vector<const Node *> en;
1582  en.push_back(side->node_ptr(0));
1583  edge_nodes.push_back(en);
1584  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1585  if (fw > 0.5)
1586  fw = 0.5;
1587  edge_face_weights.push_back(fw);
1588  }
1589  else if (xi > smooth_limit)
1590  {
1591  std::vector<const Node *> en;
1592  en.push_back(side->node_ptr(1));
1593  edge_nodes.push_back(en);
1594  Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1595  if (fw > 0.5)
1596  fw = 0.5;
1597  edge_face_weights.push_back(fw);
1598  }
1599  break;
1600  }
1601 
1602  case TRI3:
1603  case TRI6:
1604  case TRI7:
1605  {
1606  if (eta < -smooth_limit)
1607  {
1608  std::vector<const Node *> en;
1609  en.push_back(side->node_ptr(0));
1610  en.push_back(side->node_ptr(1));
1611  edge_nodes.push_back(en);
1612  Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1613  if (fw > 0.5)
1614  fw = 0.5;
1615  edge_face_weights.push_back(fw);
1616  }
1617  if ((xi + eta) > smooth_limit)
1618  {
1619  std::vector<const Node *> en;
1620  en.push_back(side->node_ptr(1));
1621  en.push_back(side->node_ptr(2));
1622  edge_nodes.push_back(en);
1623  Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
1624  if (fw > 0.5)
1625  fw = 0.5;
1626  edge_face_weights.push_back(fw);
1627  }
1628  if (xi < -smooth_limit)
1629  {
1630  std::vector<const Node *> en;
1631  en.push_back(side->node_ptr(2));
1632  en.push_back(side->node_ptr(0));
1633  edge_nodes.push_back(en);
1634  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1635  if (fw > 0.5)
1636  fw = 0.5;
1637  edge_face_weights.push_back(fw);
1638  }
1639  break;
1640  }
1641 
1642  case QUAD4:
1643  case QUAD8:
1644  case QUAD9:
1645  {
1646  if (eta < -smooth_limit)
1647  {
1648  std::vector<const Node *> en;
1649  en.push_back(side->node_ptr(0));
1650  en.push_back(side->node_ptr(1));
1651  edge_nodes.push_back(en);
1652  Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1653  if (fw > 0.5)
1654  fw = 0.5;
1655  edge_face_weights.push_back(fw);
1656  }
1657  if (xi > smooth_limit)
1658  {
1659  std::vector<const Node *> en;
1660  en.push_back(side->node_ptr(1));
1661  en.push_back(side->node_ptr(2));
1662  edge_nodes.push_back(en);
1663  Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1664  if (fw > 0.5)
1665  fw = 0.5;
1666  edge_face_weights.push_back(fw);
1667  }
1668  if (eta > smooth_limit)
1669  {
1670  std::vector<const Node *> en;
1671  en.push_back(side->node_ptr(2));
1672  en.push_back(side->node_ptr(3));
1673  edge_nodes.push_back(en);
1674  Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1675  if (fw > 0.5)
1676  fw = 0.5;
1677  edge_face_weights.push_back(fw);
1678  }
1679  if (xi < -smooth_limit)
1680  {
1681  std::vector<const Node *> en;
1682  en.push_back(side->node_ptr(3));
1683  en.push_back(side->node_ptr(0));
1684  edge_nodes.push_back(en);
1685  Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1686  if (fw > 0.5)
1687  fw = 0.5;
1688  edge_face_weights.push_back(fw);
1689  }
1690  break;
1691  }
1692 
1693  default:
1694  {
1695  mooseError("Unsupported face type: ", t);
1696  break;
1697  }
1698  }
1699 }
1700 
1701 void
1703  const Node * secondary_node,
1704  const std::set<dof_id_type> & elems_to_exclude,
1705  const std::vector<const Node *> edge_nodes,
1706  std::vector<PenetrationInfo *> & face_info_comm_edge,
1707  std::vector<PenetrationInfo *> & p_info)
1708 {
1709  // elems connected to a node on this edge, find one that has the same corners as this, and is not
1710  // the current elem
1711  auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
1712  // nodes
1713  mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
1714  const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1715 
1716  std::vector<const Elem *> elems_connected_to_edge;
1717 
1718  for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1719  {
1720  if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1721  continue;
1722  const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1723 
1724  std::vector<const Node *> nodevec;
1725  for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1726  if (elem->is_vertex(ni))
1727  nodevec.push_back(elem->node_ptr(ni));
1728 
1729  std::vector<const Node *> common_nodes;
1730  std::sort(nodevec.begin(), nodevec.end());
1731  std::set_intersection(edge_nodes.begin(),
1732  edge_nodes.end(),
1733  nodevec.begin(),
1734  nodevec.end(),
1735  std::inserter(common_nodes, common_nodes.end()));
1736 
1737  if (common_nodes.size() == edge_nodes.size())
1738  elems_connected_to_edge.push_back(elem);
1739  }
1740 
1741  if (elems_connected_to_edge.size() > 0)
1742  {
1743 
1744  // There are potentially multiple elements that share a common edge
1745  // 2D:
1746  // There can only be one element on the same surface
1747  // 3D:
1748  // If there are two edge nodes, there can only be one element on the same surface
1749  // If there is only one edge node (a corner), there could be multiple elements on the same
1750  // surface
1751  bool allowMultipleNeighbors = false;
1752 
1753  if (elems_connected_to_edge[0]->dim() == 3)
1754  {
1755  if (edge_nodes.size() == 1)
1756  {
1757  allowMultipleNeighbors = true;
1758  }
1759  }
1760 
1761  for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1762  {
1763  std::vector<PenetrationInfo *> thisElemInfo;
1764  getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1765  if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1766  {
1767  if (thisElemInfo.size() > 1)
1768  mooseError(
1769  "Found multiple neighbors to current edge/face on surface when only one is allowed");
1770  face_info_comm_edge.push_back(thisElemInfo[0]);
1771  break;
1772  }
1773 
1775  thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1776  if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1777  {
1778  if (thisElemInfo.size() > 1)
1779  mooseError(
1780  "Found multiple neighbors to current edge/face on surface when only one is allowed");
1781  face_info_comm_edge.push_back(thisElemInfo[0]);
1782  break;
1783  }
1784 
1785  for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1786  face_info_comm_edge.push_back(thisElemInfo[j]);
1787  }
1788  }
1789 }
1790 
1791 void
1792 PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1793  std::vector<PenetrationInfo *> & p_info,
1794  const Elem * elem)
1795 {
1796  for (const auto & pi : p_info)
1797  {
1798  if (!pi)
1799  continue;
1800 
1801  if (pi->_elem == elem)
1802  thisElemInfo.push_back(pi);
1803  }
1804 }
1805 
1806 void
1807 PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1808  std::vector<PenetrationInfo *> & p_info,
1809  const Node * secondary_node,
1810  const Elem * elem,
1811  const std::vector<const Node *> & nodes_that_must_be_on_side,
1812  const bool check_whether_reasonable)
1813 {
1814  const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
1815 
1816  for (auto s : elem->side_index_range())
1817  {
1818  if (!boundary_info.has_boundary_id(elem, s, _primary_boundary))
1819  continue;
1820 
1821  // Don't create info for this side if one already exists
1822  bool already_have_info_this_side = false;
1823  for (const auto & pi : thisElemInfo)
1824  if (pi->_side_num == s)
1825  {
1826  already_have_info_this_side = true;
1827  break;
1828  }
1829 
1830  if (already_have_info_this_side)
1831  break;
1832 
1833  const Elem * side = elem->build_side_ptr(s).release();
1834 
1835  // Only continue with creating info for this side if the side contains
1836  // all of the nodes in nodes_that_must_be_on_side
1837  std::vector<const Node *> nodevec;
1838  for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1839  nodevec.push_back(side->node_ptr(ni));
1840 
1841  std::sort(nodevec.begin(), nodevec.end());
1842  std::vector<const Node *> common_nodes;
1843  std::set_intersection(nodes_that_must_be_on_side.begin(),
1844  nodes_that_must_be_on_side.end(),
1845  nodevec.begin(),
1846  nodevec.end(),
1847  std::inserter(common_nodes, common_nodes.end()));
1848  if (common_nodes.size() != nodes_that_must_be_on_side.size())
1849  {
1850  delete side;
1851  break;
1852  }
1853 
1854  FEBase * fe_elem = _fes[_tid][elem->dim()];
1855  FEBase * fe_side = _fes[_tid][side->dim()];
1856 
1857  // Optionally check to see whether face is reasonable candidate based on an
1858  // estimate of how closely it is likely to project to the face
1859  if (check_whether_reasonable)
1860  if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1861  {
1862  delete side;
1863  break;
1864  }
1865 
1866  Point contact_phys;
1867  Point contact_ref;
1868  Point contact_on_face_ref;
1869  Real distance = 0.;
1870  Real tangential_distance = 0.;
1871  RealGradient normal;
1872  bool contact_point_on_side;
1873  std::vector<const Node *> off_edge_nodes;
1874  std::vector<std::vector<Real>> side_phi;
1875  std::vector<std::vector<RealGradient>> side_grad_phi;
1876  std::vector<RealGradient> dxyzdxi;
1877  std::vector<RealGradient> dxyzdeta;
1878  std::vector<RealGradient> d2xyzdxideta;
1879 
1880  std::unique_ptr<PenetrationInfo> pen_info =
1881  std::make_unique<PenetrationInfo>(elem,
1882  side,
1883  s,
1884  normal,
1885  distance,
1886  tangential_distance,
1887  contact_phys,
1888  contact_ref,
1889  contact_on_face_ref,
1890  off_edge_nodes,
1891  side_phi,
1892  side_grad_phi,
1893  dxyzdxi,
1894  dxyzdeta,
1895  d2xyzdxideta);
1896 
1897  bool search_succeeded = false;
1898  Moose::findContactPoint(*pen_info,
1899  fe_elem,
1900  fe_side,
1901  _fe_type,
1902  *secondary_node,
1903  true,
1905  contact_point_on_side,
1906  search_succeeded);
1907 
1908  // Do not add contact info from failed searches
1909  if (search_succeeded)
1910  {
1911  thisElemInfo.push_back(pen_info.get());
1912  p_info.push_back(pen_info.release());
1913  }
1914  }
1915 }
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:3134
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:162
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:3469
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:3754
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:631
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