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 "Positions.h"
      11             : #include "libmesh/parallel_algebra.h"
      12             : #include "MooseMeshUtils.h"
      13             : 
      14             : InputParameters
      15      188245 : Positions::validParams()
      16             : {
      17      188245 :   InputParameters params = GeneralReporter::validParams();
      18             : 
      19      188245 :   params.addParam<PositionsName>("initial_positions",
      20             :                                  "Positions at the beginning of the simulation");
      21             : 
      22             :   // This parameter should be set by each derived class depending on whether the generation of
      23             :   // positions is replicated or distributed. We want positions to be replicated across all ranks
      24      188245 :   params.addRequiredParam<bool>("auto_broadcast",
      25             :                                 "Wether Positions should be broadcasted across all ranks");
      26             :   // This parameter should be set by each derived class depending on whether positions should
      27             :   // be sorted after being communicated. Sorting is usually undesirable, as the user input order
      28             :   // should be respected. It may be only desirable for testing reproducibility
      29      188245 :   params.addRequiredParam<bool>("auto_sort", "Whether Positions should be sorted by coordinates");
      30             : 
      31             :   // In general, no need to refresh Positions unless the mesh moved
      32      188245 :   params.set<ExecFlagEnum>("execute_on") = EXEC_NONE;
      33             : 
      34      188245 :   params.addParamNamesToGroup("auto_broadcast auto_sort", "Advanced");
      35      188245 :   params.registerBase("Positions");
      36      188245 :   return params;
      37           0 : }
      38             : 
      39        1400 : Positions::Positions(const InputParameters & parameters)
      40             :   : GeneralReporter(parameters),
      41             :     // leverage reporter interface to keep track of consumers
      42        2800 :     _initial_positions(isParamValid("initial_positions")
      43        1416 :                            ? &getReporterValueByName<std::vector<Point>>(
      44        1416 :                                  getParam<PositionsName>("initial_positions") + "/positions_1d")
      45             :                            : nullptr),
      46             :     // Positions will be replicated on every rank so transfers may query positions from all ranks
      47        1400 :     _positions(declareValueByName<std::vector<Point>, ReporterVectorContext<Point>>(
      48             :         "positions_1d", REPORTER_MODE_REPLICATED)),
      49        1400 :     _need_broadcast(getParam<bool>("auto_broadcast")),
      50        1400 :     _need_sort(getParam<bool>("auto_sort")),
      51        2800 :     _initialized(false)
      52             : {
      53             : 
      54        1400 :   if (_initial_positions)
      55          16 :     _initial_positions_kd_tree = std::make_unique<KDTree>(*_initial_positions, 1);
      56        1400 : }
      57             : 
      58             : const Point &
      59      154319 : Positions::getPosition(unsigned int index, bool initial) const
      60             : {
      61             :   mooseAssert(initialized(initial), "Positions vector has not been initialized.");
      62             :   // Check sizes of inital positions
      63      154319 :   if (initial && _initial_positions && (*_initial_positions).size() < index)
      64           0 :     mooseError("Initial positions is not sized or initialized appropriately");
      65             : 
      66             :   // Check sizes of regular positions
      67      154319 :   if (!initial)
      68             :   {
      69      116131 :     if (_initial_positions && _positions.size() != (*_initial_positions).size())
      70           4 :       mooseError("Initial positions and current positions array length do not match");
      71      116127 :     else if (_positions.size() < index)
      72           0 :       mooseError("Positions retrieved with an out-of-bound index: '",
      73             :                  index,
      74             :                  "' when there are only ",
      75           0 :                  _positions.size(),
      76             :                  " positions.");
      77             :   }
      78             : 
      79      154315 :   if (initial && _initial_positions)
      80           0 :     return (*_initial_positions)[index];
      81      154315 :   if (_positions.size())
      82      154315 :     return _positions[index];
      83             :   else
      84           0 :     mooseError("Positions vector has not been initialized.");
      85             : }
      86             : 
      87             : const Point &
      88      178815 : Positions::getNearestPosition(const Point & target, const bool initial) const
      89             : {
      90             :   mooseAssert(initialized(initial), "Positions vector has not been initialized.");
      91      178815 :   const auto & positions = (initial && _initial_positions) ? *_initial_positions : _positions;
      92      178815 :   return positions[getNearestPositionIndex(target, initial)];
      93             : }
      94             : 
      95             : unsigned int
      96      326824 : Positions::getNearestPositionIndex(const Point & target, const bool initial) const
      97             : {
      98             :   mooseAssert(initialized(initial), "Positions vector has not been initialized.");
      99      326824 :   std::vector<std::size_t> return_index(1);
     100             : 
     101      326824 :   if (initial && _initial_positions)
     102             :   {
     103             :     mooseAssert(_initial_positions_kd_tree, "Should have an initial positions KDTree");
     104           0 :     _initial_positions_kd_tree->neighborSearch(target, 1, return_index);
     105             :   }
     106             :   else
     107             :   {
     108             :     mooseAssert(_positions_kd_tree, "Should have a positions KDTree");
     109      326824 :     _positions_kd_tree->neighborSearch(target, 1, return_index);
     110             :   }
     111             : 
     112      653648 :   return return_index[0];
     113      326824 : }
     114             : 
     115             : const std::vector<Point> &
     116         798 : Positions::getPositions(bool initial) const
     117             : {
     118         798 :   if (initial && _initial_positions)
     119          32 :     return *_initial_positions;
     120         766 :   if (_positions.size())
     121         761 :     return _positions;
     122             :   else
     123           5 :     mooseError("Positions vector has not been initialized.");
     124             : }
     125             : 
     126             : const std::vector<std::vector<Point>> &
     127           1 : Positions::getPositionsVector2D() const
     128             : {
     129           1 :   if (_positions_2d.size())
     130           0 :     return _positions_2d;
     131             :   else
     132           1 :     mooseError("2D positions vectors have not been initialized.");
     133             : }
     134             : 
     135             : const std::vector<std::vector<std::vector<Point>>> &
     136           1 : Positions::getPositionsVector3D() const
     137             : {
     138           1 :   if (_positions_3d.size())
     139           0 :     return _positions_3d;
     140             :   else
     141           1 :     mooseError("3D positions vectors have not been initialized.");
     142             : }
     143             : 
     144             : const std::vector<std::vector<std::vector<std::vector<Point>>>> &
     145           1 : Positions::getPositionsVector4D() const
     146             : {
     147           1 :   if (_positions_4d.size())
     148           0 :     return _positions_4d;
     149             :   else
     150           1 :     mooseError("4D positions vectors have not been initialized.");
     151             : }
     152             : 
     153             : void
     154        1278 : Positions::clearPositions()
     155             : {
     156        1278 :   _positions.clear();
     157        1278 :   _positions_2d.clear();
     158        1278 :   _positions_3d.clear();
     159        1278 :   _positions_4d.clear();
     160        1278 : }
     161             : 
     162             : void
     163          72 : Positions::unrollMultiDPositions()
     164             : {
     165             :   // Unroll in every dimension available
     166          72 :   std::vector<Point> temp_2d_unrolled;
     167          72 :   std::vector<Point> temp_3d_unrolled;
     168          72 :   std::vector<Point> temp_4d_unrolled;
     169          96 :   for (auto vec : _positions_2d)
     170          24 :     temp_2d_unrolled.insert(temp_2d_unrolled.end(), vec.begin(), vec.end());
     171         108 :   for (const auto & vec_vec : _positions_3d)
     172         108 :     for (const auto & vec : vec_vec)
     173          72 :       temp_3d_unrolled.insert(temp_3d_unrolled.end(), vec.begin(), vec.end());
     174          96 :   for (const auto & vec_vec_vec : _positions_4d)
     175          48 :     for (const auto & vec_vec : vec_vec_vec)
     176         120 :       for (const auto & vec : vec_vec)
     177          96 :         temp_4d_unrolled.insert(temp_4d_unrolled.end(), vec.begin(), vec.end());
     178             : 
     179             :   // for now we wont even tolerate a different ordering
     180             :   // 2D & 1D: check for conflicts
     181          72 :   if (temp_2d_unrolled.size() && _positions.size() && temp_2d_unrolled != _positions)
     182           0 :     mooseError("Inconsistency between the 2D and 1D position vectors detected. "
     183             :                "The 2D positions must unroll into the 1D positions");
     184             : 
     185             :   // 3D vs 2D & 1D
     186          72 :   if (temp_3d_unrolled.size() && temp_2d_unrolled.size() && temp_3d_unrolled != temp_3d_unrolled)
     187           0 :     mooseError("Inconsistency between the 3D and 2D position vectors detected. "
     188             :                "The 3D positions must unroll the same way as the 2D positions");
     189          72 :   if (temp_3d_unrolled.size() && _positions.size() && temp_3d_unrolled != _positions)
     190           0 :     mooseError("Inconsistency between the 3D and 1D position vectors detected. "
     191             :                "The 3D positions must unroll into the 1D positions");
     192             : 
     193             :   // 4D vs all lower dimensions
     194          72 :   if (temp_4d_unrolled.size() && temp_3d_unrolled.size() && temp_4d_unrolled != temp_3d_unrolled)
     195           0 :     mooseError("Inconsistency between the 4D and 3D position vectors detected. "
     196             :                "The 4D positions must unroll the same way as the 3D positions");
     197          72 :   if (temp_4d_unrolled.size() && temp_2d_unrolled.size() && temp_4d_unrolled != temp_2d_unrolled)
     198           0 :     mooseError("Inconsistency between the 4D and 2D position vectors detected. "
     199             :                "The 4D positions must unroll the same way as the 2D positions");
     200          72 :   if (temp_4d_unrolled.size() && _positions.size() && temp_4d_unrolled != _positions)
     201           0 :     mooseError("Inconsistency between the 4D and 1D position vectors detected. "
     202             :                "The 4D positions must unroll into the 1D positions");
     203             : 
     204             :   // Use any of the higher D to set the one D vector
     205          72 :   if (!_positions.size())
     206             :   {
     207          48 :     if (temp_2d_unrolled.size())
     208          12 :       _positions = temp_2d_unrolled;
     209          36 :     else if (temp_3d_unrolled.size())
     210          24 :       _positions = temp_3d_unrolled;
     211          12 :     else if (temp_4d_unrolled.size())
     212          12 :       _positions = temp_4d_unrolled;
     213             :     else
     214           0 :       mooseError("Positions::unrollMultiDPositions() may only be called if there is at least one "
     215             :                  "non-empty positions vector.");
     216             :   }
     217          72 : }
     218             : 
     219             : void
     220        1641 : Positions::finalize()
     221             : {
     222             :   // Gather up the positions vector on all ranks
     223             :   mooseAssert(initialized(false), "Positions vector has not been initialized.");
     224        1641 :   if (_need_broadcast)
     225             :     // The consumer/producer reporter interface can keep track of whether a reduction is needed
     226             :     // (for example if a consumer needs replicated data, but the producer is distributed) however,
     227             :     // we have currently made the decision that positions should ALWAYS be replicated
     228         183 :     _communicator.allgather(_positions, /* identical buffer lengths = */ false);
     229             : 
     230             :   // Sort positions by X then Y then Z
     231        1641 :   if (_need_sort)
     232         116 :     std::sort(_positions.begin(), _positions.end());
     233             : 
     234             :   // Make a KDTree with the positions
     235        1641 :   _positions_kd_tree = std::make_unique<KDTree>(_positions, 1);
     236             : 
     237             :   // For now at least, we expect positions to be the same on all ranks
     238             :   mooseAssert(comm().verify(_positions), "Positions should be the same across all MPI processes.");
     239        1641 : }
     240             : 
     241             : Real
     242          51 : Positions::getMinDistanceBetweenPositions() const
     243             : {
     244             :   mooseAssert(initialized(false), "Positions vector has not been initialized.");
     245             :   // Dumb nested loops. We can revisit this once we have a KDTree for nearest position searching
     246             :   // These nested loops are still faster than calling getNearestPositions on each position for now
     247          51 :   Real min_distance_sq = std::numeric_limits<Real>::max();
     248         114 :   for (const auto i1 : index_range(_positions))
     249          75 :     for (auto i2 = i1 + 1; i2 < _positions.size(); i2++)
     250          12 :       min_distance_sq = std::min(min_distance_sq, (_positions[i1] - _positions[i2]).norm_sq());
     251         102 :   return std::sqrt(min_distance_sq);
     252             : }
     253             : 
     254             : bool
     255           0 : Positions::arePositionsCoplanar() const
     256             : {
     257             :   mooseAssert(initialized(false), "Positions vector has not been initialized.");
     258           0 :   return MooseMeshUtils::isCoPlanar(_positions);
     259             : }
     260             : 
     261             : bool
     262         170 : Positions::initialized(bool initial) const
     263             : {
     264         170 :   if (initial && _initial_positions)
     265             :     // We do not forward the 'initial' status as we are not currently looking
     266             :     // to support two level nesting for initial positions
     267           0 :     return _fe_problem.getPositionsObject(getParam<PositionsName>("initial_positions"))
     268           0 :         .initialized(false);
     269             :   else
     270         170 :     return _initialized;
     271             : }