PenetrationThread.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
// Moose
#include "PenetrationThread.h"
#include "ParallelUniqueId.h"
#include "FindContactPoint.h"
#include "NearestNodeLocator.h"
#include "SubProblem.h"
#include "MooseVariableFE.h"
#include "MooseMesh.h"
#include "MooseUtils.h"
 
#include "libmesh/threads.h"
 
#include <algorithm>
 
// Anonymous namespace for helper functions that ought to be moved
// into libMesh
namespace
{
Point
closest_point_to_edge(const Point & src, const Point & p0, const Point & p1)
{
  const Point line01 = p1 - p0;
  const Real line0c_xi = ((src - p0) * line01) / line01.norm_sq();
  // The projection would be behind p0; p0 is closest
  if (line0c_xi <= 0)
    return p0;
  // The projection would be past p1; p1 is closest
  if (line0c_xi >= 1)
    return p1;
  // The projection is on the segment between p0 to p1.
  return p0 + line0c_xi * line01;
}
 
Point
closest_point_to_side(const Point & src, const Elem & side)
{
  switch (side.type())
  {
    case EDGE2:
    case EDGE3:
    case EDGE4:
      mooseAssert(side.has_affine_map(),
                  "Penetration of elements with curved sides not implemented");
      return closest_point_to_edge(src, side.point(0), side.point(1));
    case TRI3:
    case TRI6:
    {
      mooseAssert(side.has_affine_map(),
                  "Penetration of elements with curved sides not implemented");
      const Point p0 = side.point(0), p1 = side.point(1), p2 = side.point(2);
      const Point l01 = p1 - p0, l02 = p2 - p0;
      const Point tri_normal = (l01.cross(l02)).unit();
      const Point linecs = ((src - p0) * tri_normal) / tri_normal.norm_sq() * tri_normal;
      const Point in_plane = src - linecs;
      const Point planar_offset = in_plane - p0;
      // If we're outside the triangle past line 01, our closest point
      // is on that line.
      if (planar_offset.cross(l01) * tri_normal > 0)
        return closest_point_to_edge(src, p0, p1);
      // If we're outside the triangle past line 02, our closest point
      // is on that line.
      if (planar_offset.cross(l02) * tri_normal < 0)
        return closest_point_to_edge(src, p0, p2);
      // If we're outside the triangle past line 12, our closest point
      // is on that line.
      if ((in_plane - p1).cross(p2 - p1) * tri_normal > 0)
        return closest_point_to_edge(src, p1, p2);
      // We must be inside the triangle!
      return in_plane;
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    case C0POLYGON:
      mooseError("Not implemented");
    default:
      mooseError("Side type not recognized");
      break;
  }
}
 
} // anonymous namespace
 
// Mutex to use when accessing _penetration_info;
Threads::spin_mutex pinfo_mutex;
 
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,
    bool use_point_locator,
    std::vector<std::vector<FEBase *>> & fes,
    FEType & fe_type,
    NearestNodeLocator & nearest_node,
    const std::unordered_map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map)
  : _subproblem(subproblem),
    _mesh(mesh),
    _primary_boundary(primary_boundary),
    _secondary_boundary(secondary_boundary),
    _penetration_info(penetration_info),
    _check_whether_reasonable(check_whether_reasonable),
    _update_location(update_location),
    _tangential_tolerance(tangential_tolerance),
    _do_normal_smoothing(do_normal_smoothing),
    _normal_smoothing_distance(normal_smoothing_distance),
    _normal_smoothing_method(normal_smoothing_method),
    _use_point_locator(use_point_locator),
    _nodal_normal_x(NULL),
    _nodal_normal_y(NULL),
    _nodal_normal_z(NULL),
    _fes(fes),
    _fe_type(fe_type),
    _nearest_node(nearest_node),
    _node_to_elem_map(node_to_elem_map)
{
}
 
// Splitting Constructor
PenetrationThread::PenetrationThread(PenetrationThread & x, Threads::split /*split*/)
  : _subproblem(x._subproblem),
    _mesh(x._mesh),
    _primary_boundary(x._primary_boundary),
    _secondary_boundary(x._secondary_boundary),
    _penetration_info(x._penetration_info),
    _check_whether_reasonable(x._check_whether_reasonable),
    _update_location(x._update_location),
    _tangential_tolerance(x._tangential_tolerance),
    _do_normal_smoothing(x._do_normal_smoothing),
    _normal_smoothing_distance(x._normal_smoothing_distance),
    _normal_smoothing_method(x._normal_smoothing_method),
    _use_point_locator(x._use_point_locator),
    _fes(x._fes),
    _fe_type(x._fe_type),
    _nearest_node(x._nearest_node),
    _node_to_elem_map(x._node_to_elem_map)
{
}
 
void
PenetrationThread::operator()(const NodeIdRange & range)
{
  ParallelUniqueId puid;
  _tid = puid.id;
 
  // Must get the variables every time this is run because _tid can change
  if (_do_normal_smoothing &&
      _normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
  {
    _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
    _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
    _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
  }
 
  const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
  std::unique_ptr<PointLocatorBase> point_locator;
  if (_use_point_locator)
    point_locator = _mesh.getPointLocator();
 
  for (const auto & node_id : range)
  {
    const Node & node = _mesh.nodeRef(node_id);
 
    // We're going to get a reference to the pointer for the pinfo for this node
    // This will allow us to manipulate this pointer without having to go through
    // the _penetration_info map... meaning this is the only mutex we'll have to do!
    pinfo_mutex.lock();
    PenetrationInfo *& info = _penetration_info[node.id()];
    pinfo_mutex.unlock();
 
    std::vector<PenetrationInfo *> p_info;
    bool info_set(false);
 
    // See if we already have info about this node
    if (info)
    {
      FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
      FEBase * fe_side = _fes[_tid][info->_side->dim()];
 
      if (!_update_location && (info->_distance >= 0 || info->isCaptured()))
      {
        const Point contact_ref = info->_closest_point_ref;
        bool contact_point_on_side(false);
 
        // Secondary position must be the previous contact point
        // Use the previous reference coordinates
        std::vector<Point> points(1);
        points[0] = contact_ref;
        const std::vector<Point> & secondary_pos = fe_side->get_xyz();
        bool search_succeeded = false;
 
        Moose::findContactPoint(*info,
                                fe_elem,
                                fe_side,
                                _fe_type,
                                secondary_pos[0],
                                false,
                                _tangential_tolerance,
                                contact_point_on_side,
                                search_succeeded);
 
        // Restore the original reference coordinates
        info->_closest_point_ref = contact_ref;
        // Just calculated as the distance of the contact point off the surface (0).  Set to 0 to
        // avoid round-off.
        info->_distance = 0.0;
        info_set = true;
      }
      else
      {
        Real old_tangential_distance(info->_tangential_distance);
        bool contact_point_on_side(false);
        bool search_succeeded = false;
 
        Moose::findContactPoint(*info,
                                fe_elem,
                                fe_side,
                                _fe_type,
                                node,
                                false,
                                _tangential_tolerance,
                                contact_point_on_side,
                                search_succeeded);
 
        if (contact_point_on_side)
        {
          if (info->_tangential_distance <= 0.0) // on the face
          {
            info_set = true;
          }
          else if (info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
          { // off the face but within tolerance, was that way on the last step too
            if (info->_side->dim() == 2 && info->_off_edge_nodes.size() < 2)
            { // Closest point on face is on a node rather than an edge.  Another
              // face might be a better candidate.
            }
            else
            {
              info_set = true;
            }
          }
        }
      }
    }
 
    if (!info_set)
    {
      const Node * closest_node = _nearest_node.nearestNode(node.id());
 
      std::vector<dof_id_type> located_elem_ids;
      const std::vector<dof_id_type> * closest_elems;
 
      if (_use_point_locator)
      {
        std::set<const Elem *> candidate_elements;
        (*point_locator)(*closest_node, candidate_elements);
 
        if (candidate_elements.empty())
          mooseError("No proximate elements found at node ",
                     closest_node->id(),
                     " at ",
                     cast_ref<const Point &>(*closest_node),
                     " on boundary ",
                     _nearest_node._boundary1,
                     ".  This should never happen.");
 
        for (const Elem * elem : candidate_elements)
        {
          for (auto s : elem->side_index_range())
            if (boundary_info.has_boundary_id(elem, s, _primary_boundary))
            {
              located_elem_ids.push_back(elem->id());
              break;
            }
        }
 
        if (located_elem_ids.empty())
          mooseError("No proximate elements found at node ",
                     closest_node->id(),
                     " at ",
                     cast_ref<const Point &>(*closest_node),
                     " on boundary ",
                     _nearest_node._boundary1,
                     " share that boundary.  This may happen if the mesh uses the same boundary id "
                     "for a nodeset and an unrelated sideset.");
 
        closest_elems = &located_elem_ids;
      }
      else
      {
        auto node_to_elem_pair = _node_to_elem_map.find(closest_node->id());
        mooseAssert(node_to_elem_pair != _node_to_elem_map.end(),
                    "Missing entry in node to elem map");
        closest_elems = &(node_to_elem_pair->second);
      }
 
      for (const auto & elem_id : *closest_elems)
      {
        const Elem * elem = _mesh.elemPtr(elem_id);
 
        std::vector<PenetrationInfo *> thisElemInfo;
 
        std::vector<const Node *> nodesThatMustBeOnSide;
        // If we have a disconnected mesh, we might not have *any*
        // nodes that must be on a side we check; we'll rely on
        // boundary info to find valid sides, then rely on comparing
        // closest points from each to find the best.
        //
        // If we don't have a disconnected mesh, then for maximum
        // backwards compatibility we're still using the older ridge
        // and peak detection code, which depends on us ruling out
        // sides that don't touch closest_node.
        if (!_use_point_locator)
          nodesThatMustBeOnSide.push_back(closest_node);
        createInfoForElem(
            thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
      }
 
      if (_use_point_locator)
      {
        Real min_distance_sq = std::numeric_limits<Real>::max();
        Point best_point;
        unsigned int best_i = invalid_uint;
 
        // Find closest point in all p_info to the node of interest
        for (unsigned int i = 0; i < p_info.size(); ++i)
        {
          const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
          const Real distance_sq = (closest_point - node).norm_sq();
          if (distance_sq < min_distance_sq)
          {
            min_distance_sq = distance_sq;
            best_point = closest_point;
            best_i = i;
          }
        }
 
        p_info[best_i]->_closest_point = best_point;
        p_info[best_i]->_distance =
            (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) * std::sqrt(min_distance_sq);
        if (_do_normal_smoothing)
          mooseError("Normal smoothing not implemented with point locator code");
        Point normal = (best_point - node).unit();
        const Real dot = normal * p_info[best_i]->_normal;
        if (dot < 0)
          normal *= -1;
        p_info[best_i]->_normal = normal;
 
        switchInfo(info, p_info[best_i]);
        info_set = true;
      }
      else
      {
        if (p_info.size() == 1)
        {
          if (p_info[0]->_tangential_distance <= _tangential_tolerance)
          {
            switchInfo(info, p_info[0]);
            info_set = true;
          }
        }
        else if (p_info.size() > 1)
        {
          // Loop through all pairs of faces, and check for contact on ridge between each face pair
          std::vector<RidgeData> ridgeDataVec;
          for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
            for (unsigned int j = i + 1; j < p_info.size(); ++j)
            {
              Point closest_coor;
              Real tangential_distance(0.0);
              const Node * closest_node_on_ridge = NULL;
              unsigned int index = 0;
              Point closest_coor_ref;
              bool found_ridge_contact_point = findRidgeContactPoint(closest_coor,
                                                                     tangential_distance,
                                                                     closest_node_on_ridge,
                                                                     index,
                                                                     closest_coor_ref,
                                                                     p_info,
                                                                     i,
                                                                     j);
              if (found_ridge_contact_point)
              {
                RidgeData rpd;
                rpd._closest_coor = closest_coor;
                rpd._tangential_distance = tangential_distance;
                rpd._closest_node = closest_node_on_ridge;
                rpd._index = index;
                rpd._closest_coor_ref = closest_coor_ref;
                ridgeDataVec.push_back(rpd);
              }
            }
 
          if (ridgeDataVec.size() > 0) // Either find the ridge pair that is the best or find a peak
          {
            // Group together ridges for which we are off the edge of a common node.
            // Those are peaks.
            std::vector<RidgeSetData> ridgeSetDataVec;
            for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
            {
              bool foundSetWithMatchingNode = false;
              for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
              {
                if (ridgeDataVec[i]._closest_node != NULL &&
                    ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
                {
                  foundSetWithMatchingNode = true;
                  ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
                  break;
                }
              }
              if (!foundSetWithMatchingNode)
              {
                RidgeSetData rsd;
                rsd._distance = std::numeric_limits<Real>::max();
                rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
                rsd._closest_node = ridgeDataVec[i]._closest_node;
                ridgeSetDataVec.push_back(rsd);
              }
            }
            // Compute distance to each set of ridges
            for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
            {
              if (ridgeSetDataVec[i]._closest_node !=
                  NULL) // Either a peak or off the edge of single ridge
              {
                if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
                {
                  if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
                      _tangential_tolerance) // off within tolerance
                  {
                    ridgeSetDataVec[i]._closest_coor =
                        ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
                    Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
                    ridgeSetDataVec[i]._distance = contact_point_vec.norm();
                  }
                }
                else // several ridges join at common node to make a peak.  The common node is the
                     // contact point
                {
                  ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
                  Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
                  ridgeSetDataVec[i]._distance = contact_point_vec.norm();
                }
              }
              else // on a single ridge
              {
                ridgeSetDataVec[i]._closest_coor =
                    ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
                Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
                ridgeSetDataVec[i]._distance = contact_point_vec.norm();
              }
            }
            // Find the set of ridges closest to us.
            unsigned int closest_ridge_set_index(0);
            Real closest_distance(ridgeSetDataVec[0]._distance);
            Point closest_point(ridgeSetDataVec[0]._closest_coor);
            for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
            {
              if (ridgeSetDataVec[i]._distance < closest_distance)
              {
                closest_ridge_set_index = i;
                closest_distance = ridgeSetDataVec[i]._distance;
                closest_point = ridgeSetDataVec[i]._closest_coor;
              }
            }
 
            if (closest_distance <
                std::numeric_limits<Real>::max()) // contact point is on the closest ridge set
            {
              // find the face in the ridge set with the smallest index, assign that one to the
              // interaction
              // TODO: We may need to select the face with the largest projected distance
              // rather than the smallest index, similar to what is done when picking the
              // face in findRidgeContactPoint() to better condition corner-case checks for
              // sign changes. That's less likely to be a problem here, though, and the
              // code to do that would be messier.
              unsigned int face_index(std::numeric_limits<unsigned int>::max());
              for (unsigned int i = 0;
                   i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
                   ++i)
              {
                if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
                  face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
              }
 
              mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
                          "face_index invalid");
 
              p_info[face_index]->_closest_point = closest_point;
              p_info[face_index]->_distance =
                  (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
              // Calculate the normal as the vector from the ridge to the point only if we're not
              // doing normal
              // smoothing.  Normal smoothing will average out the normals on its own.
              if (!_do_normal_smoothing)
              {
                Point normal(closest_point - node);
                const Real len(normal.norm());
                if (len > 0)
                {
                  normal /= len;
                }
                const Real dot(normal * p_info[face_index]->_normal);
                if (dot < 0)
                  normal *= -1;
                p_info[face_index]->_normal = normal;
              }
              p_info[face_index]->_tangential_distance = 0.0;
 
              Point closest_point_ref;
              if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
                  1) // contact with a single ridge rather than a peak
              {
                p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
                                                               ._ridge_data_vec[0]
                                                               ._tangential_distance;
                p_info[face_index]->_closest_point_ref =
                    ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
              }
              else
              { // peak
                const Node * closest_node_on_face;
                bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
                                                      closest_node_on_face,
                                                      p_info[face_index]->_side);
                if (restricted)
                {
                  if (closest_node_on_face !=
                      ridgeSetDataVec[closest_ridge_set_index]._closest_node)
                  {
                    mooseError("Closest node when restricting point to face != closest node from "
                               "RidgeSetData");
                  }
                }
              }
 
              FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
              std::vector<Point> points(1);
              points[0] = p_info[face_index]->_closest_point_ref;
              fe->reinit(p_info[face_index]->_side, &points);
              p_info[face_index]->_side_phi = fe->get_phi();
              p_info[face_index]->_side_grad_phi = fe->get_dphi();
              p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
              p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
              p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
 
              switchInfo(info, p_info[face_index]);
              info_set = true;
            }
            else
            { // todo:remove invalid ridge cases so they don't mess up individual face
              // competition????
            }
          }
 
          if (!info_set) // contact wasn't on a ridge -- compete individual interactions
          {
            unsigned int best(0), i(1);
            do
            {
              CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
              if (CIResult == FIRST_WINS)
              {
                i++;
              }
              else if (CIResult == SECOND_WINS)
              {
                best = i;
                i++;
              }
              else if (CIResult == NEITHER_WINS)
              {
                best = i + 1;
                i += 2;
              }
            } while (i < p_info.size() && best < p_info.size());
            if (best < p_info.size())
            {
              // Ensure final info is within the tangential tolerance
              if (p_info[best]->_tangential_distance <= _tangential_tolerance)
              {
                switchInfo(info, p_info[best]);
                info_set = true;
              }
            }
          }
        }
      }
    }
 
    if (!info_set)
    {
      // If penetration is not detected within the saved patch, it is possible
      // that the secondary node has moved outside the saved patch. So, the patch
      // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
      // and the penetration detection has to be re-run on the updated patch.
 
      _recheck_secondary_nodes.push_back(node_id);
 
      delete info;
      info = NULL;
    }
    else
    {
      smoothNormal(info, p_info, node);
      FEBase * fe = _fes[_tid][info->_side->dim()];
      computeSlip(*fe, *info);
    }
 
    for (unsigned int j = 0; j < p_info.size(); ++j)
    {
      if (p_info[j])
      {
        delete p_info[j];
        p_info[j] = NULL;
      }
    }
  }
}
 
void
PenetrationThread::join(const PenetrationThread & other)
{
  _recheck_secondary_nodes.insert(_recheck_secondary_nodes.end(),
                                  other._recheck_secondary_nodes.begin(),
                                  other._recheck_secondary_nodes.end());
}
 
void
PenetrationThread::switchInfo(PenetrationInfo *& info, PenetrationInfo *& infoNew)
{
  mooseAssert(infoNew != NULL, "infoNew object is null");
  if (info)
  {
    infoNew->_starting_elem = info->_starting_elem;
    infoNew->_starting_side_num = info->_starting_side_num;
    infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
    infoNew->_incremental_slip = info->_incremental_slip;
    infoNew->_accumulated_slip = info->_accumulated_slip;
    infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
    infoNew->_frictional_energy = info->_frictional_energy;
    infoNew->_frictional_energy_old = info->_frictional_energy_old;
    infoNew->_contact_force = info->_contact_force;
    infoNew->_contact_force_old = info->_contact_force_old;
    infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
    infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
    infoNew->_locked_this_step = info->_locked_this_step;
    infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
    infoNew->_mech_status = info->_mech_status;
    infoNew->_mech_status_old = info->_mech_status_old;
  }
  else
  {
    infoNew->_starting_elem = infoNew->_elem;
    infoNew->_starting_side_num = infoNew->_side_num;
    infoNew->_starting_closest_point_ref = infoNew->_closest_point_ref;
  }
  delete info;
  info = infoNew;
  infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
}
 
PenetrationThread::CompeteInteractionResult
PenetrationThread::competeInteractions(PenetrationInfo * pi1, PenetrationInfo * pi2)
{
 
  CompeteInteractionResult result = NEITHER_WINS;
 
  if (pi1->_tangential_distance > _tangential_tolerance &&
      pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
    result = NEITHER_WINS;
 
  else if (pi1->_tangential_distance == 0.0 &&
           pi2->_tangential_distance > 0.0) // on face 1, off face 2
    result = FIRST_WINS;
 
  else if (pi2->_tangential_distance == 0.0 &&
           pi1->_tangential_distance > 0.0) // on face 2, off face 1
    result = SECOND_WINS;
 
  else if (pi1->_tangential_distance <= _tangential_tolerance &&
           pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
    result = FIRST_WINS;
 
  else if (pi2->_tangential_distance <= _tangential_tolerance &&
           pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
    result = SECOND_WINS;
 
  else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
    result = competeInteractionsBothOnFace(pi1, pi2);
 
  else if (pi1->_tangential_distance <= _tangential_tolerance &&
           pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
  {
    CommonEdgeResult cer = interactionsOffCommonEdge(pi1, pi2);
    if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
    {
      // We already checked for ridges, and it got rejected, so neither face must be valid
      result = NEITHER_WINS;
      //      mooseError("Erroneously encountered ridge case");
    }
    else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face.  Favor
                                          // the off-side face
    {
      if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
        mooseError("Invalid off_edge_nodes counts");
 
      else if (pi1->_off_edge_nodes.size() == 2)
        result = FIRST_WINS;
 
      else if (pi2->_off_edge_nodes.size() == 2)
        result = SECOND_WINS;
 
      else
        mooseError("Invalid off_edge_nodes counts");
    }
    else // The node projects to both faces within tangential tolerance.
      result = competeInteractionsBothOnFace(pi1, pi2);
  }
 
  return result;
}
 
PenetrationThread::CompeteInteractionResult
PenetrationThread::competeInteractionsBothOnFace(PenetrationInfo * pi1, PenetrationInfo * pi2)
{
  CompeteInteractionResult result = NEITHER_WINS;
 
  if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
    result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
 
  else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
    result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
 
  // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
  //       movement.  If there is some way to smooth the transition, we should do it.
  else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi1->_distance), std::abs(pi2->_distance)))
    result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
 
  else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi2->_distance), std::abs(pi1->_distance)))
    result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
 
  else // Equal within tolerance.  Favor the one with a smaller element id (for repeatibility)
  {
    if (pi1->_elem->id() < pi2->_elem->id())
      result = FIRST_WINS;
 
    else
      result = SECOND_WINS;
  }
 
  return result;
}
 
PenetrationThread::CommonEdgeResult
PenetrationThread::interactionsOffCommonEdge(PenetrationInfo * pi1, PenetrationInfo * pi2)
{
  CommonEdgeResult common_edge(NO_COMMON);
  const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
  const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
  const unsigned dim1 = pi1->_side->dim();
 
  if (dim1 == 1)
  {
    mooseAssert(pi2->_side->dim() == 1, "Incompatible dimensions.");
    mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
                "off_edge_nodes size should be <2 for 2D contact");
    if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
        off_edge_nodes1[0] == off_edge_nodes2[0])
      common_edge = COMMON_EDGE;
  }
  else
  {
    mooseAssert(dim1 == 2 && pi2->_side->dim() == 2, "Incompatible dimensions.");
    mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
                "off_edge_nodes size should be <3 for 3D contact");
    if (off_edge_nodes1.size() == 1)
    {
      if (off_edge_nodes2.size() == 1)
      {
        if (off_edge_nodes1[0] == off_edge_nodes2[0])
          common_edge = COMMON_NODE;
      }
      else if (off_edge_nodes2.size() == 2)
      {
        if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
          common_edge = EDGE_AND_COMMON_NODE;
      }
    }
    else if (off_edge_nodes1.size() == 2)
    {
      if (off_edge_nodes2.size() == 1)
      {
        if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
          common_edge = EDGE_AND_COMMON_NODE;
      }
      else if (off_edge_nodes2.size() == 2)
      {
        if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
             off_edge_nodes1[1] == off_edge_nodes2[1]) ||
            (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
          common_edge = COMMON_EDGE;
      }
    }
  }
  return common_edge;
}
 
bool
PenetrationThread::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)
{
  tangential_distance = 0.0;
  closest_node = NULL;
  PenetrationInfo * pi1 = p_info[index1];
  PenetrationInfo * pi2 = p_info[index2];
  const unsigned sidedim(pi1->_side->dim());
  mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
 
  // Nodes on faces for the two interactions
  std::vector<const Node *> side1_nodes;
  getSideCornerNodes(pi1->_side, side1_nodes);
  std::vector<const Node *> side2_nodes;
  getSideCornerNodes(pi2->_side, side2_nodes);
 
  std::sort(side1_nodes.begin(), side1_nodes.end());
  std::sort(side2_nodes.begin(), side2_nodes.end());
 
  // Find nodes shared by the two faces
  std::vector<const Node *> common_nodes;
  std::set_intersection(side1_nodes.begin(),
                        side1_nodes.end(),
                        side2_nodes.begin(),
                        side2_nodes.end(),
                        std::inserter(common_nodes, common_nodes.end()));
 
  if (common_nodes.size() != sidedim)
    return false;
 
  bool found_point1, found_point2;
  Point closest_coor_ref1(pi1->_closest_point_ref);
  const Node * closest_node1;
  found_point1 = restrictPointToSpecifiedEdgeOfFace(
      closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
 
  Point closest_coor_ref2(pi2->_closest_point_ref);
  const Node * closest_node2;
  found_point2 = restrictPointToSpecifiedEdgeOfFace(
      closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
 
  if (!found_point1 || !found_point2)
    return false;
 
  //  if (sidedim == 2)
  //  {
  // TODO:
  // We have the parametric coordinates of the closest intersection point for both faces.
  // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
  // Find that point by taking dot products of vector from contact point to secondary node point
  // with face normal vectors to see which face we're closer to.
  //  }
 
  FEBase * fe = NULL;
  std::vector<Point> points(1);
 
  // We have to pick one of the two faces to own the contact point.  Either one would
  // generally work, but to avoid some corner-case numerical roundoff issues when
  // determining signs of the normal and displacement, pick the one that the point is
  // closer to projecting onto, which is the one with the larger projected distance.
  // If that distance is the same for both faces (within numerical precision), pick the
  // face with the lowest index for repeatability.
  if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(pi1->_distance), std::abs(pi2->_distance)) ||
      (MooseUtils::absoluteFuzzyEqual(std::abs(pi1->_distance), std::abs(pi2->_distance)) &&
       index1 < index2))
  {
    fe = _fes[_tid][pi1->_side->dim()];
    contact_point_ref = closest_coor_ref1;
    points[0] = closest_coor_ref1;
    fe->reinit(pi1->_side, &points);
    index = index1;
  }
  else
  {
    fe = _fes[_tid][pi2->_side->dim()];
    contact_point_ref = closest_coor_ref2;
    points[0] = closest_coor_ref2;
    fe->reinit(pi2->_side, &points);
    index = index2;
  }
 
  contact_point = fe->get_xyz()[0];
 
  if (sidedim == 2)
  {
    if (closest_node1) // point is off the ridge between the two elements
    {
      mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
                  "If off edge of ridge, closest node must be the same on both elements");
      closest_node = closest_node1;
 
      RealGradient off_face = *closest_node1 - contact_point;
      tangential_distance = off_face.norm();
    }
  }
 
  return true;
}
 
void
PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
{
  const ElemType t(side->type());
  corner_nodes.clear();
 
  corner_nodes.push_back(side->node_ptr(0));
  corner_nodes.push_back(side->node_ptr(1));
  switch (t)
  {
    case EDGE2:
    case EDGE3:
    case EDGE4:
    {
      break;
    }
 
    case TRI3:
    case TRI6:
    case TRI7:
    {
      corner_nodes.push_back(side->node_ptr(2));
      break;
    }
 
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      corner_nodes.push_back(side->node_ptr(2));
      corner_nodes.push_back(side->node_ptr(3));
      break;
    }
 
    default:
    {
      mooseError("Unsupported face type: ", t);
      break;
    }
  }
}
 
bool
PenetrationThread::restrictPointToSpecifiedEdgeOfFace(Point & p,
                                                      const Node *& closest_node,
                                                      const Elem * side,
                                                      const std::vector<const Node *> & edge_nodes)
{
  const ElemType t = side->type();
  Real & xi = p(0);
  Real & eta = p(1);
  closest_node = NULL;
 
  std::vector<unsigned int> local_node_indices;
  for (const auto & edge_node : edge_nodes)
  {
    unsigned int local_index = side->get_node_index(edge_node);
    if (local_index == libMesh::invalid_uint)
      mooseError("Side does not contain node");
    local_node_indices.push_back(local_index);
  }
  mooseAssert(local_node_indices.size() == side->dim(),
              "Number of edge nodes must match side dimensionality");
  std::sort(local_node_indices.begin(), local_node_indices.end());
 
  bool off_of_this_edge = false;
 
  switch (t)
  {
    case EDGE2:
    case EDGE3:
    case EDGE4:
    {
      if (local_node_indices[0] == 0)
      {
        if (xi <= -1.0)
        {
          xi = -1.0;
          off_of_this_edge = true;
          closest_node = side->node_ptr(0);
        }
      }
      else if (local_node_indices[0] == 1)
      {
        if (xi >= 1.0)
        {
          xi = 1.0;
          off_of_this_edge = true;
          closest_node = side->node_ptr(1);
        }
      }
      else
      {
        mooseError("Invalid local node indices");
      }
      break;
    }
 
    case TRI3:
    case TRI6:
    case TRI7:
    {
      if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
      {
        if (eta <= 0.0)
        {
          eta = 0.0;
          off_of_this_edge = true;
          if (xi < 0.0)
            closest_node = side->node_ptr(0);
          else if (xi > 1.0)
            closest_node = side->node_ptr(1);
        }
      }
      else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
      {
        if ((xi + eta) > 1.0)
        {
          Real delta = (xi + eta - 1.0) / 2.0;
          xi -= delta;
          eta -= delta;
          off_of_this_edge = true;
          if (xi > 1.0)
            closest_node = side->node_ptr(1);
          else if (xi < 0.0)
            closest_node = side->node_ptr(2);
        }
      }
      else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
      {
        if (xi <= 0.0)
        {
          xi = 0.0;
          off_of_this_edge = true;
          if (eta > 1.0)
            closest_node = side->node_ptr(2);
          else if (eta < 0.0)
            closest_node = side->node_ptr(0);
        }
      }
      else
      {
        mooseError("Invalid local node indices");
      }
 
      break;
    }
 
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
      {
        if (eta <= -1.0)
        {
          eta = -1.0;
          off_of_this_edge = true;
          if (xi < -1.0)
            closest_node = side->node_ptr(0);
          else if (xi > 1.0)
            closest_node = side->node_ptr(1);
        }
      }
      else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
      {
        if (xi >= 1.0)
        {
          xi = 1.0;
          off_of_this_edge = true;
          if (eta < -1.0)
            closest_node = side->node_ptr(1);
          else if (eta > 1.0)
            closest_node = side->node_ptr(2);
        }
      }
      else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
      {
        if (eta >= 1.0)
        {
          eta = 1.0;
          off_of_this_edge = true;
          if (xi < -1.0)
            closest_node = side->node_ptr(3);
          else if (xi > 1.0)
            closest_node = side->node_ptr(2);
        }
      }
      else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
      {
        if (xi <= -1.0)
        {
          xi = -1.0;
          off_of_this_edge = true;
          if (eta < -1.0)
            closest_node = side->node_ptr(0);
          else if (eta > 1.0)
            closest_node = side->node_ptr(3);
        }
      }
      else
      {
        mooseError("Invalid local node indices");
      }
      break;
    }
 
    default:
    {
      mooseError("Unsupported face type: ", t);
      break;
    }
  }
  return off_of_this_edge;
}
 
bool
PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
{
  const ElemType t(side->type());
  Real & xi = p(0);
  Real & eta = p(1);
  closest_node = NULL;
 
  bool off_of_this_face(false);
 
  switch (t)
  {
    case EDGE2:
    case EDGE3:
    case EDGE4:
    {
      if (xi < -1.0)
      {
        xi = -1.0;
        off_of_this_face = true;
        closest_node = side->node_ptr(0);
      }
      else if (xi > 1.0)
      {
        xi = 1.0;
        off_of_this_face = true;
        closest_node = side->node_ptr(1);
      }
      break;
    }
 
    case TRI3:
    case TRI6:
    case TRI7:
    {
      if (eta < 0.0)
      {
        eta = 0.0;
        off_of_this_face = true;
        if (xi < 0.5)
        {
          closest_node = side->node_ptr(0);
          if (xi < 0.0)
            xi = 0.0;
        }
        else
        {
          closest_node = side->node_ptr(1);
          if (xi > 1.0)
            xi = 1.0;
        }
      }
      else if ((xi + eta) > 1.0)
      {
        Real delta = (xi + eta - 1.0) / 2.0;
        xi -= delta;
        eta -= delta;
        off_of_this_face = true;
        if (xi > 0.5)
        {
          closest_node = side->node_ptr(1);
          if (xi > 1.0)
          {
            xi = 1.0;
            eta = 0.0;
          }
        }
        else
        {
          closest_node = side->node_ptr(2);
          if (xi < 0.0)
          {
            xi = 0.0;
            eta = 1.0;
          }
        }
      }
      else if (xi < 0.0)
      {
        xi = 0.0;
        off_of_this_face = true;
        if (eta > 0.5)
        {
          closest_node = side->node_ptr(2);
          if (eta > 1.0)
            eta = 1.0;
        }
        else
        {
          closest_node = side->node_ptr(0);
          if (eta < 0.0)
            eta = 0.0;
        }
      }
      break;
    }
 
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      if (eta < -1.0)
      {
        eta = -1.0;
        off_of_this_face = true;
        if (xi < 0.0)
        {
          closest_node = side->node_ptr(0);
          if (xi < -1.0)
            xi = -1.0;
        }
        else
        {
          closest_node = side->node_ptr(1);
          if (xi > 1.0)
            xi = 1.0;
        }
      }
      else if (xi > 1.0)
      {
        xi = 1.0;
        off_of_this_face = true;
        if (eta < 0.0)
        {
          closest_node = side->node_ptr(1);
          if (eta < -1.0)
            eta = -1.0;
        }
        else
        {
          closest_node = side->node_ptr(2);
          if (eta > 1.0)
            eta = 1.0;
        }
      }
      else if (eta > 1.0)
      {
        eta = 1.0;
        off_of_this_face = true;
        if (xi < 0.0)
        {
          closest_node = side->node_ptr(3);
          if (xi < -1.0)
            xi = -1.0;
        }
        else
        {
          closest_node = side->node_ptr(2);
          if (xi > 1.0)
            xi = 1.0;
        }
      }
      else if (xi < -1.0)
      {
        xi = -1.0;
        off_of_this_face = true;
        if (eta < 0.0)
        {
          closest_node = side->node_ptr(0);
          if (eta < -1.0)
            eta = -1.0;
        }
        else
        {
          closest_node = side->node_ptr(3);
          if (eta > 1.0)
            eta = 1.0;
        }
      }
      break;
    }
 
    default:
    {
      mooseError("Unsupported face type: ", t);
      break;
    }
  }
  return off_of_this_face;
}
 
bool
PenetrationThread::isFaceReasonableCandidate(const Elem * primary_elem,
                                             const Elem * side,
                                             FEBase * fe,
                                             const Point * secondary_point,
                                             const Real tangential_tolerance)
{
  unsigned int dim = primary_elem->dim();
 
  const std::vector<Point> & phys_point = fe->get_xyz();
 
  const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
  const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
 
  Point ref_point;
 
  std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
 
  fe->reinit(side, &points);
 
  RealGradient d = *secondary_point - phys_point[0];
 
  const Real twosqrt2 = 2.8284; // way more precision than we actually need here
  Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
 
  RealVectorValue normal;
  if (dim - 1 == 2)
  {
    normal = dxyz_dxi[0].cross(dxyz_deta[0]);
  }
  else if (dim - 1 == 1)
  {
    const Node * const * elem_nodes = primary_elem->get_nodes();
    const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
    const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
 
    Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
    out_of_plane_normal /= out_of_plane_normal.norm();
 
    normal = dxyz_dxi[0].cross(out_of_plane_normal);
  }
  else
  {
    return true;
  }
  normal /= normal.norm();
 
  const Real dot(d * normal);
 
  const RealGradient normcomp = dot * normal;
  const RealGradient tangcomp = d - normcomp;
 
  const Real tangdist = tangcomp.norm();
 
  // Increase the size of the zone that we consider if the vector from the face
  // to the node has a larger normal component
  const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
 
  bool isReasonableCandidate = true;
  if (tangdist > faceExpansionFactor * max_face_length)
  {
    isReasonableCandidate = false;
  }
  return isReasonableCandidate;
}
 
void
PenetrationThread::computeSlip(FEBase & fe, PenetrationInfo & info)
{
  // Slip is current projected position of secondary node minus
  //   original projected position of secondary node
  std::vector<Point> points(1);
  points[0] = info._starting_closest_point_ref;
  const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
  fe.reinit(&side, &points);
  const std::vector<Point> & starting_point = fe.get_xyz();
  info._incremental_slip = info._closest_point - starting_point[0];
  if (info.isCaptured())
  {
    info._frictional_energy =
        info._frictional_energy_old + info._contact_force * info._incremental_slip;
    info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
  }
}
 
void
PenetrationThread::smoothNormal(PenetrationInfo * info,
                                std::vector<PenetrationInfo *> & p_info,
                                const Node & node)
{
  if (_do_normal_smoothing)
  {
    if (_normal_smoothing_method == PenetrationLocator::NSM_EDGE_BASED)
    {
      // If we are within the smoothing distance of any edges or corners, find the
      // corner nodes for those edges/corners, and weights from distance to edge/corner
      std::vector<Real> edge_face_weights;
      std::vector<PenetrationInfo *> edge_face_info;
 
      getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
 
      mooseAssert(edge_face_info.size() == edge_face_weights.size(),
                  "edge_face_info.size() != edge_face_weights.size()");
 
      if (edge_face_info.size() > 0)
      {
        // Smooth the normal using the weighting functions for all participating faces.
        RealVectorValue new_normal;
        Real this_face_weight = 1.0;
 
        for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
        {
          PenetrationInfo * npi = edge_face_info[efwi];
          if (npi)
            new_normal += npi->_normal * edge_face_weights[efwi];
 
          this_face_weight -= edge_face_weights[efwi];
        }
        mooseAssert(this_face_weight >= (0.25 - 1e-8),
                    "Sum of weights of other faces shouldn't exceed 0.75");
        new_normal += info->_normal * this_face_weight;
 
        const Real len = new_normal.norm();
        if (len > 0)
          new_normal /= len;
 
        info->_normal = new_normal;
      }
    }
    else if (_normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
    {
      // params.addParam<VariableName>("var_name","description");
      // getParam<VariableName>("var_name")
      info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
      info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
      info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
      const Real len(info->_normal.norm());
      if (len > 0)
        info->_normal /= len;
    }
  }
}
 
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)
{
  const Elem * side = info->_side;
  const Point & p = info->_closest_point_ref;
  std::set<dof_id_type> elems_to_exclude;
  elems_to_exclude.insert(info->_elem->id());
 
  std::vector<std::vector<const Node *>> edge_nodes;
 
  // Get the pairs of nodes along every edge that we are close enough to smooth with
  getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
  std::vector<Elem *> edge_neighbor_elems;
  edge_face_info.resize(edge_nodes.size(), NULL);
 
  std::vector<unsigned int> edges_without_neighbors;
 
  for (unsigned int i = 0; i < edge_nodes.size(); ++i)
  {
    // Sort all sets of edge nodes (needed for comparing edges)
    std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
 
    std::vector<PenetrationInfo *> face_info_comm_edge;
    getInfoForFacesWithCommonNodes(
        &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
 
    if (face_info_comm_edge.size() == 0)
      edges_without_neighbors.push_back(i);
    else if (face_info_comm_edge.size() > 1)
      mooseError("Only one neighbor allowed per edge");
    else
      edge_face_info[i] = face_info_comm_edge[0];
  }
 
  // Remove edges without neighbors from the vector, starting from end
  std::vector<unsigned int>::reverse_iterator rit;
  for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
  {
    unsigned int index = *rit;
    edge_nodes.erase(edge_nodes.begin() + index);
    edge_face_weights.erase(edge_face_weights.begin() + index);
    edge_face_info.erase(edge_face_info.begin() + index);
  }
 
  // Handle corner case
  if (edge_nodes.size() > 1)
  {
    if (edge_nodes.size() != 2)
      mooseError("Invalid number of smoothing edges");
 
    // find common node
    std::vector<const Node *> common_nodes;
    std::set_intersection(edge_nodes[0].begin(),
                          edge_nodes[0].end(),
                          edge_nodes[1].begin(),
                          edge_nodes[1].end(),
                          std::inserter(common_nodes, common_nodes.end()));
 
    if (common_nodes.size() != 1)
      mooseError("Invalid number of common nodes");
 
    for (const auto & pinfo : edge_face_info)
      elems_to_exclude.insert(pinfo->_elem->id());
 
    std::vector<PenetrationInfo *> face_info_comm_edge;
    getInfoForFacesWithCommonNodes(
        &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
 
    unsigned int num_corner_neighbors = face_info_comm_edge.size();
 
    if (num_corner_neighbors > 0)
    {
      Real fw0 = edge_face_weights[0];
      Real fw1 = edge_face_weights[1];
 
      // Corner weight is product of edge weights.  Spread out over multiple neighbors.
      Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
 
      // Adjust original edge weights
      edge_face_weights[0] *= (1.0 - fw1);
      edge_face_weights[1] *= (1.0 - fw0);
 
      for (unsigned int i = 0; i < num_corner_neighbors; ++i)
      {
        edge_face_weights.push_back(fw_corner);
        edge_face_info.push_back(face_info_comm_edge[i]);
      }
    }
  }
}
 
void
PenetrationThread::getSmoothingEdgeNodesAndWeights(
    const Point & p,
    const Elem * side,
    std::vector<std::vector<const Node *>> & edge_nodes,
    std::vector<Real> & edge_face_weights)
{
  const ElemType t(side->type());
  const Real & xi = p(0);
  const Real & eta = p(1);
 
  Real smooth_limit = 1.0 - _normal_smoothing_distance;
 
  switch (t)
  {
    case EDGE2:
    case EDGE3:
    case EDGE4:
    {
      if (xi < -smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(0));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      else if (xi > smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(1));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      break;
    }
 
    case TRI3:
    case TRI6:
    case TRI7:
    {
      if (eta < -smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(0));
        en.push_back(side->node_ptr(1));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      if ((xi + eta) > smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(1));
        en.push_back(side->node_ptr(2));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      if (xi < -smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(2));
        en.push_back(side->node_ptr(0));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      break;
    }
 
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      if (eta < -smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(0));
        en.push_back(side->node_ptr(1));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      if (xi > smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(1));
        en.push_back(side->node_ptr(2));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      if (eta > smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(2));
        en.push_back(side->node_ptr(3));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      if (xi < -smooth_limit)
      {
        std::vector<const Node *> en;
        en.push_back(side->node_ptr(3));
        en.push_back(side->node_ptr(0));
        edge_nodes.push_back(en);
        Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
        if (fw > 0.5)
          fw = 0.5;
        edge_face_weights.push_back(fw);
      }
      break;
    }
 
    default:
    {
      mooseError("Unsupported face type: ", t);
      break;
    }
  }
}
 
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)
{
  // elems connected to a node on this edge, find one that has the same corners as this, and is not
  // the current elem
  auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
                                                                        // nodes
  mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
  const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
 
  std::vector<const Elem *> elems_connected_to_edge;
 
  for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
  {
    if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
      continue;
    const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
 
    std::vector<const Node *> nodevec;
    for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
      if (elem->is_vertex(ni))
        nodevec.push_back(elem->node_ptr(ni));
 
    std::vector<const Node *> common_nodes;
    std::sort(nodevec.begin(), nodevec.end());
    std::set_intersection(edge_nodes.begin(),
                          edge_nodes.end(),
                          nodevec.begin(),
                          nodevec.end(),
                          std::inserter(common_nodes, common_nodes.end()));
 
    if (common_nodes.size() == edge_nodes.size())
      elems_connected_to_edge.push_back(elem);
  }
 
  if (elems_connected_to_edge.size() > 0)
  {
 
    // There are potentially multiple elements that share a common edge
    // 2D:
    // There can only be one element on the same surface
    // 3D:
    // If there are two edge nodes, there can only be one element on the same surface
    // If there is only one edge node (a corner), there could be multiple elements on the same
    // surface
    bool allowMultipleNeighbors = false;
 
    if (elems_connected_to_edge[0]->dim() == 3)
    {
      if (edge_nodes.size() == 1)
      {
        allowMultipleNeighbors = true;
      }
    }
 
    for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
    {
      std::vector<PenetrationInfo *> thisElemInfo;
      getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
      if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
      {
        if (thisElemInfo.size() > 1)
          mooseError(
              "Found multiple neighbors to current edge/face on surface when only one is allowed");
        face_info_comm_edge.push_back(thisElemInfo[0]);
        break;
      }
 
      createInfoForElem(
          thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
      if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
      {
        if (thisElemInfo.size() > 1)
          mooseError(
              "Found multiple neighbors to current edge/face on surface when only one is allowed");
        face_info_comm_edge.push_back(thisElemInfo[0]);
        break;
      }
 
      for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
        face_info_comm_edge.push_back(thisElemInfo[j]);
    }
  }
}
 
void
PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
                                  std::vector<PenetrationInfo *> & p_info,
                                  const Elem * elem)
{
  for (const auto & pi : p_info)
  {
    if (!pi)
      continue;
 
    if (pi->_elem == elem)
      thisElemInfo.push_back(pi);
  }
}
 
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)
{
  const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
 
  for (auto s : elem->side_index_range())
  {
    if (!boundary_info.has_boundary_id(elem, s, _primary_boundary))
      continue;
 
    // Don't create info for this side if one already exists
    bool already_have_info_this_side = false;
    for (const auto & pi : thisElemInfo)
      if (pi->_side_num == s)
      {
        already_have_info_this_side = true;
        break;
      }
 
    if (already_have_info_this_side)
      break;
 
    const Elem * side = elem->build_side_ptr(s).release();
 
    // Only continue with creating info for this side if the side contains
    // all of the nodes in nodes_that_must_be_on_side
    std::vector<const Node *> nodevec;
    for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
      nodevec.push_back(side->node_ptr(ni));
 
    std::sort(nodevec.begin(), nodevec.end());
    std::vector<const Node *> common_nodes;
    std::set_intersection(nodes_that_must_be_on_side.begin(),
                          nodes_that_must_be_on_side.end(),
                          nodevec.begin(),
                          nodevec.end(),
                          std::inserter(common_nodes, common_nodes.end()));
    if (common_nodes.size() != nodes_that_must_be_on_side.size())
    {
      delete side;
      break;
    }
 
    FEBase * fe_elem = _fes[_tid][elem->dim()];
    FEBase * fe_side = _fes[_tid][side->dim()];
 
    // Optionally check to see whether face is reasonable candidate based on an
    // estimate of how closely it is likely to project to the face
    if (check_whether_reasonable)
      if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
      {
        delete side;
        break;
      }
 
    Point contact_phys;
    Point contact_ref;
    Point contact_on_face_ref;
    Real distance = 0.;
    Real tangential_distance = 0.;
    RealGradient normal;
    bool contact_point_on_side;
    std::vector<const Node *> off_edge_nodes;
    std::vector<std::vector<Real>> side_phi;
    std::vector<std::vector<RealGradient>> side_grad_phi;
    std::vector<RealGradient> dxyzdxi;
    std::vector<RealGradient> dxyzdeta;
    std::vector<RealGradient> d2xyzdxideta;
 
    std::unique_ptr<PenetrationInfo> pen_info =
        std::make_unique<PenetrationInfo>(elem,
                                          side,
                                          s,
                                          normal,
                                          distance,
                                          tangential_distance,
                                          contact_phys,
                                          contact_ref,
                                          contact_on_face_ref,
                                          off_edge_nodes,
                                          side_phi,
                                          side_grad_phi,
                                          dxyzdxi,
                                          dxyzdeta,
                                          d2xyzdxideta);
 
    bool search_succeeded = false;
    Moose::findContactPoint(*pen_info,
                            fe_elem,
                            fe_side,
                            _fe_type,
                            *secondary_node,
                            true,
                            _tangential_tolerance,
                            contact_point_on_side,
                            search_succeeded);
 
    // Do not add contact info from failed searches
    if (search_succeeded)
    {
      thisElemInfo.push_back(pen_info.get());
      p_info.push_back(pen_info.release());
    }
  }
}