Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "MeshAlignment.h"
      11             : #include "Assembly.h"
      12             : #include "KDTree.h"
      13             : 
      14             : #include "libmesh/elem.h"
      15             : 
      16         454 : MeshAlignment::MeshAlignment(const MooseMesh & mesh, bool require_same_translation)
      17             :   : MeshAlignmentBase(mesh),
      18         454 :     _require_same_translation(require_same_translation),
      19         454 :     _meshes_are_coincident(false),
      20         454 :     _meshes_have_same_translation(false)
      21             : {
      22         454 : }
      23             : 
      24             : void
      25         296 : MeshAlignment::initialize(
      26             :     const std::vector<dof_id_type> & primary_elem_ids,
      27             :     const std::vector<std::tuple<dof_id_type, unsigned short int>> & secondary_boundary_info)
      28             : {
      29         296 :   _primary_elem_ids = primary_elem_ids;
      30         296 :   extractFrom1DElements(
      31         296 :       _primary_elem_ids, _primary_elem_points, _primary_node_ids, _primary_node_points);
      32             : 
      33         296 :   extractFromBoundaryInfo(secondary_boundary_info,
      34         296 :                           _secondary_elem_ids,
      35         296 :                           _secondary_side_ids,
      36         296 :                           _secondary_elem_points,
      37         296 :                           _secondary_node_ids,
      38         296 :                           _secondary_node_points);
      39             : 
      40         296 :   buildMapping();
      41         296 : }
      42             : 
      43             : void
      44         146 : MeshAlignment::initialize(
      45             :     const std::vector<std::tuple<dof_id_type, unsigned short int>> & primary_boundary_info,
      46             :     const std::vector<std::tuple<dof_id_type, unsigned short int>> & secondary_boundary_info)
      47             : {
      48         146 :   extractFromBoundaryInfo(primary_boundary_info,
      49         146 :                           _primary_elem_ids,
      50         146 :                           _primary_side_ids,
      51         146 :                           _primary_elem_points,
      52         146 :                           _primary_node_ids,
      53         146 :                           _primary_node_points);
      54             : 
      55         146 :   extractFromBoundaryInfo(secondary_boundary_info,
      56         146 :                           _secondary_elem_ids,
      57         146 :                           _secondary_side_ids,
      58         146 :                           _secondary_elem_points,
      59         146 :                           _secondary_node_ids,
      60         146 :                           _secondary_node_points);
      61             : 
      62         146 :   buildMapping();
      63         146 : }
      64             : 
      65             : void
      66         442 : MeshAlignment::buildMapping()
      67             : {
      68         442 :   _meshes_are_coincident = true;
      69         442 :   std::vector<unsigned int> primary_elem_pairing_count(_primary_elem_ids.size(), 0);
      70             : 
      71             :   // Build the element mapping
      72         442 :   if (_primary_elem_points.size() > 0 && _secondary_elem_points.size() > 0)
      73             :   {
      74             :     // find the primary elements that are nearest to the secondary elements
      75         442 :     KDTree kd_tree(_primary_elem_points, _mesh.getMaxLeafSize());
      76        6232 :     for (std::size_t i_secondary = 0; i_secondary < _secondary_elem_points.size(); i_secondary++)
      77             :     {
      78             :       unsigned int patch_size = 1;
      79        5790 :       std::vector<std::size_t> return_index(patch_size);
      80        5790 :       kd_tree.neighborSearch(_secondary_elem_points[i_secondary], patch_size, return_index);
      81        5790 :       const std::size_t i_primary = return_index[0];
      82             : 
      83             :       // Flip flag if any pair of points are not coincident
      84        5790 :       if (!_secondary_elem_points[i_secondary].absolute_fuzzy_equals(
      85             :               _primary_elem_points[i_primary]))
      86        3173 :         _meshes_are_coincident = false;
      87             : 
      88        5790 :       const auto primary_elem_id = _primary_elem_ids[i_primary];
      89        5790 :       const auto secondary_elem_id = _secondary_elem_ids[i_secondary];
      90             : 
      91        5790 :       _coupled_elem_ids.insert({primary_elem_id, secondary_elem_id});
      92        5790 :       _coupled_elem_ids.insert({secondary_elem_id, primary_elem_id});
      93             : 
      94        5790 :       primary_elem_pairing_count[i_primary]++;
      95             :     }
      96             :   }
      97             : 
      98             :   // Build the node mapping
      99             :   Point reference_translation_vector;
     100         442 :   _meshes_have_same_translation = true;
     101         442 :   if (_primary_node_points.size() > 0 && _secondary_node_points.size() > 0)
     102             :   {
     103             :     // find the primary nodes that are nearest to the secondary nodes
     104         442 :     KDTree kd_tree(_primary_node_points, _mesh.getMaxLeafSize());
     105        6674 :     for (std::size_t i_secondary = 0; i_secondary < _secondary_node_points.size(); i_secondary++)
     106             :     {
     107             :       unsigned int patch_size = 1;
     108        6232 :       std::vector<std::size_t> return_index(patch_size);
     109        6232 :       kd_tree.neighborSearch(_secondary_node_points[i_secondary], patch_size, return_index);
     110        6232 :       const std::size_t i_primary = return_index[0];
     111             : 
     112             :       const Point translation_vector =
     113             :           _primary_node_points[i_primary] - _secondary_node_points[i_secondary];
     114             : 
     115        6232 :       if (i_secondary == 0)
     116         442 :         reference_translation_vector = translation_vector;
     117             :       else
     118             :       {
     119        5790 :         if (!MooseUtils::absoluteFuzzyEqual(
     120        5790 :                 (reference_translation_vector - translation_vector).norm(), 0))
     121          66 :           _meshes_have_same_translation = false;
     122             :       }
     123             : 
     124        6232 :       const auto primary_node_id = _primary_node_ids[i_primary];
     125        6232 :       const auto secondary_node_id = _secondary_node_ids[i_secondary];
     126             : 
     127        6232 :       _coupled_node_ids.insert({primary_node_id, secondary_node_id});
     128        6232 :       _coupled_node_ids.insert({secondary_node_id, primary_node_id});
     129             :     }
     130             : 
     131         876 :     if (_meshes_have_same_translation &&
     132         434 :         MooseUtils::absoluteFuzzyEqual(reference_translation_vector.norm(), 0))
     133         207 :       _meshes_are_coincident = true;
     134             :     else
     135         235 :       _meshes_are_coincident = false;
     136             : 
     137             :     // Check if meshes are aligned
     138         442 :     if (_require_same_translation)
     139         296 :       _meshes_are_aligned = _meshes_have_same_translation;
     140             :     else
     141             :     {
     142             :       // Else require that all primary boundary elements have exactly one pairing
     143         146 :       _meshes_are_aligned = true;
     144        1426 :       for (std::size_t i_primary = 0; i_primary < _primary_elem_ids.size(); i_primary++)
     145        1280 :         if (primary_elem_pairing_count[i_primary] != 1)
     146          24 :           _meshes_are_aligned = false;
     147             :     }
     148             :   }
     149         442 : }
     150             : 
     151             : void
     152         374 : MeshAlignment::buildCoupledElemQpIndexMap(Assembly & assembly)
     153             : {
     154             :   // get local quadrature point maps on each side
     155             :   std::map<dof_id_type, std::vector<Point>> primary_qp_map =
     156         374 :       getLocalQuadraturePointMap(assembly, _primary_elem_ids, _primary_side_ids);
     157             :   std::map<dof_id_type, std::vector<Point>> secondary_qp_map =
     158         374 :       getLocalQuadraturePointMap(assembly, _secondary_elem_ids, _secondary_side_ids);
     159             : 
     160             :   // build mapping
     161        6214 :   for (const auto & primary_elem_id : _primary_elem_ids)
     162             :   {
     163        5840 :     const Elem * primary_elem = _mesh.queryElemPtr(primary_elem_id);
     164             :     // The PID check is needed to exclude ghost elements, since those don't
     165             :     // necessarily have a coupled element on this processor:
     166        5840 :     if (primary_elem && primary_elem->processor_id() == _mesh.processor_id())
     167             :     {
     168        4012 :       std::vector<Point> & primary_qps = primary_qp_map[primary_elem_id];
     169             : 
     170        4012 :       const dof_id_type secondary_elem_id = getCoupledElemID(primary_elem_id);
     171        4012 :       std::vector<Point> & secondary_qps = secondary_qp_map[secondary_elem_id];
     172             : 
     173             :       mooseAssert(primary_qps.size() == secondary_qps.size(),
     174             :                   "The numbers of quadrature points for each element must be the same");
     175             : 
     176        4012 :       _coupled_elem_qp_indices[primary_elem_id].resize(primary_qps.size());
     177        4012 :       KDTree kd_tree_qp(secondary_qps, _mesh.getMaxLeafSize());
     178       12036 :       for (std::size_t i = 0; i < primary_qps.size(); i++)
     179             :       {
     180             :         unsigned int patch_size = 1;
     181        8024 :         std::vector<std::size_t> return_index(patch_size);
     182        8024 :         kd_tree_qp.neighborSearch(primary_qps[i], patch_size, return_index);
     183        8024 :         _coupled_elem_qp_indices[primary_elem_id][i] = return_index[0];
     184             :       }
     185             :     }
     186             :   }
     187         374 : }
     188             : 
     189             : std::map<dof_id_type, std::vector<Point>>
     190         748 : MeshAlignment::getLocalQuadraturePointMap(Assembly & assembly,
     191             :                                           const std::vector<dof_id_type> & elem_ids,
     192             :                                           const std::vector<unsigned short int> & side_ids) const
     193             : {
     194             :   std::map<dof_id_type, std::vector<Point>> elem_id_to_qps;
     195       12428 :   for (unsigned int i = 0; i < elem_ids.size(); i++)
     196             :   {
     197       11680 :     const auto elem_id = elem_ids[i];
     198       11680 :     const Elem * elem = _mesh.queryElemPtr(elem_id);
     199       11680 :     if (elem)
     200             :     {
     201             :       assembly.setCurrentSubdomainID(elem->subdomain_id());
     202             : 
     203             :       MooseArray<Point> q_points;
     204       10673 :       if (side_ids.size() > 0)
     205             :       {
     206        5387 :         assembly.reinit(elem, side_ids[i]);
     207        5387 :         q_points = assembly.qPointsFace();
     208             :       }
     209             :       else
     210             :       {
     211        5286 :         assembly.reinit(elem);
     212        5286 :         q_points = assembly.qPoints();
     213             :       }
     214             : 
     215       32019 :       for (std::size_t i = 0; i < q_points.size(); i++)
     216       21346 :         elem_id_to_qps[elem_id].push_back(q_points[i]);
     217             :     }
     218             :   }
     219             : 
     220         748 :   return elem_id_to_qps;
     221             : }
     222             : 
     223             : bool
     224        5724 : MeshAlignment::hasCoupledElemID(const dof_id_type & elem_id) const
     225             : {
     226        5724 :   return _coupled_elem_ids.find(elem_id) != _coupled_elem_ids.end();
     227             : }
     228             : 
     229             : const dof_id_type &
     230      171731 : MeshAlignment::getCoupledElemID(const dof_id_type & elem_id) const
     231             : {
     232             :   mooseAssert(hasCoupledElemID(elem_id), "The element ID has no coupled element.");
     233      171731 :   return _coupled_elem_ids.find(elem_id)->second;
     234             : }
     235             : 
     236             : bool
     237     1536000 : MeshAlignment::hasCoupledNodeID(const dof_id_type & node_id) const
     238             : {
     239     1536000 :   return _coupled_node_ids.find(node_id) != _coupled_node_ids.end();
     240             : }
     241             : 
     242             : const dof_id_type &
     243     2389316 : MeshAlignment::getCoupledNodeID(const dof_id_type & node_id) const
     244             : {
     245             :   mooseAssert(hasCoupledNodeID(node_id), "The node ID has no coupled node.");
     246     2389316 :   return _coupled_node_ids.find(node_id)->second;
     247             : }
     248             : 
     249             : unsigned int
     250      323990 : MeshAlignment::getCoupledElemQpIndex(const dof_id_type & elem_id, const unsigned int & qp) const
     251             : {
     252             :   auto it = _coupled_elem_qp_indices.find(elem_id);
     253             :   mooseAssert(it != _coupled_elem_qp_indices.end(),
     254             :               "The element ID has no coupled quadrature point indices.");
     255             :   mooseAssert(qp < it->second.size(), "The quadrature index does not exist in the map.");
     256      323990 :   return it->second[qp];
     257             : }