AutomaticMortarGeneration.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "MortarSegmentInfo.h"
#include "Mortar3DSubpatchPlane.h"
#include "MooseHashing.h"
#include "ConsoleStreamInterface.h"
#include "MooseError.h"
#include "MooseUtils.h"

// libMesh includes
#include "libmesh/id_types.h"
#include "libmesh/equation_systems.h"
#include "libmesh/elem.h"
#include "libmesh/int_range.h"

// C++ includes
#include <array>
#include <optional>
#include <set>
#include <memory>
#include <vector>
#include <unordered_map>

// Forward declarations
namespace libMesh
{
class MeshBase;
class System;
}
class GetPot;

// Using statements
using libMesh::boundary_id_type;
using libMesh::CompareDofObjectsByID;
using libMesh::dof_id_type;
using libMesh::Elem;
using libMesh::MeshBase;
using libMesh::Node;
using libMesh::Point;
using libMesh::Real;
using libMesh::subdomain_id_type;

typedef boundary_id_type BoundaryID;
typedef subdomain_id_type SubdomainID;

struct MortarSegmentReferencePoints
{
  std::array<Point, 3> secondary_reference_points;
  std::array<Point, 3> primary_reference_points;
};

class AutomaticMortarGeneration : public ConsoleStreamInterface
{
public:
  const static std::string system_name;

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

  void buildNodeToElemMaps();

  void computeNodalGeometry();

  void projectSecondaryNodes();

  void projectPrimaryNodes();

  void buildMortarSegmentMesh();

  void buildMortarSegmentMesh3d();

  struct MsmSubdomainStats
  {
    SubdomainID primary_subd_id;
    SubdomainID secondary_subd_id;
    std::size_t secondary_lower_n_elems;
    Real secondary_lower_max_volume;
    Real secondary_lower_min_volume;
    Real secondary_lower_median_volume;
    std::size_t primary_lower_n_elems;
    Real primary_lower_max_volume;
    Real primary_lower_min_volume;
    Real primary_lower_median_volume;
    std::size_t msm_n_elems;
    Real msm_max_volume;
    Real msm_min_volume;
    Real msm_median_volume;
  };

  std::vector<MsmSubdomainStats> computeMsmStatistics();

  void msmStatistics();

  void clear();

  void meshChanged() { _msm_node_id_start = std::nullopt; }

  bool onDisplaced() const { return _on_displaced; }

  std::vector<Point> getNodalNormals(const Elem & secondary_elem) const;

  std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
  getNodalTangents(const Elem & secondary_elem) const;

  std::map<unsigned int, unsigned int>
  getSecondaryIpToLowerElementMap(const Elem & lower_secondary_elem) const;

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

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

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

  const Elem * getSecondaryLowerdElemFromSecondaryElem(dof_id_type secondary_elem_id) const;

  void computeInactiveLMNodes();

  void computeIncorrectEdgeDroppingInactiveLMNodes();

  void computeInactiveLMElems();

  const std::unordered_map<dof_id_type, std::unordered_set<dof_id_type>> &
  mortarInterfaceCoupling() const
  {
    return _mortar_interface_coupling;
  }

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

  const MeshBase & mortarSegmentMesh() const { return *_mortar_segment_mesh; }

  const std::unordered_map<const Elem *, MortarSegmentInfo> & mortarSegmentMeshElemToInfo() const
  {
    return _msm_elem_to_info;
  }

  int dim() const { return _mesh.mesh_dimension(); }

  Mortar3DQuadraturePointMapping mortar3DQpMapping() const { return _mortar_3d_qp_mapping; }

  const MortarSegmentReferencePoints &
  mortarSegmentReferencePoints(const Elem & mortar_segment_elem) const;

  const std::unordered_set<const Node *> & getInactiveLMNodes() const
  {
    return _inactive_local_lm_nodes;
  }

  const std::unordered_set<const Elem *> & getInactiveLMElems() const
  {
    return _inactive_local_lm_elems;
  }

  bool incorrectEdgeDropping() const { return !_correct_edge_dropping; }

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

  std::vector<MortarFilterIter> secondariesToMortarSegments(const Node & node) const;

  const std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID>> &
  secondariesToMortarSegments() const
  {
    return _secondary_elems_to_mortar_segments;
  }

  const std::set<SubdomainID> & secondaryIPSubIDs() const { return _secondary_ip_sub_ids; }

  const std::set<SubdomainID> & primaryIPSubIDs() const { return _primary_ip_sub_ids; }

  const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodesToSecondaryElem() const
  {
    return _nodes_to_secondary_elem_map;
  }

  void initOutput();

private:
  void outputMortarMesh();

  std::string mortarInterfaceName() const;

  void buildCouplingInformation();

  MooseApp & _app;

  MeshBase & _mesh;

  std::set<BoundaryID> _secondary_requested_boundary_ids;

  std::set<BoundaryID> _primary_requested_boundary_ids;

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

  std::unordered_map<dof_id_type, std::vector<const Elem *>> _nodes_to_secondary_elem_map;
  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;

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

  std::unique_ptr<MeshBase> _mortar_segment_mesh;

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

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

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

  std::set<SubdomainID> _secondary_boundary_subdomain_ids;
  std::set<SubdomainID> _primary_boundary_subdomain_ids;

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

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

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

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

  // List of inactive lagrange multiplier nodes (for nodal variables)
  std::unordered_set<const Node *> _inactive_local_lm_nodes;

  std::unordered_set<const Elem *> _inactive_local_lm_elems;

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

  std::set<SubdomainID> _secondary_ip_sub_ids;

  std::set<SubdomainID> _primary_ip_sub_ids;

  void projectSecondaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id,
                                       SubdomainID lower_dimensional_secondary_subdomain_id);

  void projectPrimaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id,
                                     SubdomainID lower_dimensional_secondary_subdomain_id);

  void
  householderOrthogolization(const Point & normal, Point & tangent_one, Point & tangent_two) const;

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

  const bool _debug;

  const bool _on_displaced;

  const bool _periodic;

  const bool _distributed;

  Real _newton_tolerance = 1e-12;

  Real _xi_tolerance = 1e-6;

  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;

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

  std::unique_ptr<InputParameters> _output_params;

  std::optional<dof_id_type> _msm_node_id_start;

  std::unordered_set<dof_id_type> _projected_secondary_nodes;

  std::unordered_set<dof_id_type> _failed_secondary_node_projections;

  friend class MortarNodalGeometryOutput;
  friend class AugmentSparsityOnInterface;
};

inline const std::pair<BoundaryID, BoundaryID> &
AutomaticMortarGeneration::primarySecondaryBoundaryIDPair() const
{
  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();
}