Add a line search in newton solve to find contact point

framework/include/geomsearch/FindContactPoint.h

@@ -34,7 +34,8 @@ void findContactPoint(PenetrationInfo & p_info,
                       const Point & secondary_point,
                       bool start_with_centroid,
                       const Real tangential_tolerance,
-                      bool & contact_point_on_side);
+                      bool & contact_point_on_side,
+                      bool & search_succeeded);
 
 void restrictPointToFace(Point & p, const Elem * side, std::vector<const Node *> & off_edge_nodes);
 }

framework/src/constraints/MortarConstraintBase.C

@@ -51,7 +51,7 @@ MortarConstraintBase::validParams()
   // Neither is guaranteed to be a superset of the other. For instance ghosting of lower-d point
   // neighbors (AugmentSparsityOnInterface with ghost_point_neighbors = true) is only guaranteed to
   // ghost those lower-d point neighbors on *processes that own lower-d elements*. And you may have
-  // a process that only owns higher-dimensionsional elements
+  // a process that only owns higher-dimensional elements
   //
   // Note that in my experience it is only important for the higher-d lower-d point neighbors to be
   // ghosted when forming sparsity patterns and so I'm putting this here instead of at the

framework/src/geomsearch/FindContactPoint.C

@@ -40,6 +40,8 @@ namespace Moose
  * @param secondary_point The physical space coordinates of the secondary node
  * @param contact_point_on_side whether or not the contact_point actually lies on _that_ side of the
  * element.
+ * @param search_succeeded whether or not the search for the contact point succeeded. If not it
+ * should likely be discarded
  */
 void
 findContactPoint(PenetrationInfo & p_info,
@@ -49,8 +51,12 @@ findContactPoint(PenetrationInfo & p_info,
                  const Point & secondary_point,
                  bool start_with_centroid,
                  const Real tangential_tolerance,
-                 bool & contact_point_on_side)
+                 bool & contact_point_on_side,
+                 bool & search_succeeded)
 {
+  // Default to true and we'll switch on failures
+  search_succeeded = true;
+
   const Elem * primary_elem = p_info._elem;
 
   unsigned int dim = primary_elem->dim();
@@ -116,7 +122,6 @@ findContactPoint(PenetrationInfo & p_info,
   for (unsigned int it = 0; it < 3 && update_size > TOLERANCE * 1e3; ++it)
   {
     DenseMatrix<Real> jac(dim - 1, dim - 1);
-
     jac(0, 0) = -(dxyz_dxi[0] * dxyz_dxi[0]);
 
     if (dim - 1 == 2)
@@ -127,14 +132,12 @@ findContactPoint(PenetrationInfo & p_info,
     }
 
     DenseVector<Real> rhs(dim - 1);
-
     rhs(0) = dxyz_dxi[0] * d;
 
     if (dim - 1 == 2)
       rhs(1) = dxyz_deta[0] * d;
 
     DenseVector<Real> update(dim - 1);
-
     jac.lu_solve(rhs, update);
 
     ref_point(0) -= update(0);
@@ -154,12 +157,13 @@ findContactPoint(PenetrationInfo & p_info,
   unsigned nit = 0;
 
   // Newton Loop
-  for (; nit < 12 && update_size > TOLERANCE * TOLERANCE; nit++)
+  const auto max_newton_its = 25;
+  const auto tolerance_newton = 1e3 * TOLERANCE * TOLERANCE;
+  for (; nit < max_newton_its && update_size > tolerance_newton; nit++)
   {
     d = secondary_point - phys_point[0];
 
     DenseMatrix<Real> jac(dim - 1, dim - 1);
-
     jac(0, 0) = (d2xyz_dxi2[0] * d) - (dxyz_dxi[0] * dxyz_dxi[0]);
 
     if (dim - 1 == 2)
@@ -171,32 +175,75 @@ findContactPoint(PenetrationInfo & p_info,
     }
 
     DenseVector<Real> rhs(dim - 1);
-
     rhs(0) = -dxyz_dxi[0] * d;
 
     if (dim - 1 == 2)
       rhs(1) = -dxyz_deta[0] * d;
 
     DenseVector<Real> update(dim - 1);
-
     jac.lu_solve(rhs, update);
 
-    ref_point(0) += update(0);
+    // Improvised line search in case the update is too large and gets out of the element so bad
+    // that we cannot reinit at the new point
+    Real mult = 1;
+    while (true)
+    {
+      try
+      {
+        ref_point(0) += mult * update(0);
 
-    if (dim - 1 == 2)
-      ref_point(1) += update(1);
+        if (dim - 1 == 2)
+          ref_point(1) += mult * update(1);
 
-    points[0] = ref_point;
-    fe_side->reinit(side, &points);
-    d = secondary_point - phys_point[0];
+        points[0] = ref_point;
+        fe_side->reinit(side, &points);
+        d = secondary_point - phys_point[0];
 
-    update_size = update.l2_norm();
+        // we don't multiply by 'mult' because it is used for convergence
+        update_size = update.l2_norm();
+        break;
+      }
+      catch (libMesh::LogicError & e)
+      {
+        ref_point(0) -= mult * update(0);
+        if (dim - 1 == 2)
+          ref_point(1) -= mult * update(1);
+
+        mult *= 0.5;
+        if (mult < 1e-6)
+        {
+#ifndef NDEBUG
+          mooseWarning("We could not solve for the contact point.", e.what());
+#endif
+          update_size = update.l2_norm();
+          d = (secondary_point - phys_point[0]) * mult;
+          break;
+        }
+      }
+    }
+    // We failed the line search, make sure to trigger the error
+    if (mult < 1e-6)
+    {
+      nit = max_newton_its;
+      update_size = 1;
+      break;
+    }
   }
 
-  /*
-    if (nit == 12 && update_size > TOLERANCE*TOLERANCE)
-      Moose::err<<"Warning!  Newton solve for contact point failed to converge!"<<std::endl;
-  */
+  if (nit == max_newton_its && update_size > tolerance_newton)
+  {
+    search_succeeded = false;
+#ifndef NDEBUG
+    const auto initial_point =
+        start_with_centroid ? side->vertex_average() : ref_point = p_info._closest_point_ref;
+    Moose::err << "Warning!  Newton solve for contact point failed to converge!\nLast update "
+                  "distance was: "
+               << update_size << "\nInitial point guess:   " << initial_point
+               << "\nLast considered point: " << phys_point[0]
+               << "\nThis potential contact pair (face, point) will be discarded." << std::endl;
+#endif
+    return;
+  }
 
   p_info._closest_point_ref = ref_point;
   p_info._closest_point = phys_point[0];
@@ -205,7 +252,8 @@ findContactPoint(PenetrationInfo & p_info,
   if (dim - 1 == 2)
   {
     p_info._normal = dxyz_dxi[0].cross(dxyz_deta[0]);
-    p_info._normal /= p_info._normal.norm();
+    if (!MooseUtils::absoluteFuzzyEqual(p_info._normal.norm(), 0))
+      p_info._normal /= p_info._normal.norm();
   }
   else
   {

framework/src/geomsearch/PenetrationThread.C

@@ -129,6 +129,7 @@ PenetrationThread::operator()(const NodeIdRange & range)
         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,
@@ -137,7 +138,8 @@ PenetrationThread::operator()(const NodeIdRange & range)
                                 secondary_pos[0],
                                 false,
                                 _tangential_tolerance,
-                                contact_point_on_side);
+                                contact_point_on_side,
+                                search_succeeded);
 
         // Restore the original reference coordinates
         info->_closest_point_ref = contact_ref;
@@ -150,6 +152,7 @@ PenetrationThread::operator()(const NodeIdRange & range)
       {
         Real old_tangential_distance(info->_tangential_distance);
         bool contact_point_on_side(false);
+        bool search_succeeded = false;
 
         Moose::findContactPoint(*info,
                                 fe_elem,
@@ -158,7 +161,8 @@ PenetrationThread::operator()(const NodeIdRange & range)
                                 node,
                                 false,
                                 _tangential_tolerance,
-                                contact_point_on_side);
+                                contact_point_on_side,
+                                search_succeeded);
 
         if (contact_point_on_side)
         {
@@ -211,7 +215,7 @@ PenetrationThread::operator()(const NodeIdRange & range)
       else if (p_info.size() > 1)
       {
 
-        // Loop through all pairs of faces, and check for contact on ridge betweeen each face pair
+        // 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)
@@ -1723,18 +1727,23 @@ PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemIn
                                                      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);
+                            contact_point_on_side,
+                            search_succeeded);
 
-    thisElemInfo.push_back(pen_info);
-
-    p_info.push_back(pen_info);
+    // Do not add contact info from failed searches
+    if (search_succeeded)
+    {
+      thisElemInfo.push_back(pen_info);
+      p_info.push_back(pen_info);
+    }
   }
 }
 

test/tests/geomsearch/3d_penetration_locator/3d_rings.i

@@ -0,0 +1,45 @@
+[Mesh]
+  [fmg]
+    type = FileMeshGenerator
+    file = rings_facing.e
+  []
+  allow_renumbering = false
+[]
+
+[AuxVariables]
+  [penet]
+    family = MONOMIAL
+    order = CONSTANT
+    block = structure
+  []
+[]
+
+[AuxKernels]
+  [penet]
+    type = PenetrationAux
+    variable = penet
+    boundary = 'bottom_structure'
+    paired_boundary = 'top_heat_source'
+  []
+[]
+
+[VectorPostprocessors]
+  [values]
+    type = ElementValueSampler
+    variable = penet
+    sort_by = 'id'
+    block = 'structure'
+  []
+[]
+
+[Executioner]
+  type = Steady
+[]
+
+[Problem]
+  solve = false
+[]
+
+[Outputs]
+  csv = true
+[]