return;
  }

  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 & 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());
  }

  // If there's no derivatives then there is nothing to do. Moreover, if we pass zero size column

      // 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>

 */
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))}};
}
}
}

  _mortar_interface_coupling.clear();
  _secondary_node_to_nodal_normal.clear();
  _secondary_node_to_hh_nodal_tangents.clear();
  _nodal_geometry_coordinate_snapshot.clear();
  _secondary_element_to_secondary_lowerd_element.clear();
  _secondary_elems_to_mortar_segments.clear();
  _secondary_ip_sub_ids.clear();

    // Rows attached to a secondary node or either volume parent may need the complete secondary
    // face-star coupling, even when this lower-dimensional face is owned by another processor.
    std::set<processor_id_type> secondary_face_coupling_processors{
        secondary_elem->processor_id(),
        secondary_elem->interior_parent()->processor_id(),
        primary_elem->interior_parent()->processor_id()};
    for (const auto n : secondary_elem->node_index_range())
      secondary_face_coupling_processors.insert(secondary_elem->node_ptr(n)->processor_id());

    // LowerSecondary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(),
                                            secondary_elem->interior_parent()->id());

    // LowerPrimary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(),
                                            primary_elem->interior_parent()->id());

    // Lower-LowerDimensionalPrimary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(), primary_elem->id());

    // SecondaryLower
    coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(

    // PrimarySecondary
    coupling_info[primary_elem->interior_parent()->processor_id()].emplace_back(
        primary_elem->interior_parent()->id(), secondary_elem->interior_parent()->id());
  }

  // A face pair may produce several mortar segments. Remove repeated entries before communication.
  for (auto & [_, couplings] : coupling_info)
  {
    std::sort(couplings.begin(), couplings.end());
    couplings.erase(std::unique(couplings.begin(), couplings.end()), couplings.end());
  }

  // Push the coupling information

  /// The _periodic flag tells us whether we want to inward vs outward facing normals
  Real sign = _periodic ? -1 : 1;

  _nodal_geometry_coordinate_snapshot.clear();

  // First loop over lower-dimensional secondary side elements and compute/save the outward normal
  // for each one. We loop over all active elements currently, but this procedure could be

    // Save the coordinates used for the Real nodal geometry. The AD path later adds displacement
    // derivatives to these same values so both paths use the same geometry state.
    auto parent_side_elem = interior_parent->build_side_ptr(s);
    for (const auto & node : parent_side_elem->node_ref_range())
      _nodal_geometry_coordinate_snapshot.emplace(&node, node);

    // Reinit the face FE object on side s.
    nnx_fe_face->reinit(interior_parent, s);

    nodal_normal = nodal_normal.unit();

    _secondary_node_to_nodal_normal[_mesh.node_ptr(node_id)] = nodal_normal;
    _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)] =
        Moose::Mortar::householderTangents(nodal_normal);
  }
}

void
AutomaticMortarGeneration::computeADNodalNormals(
    const std::function<ADPoint(const Node &, const Point &)> & coordinate,
    std::unordered_map<const Node *, ADRealVectorValue> & nodal_normals) const
{
  const auto dim = _mesh.mesh_dimension();
  mooseAssert(dim == 2 || dim == 3, "AD nodal normals are only valid for 2D or 3D mortar.");

  QNodal qface(dim - 1);
  const auto & qweights = qface.get_weights();
  const Real sign = _periodic ? -1 : 1;

  std::unordered_map<const Node *, ADRealVectorValue> weighted_area_vectors;
  std::unordered_map<const Node *, Real> weighted_area_magnitudes;

  // AD propagates coordinate derivatives through the parent-side area vectors and normalization.
  for (const auto secondary_subdomain_id : _secondary_boundary_subdomain_ids)
    for (const Elem * const secondary_elem :
         as_range(_mesh.active_subdomain_elements_begin(secondary_subdomain_id),
                  _mesh.active_subdomain_elements_end(secondary_subdomain_id)))
    {
      const Elem * const interior_parent = secondary_elem->interior_parent();
      mooseAssert(interior_parent,
                  "No interior parent exists for element "
                      << secondary_elem->id()
                      << ". There may be a problem with your sideset set-up.");

      const Real element_size = secondary_elem->hmax();
      if (element_size <= 0)
        mooseError("Cannot compute an AD nodal normal for secondary mortar element ",
                   secondary_elem->id(),
                   " because it has a non-positive characteristic length.");

      const auto side = interior_parent->which_side_am_i(secondary_elem);
      FEType fe_type(secondary_elem->default_order(), LAGRANGE);
      std::unique_ptr<FEBase> fe_face(FEBase::build(dim, fe_type));
      fe_face->attach_quadrature_rule(&qface);

      const auto & face_normals = fe_face->get_normals();
      const auto & face_points = fe_face->get_xyz();
      const auto & map_dpsidxi = fe_face->get_fe_map().get_dpsidxi();
      const auto * const map_dpsideta = dim == 3 ? &fe_face->get_fe_map().get_dpsideta() : nullptr;

      fe_face->reinit(interior_parent, side);

      const auto qpoint_to_secondary_node =
          nodalQuadraturePointToSecondaryNodeMap(*secondary_elem, face_points);

      mooseAssert(face_normals.size() == face_points.size() &&
                      qweights.size() == face_points.size(),

      // Use the parent-side map rather than a standalone lower-dimensional basis. This preserves
      // enriched mappings such as the TET14-to-TRI7 face.
      auto parent_side_elem = interior_parent->build_side_ptr(side);
      const unsigned int n_mapping_shape_functions =
          dim == 2
              ? FE<2, LAGRANGE>::n_dofs(parent_side_elem.get(), parent_side_elem->default_order())
              : FE<3, LAGRANGE>::n_dofs(parent_side_elem.get(), parent_side_elem->default_order());
      if (map_dpsidxi.size() != n_mapping_shape_functions ||
          (map_dpsideta && map_dpsideta->size() != n_mapping_shape_functions))
        mooseError("The parent-side geometric map for secondary mortar element ",
                   secondary_elem->id(),
                   " of type ",
                   libMesh::Utility::enum_to_string<ElemType>(secondary_elem->type()),
                   " does not have one mapping shape function per element node.");

      std::vector<ADPoint> side_points;
      side_points.reserve(n_mapping_shape_functions);
      for (const auto side_node_index : make_range(n_mapping_shape_functions))
      {
        const Node & node = parent_side_elem->node_ref(side_node_index);
        side_points.push_back(
            coordinate(node, libmesh_map_find(_nodal_geometry_coordinate_snapshot, &node)));
      }

      for (const auto qp : make_range(face_points.size()))
      {
        const auto secondary_node_index = qpoint_to_secondary_node[qp];
        const Node * const secondary_node = secondary_elem->node_ptr(secondary_node_index);

        ADRealVectorValue tangent_xi;
        ADRealVectorValue tangent_eta;
        for (const auto side_node_index : make_range(n_mapping_shape_functions))
        {
          tangent_xi.add_scaled(side_points[side_node_index], map_dpsidxi[side_node_index][qp]);
          if (dim == 3)
            tangent_eta.add_scaled(side_points[side_node_index],
                                   (*map_dpsideta)[side_node_index][qp]);
        }

        auto area_vector = dim == 2 ? ADRealVectorValue(tangent_xi(1), -tangent_xi(0), 0)
                                    : tangent_xi.cross(tangent_eta);
        const Real area = MetaPhysicL::raw_value(area_vector.norm());
        const Real degeneracy_tolerance =
            dim == 2 ? TOLERANCE * element_size : TOLERANCE * element_size * element_size;
        if (area <= degeneracy_tolerance)
          mooseError("Cannot compute an AD nodal normal on degenerate secondary mortar element ",
                     secondary_elem->id(),
                     ".");

        // qweight * area_vector is equivalent to JxW * face_normal. Fix its orientation to the
        // Real-valued face normal because the incident-face set and orientation are not
        // differentiated.
        const Real orientation =
            MetaPhysicL::raw_value(area_vector) * face_normals[qp] < 0 ? -1 : 1;

        weighted_area_vectors[secondary_node] += sign * orientation * qweights[qp] * area_vector;
        weighted_area_magnitudes[secondary_node] += std::abs(qweights[qp]) * area;
      }
    }

  nodal_normals.clear();
  nodal_normals.reserve(weighted_area_vectors.size());
  for (auto & [secondary_node, weighted_area_vector] : weighted_area_vectors)
  {
    const Real weighted_area_vector_norm = MetaPhysicL::raw_value(weighted_area_vector.norm());
    const Real weighted_area_magnitude = libmesh_map_find(weighted_area_magnitudes, secondary_node);
    if (weighted_area_vector_norm <= TOLERANCE * weighted_area_magnitude)
      mooseError("Cannot compute a normalized AD nodal normal for secondary node ",
                 secondary_node->id(),
                 " because its incident face contributions cancel to a near-zero weighted normal.");

    auto nodal_normal = weighted_area_vector.unit();
    const auto & stored_normal = libmesh_map_find(_secondary_node_to_nodal_normal, secondary_node);
    mooseAssert((MetaPhysicL::raw_value(nodal_normal) - stored_normal).norm() < 100 * TOLERANCE,
                "The stored and AD secondary nodal normals must use the same geometry state.");
    // Set the raw value from the stored contact direction while retaining the derivatives produced
    // by the equivalent AD area-vector calculation.
    for (const auto component : make_range(3u))
      nodal_normal(component).value() = stored_normal(component);
    nodal_normals.emplace(secondary_node, std::move(nodal_normal));
  }
}

// Project secondary nodes onto their corresponding primary elements for each primary/secondary

                                      obj_params.get<SubdomainName>("secondary_subdomain");
                                  rm_params.set<SubdomainName>("primary_subdomain") =
                                      obj_params.get<SubdomainName>("primary_subdomain");
                                  rm_params.set<bool>("ghost_point_neighbors") =
                                      obj_params.get<bool>("ghost_point_neighbors");
                                  rm_params.set<bool>("ghost_higher_d_neighbors") =
                                      obj_params.get<bool>("ghost_higher_d_neighbors");
                                });

  for (const auto coupled_elem_id : coupled_set)
  {
    // Coupling data can precede a mesh rebuild, so IDs absent from the current mesh are ignored.
    const Elem * const coupled_elem = _mesh->query_elem_ptr(coupled_elem_id);
    if (coupled_elem && coupled_elem->processor_id() != p)
      coupled_elements.emplace(coupled_elem, _null_mat);
  }
}

  // Node ownership can follow a higher-dimensional parent, so the query may be a secondary face,
  // its interior parent, or a paired primary parent. Start from every associated secondary face so
  // the ghosting graph contains the complete remote nodal-normal support in either row direction.
  std::vector<const Elem *> secondary_lower_elem_candidates;
  if (query_elem->subdomain_id() == secondary_subdomain_id)
    secondary_lower_elem_candidates.push_back(query_elem);

  // Coupling functors can run before mortar data is rebuilt after mesh refinement. Use the
  // relationship manager's current mesh for face topology instead of AMG-owned element pointers.
  for (const auto side : query_elem->side_index_range())
    if (const Elem * const lower_elem = _moose_mesh->getLowerDElem(query_elem, side);
        lower_elem && lower_elem->subdomain_id() == secondary_subdomain_id)
      secondary_lower_elem_candidates.push_back(lower_elem);

  const auto & mic = amg.mortarInterfaceCoupling();
  if (const auto find_it = mic.find(query_elem->id()); find_it != mic.end())
    for (const auto coupled_elem_id : find_it->second)
      if (const Elem * const coupled_elem = _mesh->query_elem_ptr(coupled_elem_id);
          coupled_elem && coupled_elem->subdomain_id() == secondary_subdomain_id)
        secondary_lower_elem_candidates.push_back(coupled_elem);

  std::set<dof_id_type> secondary_lower_elems_handled;
  for (const Elem * const candidate : secondary_lower_elem_candidates)
  {
    if (candidate->subactive())
      continue;

    // Resolve inactive ancestors to the active faces used by the current coupling graph.
    std::vector<const Elem *> active_secondary_lower_elems;
    candidate->active_family_tree(active_secondary_lower_elems);
    for (const Elem * const secondary_lower_elem : active_secondary_lower_elems)
    {
      if (secondary_lower_elem->subdomain_id() != secondary_subdomain_id ||
          !secondary_lower_elems_handled.insert(secondary_lower_elem->id()).second)
        continue;

      std::set<const Elem *> secondary_lower_elem_point_neighbors;
      secondary_lower_elem->find_point_neighbors(secondary_lower_elem_point_neighbors);
      secondary_lower_elem_point_neighbors.insert(secondary_lower_elem);

      for (const Elem * const neigh : secondary_lower_elem_point_neighbors)
      {
        if (!neigh->active() || neigh->subdomain_id() != secondary_subdomain_id)
          continue;

        // Ghosting functors return only elements not owned by the requesting processor.
        if (neigh->processor_id() != p)
          coupled_elements.emplace(neigh, _null_mat);

        // Every point-neighbor face contributes to the smoothed nodal normal, even when it has no
        // mortar segment. Its interior parent therefore contributes displacement derivatives.
        const Elem * const interior_parent = neigh->interior_parent();
        mooseAssert(interior_parent,
                    "A lower-dimensional secondary mortar element must have an interior parent.");
        if (interior_parent->processor_id() != p)
          coupled_elements.emplace(interior_parent, _null_mat);

        ghostMortarInterfaceCouplings(p, neigh, coupled_elements, amg);
      }
    }
  }
}

void