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