Include nodal normal derivatives in mortar contact Jacobians

framework/doc/content/automatic_differentiation/index.md

@@ -191,6 +191,13 @@ i.e. the "sparse size" (stored as a `_dynamic_N` data member) of its
 data containers will never exceed what is required for the run-time problem,
 e.g. 18 for the 2D second-order solid mechanics example.
 
+When a batch of AD residual rows is constrained and cached, different rows may
+depend on different sparse sets of degrees of freedom. MOOSE first forms the
+union of those derivative indices so that every row uses one local matrix
+column layout. Missing row entries are zero while libMesh applies constraints;
+only the resulting nonzero entries are then cached. Residual batches with the
+same derivative support retain the direct common-layout path.
+
 ## AD in MOOSE
 
 As mentioned in above, MetaPhysicL is a forward-mode [!ac](AD)

framework/doc/content/source/constraints/AutomaticMortarGeneration.md

@@ -16,6 +16,26 @@ that have support from the shape function associated with a given node. In the f
 ghosting functors to the mortar segment mesh, which would allow us to delete unnecessary non-local
 elements.
 
+## Nodal normal derivatives
+
+`AutomaticMortarGeneration` constructs a unit normal at each secondary mortar node from the
+nodal-quadrature normals and weights of the incident secondary faces,
+
+!equation
+\boldsymbol{a}_A = \sum_{e \in \operatorname{star}(A)}
+  J_{eA}\boldsymbol{n}_{eA}, \qquad
+\boldsymbol{n}_A = \frac{\boldsymbol{a}_A}{\lVert\boldsymbol{a}_A\rVert}.
+
+For supported mechanical-contact Jacobian evaluations, the face tangents are formed from AD nodal
+coordinates. The equivalent quadrature-weighted oriented area vector is formed from an edge
+rotation in 2D or a tangent cross product in 3D. AD then differentiates the area weighting,
+summation over the secondary face star, and final normalization. The stored nodal-normal value is
+retained, so residual and Jacobian evaluations use the same contact direction.
+
+The incident-face set and orientation signs are fixed during an evaluation. Mortar segment
+topology, projections, integration measures, moving overlap boundaries, and active-set changes are
+not differentiated.
+
 ## 2D
 
 Generation of the 2D mortar segment mesh is outlined in [!cite](osti_1468630) and follows the

framework/include/base/Assembly.h

@@ -24,6 +24,7 @@
 #include "libmesh/numeric_vector.h"
 #include "libmesh/elem_side_builder.h"
 
+#include <algorithm>
 #include <unordered_map>
 
 // libMesh forward declarations
@@ -3131,23 +3132,43 @@ Assembly::cacheJacobian(const Residuals & residuals,
     return;
   }
 
-  const auto & compare_dofs = residuals[0].derivatives().nude_indices();
-#ifndef NDEBUG
-  auto compare_dofs_set = std::set<dof_id_type>(compare_dofs.begin(), compare_dofs.end());
-
+  const auto & first_dofs = residuals[0].derivatives().nude_indices();
+  bool supports_match = true;
   for (const auto i : make_range(decltype(residuals.size())(1), residuals.size()))
   {
-    const auto & residual = residuals[i];
-    auto current_dofs_set = std::set<dof_id_type>(residual.derivatives().nude_indices().begin(),
-                                                  residual.derivatives().nude_indices().end());
-    mooseAssert(compare_dofs_set == current_dofs_set,
-                "We're going to see whether the dof sets are the same. IIRC the degree of freedom "
-                "dependence (as indicated by the dof index set held by the ADReal) has to be the "
-                "same for every residual passed to this method otherwise constrain_element_matrix "
-                "will not work.");
+    const auto & current_dofs = residuals[i].derivatives().nude_indices();
+    // MetaPhysicL stores sparse derivative indices in sorted order, so equal index arrays identify
+    // rows with the same derivative support regardless of insertion order.
+    if (current_dofs.size() != first_dofs.size() ||
+        !std::equal(first_dofs.begin(), first_dofs.end(), current_dofs.begin()))
+    {
+      supports_match = false;
+      break;
+    }
+  }
+
+  // Keep the common-layout path whenever every row has the same derivative support.
+  if (supports_match)
+    _column_indices.assign(first_dofs.begin(), first_dofs.end());
+  else
+  {
+    // Constraining the local matrix requires one column layout shared by every residual row.
+    // Gather the union of the sparse AD supports, then leave entries absent from a row as zero.
+    std::size_t combined_support_size = 0;
+    for (const auto i : index_range(residuals))
+      combined_support_size += residuals[i].derivatives().nude_indices().size();
+
+    _column_indices.clear();
+    _column_indices.reserve(combined_support_size);
+    for (const auto i : index_range(residuals))
+    {
+      const auto & current_dofs = residuals[i].derivatives().nude_indices();
+      _column_indices.insert(_column_indices.end(), current_dofs.begin(), current_dofs.end());
+    }
+    std::sort(_column_indices.begin(), _column_indices.end());
+    _column_indices.erase(std::unique(_column_indices.begin(), _column_indices.end()),
+                          _column_indices.end());
   }
-#endif
-  _column_indices.assign(compare_dofs.begin(), compare_dofs.end());
 
   // If there's no derivatives then there is nothing to do. Moreover, if we pass zero size column
   // indices to constrain_element_matrix then we will potentially get errors out of BLAS
@@ -3170,7 +3191,11 @@ Assembly::cacheJacobian(const Residuals & residuals,
 
   for (const auto i : index_range(_row_indices))
     for (const auto j : index_range(_column_indices))
-      cacheJacobian(_row_indices[i], _column_indices[j], _element_matrix(i, j), {}, matrix_tags);
+      // Constraints may turn entries that were absent from an original AD row into nonzero
+      // contributions. Cache those entries while avoiding structural zeros introduced by the
+      // union above.
+      if (supports_match || _element_matrix(i, j) != 0.0)
+        cacheJacobian(_row_indices[i], _column_indices[j], _element_matrix(i, j), {}, matrix_tags);
 }
 
 template <typename Residuals, typename Indices>

framework/include/constraints/AutomaticMortarGeneration.h

@@ -21,9 +21,14 @@
 #include "libmesh/equation_systems.h"
 #include "libmesh/elem.h"
 #include "libmesh/int_range.h"
+#include "libmesh/vector_value.h"
+
+#include "metaphysicl/raw_type.h"
 
 // C++ includes
 #include <array>
+#include <cmath>
+#include <functional>
 #include <optional>
 #include <set>
 #include <memory>
@@ -48,10 +53,46 @@ using libMesh::Node;
 using libMesh::Point;
 using libMesh::Real;
 using libMesh::subdomain_id_type;
+using libMesh::VectorValue;
 
 typedef boundary_id_type BoundaryID;
 typedef subdomain_id_type SubdomainID;
 
+namespace Moose
+{
+namespace Mortar
+{
+/**
+ * Construct the two tangent vectors used by mortar from a unit normal. The template keeps the
+ * Real-valued and AD contact geometry on the same Householder chart and therefore gives them
+ * identical values.
+ * See Lopes, Silva, and Ambrosio, Computer-Aided Design 45(3), 2013, pp. 683-694.
+ */
+template <typename Vector>
+std::array<Vector, 2>
+householderTangents(const Vector & normal)
+{
+  mooseAssert(MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(normal.norm()), 1),
+              "The input nodal normal should have unity norm");
+
+  const Vector h_vector(normal(0) + 1.0, normal(1), normal(2));
+
+  // This chart is singular at (-1, 0, 0), where a constant orthonormal basis with zero
+  // derivatives is used.
+  if (std::abs(MetaPhysicL::raw_value(h_vector(0))) < TOLERANCE)
+    return {{Vector(0, 1, 0), Vector(0, 0, -1)}};
+
+  const auto h = h_vector.norm();
+  return {{Vector(-2.0 * h_vector(0) * h_vector(1) / (h * h),
+                  1.0 - 2.0 * h_vector(1) * h_vector(1) / (h * h),
+                  -2.0 * h_vector(1) * h_vector(2) / (h * h)),
+           Vector(-2.0 * h_vector(0) * h_vector(2) / (h * h),
+                  -2.0 * h_vector(1) * h_vector(2) / (h * h),
+                  1.0 - 2.0 * h_vector(2) * h_vector(2) / (h * h))}};
+}
+}
+}
+
 /**
  * Parent-face reference coordinates associated with the vertices of one triangular mortar segment.
  */
@@ -233,6 +274,17 @@ class AutomaticMortarGeneration : public ConsoleStreamInterface
   std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
   getNodalTangents(const Elem & secondary_elem) const;
 
+  /**
+   * Build the normalized JxW-weighted secondary nodal normals from AD nodal coordinates.
+   *
+   * The coordinate functor combines the node's displacement degrees of freedom with the coordinate
+   * snapshot used to compute the stored nodal geometry. This keeps residual values and their
+   * derivatives evaluated at the same geometry state.
+   */
+  void
+  computeADNodalNormals(const std::function<ADPoint(const Node &, const Point &)> & coordinate,
+                        std::unordered_map<const Node *, ADRealVectorValue> & nodal_normals) const;
+
   /**
    * Compute on-the-fly mapping from secondary interior parent nodes to lower dimensional nodes
    * @return The map from secondary interior parent nodes to lower dimensional nodes
@@ -501,6 +553,9 @@ class AutomaticMortarGeneration : public ConsoleStreamInterface
   /// (Householder approach).
   std::unordered_map<const Node *, std::array<Point, 2>> _secondary_node_to_hh_nodal_tangents;
 
+  /// Coordinates used to construct the stored nodal normals
+  std::unordered_map<const Node *, Point> _nodal_geometry_coordinate_snapshot;
+
   /// Map from full dimensional secondary element id to lower dimensional secondary element
   std::unordered_map<dof_id_type, const Elem *> _secondary_element_to_secondary_lowerd_element;
 
@@ -536,12 +591,6 @@ class AutomaticMortarGeneration : public ConsoleStreamInterface
   void projectPrimaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id,
                                      SubdomainID lower_dimensional_secondary_subdomain_id);
 
-  /**
-   * Householder orthogonalization procedure to obtain proper basis for tangent and binormal vectors
-   */
-  void
-  householderOrthogolization(const Point & normal, Point & tangent_one, Point & tangent_two) const;
-
   /**
    * Process aligned nodes
    * @returns whether mortar segment(s) were created

framework/include/relationshipmanagers/AugmentSparsityOnInterface.h

@@ -72,13 +72,13 @@ class AugmentSparsityOnInterface : public RelationshipManager
   /**
    * Query the mortar interface couplings of the query element. If a lower dimensional secondary
    * element is found, then we search for its point neighbors, which we ghost, as well as all of the
-   * mortar interface couplings of the point neighbors. This kind of ghosting is required for mortar
-   * nodal auxiliary kernels
+   * point neighbors, their interior parents, and the point neighbors' mortar interface couplings.
+   * This kind of ghosting is required for mortar nodal auxiliary kernels and nodal-normal contact
+   * Jacobians.
    */
   void ghostLowerDSecondaryElemPointNeighbors(const processor_id_type p,
                                               const Elem * const query_elem,
                                               map_type & coupled_elements,
-                                              BoundaryID secondary_boundary_id,
                                               SubdomainID secondary_subdomain_id,
                                               const AutomaticMortarGeneration & amg) const;
 
@@ -98,8 +98,8 @@ class AugmentSparsityOnInterface : public RelationshipManager
   /// the matrix sparsity pattern
   const bool _is_coupling_functor;
 
-  /// Whether to ghost point neighbors of secondary lower subdomain elements and their
-  /// cross mortar interface counterparts for applications such as mortar nodal auxiliary kernels
+  /// Whether to ghost secondary face point neighbors, their interior parents, and their
+  /// cross-interface counterparts
   const bool _ghost_point_neighbors;
 
   /// Whether we should ghost higher-dimensional neighbors. This is necessary when we are doing