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"
#include "libmesh/vector_value.h"
 
#include "metaphysicl/raw_type.h"
 
// C++ includes
#include <array>
#include <cmath>
#include <functional>
#include <optional>
#include <set>
#include <memory>
#include <vector>
#include <unordered_map>
 
// Forward declarations
class GetPot;
 
// Using statements
using libMesh::CompareDofObjectsByID;
 
typedef boundary_id_type BoundaryID;
typedef subdomain_id_type SubdomainID;
 
namespace Moose
{
namespace Mortar
{
template <typename Vector>
std::array<Vector, 2>
householderTangents(const Vector & normal)
{
  mooseAssert(MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(normal.norm()), 1),
              "The input nodal normal should have unity norm");
 
  const Vector h_vector(normal(0) + 1.0, normal(1), normal(2));
 
  // The h = n + e_x chart is singular at (-1, 0, 0), where the complementary h2 = n - e_x chart
  // is regular. The h2 frame has determinant -1, so negating its second tangent gives the same
  // positive orientation as the primary chart. At the pole, the frame is
  // {(0, 1, 0), (0, 0, -1)} with nonzero derivatives.
  if (std::abs(MetaPhysicL::raw_value(h_vector(0))) < TOLERANCE)
  {
    const Vector h2_vector(normal(0) - 1.0, normal(1), normal(2));
    const auto h2_squared = h2_vector.norm_sq();
    return {{Vector(-2.0 * h2_vector(0) * h2_vector(1) / h2_squared,
                    1.0 - 2.0 * h2_vector(1) * h2_vector(1) / h2_squared,
                    -2.0 * h2_vector(1) * h2_vector(2) / h2_squared),
             Vector(2.0 * h2_vector(0) * h2_vector(2) / h2_squared,
                    2.0 * h2_vector(1) * h2_vector(2) / h2_squared,
                    2.0 * h2_vector(2) * h2_vector(2) / h2_squared - 1.0)}};
  }
 
  // The squared norm avoids propagating derivatives through a square root.
  const auto h_squared = h_vector.norm_sq();
  return {{Vector(-2.0 * h_vector(0) * h_vector(1) / h_squared,
                  1.0 - 2.0 * h_vector(1) * h_vector(1) / h_squared,
                  -2.0 * h_vector(1) * h_vector(2) / h_squared),
           Vector(-2.0 * h_vector(0) * h_vector(2) / h_squared,
                  -2.0 * h_vector(1) * h_vector(2) / h_squared,
                  1.0 - 2.0 * h_vector(2) * h_vector(2) / h_squared)}};
}
}
}
 
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;
 
  void
  computeADNodalNormals(const std::function<ADPoint(const Node &)> & coordinate,
                        std::unordered_map<const Node *, ADRealVectorValue> & nodal_normals) 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;
 
  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);
 
  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();
}