AutomaticMortarGeneration.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

#include "AutomaticMortarGeneration.h"
#include "MortarSegmentInfo.h"
#include "NanoflannMeshAdaptor.h"
#include "MooseError.h"
#include "MooseTypes.h"
#include "MooseLagrangeHelpers.h"
#include "MortarSegmentHelper.h"
#include "MortarUtils.h"
#include "FormattedTable.h"
#include "FEProblemBase.h"
#include "DisplacedProblem.h"
#include "Output.h"

#include "libmesh/mesh_tools.h"
#include "libmesh/explicit_system.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/elem.h"
#include "libmesh/node.h"
#include "libmesh/dof_map.h"
#include "libmesh/edge_edge2.h"
#include "libmesh/edge_edge3.h"
#include "libmesh/face_tri3.h"
#include "libmesh/face_tri6.h"
#include "libmesh/face_tri7.h"
#include "libmesh/face_quad4.h"
#include "libmesh/face_quad8.h"
#include "libmesh/face_quad9.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/quadrature_nodal.h"
#include "libmesh/distributed_mesh.h"
#include "libmesh/replicated_mesh.h"
#include "libmesh/enum_to_string.h"
#include "libmesh/statistics.h"
#include "libmesh/equation_systems.h"

#include "metaphysicl/dualnumber.h"

#include "timpi/communicator.h"
#include "timpi/parallel_sync.h"

#include <array>
#include <algorithm>
#include <cmath>
#include <limits>

using namespace libMesh;
using MetaPhysicL::DualNumber;

// Make newer nanoflann API spelling compatible with older nanoflann
// versions
#if NANOFLANN_VERSION < 0x150
namespace nanoflann
{
typedef SearchParams SearchParameters;
}
#endif

namespace
{
// QNodal on a parent side returns normals, weights, and physical points in the
// parent-side quadrature ordering. That ordering is not guaranteed to match the
// node ordering of the generated lower-dimensional secondary element, especially
// for higher-order faces. Build the association geometrically so each
// quadrature value is attached to the secondary node at the same physical point.
std::vector<unsigned int>
nodalQuadraturePointToSecondaryNodeMap(const Elem & secondary_elem,
                                       const std::vector<Point> & q_points)
{
  const auto n_nodes = secondary_elem.n_nodes();
  if (q_points.size() != n_nodes)
    mooseError("Nodal quadrature produced ",
               q_points.size(),
               " points for secondary mortar element ",
               secondary_elem.id(),
               " of type ",
               libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
               ", but the element has ",
               n_nodes,
               " nodes.");

  const auto invalid_node = std::numeric_limits<unsigned int>::max();
  std::vector<unsigned int> qpoint_to_node(n_nodes, invalid_node);
  std::vector<bool> node_used(n_nodes, false);

  const Real element_size = secondary_elem.hmax();
  mooseAssert(element_size > 0,
              "Secondary mortar element "
                  << secondary_elem.id() << " of type "
                  << libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type())
                  << " has a non-positive hmax and cannot be used for nodal quadrature point "
                     "matching.");

  // The nodal quadrature locations and the generated secondary nodes are two floating-point
  // reconstructions of the same physical points. Scale the tolerance by element size so the
  // matching is insensitive to coordinate magnitude; the 100*TOLERANCE factor allows roundoff
  // from FE reinitialization and mesh generation while remaining far below a valid node spacing.
  const Real matching_tol = 100 * TOLERANCE * element_size;
  const Real matching_tol_sq = matching_tol * matching_tol;

  // Each nodal quadrature point should coincide with exactly one still-unused
  // secondary node. The unused-node search makes the mapping one-to-one and
  // avoids silently assigning two quadrature entries to the same node.
  for (const auto qp : make_range(q_points.size()))
  {
    unsigned int closest_node = invalid_node;
    Real closest_dist_sq = std::numeric_limits<Real>::max();
    Real second_closest_dist_sq = std::numeric_limits<Real>::max();

    for (const auto n : make_range(n_nodes))
    {
      if (node_used[n])
        continue;

      const Real dist_sq = (q_points[qp] - secondary_elem.point(n)).norm_sq();
      if (dist_sq < closest_dist_sq)
      {
        second_closest_dist_sq = closest_dist_sq;
        closest_dist_sq = dist_sq;
        closest_node = n;
      }
      else if (dist_sq < second_closest_dist_sq)
        second_closest_dist_sq = dist_sq;
    }

    if (closest_node == invalid_node || closest_dist_sq > matching_tol_sq)
      mooseError("Could not match nodal quadrature point ",
                 qp,
                 " at ",
                 q_points[qp],
                 " to a node on secondary mortar element ",
                 secondary_elem.id(),
                 " of type ",
                 libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
                 ". The nearest unmatched node distance is ",
                 std::sqrt(closest_dist_sq),
                 ", which exceeds the tolerance ",
                 matching_tol,
                 ".");

    if (second_closest_dist_sq <= matching_tol_sq)
      mooseError("Nodal quadrature point ",
                 qp,
                 " at ",
                 q_points[qp],
                 " does not map uniquely to secondary mortar element ",
                 secondary_elem.id(),
                 " of type ",
                 libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
                 ". Two unmatched nodes are within the matching tolerance ",
                 matching_tol,
                 ".");

    qpoint_to_node[qp] = closest_node;
    node_used[closest_node] = true;
  }

#ifdef DEBUG
  // In optimized builds the mapping above skips already matched nodes for speed. In debug builds,
  // audit the full candidate set to catch ambiguous geometry or accidental many-to-one matches.
  std::vector<unsigned int> node_to_qpoint(n_nodes, invalid_node);
  for (const auto qp : make_range(q_points.size()))
  {
    const auto mapped_node = qpoint_to_node[qp];
    mooseAssert(mapped_node != invalid_node && mapped_node < n_nodes,
                "Invalid secondary node mapping for nodal quadrature point " << qp << ".");
    mooseAssert(node_to_qpoint[mapped_node] == invalid_node,
                "Secondary node " << mapped_node << " on mortar element " << secondary_elem.id()
                                  << " was matched to both nodal quadrature point "
                                  << node_to_qpoint[mapped_node] << " and " << qp << ".");
    node_to_qpoint[mapped_node] = qp;

    // Check the qp -> node direction without excluding nodes already matched by previous qps.
    unsigned int candidate_count = 0;
    unsigned int candidate_node = invalid_node;
    for (const auto n : make_range(n_nodes))
      if ((q_points[qp] - secondary_elem.point(n)).norm_sq() <= matching_tol_sq)
      {
        ++candidate_count;
        candidate_node = n;
      }

    mooseAssert(candidate_count == 1,
                "Nodal quadrature point " << qp << " on mortar element " << secondary_elem.id()
                                          << " has " << candidate_count
                                          << " secondary node candidates within tolerance "
                                          << matching_tol << ".");
    mooseAssert(candidate_node == mapped_node,
                "Nodal quadrature point " << qp << " on mortar element " << secondary_elem.id()
                                          << " was matched to node " << mapped_node
                                          << ", but the full candidate search found node "
                                          << candidate_node << ".");
  }

  for (const auto n : make_range(n_nodes))
  {
    mooseAssert(node_to_qpoint[n] != invalid_node,
                "Secondary node " << n << " on mortar element " << secondary_elem.id()
                                  << " was not matched to a nodal quadrature point.");

    // Check the node -> qp direction so every secondary node is also uniquely represented.
    unsigned int candidate_count = 0;
    unsigned int candidate_qp = invalid_node;
    for (const auto qp : make_range(q_points.size()))
      if ((q_points[qp] - secondary_elem.point(n)).norm_sq() <= matching_tol_sq)
      {
        ++candidate_count;
        candidate_qp = qp;
      }

    mooseAssert(candidate_count == 1,
                "Secondary node " << n << " on mortar element " << secondary_elem.id() << " has "
                                  << candidate_count
                                  << " nodal quadrature point candidates within tolerance "
                                  << matching_tol << ".");
    mooseAssert(candidate_qp == node_to_qpoint[n],
                "Secondary node " << n << " on mortar element " << secondary_elem.id()
                                  << " was matched to nodal quadrature point " << node_to_qpoint[n]
                                  << ", but the full candidate search found point " << candidate_qp
                                  << ".");
  }
#endif

  return qpoint_to_node;
}
}

class MortarNodalGeometryOutput : public Output
{
public:
  static InputParameters validParams()
  {
    auto params = Output::validParams();
    params.addPrivateParam<AutomaticMortarGeneration *>("_amg", nullptr);
    params.addPrivateParam<MooseApp *>(MooseBase::app_param, nullptr);
    params.set<std::string>(MooseBase::type_param) = "MortarNodalGeometryOutput";
    return params;
  };

  MortarNodalGeometryOutput(const InputParameters & params)
    : Output(params), _amg(*getCheckedPointerParam<AutomaticMortarGeneration *>("_amg"))
  {
  }

  void output() override
  {
    // Must call compute_nodal_geometry first!
    if (_amg._secondary_node_to_nodal_normal.empty() ||
        _amg._secondary_node_to_hh_nodal_tangents.empty())
      mooseError("No entries found in the secondary node -> nodal geometry map.");

    auto & problem = _app.feProblem();
    auto & subproblem = _amg._on_displaced
                            ? static_cast<SubProblem &>(*problem.getDisplacedProblem())
                            : static_cast<SubProblem &>(problem);
    auto & nodal_normals_es = subproblem.es();

    const std::string nodal_normals_sys_name = "nodal_normals";

    if (!_nodal_normals_system)
    {
      for (const auto s : make_range(nodal_normals_es.n_systems()))
        if (!nodal_normals_es.get_system(s).is_initialized())
          // This is really early on in the simulation and the systems have not been initialized. We
          // thus need to avoid calling reinit on systems that haven't even had their first init yet
          return;

      _nodal_normals_system =
          &nodal_normals_es.template add_system<ExplicitSystem>(nodal_normals_sys_name);
      _nnx_var_num = _nodal_normals_system->add_variable("nodal_normal_x", FEType(FIRST, LAGRANGE)),
      _nny_var_num = _nodal_normals_system->add_variable("nodal_normal_y", FEType(FIRST, LAGRANGE));
      _nnz_var_num = _nodal_normals_system->add_variable("nodal_normal_z", FEType(FIRST, LAGRANGE));

      _t1x_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_1_x", FEType(FIRST, LAGRANGE)),
      _t1y_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_1_y", FEType(FIRST, LAGRANGE));
      _t1z_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_1_z", FEType(FIRST, LAGRANGE));

      _t2x_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_2_x", FEType(FIRST, LAGRANGE)),
      _t2y_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_2_y", FEType(FIRST, LAGRANGE));
      _t2z_var_num =
          _nodal_normals_system->add_variable("nodal_tangent_2_z", FEType(FIRST, LAGRANGE));
      nodal_normals_es.reinit();
    }

    const DofMap & dof_map = _nodal_normals_system->get_dof_map();
    std::vector<dof_id_type> dof_indices_nnx, dof_indices_nny, dof_indices_nnz;
    std::vector<dof_id_type> dof_indices_t1x, dof_indices_t1y, dof_indices_t1z;
    std::vector<dof_id_type> dof_indices_t2x, dof_indices_t2y, dof_indices_t2z;

    for (MeshBase::const_element_iterator el = _amg._mesh.elements_begin(),
                                          end_el = _amg._mesh.elements_end();
         el != end_el;
         ++el)
    {
      const Elem * elem = *el;

      // Get the nodal dofs for this Elem.
      dof_map.dof_indices(elem, dof_indices_nnx, _nnx_var_num);
      dof_map.dof_indices(elem, dof_indices_nny, _nny_var_num);
      dof_map.dof_indices(elem, dof_indices_nnz, _nnz_var_num);

      dof_map.dof_indices(elem, dof_indices_t1x, _t1x_var_num);
      dof_map.dof_indices(elem, dof_indices_t1y, _t1y_var_num);
      dof_map.dof_indices(elem, dof_indices_t1z, _t1z_var_num);

      dof_map.dof_indices(elem, dof_indices_t2x, _t2x_var_num);
      dof_map.dof_indices(elem, dof_indices_t2y, _t2y_var_num);
      dof_map.dof_indices(elem, dof_indices_t2z, _t2z_var_num);

      //

      // For each node of the Elem, if it is in the secondary_node_to_nodal_normal
      // container, set the corresponding nodal normal dof values.
      for (MooseIndex(elem->n_vertices()) n = 0; n < elem->n_vertices(); ++n)
      {
        auto it = _amg._secondary_node_to_nodal_normal.find(elem->node_ptr(n));
        if (it != _amg._secondary_node_to_nodal_normal.end())
        {
          _nodal_normals_system->solution->set(dof_indices_nnx[n], it->second(0));
          _nodal_normals_system->solution->set(dof_indices_nny[n], it->second(1));
          _nodal_normals_system->solution->set(dof_indices_nnz[n], it->second(2));
        }

        auto it_tangent = _amg._secondary_node_to_hh_nodal_tangents.find(elem->node_ptr(n));
        if (it_tangent != _amg._secondary_node_to_hh_nodal_tangents.end())
        {
          _nodal_normals_system->solution->set(dof_indices_t1x[n], it_tangent->second[0](0));
          _nodal_normals_system->solution->set(dof_indices_t1y[n], it_tangent->second[0](1));
          _nodal_normals_system->solution->set(dof_indices_t1z[n], it_tangent->second[0](2));

          _nodal_normals_system->solution->set(dof_indices_t2x[n], it_tangent->second[1](0));
          _nodal_normals_system->solution->set(dof_indices_t2y[n], it_tangent->second[1](1));
          _nodal_normals_system->solution->set(dof_indices_t2z[n], it_tangent->second[1](2));
        }

      } // end loop over nodes
    } // end loop over elems

    // Finish assembly.
    _nodal_normals_system->solution->close();

    std::set<std::string> sys_names = {nodal_normals_sys_name};

    // Write the nodal normals to file
    ExodusII_IO nodal_normals_writer(_amg._mesh);

    // Default to non-HDF5 output for wider compatibility
    nodal_normals_writer.set_hdf5_writing(false);

    nodal_normals_writer.write_equation_systems(
        "nodal_geometry_only.e", nodal_normals_es, &sys_names);
  }

private:
  AutomaticMortarGeneration & _amg;


  libMesh::System * _nodal_normals_system = nullptr;
  unsigned int _nnx_var_num;
  unsigned int _nny_var_num;
  unsigned int _nnz_var_num;

  unsigned int _t1x_var_num;
  unsigned int _t1y_var_num;
  unsigned int _t1z_var_num;

  unsigned int _t2x_var_num;
  unsigned int _t2y_var_num;
  unsigned int _t2z_var_num;
};

AutomaticMortarGeneration::AutomaticMortarGeneration(
    MooseApp & app,
    MeshBase & mesh_in,
    const std::pair<BoundaryID, BoundaryID> & boundary_key,
    const std::pair<SubdomainID, SubdomainID> & subdomain_key,
    bool on_displaced,
    bool periodic,
    const bool debug,
    const bool correct_edge_dropping,
    const Real minimum_projection_angle,
    const Mortar3DSubpatchPlane mortar_3d_subpatch_plane,
    const MortarSegmentTriangulationMode triangulation_mode,
    const bool triangulate_triangles,
    const Mortar3DQuadraturePointMapping mortar_3d_qp_mapping)
  : ConsoleStreamInterface(app),
    _app(app),
    _mesh(mesh_in),
    _debug(debug),
    _on_displaced(on_displaced),
    _periodic(periodic),
    // 3D mortar always builds the mortar segment mesh distributedly (each rank adds only its local
    // secondary elements). For 2D, we ghost the entire mortar interface when displaced, so
    // displaced meshes are always replicated; otherwise follow the parent mesh.
    _distributed(_mesh.mesh_dimension() == 3 ? true : (!_on_displaced && !_mesh.is_replicated())),
    _correct_edge_dropping(correct_edge_dropping),
    _minimum_projection_angle(minimum_projection_angle),
    _mortar_3d_subpatch_plane(mortar_3d_subpatch_plane),
    _triangulation_mode(triangulation_mode),
    _triangulate_triangles(triangulate_triangles),
    _mortar_3d_qp_mapping(mortar_3d_qp_mapping)
{
  _primary_secondary_boundary_id_pairs.push_back(boundary_key);
  _primary_requested_boundary_ids.insert(boundary_key.first);
  _secondary_requested_boundary_ids.insert(boundary_key.second);
  _primary_secondary_subdomain_id_pairs.push_back(subdomain_key);
  _primary_boundary_subdomain_ids.insert(subdomain_key.first);
  _secondary_boundary_subdomain_ids.insert(subdomain_key.second);

  if (_distributed)
    _mortar_segment_mesh =
        std::make_unique<DistributedMesh>(_mesh.comm(), _mesh.spatial_dimension());
  else
    _mortar_segment_mesh =
        std::make_unique<ReplicatedMesh>(_mesh.comm(), _mesh.spatial_dimension());
}

std::string
AutomaticMortarGeneration::mortarInterfaceName() const
{
  std::vector<std::string> string_vec(_primary_secondary_boundary_id_pairs.size() * 2 + 1);
  for (const auto i : index_range(_primary_secondary_boundary_id_pairs))
  {
    const auto [primary_bnd_id, secondary_bnd_id] = _primary_secondary_boundary_id_pairs[i];
    string_vec[2 * i] = std::to_string(primary_bnd_id);
    string_vec[2 * i + 1] = std::to_string(secondary_bnd_id);
  }
  string_vec.back() = _on_displaced ? "displaced" : "undisplaced";
  return MooseUtils::join(string_vec, "_");
}

void
AutomaticMortarGeneration::initOutput()
{
  if (!_debug)
    return;

  _output_params = std::make_unique<InputParameters>(MortarNodalGeometryOutput::validParams());
  _output_params->set<AutomaticMortarGeneration *>("_amg") = this;
  _output_params->set<FEProblemBase *>("_fe_problem_base") = &_app.feProblem();
  _output_params->set<MooseApp *>(MooseBase::app_param) = &_app;
  _output_params->set<std::string>(MooseBase::name_param) =
      "mortar_nodal_geometry_" + mortarInterfaceName();
  _output_params->finalize("MortarNodalGeometryOutput");
  _app.getOutputWarehouse().addOutput(std::make_shared<MortarNodalGeometryOutput>(*_output_params));
}

void
AutomaticMortarGeneration::clear()
{
  _msm_elem_to_reference_points.clear();
  _mortar_segment_mesh->clear();
  _nodes_to_secondary_elem_map.clear();
  _nodes_to_primary_elem_map.clear();
  _secondary_node_and_elem_to_xi2_primary_elem.clear();
  _primary_node_and_elem_to_xi1_secondary_elem.clear();
  _msm_elem_to_info.clear();
  _lower_elem_to_side_id.clear();
  _mortar_interface_coupling.clear();
  _secondary_node_to_nodal_normal.clear();
  _secondary_node_to_hh_nodal_tangents.clear();
  _secondary_element_to_secondary_lowerd_element.clear();
  _secondary_elems_to_mortar_segments.clear();
  _secondary_ip_sub_ids.clear();
  _primary_ip_sub_ids.clear();
  _projected_secondary_nodes.clear();
  _failed_secondary_node_projections.clear();
}

const MortarSegmentReferencePoints &
AutomaticMortarGeneration::mortarSegmentReferencePoints(const Elem & mortar_segment_elem) const
{
  if (_mortar_3d_qp_mapping != Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION)
    mooseError("Mortar segment reference points were requested for mortar segment element ",
               mortar_segment_elem.id(),
               ", but the reference-interpolation mapping mode is not enabled.");

  const auto reference_points_it = _msm_elem_to_reference_points.find(&mortar_segment_elem);
  if (reference_points_it == _msm_elem_to_reference_points.end())
    mooseError("No reference-point record was found for mortar segment element ",
               mortar_segment_elem.id(),
               ". The mortar segment info and reference-point maps are not aligned.");

  return reference_points_it->second;
}

void
AutomaticMortarGeneration::buildNodeToElemMaps()
{
  if (_secondary_requested_boundary_ids.empty() || _primary_requested_boundary_ids.empty())
    mooseError(
        "Must specify secondary and primary boundary ids before building node-to-elem maps.");

  // Construct nodes_to_secondary_elem_map
  for (const auto & secondary_elem :
       as_range(_mesh.active_elements_begin(), _mesh.active_elements_end()))
  {
    // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
    if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
      continue;

    for (const auto & nd : secondary_elem->node_ref_range())
    {
      std::vector<const Elem *> & vec = _nodes_to_secondary_elem_map[nd.id()];
      vec.push_back(secondary_elem);
    }
  }

  // Construct nodes_to_primary_elem_map
  for (const auto & primary_elem :
       as_range(_mesh.active_elements_begin(), _mesh.active_elements_end()))
  {
    // If this is not one of the lower-dimensional primary side elements, go on to the next one.
    if (!this->_primary_boundary_subdomain_ids.count(primary_elem->subdomain_id()))
      continue;

    for (const auto & nd : primary_elem->node_ref_range())
    {
      std::vector<const Elem *> & vec = _nodes_to_primary_elem_map[nd.id()];
      vec.push_back(primary_elem);
    }
  }
}

std::vector<Point>
AutomaticMortarGeneration::getNodalNormals(const Elem & secondary_elem) const
{
  std::vector<Point> nodal_normals(secondary_elem.n_nodes());
  for (const auto n : make_range(secondary_elem.n_nodes()))
    nodal_normals[n] = _secondary_node_to_nodal_normal.at(secondary_elem.node_ptr(n));

  return nodal_normals;
}

const Elem *
AutomaticMortarGeneration::getSecondaryLowerdElemFromSecondaryElem(
    dof_id_type secondary_elem_id) const
{
  mooseAssert(_secondary_element_to_secondary_lowerd_element.count(secondary_elem_id),
              "Map should locate secondary element");

  return _secondary_element_to_secondary_lowerd_element.at(secondary_elem_id);
}

std::map<unsigned int, unsigned int>
AutomaticMortarGeneration::getSecondaryIpToLowerElementMap(const Elem & lower_secondary_elem) const
{
  std::map<unsigned int, unsigned int> secondary_ip_i_to_lower_secondary_i;
  const Elem * const secondary_ip = lower_secondary_elem.interior_parent();
  mooseAssert(secondary_ip, "This should be non-null");

  for (const auto i : make_range(lower_secondary_elem.n_nodes()))
  {
    const auto & nd = lower_secondary_elem.node_ref(i);
    secondary_ip_i_to_lower_secondary_i[secondary_ip->get_node_index(&nd)] = i;
  }

  return secondary_ip_i_to_lower_secondary_i;
}

std::map<unsigned int, unsigned int>
AutomaticMortarGeneration::getPrimaryIpToLowerElementMap(
    const Elem & lower_primary_elem,
    const Elem & primary_elem,
    const Elem & /*lower_secondary_elem*/) const
{
  std::map<unsigned int, unsigned int> primary_ip_i_to_lower_primary_i;

  for (const auto i : make_range(lower_primary_elem.n_nodes()))
  {
    const auto & nd = lower_primary_elem.node_ref(i);
    primary_ip_i_to_lower_primary_i[primary_elem.get_node_index(&nd)] = i;
  }

  return primary_ip_i_to_lower_primary_i;
}

std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
AutomaticMortarGeneration::getNodalTangents(const Elem & secondary_elem) const
{
  // MetaPhysicL will check if we ran out of allocated space.
  MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_one(0);
  MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_two(0);

  for (const auto n : make_range(secondary_elem.n_nodes()))
  {
    const auto & tangent_vectors =
        libmesh_map_find(_secondary_node_to_hh_nodal_tangents, secondary_elem.node_ptr(n));
    nodal_tangents_one.push_back(tangent_vectors[0]);
    nodal_tangents_two.push_back(tangent_vectors[1]);
  }

  return {{nodal_tangents_one, nodal_tangents_two}};
}

std::vector<Point>
AutomaticMortarGeneration::getNormals(const Elem & secondary_elem,
                                      const std::vector<Real> & oned_xi1_pts) const
{
  std::vector<Point> xi1_pts(oned_xi1_pts.size());
  for (const auto qp : index_range(oned_xi1_pts))
    xi1_pts[qp] = oned_xi1_pts[qp];

  return getNormals(secondary_elem, xi1_pts);
}

std::vector<Point>
AutomaticMortarGeneration::getNormals(const Elem & secondary_elem,
                                      const std::vector<Point> & xi1_pts) const
{
  const auto mortar_dim = _mesh.mesh_dimension() - 1;
  const auto num_qps = xi1_pts.size();
  const auto nodal_normals = getNodalNormals(secondary_elem);
  std::vector<Point> normals(num_qps);

  for (const auto n : make_range(secondary_elem.n_nodes()))
    for (const auto qp : make_range(num_qps))
    {
      const auto phi =
          (mortar_dim == 1)
              ? Moose::fe_lagrange_1D_shape(secondary_elem.default_order(), n, xi1_pts[qp](0))
              : Moose::fe_lagrange_2D_shape(secondary_elem.type(),
                                            secondary_elem.default_order(),
                                            n,
                                            static_cast<const TypeVector<Real> &>(xi1_pts[qp]));
      normals[qp] += phi * nodal_normals[n];
    }

  if (_periodic)
    for (auto & normal : normals)
      normal *= -1;

  return normals;
}

void
AutomaticMortarGeneration::buildMortarSegmentMesh()
{
  using std::abs;

  dof_id_type local_id_index = 0;
  std::size_t node_unique_id_offset = 0;

  // Create an offset by the maximum number of mortar segment elements that can be created *plus*
  // the number of lower-dimensional secondary subdomain elements. Recall that the number of mortar
  // segments created is a function of node projection, *and* that if we split elems we will delete
  // that elem which has already taken a unique id
  for (const auto & pr : _primary_secondary_boundary_id_pairs)
  {
    const auto primary_bnd_id = pr.first;
    const auto secondary_bnd_id = pr.second;
    const auto num_primary_nodes =
        std::distance(_mesh.bid_nodes_begin(primary_bnd_id), _mesh.bid_nodes_end(primary_bnd_id));
    const auto num_secondary_nodes = std::distance(_mesh.bid_nodes_begin(secondary_bnd_id),
                                                   _mesh.bid_nodes_end(secondary_bnd_id));
    mooseAssert(num_primary_nodes,
                "There are no primary nodes on boundary ID "
                    << primary_bnd_id << ". Does that bondary ID even exist on the mesh?");
    mooseAssert(num_secondary_nodes,
                "There are no secondary nodes on boundary ID "
                    << secondary_bnd_id << ". Does that bondary ID even exist on the mesh?");

    node_unique_id_offset += num_primary_nodes + 2 * num_secondary_nodes;
  }

  // 1.) Add all lower-dimensional secondary side elements as the "initial" mortar segments.
  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
                                        end_el = _mesh.active_elements_end();
       el != end_el;
       ++el)
  {
    const Elem * secondary_elem = *el;

    // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
    if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
      continue;

    std::vector<Node *> new_nodes;
    for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
    {
      new_nodes.push_back(_mortar_segment_mesh->add_point(
          secondary_elem->point(n), secondary_elem->node_id(n), secondary_elem->processor_id()));
      Node * const new_node = new_nodes.back();
      new_node->set_unique_id(new_node->id() + node_unique_id_offset);
    }

    std::unique_ptr<Elem> new_elem;
    if (secondary_elem->default_order() == SECOND)
      new_elem = std::make_unique<Edge3>();
    else
      new_elem = std::make_unique<Edge2>();

    new_elem->processor_id() = secondary_elem->processor_id();
    new_elem->subdomain_id() = secondary_elem->subdomain_id();
    new_elem->set_id(local_id_index++);
    new_elem->set_unique_id(new_elem->id());

    for (MooseIndex(new_elem->n_nodes()) n = 0; n < new_elem->n_nodes(); ++n)
      new_elem->set_node(n, new_nodes[n]);

    Elem * new_elem_ptr = _mortar_segment_mesh->add_elem(new_elem.release());

    // The xi^(1) values for this mortar segment are initially -1 and 1.
    MortarSegmentInfo msinfo;
    msinfo.xi1_a = -1;
    msinfo.xi1_b = +1;
    msinfo.secondary_elem = secondary_elem;

    auto new_container_it0 = _secondary_node_and_elem_to_xi2_primary_elem.find(
             std::make_pair(secondary_elem->node_ptr(0), secondary_elem)),
         new_container_it1 = _secondary_node_and_elem_to_xi2_primary_elem.find(
             std::make_pair(secondary_elem->node_ptr(1), secondary_elem));

    bool new_container_node0_found =
             (new_container_it0 != _secondary_node_and_elem_to_xi2_primary_elem.end()),
         new_container_node1_found =
             (new_container_it1 != _secondary_node_and_elem_to_xi2_primary_elem.end());

    const Elem * node0_primary_candidate = nullptr;
    const Elem * node1_primary_candidate = nullptr;

    if (new_container_node0_found)
    {
      const auto & xi2_primary_elem_pair = new_container_it0->second;
      msinfo.xi2_a = xi2_primary_elem_pair.first;
      node0_primary_candidate = xi2_primary_elem_pair.second;
    }

    if (new_container_node1_found)
    {
      const auto & xi2_primary_elem_pair = new_container_it1->second;
      msinfo.xi2_b = xi2_primary_elem_pair.first;
      node1_primary_candidate = xi2_primary_elem_pair.second;
    }

    // If both node0 and node1 agree on the primary element they are
    // projected into, then this mortar segment fits entirely within
    // a single primary element, and we can go ahead and set the
    // msinfo.primary_elem pointer now.
    if (node0_primary_candidate == node1_primary_candidate)
      msinfo.primary_elem = node0_primary_candidate;

    // Associate this MSM elem with the MortarSegmentInfo.
    _msm_elem_to_info.emplace(new_elem_ptr, msinfo);

    // Maintain the mapping between secondary elems and mortar segment elems contained within them.
    // Initially, only the original secondary_elem is present.
    _secondary_elems_to_mortar_segments[secondary_elem->id()].insert(new_elem_ptr);
  }

  // 2.) Insert new nodes from primary side and split mortar segments as necessary.
  for (const auto & pr : _primary_node_and_elem_to_xi1_secondary_elem)
  {
    auto key = pr.first;
    auto val = pr.second;

    const Node * primary_node = std::get<1>(key);
    Real xi1 = val.first;
    const Elem * secondary_elem = val.second;

    // If this is an aligned node, we don't need to do anything.
    if (abs(abs(xi1) - 1.) < _xi_tolerance)
      continue;

    auto && order = secondary_elem->default_order();

    // Determine physical location of new point to be inserted.
    Point new_pt(0);
    for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
      new_pt += Moose::fe_lagrange_1D_shape(order, n, xi1) * secondary_elem->point(n);

    // Find the current mortar segment that will have to be split.
    auto & mortar_segment_set = _secondary_elems_to_mortar_segments[secondary_elem->id()];
    Elem * current_mortar_segment = nullptr;
    MortarSegmentInfo * info = nullptr;

    for (const auto & mortar_segment_candidate : mortar_segment_set)
    {
      try
      {
        info = &_msm_elem_to_info.at(mortar_segment_candidate);
      }
      catch (std::out_of_range &)
      {
        mooseError("MortarSegmentInfo not found for the mortar segment candidate");
      }
      if (info->xi1_a <= xi1 && xi1 <= info->xi1_b)
      {
        current_mortar_segment = mortar_segment_candidate;
        break;
      }
    }

    // Make sure we found one.
    if (current_mortar_segment == nullptr)
      mooseError("Unable to find appropriate mortar segment during linear search!");

    // If node lands on endpoint of segment, don't split.
    // Jacob: This condition was getting missed by the < comparison a few lines above. To fix it I
    // just made it <= and put this condition in to handle equality different. It probably could be
    // done with a tolerance but the the toleranced equality is already handled later when we drop
    // segments with small volume.
    if (info->xi1_a == xi1 || xi1 == info->xi1_b)
      continue;

    const auto new_id = _mortar_segment_mesh->max_node_id();
    mooseAssert(_mortar_segment_mesh->comm().verify(new_id),
                "new_id must be the same on all processes");
    Node * const new_node =
        _mortar_segment_mesh->add_point(new_pt, new_id, secondary_elem->processor_id());
    new_node->set_unique_id(new_id + node_unique_id_offset);

    // Reconstruct the nodal normal at xi1. This will help us
    // determine the orientation of the primary elems relative to the
    // new mortar segments.
    const Point normal = getNormals(*secondary_elem, std::vector<Real>({xi1}))[0];

    // Get the set of primary_node neighbors.
    if (this->_nodes_to_primary_elem_map.find(primary_node->id()) ==
        this->_nodes_to_primary_elem_map.end())
      mooseError("We should already have built this primary node to elem pair!");
    const std::vector<const Elem *> & primary_node_neighbors =
        this->_nodes_to_primary_elem_map[primary_node->id()];

    // Sanity check
    if (primary_node_neighbors.size() == 0 || primary_node_neighbors.size() > 2)
      mooseError("We must have either 1 or 2 primary side nodal neighbors, but we had ",
                 primary_node_neighbors.size());

    // Primary Elem pointers which we will eventually assign to the
    // mortar segments being created.  We start by assuming
    // primary_node_neighbor[0] is on the "left" and
    // primary_node_neighbor[1]/"nothing" is on the "right" and then
    // swap them if that's not the case.
    const Elem * left_primary_elem = primary_node_neighbors[0];
    const Elem * right_primary_elem =
        (primary_node_neighbors.size() == 2) ? primary_node_neighbors[1] : nullptr;

    Real left_xi2 = MortarSegmentInfo::invalid_xi, right_xi2 = MortarSegmentInfo::invalid_xi;

    // Storage for z-component of cross products for determining
    // orientation.
    std::array<Real, 2> secondary_node_cps;
    std::vector<Real> primary_node_cps(primary_node_neighbors.size());

    // Store z-component of left and right secondary node cross products with the nodal normal.
    for (unsigned int nid = 0; nid < 2; ++nid)
      secondary_node_cps[nid] = normal.cross(secondary_elem->point(nid) - new_pt)(2);

    for (MooseIndex(primary_node_neighbors) mnn = 0; mnn < primary_node_neighbors.size(); ++mnn)
    {
      const Elem * primary_neigh = primary_node_neighbors[mnn];
      Point opposite = (primary_neigh->node_ptr(0) == primary_node) ? primary_neigh->point(1)
                                                                    : primary_neigh->point(0);
      Point cp = normal.cross(opposite - new_pt);
      primary_node_cps[mnn] = cp(2);
    }

    // We will verify that only 1 orientation is actually valid.
    bool orientation1_valid = false, orientation2_valid = false;

    if (primary_node_neighbors.size() == 2)
    {
      // 2 primary neighbor case
      orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.) &&
                           (secondary_node_cps[1] * primary_node_cps[1] > 0.);

      orientation2_valid = (secondary_node_cps[0] * primary_node_cps[1] > 0.) &&
                           (secondary_node_cps[1] * primary_node_cps[0] > 0.);
    }
    else if (primary_node_neighbors.size() == 1)
    {
      // 1 primary neighbor case
      orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.);
      orientation2_valid = (secondary_node_cps[1] * primary_node_cps[0] > 0.);
    }
    else
      mooseError("Invalid primary node neighbors size ", primary_node_neighbors.size());

    // Verify that both orientations are not simultaneously valid/invalid. If they are not, then we
    // are going to throw an exception instead of erroring out since we can easily reach this point
    // if we have one bad linear solve. It's better in general to catch the error and then try a
    // smaller time-step
    if (orientation1_valid && orientation2_valid)
      throw MooseException(
          "AutomaticMortarGeneration: Both orientations cannot simultaneously be valid.");

    // We are going to treat the case where both orientations are invalid as a case in which we
    // should not be splitting the mortar mesh to incorporate primary mesh elements.
    // In practice, this case has appeared for very oblique projections, so we assume these cases
    // will not be considered in mortar thermomechanical contact.
    if (!orientation1_valid && !orientation2_valid)
    {
      mooseDoOnce(mooseWarning(
          "AutomaticMortarGeneration: Unable to determine valid secondary-primary orientation. "
          "Consequently we will consider projection of the primary node invalid and not split the "
          "mortar segment. "
          "This situation can indicate there are very oblique projections between primary (mortar) "
          "and secondary (non-mortar) surfaces for a good problem set up. It can also mean your "
          "time step is too large. This message is only printed once."));
      continue;
    }

    // Make an Elem on the left
    std::unique_ptr<Elem> new_elem_left;
    if (order == SECOND)
      new_elem_left = std::make_unique<Edge3>();
    else
      new_elem_left = std::make_unique<Edge2>();

    new_elem_left->processor_id() = current_mortar_segment->processor_id();
    new_elem_left->subdomain_id() = current_mortar_segment->subdomain_id();
    new_elem_left->set_id(local_id_index++);
    new_elem_left->set_unique_id(new_elem_left->id());
    new_elem_left->set_node(0, current_mortar_segment->node_ptr(0));
    new_elem_left->set_node(1, new_node);

    // Make an Elem on the right
    std::unique_ptr<Elem> new_elem_right;
    if (order == SECOND)
      new_elem_right = std::make_unique<Edge3>();
    else
      new_elem_right = std::make_unique<Edge2>();

    new_elem_right->processor_id() = current_mortar_segment->processor_id();
    new_elem_right->subdomain_id() = current_mortar_segment->subdomain_id();
    new_elem_right->set_id(local_id_index++);
    new_elem_right->set_unique_id(new_elem_right->id());
    new_elem_right->set_node(0, new_node);
    new_elem_right->set_node(1, current_mortar_segment->node_ptr(1));

    if (order == SECOND)
    {
      // left
      Point left_interior_point(0);
      Real left_interior_xi = (xi1 + info->xi1_a) / 2;

      // This is eta for the current mortar segment that we're splitting
      Real current_left_interior_eta =
          (2. * left_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);

      for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
           n < current_mortar_segment->n_nodes();
           ++n)
        left_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_left_interior_eta) *
                               current_mortar_segment->point(n);

      const auto new_interior_left_id = _mortar_segment_mesh->max_node_id();
      mooseAssert(_mortar_segment_mesh->comm().verify(new_interior_left_id),
                  "new_id must be the same on all processes");
      Node * const new_interior_node_left = _mortar_segment_mesh->add_point(
          left_interior_point, new_interior_left_id, new_elem_left->processor_id());
      new_elem_left->set_node(2, new_interior_node_left);
      new_interior_node_left->set_unique_id(new_interior_left_id + node_unique_id_offset);

      // right
      Point right_interior_point(0);
      Real right_interior_xi = (xi1 + info->xi1_b) / 2;
      // This is eta for the current mortar segment that we're splitting
      Real current_right_interior_eta =
          (2. * right_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);

      for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
           n < current_mortar_segment->n_nodes();
           ++n)
        right_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_right_interior_eta) *
                                current_mortar_segment->point(n);

      const auto new_interior_id_right = _mortar_segment_mesh->max_node_id();
      mooseAssert(_mortar_segment_mesh->comm().verify(new_interior_id_right),
                  "new_id must be the same on all processes");
      Node * const new_interior_node_right = _mortar_segment_mesh->add_point(
          right_interior_point, new_interior_id_right, new_elem_right->processor_id());
      new_elem_right->set_node(2, new_interior_node_right);
      new_interior_node_right->set_unique_id(new_interior_id_right + node_unique_id_offset);
    }

    // If orientation 2 was valid, swap the left and right primaries.
    if (orientation2_valid)
      std::swap(left_primary_elem, right_primary_elem);

    // Now that we know left_primary_elem and right_primary_elem, we can determine left_xi2 and
    // right_xi2.
    if (left_primary_elem)
      left_xi2 = (primary_node == left_primary_elem->node_ptr(0)) ? -1 : +1;
    if (right_primary_elem)
      right_xi2 = (primary_node == right_primary_elem->node_ptr(0)) ? -1 : +1;

    // Grab the MortarSegmentInfo object associated with this
    // segment. We can use "at()" here since we want this to fail if
    // current_mortar_segment is not found... Since we're going to
    // erase this entry from the map momentarily, we make an actual
    // copy rather than grabbing a reference.
    auto msm_it = _msm_elem_to_info.find(current_mortar_segment);
    if (msm_it == _msm_elem_to_info.end())
      mooseError("MortarSegmentInfo not found for current_mortar_segment.");
    MortarSegmentInfo current_msinfo = msm_it->second;

    // add_left
    {
      Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_left.release());

      // Create new MortarSegmentInfo objects for new_elem_left
      MortarSegmentInfo new_msinfo_left;

      // The new MortarSegmentInfo info objects inherit their "outer"
      // information from current_msinfo and the rest is determined by
      // the Node being inserted.
      new_msinfo_left.xi1_a = current_msinfo.xi1_a;
      new_msinfo_left.xi2_a = current_msinfo.xi2_a;
      new_msinfo_left.secondary_elem = secondary_elem;
      new_msinfo_left.xi1_b = xi1;
      new_msinfo_left.xi2_b = left_xi2;
      new_msinfo_left.primary_elem = left_primary_elem;

      // Add new msinfo objects to the map.
      _msm_elem_to_info.emplace(msm_new_elem, new_msinfo_left);

      // We need to insert new_elem_left in
      // the mortar_segment_set for this secondary_elem.
      mortar_segment_set.insert(msm_new_elem);
    }

    // add_right
    {
      Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_right.release());

      // Create new MortarSegmentInfo objects for new_elem_right
      MortarSegmentInfo new_msinfo_right;

      new_msinfo_right.xi1_b = current_msinfo.xi1_b;
      new_msinfo_right.xi2_b = current_msinfo.xi2_b;
      new_msinfo_right.secondary_elem = secondary_elem;
      new_msinfo_right.xi1_a = xi1;
      new_msinfo_right.xi2_a = right_xi2;
      new_msinfo_right.primary_elem = right_primary_elem;

      _msm_elem_to_info.emplace(msm_new_elem, new_msinfo_right);

      mortar_segment_set.insert(msm_new_elem);
    }

    // Erase the MortarSegmentInfo object for current_mortar_segment from the map.
    _msm_elem_to_info.erase(msm_it);

    // current_mortar_segment must be erased from the
    // mortar_segment_set since it has now been split.
    mortar_segment_set.erase(current_mortar_segment);

    // The original mortar segment has been split, so erase it from
    // the mortar segment mesh.
    _mortar_segment_mesh->delete_elem(current_mortar_segment);
  }

  // Remove all MSM elements without a primary contribution
  for (auto msm_elem : _mortar_segment_mesh->active_element_ptr_range())
  {
    MortarSegmentInfo & msinfo = libmesh_map_find(_msm_elem_to_info, msm_elem);
    Elem * primary_elem = const_cast<Elem *>(msinfo.primary_elem);
    if (primary_elem == nullptr || abs(msinfo.xi2_a) > 1.0 + TOLERANCE ||
        abs(msinfo.xi2_b) > 1.0 + TOLERANCE)
    {
      // Erase from secondary to msms map
      auto it = _secondary_elems_to_mortar_segments.find(msinfo.secondary_elem->id());
      mooseAssert(it != _secondary_elems_to_mortar_segments.end(),
                  "We should have found the element");
      auto & msm_set = it->second;
      msm_set.erase(msm_elem);
      // We may be creating nodes with only one element neighbor where before this removal there
      // were two. But the nodal normal used in computations will reflect the two-neighbor geometry.
      // For a lower-d secondary mesh corner, that will imply the corner node will have a tilted
      // normal vector (same for tangents) despite the mortar segment mesh not including its
      // vertical neighboring element. It is the secondary element neighbors (not mortar segment
      // mesh neighbors) that determine the nodal normal field.
      if (msm_set.empty())
        _secondary_elems_to_mortar_segments.erase(it);

      // Erase msinfo
      _msm_elem_to_info.erase(msm_elem);

      // Remove element from mortar segment mesh
      _mortar_segment_mesh->delete_elem(msm_elem);
    }
    else
    {
      _secondary_ip_sub_ids.insert(msinfo.secondary_elem->interior_parent()->subdomain_id());
      _primary_ip_sub_ids.insert(msinfo.primary_elem->interior_parent()->subdomain_id());
    }
  }

  std::unordered_set<Node *> msm_connected_nodes;

  // Deleting elements may produce isolated nodes.
  // Loops for identifying and removing such nodes from mortar segment mesh.
  for (const auto & element : _mortar_segment_mesh->element_ptr_range())
    for (auto & n : element->node_ref_range())
      msm_connected_nodes.insert(&n);

  for (const auto & node : _mortar_segment_mesh->node_ptr_range())
    if (!msm_connected_nodes.count(node))
      _mortar_segment_mesh->delete_node(node);

#ifdef DEBUG
  // Verify that all segments without primary contribution have been deleted
  for (auto msm_elem : _mortar_segment_mesh->active_element_ptr_range())
  {
    const MortarSegmentInfo & msinfo = libmesh_map_find(_msm_elem_to_info, msm_elem);
    mooseAssert(msinfo.primary_elem != nullptr,
                "All mortar segment elements should have valid "
                "primary element.");
  }
#endif

  _mortar_segment_mesh->cache_elem_data();

  // (Optionally) Write the mortar segment mesh to file for inspection
  if (_debug)
    outputMortarMesh();

  buildCouplingInformation();
}

void
AutomaticMortarGeneration::outputMortarMesh()
{
  ExodusII_IO mortar_segment_mesh_writer(*_mortar_segment_mesh);

  // Default to non-HDF5 output for wider compatibility
  mortar_segment_mesh_writer.set_hdf5_writing(false);

  std::array<std::string, 3> file_pieces = {
      _app.getOutputFileBase(/*for_non_moose_build_output=*/true),
      mortarInterfaceName(),
      "mortar_segment_mesh.e"};
  mortar_segment_mesh_writer.write(MooseUtils::join(file_pieces, "_"));
}

void
AutomaticMortarGeneration::buildMortarSegmentMesh3d()
{
  const bool use_reference_interpolation =
      _mortar_3d_qp_mapping == Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION;

  // Add an integer flag to mortar segment mesh to keep track of which subelem
  // of second order primal elements mortar segments correspond to
  auto secondary_sub_elem = _mortar_segment_mesh->add_elem_integer("secondary_sub_elem");
  auto primary_sub_elem = _mortar_segment_mesh->add_elem_integer("primary_sub_elem");

  // Assign globally unique node/element IDs via an exclusive prefix scan: each rank's bound is
  // local_secondary_sub_elems * visible_primary_sub_elems * 9, where 9 is the maximum nodes a
  // single secondary/primary sub-element pair can produce (8-vertex clipped polygon + center).
  // The result is cached and invalidated by meshChanged(), so the allgather only runs on topology
  // changes, not on every displaced-mesh residual update.
  if (!_msm_node_id_start.has_value())
  {
    dof_id_type local_secondary_sub_elems = 0, visible_primary_sub_elems = 0;
    for (const auto & [primary_sub_id, secondary_sub_id] : _primary_secondary_subdomain_id_pairs)
    {
      for (const auto * const el :
           _mesh.active_local_subdomain_elements_ptr_range(secondary_sub_id))
        local_secondary_sub_elems += el->n_sub_elem();
      for (const auto * const el : _mesh.active_subdomain_elements_ptr_range(primary_sub_id))
        visible_primary_sub_elems += el->n_sub_elem();
    }
    const dof_id_type per_rank_bound = local_secondary_sub_elems * visible_primary_sub_elems * 9;
    std::vector<dof_id_type> per_rank_bounds;
    _mesh.comm().allgather(per_rank_bound, per_rank_bounds);
    dof_id_type start = 0;
    for (const auto r : make_range(_mesh.processor_id()))
      start += per_rank_bounds[r];
    _msm_node_id_start = start;
  }
  dof_id_type next_node_id = *_msm_node_id_start;
  // Element IDs use the same starting offset: node and element IDs are separately numbered, and
  // element count per clip (n triangles) is always <= node count (n+1), so per_rank_bound covers
  // both.
  dof_id_type next_elem_id = next_node_id;

  // Loop through mortar secondary and primary pairs to create mortar segment mesh between each
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
  {
    const auto primary_subd_id = pr.first;
    const auto secondary_subd_id = pr.second;

    // Build k-d tree for use in Step 1.2 for primary interface coarse screening
    NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, primary_subd_id);
    subdomain_kd_tree_t kd_tree(
        3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));

    // Construct the KD tree.
    kd_tree.buildIndex();

    // Return the unoriented geometric normal of a linearized subpatch. These expressions are the
    // TRI3 and QUAD4 mapping tangents evaluated at the reference center, equivalent to evaluating
    // the first-order finite-element normal there without constructing a temporary element.
    auto get_sub_elem_geometric_normal = [](const std::vector<Point> & nodes)
    {
      Point dxdxi;
      Point dxdeta;
      if (nodes.size() == 3)
      {
        dxdxi = nodes[1] - nodes[0];
        dxdeta = nodes[2] - nodes[0];
      }
      else if (nodes.size() == 4)
      {
        // Bilinear center tangents define one normal for the full quad instead of selecting one of
        // the two diagonal triangle normals.
        dxdxi = 0.25 * (nodes[1] + nodes[2] - nodes[0] - nodes[3]);
        dxdeta = 0.25 * (nodes[2] + nodes[3] - nodes[0] - nodes[1]);
      }
      else
        mooseError("GEOMETRIC_NORMAL 3D mortar subpatch plane construction only supports "
                   "triangular and quadrilateral subpatches, but received ",
                   nodes.size(),
                   " nodes.");

      Point geometric_normal = dxdxi.cross(dxdeta);
      const auto normal_norm = geometric_normal.norm();
      // The cross product has units of area, so compare it with the product of tangent lengths.
      // Their ratio is the sine of the included angle and is independent of the mesh length scale.
      if (normal_norm <= TOLERANCE * dxdxi.norm() * dxdeta.norm())
        mooseError("GEOMETRIC_NORMAL 3D mortar subpatch plane construction encountered a "
                   "degenerate subpatch.");

      geometric_normal /= normal_norm;
      return geometric_normal;
    };

    for (MeshBase::const_element_iterator el = _mesh.active_local_elements_begin(),
                                          end_el = _mesh.active_local_elements_end();
         el != end_el;
         ++el)
    {
      const Elem * secondary_side_elem = *el;

      const Real secondary_volume = secondary_side_elem->volume();

      // If this Elem is not in the current secondary subdomain, go on to the next one.
      if (secondary_side_elem->subdomain_id() != secondary_subd_id)
        continue;

      auto [secondary_elem_to_msm_map_it, insertion_happened] =
          _secondary_elems_to_mortar_segments.emplace(secondary_side_elem->id(),
                                                      std::set<Elem *, CompareDofObjectsByID>{});
      libmesh_ignore(insertion_happened);
      auto & secondary_to_msm_element_set = secondary_elem_to_msm_map_it->second;

      std::vector<std::unique_ptr<MortarSegmentHelper>> mortar_segment_helper(
          secondary_side_elem->n_sub_elem());
      const auto nodal_normals = getNodalNormals(*secondary_side_elem);

      for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
      {
        // Get indices of sub-element nodes in element
        const auto sub_elem_nodes =
            Moose::Mortar::getMortarSubElementNodeIndices(*secondary_side_elem, sel);

        // Secondary sub-element center, normal, and nodes
        Point center;
        Point normal;
        std::vector<Point> nodes(sub_elem_nodes.size());

        // Collect the sub-element points and evaluate its center and averaged nodal normal.
        for (auto iv : make_range(sub_elem_nodes.size()))
        {
          const auto n = sub_elem_nodes[iv];
          nodes[iv] = secondary_side_elem->point(n);
          center += secondary_side_elem->point(n);
          normal += nodal_normals[n];
        }
        center /= sub_elem_nodes.size();
        normal = normal.unit();

        if (_mortar_3d_subpatch_plane == Mortar3DSubpatchPlane::GEOMETRIC_NORMAL)
        {
          const Point averaged_normal = normal;
          normal = get_sub_elem_geometric_normal(nodes);
          if (normal * averaged_normal < 0)
            normal *= -1;
        }

        if (use_reference_interpolation)
        {
          std::vector<Point> sub_elem_reference_points;
          sub_elem_reference_points.reserve(sub_elem_nodes.size());
          for (const auto node_index : sub_elem_nodes)
            sub_elem_reference_points.push_back(secondary_side_elem->master_point(node_index));

          mortar_segment_helper[sel] =
              std::make_unique<MortarSegmentHelper>(std::move(nodes),
                                                    std::move(sub_elem_reference_points),
                                                    center,
                                                    normal,
                                                    _triangulation_mode,
                                                    _triangulate_triangles);
        }
        else
          mortar_segment_helper[sel] = std::make_unique<MortarSegmentHelper>(
              std::move(nodes), center, normal, _triangulation_mode, _triangulate_triangles);
      }

      // Search point for performing Nanoflann (k-d tree) searches.
      // In each case we use the center point of the original element (not sub-elements for second
      // order elements). This is to do search for all sub-elements simultaneously
      std::array<Real, 3> query_pt;
      Point center_point;
      switch (secondary_side_elem->type())
      {
        case TRI3:
        case QUAD4:
          center_point = mortar_segment_helper[0]->center();
          query_pt = {{center_point(0), center_point(1), center_point(2)}};
          break;
        case TRI6:
        case TRI7:
          center_point = mortar_segment_helper[1]->center();
          query_pt = {{center_point(0), center_point(1), center_point(2)}};
          break;
        case QUAD8:
          center_point = mortar_segment_helper[4]->center();
          query_pt = {{center_point(0), center_point(1), center_point(2)}};
          break;
        case QUAD9:
          center_point = secondary_side_elem->point(8);
          query_pt = {{center_point(0), center_point(1), center_point(2)}};
          break;
        default:
          mooseError(
              "Face element type: ", secondary_side_elem->type(), "not supported for 3D mortar");
      }

      // The number of results we want to get. These results will only be used to find
      // a single element with non-trivial overlap, after an element is identified a breadth
      // first search is done on neighbors
      const std::size_t num_results = 3;

      // Initialize result_set and do the search.
      std::vector<size_t> ret_index(num_results);
      std::vector<Real> out_dist_sqr(num_results);
      nanoflann::KNNResultSet<Real> result_set(num_results);
      result_set.init(&ret_index[0], &out_dist_sqr[0]);
      kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());

      // Initialize list of processed primary elements, we don't want to revisit processed elements
      std::set<const Elem *, CompareDofObjectsByID> processed_primary_elems;

      // Initialize candidate set and flag for switching between coarse screening and breadth-first
      // search
      bool primary_elem_found = false;
      std::set<const Elem *, CompareDofObjectsByID> primary_elem_candidates;
      const bool use_geometric_subpatch_normals =
          _mortar_3d_subpatch_plane == Mortar3DSubpatchPlane::GEOMETRIC_NORMAL;
      // In geometric mode the projection-angle cutoff also rejects near-orthogonal subpatch pairs.
      // The absolute dot product below keeps opposing primary/secondary orientations admissible.
      const Real minimum_subpatch_normal_alignment =
          use_geometric_subpatch_normals ? std::sin(_minimum_projection_angle * libMesh::pi / 180.0)
                                         : 0.0;

      // Loop candidate nodes (returned by Nanoflann) and add all adjoining elems to candidate set
      for (auto r : make_range(result_set.size()))
      {
        // Verify that the squared distance we compute is the same as nanoflann's
        mooseAssert(abs((_mesh.point(ret_index[r]) - center_point).norm_sq() - out_dist_sqr[r]) <=
                        TOLERANCE,
                    "Lower-dimensional element squared distance verification failed.");

        // Get list of elems connected to node
        std::vector<const Elem *> & node_elems =
            this->_nodes_to_primary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));

        // Uniquely add elems to candidate set
        for (auto elem : node_elems)
          primary_elem_candidates.insert(elem);
      }

      while (!primary_elem_candidates.empty())
      {
        const Elem * primary_elem_candidate = *primary_elem_candidates.begin();

        // If we've already processed this candidate, we don't need to check it again.
        if (processed_primary_elems.count(primary_elem_candidate))
        {
          primary_elem_candidates.erase(primary_elem_candidate);
          continue;
        }

        // Initialize set of nodes used to construct mortar segment elements
        std::vector<Point> nodal_points;

        // Initialize map from mortar segment elements to nodes
        std::vector<std::vector<unsigned int>> elem_to_node_map;

        // Initialize list of secondary and primary sub-elements that formed each mortar segment
        std::vector<std::pair<unsigned int, unsigned int>> sub_elem_map;
        std::vector<std::array<Point, 3>> elem_to_secondary_reference_points;
        std::vector<std::array<Point, 3>> elem_to_primary_reference_points;

        for (auto p_el : make_range(primary_elem_candidate->n_sub_elem()))
        {
          // Get nodes of primary sub-elements
          const auto sub_elem_nodes =
              Moose::Mortar::getMortarSubElementNodeIndices(*primary_elem_candidate, p_el);

          // Get list of primary sub-element vertex nodes
          std::vector<Point> primary_sub_elem(sub_elem_nodes.size());
          for (auto iv : make_range(sub_elem_nodes.size()))
          {
            const auto n = sub_elem_nodes[iv];
            primary_sub_elem[iv] = primary_elem_candidate->point(n);
          }
          Point primary_sub_elem_normal;
          if (use_geometric_subpatch_normals)
            primary_sub_elem_normal = get_sub_elem_geometric_normal(primary_sub_elem);

          std::vector<Point> sub_elem_reference_points;
          if (use_reference_interpolation)
          {
            sub_elem_reference_points.reserve(sub_elem_nodes.size());
            for (const auto node_index : sub_elem_nodes)
              sub_elem_reference_points.push_back(primary_elem_candidate->master_point(node_index));
          }

          // Loop through secondary sub-elements
          for (auto s_el : make_range(secondary_side_elem->n_sub_elem()))
          {
            // Nearby primary candidates can include adjacent corner faces. Those faces may clip to
            // numerical slivers, which we do not consider valid face-to-face mortar pairs for this
            // search.
            if (use_geometric_subpatch_normals &&
                std::abs(primary_sub_elem_normal * mortar_segment_helper[s_el]->normal()) <
                    minimum_subpatch_normal_alignment)
              continue;

            // Mortar segment helpers were defined for each secondary sub-element, they will:
            //  1. Project primary sub-element onto linearized secondary sub-element
            //  2. Clip projected primary sub-element against secondary sub-element
            //  3. Triangulate clipped polygon to form mortar segments
            //
            // Mortar segment helpers append a list of mortar segment nodes and connectivities that
            // can be directly used to build mortar segments
            const auto segments_before_helper = elem_to_node_map.size();
            if (use_reference_interpolation)
              mortar_segment_helper[s_el]->getMortarSegments(primary_sub_elem,
                                                             sub_elem_reference_points,
                                                             nodal_points,
                                                             elem_to_node_map,
                                                             elem_to_secondary_reference_points,
                                                             elem_to_primary_reference_points,
                                                             TOLERANCE * secondary_volume);
            else
              mortar_segment_helper[s_el]->getMortarSegments(
                  primary_sub_elem, nodal_points, elem_to_node_map);

            // Keep track of which secondary and primary sub-elements created segment
            for (auto i = segments_before_helper; i < elem_to_node_map.size(); ++i)
              sub_elem_map.push_back(std::make_pair(s_el, p_el));
          }
        }

        // Mark primary element as processed and remove from candidate list
        processed_primary_elems.insert(primary_elem_candidate);
        primary_elem_candidates.erase(primary_elem_candidate);

        // If overlap of polygons was non-trivial (created mortar segment elements)
        if (!elem_to_node_map.empty())
        {
          if (sub_elem_map.size() != elem_to_node_map.size())
            mooseError("The mortar segment subpatch map is not aligned with the mortar segment "
                       "connectivity map.");
          if (use_reference_interpolation &&
              (elem_to_secondary_reference_points.size() != elem_to_node_map.size() ||
               elem_to_primary_reference_points.size() != elem_to_node_map.size()))
            mooseError("The mortar segment reference-point maps are not aligned with the mortar "
                       "segment connectivity map.");

          // Only overlap polygons large enough to become mortar segments may switch the candidate
          // search to breadth first.
          bool seed_breadth_first_search = false;
          std::vector<bool> retained_mortar_segments(elem_to_node_map.size(), false);
          for (const auto el : index_range(elem_to_node_map))
          {
            const auto & node_map = elem_to_node_map[el];
            if (node_map.size() != 3)
              mooseError(
                  "Active mortar segments only supports TRI elements, 3 nodes expected but: ",
                  node_map.size(),
                  " provided.");

            const Point e1 = nodal_points[node_map[1]] - nodal_points[node_map[0]];
            const Point e2 = nodal_points[node_map[2]] - nodal_points[node_map[0]];
            retained_mortar_segments[el] =
                0.5 * e1.cross(e2).norm() / secondary_volume >= TOLERANCE;
            seed_breadth_first_search = seed_breadth_first_search || retained_mortar_segments[el];
          }

          if (seed_breadth_first_search)
          {
            // If this is the first element with a qualifying overlap, set flag. Candidates will
            // now be neighbors of elements that had qualifying overlap.
            if (!primary_elem_found)
            {
              primary_elem_found = true;
              primary_elem_candidates.clear();
            }

            // Add neighbors to candidate list
            for (auto neighbor : primary_elem_candidate->neighbor_ptr_range())
            {
              // If not valid or not on lower dimensional secondary subdomain, skip
              if (neighbor == nullptr || neighbor->subdomain_id() != primary_subd_id)
                continue;
              // If already processed, skip
              if (processed_primary_elems.count(neighbor))
                continue;
              // Otherwise, add to candidates
              primary_elem_candidates.insert(neighbor);
            }
          }

          std::vector<Node *> new_nodes;
          // Clipping can append points for triangles later rejected by the area tolerance. Add only
          // points referenced by retained triangles so the mortar mesh has no orphan nodes.
          std::vector<bool> retained_nodes(nodal_points.size(), false);
          for (const auto el : index_range(elem_to_node_map))
            if (retained_mortar_segments[el])
              for (const auto node : elem_to_node_map[el])
                retained_nodes[node] = true;

          new_nodes.resize(nodal_points.size(), nullptr);
          for (const auto node : index_range(nodal_points))
            if (retained_nodes[node])
              new_nodes[node] = _mortar_segment_mesh->add_point(
                  nodal_points[node], next_node_id++, secondary_side_elem->processor_id());

          // Loop through triangular elements in map
          for (auto el : index_range(elem_to_node_map))
          {
            if (!retained_mortar_segments[el])
              continue;

            std::unique_ptr<Elem> new_elem;
            if (elem_to_node_map[el].size() == 3)
              new_elem = std::make_unique<Tri3>();
            else
              mooseError("Active mortar segments only supports TRI elements, 3 nodes expected "
                         "but: ",
                         elem_to_node_map[el].size(),
                         " provided.");

            new_elem->processor_id() = secondary_side_elem->processor_id();
            new_elem->subdomain_id() = secondary_side_elem->subdomain_id();
            new_elem->set_id(next_elem_id++);

            // Attach newly created nodes
            for (auto i : index_range(elem_to_node_map[el]))
              new_elem->set_node(i, new_nodes[elem_to_node_map[el][i]]);

            // If element is smaller than tolerance, don't add to msm
            if (new_elem->volume() / secondary_volume < TOLERANCE)
              continue;

            // Add elements to mortar segment mesh
            Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem.release());

            msm_new_elem->set_extra_integer(secondary_sub_elem, sub_elem_map[el].first);
            msm_new_elem->set_extra_integer(primary_sub_elem, sub_elem_map[el].second);

            // Fill out mortar segment info
            MortarSegmentInfo msinfo;
            msinfo.secondary_elem = secondary_side_elem;
            msinfo.primary_elem = primary_elem_candidate;

            // Associate this MSM elem with the MortarSegmentInfo.
            _msm_elem_to_info.emplace(msm_new_elem, msinfo);

            // Store reference data only for retained segments.
            if (use_reference_interpolation)
            {
              MortarSegmentReferencePoints reference_points{elem_to_secondary_reference_points[el],
                                                            elem_to_primary_reference_points[el]};
              _msm_elem_to_reference_points.emplace(msm_new_elem, reference_points);
            }

            // Add this mortar segment to the secondary elem to mortar segment map
            secondary_to_msm_element_set.insert(msm_new_elem);

            _secondary_ip_sub_ids.insert(msinfo.secondary_elem->interior_parent()->subdomain_id());
            // Unlike for 2D, we always have a primary when building the mortar mesh so we don't
            // have to check for null
            _primary_ip_sub_ids.insert(msinfo.primary_elem->interior_parent()->subdomain_id());
          }
        }
        // End loop through primary element candidates
      }

      if (use_geometric_subpatch_normals)
      {
        // A geometric corner filter may intentionally leave individual subpatches uncovered. Warn
        // only when the complete secondary element failed to produce a retained segment.
        if (secondary_to_msm_element_set.empty())
          mooseDoOnce(
              mooseWarning("Some secondary elements on mortar interface were unable to identify"
                           " a corresponding primary element; this may be expected depending on"
                           " problem geometry but may indicate a failure of the element search"
                           " or projection"));
      }
      else
        for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
          if (mortar_segment_helper[sel]->remainder() == 1.0)
            mooseDoOnce(
                mooseWarning("Some secondary elements on mortar interface were unable to identify"
                             " a corresponding primary element; this may be expected depending on"
                             " problem geometry but may indicate a failure of the element search"
                             " or projection"));

      if (secondary_to_msm_element_set.empty())
        _secondary_elems_to_mortar_segments.erase(secondary_elem_to_msm_map_it);
    } // End loop through secondary elements
  } // End loop through mortar constraint pairs

  mooseAssert(!use_reference_interpolation ||
                  _msm_elem_to_reference_points.size() == _msm_elem_to_info.size(),
              "Mortar segment info and reference-point maps must remain aligned.");

  _mortar_segment_mesh->cache_elem_data();

  // The mesh was built distributedly (each rank owns only its local elements), so mark it
  // as such so MeshSerializer correctly gathers it to proc 0 for Exodus output.
  _mortar_segment_mesh->set_distributed();

  // Output mortar segment mesh
  if (_debug)
  {
    // If element is not triangular, increment subdomain id
    // (ExodusII does not support mixed element types in a single subdomain)
    for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
      if (msm_el->type() != TRI3)
        msm_el->subdomain_id()++;

    outputMortarMesh();

    // Undo increment
    for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
      if (msm_el->type() != TRI3)
        msm_el->subdomain_id()--;
  }

  buildCouplingInformation();

  // Print mortar segment mesh statistics
  if (_debug)
  {
    msmStatistics();
  }
}

void
AutomaticMortarGeneration::buildCouplingInformation()
{
  std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type>>>
      coupling_info;

  // Loop over the msm_elem_to_info object and build a bi-directional
  // multimap from secondary elements to the primary Elems which they are
  // coupled to and vice-versa. This is used in the
  // AugmentSparsityOnInterface functor to determine whether a given
  // secondary Elem is coupled across the mortar interface to a primary
  // element.
  for (const auto & pr : _msm_elem_to_info)
  {
    const Elem * secondary_elem = pr.second.secondary_elem;
    const Elem * primary_elem = pr.second.primary_elem;

    // LowerSecondary
    coupling_info[secondary_elem->processor_id()].emplace_back(
        secondary_elem->id(), secondary_elem->interior_parent()->id());
    if (secondary_elem->processor_id() != _mesh.processor_id())
      // We want to keep information for nonlocal lower-dimensional secondary element point
      // neighbors for mortar nodal aux kernels
      _mortar_interface_coupling[secondary_elem->id()].insert(
          secondary_elem->interior_parent()->id());

    // LowerPrimary
    coupling_info[secondary_elem->processor_id()].emplace_back(
        secondary_elem->id(), primary_elem->interior_parent()->id());
    if (secondary_elem->processor_id() != _mesh.processor_id())
      // We want to keep information for nonlocal lower-dimensional secondary element point
      // neighbors for mortar nodal aux kernels
      _mortar_interface_coupling[secondary_elem->id()].insert(
          primary_elem->interior_parent()->id());

    // Lower-LowerDimensionalPrimary
    coupling_info[secondary_elem->processor_id()].emplace_back(secondary_elem->id(),
                                                               primary_elem->id());
    if (secondary_elem->processor_id() != _mesh.processor_id())
      // We want to keep information for nonlocal lower-dimensional secondary element point
      // neighbors for mortar nodal aux kernels
      _mortar_interface_coupling[secondary_elem->id()].insert(primary_elem->id());

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

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

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

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

  // Push the coupling information
  auto action_functor =
      [this](processor_id_type,
             const std::vector<std::pair<dof_id_type, dof_id_type>> & coupling_info)
  {
    for (auto [i, j] : coupling_info)
      _mortar_interface_coupling[i].insert(j);
  };
  TIMPI::push_parallel_vector_data(_mesh.comm(), coupling_info, action_functor);
}

std::vector<AutomaticMortarGeneration::MsmSubdomainStats>
AutomaticMortarGeneration::computeMsmStatistics()
{
  std::vector<MsmSubdomainStats> result;
  StatisticsVector<Real> primary;
  StatisticsVector<Real> secondary;
  StatisticsVector<Real> msm;
  std::unordered_map<dof_id_type, Real> primary_elems_to_volume;

  for (const auto & [primary_subd_id, secondary_subd_id] : _primary_secondary_subdomain_id_pairs)
  {
    for (const auto * const secondary_el :
         _mesh.active_local_subdomain_element_ptr_range(secondary_subd_id))
    {
      secondary.push_back(secondary_el->volume());
      // We may not have projected onto a primary face in which case we may not have created mortar
      // segments
      if (auto it = _secondary_elems_to_mortar_segments.find(secondary_el->id());
          it != _secondary_elems_to_mortar_segments.end())
        for (const auto * const msm_elem : it->second)
        {
          msm.push_back(msm_elem->volume());
          const auto & msm_info = libmesh_map_find(_msm_elem_to_info, msm_elem);
          // Now it's also possible that we didn't project onto a primary face and we *did* create
          // mortar segments
          if (msm_info.primary_elem)
          {
            if (msm_info.primary_elem->subdomain_id() != primary_subd_id)
              mooseError("Unhandled primary-secondary pairing when computing mortar segment "
                         "statistics. This could happen if you have the same secondary "
                         "lower-dimensional subdomain ID paired with multiple lower-dimensional "
                         "primary subdomain IDs. Contact a MOOSE developer for help.");
            if (const auto [it, inserted] =
                    primary_elems_to_volume.emplace(msm_info.primary_elem->id(), Real{});
                inserted)
              it->second = msm_info.primary_elem->volume();
            else
              mooseAssert(
                  MooseUtils::absoluteFuzzyEqual(it->second, msm_info.primary_elem->volume()),
                  "Volumes should be consistent");
          }
        }
    }

    _mesh.comm().set_union(primary_elems_to_volume);
    _mesh.comm().allgather(static_cast<std::vector<Real> &>(secondary));
    _mesh.comm().allgather(static_cast<std::vector<Real> &>(msm));
    primary.reserve(primary_elems_to_volume.size());
    for (const auto [_, volume] : primary_elems_to_volume)
      primary.push_back(volume);

    MsmSubdomainStats stats;
    stats.primary_subd_id = primary_subd_id;
    stats.secondary_subd_id = secondary_subd_id;
    stats.secondary_lower_n_elems = secondary.size();
    stats.secondary_lower_max_volume = secondary.maximum();
    stats.secondary_lower_min_volume = secondary.minimum();
    stats.secondary_lower_median_volume = secondary.median();
    stats.primary_lower_n_elems = primary.size();
    stats.primary_lower_max_volume = primary.maximum();
    stats.primary_lower_min_volume = primary.minimum();
    stats.primary_lower_median_volume = primary.median();
    stats.msm_n_elems = msm.size();
    stats.msm_max_volume = msm.maximum();
    stats.msm_min_volume = msm.minimum();
    stats.msm_median_volume = msm.median();
    result.push_back(stats);

    primary.clear();
    secondary.clear();
    msm.clear();
    primary_elems_to_volume.clear();
  }

  return result;
}

void
AutomaticMortarGeneration::msmStatistics()
{
  const auto all_stats = computeMsmStatistics();

  if (_mesh.processor_id() != 0)
    return;

  Moose::out << "Mortar Interface Statistics:" << std::endl;
  for (const auto & stats : all_stats)
  {
    std::vector<std::string> col_names = {"mesh", "n_elems", "max", "min", "median"};
    std::vector<std::string> subds = {"secondary_lower", "primary_lower", "mortar_segment"};
    std::vector<size_t> n_elems = {
        stats.secondary_lower_n_elems, stats.primary_lower_n_elems, stats.msm_n_elems};
    std::vector<Real> maxs = {
        stats.secondary_lower_max_volume, stats.primary_lower_max_volume, stats.msm_max_volume};
    std::vector<Real> mins = {
        stats.secondary_lower_min_volume, stats.primary_lower_min_volume, stats.msm_min_volume};
    std::vector<Real> medians = {stats.secondary_lower_median_volume,
                                 stats.primary_lower_median_volume,
                                 stats.msm_median_volume};

    FormattedTable table;
    table.clear();
    for (auto i : index_range(subds))
    {
      table.addRow(i);
      table.addData<std::string>(col_names[0], subds[i]);
      table.addData<size_t>(col_names[1], n_elems[i]);
      table.addData<Real>(col_names[2], maxs[i]);
      table.addData<Real>(col_names[3], mins[i]);
      table.addData<Real>(col_names[4], medians[i]);
    }

    Moose::out << "secondary subdomain: " << stats.secondary_subd_id
               << " \tprimary subdomain: " << stats.primary_subd_id << std::endl;
    table.printTable(Moose::out, subds.size());
  }
}

// The blocks marked with **** are for regressing edge dropping treatment and should be removed
// eventually.
//****
// Compute inactve nodes when the old (incorrect) edge dropping treatemnt is enabled
void
AutomaticMortarGeneration::computeIncorrectEdgeDroppingInactiveLMNodes()
{
  using std::abs;

  // Note that in 3D our trick to check whether an element has edge dropping needs loose tolerances
  // since the mortar segments are on the linearized element and comparing the volume of the
  // linearized element does not have the same volume as the warped element
  const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;

  std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_inactive_nodes_set;
  const auto my_pid = _mesh.processor_id();

  // List of inactive nodes on local secondary elements
  std::unordered_set<dof_id_type> inactive_node_ids;

  std::unordered_map<const Elem *, Real> active_volume{};

  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    for (const auto el : _mesh.active_subdomain_elements_ptr_range(pr.second))
      active_volume[el] = 0.;

  // Compute fraction of elements with corresponding primary elements
  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
  {
    const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
    const Elem * secondary_elem = msinfo.secondary_elem;

    active_volume[secondary_elem] += msm_elem->volume();
  }

  // Mark all inactive local nodes
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    // Loop through all elements on my processor
    for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(pr.second))
      // If elem fully or partially dropped
      if (abs(active_volume[el] / el->volume() - 1.0) > tol)
      {
        // Add all nodes to list of inactive
        for (auto n : make_range(el->n_nodes()))
          inactive_node_ids.insert(el->node_id(n));
      }

  // Assemble list of procs that nodes contribute to
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
  {
    const auto secondary_subd_id = pr.second;

    // Loop through all elements not on my processor
    for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
    {
      // Get processor_id
      const auto pid = el->processor_id();

      // If element is in my subdomain, skip
      if (pid == my_pid)
        continue;

      // If element on proc pid shares any of my inactive nodes, mark to send
      for (const auto n : make_range(el->n_nodes()))
      {
        const auto node_id = el->node_id(n);
        if (inactive_node_ids.find(node_id) != inactive_node_ids.end())
          proc_to_inactive_nodes_set[pid].insert(node_id);
      }
    }
  }

  // Send list of inactive nodes
  {
    // Pack set into vector for sending (push_parallel_vector_data doesn't like sets)
    std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_inactive_nodes_vector;
    for (const auto & proc_set : proc_to_inactive_nodes_set)
      proc_to_inactive_nodes_vector[proc_set.first].insert(
          proc_to_inactive_nodes_vector[proc_set.first].end(),
          proc_set.second.begin(),
          proc_set.second.end());

    // First push data
    auto action_functor = [this, &inactive_node_ids](const processor_id_type pid,
                                                     const std::vector<dof_id_type> & sent_data)
    {
      if (pid == _mesh.processor_id())
        mooseError("Should not be communicating with self.");
      for (const auto pr : sent_data)
        inactive_node_ids.insert(pr);
    };
    TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_inactive_nodes_vector, action_functor);
  }
  _inactive_local_lm_nodes.clear();
  for (const auto node_id : inactive_node_ids)
    _inactive_local_lm_nodes.insert(_mesh.node_ptr(node_id));
}

void
AutomaticMortarGeneration::computeInactiveLMNodes()
{
  if (!_correct_edge_dropping)
  {
    computeIncorrectEdgeDroppingInactiveLMNodes();
    return;
  }

  std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_active_nodes_set;
  const auto my_pid = _mesh.processor_id();

  // List of active nodes on local secondary elements
  std::unordered_set<dof_id_type> active_local_nodes;

  // Mark all active local nodes
  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
  {
    const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
    const Elem * secondary_elem = msinfo.secondary_elem;

    for (auto n : make_range(secondary_elem->n_nodes()))
      active_local_nodes.insert(secondary_elem->node_id(n));
  }

  // Assemble list of procs that nodes contribute to
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
  {
    const auto secondary_subd_id = pr.second;

    // Loop through all elements not on my processor
    for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
    {
      // Get processor_id
      const auto pid = el->processor_id();

      // If element is in my subdomain, skip
      if (pid == my_pid)
        continue;

      // If element on proc pid shares any of my active nodes, mark to send
      for (const auto n : make_range(el->n_nodes()))
      {
        const auto node_id = el->node_id(n);
        if (active_local_nodes.find(node_id) != active_local_nodes.end())
          proc_to_active_nodes_set[pid].insert(node_id);
      }
    }
  }

  // Send list of active nodes
  {
    // Pack set into vector for sending (push_parallel_vector_data doesn't like sets)
    std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_active_nodes_vector;
    for (const auto & proc_set : proc_to_active_nodes_set)
    {
      proc_to_active_nodes_vector[proc_set.first].reserve(proc_to_active_nodes_set.size());
      for (const auto node_id : proc_set.second)
        proc_to_active_nodes_vector[proc_set.first].push_back(node_id);
    }

    // First push data
    auto action_functor = [this, &active_local_nodes](const processor_id_type pid,
                                                      const std::vector<dof_id_type> & sent_data)
    {
      if (pid == _mesh.processor_id())
        mooseError("Should not be communicating with self.");
      active_local_nodes.insert(sent_data.begin(), sent_data.end());
    };
    TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_active_nodes_vector, action_functor);
  }

  // Every proc has correct list of active local nodes, now take complement (list of inactive nodes)
  // and store to use later to zero LM DoFs on inactive nodes
  _inactive_local_lm_nodes.clear();
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
             /*secondary_subd_id*/ pr.second))
      for (const auto n : make_range(el->n_nodes()))
        if (active_local_nodes.find(el->node_id(n)) == active_local_nodes.end())
          _inactive_local_lm_nodes.insert(el->node_ptr(n));
}

// Note: could be combined with previous routine, keeping separate for clarity (for now)
void
AutomaticMortarGeneration::computeInactiveLMElems()
{
  // Mark all active secondary elements
  std::unordered_set<const Elem *> active_local_elems;

  //****
  // Note that in 3D our trick to check whether an element has edge dropping needs loose tolerances
  // since the mortar segments are on the linearized element and comparing the volume of the
  // linearized element does not have the same volume as the warped element
  const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;

  std::unordered_map<const Elem *, Real> active_volume;

  // Compute fraction of elements with corresponding primary elements
  if (!_correct_edge_dropping)
    for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
    {
      const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
      const Elem * secondary_elem = msinfo.secondary_elem;

      active_volume[secondary_elem] += msm_elem->volume();
    }
  //****

  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
  {
    const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
    const Elem * secondary_elem = msinfo.secondary_elem;

    //****
    if (!_correct_edge_dropping)
      if (abs(active_volume[secondary_elem] / secondary_elem->volume() - 1.0) > tol)
        continue;
    //****

    active_local_elems.insert(secondary_elem);
  }

  // Take complement of active elements in active local subdomain to get inactive local elements
  _inactive_local_lm_elems.clear();
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
             /*secondary_subd_id*/ pr.second))
      if (active_local_elems.find(el) == active_local_elems.end())
        _inactive_local_lm_elems.insert(el);
}

void
AutomaticMortarGeneration::computeNodalGeometry()
{
  // The dimension according to Mesh::mesh_dimension().
  const auto dim = _mesh.mesh_dimension();

  mooseAssert(dim == 2 || dim == 3,
              "AutomaticMortarGeneration::computeNodalGeometry() is only valid for "
              "mortar constraints on 2D or 3D meshes.");
  // A nodal lower-dimensional nodal quadrature rule to be used on faces.
  QNodal qface(dim - 1);

  // A map from the node id to the attached elemental normals/weights evaluated at the node. Th
  // length of the vector will correspond to the number of elements attached to the node. If it is a
  // vertex node, for a 1D mortar mesh, the vector length will be two. If it is an interior node,
  // the vector will be length 1. The first member of the pair is that element's normal at the node.
  // The second member is that element's JxW at the node
  std::map<dof_id_type, std::vector<std::pair<Point, Real>>> node_to_normals_map;

  Real sign = _periodic ? -1 : 1;

  // 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
  // parallelized as well.
  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
                                        end_el = _mesh.active_elements_end();
       el != end_el;
       ++el)
  {
    const Elem * secondary_elem = *el;

    // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
    if (!_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
      continue;

    // We will create an FE object and attach the nodal quadrature rule such that we can get out the
    // normals at the element nodes
    FEType nnx_fe_type(secondary_elem->default_order(), LAGRANGE);
    std::unique_ptr<FEBase> nnx_fe_face(FEBase::build(dim, nnx_fe_type));
    nnx_fe_face->attach_quadrature_rule(&qface);
    const auto & face_normals = nnx_fe_face->get_normals();
    const auto & face_points = nnx_fe_face->get_xyz();

    const auto & JxW = nnx_fe_face->get_JxW();

    // Which side of the parent are we? We need to know this to know
    // which side to reinit.
    const Elem * 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.");

    // Map to get lower dimensional element from interior parent on secondary surface
    // This map can be used to provide a handle to methods in this class that need to
    // operate on lower dimensional elements.
    _secondary_element_to_secondary_lowerd_element.emplace(interior_parent->id(), secondary_elem);

    // Look up which side of the interior parent secondary_elem is.
    auto s = interior_parent->which_side_am_i(secondary_elem);

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

    // Match by physical location instead of assuming that parent-side nodal
    // quadrature ordering and lower-dimensional side-element node ordering are
    // identical.
    const auto qpoint_to_secondary_node =
        nodalQuadraturePointToSecondaryNodeMap(*secondary_elem, face_points);

    mooseAssert(face_normals.size() == face_points.size() && JxW.size() == face_points.size(),
                "Face nodal geometry vectors must have the same size.");

    for (const auto qp : make_range(face_points.size()))
    {
      const auto n = qpoint_to_secondary_node[qp];
      auto & normals_and_weights_vec = node_to_normals_map[secondary_elem->node_id(n)];
      normals_and_weights_vec.push_back(std::make_pair(sign * face_normals[qp], JxW[qp]));
    }
  }

  for (const auto & pr : node_to_normals_map)
  {
    // Compute normal vector
    const auto & node_id = pr.first;
    const auto & normals_and_weights_vec = pr.second;

    Point nodal_normal;
    for (const auto & norm_and_weight : normals_and_weights_vec)
      nodal_normal += norm_and_weight.first * norm_and_weight.second;
    nodal_normal = nodal_normal.unit();

    _secondary_node_to_nodal_normal[_mesh.node_ptr(node_id)] = nodal_normal;

    Point nodal_tangent_one;
    Point nodal_tangent_two;
    householderOrthogolization(nodal_normal, nodal_tangent_one, nodal_tangent_two);

    _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][0] = nodal_tangent_one;
    _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][1] = nodal_tangent_two;
  }
}

void
AutomaticMortarGeneration::householderOrthogolization(const Point & nodal_normal,
                                                      Point & nodal_tangent_one,
                                                      Point & nodal_tangent_two) const
{
  using std::abs;

  mooseAssert(MooseUtils::absoluteFuzzyEqual(nodal_normal.norm(), 1),
              "The input nodal normal should have unity norm");

  const Real nx = nodal_normal(0);
  const Real ny = nodal_normal(1);
  const Real nz = nodal_normal(2);

  // See Lopes DS, Silva MT, Ambrosio JA. Tangent vectors to a 3-D surface normal: A geometric tool
  // to find orthogonal vectors based on the Householder transformation. Computer-Aided Design. 2013
  // Mar 1;45(3):683-94. We choose one definition of h_vector and deal with special case.
  const Point h_vector(nx + 1.0, ny, nz);

  // Avoid singularity of the equations at the end of routine by providing the solution to
  // (nx,ny,nz)=(-1,0,0) Normal/tangent fields can be visualized by outputting nodal geometry mesh
  // on a spherical problem.
  if (abs(h_vector(0)) < TOLERANCE)
  {
    nodal_tangent_one(0) = 0;
    nodal_tangent_one(1) = 1;
    nodal_tangent_one(2) = 0;

    nodal_tangent_two(0) = 0;
    nodal_tangent_two(1) = 0;
    nodal_tangent_two(2) = -1;

    return;
  }

  const Real h = h_vector.norm();

  nodal_tangent_one(0) = -2.0 * h_vector(0) * h_vector(1) / (h * h);
  nodal_tangent_one(1) = 1.0 - 2.0 * h_vector(1) * h_vector(1) / (h * h);
  nodal_tangent_one(2) = -2.0 * h_vector(1) * h_vector(2) / (h * h);

  nodal_tangent_two(0) = -2.0 * h_vector(0) * h_vector(2) / (h * h);
  nodal_tangent_two(1) = -2.0 * h_vector(1) * h_vector(2) / (h * h);
  nodal_tangent_two(2) = 1.0 - 2.0 * h_vector(2) * h_vector(2) / (h * h);
}

// Project secondary nodes onto their corresponding primary elements for each primary/secondary
// pair.
void
AutomaticMortarGeneration::projectSecondaryNodes()
{
  // For each primary/secondary boundary id pair, call the
  // project_secondary_nodes_single_pair() helper function.
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    projectSecondaryNodesSinglePair(pr.first, pr.second);
}

bool
AutomaticMortarGeneration::processAlignedNodes(
    const Node & secondary_node,
    const Node & primary_node,
    const std::vector<const Elem *> * secondary_node_neighbors,
    const std::vector<const Elem *> * primary_node_neighbors,
    const VectorValue<Real> & nodal_normal,
    const Elem & candidate_element,
    std::set<const Elem *> & rejected_elem_candidates)
{
  if (!secondary_node_neighbors)
    secondary_node_neighbors = &libmesh_map_find(_nodes_to_secondary_elem_map, secondary_node.id());
  if (!primary_node_neighbors)
    primary_node_neighbors = &libmesh_map_find(_nodes_to_primary_elem_map, primary_node.id());

  std::vector<bool> primary_elems_mapped(primary_node_neighbors->size(), false);

  // Add entries to secondary_node_and_elem_to_xi2_primary_elem container.
  //
  // First, determine "on left" vs. "on right" orientation of the nodal neighbors.
  // There can be a max of 2 nodal neighbors, and we want to make sure that the
  // secondary nodal neighbor on the "left" is associated with the primary nodal
  // neighbor on the "left" and similarly for the "right". We use cross products to determine
  // alignment. In the below diagram, 'x' denotes a node, and connected '|' are lower dimensional
  // elements.
  //                   x
  //           x       |
  //           |       |
  // secondary x ----> x primary
  //           |       |
  //           |       x
  //           x
  //
  //  Looking at the aligned nodes, the secondary node first, if we pick the top secondary lower
  //  dimensional element, then the cross product as written a few lines below points out of the
  //  screen towards you. (Point in the direction of the secondary nodal normal, and then curl your
  //  hand towards the secondary element's opposite node, then the thumb points in the direction of
  //  the cross product). Doing the same with the aligned primary node, if we pick the top primary
  //  element, then the cross product also points out of the screen. Because the cross products
  //  point in the same direction (positive dot product), then we know to associate the
  //  secondary-primary element pair. If we had picked the bottom primary element whose cross
  //  product points into the screen, then clearly the cross products point in the opposite
  //  direction and we don't have a match
  std::array<Real, 2> secondary_node_neighbor_cps, primary_node_neighbor_cps;

  for (const auto nn : index_range(*secondary_node_neighbors))
  {
    const Elem * const secondary_neigh = (*secondary_node_neighbors)[nn];
    const Point opposite = (secondary_neigh->node_ptr(0) == &secondary_node)
                               ? secondary_neigh->point(1)
                               : secondary_neigh->point(0);
    const Point cp = nodal_normal.cross(opposite - secondary_node);
    secondary_node_neighbor_cps[nn] = cp(2);
  }

  for (const auto nn : index_range(*primary_node_neighbors))
  {
    const Elem * const primary_neigh = (*primary_node_neighbors)[nn];
    const Point opposite = (primary_neigh->node_ptr(0) == &primary_node) ? primary_neigh->point(1)
                                                                         : primary_neigh->point(0);
    const Point cp = nodal_normal.cross(opposite - primary_node);
    primary_node_neighbor_cps[nn] = cp(2);
  }

  // Associate secondary/primary elems on matching sides.
  bool found_match = false;
  for (const auto snn : index_range(*secondary_node_neighbors))
    for (const auto mnn : index_range(*primary_node_neighbors))
      if (secondary_node_neighbor_cps[snn] * primary_node_neighbor_cps[mnn] > 0)
      {
        found_match = true;
        if (primary_elems_mapped[mnn])
          continue;
        primary_elems_mapped[mnn] = true;

        // Figure out xi^(2) value by looking at which node primary_node is
        // of the current primary node neighbor.
        const Real xi2 = (&primary_node == (*primary_node_neighbors)[mnn]->node_ptr(0)) ? -1 : +1;
        const auto secondary_key =
            std::make_pair(&secondary_node, (*secondary_node_neighbors)[snn]);
        const auto primary_val = std::make_pair(xi2, (*primary_node_neighbors)[mnn]);
        _secondary_node_and_elem_to_xi2_primary_elem.emplace(secondary_key, primary_val);

        // Also map in the other direction.
        const Real xi1 =
            (&secondary_node == (*secondary_node_neighbors)[snn]->node_ptr(0)) ? -1 : +1;

        const auto primary_key =
            std::make_tuple(primary_node.id(), &primary_node, (*primary_node_neighbors)[mnn]);
        const auto secondary_val = std::make_pair(xi1, (*secondary_node_neighbors)[snn]);
        _primary_node_and_elem_to_xi1_secondary_elem.emplace(primary_key, secondary_val);
      }

  if (!found_match)
  {
    // There could be coincident nodes and this might be a bad primary candidate (see
    // issue #21680). Instead of giving up, let's try continuing
    rejected_elem_candidates.insert(&candidate_element);
    return false;
  }

  // We need to handle the case where we've exactly projected a secondary node onto a
  // primary node, but our secondary node is at one of the secondary boundary face endpoints and
  // our primary node is not.
  if (secondary_node_neighbors->size() == 1 && primary_node_neighbors->size() == 2)
    for (const auto i : index_range(primary_elems_mapped))
      if (!primary_elems_mapped[i])
      {
        _primary_node_and_elem_to_xi1_secondary_elem.emplace(
            std::make_tuple(primary_node.id(), &primary_node, (*primary_node_neighbors)[i]),
            std::make_pair(1, nullptr));
      }

  return found_match;
}

void
AutomaticMortarGeneration::projectSecondaryNodesSinglePair(
    SubdomainID lower_dimensional_primary_subdomain_id,
    SubdomainID lower_dimensional_secondary_subdomain_id)
{
  using std::abs;

  // Build the "subdomain" adaptor based KD Tree.
  NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_primary_subdomain_id);
  subdomain_kd_tree_t kd_tree(
      3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));

  // Construct the KD tree.
  kd_tree.buildIndex();

  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
                                        end_el = _mesh.active_elements_end();
       el != end_el;
       ++el)
  {
    const Elem * secondary_side_elem = *el;

    // If this Elem is not in the current secondary subdomain, go on to the next one.
    if (secondary_side_elem->subdomain_id() != lower_dimensional_secondary_subdomain_id)
      continue;

    // For each node on the lower-dimensional element, find the nearest
    // node on the primary side using the KDTree, then
    // search in nearby elements for where it projects
    // along the nodal normal direction.
    for (MooseIndex(secondary_side_elem->n_vertices()) n = 0; n < secondary_side_elem->n_vertices();
         ++n)
    {
      const Node * secondary_node = secondary_side_elem->node_ptr(n);

      // Get the nodal neighbors for secondary_node, so we can check whether we've
      // already successfully projected it.
      const std::vector<const Elem *> & secondary_node_neighbors =
          this->_nodes_to_secondary_elem_map.at(secondary_node->id());

      // Check whether we've already mapped this secondary node
      // successfully for all of its nodal neighbors.
      bool is_mapped = true;
      for (MooseIndex(secondary_node_neighbors) snn = 0; snn < secondary_node_neighbors.size();
           ++snn)
      {
        auto secondary_key = std::make_pair(secondary_node, secondary_node_neighbors[snn]);
        if (!_secondary_node_and_elem_to_xi2_primary_elem.count(secondary_key))
        {
          is_mapped = false;
          break;
        }
      }

      // Go to the next node if this one has already been mapped.
      if (is_mapped)
        continue;

      // Look up the new nodal normal value in the local storage, error if not found.
      Point nodal_normal = _secondary_node_to_nodal_normal.at(secondary_node);

      // Data structure for performing Nanoflann searches.
      std::array<Real, 3> query_pt = {
          {(*secondary_node)(0), (*secondary_node)(1), (*secondary_node)(2)}};

      // The number of results we want to get.  We'll look for a
      // "few" nearest nodes, hopefully that is enough to let us
      // figure out which lower-dimensional Elem on the primary
      // side we are across from.
      const std::size_t num_results = 3;

      // Initialize result_set and do the search.
      std::vector<size_t> ret_index(num_results);
      std::vector<Real> out_dist_sqr(num_results);
      nanoflann::KNNResultSet<Real> result_set(num_results);
      result_set.init(&ret_index[0], &out_dist_sqr[0]);
      kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());

      // If this flag gets set in the loop below, we can break out of the outer r-loop as well.
      bool projection_succeeded = false;

      // Once we've rejected a candidate for a given secondary_node,
      // there's no reason to check it again.
      std::set<const Elem *> rejected_primary_elem_candidates;

      // Loop over the closest nodes, check whether
      // the secondary node successfully projects into
      // either of the closest neighbors, stop when
      // the projection succeeds.
      for (MooseIndex(result_set) r = 0; r < result_set.size(); ++r)
      {
        // Verify that the squared distance we compute is the same as nanoflann'sFss
        mooseAssert(abs((_mesh.point(ret_index[r]) - *secondary_node).norm_sq() -
                        out_dist_sqr[r]) <= TOLERANCE,
                    "Lower-dimensional element squared distance verification failed.");

        // Get a reference to the vector of lower dimensional elements from the
        // nodes_to_primary_elem_map.
        std::vector<const Elem *> & primary_elem_candidates =
            this->_nodes_to_primary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));

        // Search the Elems connected to this node on the primary mesh side.
        for (MooseIndex(primary_elem_candidates) e = 0; e < primary_elem_candidates.size(); ++e)
        {
          const Elem * primary_elem_candidate = primary_elem_candidates[e];

          // If we've already rejected this candidate, we don't need to check it again.
          if (rejected_primary_elem_candidates.count(primary_elem_candidate))
            continue;

          // Now generically solve for xi2
          const auto order = primary_elem_candidate->default_order();
          DualNumber<Real> xi2_dn{0, 1};
          unsigned int current_iterate = 0, max_iterates = 10;

          // Newton loop
          do
          {
            VectorValue<DualNumber<Real>> x2(0);
            for (MooseIndex(primary_elem_candidate->n_nodes()) n = 0;
                 n < primary_elem_candidate->n_nodes();
                 ++n)
              x2 +=
                  Moose::fe_lagrange_1D_shape(order, n, xi2_dn) * primary_elem_candidate->point(n);
            const auto u = x2 - (*secondary_node);
            const auto F = u(0) * nodal_normal(1) - u(1) * nodal_normal(0);

            if (abs(F) < _newton_tolerance)
              break;

            if (F.derivatives())
            {
              Real dxi2 = -F.value() / F.derivatives();

              xi2_dn += dxi2;
            }
            else
              // It's possible that the secondary surface nodal normal is completely orthogonal to
              // the primary surface normal, in which case the derivative is 0. We know in this case
              // that the projection should be a failure
              current_iterate = max_iterates;
          } while (++current_iterate < max_iterates);

          Real xi2 = xi2_dn.value();

          // Check whether the projection worked. The last condition checks for obliqueness of the
          // projection
          //
          // We are projecting on one side first and the other side second. If we make the
          // tolerance bigger and remove the (5) factor we are going to continue to miss the
          // second projection and fall into the exception message in
          // projectPrimaryNodesSinglePair. What makes this modification to not fall in the
          // exception is that we are projecting on one side more xi than in the other. There
          // should be a better way of doing this by using actual distances and not parametric
          // coordinates. But I believe making the tolerance uniformly larger or smaller won't do
          // the trick here.
          if ((current_iterate < max_iterates) && (std::abs(xi2) <= 1. + 5 * _xi_tolerance) &&
              (abs((primary_elem_candidate->point(0) - primary_elem_candidate->point(1)).unit() *
                   nodal_normal) < std::cos(_minimum_projection_angle * libMesh::pi / 180)))
          {
            // If xi2 == +1 or -1 then this secondary node mapped directly to a node on the primary
            // surface. This isn't as unlikely as you might think, it will happen if the meshes
            // on the interface start off being perfectly aligned. In this situation, we need to
            // associate the secondary node with two different elements (and two corresponding
            // xi^(2) values.
            if (abs(abs(xi2) - 1.) <= _xi_tolerance * 5.0)
            {
              const Node * primary_node = (xi2 < 0) ? primary_elem_candidate->node_ptr(0)
                                                    : primary_elem_candidate->node_ptr(1);
              const bool created_mortar_segment =
                  processAlignedNodes(*secondary_node,
                                      *primary_node,
                                      &secondary_node_neighbors,
                                      nullptr,
                                      nodal_normal,
                                      *primary_elem_candidate,
                                      rejected_primary_elem_candidates);

              if (!created_mortar_segment)
                continue;
            }
            else // Point falls somewhere in the middle of the Elem.
            {
              // Add two entries to secondary_node_and_elem_to_xi2_primary_elem.
              for (MooseIndex(secondary_node_neighbors) nn = 0;
                   nn < secondary_node_neighbors.size();
                   ++nn)
              {
                const Elem * neigh = secondary_node_neighbors[nn];
                for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
                {
                  const Node * neigh_node = neigh->node_ptr(nid);
                  if (secondary_node == neigh_node)
                  {
                    auto key = std::make_pair(neigh_node, neigh);
                    auto val = std::make_pair(xi2, primary_elem_candidate);
                    _secondary_node_and_elem_to_xi2_primary_elem.emplace(key, val);
                  }
                }
              }
            }

            projection_succeeded = true;
            break; // out of e-loop
          }
          else
            // The current secondary_node is not in this Elem, so keep track of the rejects.
            rejected_primary_elem_candidates.insert(primary_elem_candidate);
        }

        if (projection_succeeded)
          break; // out of r-loop
      } // r-loop

      if (!projection_succeeded)
      {
        _failed_secondary_node_projections.insert(secondary_node->id());
        if (_debug)
          _console << "Failed to find primary Elem into which secondary node "
                   << static_cast<const Point &>(*secondary_node) << ", id '"
                   << secondary_node->id() << "', projects onto\n"
                   << std::endl;
      }
      else if (_debug)
        _projected_secondary_nodes.insert(secondary_node->id());
    } // loop over side nodes
  } // end loop over lower-dimensional elements

  if (_distributed)
  {
    if (_debug)
      _mesh.comm().set_union(_projected_secondary_nodes);
    _mesh.comm().set_union(_failed_secondary_node_projections);
  }

  if (_debug)
    _console << "\n"
             << _projected_secondary_nodes.size() << " out of "
             << _projected_secondary_nodes.size() + _failed_secondary_node_projections.size()
             << " secondary nodes were successfully projected\n"
             << std::endl;
}

// Inverse map primary nodes onto their corresponding secondary elements for each primary/secondary
// pair.
void
AutomaticMortarGeneration::projectPrimaryNodes()
{
  // For each primary/secondary boundary id pair, call the
  // project_primary_nodes_single_pair() helper function.
  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
    projectPrimaryNodesSinglePair(pr.first, pr.second);
}

void
AutomaticMortarGeneration::projectPrimaryNodesSinglePair(
    SubdomainID lower_dimensional_primary_subdomain_id,
    SubdomainID lower_dimensional_secondary_subdomain_id)
{
  using std::abs;

  // Build a Nanoflann object on the lower-dimensional secondary elements of the Mesh.
  NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_secondary_subdomain_id);
  subdomain_kd_tree_t kd_tree(
      3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));

  // Construct the KD tree for lower-dimensional elements in the volume mesh.
  kd_tree.buildIndex();

  std::unordered_set<dof_id_type> primary_nodes_visited;

  for (const auto & primary_side_elem : _mesh.active_element_ptr_range())
  {
    // If this is not one of the lower-dimensional primary side elements, go on to the next one.
    if (primary_side_elem->subdomain_id() != lower_dimensional_primary_subdomain_id)
      continue;

    // For each node on this side, find the nearest node on the secondary side using the KDTree,
    // then search in nearby elements for where it projects along the nodal normal direction.
    for (MooseIndex(primary_side_elem->n_vertices()) n = 0; n < primary_side_elem->n_vertices();
         ++n)
    {
      // Get a pointer to this node.
      const Node * primary_node = primary_side_elem->node_ptr(n);

      // Get the nodal neighbors connected to this primary node.
      const std::vector<const Elem *> & primary_node_neighbors =
          _nodes_to_primary_elem_map.at(primary_node->id());

      // Check whether we have already successfully inverse mapped this primary node (whether during
      // secondary node projection or now during primary node projection) or we have already failed
      // to inverse map this primary node (now during primary node projection), and then skip if
      // either of those things is true
      auto primary_key =
          std::make_tuple(primary_node->id(), primary_node, primary_node_neighbors[0]);
      if (!primary_nodes_visited.insert(primary_node->id()).second ||
          _primary_node_and_elem_to_xi1_secondary_elem.count(primary_key))
        continue;

      // Data structure for performing Nanoflann searches.
      Real query_pt[3] = {(*primary_node)(0), (*primary_node)(1), (*primary_node)(2)};

      // The number of results we want to get.  We'll look for a
      // "few" nearest nodes, hopefully that is enough to let us
      // figure out which lower-dimensional Elem on the secondary side
      // we are across from.
      const size_t num_results = 3;

      // Initialize result_set and do the search.
      std::vector<size_t> ret_index(num_results);
      std::vector<Real> out_dist_sqr(num_results);
      nanoflann::KNNResultSet<Real> result_set(num_results);
      result_set.init(&ret_index[0], &out_dist_sqr[0]);
      kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());

      // If this flag gets set in the loop below, we can break out of the outer r-loop as well.
      bool projection_succeeded = false;

      // Once we've rejected a candidate for a given
      // primary_node, there's no reason to check it
      // again.
      std::set<const Elem *> rejected_secondary_elem_candidates;

      // Loop over the closest nodes, check whether the secondary node successfully projects into
      // either of the closest neighbors, stop when the projection succeeds.
      for (MooseIndex(result_set) r = 0; r < result_set.size(); ++r)
      {
        // Verify that the squared distance we compute is the same as nanoflann's
        mooseAssert(abs((_mesh.point(ret_index[r]) - *primary_node).norm_sq() - out_dist_sqr[r]) <=
                        TOLERANCE,
                    "Lower-dimensional element squared distance verification failed.");

        // Get a reference to the vector of lower dimensional elements from the
        // nodes_to_secondary_elem_map.
        const std::vector<const Elem *> & secondary_elem_candidates =
            _nodes_to_secondary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));

        // Print the Elems connected to this node on the secondary mesh side.
        for (MooseIndex(secondary_elem_candidates) e = 0; e < secondary_elem_candidates.size(); ++e)
        {
          const Elem * secondary_elem_candidate = secondary_elem_candidates[e];

          // If we've already rejected this candidate, we don't need to check it again.
          if (rejected_secondary_elem_candidates.count(secondary_elem_candidate))
            continue;

          std::vector<Point> nodal_normals(secondary_elem_candidate->n_nodes());
          for (const auto n : make_range(secondary_elem_candidate->n_nodes()))
            nodal_normals[n] =
                _secondary_node_to_nodal_normal.at(secondary_elem_candidate->node_ptr(n));

          // Use equation 2.4.6 from Bin Yang's dissertation to try and solve for
          // the position on the secondary element where this primary came from.  This
          // requires a Newton iteration in general.
          DualNumber<Real> xi1_dn{0, 1}; // initial guess
          auto && order = secondary_elem_candidate->default_order();
          unsigned int current_iterate = 0, max_iterates = 10;

          VectorValue<DualNumber<Real>> normals(0);

          // Newton iteration loop - this to converge in 1 iteration when it
          // succeeds, and possibly two iterations when it converges to a
          // xi outside the reference element. I don't know any reason why it should
          // only take 1 iteration -- the Jacobian is not constant in general...
          do
          {
            VectorValue<DualNumber<Real>> x1(0);
            for (MooseIndex(secondary_elem_candidate->n_nodes()) n = 0;
                 n < secondary_elem_candidate->n_nodes();
                 ++n)
            {
              const auto phi = Moose::fe_lagrange_1D_shape(order, n, xi1_dn);
              x1 += phi * secondary_elem_candidate->point(n);
              normals += phi * nodal_normals[n];
            }

            const auto u = x1 - (*primary_node);

            const auto F = u(0) * normals(1) - u(1) * normals(0);

            if (abs(F) < _newton_tolerance)
              break;

            // Unlike for projection of nodal normals onto primary surfaces, we should never have a
            // case where the nodal normal is completely orthogonal to the secondary surface, so we
            // do not have to guard against F.derivatives() == 0 here
            Real dxi1 = -F.value() / F.derivatives();

            xi1_dn += dxi1;

            normals = 0;
          } while (++current_iterate < max_iterates);

          Real xi1 = xi1_dn.value();

          // Check for convergence to a valid solution... The last condition checks for obliqueness
          // of the projection
          if ((current_iterate < max_iterates) && (abs(xi1) <= 1. + _xi_tolerance) &&
              (abs((primary_side_elem->point(0) - primary_side_elem->point(1)).unit() *
                   MetaPhysicL::raw_value(normals).unit()) <
               std::cos(_minimum_projection_angle * libMesh::pi / 180.0)))
          {
            if (abs(abs(xi1) - 1.) < _xi_tolerance)
            {
              // Special case: xi1=+/-1.
              // It is unlikely that we get here, because this primary node should already
              // have been mapped during the project_secondary_nodes() routine, but
              // there is still a chance since the tolerances are applied to
              // the xi coordinate and that value may be different on a primary element and a
              // secondary element since they may have different sizes. It's also possible that we
              // may reach this point if the solve has yielded a non-physical configuration such as
              // one block being pushed way out into space
              const Node & secondary_node = (xi1 < 0) ? secondary_elem_candidate->node_ref(0)
                                                      : secondary_elem_candidate->node_ref(1);
              bool created_mortar_segment = false;

              // If we have failed to project this secondary node, let's try again now
              if (_failed_secondary_node_projections.count(secondary_node.id()))
                created_mortar_segment = processAlignedNodes(secondary_node,
                                                             *primary_node,
                                                             nullptr,
                                                             &primary_node_neighbors,
                                                             MetaPhysicL::raw_value(normals),
                                                             *secondary_elem_candidate,
                                                             rejected_secondary_elem_candidates);
              else
                rejected_secondary_elem_candidates.insert(secondary_elem_candidate);

              if (!created_mortar_segment)
                // We used to throw an exception in this scope but now that we support processing
                // aligned nodes within this primary node projection method, I don't see any harm in
                // simply rejecting the secondary element candidate in the case of failure and
                // continuing just as we do when projecting secondary nodes
                continue;
            }
            else // somewhere in the middle of the Elem
            {
              // Add entry to primary_node_and_elem_to_xi1_secondary_elem
              //
              // Note: we originally duplicated the map values for the keys (node, left_neighbor)
              // and (node, right_neighbor) but I don't think that should be necessary. Instead we
              // just do it for neighbor 0, but really maybe we don't even need to do that since
              // we can always look up the neighbors later given the Node... keeping it like this
              // helps to maintain the "symmetry" of the two containers.
              const Elem * neigh = primary_node_neighbors[0];
              for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
              {
                const Node * neigh_node = neigh->node_ptr(nid);
                if (primary_node == neigh_node)
                {
                  auto key = std::make_tuple(neigh_node->id(), neigh_node, neigh);
                  auto val = std::make_pair(xi1, secondary_elem_candidate);
                  _primary_node_and_elem_to_xi1_secondary_elem.emplace(key, val);
                }
              }
            }

            projection_succeeded = true;
            break; // out of e-loop
          }
          else
          {
            // The current primary_point is not in this Elem, so keep track of the rejects.
            rejected_secondary_elem_candidates.insert(secondary_elem_candidate);
          }
        } // end e-loop over candidate elems

        if (projection_succeeded)
          break; // out of r-loop
      } // r-loop

      if (!projection_succeeded && _debug)
      {
        _console << "\nFailed to find point from which primary node "
                 << static_cast<const Point &>(*primary_node) << " was projected." << std::endl
                 << std::endl;
      }
    } // loop over side nodes
  } // end loop over elements for finding where primary points would have projected from.
}

std::vector<AutomaticMortarGeneration::MortarFilterIter>
AutomaticMortarGeneration::secondariesToMortarSegments(const Node & node) const
{
  auto secondary_it = _nodes_to_secondary_elem_map.find(node.id());
  if (secondary_it == _nodes_to_secondary_elem_map.end())
    return {};

  const auto & secondary_elems = secondary_it->second;
  std::vector<MortarFilterIter> ret;
  ret.reserve(secondary_elems.size());

  for (const auto i : index_range(secondary_elems))
  {
    auto * const secondary_elem = secondary_elems[i];
    auto msm_it = _secondary_elems_to_mortar_segments.find(secondary_elem->id());
    if (msm_it == _secondary_elems_to_mortar_segments.end())
      // We may have removed this element key from this map
      continue;

    mooseAssert(secondary_elem->active(),
                "We loop over active elements when building the mortar segment mesh, so we golly "
                "well hope this is active.");
    mooseAssert(!msm_it->second.empty(),
                "We should have removed all secondaries from this map if they do not have any "
                "mortar segments associated with them.");
    ret.push_back(msm_it);
  }

  return ret;
}