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 "MeshAlignment2D3D.h"
      11             : #include "Assembly.h"
      12             : #include "KDTree.h"
      13             : #include "MooseUtils.h"
      14             : 
      15             : #include "libmesh/elem.h"
      16             : 
      17          69 : MeshAlignment2D3D::MeshAlignment2D3D(const MooseMesh & mesh) : MeshAlignmentOneToMany(mesh) {}
      18             : 
      19             : void
      20          61 : MeshAlignment2D3D::initialize(
      21             :     const std::vector<std::tuple<dof_id_type, unsigned short int>> & primary_boundary_info,
      22             :     const std::vector<std::tuple<dof_id_type, unsigned short int>> & secondary_boundary_info,
      23             :     const Point & axis_point,
      24             :     const RealVectorValue & axis_direction)
      25             : {
      26          61 :   extractFromBoundaryInfo(primary_boundary_info,
      27          61 :                           _primary_elem_ids,
      28          61 :                           _primary_side_ids,
      29          61 :                           _primary_elem_points,
      30          61 :                           _primary_node_ids,
      31          61 :                           _primary_node_points);
      32             : 
      33          61 :   extractFromBoundaryInfo(secondary_boundary_info,
      34          61 :                           _secondary_elem_ids,
      35          61 :                           _secondary_side_ids,
      36          61 :                           _secondary_elem_points,
      37          61 :                           _secondary_node_ids,
      38          61 :                           _secondary_node_points);
      39             : 
      40          61 :   buildMapping();
      41          61 :   checkAlignment(axis_point, axis_direction);
      42          61 : }
      43             : 
      44             : void
      45          64 : MeshAlignment2D3D::buildCoupledElemQpIndexMapSecondary(Assembly & assembly)
      46             : {
      47        6464 :   for (const auto i_secondary : index_range(_secondary_elem_ids))
      48             :   {
      49        6400 :     const auto secondary_elem_id = _secondary_elem_ids[i_secondary];
      50        6400 :     const Elem * secondary_elem = _mesh.queryElemPtr(secondary_elem_id);
      51             :     // The PID check is needed to exclude ghost elements, since those don't
      52             :     // necessarily have a coupled element on this processor:
      53        6400 :     if (secondary_elem && secondary_elem->processor_id() == _mesh.processor_id())
      54             :     {
      55        4000 :       const auto secondary_side_id = _secondary_side_ids[i_secondary];
      56             :       assembly.setCurrentSubdomainID(secondary_elem->subdomain_id());
      57        4000 :       assembly.reinit(secondary_elem, secondary_side_id);
      58        4000 :       const auto secondary_qps = assembly.qPointsFace().stdVector();
      59        4000 :       _n_qp_secondary = secondary_qps.size();
      60             : 
      61        4000 :       const auto primary_elem_id = getCoupledPrimaryElemID(secondary_elem_id);
      62        4000 :       auto it = std::find(_primary_elem_ids.begin(), _primary_elem_ids.end(), primary_elem_id);
      63             :       const auto i_primary = std::distance(_primary_elem_ids.begin(), it);
      64        4000 :       const auto primary_side_id = _primary_side_ids[i_primary];
      65             : 
      66             :       // This element should be ghosted, so this should be safe
      67        4000 :       const Elem * primary_elem = _mesh.elemPtr(primary_elem_id);
      68             :       assembly.setCurrentSubdomainID(primary_elem->subdomain_id());
      69        4000 :       assembly.reinit(primary_elem, primary_side_id);
      70        4000 :       auto primary_qps = assembly.qPointsFace().stdVector();
      71        4000 :       _n_qp_primary = primary_qps.size();
      72             : 
      73             :       // compute the area for each quadrature point on the primary side
      74        4000 :       if (_primary_elem_id_to_area.find(primary_elem_id) == _primary_elem_id_to_area.end())
      75             :       {
      76             :         const auto & JxW = assembly.JxWFace();
      77             :         const auto & coord = assembly.coordTransformation();
      78         469 :         std::vector<Real> area(_n_qp_primary, 0.0);
      79        1407 :         for (unsigned int qp = 0; qp < _n_qp_primary; qp++)
      80         938 :           area[qp] = JxW[qp] * coord[qp];
      81         469 :         _primary_elem_id_to_area[primary_elem_id] = area;
      82             :       }
      83             : 
      84        4000 :       _secondary_elem_id_to_qp_indices[secondary_elem_id].resize(secondary_qps.size());
      85        4000 :       KDTree kd_tree_qp(primary_qps, _mesh.getMaxLeafSize());
      86       20000 :       for (const auto i : index_range(secondary_qps))
      87             :       {
      88             :         unsigned int patch_size = 1;
      89       16000 :         std::vector<std::size_t> return_index(patch_size);
      90       16000 :         kd_tree_qp.neighborSearch(secondary_qps[i], patch_size, return_index);
      91       16000 :         _secondary_elem_id_to_qp_indices[secondary_elem_id][i] = return_index[0];
      92             :       }
      93             :     }
      94             :   }
      95          64 : }
      96             : 
      97             : const std::vector<Real> &
      98       90401 : MeshAlignment2D3D::getPrimaryArea(const dof_id_type primary_elem_id) const
      99             : {
     100             :   auto it = _primary_elem_id_to_area.find(primary_elem_id);
     101             :   mooseAssert(it != _primary_elem_id_to_area.end(), "The element ID has no area stored.");
     102       90401 :   return it->second;
     103             : }