https://mooseframework.inl.gov
Classes | Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Protected Attributes | List of all members
PenetrationThread Class Reference

#include <PenetrationThread.h>

Classes

struct  RidgeData
 
struct  RidgeSetData
 

Public Member Functions

 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)
 
 PenetrationThread (PenetrationThread &x, Threads::split split)
 
void operator() (const NodeIdRange &range)
 
void join (const PenetrationThread &other)
 

Public Attributes

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 patch has to be updated. More...
 

Protected Types

enum  CompeteInteractionResult { FIRST_WINS, SECOND_WINS, NEITHER_WINS }
 
enum  CommonEdgeResult { NO_COMMON, COMMON_EDGE, COMMON_NODE, EDGE_AND_COMMON_NODE }
 

Protected Member Functions

CompeteInteractionResult competeInteractions (PenetrationInfo *pi1, PenetrationInfo *pi2)
 When interactions are identified between a node and two faces, compete between the faces to determine whether first (pi1) or second (pi2) interaction is stronger. More...
 
CompeteInteractionResult competeInteractionsBothOnFace (PenetrationInfo *pi1, PenetrationInfo *pi2)
 Determine whether first (pi1) or second (pi2) interaction is stronger when it is known that the node projects to both of the two competing faces. More...
 
CommonEdgeResult interactionsOffCommonEdge (PenetrationInfo *pi1, PenetrationInfo *pi2)
 
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)
 
void getSideCornerNodes (const Elem *side, std::vector< const Node *> &corner_nodes)
 
bool restrictPointToSpecifiedEdgeOfFace (libMesh::Point &p, const Node *&closest_node, const Elem *side, const std::vector< const Node *> &edge_nodes)
 
bool restrictPointToFace (libMesh::Point &p, const Node *&closest_node, const Elem *side)
 
bool isFaceReasonableCandidate (const Elem *primary_elem, const Elem *side, libMesh::FEBase *fe, const libMesh::Point *secondary_point, const Real tangential_tolerance)
 
void smoothNormal (PenetrationInfo *info, std::vector< PenetrationInfo *> &p_info, const Node &node)
 
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)
 
void getSmoothingEdgeNodesAndWeights (const libMesh::Point &p, const Elem *side, std::vector< std::vector< const Node *>> &edge_nodes, std::vector< Real > &edge_face_weights)
 
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)
 
void getInfoForElem (std::vector< PenetrationInfo *> &thisElemInfo, std::vector< PenetrationInfo *> &p_info, const Elem *elem)
 
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)
 
void getSidesOnPrimaryBoundary (std::vector< unsigned int > &sides, const Elem *const elem)
 
void computeSlip (libMesh::FEBase &fe, PenetrationInfo &info)
 
void switchInfo (PenetrationInfo *&info, PenetrationInfo *&infoNew)
 

Protected Attributes

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
 
MooseVariable_nodal_normal_x
 
MooseVariable_nodal_normal_y
 
MooseVariable_nodal_normal_z
 
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
 
THREAD_ID _tid
 
libMesh::ElemSideBuilder _elem_side_builder
 Helper for building element sides without extraneous allocation. More...
 

Detailed Description

Definition at line 24 of file PenetrationThread.h.

Member Enumeration Documentation

◆ CommonEdgeResult

Enumerator
NO_COMMON 
COMMON_EDGE 
COMMON_NODE 
EDGE_AND_COMMON_NODE 

Definition at line 98 of file PenetrationThread.h.

◆ CompeteInteractionResult

Enumerator
FIRST_WINS 
SECOND_WINS 
NEITHER_WINS 

Definition at line 91 of file PenetrationThread.h.

Constructor & Destructor Documentation

◆ PenetrationThread() [1/2]

PenetrationThread::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 
)

Definition at line 29 of file PenetrationThread.C.

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 }
MooseVariable * _nodal_normal_z
NearestNodeLocator & _nearest_node
const std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > & _bc_tuples
libMesh::FEType & _fe_type
BoundaryID _primary_boundary
SubProblem & _subproblem
const MooseMesh & _mesh
BoundaryID _secondary_boundary
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
MooseVariable * _nodal_normal_y
const std::map< dof_id_type, std::vector< dof_id_type > > & _node_to_elem_map
MooseVariable * _nodal_normal_x
std::vector< std::vector< libMesh::FEBase * > > & _fes
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method

◆ PenetrationThread() [2/2]

PenetrationThread::PenetrationThread ( PenetrationThread x,
Threads::split  split 
)

Definition at line 69 of file PenetrationThread.C.

71  _mesh(x._mesh),
81  _fes(x._fes),
82  _fe_type(x._fe_type),
86 {
87 }
NearestNodeLocator & _nearest_node
const std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > & _bc_tuples
libMesh::FEType & _fe_type
BoundaryID _primary_boundary
SubProblem & _subproblem
const MooseMesh & _mesh
BoundaryID _secondary_boundary
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
const std::map< dof_id_type, std::vector< dof_id_type > > & _node_to_elem_map
std::vector< std::vector< libMesh::FEBase * > > & _fes
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method

Member Function Documentation

◆ competeInteractions()

PenetrationThread::CompeteInteractionResult PenetrationThread::competeInteractions ( PenetrationInfo pi1,
PenetrationInfo pi2 
)
protected

When interactions are identified between a node and two faces, compete between the faces to determine whether first (pi1) or second (pi2) interaction is stronger.

Parameters
pi1Pointer to the PenetrationInfo object for the first face
pi2Pointer to the PenetrationInfo object for the second face
Returns
Appropriate ComputeInterationResult enum entry identifying which face is a better match

Definition at line 514 of file PenetrationThread.C.

Referenced by operator()().

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 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
std::vector< const Node * > _off_edge_nodes
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)

◆ competeInteractionsBothOnFace()

PenetrationThread::CompeteInteractionResult PenetrationThread::competeInteractionsBothOnFace ( PenetrationInfo pi1,
PenetrationInfo pi2 
)
protected

Determine whether first (pi1) or second (pi2) interaction is stronger when it is known that the node projects to both of the two competing faces.

Parameters
pi1Pointer to the PenetrationInfo object for the first face
pi2Pointer to the PenetrationInfo object for the second face
Returns
Appropriate ComputeInterationResult enum entry identifying which face is a better match

Definition at line 575 of file PenetrationThread.C.

Referenced by competeInteractions().

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 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:42
dof_id_type id() const
const Elem * _elem
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

◆ computeSlip()

void PenetrationThread::computeSlip ( libMesh::FEBase fe,
PenetrationInfo info 
)
protected

Definition at line 1223 of file PenetrationThread.C.

Referenced by operator()().

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 }
MPI_Info info
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
virtual_for_inffe const std::vector< Point > & get_xyz() const
libMesh::ElemSideBuilder _elem_side_builder
Helper for building element sides without extraneous allocation.

◆ createInfoForElem()

void PenetrationThread::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 
)
protected

Definition at line 1641 of file PenetrationThread.C.

Referenced by getInfoForFacesWithCommonNodes(), and operator()().

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 }
void getSidesOnPrimaryBoundary(std::vector< unsigned int > &sides, const Elem *const elem)
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
Real distance(const Point &p)
libMesh::FEType & _fe_type
virtual unsigned int n_nodes() const=0
char ** sides
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...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const=0
const Node * node_ptr(const unsigned int i) const
std::vector< std::vector< libMesh::FEBase * > > & _fes
bool isFaceReasonableCandidate(const Elem *primary_elem, const Elem *side, libMesh::FEBase *fe, const libMesh::Point *secondary_point, const Real tangential_tolerance)
const Real pi

◆ findRidgeContactPoint()

bool PenetrationThread::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 
)
protected

Definition at line 660 of file PenetrationThread.C.

Referenced by operator()().

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 }
auto norm() const -> decltype(std::norm(Real()))
Data structure used to hold penetration information.
void getSideCornerNodes(const Elem *side, std::vector< const Node *> &corner_nodes)
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
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
virtual unsigned short dim() const=0
std::vector< std::vector< libMesh::FEBase * > > & _fes

◆ getInfoForElem()

void PenetrationThread::getInfoForElem ( std::vector< PenetrationInfo *> &  thisElemInfo,
std::vector< PenetrationInfo *> &  p_info,
const Elem elem 
)
protected

Definition at line 1626 of file PenetrationThread.C.

Referenced by getInfoForFacesWithCommonNodes().

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 }
const Real pi

◆ getInfoForFacesWithCommonNodes()

void PenetrationThread::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 
)
protected

Definition at line 1536 of file PenetrationThread.C.

Referenced by getSmoothingFacesAndWeights().

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 }
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:3108
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
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)
const MooseMesh & _mesh
virtual unsigned int n_nodes() const=0
const Node * node_ptr(const unsigned int i) const
virtual bool is_vertex(const unsigned int i) const=0
void getInfoForElem(std::vector< PenetrationInfo *> &thisElemInfo, std::vector< PenetrationInfo *> &p_info, const Elem *elem)
const std::map< dof_id_type, std::vector< dof_id_type > > & _node_to_elem_map

◆ getSideCornerNodes()

void PenetrationThread::getSideCornerNodes ( const Elem side,
std::vector< const Node *> &  corner_nodes 
)
protected

Definition at line 760 of file PenetrationThread.C.

Referenced by findRidgeContactPoint().

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 }
ElemType
QUAD8
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
EDGE4
TRI3
QUAD4
TRI6
EDGE2
const Node * node_ptr(const unsigned int i) const
TRI7
QUAD9
EDGE3
virtual ElemType type() const=0

◆ getSidesOnPrimaryBoundary()

void PenetrationThread::getSidesOnPrimaryBoundary ( std::vector< unsigned int > &  sides,
const Elem *const  elem 
)
protected

Definition at line 1760 of file PenetrationThread.C.

Referenced by createInfoForElem().

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 }
std::tuple< dof_id_type, unsigned short int, boundary_id_type > BCTuple
const std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > & _bc_tuples
BoundaryID _primary_boundary
int8_t boundary_id_type
dof_id_type id() const
char ** sides
void operator()(const NodeIdRange &range)
uint8_t dof_id_type

◆ getSmoothingEdgeNodesAndWeights()

void PenetrationThread::getSmoothingEdgeNodesAndWeights ( const libMesh::Point p,
const Elem side,
std::vector< std::vector< const Node *>> &  edge_nodes,
std::vector< Real > &  edge_face_weights 
)
protected

Definition at line 1395 of file PenetrationThread.C.

Referenced by getSmoothingFacesAndWeights().

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 }
ElemType
QUAD8
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
EDGE4
TRI3
QUAD4
TRI6
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EDGE2
const Node * node_ptr(const unsigned int i) const
TRI7
QUAD9
EDGE3
virtual ElemType type() const=0

◆ getSmoothingFacesAndWeights()

void PenetrationThread::getSmoothingFacesAndWeights ( PenetrationInfo info,
std::vector< PenetrationInfo *> &  edge_face_info,
std::vector< Real > &  edge_face_weights,
std::vector< PenetrationInfo *> &  p_info,
const Node secondary_node 
)
protected

Definition at line 1300 of file PenetrationThread.C.

Referenced by smoothNormal().

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 }
void getSmoothingEdgeNodesAndWeights(const libMesh::Point &p, const Elem *side, std::vector< std::vector< const Node *>> &edge_nodes, std::vector< Real > &edge_face_weights)
MPI_Info info
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
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)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ interactionsOffCommonEdge()

PenetrationThread::CommonEdgeResult PenetrationThread::interactionsOffCommonEdge ( PenetrationInfo pi1,
PenetrationInfo pi2 
)
protected

Definition at line 606 of file PenetrationThread.C.

Referenced by competeInteractions().

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 }
const Elem * _side
std::vector< const Node * > _off_edge_nodes
virtual unsigned short dim() const=0

◆ isFaceReasonableCandidate()

bool PenetrationThread::isFaceReasonableCandidate ( const Elem primary_elem,
const Elem side,
libMesh::FEBase fe,
const libMesh::Point secondary_point,
const Real  tangential_tolerance 
)
protected

Definition at line 1157 of file PenetrationThread.C.

Referenced by createInfoForElem().

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 }
auto norm() const -> decltype(std::norm(Real()))
virtual_for_inffe const std::vector< RealGradient > & get_dxyzdeta() const
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
virtual Real hmax() const
const Node *const * get_nodes() const
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
virtual_for_inffe const std::vector< Point > & get_xyz() const
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const=0
virtual_for_inffe const std::vector< RealGradient > & get_dxyzdxi() const

◆ join()

void PenetrationThread::join ( const PenetrationThread other)

Definition at line 472 of file PenetrationThread.C.

473 {
475  other._recheck_secondary_nodes.begin(),
476  other._recheck_secondary_nodes.end());
477 }
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...

◆ operator()()

void PenetrationThread::operator() ( const NodeIdRange range)

Definition at line 90 of file PenetrationThread.C.

Referenced by getSidesOnPrimaryBoundary().

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;
271  rsd._distance = std::numeric_limits<Real>::max();
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 }
MooseVariable * _nodal_normal_z
auto norm() const -> decltype(std::norm(Real()))
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)
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.
NearestNodeLocator & _nearest_node
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)
Threads::spin_mutex pinfo_mutex
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
auto max(const L &left, const R &right)
SubProblem & _subproblem
const std::vector< std::vector< OutputGradient > > & get_dphi() const
const MooseMesh & _mesh
dof_id_type id() const
CompeteInteractionResult competeInteractions(PenetrationInfo *pi1, PenetrationInfo *pi2)
When interactions are identified between a node and two faces, compete between the faces to determine...
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...
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 restrictPointToFace(libMesh::Point &p, const Node *&closest_node, const Elem *side)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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...
MooseVariable * _nodal_normal_y
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
MooseVariable * _nodal_normal_x
std::vector< std::vector< libMesh::FEBase * > > & _fes
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method
const std::vector< std::vector< OutputShape > > & get_phi() const

◆ restrictPointToFace()

bool PenetrationThread::restrictPointToFace ( libMesh::Point p,
const Node *&  closest_node,
const Elem side 
)
protected

Definition at line 976 of file PenetrationThread.C.

Referenced by operator()().

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 }
ElemType
QUAD8
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
EDGE4
TRI3
QUAD4
TRI6
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EDGE2
const Node * node_ptr(const unsigned int i) const
TRI7
QUAD9
EDGE3
virtual ElemType type() const=0

◆ restrictPointToSpecifiedEdgeOfFace()

bool PenetrationThread::restrictPointToSpecifiedEdgeOfFace ( libMesh::Point p,
const Node *&  closest_node,
const Elem side,
const std::vector< const Node *> &  edge_nodes 
)
protected

Definition at line 802 of file PenetrationThread.C.

Referenced by findRidgeContactPoint().

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 }
ElemType
QUAD8
const unsigned int invalid_uint
unsigned int get_node_index(const Node *node_ptr) const
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
EDGE4
TRI3
QUAD4
TRI6
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EDGE2
virtual unsigned short dim() const=0
const Node * node_ptr(const unsigned int i) const
TRI7
QUAD9
EDGE3
virtual ElemType type() const=0

◆ smoothNormal()

void PenetrationThread::smoothNormal ( PenetrationInfo info,
std::vector< PenetrationInfo *> &  p_info,
const Node node 
)
protected

Definition at line 1242 of file PenetrationThread.C.

Referenced by operator()().

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 }
MooseVariable * _nodal_normal_z
auto norm() const -> decltype(std::norm(Real()))
RealVectorValue _normal
MPI_Info info
OutputType getValue(const Elem *elem, const std::vector< std::vector< OutputShape >> &phi) const
Compute the variable value at a point on an element.
Data structure used to hold penetration information.
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)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseVariable * _nodal_normal_y
MooseVariable * _nodal_normal_x
PenetrationLocator::NORMAL_SMOOTHING_METHOD _normal_smoothing_method

◆ switchInfo()

void PenetrationThread::switchInfo ( PenetrationInfo *&  info,
PenetrationInfo *&  infoNew 
)
protected

Definition at line 480 of file PenetrationThread.C.

Referenced by operator()().

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 }
MPI_Info info
const Elem * _starting_elem
unsigned int _stick_locked_this_step
RealVectorValue _contact_force_old
unsigned int _locked_this_step
const Elem * _elem
unsigned int _starting_side_num
unsigned int _side_num
RealVectorValue _contact_force
MECH_STATUS_ENUM _mech_status
Point _starting_closest_point_ref
MECH_STATUS_ENUM _mech_status_old
RealVectorValue _lagrange_multiplier_slip

Member Data Documentation

◆ _bc_tuples

const std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type> >& PenetrationThread::_bc_tuples
protected

Definition at line 84 of file PenetrationThread.h.

Referenced by getSidesOnPrimaryBoundary().

◆ _check_whether_reasonable

bool PenetrationThread::_check_whether_reasonable
protected

Definition at line 65 of file PenetrationThread.h.

Referenced by operator()().

◆ _do_normal_smoothing

bool PenetrationThread::_do_normal_smoothing
protected

Definition at line 68 of file PenetrationThread.h.

Referenced by operator()(), and smoothNormal().

◆ _elem_side_builder

libMesh::ElemSideBuilder PenetrationThread::_elem_side_builder
protected

Helper for building element sides without extraneous allocation.

Definition at line 89 of file PenetrationThread.h.

Referenced by computeSlip().

◆ _fe_type

libMesh::FEType& PenetrationThread::_fe_type
protected

Definition at line 77 of file PenetrationThread.h.

Referenced by createInfoForElem(), and operator()().

◆ _fes

std::vector<std::vector<libMesh::FEBase *> >& PenetrationThread::_fes
protected

Definition at line 75 of file PenetrationThread.h.

Referenced by createInfoForElem(), findRidgeContactPoint(), and operator()().

◆ _mesh

const MooseMesh& PenetrationThread::_mesh
protected

Definition at line 58 of file PenetrationThread.h.

Referenced by getInfoForFacesWithCommonNodes(), and operator()().

◆ _nearest_node

NearestNodeLocator& PenetrationThread::_nearest_node
protected

Definition at line 79 of file PenetrationThread.h.

Referenced by operator()().

◆ _nodal_normal_x

MooseVariable* PenetrationThread::_nodal_normal_x
protected

Definition at line 71 of file PenetrationThread.h.

Referenced by operator()(), and smoothNormal().

◆ _nodal_normal_y

MooseVariable* PenetrationThread::_nodal_normal_y
protected

Definition at line 72 of file PenetrationThread.h.

Referenced by operator()(), and smoothNormal().

◆ _nodal_normal_z

MooseVariable* PenetrationThread::_nodal_normal_z
protected

Definition at line 73 of file PenetrationThread.h.

Referenced by operator()(), and smoothNormal().

◆ _node_to_elem_map

const std::map<dof_id_type, std::vector<dof_id_type> >& PenetrationThread::_node_to_elem_map
protected

Definition at line 81 of file PenetrationThread.h.

Referenced by getInfoForFacesWithCommonNodes(), and operator()().

◆ _normal_smoothing_distance

Real PenetrationThread::_normal_smoothing_distance
protected

Definition at line 69 of file PenetrationThread.h.

Referenced by getSmoothingEdgeNodesAndWeights().

◆ _normal_smoothing_method

PenetrationLocator::NORMAL_SMOOTHING_METHOD PenetrationThread::_normal_smoothing_method
protected

Definition at line 70 of file PenetrationThread.h.

Referenced by operator()(), and smoothNormal().

◆ _penetration_info

std::map<dof_id_type, PenetrationInfo *>& PenetrationThread::_penetration_info
protected

Definition at line 63 of file PenetrationThread.h.

Referenced by operator()().

◆ _primary_boundary

BoundaryID PenetrationThread::_primary_boundary
protected

Definition at line 59 of file PenetrationThread.h.

Referenced by getSidesOnPrimaryBoundary().

◆ _recheck_secondary_nodes

std::vector<dof_id_type> PenetrationThread::_recheck_secondary_nodes

List of secondary nodes for which penetration was not detected in the current patch and for which patch has to be updated.

Definition at line 53 of file PenetrationThread.h.

Referenced by join(), and operator()().

◆ _secondary_boundary

BoundaryID PenetrationThread::_secondary_boundary
protected

Definition at line 60 of file PenetrationThread.h.

◆ _subproblem

SubProblem& PenetrationThread::_subproblem
protected

Definition at line 56 of file PenetrationThread.h.

Referenced by operator()().

◆ _tangential_tolerance

Real PenetrationThread::_tangential_tolerance
protected

Definition at line 67 of file PenetrationThread.h.

Referenced by competeInteractions(), createInfoForElem(), and operator()().

◆ _tid

THREAD_ID PenetrationThread::_tid
protected

Definition at line 86 of file PenetrationThread.h.

Referenced by createInfoForElem(), findRidgeContactPoint(), and operator()().

◆ _update_location

bool PenetrationThread::_update_location
protected

Definition at line 66 of file PenetrationThread.h.

Referenced by operator()().


The documentation for this class was generated from the following files: