AutomaticMortarGeneration Class Reference

This class is a container/interface for the objects involved in automatic generation of mortar spaces. More...

#include <AutomaticMortarGeneration.h>

Inheritance diagram for AutomaticMortarGeneration:

Classes
struct	MsmSubdomainStats
	Statistics for one primary-secondary subdomain pair. More...

Public Types
using	MortarFilterIter = std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > >::const_iterator

Public Member Functions
	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=Mortar3DQuadraturePointMapping::NORMAL_PROJECTION)
	Must be constructed with a reference to the Mesh we are generating mortar spaces for.
void	buildNodeToElemMaps ()
	Once the secondary_requested_boundary_ids and primary_requested_boundary_ids containers have been filled in, call this function to build node-to-Elem maps for the lower-dimensional elements.
void	computeNodalGeometry ()
	Computes and stores the nodal normal/tangent vectors in a local data structure instead of using the ExplicitSystem/NumericVector approach.
void	projectSecondaryNodes ()
	Project secondary nodes (find xi^(2) values) to the closest points on the primary surface.
void	projectPrimaryNodes ()
	(Inverse) project primary nodes to the points on the secondary surface where they would have come from (find (xi^(1) values)).
void	buildMortarSegmentMesh ()
	Builds the mortar segment mesh once the secondary and primary node projections have been completed.
void	buildMortarSegmentMesh3d ()
	Builds the mortar segment mesh once the secondary and primary node projections have been completed.
std::vector< MsmSubdomainStats >	computeMsmStatistics ()
	Computes mortar segment mesh statistics and returns one entry per subdomain pair.
void	msmStatistics ()
	Prints mortar segment mesh statistics to console (calls computeMsmStatistics internally)
void	clear ()
	Clears the mortar segment mesh and accompanying data structures.
void	meshChanged ()
	Invalidates the cached MSM node/element ID starting offset so that the next call to buildMortarSegmentMesh3d() recomputes it via allgather.
bool	onDisplaced () const
	returns whether this object is on the displaced mesh
std::vector< Point >	getNodalNormals (const Elem &secondary_elem) const
std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 >	getNodalTangents (const Elem &secondary_elem) const
	Compute the two nodal tangents, which are built on-the-fly.
void	computeADNodalNormals (const std::function< ADPoint(const Node &)> &coordinate, std::unordered_map< const Node *, ADRealVectorValue > &nodal_normals) const
	Build the normalized JxW-weighted secondary nodal normals from AD nodal coordinates.
std::map< unsigned int, unsigned int >	getSecondaryIpToLowerElementMap (const Elem &lower_secondary_elem) const
	Compute on-the-fly mapping from secondary interior parent nodes to lower dimensional nodes.
std::map< unsigned int, unsigned int >	getPrimaryIpToLowerElementMap (const Elem &primary_elem, const Elem &primary_elem_ip, const Elem &lower_secondary_elem) const
	Compute on-the-fly mapping from primary interior parent nodes to its corresponding lower dimensional nodes.
std::vector< Point >	getNormals (const Elem &secondary_elem, const std::vector< Point > &xi1_pts) const
	Compute the normals at given reference points on a secondary element.
std::vector< Point >	getNormals (const Elem &secondary_elem, const std::vector< Real > &oned_xi1_pts) const
	Compute the normals at given reference points on a secondary element.
const Elem *	getSecondaryLowerdElemFromSecondaryElem (dof_id_type secondary_elem_id) const
	Return lower dimensional secondary element given its interior parent.
void	computeInactiveLMNodes ()
	Get list of secondary nodes that don't contribute to interaction with any primary element.
void	computeIncorrectEdgeDroppingInactiveLMNodes ()
	Computes inactive secondary nodes when incorrect edge dropping behavior is enabled (any node touching a partially or fully dropped element is dropped)
void	computeInactiveLMElems ()
	Get list of secondary elems without any corresponding primary elements.
const std::unordered_map< dof_id_type, std::unordered_set< dof_id_type > > &	mortarInterfaceCoupling () const
void	requestNonlocalSecondaryCouplings ()
	Request that every process receive the mortar interface couplings of the nonlocal secondary lower-dimensional elements that point-neighbor ghosting reads for its local elements: the point neighbors of its local secondary elements and of the secondary elements coupled to its local elements.
const std::pair< BoundaryID, BoundaryID > &	primarySecondaryBoundaryIDPair () const
const MeshBase &	mortarSegmentMesh () const
const std::unordered_map< const Elem *, MortarSegmentInfo > &	mortarSegmentMeshElemToInfo () const
int	dim () const
Mortar3DQuadraturePointMapping	mortar3DQpMapping () const
	Return the 3D mortar quadrature-point mapping method.
const MortarSegmentReferencePoints &	mortarSegmentReferencePoints (const Elem &mortar_segment_elem) const
	Return the parent-face reference coordinates for a mortar segment.
const std::unordered_set< const Node * > &	getInactiveLMNodes () const
const std::unordered_set< const Elem * > &	getInactiveLMElems () const
bool	incorrectEdgeDropping () const
std::vector< MortarFilterIter >	secondariesToMortarSegments (const Node &node) const
const std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > > &	secondariesToMortarSegments () const
const std::set< SubdomainID > &	secondaryIPSubIDs () const
const std::set< SubdomainID > &	primaryIPSubIDs () const
const std::unordered_map< dof_id_type, std::vector< const Elem * > > &	nodesToSecondaryElem () const
void	initOutput ()
	initialize mortar-mesh based output

Public Attributes
const ConsoleStream	_console
	An instance of helper class to write streams to the Console objects.

Static Public Attributes
static const std::string	system_name
	The name of the nodal normals system.

Private Member Functions
void	outputMortarMesh ()
	Write the mortar segment mesh to exodus.
std::string	mortarInterfaceName () const
void	buildCouplingInformation ()
	build the _mortar_interface_coupling data
void	projectSecondaryNodesSinglePair (SubdomainID lower_dimensional_primary_subdomain_id, SubdomainID lower_dimensional_secondary_subdomain_id)
	Helper function responsible for projecting secondary nodes onto primary elements for a single primary/secondary pair.
void	projectPrimaryNodesSinglePair (SubdomainID lower_dimensional_primary_subdomain_id, SubdomainID lower_dimensional_secondary_subdomain_id)
	Helper function used internally by AutomaticMortarGeneration::project_primary_nodes().
bool	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_element_candidates)
	Process aligned nodes.

Private Attributes
MooseApp &	_app
	The Moose app.
MeshBase &	_mesh
	Reference to the mesh stored in equation_systems.
std::set< BoundaryID >	_secondary_requested_boundary_ids
	The boundary ids corresponding to all the secondary surfaces.
std::set< BoundaryID >	_primary_requested_boundary_ids
	The boundary ids corresponding to all the primary surfaces.
std::vector< std::pair< BoundaryID, BoundaryID > >	_primary_secondary_boundary_id_pairs
	A list of primary/secondary boundary id pairs corresponding to each side of the mortar interface.
std::unordered_map< dof_id_type, std::vector< const Elem * > >	_nodes_to_secondary_elem_map
	Map from nodes to connected lower-dimensional elements on the secondary/primary subdomains.
std::unordered_map< dof_id_type, std::vector< const Elem * > >	_nodes_to_primary_elem_map
std::unordered_map< std::pair< const Node *, const Elem * >, std::pair< Real, const Elem * > >	_secondary_node_and_elem_to_xi2_primary_elem
	Similar to the map above, but associates a (Secondary Node, Secondary Elem) pair to a (xi^(2), primary Elem) pair.
std::map< std::tuple< dof_id_type, const Node *, const Elem * >, std::pair< Real, const Elem * > >	_primary_node_and_elem_to_xi1_secondary_elem
	Same type of container, but for mapping (Primary Node ID, Primary Node, Primary Elem) -> (xi^(1), Secondary Elem) where they are inverse-projected along the nodal normal direction.
std::unique_ptr< MeshBase >	_mortar_segment_mesh
	1D Mesh of mortar segment elements which gets built by the call to build_mortar_segment_mesh().
std::unordered_map< const Elem *, MortarSegmentInfo >	_msm_elem_to_info
	Map between Elems in the mortar segment mesh and their info structs.
std::unordered_map< const Elem *, unsigned int >	_lower_elem_to_side_id
	Keeps track of the mapping between lower-dimensional elements and the side_id of the interior_parent which they are.
std::vector< std::pair< SubdomainID, SubdomainID > >	_primary_secondary_subdomain_id_pairs
	A list of primary/secondary subdomain id pairs corresponding to each side of the mortar interface.
std::set< SubdomainID >	_secondary_boundary_subdomain_ids
	The secondary/primary lower-dimensional boundary subdomain ids are the secondary/primary boundary ids.
std::set< SubdomainID >	_primary_boundary_subdomain_ids
std::unordered_map< dof_id_type, std::unordered_set< dof_id_type > >	_mortar_interface_coupling
	Used by the AugmentSparsityOnInterface functor to determine whether a given Elem is coupled to any others across the gap, and to explicitly set up the dependence between interior_parent() elements on the secondary side and their lower-dimensional sides which are on the interface.
std::unordered_map< const Node *, Point >	_secondary_node_to_nodal_normal
	Container for storing the nodal normal vector associated with each secondary node.
std::unordered_map< const Node *, std::array< Point, 2 > >	_secondary_node_to_hh_nodal_tangents
	Container for storing the nodal tangent/binormal vectors associated with each secondary node (Householder approach).
std::unordered_map< dof_id_type, const Elem * >	_secondary_element_to_secondary_lowerd_element
	Map from full dimensional secondary element id to lower dimensional secondary element.
std::unordered_set< const Node * >	_inactive_local_lm_nodes
	List of inactive lagrange multiplier nodes owned by this processor (for nodal variables)
std::unordered_set< const Elem * >	_inactive_local_lm_elems
	List of inactive lagrange multiplier nodes (for elemental variables)
std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > >	_secondary_elems_to_mortar_segments
	We maintain a mapping from lower-dimensional secondary elements in the original mesh to (sets of) elements in mortar_segment_mesh.
std::set< SubdomainID >	_secondary_ip_sub_ids
	All the secondary interior parent subdomain IDs associated with the mortar mesh.
std::set< SubdomainID >	_primary_ip_sub_ids
	All the primary interior parent subdomain IDs associated with the mortar mesh.
const bool	_debug
	Whether to print debug output.
const bool	_on_displaced
	Whether this object is on the displaced mesh.
const bool	_periodic
	Whether this object will be generating a mortar segment mesh for periodic constraints.
const bool	_distributed
	Whether the mortar segment mesh is distributed.
bool	_nonlocal_secondary_couplings = false
	Whether buildCouplingInformation() also gathers the couplings of nonlocal secondary lower-dimensional elements.
Real	_newton_tolerance = 1e-12
	Newton solve tolerance for node projections.
Real	_xi_tolerance = 1e-6
	Tolerance for checking projection xi values.
const bool	_correct_edge_dropping
	Flag to enable regressed treatment of edge dropping where all LM DoFs on edge dropping element are strongly set to 0.
const Real	_minimum_projection_angle
	Parameter to control which angle (in degrees) is admissible for the creation of mortar segments.
const Mortar3DSubpatchPlane	_mortar_3d_subpatch_plane
	Method used to define the local projection planes for 3D secondary subpatches.
const MortarSegmentTriangulationMode	_triangulation_mode
	Triangulation mode used for clipped 3D mortar polygons.
const bool	_triangulate_triangles
	Whether already-triangular clipped polygons should still be centroid-subdivided.
const Mortar3DQuadraturePointMapping	_mortar_3d_qp_mapping
	Method used to map 3D mortar segment quadrature points to primary and secondary faces.
std::unordered_map< const Elem *, MortarSegmentReferencePoints >	_msm_elem_to_reference_points
	Reference-coordinate data used only by the reference-interpolation mapping mode.
std::unique_ptr< InputParameters >	_output_params
	Storage for the input parameters used by the mortar nodal geometry output.
std::optional< dof_id_type >	_msm_node_id_start
	Cached per-rank starting ID for 3D MSM nodes/elements.
std::unordered_set< dof_id_type >	_projected_secondary_nodes
	Debugging container for printing information about fraction of successful projections for secondary nodes.
std::unordered_set< dof_id_type >	_failed_secondary_node_projections
	Secondary nodes that failed to project.

Friends
class	MortarNodalGeometryOutput
class	AugmentSparsityOnInterface

Detailed Description

This class is a container/interface for the objects involved in automatic generation of mortar spaces.

Definition at line 133 of file AutomaticMortarGeneration.h.

Member Typedef Documentation

MortarFilterIter

using AutomaticMortarGeneration::MortarFilterIter = std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID> >::const_iterator

Definition at line 442 of file AutomaticMortarGeneration.h.

Constructor & Destructor Documentation

AutomaticMortarGeneration()

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 = Mortar3DQuadraturePointMapping::NORMAL_PROJECTION
	)

Must be constructed with a reference to the Mesh we are generating mortar spaces for.

Definition at line 221 of file AutomaticMortarGeneration.C.

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

Member Function Documentation

buildCouplingInformation()

void AutomaticMortarGeneration::buildCouplingInformation ( )

private

build the _mortar_interface_coupling data

Definition at line 1548 of file AutomaticMortarGeneration.C.

{
  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;
 
    // Every process owning a row on a secondary node or either volume parent needs the complete
    // secondary face-star coupling, including when that process also owns the lower-dimensional
    // face.
    std::set<processor_id_type> secondary_face_coupling_processors{
        secondary_elem->processor_id(),
        secondary_elem->interior_parent()->processor_id(),
        primary_elem->interior_parent()->processor_id()};
    for (const auto n : secondary_elem->node_index_range())
      secondary_face_coupling_processors.insert(secondary_elem->node_ptr(n)->processor_id());
 
    // LowerSecondary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(),
                                            secondary_elem->interior_parent()->id());
 
    // LowerPrimary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(),
                                            primary_elem->interior_parent()->id());
 
    // Lower-LowerDimensionalPrimary
    for (const auto processor : secondary_face_coupling_processors)
      coupling_info[processor].emplace_back(secondary_elem->id(), primary_elem->id());
 
    // SecondaryLower
    coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(
        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());
  }
 
  // A face pair may produce several mortar segments. Remove repeated entries before communication.
  for (auto & [_, couplings] : coupling_info)
  {
    std::sort(couplings.begin(), couplings.end());
    couplings.erase(std::unique(couplings.begin(), couplings.end()), couplings.end());
  }
 
  // Push the coupling information
  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);
 
  if (!_nonlocal_secondary_couplings)
    return;
 
  // The push above sends a secondary lower-dimensional element's couplings only to the processes
  // owning it, its nodes, or its segments' volume parents. Some consumers (e.g. penalty contact)
  // apply a secondary node's constraint on every process that contributes to it, not only on the
  // node's owner. The weighted gap at a node shared by a local element E and a nonlocal element F
  // carries the derivatives of the primary elements coupled to F, and the process owning E (or a
  // volume parent coupled to E) applies that gap to its own rows. When AugmentSparsityOnInterface
  // ghosts point neighbors, it reads the couplings of F to add those primary columns to the rows.
  // For a local element, it reads the couplings of the point neighbors of the element's own
  // secondary lower-dimensional elements and of the secondary elements coupled to it (e.g. those
  // paired with a local primary element), so pull the couplings of those point neighbors from
  // their owners, which hold the complete sets after the push
  const auto is_secondary = [this](const Elem * const elem)
  { return _secondary_boundary_subdomain_ids.count(elem->subdomain_id()); };
  const auto pid = _mesh.processor_id();
  std::unordered_set<const Elem *> secondary_elems_read;
  for (const auto secondary_subdomain_id : _secondary_boundary_subdomain_ids)
    for (const Elem * const secondary_elem :
         as_range(_mesh.active_subdomain_elements_begin(secondary_subdomain_id),
                  _mesh.active_subdomain_elements_end(secondary_subdomain_id)))
      if (secondary_elem->processor_id() == pid ||
          secondary_elem->interior_parent()->processor_id() == pid)
        secondary_elems_read.insert(secondary_elem);
  for (const auto & [elem_id, coupled_elem_ids] : _mortar_interface_coupling)
    if (const Elem * const elem = _mesh.query_elem_ptr(elem_id);
        elem && elem->processor_id() == pid)
      for (const auto coupled_elem_id : coupled_elem_ids)
        if (const Elem * const coupled_elem = _mesh.query_elem_ptr(coupled_elem_id);
            coupled_elem && is_secondary(coupled_elem))
          secondary_elems_read.insert(coupled_elem);
 
  std::unordered_set<const Elem *> secondary_elems_queried;
  std::unordered_map<processor_id_type, std::vector<dof_id_type>> secondary_elems_to_query;
  for (const Elem * const secondary_elem : secondary_elems_read)
  {
    std::set<const Elem *> point_neighbors;
    secondary_elem->find_point_neighbors(point_neighbors);
    point_neighbors.insert(secondary_elem);
    for (const Elem * const neigh : point_neighbors)
      if (neigh->active() && is_secondary(neigh) && neigh->processor_id() != pid &&
          secondary_elems_queried.insert(neigh).second)
        secondary_elems_to_query[neigh->processor_id()].push_back(neigh->id());
  }
 
  auto gather_functor = [this](processor_id_type,
                               const std::vector<dof_id_type> & secondary_elem_ids,
                               std::vector<std::vector<dof_id_type>> & coupled_elem_ids)
  {
    coupled_elem_ids.resize(secondary_elem_ids.size());
    for (const auto i : index_range(secondary_elem_ids))
      if (const auto it = _mortar_interface_coupling.find(secondary_elem_ids[i]);
          it != _mortar_interface_coupling.end())
        coupled_elem_ids[i].assign(it->second.begin(), it->second.end());
  };
  auto act_on_data_functor = [this](processor_id_type,
                                    const std::vector<dof_id_type> & secondary_elem_ids,
                                    const std::vector<std::vector<dof_id_type>> & coupled_elem_ids)
  {
    for (const auto i : index_range(secondary_elem_ids))
      _mortar_interface_coupling[secondary_elem_ids[i]].insert(coupled_elem_ids[i].begin(),
                                                               coupled_elem_ids[i].end());
  };
  std::vector<dof_id_type> * example = nullptr;
  TIMPI::pull_parallel_vector_data(
      _mesh.comm(), secondary_elems_to_query, gather_functor, act_on_data_functor, example);
}

Referenced by buildMortarSegmentMesh(), and buildMortarSegmentMesh3d().

buildMortarSegmentMesh()

void AutomaticMortarGeneration::buildMortarSegmentMesh

(

)

Builds the mortar segment mesh once the secondary and primary node projections have been completed.

Inputs:

• mesh
• primary_node_and_elem_to_xi1_secondary_elem
• secondary_node_and_elem_to_xi2_primary_elem
• nodes_to_primary_elem_map

Outputs:

• mortar_segment_mesh
• msm_elem_to_info

Defined in the file build_mortar_segment_mesh.C.

This was a change to how inactive LM DoFs are handled. Now mortar segment elements are not used in assembly if there is no corresponding primary element and inactive LM DoFs (those with no contribution to an active primary element) are zeroed.

Definition at line 486 of file AutomaticMortarGeneration.C.

{
  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)
    {
      static const auto invalid_orientation_id =
          moose::internal::getSolutionInvalidityRegistry().registerInvalidity(
              "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.",
              true);
      _app.solutionInvalidity().flagInvalidSolutionInternal(invalid_orientation_id);
      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();
}

Referenced by MortarInterfaceWarehouse::update().

buildMortarSegmentMesh3d()

void AutomaticMortarGeneration::buildMortarSegmentMesh3d

(

)

Builds the mortar segment mesh once the secondary and primary node projections have been completed.

Inputs:

• mesh

Outputs:

• mortar_segment_mesh
• msm_elem_to_info

Step 1: Build mortar segments for all secondary elements

Step 1.1: Linearize secondary face elements

For first order face elements (Tri3 and Quad4) elements are simply linearized around center For second order (Tri6 and Quad9) and third order (Tri7) face elements, elements are sub-divided into four first order elements then each of the sub-elements is linearized around their respective centers For Quad8 elements, they are sub-divided into one quad and four triangle elements and each sub-element is linearized around their respective centers

Step 1.2: Coarse screening using a k-d tree to find nodes on the primary interface that are 'close to' a center point of the secondary element.

Step 1.3: Loop through primary candidate nodes, create mortar segments

Once an element with non-trivial projection onto secondary element identified, switch to breadth-first search (drop all current candidates and add only neighbors of elements with non-trivial overlap)

Step 1.3.2: Sub-divide primary element candidate, then project onto secondary sub-elements, perform polygon clipping, and triangulate to form mortar segments

Step 1.3.3: Create mortar segments and add to mortar segment mesh

Definition at line 996 of file AutomaticMortarGeneration.C.

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

Referenced by MortarInterfaceWarehouse::update().

buildNodeToElemMaps()

void AutomaticMortarGeneration::buildNodeToElemMaps

(

)

Once the secondary_requested_boundary_ids and primary_requested_boundary_ids containers have been filled in, call this function to build node-to-Elem maps for the lower-dimensional elements.

Definition at line 337 of file AutomaticMortarGeneration.C.

{
  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);
    }
  }
}

Referenced by MortarInterfaceWarehouse::update().

clear()

void AutomaticMortarGeneration::clear

(

)

Clears the mortar segment mesh and accompanying data structures.

Definition at line 298 of file AutomaticMortarGeneration.C.

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

Referenced by MortarInterfaceWarehouse::update().

computeADNodalNormals()

void AutomaticMortarGeneration::computeADNodalNormals	(	const std::function< ADPoint(const Node &)> &	coordinate,
		std::unordered_map< const Node *, ADRealVectorValue > &	nodal_normals
	)		const

Build the normalized JxW-weighted secondary nodal normals from AD nodal coordinates.

The coordinate functor adds displacement derivatives to each node's current coordinates, so the stored normal values and their derivatives use the same geometry state.

Definition at line 2089 of file AutomaticMortarGeneration.C.

{
  const auto dim = _mesh.mesh_dimension();
  mooseAssert(dim == 2 || dim == 3, "AD nodal normals are only valid for 2D or 3D mortar.");
 
  libMesh::QNodal qface(dim - 1);
  const auto & qweights = qface.get_weights();
  const Real sign = _periodic ? -1 : 1;
 
  std::unordered_map<const Node *, ADRealVectorValue> weighted_area_vectors;
  std::unordered_map<const Node *, Real> weighted_area_magnitudes;
 
  // One face FE object per secondary element order, reused across elements of that order
  std::map<Order, std::unique_ptr<FEBase>> fe_faces;
 
  // AD propagates coordinate derivatives through the parent-side area vectors and normalization.
  for (const auto secondary_subdomain_id : _secondary_boundary_subdomain_ids)
    for (const Elem * const secondary_elem :
         as_range(_mesh.active_subdomain_elements_begin(secondary_subdomain_id),
                  _mesh.active_subdomain_elements_end(secondary_subdomain_id)))
    {
      const Elem * const interior_parent = secondary_elem->interior_parent();
      mooseAssert(interior_parent,
                  "No interior parent exists for element "
                      << secondary_elem->id()
                      << ". There may be a problem with your sideset set-up.");
 
      const Real element_size = secondary_elem->hmax();
      if (element_size <= 0)
        mooseError("Cannot compute an AD nodal normal for secondary mortar element ",
                   secondary_elem->id(),
                   " because it has a non-positive characteristic length.");
 
      const auto side = interior_parent->which_side_am_i(secondary_elem);
      auto & fe_face = fe_faces[secondary_elem->default_order()];
      if (!fe_face)
      {
        fe_face = FEBase::build(dim, FEType(secondary_elem->default_order(), LAGRANGE));
        fe_face->attach_quadrature_rule(&qface);
      }
 
      const auto & face_normals = fe_face->get_normals();
      const auto & face_points = fe_face->get_xyz();
      const auto & map_dpsidxi = fe_face->get_fe_map().get_dpsidxi();
      const auto * const map_dpsideta = dim == 3 ? &fe_face->get_fe_map().get_dpsideta() : nullptr;
 
      fe_face->reinit(interior_parent, side);
 
      const auto qpoint_to_secondary_node =
          MooseMeshUtils::nodalQuadraturePointToNodeMap(*secondary_elem, face_points);
 
      mooseAssert(face_normals.size() == face_points.size() &&
                      qweights.size() == face_points.size(),
                  "Face nodal geometry vectors and quadrature weights must have the same size.");
 
      // The parent-side map includes enriched face mappings such as TET14-to-TRI7.
      auto parent_side_elem = interior_parent->build_side_ptr(side);
      const unsigned int n_mapping_shape_functions =
          dim == 2 ? libMesh::FE<2, LAGRANGE>::n_dofs(parent_side_elem.get(),
                                                      parent_side_elem->default_order())
                   : libMesh::FE<3, LAGRANGE>::n_dofs(parent_side_elem.get(),
                                                      parent_side_elem->default_order());
      if (map_dpsidxi.size() != n_mapping_shape_functions ||
          (map_dpsideta && map_dpsideta->size() != n_mapping_shape_functions))
        mooseError("The parent-side geometric map for secondary mortar element ",
                   secondary_elem->id(),
                   " of type ",
                   libMesh::Utility::enum_to_string<ElemType>(secondary_elem->type()),
                   " does not have one mapping shape function per element node.");
      mooseAssert(n_mapping_shape_functions <= parent_side_elem->n_nodes(),
                  "The mapping basis must not index past the side element's nodes.");
 
      std::vector<ADPoint> side_points;
      side_points.reserve(n_mapping_shape_functions);
      for (const auto side_node_index : make_range(n_mapping_shape_functions))
        side_points.push_back(coordinate(parent_side_elem->node_ref(side_node_index)));
 
      for (const auto qp : index_range(face_points))
      {
        const auto secondary_node_index = qpoint_to_secondary_node[qp];
        const Node * const secondary_node = secondary_elem->node_ptr(secondary_node_index);
 
        ADRealVectorValue tangent_xi;
        ADRealVectorValue tangent_eta;
        for (const auto side_node_index : make_range(n_mapping_shape_functions))
        {
          tangent_xi.add_scaled(side_points[side_node_index], map_dpsidxi[side_node_index][qp]);
          if (dim == 3)
            tangent_eta.add_scaled(side_points[side_node_index],
                                   (*map_dpsideta)[side_node_index][qp]);
        }
 
        auto area_vector = dim == 2 ? ADRealVectorValue(tangent_xi(1), -tangent_xi(0), 0)
                                    : tangent_xi.cross(tangent_eta);
        const Real area = MetaPhysicL::raw_value(area_vector.norm());
        const Real degeneracy_tolerance =
            dim == 2 ? TOLERANCE * element_size : TOLERANCE * element_size * element_size;
        if (area <= degeneracy_tolerance)
          mooseError("Cannot compute an AD nodal normal on degenerate secondary mortar element ",
                     secondary_elem->id(),
                     ".");
 
        // qweight * area_vector is equivalent to JxW * face_normal. Fix its orientation to the
        // Real-valued face normal because the incident-face set and orientation are not
        // differentiated.
        const Real orientation =
            MetaPhysicL::raw_value(area_vector) * face_normals[qp] < 0 ? -1 : 1;
 
        weighted_area_vectors[secondary_node] += sign * orientation * qweights[qp] * area_vector;
        weighted_area_magnitudes[secondary_node] += std::abs(qweights[qp]) * area;
      }
    }
 
  nodal_normals.clear();
  nodal_normals.reserve(weighted_area_vectors.size());
  for (auto & [secondary_node, weighted_area_vector] : weighted_area_vectors)
  {
    const Real weighted_area_vector_norm = MetaPhysicL::raw_value(weighted_area_vector.norm());
    const Real weighted_area_magnitude = libmesh_map_find(weighted_area_magnitudes, secondary_node);
    if (weighted_area_vector_norm <= TOLERANCE * weighted_area_magnitude)
      mooseError("Cannot compute a normalized AD nodal normal for secondary node ",
                 secondary_node->id(),
                 " because its incident face contributions cancel to a near-zero weighted normal.");
 
    auto nodal_normal = weighted_area_vector.unit();
    const auto & stored_normal = libmesh_map_find(_secondary_node_to_nodal_normal, secondary_node);
    // Equivalent FE and AD evaluations of the same unit normal may differ through operation
    // ordering. The tolerance admits roundoff while detecting inconsistent geometry states.
    mooseAssert((MetaPhysicL::raw_value(nodal_normal) - stored_normal).norm() < 1e-10,
                "The stored and AD secondary nodal normals must use the same geometry state.");
    // Set the raw value from the stored contact direction while retaining the derivatives produced
    // by the equivalent AD area-vector calculation.
    for (const auto component : make_range(Moose::dim))
      nodal_normal(component).value() = stored_normal(component);
    nodal_normals.emplace(secondary_node, std::move(nodal_normal));
  }
}

computeInactiveLMElems()

void AutomaticMortarGeneration::computeInactiveLMElems

(

)

Get list of secondary elems without any corresponding primary elements.

Used to enforce zero values on inactive DoFs of elemental variables.

Definition at line 1941 of file AutomaticMortarGeneration.C.

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

Referenced by MortarInterfaceWarehouse::update().

computeInactiveLMNodes()

void AutomaticMortarGeneration::computeInactiveLMNodes

(

)

Get list of secondary nodes that don't contribute to interaction with any primary element.

Used to enforce zero values on inactive DoFs of nodal variables.

Definition at line 1877 of file AutomaticMortarGeneration.C.

{
  if (!_correct_edge_dropping)
  {
    computeIncorrectEdgeDroppingInactiveLMNodes();
    return;
  }
 
  // Only the owner of a node zeroes its LM DoFs, but the owner need not own (or even ghost) any
  // secondary lower-dimensional element touching it. Each processor therefore sends every node of
  // its local secondary elements to the node's owner, flagged as active if a local mortar segment
  // touches it. A node is inactive if no processor flags it as active.
 
  // 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));
  }
 
  // List of all nodes on local secondary elements
  std::unordered_set<dof_id_type> local_nodes;
  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()))
        local_nodes.insert(el->node_id(n));
 
  // Pairs of (node id, 1 if active on the sending processor and 0 otherwise)
  std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned char>>>
      proc_to_nodes;
  for (const auto node_id : local_nodes)
    proc_to_nodes[_mesh.node_ptr(node_id)->processor_id()].emplace_back(
        node_id, active_local_nodes.count(node_id));
 
  std::unordered_set<dof_id_type> candidate_nodes, active_nodes;
  auto action_functor = [&candidate_nodes, &active_nodes](
                            const processor_id_type,
                            const std::vector<std::pair<dof_id_type, unsigned char>> & sent_data)
  {
    for (const auto & [node_id, active] : sent_data)
    {
      candidate_nodes.insert(node_id);
      if (active)
        active_nodes.insert(node_id);
    }
  };
  TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_nodes, action_functor);
 
  // Store the owned inactive nodes to use later to zero LM DoFs on inactive nodes
  _inactive_local_lm_nodes.clear();
  for (const auto node_id : candidate_nodes)
    if (!active_nodes.count(node_id))
      _inactive_local_lm_nodes.insert(_mesh.node_ptr(node_id));
}

Referenced by MortarInterfaceWarehouse::update().

computeIncorrectEdgeDroppingInactiveLMNodes()

void AutomaticMortarGeneration::computeIncorrectEdgeDroppingInactiveLMNodes

(

)

Computes inactive secondary nodes when incorrect edge dropping behavior is enabled (any node touching a partially or fully dropped element is dropped)

Definition at line 1821 of file AutomaticMortarGeneration.C.

{
  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;
 
  // 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));
      }
 
  // Only the owner of a node zeroes its LM DoFs, but the owner need not own (or even ghost) any
  // secondary lower-dimensional element touching it, so send each inactive node to its owner,
  // which takes the union
  std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_inactive_nodes;
  for (const auto node_id : inactive_node_ids)
    proc_to_inactive_nodes[_mesh.node_ptr(node_id)->processor_id()].push_back(node_id);
 
  _inactive_local_lm_nodes.clear();
  auto action_functor = [this](const processor_id_type, const std::vector<dof_id_type> & sent_data)
  {
    for (const auto node_id : sent_data)
      _inactive_local_lm_nodes.insert(_mesh.node_ptr(node_id));
  };
  TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_inactive_nodes, action_functor);
}

Referenced by computeInactiveLMNodes().

computeMsmStatistics()

std::vector< AutomaticMortarGeneration::MsmSubdomainStats > AutomaticMortarGeneration::computeMsmStatistics

(

)

Computes mortar segment mesh statistics and returns one entry per subdomain pair.

Must be called collectively on all ranks.

Definition at line 1692 of file AutomaticMortarGeneration.C.

{
  std::vector<MsmSubdomainStats> result;
  libMesh::StatisticsVector<Real> primary;
  libMesh::StatisticsVector<Real> secondary;
  libMesh::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(cast_ref<std::vector<Real> &>(secondary));
    _mesh.comm().allgather(cast_ref<std::vector<Real> &>(msm));
    primary.reserve(primary_elems_to_volume.size());
    for (const auto [_, volume] : primary_elems_to_volume)
      primary.push_back(volume);
 
    // An interface on which no mortar segments formed has empty primary and segment lists, for
    // which libMesh's maximum() and minimum() are undefined. Report 0, as median() does.
    const auto max_volume = [](const auto & volumes)
    { return volumes.empty() ? Real(0) : volumes.maximum(); };
    const auto min_volume = [](const auto & volumes)
    { return volumes.empty() ? Real(0) : volumes.minimum(); };
 
    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 = max_volume(secondary);
    stats.secondary_lower_min_volume = min_volume(secondary);
    stats.secondary_lower_median_volume = secondary.median();
    stats.primary_lower_n_elems = primary.size();
    stats.primary_lower_max_volume = max_volume(primary);
    stats.primary_lower_min_volume = min_volume(primary);
    stats.primary_lower_median_volume = primary.median();
    stats.msm_n_elems = msm.size();
    stats.msm_max_volume = max_volume(msm);
    stats.msm_min_volume = min_volume(msm);
    stats.msm_median_volume = msm.median();
    result.push_back(stats);
 
    primary.clear();
    secondary.clear();
    msm.clear();
    primary_elems_to_volume.clear();
  }
 
  return result;
}

Referenced by msmStatistics().

computeNodalGeometry()

void AutomaticMortarGeneration::computeNodalGeometry

(

)

Computes and stores the nodal normal/tangent vectors in a local data structure instead of using the ExplicitSystem/NumericVector approach.

This design was triggered by the way that the GhostingFunctor operates, but I think it is a better/more efficient way to do it anyway.

The _periodic flag tells us whether we want to inward vs outward facing normals

Definition at line 1989 of file AutomaticMortarGeneration.C.

{
  // 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.
  libMesh::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 =
        MooseMeshUtils::nodalQuadraturePointToNodeMap(*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]));
    }
  }
 
  // Normalize the JxW-weighted sum of incident face normals at each secondary node.
  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;
    _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)] =
        Moose::Mortar::nodalTangents(nodal_normal, dim);
  }
}

Referenced by MortarInterfaceWarehouse::update().

dim()

int AutomaticMortarGeneration::dim ( ) const

inline

Definition at line 415 of file AutomaticMortarGeneration.h.

{ return _mesh.mesh_dimension(); }

Referenced by computeADNodalNormals(), computeInactiveLMElems(), computeIncorrectEdgeDroppingInactiveLMNodes(), computeNodalGeometry(), Moose::Mortar::loopOverMortarSegments(), and MortarInterfaceWarehouse::update().

getInactiveLMElems()

const std::unordered_set< const Elem * > & AutomaticMortarGeneration::getInactiveLMElems ( ) const

inline

Returns

The list of secondary elems on which mortar constraint is not active

Definition at line 435 of file AutomaticMortarGeneration.h.

  {
    return _inactive_local_lm_elems;
  }

Referenced by ComputeMortarFunctor::operator()().

getInactiveLMNodes()

const std::unordered_set< const Node * > & AutomaticMortarGeneration::getInactiveLMNodes ( ) const

inline

Returns

The set of nodes on which mortar constraints are not active

Definition at line 427 of file AutomaticMortarGeneration.h.

  {
    return _inactive_local_lm_nodes;
  }

Referenced by ComputeMortarFunctor::operator()().

getNodalNormals()

std::vector< Point > AutomaticMortarGeneration::getNodalNormals

(

const Elem &

secondary_elem

)

const

Returns

The nodal normals associated with the provided secondary_elem

Definition at line 375 of file AutomaticMortarGeneration.C.

{
  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;
}

Referenced by buildMortarSegmentMesh3d(), getNormals(), and MortarConsumerInterface::setNormals().

getNodalTangents()

std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 > AutomaticMortarGeneration::getNodalTangents

(

const Elem &

secondary_elem

)

const

Compute the two nodal tangents, which are built on-the-fly.

Returns

The nodal tangents associated with the provided secondary_elem

Definition at line 428 of file AutomaticMortarGeneration.C.

{
  // 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}};
}

getNormals() [1/2]

std::vector< Point > AutomaticMortarGeneration::getNormals	(	const Elem &	secondary_elem,
		const std::vector< Point > &	xi1_pts
	)		const

Compute the normals at given reference points on a secondary element.

Parameters

secondary_elem	The secondary element used to query for associated nodal normals
xi1_pts	The reference points on the secondary element to evaluate the normals at. The points should only be non-zero in the zeroth entry because right now our mortar mesh elements are always 1D

Returns

The normals

Definition at line 457 of file AutomaticMortarGeneration.C.

{
  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,
                                            cast_ref<const TypeVector<Real> &>(xi1_pts[qp]));
      normals[qp] += phi * nodal_normals[n];
    }
 
  if (_periodic)
    for (auto & normal : normals)
      normal *= -1;
 
  return normals;
}

Referenced by buildMortarSegmentMesh(), getNormals(), and MortarConsumerInterface::setNormals().

getNormals() [2/2]

std::vector< Point > AutomaticMortarGeneration::getNormals	(	const Elem &	secondary_elem,
		const std::vector< Real > &	oned_xi1_pts
	)		const

Compute the normals at given reference points on a secondary element.

Parameters

secondary_elem	The secondary element used to query for associated nodal normals
1d_xi1_pts	The reference points on the secondary element to evaluate the normals at. The "points" are single reals corresponding to xi because right now our mortar mesh elements are always 1D

Returns

The normals

Definition at line 446 of file AutomaticMortarGeneration.C.

{
  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);
}

getPrimaryIpToLowerElementMap()

std::map< unsigned int, unsigned int > AutomaticMortarGeneration::getPrimaryIpToLowerElementMap	(	const Elem &	primary_elem,
		const Elem &	primary_elem_ip,
		const Elem &	lower_secondary_elem
	)		const

Compute on-the-fly mapping from primary interior parent nodes to its corresponding lower dimensional nodes.

Returns

The map from primary interior parent nodes to its corresponding lower dimensional nodes

Definition at line 411 of file AutomaticMortarGeneration.C.

{
  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;
}

getSecondaryIpToLowerElementMap()

std::map< unsigned int, unsigned int > AutomaticMortarGeneration::getSecondaryIpToLowerElementMap

(

const Elem &

lower_secondary_elem

)

const

Compute on-the-fly mapping from secondary interior parent nodes to lower dimensional nodes.

Returns

The map from secondary interior parent nodes to lower dimensional nodes

Definition at line 395 of file AutomaticMortarGeneration.C.

{
  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;
}

getSecondaryLowerdElemFromSecondaryElem()

const Elem * AutomaticMortarGeneration::getSecondaryLowerdElemFromSecondaryElem

(

dof_id_type

secondary_elem_id

)

const

Return lower dimensional secondary element given its interior parent.

Helpful outside the mortar generation to locate mortar-related quantities.

Parameters

secondary_elem_id

The secondary interior parent element id used to query for associated lower dimensional element

Returns

The corresponding lower dimensional secondary element

Definition at line 385 of file AutomaticMortarGeneration.C.

{
  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);
}

incorrectEdgeDropping()

bool AutomaticMortarGeneration::incorrectEdgeDropping ( ) const

inline

Definition at line 440 of file AutomaticMortarGeneration.h.

{ return !_correct_edge_dropping; }

Referenced by ComputeMortarFunctor::operator()().

initOutput()

void AutomaticMortarGeneration::initOutput

(

)

initialize mortar-mesh based output

Definition at line 282 of file AutomaticMortarGeneration.C.

{
  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));
}

meshChanged()

void AutomaticMortarGeneration::meshChanged ( )

inline

Invalidates the cached MSM node/element ID starting offset so that the next call to buildMortarSegmentMesh3d() recomputes it via allgather.

Call this when mesh topology changes.

Definition at line 282 of file AutomaticMortarGeneration.h.

{ _msm_node_id_start = std::nullopt; }

Referenced by MortarInterfaceWarehouse::meshChanged().

mortar3DQpMapping()

Mortar3DQuadraturePointMapping AutomaticMortarGeneration::mortar3DQpMapping ( ) const

inline

Return the 3D mortar quadrature-point mapping method.

Definition at line 418 of file AutomaticMortarGeneration.h.

{ return _mortar_3d_qp_mapping; }

Referenced by Moose::Mortar::loopOverMortarSegments().

mortarInterfaceCoupling()

const std::unordered_map< dof_id_type, std::unordered_set< dof_id_type > > & AutomaticMortarGeneration::mortarInterfaceCoupling ( ) const

inline

Returns

The mortar interface coupling

Definition at line 382 of file AutomaticMortarGeneration.h.

  {
    return _mortar_interface_coupling;
  }

mortarInterfaceName()

std::string AutomaticMortarGeneration::mortarInterfaceName ( ) const

private

Returns

A string uniquely identifying this mortar interface

Definition at line 268 of file AutomaticMortarGeneration.C.

{
  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, "_");
}

Referenced by initOutput(), and outputMortarMesh().

mortarSegmentMesh()

const MeshBase & AutomaticMortarGeneration::mortarSegmentMesh ( ) const

inline

Returns

The mortar segment mesh

Definition at line 405 of file AutomaticMortarGeneration.h.

{ return *_mortar_segment_mesh; }

Referenced by Moose::Mortar::loopOverMortarSegments().

mortarSegmentMeshElemToInfo()

const std::unordered_map< const Elem *, MortarSegmentInfo > & AutomaticMortarGeneration::mortarSegmentMeshElemToInfo ( ) const

inline

Returns

The mortar segment element to corresponding information

Definition at line 410 of file AutomaticMortarGeneration.h.

  {
    return _msm_elem_to_info;
  }

Referenced by Moose::Mortar::loopOverMortarSegments().

mortarSegmentReferencePoints()

const MortarSegmentReferencePoints & AutomaticMortarGeneration::mortarSegmentReferencePoints

(

const Elem &

mortar_segment_elem

)

const

Return the parent-face reference coordinates for a mortar segment.

Definition at line 320 of file AutomaticMortarGeneration.C.

{
  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;
}

Referenced by Moose::Mortar::loopOverMortarSegments().

msmStatistics()

void AutomaticMortarGeneration::msmStatistics

(

)

Prints mortar segment mesh statistics to console (calls computeMsmStatistics internally)

Definition at line 1776 of file AutomaticMortarGeneration.C.

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

Referenced by buildMortarSegmentMesh3d().

nodesToSecondaryElem()

const std::unordered_map< dof_id_type, std::vector< const Elem * > > & AutomaticMortarGeneration::nodesToSecondaryElem ( ) const

inline

Returns

Map from node id to secondary lower-d element pointer

Definition at line 476 of file AutomaticMortarGeneration.h.

  {
    return _nodes_to_secondary_elem_map;
  }

onDisplaced()

bool AutomaticMortarGeneration::onDisplaced ( ) const

inline

returns whether this object is on the displaced mesh

Definition at line 287 of file AutomaticMortarGeneration.h.

{ return _on_displaced; }

outputMortarMesh()

void AutomaticMortarGeneration::outputMortarMesh ( )

private

Write the mortar segment mesh to exodus.

Definition at line 981 of file AutomaticMortarGeneration.C.

{
  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, "_"));
}

Referenced by buildMortarSegmentMesh(), and buildMortarSegmentMesh3d().

primaryIPSubIDs()

const std::set< SubdomainID > & AutomaticMortarGeneration::primaryIPSubIDs ( ) const

inline

Returns

All the primary interior parent subdomain IDs associated with the mortar mesh

Definition at line 471 of file AutomaticMortarGeneration.h.

{ return _primary_ip_sub_ids; }

Referenced by Moose::Mortar::setupMortarMaterials(), and MortarInterfaceWarehouse::update().

primarySecondaryBoundaryIDPair()

const std::pair< BoundaryID, BoundaryID > & AutomaticMortarGeneration::primarySecondaryBoundaryIDPair ( ) const

inline

Returns

The primary-secondary boundary ID pair

Definition at line 705 of file AutomaticMortarGeneration.h.

{
  mooseAssert(_primary_secondary_boundary_id_pairs.size() == 1,
              "We currently only support a single boundary pair per mortar generation object");
 
  return _primary_secondary_boundary_id_pairs.front();
}

Referenced by Moose::Mortar::loopOverMortarSegments(), and Moose::Mortar::setupMortarMaterials().

processAlignedNodes()

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_element_candidates  )

private

Process aligned nodes.

Returns

whether mortar segment(s) were created

Definition at line 2241 of file AutomaticMortarGeneration.C.

{
  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;
}

Referenced by projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

projectPrimaryNodes()

void AutomaticMortarGeneration::projectPrimaryNodes

(

)

(Inverse) project primary nodes to the points on the secondary surface where they would have come from (find (xi^(1) values)).

Inputs:

• The nodal normals values
• mesh
• nodes_to_secondary_elem_map

Outputs:

• primary_node_and_elem_to_xi1_secondary_elem

Defined in the file project_primary_nodes.C.

Definition at line 2602 of file AutomaticMortarGeneration.C.

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

Referenced by MortarInterfaceWarehouse::update().

projectPrimaryNodesSinglePair()

void AutomaticMortarGeneration::projectPrimaryNodesSinglePair ( SubdomainID  lower_dimensional_primary_subdomain_id, SubdomainID  lower_dimensional_secondary_subdomain_id  )

private

Helper function used internally by AutomaticMortarGeneration::project_primary_nodes().

Definition at line 2611 of file AutomaticMortarGeneration.C.

{
  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 "
                 << cast_ref<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.
}

Referenced by projectPrimaryNodes().

projectSecondaryNodes()

void AutomaticMortarGeneration::projectSecondaryNodes

(

)

Project secondary nodes (find xi^(2) values) to the closest points on the primary surface.

Inputs:

• The nodal normals values
• mesh
• nodes_to_primary_elem_map

Outputs:

• secondary_node_and_elem_to_xi2_primary_elem

Defined in the file project_secondary_nodes.C.

Definition at line 2232 of file AutomaticMortarGeneration.C.

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

Referenced by MortarInterfaceWarehouse::update().

projectSecondaryNodesSinglePair()

void AutomaticMortarGeneration::projectSecondaryNodesSinglePair ( SubdomainID  lower_dimensional_primary_subdomain_id, SubdomainID  lower_dimensional_secondary_subdomain_id  )

private

Helper function responsible for projecting secondary nodes onto primary elements for a single primary/secondary pair.

Called by the class member AutomaticMortarGeneration::project_secondary_nodes().

Definition at line 2357 of file AutomaticMortarGeneration.C.

{
  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 "
                   << cast_ref<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;
}

Referenced by projectSecondaryNodes().

requestNonlocalSecondaryCouplings()

void AutomaticMortarGeneration::requestNonlocalSecondaryCouplings ( )

inline

Request that every process receive the mortar interface couplings of the nonlocal secondary lower-dimensional elements that point-neighbor ghosting reads for its local elements: the point neighbors of its local secondary elements and of the secondary elements coupled to its local elements.

Consumers that apply a node's constraint on every process contributing to it, rather than only on the node's owner, need these couplings to preallocate the rows of their local elements.

Definition at line 395 of file AutomaticMortarGeneration.h.

{ _nonlocal_secondary_couplings = true; }

secondariesToMortarSegments() [1/2]

const std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > > & AutomaticMortarGeneration::secondariesToMortarSegments ( ) const

inline

Returns

the lower dimensional secondary element ids and their associated mortar segment elements

Definition at line 458 of file AutomaticMortarGeneration.h.

  {
    return _secondary_elems_to_mortar_segments;
  }

secondariesToMortarSegments() [2/2]

std::vector< AutomaticMortarGeneration::MortarFilterIter > AutomaticMortarGeneration::secondariesToMortarSegments

(

const Node &

node

)

const

Returns

A vector of iterators that point to the lower dimensional secondary elements and their associated mortar segment elements that would have nonzero values for a Lagrange shape function associated with the provided node. This method may return an empty container if the node is away from the mortar mesh

Definition at line 2838 of file AutomaticMortarGeneration.C.

{
  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;
}

Referenced by MortarUserObjectThread::operator()(), and ComputeMortarFunctor::operator()().

secondaryIPSubIDs()

const std::set< SubdomainID > & AutomaticMortarGeneration::secondaryIPSubIDs ( ) const

inline

Returns

All the secondary interior parent subdomain IDs associated with the mortar mesh

Definition at line 466 of file AutomaticMortarGeneration.h.

{ return _secondary_ip_sub_ids; }

Referenced by Moose::Mortar::setupMortarMaterials(), and MortarInterfaceWarehouse::update().

Friends And Related Symbol Documentation

AugmentSparsityOnInterface

friend class AugmentSparsityOnInterface

friend

Definition at line 701 of file AutomaticMortarGeneration.h.

MortarNodalGeometryOutput

friend class MortarNodalGeometryOutput

friend

Definition at line 700 of file AutomaticMortarGeneration.h.

Member Data Documentation

_app

MooseApp& AutomaticMortarGeneration::_app

private

The Moose app.

Definition at line 503 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), initOutput(), and outputMortarMesh().

_console

const ConsoleStream ConsoleStreamInterface::_console

inherited

An instance of helper class to write streams to the Console objects.

Definition at line 31 of file ConsoleStreamInterface.h.

Referenced by IterationAdaptiveDT::acceptStep(), MaterialOutputAction::act(), MeshOnlyAction::act(), SetupDebugAction::act(), FEProblemBase::adaptMesh(), Adaptivity::adaptMesh(), PerfGraph::addToExecutionList(), SimplePredictor::apply(), SystemBase::applyScalingFactors(), MultiApp::backup(), FEProblemBase::backupMultiApps(), CoarsenedPiecewiseLinear::buildCoarsenedGrid(), DefaultSteadyStateConvergence::checkConvergence(), MeshDiagnosticsGenerator::checkElementOverlap(), MeshDiagnosticsGenerator::checkElementTypes(), MeshDiagnosticsGenerator::checkElementVolumes(), FEProblemBase::checkExceptionAndStopSolve(), SolverSystem::checkInvalidSolution(), MeshDiagnosticsGenerator::checkLocalJacobians(), MeshDiagnosticsGenerator::checkNonConformalMesh(), MeshDiagnosticsGenerator::checkNonConformalMeshFromAdaptivity(), MeshDiagnosticsGenerator::checkNonConformingFaces(), MeshDiagnosticsGenerator::checkNonMatchingEdges(), MeshDiagnosticsGenerator::checkNonPlanarSides(), MeshDiagnosticsGenerator::checkPolygons(), FEProblemBase::checkProblemIntegrity(), ReferenceResidualConvergence::checkResidualConvergence(), MeshDiagnosticsGenerator::checkSidesetsOrientation(), MeshDiagnosticsGenerator::checkWatertightNodesets(), MeshDiagnosticsGenerator::checkWatertightSidesets(), CoarsenSurfaceMeshAlongSidesetGenerator::coarsenAlongSidesets(), IterationAdaptiveDT::computeAdaptiveDT(), TransientBase::computeConstrainedDT(), DefaultMultiAppFixedPointConvergence::computeCustomConvergencePostprocessor(), NonlinearSystemBase::computeDamping(), FixedPointIterationAdaptiveDT::computeDT(), IterationAdaptiveDT::computeDT(), IterationAdaptiveDT::computeFailedDT(), IterationAdaptiveDT::computeInitialDT(), IterationAdaptiveDT::computeInterpolationDT(), FEProblemBase::computeLinearSystemTags(), LinearSystem::computeLinearSystemTags(), NonlinearSystemBase::computeScaling(), Problem::console(), TimeStepper::constrainStep(), IterationAdaptiveDT::constrainStep(), MultiApp::createApp(), FEProblemBase::execMultiApps(), FEProblemBase::execMultiAppTransfers(), Eigenvalue::execute(), SteadyBase::execute(), MFEMSteady::execute(), MessageFromInput::execute(), ActionWarehouse::executeActionsWithAction(), ActionWarehouse::executeAllActions(), MeshGeneratorSystem::executeMeshGenerators(), SidesetAroundSubdomainUpdater::finalize(), ElementQualityChecker::finalize(), FEProblemBase::finishMultiAppStep(), MeshRepairGenerator::fixOverlappingNodes(), SurfaceSubdomainsDelaunayRemesher::General2DDelaunay(), CoarsenBlockGenerator::generate(), CoarsenSurfaceMeshAlongSidesetGenerator::generate(), OrientSurfaceMeshGenerator::generate(), PolyLineMeshFollowingNodeSetGenerator::generate(), StitchMeshGenerator::generate(), MeshGenerator::generateInternal(), VariableCondensationPreconditioner::getDofToCondense(), InversePowerMethod::init(), NonlinearEigen::init(), FEProblemBase::initialAdaptMesh(), SubProblem::initialSetup(), EigenExecutionerBase::inversePowerIteration(), FEProblemBase::joinAndFinalize(), TransientBase::keepGoing(), IterationAdaptiveDT::limitDTByFunction(), IterationAdaptiveDT::limitDTToPostprocessorValue(), FEProblemBase::logAdd(), EigenExecutionerBase::makeBXConsistent(), Console::meshChanged(), SurfaceDelaunayGeneratorBase::meshNormalDeviation2D(), MooseBase::mooseDeprecated(), MooseBase::mooseDeprecatedNoTrace(), MooseBase::mooseInfo(), MooseBase::mooseWarning(), MooseBase::mooseWarningNonPrefixed(), NEML2Action::NEML2Action(), ReferenceResidualConvergence::nonlinearConvergenceSetup(), Console::output(), DOFMapOutput::output(), MaterialPropertyDebugOutput::output(), PerfGraphOutput::output(), ReporterDebugOutput::output(), SolutionInvalidityOutput::output(), VariableResidualNormsDebugOutput::output(), ControlOutput::outputActiveObjects(), ControlOutput::outputChangedControls(), ControlOutput::outputControls(), Console::outputInput(), WebServerControl::outputMessage(), Console::outputPostprocessors(), PseudoTimestep::outputPseudoTimestep(), Console::outputReporters(), DefaultMultiAppFixedPointConvergence::outputResidualNorm(), Console::outputScalarVariables(), Console::outputSystemInformation(), FEProblemBase::possiblyRebuildGeomSearchPatches(), EigenExecutionerBase::postExecute(), AB2PredictorCorrector::postSolve(), DefaultMultiAppFixedPointConvergence::preLoop(), ActionWarehouse::printActionDependencySets(), BlockRestrictionDebugOutput::printBlockRestrictionGroups(), BlockRestrictionDebugOutput::printBlockRestrictionMap(), BlockRestrictionDebugOutput::printBoundaryRestrictionGroups(), SolutionInvalidity::printDebug(), EigenExecutionerBase::printEigenvalue(), PicardSolve::printFixedPointConvergenceHistory(), SecantSolve::printFixedPointConvergenceHistory(), SteffensenSolve::printFixedPointConvergenceHistory(), FixedPointSolve::printFixedPointConvergenceReason(), PerfGraphLivePrint::printLiveMessage(), MaterialPropertyDebugOutput::printMaterialMap(), PerfGraphLivePrint::printStats(), NEML2Action::printSummary(), projectPrimaryNodesSinglePair(), projectSecondaryNodesSinglePair(), CoarsenBlockGenerator::recursiveCoarsen(), SolutionTimeAdaptiveDT::rejectStep(), MultiApp::restore(), FEProblemBase::restoreMultiApps(), FEProblemBase::restoreSolutions(), NonlinearSystemBase::setInitialSolution(), MooseApp::setupOptions(), Checkpoint::shouldOutput(), SubProblem::showFunctorRequestors(), SubProblem::showFunctors(), FullSolveMultiApp::showStatusMessage(), FEProblemSolve::solve(), FixedPointSolve::solve(), LinearSystem::solve(), NonlinearSystem::solve(), AStableDirk4::solve(), ExplicitRK2::solve(), ExplicitTVDRK2::solve(), ImplicitMidpoint::solve(), LStableDirk2::solve(), LStableDirk3::solve(), LStableDirk4::solve(), EigenProblem::solve(), FixedPointSolve::solveStep(), TransientMultiApp::solveStep(), MeshRepairGenerator::splitNonConvexPolygons(), PerfGraphLivePrint::start(), WebServerControl::startServer(), AB2PredictorCorrector::step(), NonlinearEigen::takeStep(), TransientBase::takeStep(), MFEMTransient::takeStep(), TerminateChainControl::terminate(), SubProblem::timestepSetup(), FEProblemBase::updateMeshXFEM(), Convergence::verboseOutput(), Console::writeTimestepInformation(), Console::writeVariableNorms(), and FEProblemBase::~FEProblemBase().

_correct_edge_dropping

const bool AutomaticMortarGeneration::_correct_edge_dropping

private

Flag to enable regressed treatment of edge dropping where all LM DoFs on edge dropping element are strongly set to 0.

Definition at line 662 of file AutomaticMortarGeneration.h.

Referenced by computeInactiveLMElems(), computeInactiveLMNodes(), and incorrectEdgeDropping().

_debug

const bool AutomaticMortarGeneration::_debug

private

Whether to print debug output.

Definition at line 637 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), initOutput(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_distributed

const bool AutomaticMortarGeneration::_distributed

private

Whether the mortar segment mesh is distributed.

Definition at line 646 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), and projectSecondaryNodesSinglePair().

_failed_secondary_node_projections

std::unordered_set<dof_id_type> AutomaticMortarGeneration::_failed_secondary_node_projections

private

Secondary nodes that failed to project.

Definition at line 698 of file AutomaticMortarGeneration.h.

Referenced by clear(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_inactive_local_lm_elems

std::unordered_set<const Elem *> AutomaticMortarGeneration::_inactive_local_lm_elems

private

List of inactive lagrange multiplier nodes (for elemental variables)

Definition at line 596 of file AutomaticMortarGeneration.h.

Referenced by computeInactiveLMElems(), and getInactiveLMElems().

_inactive_local_lm_nodes

std::unordered_set<const Node *> AutomaticMortarGeneration::_inactive_local_lm_nodes

private

List of inactive lagrange multiplier nodes owned by this processor (for nodal variables)

Definition at line 593 of file AutomaticMortarGeneration.h.

Referenced by computeInactiveLMNodes(), computeIncorrectEdgeDroppingInactiveLMNodes(), and getInactiveLMNodes().

_lower_elem_to_side_id

std::unordered_map<const Elem *, unsigned int> AutomaticMortarGeneration::_lower_elem_to_side_id

private

Keeps track of the mapping between lower-dimensional elements and the side_id of the interior_parent which they are.

Definition at line 561 of file AutomaticMortarGeneration.h.

Referenced by clear().

_mesh

MeshBase& AutomaticMortarGeneration::_mesh

private

Reference to the mesh stored in equation_systems.

Definition at line 506 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), buildCouplingInformation(), buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), buildNodeToElemMaps(), computeADNodalNormals(), computeInactiveLMElems(), computeInactiveLMNodes(), computeIncorrectEdgeDroppingInactiveLMNodes(), computeMsmStatistics(), computeNodalGeometry(), dim(), getNormals(), msmStatistics(), MortarNodalGeometryOutput::output(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_minimum_projection_angle

const Real AutomaticMortarGeneration::_minimum_projection_angle

private

Parameter to control which angle (in degrees) is admissible for the creation of mortar segments.

If set to a value close to zero, very oblique projections are allowed, which can result in mortar segments solving physics not meaningfully and overprojection of primary nodes onto the mortar segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in 2D.

Definition at line 668 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_mortar_3d_qp_mapping

const Mortar3DQuadraturePointMapping AutomaticMortarGeneration::_mortar_3d_qp_mapping

private

Method used to map 3D mortar segment quadrature points to primary and secondary faces.

Definition at line 680 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d(), mortar3DQpMapping(), and mortarSegmentReferencePoints().

_mortar_3d_subpatch_plane

const Mortar3DSubpatchPlane AutomaticMortarGeneration::_mortar_3d_subpatch_plane

private

Method used to define the local projection planes for 3D secondary subpatches.

Definition at line 671 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d().

_mortar_interface_coupling

std::unordered_map<dof_id_type, std::unordered_set<dof_id_type> > AutomaticMortarGeneration::_mortar_interface_coupling

private

Used by the AugmentSparsityOnInterface functor to determine whether a given Elem is coupled to any others across the gap, and to explicitly set up the dependence between interior_parent() elements on the secondary side and their lower-dimensional sides which are on the interface.

This latter type of coupling must be explicitly declared when there is no primary_elem for a given mortar segment and you are using e.g. a P^1-P^0 discretization which does not induce the coupling automatically.

Definition at line 580 of file AutomaticMortarGeneration.h.

Referenced by buildCouplingInformation(), clear(), and mortarInterfaceCoupling().

_mortar_segment_mesh

std::unique_ptr<MeshBase> AutomaticMortarGeneration::_mortar_segment_mesh

private

1D Mesh of mortar segment elements which gets built by the call to build_mortar_segment_mesh().

Definition at line 552 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), clear(), computeInactiveLMElems(), computeInactiveLMNodes(), computeIncorrectEdgeDroppingInactiveLMNodes(), mortarSegmentMesh(), and outputMortarMesh().

_msm_elem_to_info

std::unordered_map<const Elem *, MortarSegmentInfo> AutomaticMortarGeneration::_msm_elem_to_info

private

Map between Elems in the mortar segment mesh and their info structs.

This gets filled in by the call to build_mortar_segment_mesh().

Definition at line 557 of file AutomaticMortarGeneration.h.

Referenced by buildCouplingInformation(), buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), clear(), computeInactiveLMElems(), computeInactiveLMNodes(), computeIncorrectEdgeDroppingInactiveLMNodes(), computeMsmStatistics(), and mortarSegmentMeshElemToInfo().

_msm_elem_to_reference_points

std::unordered_map<const Elem *, MortarSegmentReferencePoints> AutomaticMortarGeneration::_msm_elem_to_reference_points

private

Reference-coordinate data used only by the reference-interpolation mapping mode.

Definition at line 683 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d(), clear(), and mortarSegmentReferencePoints().

_msm_node_id_start

std::optional<dof_id_type> AutomaticMortarGeneration::_msm_node_id_start

private

Cached per-rank starting ID for 3D MSM nodes/elements.

nullopt forces recomputation on next buildMortarSegmentMesh3d() call. Reset by meshChanged() on topology change so the allgather is skipped for displaced-mesh residual updates that only move nodes.

Definition at line 691 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d(), and meshChanged().

_newton_tolerance

Real AutomaticMortarGeneration::_newton_tolerance = 1e-12

private

Newton solve tolerance for node projections.

Definition at line 653 of file AutomaticMortarGeneration.h.

Referenced by projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_nodes_to_primary_elem_map

std::unordered_map<dof_id_type, std::vector<const Elem *> > AutomaticMortarGeneration::_nodes_to_primary_elem_map

private

Definition at line 520 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), buildNodeToElemMaps(), clear(), processAlignedNodes(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_nodes_to_secondary_elem_map

std::unordered_map<dof_id_type, std::vector<const Elem *> > AutomaticMortarGeneration::_nodes_to_secondary_elem_map

private

Map from nodes to connected lower-dimensional elements on the secondary/primary subdomains.

Definition at line 519 of file AutomaticMortarGeneration.h.

Referenced by buildNodeToElemMaps(), clear(), nodesToSecondaryElem(), processAlignedNodes(), projectPrimaryNodesSinglePair(), projectSecondaryNodesSinglePair(), and secondariesToMortarSegments().

_nonlocal_secondary_couplings

bool AutomaticMortarGeneration::_nonlocal_secondary_couplings = false

private

Whether buildCouplingInformation() also gathers the couplings of nonlocal secondary lower-dimensional elements.

Definition at line 650 of file AutomaticMortarGeneration.h.

Referenced by buildCouplingInformation(), and requestNonlocalSecondaryCouplings().

_on_displaced

const bool AutomaticMortarGeneration::_on_displaced

private

Whether this object is on the displaced mesh.

Definition at line 640 of file AutomaticMortarGeneration.h.

Referenced by mortarInterfaceName(), onDisplaced(), and MortarNodalGeometryOutput::output().

_output_params

std::unique_ptr<InputParameters> AutomaticMortarGeneration::_output_params

private

Storage for the input parameters used by the mortar nodal geometry output.

Definition at line 686 of file AutomaticMortarGeneration.h.

Referenced by initOutput().

_periodic

const bool AutomaticMortarGeneration::_periodic

private

Whether this object will be generating a mortar segment mesh for periodic constraints.

Definition at line 643 of file AutomaticMortarGeneration.h.

Referenced by computeADNodalNormals(), computeNodalGeometry(), and getNormals().

_primary_boundary_subdomain_ids

std::set<SubdomainID> AutomaticMortarGeneration::_primary_boundary_subdomain_ids

private

Definition at line 570 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), and buildNodeToElemMaps().

_primary_ip_sub_ids

std::set<SubdomainID> AutomaticMortarGeneration::_primary_ip_sub_ids

private

All the primary interior parent subdomain IDs associated with the mortar mesh.

Definition at line 608 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), clear(), and primaryIPSubIDs().

_primary_node_and_elem_to_xi1_secondary_elem

std::map<std::tuple<dof_id_type, const Node *, const Elem *>, std::pair<Real, const Elem *> > AutomaticMortarGeneration::_primary_node_and_elem_to_xi1_secondary_elem

private

Same type of container, but for mapping (Primary Node ID, Primary Node, Primary Elem) -> (xi^(1), Secondary Elem) where they are inverse-projected along the nodal normal direction.

Note that the first item of the key, the primary node ID, is important for storing the key-value pairs in a consistent order across processes, e.g. this container has to be ordered!

Definition at line 548 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), clear(), processAlignedNodes(), and projectPrimaryNodesSinglePair().

_primary_requested_boundary_ids

std::set<BoundaryID> AutomaticMortarGeneration::_primary_requested_boundary_ids

private

The boundary ids corresponding to all the primary surfaces.

Definition at line 512 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), and buildNodeToElemMaps().

_primary_secondary_boundary_id_pairs

std::vector<std::pair<BoundaryID, BoundaryID> > AutomaticMortarGeneration::_primary_secondary_boundary_id_pairs

private

A list of primary/secondary boundary id pairs corresponding to each side of the mortar interface.

Definition at line 516 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), buildMortarSegmentMesh(), mortarInterfaceName(), and primarySecondaryBoundaryIDPair().

_primary_secondary_subdomain_id_pairs

std::vector<std::pair<SubdomainID, SubdomainID> > AutomaticMortarGeneration::_primary_secondary_subdomain_id_pairs

private

A list of primary/secondary subdomain id pairs corresponding to each side of the mortar interface.

Definition at line 565 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), buildMortarSegmentMesh3d(), computeInactiveLMElems(), computeInactiveLMNodes(), computeIncorrectEdgeDroppingInactiveLMNodes(), computeMsmStatistics(), projectPrimaryNodes(), and projectSecondaryNodes().

_projected_secondary_nodes

std::unordered_set<dof_id_type> AutomaticMortarGeneration::_projected_secondary_nodes

private

Debugging container for printing information about fraction of successful projections for secondary nodes.

If !_debug then this should always be empty

Definition at line 695 of file AutomaticMortarGeneration.h.

Referenced by clear(), and projectSecondaryNodesSinglePair().

_secondary_boundary_subdomain_ids

std::set<SubdomainID> AutomaticMortarGeneration::_secondary_boundary_subdomain_ids

private

The secondary/primary lower-dimensional boundary subdomain ids are the secondary/primary boundary ids.

Definition at line 569 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), buildCouplingInformation(), buildMortarSegmentMesh(), buildNodeToElemMaps(), computeADNodalNormals(), and computeNodalGeometry().

_secondary_element_to_secondary_lowerd_element

std::unordered_map<dof_id_type, const Elem *> AutomaticMortarGeneration::_secondary_element_to_secondary_lowerd_element

private

Map from full dimensional secondary element id to lower dimensional secondary element.

Definition at line 590 of file AutomaticMortarGeneration.h.

Referenced by clear(), computeNodalGeometry(), and getSecondaryLowerdElemFromSecondaryElem().

_secondary_elems_to_mortar_segments

std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID> > AutomaticMortarGeneration::_secondary_elems_to_mortar_segments

private

We maintain a mapping from lower-dimensional secondary elements in the original mesh to (sets of) elements in mortar_segment_mesh.

This allows us to quickly determine which elements need to be split.

Definition at line 602 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), clear(), computeMsmStatistics(), secondariesToMortarSegments(), and secondariesToMortarSegments().

_secondary_ip_sub_ids

std::set<SubdomainID> AutomaticMortarGeneration::_secondary_ip_sub_ids

private

All the secondary interior parent subdomain IDs associated with the mortar mesh.

Definition at line 605 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), buildMortarSegmentMesh3d(), clear(), and secondaryIPSubIDs().

_secondary_node_and_elem_to_xi2_primary_elem

std::unordered_map<std::pair<const Node *, const Elem *>, std::pair<Real, const Elem *> > AutomaticMortarGeneration::_secondary_node_and_elem_to_xi2_primary_elem

private

Similar to the map above, but associates a (Secondary Node, Secondary Elem) pair to a (xi^(2), primary Elem) pair.

This allows a single secondary node, which is potentially connected to two elements on the secondary side, to be associated with multiple primary Elem/xi^(2) values to handle the case where the primary and secondary nodes are "matching". In this configuration:

A B o--—o--—o (secondary orientation ->) | v ---—x---— (primary orientation <-) C D

The entries in the map should be: (Elem A, Node 1) -> (Elem C, xi^(2)=-1) (Elem B, Node 0) -> (Elem D, xi^(2)=+1)

Definition at line 540 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), clear(), processAlignedNodes(), and projectSecondaryNodesSinglePair().

_secondary_node_to_hh_nodal_tangents

std::unordered_map<const Node *, std::array<Point, 2> > AutomaticMortarGeneration::_secondary_node_to_hh_nodal_tangents

private

Container for storing the nodal tangent/binormal vectors associated with each secondary node (Householder approach).

Definition at line 587 of file AutomaticMortarGeneration.h.

Referenced by clear(), computeNodalGeometry(), getNodalTangents(), and MortarNodalGeometryOutput::output().

_secondary_node_to_nodal_normal

std::unordered_map<const Node *, Point> AutomaticMortarGeneration::_secondary_node_to_nodal_normal

private

Container for storing the nodal normal vector associated with each secondary node.

Definition at line 583 of file AutomaticMortarGeneration.h.

Referenced by clear(), computeADNodalNormals(), computeNodalGeometry(), getNodalNormals(), MortarNodalGeometryOutput::output(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

_secondary_requested_boundary_ids

std::set<BoundaryID> AutomaticMortarGeneration::_secondary_requested_boundary_ids

private

The boundary ids corresponding to all the secondary surfaces.

Definition at line 509 of file AutomaticMortarGeneration.h.

Referenced by AutomaticMortarGeneration(), and buildNodeToElemMaps().

_triangulate_triangles

const bool AutomaticMortarGeneration::_triangulate_triangles

private

Whether already-triangular clipped polygons should still be centroid-subdivided.

Definition at line 677 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d().

_triangulation_mode

const MortarSegmentTriangulationMode AutomaticMortarGeneration::_triangulation_mode

private

Triangulation mode used for clipped 3D mortar polygons.

Definition at line 674 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh3d().

_xi_tolerance

Real AutomaticMortarGeneration::_xi_tolerance = 1e-6

private

Tolerance for checking projection xi values.

Usually we are checking whether we projected onto a certain element (in which case -1 <= xi <= 1) or whether we should have already projected a primary node (in which case we error if abs(xi) is sufficiently close to 1)

Definition at line 658 of file AutomaticMortarGeneration.h.

Referenced by buildMortarSegmentMesh(), projectPrimaryNodesSinglePair(), and projectSecondaryNodesSinglePair().

system_name

const std::string AutomaticMortarGeneration::system_name

static

The name of the nodal normals system.

We store this in one place so it's easy to change later.

Definition at line 140 of file AutomaticMortarGeneration.h.

---

The documentation for this class was generated from the following files:

• include/constraints/AutomaticMortarGeneration.h
• src/constraints/AutomaticMortarGeneration.C