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