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       66959 : GeometricSearchData::GeometricSearchData(SubProblem & subproblem, MooseMesh & mesh)
      24       66959 :   : _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       66959 :   const auto & nodeset_ids = _mesh.meshNodesetIds();
      28      324738 :   for (const auto id : nodeset_ids)
      29      257779 :     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       66959 : }
      33             : 
      34       62455 : GeometricSearchData::~GeometricSearchData()
      35             : {
      36       64325 :   for (auto & it : _penetration_locators)
      37        1870 :     delete it.second;
      38             : 
      39       64416 :   for (auto & it : _nearest_node_locators)
      40        1961 :     delete it.second;
      41       62455 : }
      42             : 
      43             : void
      44      292530 : GeometricSearchData::update(GeometricSearchType type)
      45             : {
      46      292530 :   if (type == ALL || type == QUADRATURE || type == NEAREST_NODE)
      47             :   {
      48      292530 :     if (_first) // Only do this once
      49             :     {
      50       63988 :       _first = false;
      51             : 
      52       64133 :       for (const auto & it : _secondary_to_qsecondary)
      53         145 :         generateQuadratureNodes(it.first, it.second);
      54             : 
      55             :       // reinit on displaced mesh before update
      56       63988 :       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      314135 :     for (const auto & qbnd : _quadrature_boundaries)
      65       21605 :       updateQuadratureNodes(qbnd);
      66             :   }
      67             : 
      68      292530 :   if (type == ALL || type == NEAREST_NODE)
      69             :   {
      70      480765 :     for (const auto & nnl_it : _nearest_node_locators)
      71             :     {
      72      188235 :       NearestNodeLocator * nnl = nnl_it.second;
      73      188235 :       nnl->findNodes();
      74             :     }
      75             :   }
      76             : 
      77      292530 :   if (type == ALL || type == PENETRATION)
      78             :   {
      79      414476 :     for (const auto & pl_it : _penetration_locators)
      80             :     {
      81      185938 :       PenetrationLocator * pl = pl_it.second;
      82      185938 :       pl->detectPenetration();
      83             :     }
      84             :   }
      85             : 
      86      292526 :   if (type == ALL || type == PENETRATION)
      87             :   {
      88      228538 :     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      292526 : }
      95             : 
      96             : void
      97        8220 : GeometricSearchData::reinit()
      98             : {
      99        8220 :   _mesh.clearQuadratureNodes();
     100             :   // Update the position of quadrature nodes first
     101        8282 :   for (const auto & qbnd : _quadrature_boundaries)
     102          62 :     reinitQuadratureNodes(qbnd);
     103             : 
     104        8370 :   for (const auto & nnl_it : _nearest_node_locators)
     105             :   {
     106         150 :     NearestNodeLocator * nnl = nnl_it.second;
     107         150 :     nnl->reinit();
     108             :   }
     109             : 
     110        8346 :   for (const auto & pl_it : _penetration_locators)
     111             :   {
     112         126 :     PenetrationLocator * pl = pl_it.second;
     113         126 :     pl->reinit();
     114             :   }
     115             : 
     116        8220 :   for (const auto & epl_it : _element_pair_locators)
     117             :   {
     118           0 :     ElementPairLocator & epl = *(epl_it.second);
     119           0 :     epl.reinit();
     120             :   }
     121        8220 : }
     122             : 
     123             : void
     124         942 : GeometricSearchData::clearNearestNodeLocators()
     125             : {
     126        1884 :   for (const auto & nnl_it : _nearest_node_locators)
     127             :   {
     128         942 :     NearestNodeLocator * nnl = nnl_it.second;
     129         942 :     nnl->reinit();
     130             :   }
     131         942 : }
     132             : 
     133             : Real
     134         361 : GeometricSearchData::maxPatchPercentage()
     135             : {
     136         361 :   Real max = 0.0;
     137             : 
     138        1083 :   for (const auto & nnl_it : _nearest_node_locators)
     139             :   {
     140         722 :     NearestNodeLocator * nnl = nnl_it.second;
     141             : 
     142         722 :     if (nnl->_max_patch_percentage > max)
     143         349 :       max = nnl->_max_patch_percentage;
     144             :   }
     145             : 
     146         361 :   return max;
     147             : }
     148             : 
     149             : PenetrationLocator &
     150       13619 : GeometricSearchData::getPenetrationLocator(const BoundaryName & primary,
     151             :                                            const BoundaryName & secondary,
     152             :                                            Order order)
     153             : {
     154       13619 :   auto primary_id = _mesh.getBoundaryID(primary);
     155       13619 :   auto secondary_id = _mesh.getBoundaryID(secondary);
     156             : 
     157       13619 :   _subproblem.addGhostedBoundary(primary_id);
     158       13619 :   _subproblem.addGhostedBoundary(secondary_id);
     159             : 
     160             :   PenetrationLocator * pl =
     161       13619 :       _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)];
     162             : 
     163       13619 :   if (!pl)
     164             :   {
     165           0 :     pl = new PenetrationLocator(_subproblem,
     166             :                                 *this,
     167             :                                 _mesh,
     168             :                                 primary_id,
     169             :                                 secondary_id,
     170             :                                 order,
     171        1759 :                                 getNearestNodeLocator(primary_id, secondary_id));
     172             : 
     173        1759 :     if (_search_using_point_locator)
     174          39 :       pl->setUsePointLocator(true);
     175             : 
     176        1759 :     _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)] = pl;
     177             :   }
     178             : 
     179       13619 :   return *pl;
     180             : }
     181             : 
     182             : PenetrationLocator &
     183         129 : GeometricSearchData::getQuadraturePenetrationLocator(const BoundaryName & primary,
     184             :                                                      const BoundaryName & secondary,
     185             :                                                      Order order)
     186             : {
     187         129 :   const auto primary_id = _mesh.getBoundaryID(primary);
     188         129 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     189             : 
     190         129 :   _subproblem.addGhostedBoundary(primary_id);
     191         129 :   _subproblem.addGhostedBoundary(secondary_id);
     192             : 
     193             :   // Generate a new boundary id (we use the negative numbers and subtract 1 to disambiguate id=0)
     194         129 :   const auto qsecondary_id = -secondary_id - 1;
     195             : 
     196         129 :   _secondary_to_qsecondary[secondary_id] = qsecondary_id;
     197             : 
     198             :   PenetrationLocator * pl =
     199         129 :       _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, qsecondary_id)];
     200             : 
     201         129 :   if (!pl)
     202             :   {
     203           0 :     pl = new PenetrationLocator(_subproblem,
     204             :                                 *this,
     205             :                                 _mesh,
     206             :                                 primary_id,
     207             :                                 qsecondary_id,
     208             :                                 order,
     209         119 :                                 getQuadratureNearestNodeLocator(primary_id, secondary_id));
     210             : 
     211         119 :     if (_search_using_point_locator)
     212           0 :       pl->setUsePointLocator(true);
     213             : 
     214         119 :     _penetration_locators[std::pair<BoundaryID, BoundaryID>(primary_id, qsecondary_id)] = pl;
     215             :   }
     216             : 
     217         129 :   return *pl;
     218             : }
     219             : 
     220             : NearestNodeLocator &
     221          84 : GeometricSearchData::getNearestNodeLocator(const BoundaryName & primary,
     222             :                                            const BoundaryName & secondary)
     223             : {
     224          84 :   const auto primary_id = _mesh.getBoundaryID(primary);
     225          84 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     226             : 
     227          84 :   _subproblem.addGhostedBoundary(primary_id);
     228          84 :   _subproblem.addGhostedBoundary(secondary_id);
     229             : 
     230          84 :   return getNearestNodeLocator(primary_id, secondary_id);
     231             : }
     232             : 
     233             : NearestNodeLocator &
     234        1990 : GeometricSearchData::getNearestNodeLocator(const BoundaryID primary_id,
     235             :                                            const BoundaryID secondary_id)
     236             : {
     237             :   NearestNodeLocator * nnl =
     238        1990 :       _nearest_node_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)];
     239             : 
     240        1990 :   _subproblem.addGhostedBoundary(primary_id);
     241        1990 :   _subproblem.addGhostedBoundary(secondary_id);
     242             : 
     243        1990 :   if (!nnl)
     244             :   {
     245        1969 :     nnl = new NearestNodeLocator(_subproblem, _mesh, primary_id, secondary_id);
     246        1969 :     _nearest_node_locators[std::pair<BoundaryID, BoundaryID>(primary_id, secondary_id)] = nnl;
     247             :   }
     248             : 
     249        1990 :   return *nnl;
     250             : }
     251             : 
     252             : NearestNodeLocator &
     253          28 : GeometricSearchData::getQuadratureNearestNodeLocator(const BoundaryName & primary,
     254             :                                                      const BoundaryName & secondary)
     255             : {
     256          28 :   const auto primary_id = _mesh.getBoundaryID(primary);
     257          28 :   const auto secondary_id = _mesh.getBoundaryID(secondary);
     258             : 
     259          28 :   _subproblem.addGhostedBoundary(primary_id);
     260          28 :   _subproblem.addGhostedBoundary(secondary_id);
     261             : 
     262          28 :   return getQuadratureNearestNodeLocator(primary_id, secondary_id);
     263             : }
     264             : 
     265             : NearestNodeLocator &
     266         147 : 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         147 :   const auto qsecondary_id = -secondary_id - 1;
     271             : 
     272         147 :   _secondary_to_qsecondary[secondary_id] = qsecondary_id;
     273         147 :   return getNearestNodeLocator(primary_id, qsecondary_id);
     274             : }
     275             : 
     276             : void
     277         207 : 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         207 :   if (_quadrature_boundaries.find(secondary_id) != _quadrature_boundaries.end())
     283             :   {
     284          62 :     if (!reiniting)
     285           0 :       return;
     286             :   }
     287             :   else
     288         145 :     _quadrature_boundaries.insert(secondary_id);
     289             : 
     290         207 :   const MooseArray<Point> & points_face = _subproblem.assembly(0, 0).qPointsFace();
     291             : 
     292         207 :   ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
     293       37763 :   for (const auto & belem : range)
     294             :   {
     295       37556 :     const Elem * elem = belem->_elem;
     296       37556 :     const auto side = belem->_side;
     297       37556 :     const auto boundary_id = belem->_bnd_id;
     298             : 
     299       37556 :     if (elem->processor_id() == _subproblem.processor_id())
     300             :     {
     301       28844 :       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        5466 :         _subproblem.setCurrentSubdomainID(elem, 0);
     306        5466 :         _subproblem.assembly(0, 0).reinit(elem, side);
     307             : 
     308       11556 :         for (const auto qp : index_range(points_face))
     309        6090 :           _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          42 : GeometricSearchData::setSearchUsingPointLocator(bool state)
     324             : {
     325          42 :   if (state != _search_using_point_locator)
     326             :   {
     327          39 :     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          42 :   _search_using_point_locator = state;
     336          42 : }
     337             : 
     338             : void
     339       21605 : GeometricSearchData::updateQuadratureNodes(const BoundaryID secondary_id)
     340             : {
     341       21605 :   const MooseArray<Point> & points_face = _subproblem.assembly(0, 0).qPointsFace();
     342             : 
     343       21605 :   ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
     344     4291442 :   for (const auto & belem : range)
     345             :   {
     346     4269837 :     const Elem * elem = belem->_elem;
     347     4269837 :     const auto side = belem->_side;
     348     4269837 :     const auto boundary_id = belem->_bnd_id;
     349             : 
     350     4269837 :     if (elem->processor_id() == _subproblem.processor_id())
     351             :     {
     352     3403868 :       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       45178 :         _subproblem.setCurrentSubdomainID(elem, 0);
     357       45178 :         _subproblem.assembly(0, 0).reinit(elem, side);
     358             : 
     359      125850 :         for (const auto qp : index_range(points_face))
     360       80672 :           (*_mesh.getQuadratureNode(elem, side, qp)) = points_face[qp];
     361             :       }
     362             :     }
     363             :   }
     364       21605 : }
     365             : 
     366             : void
     367          62 : GeometricSearchData::reinitQuadratureNodes(const BoundaryID /*secondary_id*/)
     368             : {
     369             :   // Regenerate the quadrature nodes
     370         124 :   for (const auto & it : _secondary_to_qsecondary)
     371          62 :     generateQuadratureNodes(it.first, it.second, /*reiniting=*/true);
     372          62 : }
     373             : 
     374             : void
     375         796 : GeometricSearchData::updateGhostedElems()
     376             : {
     377        1518 :   for (const auto & nnl_it : _nearest_node_locators)
     378             :   {
     379         722 :     NearestNodeLocator * nnl = nnl_it.second;
     380         722 :     nnl->updateGhostedElems();
     381             :   }
     382         796 : }