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 "GeometricSearchData.h"
      11             : // Moose includes
      12             : #include "NearestNodeLocator.h"
      13             : #include "PenetrationLocator.h"
      14             : #include "ElementPairLocator.h"
      15             : #include "SubProblem.h"
      16             : #include "MooseMesh.h"
      17             : #include "Assembly.h"
      18             : 
      19             : #include "libmesh/elem.h"
      20             : #include "libmesh/node.h"
      21             : #include "libmesh/int_range.h"
      22             : 
      23       63665 : GeometricSearchData::GeometricSearchData(SubProblem & subproblem, MooseMesh & mesh)
      24       63665 :   : _subproblem(subproblem), _mesh(mesh), _first(true), _search_using_point_locator(false)
      25             : {
      26             :   // do a check on the boundary IDs to see if any conflict will rise from computing QP node set IDs
      27       63665 :   const auto & nodeset_ids = _mesh.meshNodesetIds();
      28      306706 :   for (const auto id : nodeset_ids)
      29      243041 :     if (nodeset_ids.find(-id - 1) != nodeset_ids.end())
      30           0 :       mooseError("Your mesh contains nodesets with negative IDs that interfere with QP nodeset IDs "
      31             :                  "potentially generated by the GeometricSearch system.");
      32       63665 : }
      33             : 
      34       60555 : GeometricSearchData::~GeometricSearchData()
      35             : {
      36       62277 :   for (auto & it : _penetration_locators)
      37        1722 :     delete it.second;
      38             : 
      39       62361 :   for (auto & it : _nearest_node_locators)
      40        1806 :     delete it.second;
      41       60555 : }
      42             : 
      43             : void
      44      271346 : GeometricSearchData::update(GeometricSearchType type)
      45             : {
      46      271346 :   if (type == ALL || type == QUADRATURE || type == NEAREST_NODE)
      47             :   {
      48      271346 :     if (_first) // Only do this once
      49             :     {
      50       61179 :       _first = false;
      51             : 
      52       61312 :       for (const auto & it : _secondary_to_qsecondary)
      53         133 :         generateQuadratureNodes(it.first, it.second);
      54             : 
      55             :       // reinit on displaced mesh before update
      56       61179 :       for (const auto & epl_it : _element_pair_locators)
      57             :       {
      58           0 :         ElementPairLocator & epl = *(epl_it.second);
      59           0 :         epl.reinit();
      60             :       }
      61             :     }
      62             : 
      63             :     // Update the position of quadrature nodes first
      64      290919 :     for (const auto & qbnd : _quadrature_boundaries)
      65       19573 :       updateQuadratureNodes(qbnd);
      66             :   }
      67             : 
      68      271346 :   if (type == ALL || type == NEAREST_NODE)
      69             :   {
      70      443313 :     for (const auto & nnl_it : _nearest_node_locators)
      71             :     {
      72      171967 :       NearestNodeLocator * nnl = nnl_it.second;
      73      171967 :       nnl->findNodes();
      74             :     }
      75             :   }
      76             : 
      77      271346 :   if (type == ALL || type == PENETRATION)
      78             :   {
      79      380016 :     for (const auto & pl_it : _penetration_locators)
      80             :     {
      81      169852 :       PenetrationLocator * pl = pl_it.second;
      82      169852 :       pl->detectPenetration();
      83             :     }
      84             :   }
      85             : 
      86      271343 :   if (type == ALL || type == PENETRATION)
      87             :   {
      88      210164 :     for (auto & elem_pair_locator_pair : _element_pair_locators)
      89             :     {
      90           0 :       ElementPairLocator & epl = (*elem_pair_locator_pair.second);
      91           0 :       epl.update();
      92             :     }
      93             :   }
      94      271343 : }
      95             : 
      96             : void
      97        7567 : GeometricSearchData::reinit()
      98             : {
      99        7567 :   _mesh.clearQuadratureNodes();
     100             :   // Update the position of quadrature nodes first
     101        7624 :   for (const auto & qbnd : _quadrature_boundaries)
     102          57 :     reinitQuadratureNodes(qbnd);
     103             : 
     104        7705 :   for (const auto & nnl_it : _nearest_node_locators)
     105             :   {
     106         138 :     NearestNodeLocator * nnl = nnl_it.second;
     107         138 :     nnl->reinit();
     108             :   }
     109             : 
     110        7683 :   for (const auto & pl_it : _penetration_locators)
     111             :   {
     112         116 :     PenetrationLocator * pl = pl_it.second;
     113         116 :     pl->reinit();
     114             :   }
     115             : 
     116        7567 :   for (const auto & epl_it : _element_pair_locators)
     117             :   {
     118           0 :     ElementPairLocator & epl = *(epl_it.second);
     119           0 :     epl.reinit();
     120             :   }
     121        7567 : }
     122             : 
     123             : void
     124         864 : GeometricSearchData::clearNearestNodeLocators()
     125             : {
     126        1728 :   for (const auto & nnl_it : _nearest_node_locators)
     127             :   {
     128         864 :     NearestNodeLocator * nnl = nnl_it.second;
     129         864 :     nnl->reinit();
     130             :   }
     131         864 : }
     132             : 
     133             : Real
     134         331 : GeometricSearchData::maxPatchPercentage()
     135             : {
     136         331 :   Real max = 0.0;
     137             : 
     138         993 :   for (const auto & nnl_it : _nearest_node_locators)
     139             :   {
     140         662 :     NearestNodeLocator * nnl = nnl_it.second;
     141             : 
     142         662 :     if (nnl->_max_patch_percentage > max)
     143         320 :       max = nnl->_max_patch_percentage;
     144             :   }
     145             : 
     146         331 :   return max;
     147             : }
     148             : 
     149             : PenetrationLocator &
     150       12634 : GeometricSearchData::getPenetrationLocator(const BoundaryName & primary,
     151             :                                            const BoundaryName & secondary,
     152             :                                            Order order)
     153             : {
     154       12634 :   auto primary_id = _mesh.getBoundaryID(primary);
     155       12634 :   auto secondary_id = _mesh.getBoundaryID(secondary);
     156             : 
     157       12634 :   _subproblem.addGhostedBoundary(primary_id);
     158       12634 :   _subproblem.addGhostedBoundary(secondary_id);
     159             : 
     160             :   PenetrationLocator * pl =
     161       12634 :       _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)];
     162             : 
     163       12634 :   if (!pl)
     164             :   {
     165           0 :     pl = new PenetrationLocator(_subproblem,
     166             :                                 *this,
     167             :                                 _mesh,
     168             :                                 primary_id,
     169             :                                 secondary_id,
     170             :                                 order,
     171        1619 :                                 getNearestNodeLocator(primary_id, secondary_id));
     172             : 
     173        1619 :     if (_search_using_point_locator)
     174          36 :       pl->setUsePointLocator(true);
     175             : 
     176        1619 :     _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)] = pl;
     177             :   }
     178             : 
     179       12634 :   return *pl;
     180             : }
     181             : 
     182             : PenetrationLocator &
     183         119 : GeometricSearchData::getQuadraturePenetrationLocator(const BoundaryName & primary,
     184             :                                                      const BoundaryName & secondary,
     185             :                                                      Order order)
     186             : {
     187         119 :   const auto primary_id = _mesh.getBoundaryID(primary);
     188         119 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     189             : 
     190         119 :   _subproblem.addGhostedBoundary(primary_id);
     191         119 :   _subproblem.addGhostedBoundary(secondary_id);
     192             : 
     193             :   // Generate a new boundary id (we use the negative numbers and subtract 1 to disambiguate id=0)
     194         119 :   const auto qsecondary_id = -secondary_id - 1;
     195             : 
     196         119 :   _secondary_to_qsecondary[secondary_id] = qsecondary_id;
     197             : 
     198             :   PenetrationLocator * pl =
     199         119 :       _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, qsecondary_id)];
     200             : 
     201         119 :   if (!pl)
     202             :   {
     203           0 :     pl = new PenetrationLocator(_subproblem,
     204             :                                 *this,
     205             :                                 _mesh,
     206             :                                 primary_id,
     207             :                                 qsecondary_id,
     208             :                                 order,
     209         109 :                                 getQuadratureNearestNodeLocator(primary_id, secondary_id));
     210             : 
     211         109 :     if (_search_using_point_locator)
     212           0 :       pl->setUsePointLocator(true);
     213             : 
     214         109 :     _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, qsecondary_id)] = pl;
     215             :   }
     216             : 
     217         119 :   return *pl;
     218             : }
     219             : 
     220             : NearestNodeLocator &
     221          78 : GeometricSearchData::getNearestNodeLocator(const BoundaryName & primary,
     222             :                                            const BoundaryName & secondary)
     223             : {
     224          78 :   const auto primary_id = _mesh.getBoundaryID(primary);
     225          78 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     226             : 
     227          78 :   _subproblem.addGhostedBoundary(primary_id);
     228          78 :   _subproblem.addGhostedBoundary(secondary_id);
     229             : 
     230          78 :   return getNearestNodeLocator(primary_id, secondary_id);
     231             : }
     232             : 
     233             : NearestNodeLocator &
     234        1832 : GeometricSearchData::getNearestNodeLocator(const BoundaryID primary_id,
     235             :                                            const BoundaryID secondary_id)
     236             : {
     237             :   NearestNodeLocator * nnl =
     238        1832 :       _nearest_node_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)];
     239             : 
     240        1832 :   _subproblem.addGhostedBoundary(primary_id);
     241        1832 :   _subproblem.addGhostedBoundary(secondary_id);
     242             : 
     243        1832 :   if (!nnl)
     244             :   {
     245        1812 :     nnl = new NearestNodeLocator(_subproblem, _mesh, primary_id, secondary_id);
     246        1812 :     _nearest_node_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)] = nnl;
     247             :   }
     248             : 
     249        1832 :   return *nnl;
     250             : }
     251             : 
     252             : NearestNodeLocator &
     253          26 : GeometricSearchData::getQuadratureNearestNodeLocator(const BoundaryName & primary,
     254             :                                                      const BoundaryName & secondary)
     255             : {
     256          26 :   const auto primary_id = _mesh.getBoundaryID(primary);
     257          26 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     258             : 
     259          26 :   _subproblem.addGhostedBoundary(primary_id);
     260          26 :   _subproblem.addGhostedBoundary(secondary_id);
     261             : 
     262          26 :   return getQuadratureNearestNodeLocator(primary_id, secondary_id);
     263             : }
     264             : 
     265             : NearestNodeLocator &
     266         135 : GeometricSearchData::getQuadratureNearestNodeLocator(const BoundaryID primary_id,
     267             :                                                      const BoundaryID secondary_id)
     268             : {
     269             :   // Generate a new boundary id (we use the negative numbers and subtract 1 to disambiguate id=0)
     270         135 :   const auto qsecondary_id = -secondary_id - 1;
     271             : 
     272         135 :   _secondary_to_qsecondary[secondary_id] = qsecondary_id;
     273         135 :   return getNearestNodeLocator(primary_id, qsecondary_id);
     274             : }
     275             : 
     276             : void
     277         190 : GeometricSearchData::generateQuadratureNodes(const BoundaryID secondary_id,
     278             :                                              const BoundaryID qsecondary_id,
     279             :                                              bool reiniting)
     280             : {
     281             :   // Have we already generated quadrature nodes for this boundary id?
     282         190 :   if (_quadrature_boundaries.find(secondary_id) != _quadrature_boundaries.end())
     283             :   {
     284          57 :     if (!reiniting)
     285           0 :       return;
     286             :   }
     287             :   else
     288         133 :     _quadrature_boundaries.insert(secondary_id);
     289             : 
     290         190 :   const MooseArray<Point> & points_face = _subproblem.assembly(0, 0).qPointsFace();
     291             : 
     292         190 :   ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
     293       34473 :   for (const auto & belem : range)
     294             :   {
     295       34283 :     const Elem * elem = belem->_elem;
     296       34283 :     const auto side = belem->_side;
     297       34283 :     const auto boundary_id = belem->_bnd_id;
     298             : 
     299       34283 :     if (elem->processor_id() == _subproblem.processor_id())
     300             :     {
     301       25571 :       if (boundary_id == (BoundaryID)secondary_id)
     302             :       {
     303             :         // All we should need to do here is reinit the underlying libMesh::FE object because that
     304             :         // will get us the correct points for the current element and side
     305        4797 :         _subproblem.setCurrentSubdomainID(elem, 0);
     306        4797 :         _subproblem.assembly(0, 0).reinit(elem, side);
     307             : 
     308       10154 :         for (const auto qp : index_range(points_face))
     309        5357 :           _mesh.addQuadratureNode(elem, side, qp, qsecondary_id, points_face[qp]);
     310             :       }
     311             :     }
     312             :   }
     313             : }
     314             : 
     315             : void
     316           0 : GeometricSearchData::addElementPairLocator(const BoundaryID interface_id,
     317             :                                            std::shared_ptr<ElementPairLocator> epl)
     318             : {
     319           0 :   _element_pair_locators[interface_id] = epl;
     320           0 : }
     321             : 
     322             : void
     323          39 : GeometricSearchData::setSearchUsingPointLocator(bool state)
     324             : {
     325          39 :   if (state != _search_using_point_locator)
     326             :   {
     327          36 :     for (auto & [key, val] : _penetration_locators)
     328             :     {
     329           0 :       libmesh_ignore(key);
     330           0 :       if (val)
     331           0 :         val->setUsePointLocator(state);
     332             :     }
     333             :   }
     334             : 
     335          39 :   _search_using_point_locator = state;
     336          39 : }
     337             : 
     338             : void
     339       19573 : GeometricSearchData::updateQuadratureNodes(const BoundaryID secondary_id)
     340             : {
     341       19573 :   const MooseArray<Point> & points_face = _subproblem.assembly(0, 0).qPointsFace();
     342             : 
     343       19573 :   ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
     344     3870502 :   for (const auto & belem : range)
     345             :   {
     346     3850929 :     const Elem * elem = belem->_elem;
     347     3850929 :     const auto side = belem->_side;
     348     3850929 :     const auto boundary_id = belem->_bnd_id;
     349             : 
     350     3850929 :     if (elem->processor_id() == _subproblem.processor_id())
     351             :     {
     352     2984960 :       if (boundary_id == (BoundaryID)secondary_id)
     353             :       {
     354             :         // All we should need to do here is reinit the underlying libMesh::FE object because that
     355             :         // will get us the correct points for the current element and side
     356       39672 :         _subproblem.setCurrentSubdomainID(elem, 0);
     357       39672 :         _subproblem.assembly(0, 0).reinit(elem, side);
     358             : 
     359      110542 :         for (const auto qp : index_range(points_face))
     360       70870 :           (*_mesh.getQuadratureNode(elem, side, qp)) = points_face[qp];
     361             :       }
     362             :     }
     363             :   }
     364       19573 : }
     365             : 
     366             : void
     367          57 : GeometricSearchData::reinitQuadratureNodes(const BoundaryID /*secondary_id*/)
     368             : {
     369             :   // Regenerate the quadrature nodes
     370         114 :   for (const auto & it : _secondary_to_qsecondary)
     371          57 :     generateQuadratureNodes(it.first, it.second, /*reiniting=*/true);
     372          57 : }
     373             : 
     374             : void
     375         730 : GeometricSearchData::updateGhostedElems()
     376             : {
     377        1392 :   for (const auto & nnl_it : _nearest_node_locators)
     378             :   {
     379         662 :     NearestNodeLocator * nnl = nnl_it.second;
     380         662 :     nnl->updateGhostedElems();
     381             :   }
     382         730 : }