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 "SphericalGridDivision.h"
      11             : #include "MooseMesh.h"
      12             : #include "FEProblemBase.h"
      13             : #include "Positions.h"
      14             : 
      15             : #include "libmesh/elem.h"
      16             : #include <math.h>
      17             : 
      18             : registerMooseObject("MooseApp", SphericalGridDivision);
      19             : 
      20             : InputParameters
      21       14319 : SphericalGridDivision::validParams()
      22             : {
      23       14319 :   InputParameters params = MeshDivision::validParams();
      24       14319 :   params.addClassDescription("Divide the mesh along a spherical grid.");
      25             : 
      26             :   // Definition of the sphere(s)
      27       14319 :   params.addParam<Point>("center", "Center of the sphere");
      28       14319 :   params.addParam<PositionsName>("center_positions", "Positions of the centers of the spheres");
      29             : 
      30             :   // Spatial bounds of the sphere
      31       14319 :   params.addRangeCheckedParam<Real>(
      32             :       "r_min", 0, "r_min>=0", "Minimum radial coordinate (for a hollow sphere)");
      33       14319 :   params.addRequiredRangeCheckedParam<Real>("r_max", "r_max>0", "Maximum radial coordinate");
      34             : 
      35             :   // Number of bins
      36       14319 :   params.addRequiredRangeCheckedParam<unsigned int>(
      37             :       "n_radial", "n_radial>0", "Number of divisions in the sphere radial direction");
      38             : 
      39       42957 :   params.addParam<bool>("assign_domain_outside_grid_to_border",
      40       28638 :                         false,
      41             :                         "Whether to map the domain outside the grid back to the border of the grid "
      42             :                         "(radially or axially)");
      43             : 
      44       14319 :   return params;
      45           0 : }
      46             : 
      47          30 : SphericalGridDivision::SphericalGridDivision(const InputParameters & parameters)
      48             :   : MeshDivision(parameters),
      49          30 :     _center(isParamValid("center") ? &getParam<Point>("center") : nullptr),
      50          30 :     _center_positions(
      51          60 :         isParamValid("center_positions")
      52          30 :             ? &_fe_problem->getPositionsObject(getParam<PositionsName>("center_positions"))
      53             :             : nullptr),
      54          30 :     _min_r(getParam<Real>("r_min")),
      55          30 :     _max_r(getParam<Real>("r_max")),
      56          30 :     _n_radial(getParam<unsigned int>("n_radial")),
      57          60 :     _outside_grid_counts_as_border(getParam<bool>("assign_domain_outside_grid_to_border"))
      58             : {
      59          30 :   SphericalGridDivision::initialize();
      60             : 
      61             :   // Check that we know the centers
      62          26 :   if (!_center && !_center_positions)
      63           0 :     paramError("center", "You must pass a parameter for the center of the spherical frame");
      64             : 
      65             :   // Check non-negative size
      66          26 :   if (_max_r < _min_r)
      67           0 :     paramError("r_min", "Maximum radius must be larger than minimum radius");
      68             : 
      69             :   // Check non-zero size if subdivided
      70          26 :   if (_n_radial > 1 && MooseUtils::absoluteFuzzyEqual(_min_r, _max_r))
      71           0 :     paramError("n_radial", "Zero-thickness sphere cannot be subdivided radially");
      72          26 : }
      73             : 
      74             : void
      75          30 : SphericalGridDivision::initialize()
      76             : {
      77          30 :   if (!_center_positions)
      78          13 :     setNumDivisions(_n_radial);
      79             :   else
      80          17 :     setNumDivisions(_center_positions->getNumPositions() * _n_radial);
      81             : 
      82             :   // Check that the grid is well-defined
      83          30 :   if (_center_positions)
      84             :   {
      85          17 :     Real min_dist = 2 * _max_r;
      86          17 :     Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
      87             :     // Note that if the positions are not co-planar, the distance reported would be bigger but there
      88             :     // could still be an overlap. Looking at min_center_dist is not enough
      89          17 :     if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
      90           4 :       mooseWarning(
      91             :           "Spherical grids centered on the positions are too close to each other (min distance: ",
      92             :           min_center_dist,
      93             :           "), closer than the radial extent of each grid. Mesh division is ill-defined");
      94             :   }
      95             : 
      96             :   // We could alternatively check every point in the mesh but it seems expensive
      97             :   // A bounding box check could also suffice but the sphere would need to be excessively large to
      98             :   // enclose the entire mesh
      99          26 :   _mesh_fully_indexed = true;
     100          26 :   if (!_outside_grid_counts_as_border || _center_positions)
     101          26 :     _mesh_fully_indexed = false;
     102          26 : }
     103             : 
     104             : unsigned int
     105        1280 : SphericalGridDivision::divisionIndex(const Elem & elem) const
     106             : {
     107        1280 :   return divisionIndex(elem.vertex_average());
     108             : }
     109             : 
     110             : unsigned int
     111        1280 : SphericalGridDivision::divisionIndex(const Point & pt) const
     112             : {
     113             :   // Compute coordinates of the point in the spherical coordinates
     114        1280 :   Point pc;
     115        1280 :   unsigned int offset = 0;
     116        1280 :   if (_center)
     117         640 :     pc(0) = (pt - *_center).norm();
     118             :   else
     119             :   {
     120             :     // If distributing using positions, find the closest position
     121         640 :     const bool initial = _fe_problem->getCurrentExecuteOnFlag() == EXEC_INITIAL;
     122         640 :     const auto nearest_center_index = _center_positions->getNearestPositionIndex(pt, initial);
     123         640 :     offset = nearest_center_index * _n_radial;
     124         640 :     pc(0) = (pt - _center_positions->getPosition(nearest_center_index, initial)).norm();
     125             :   }
     126             : 
     127        1280 :   if (!_outside_grid_counts_as_border)
     128             :   {
     129        2544 :     if (MooseUtils::absoluteFuzzyLessThan(pc(0), _min_r) ||
     130        1264 :         MooseUtils::absoluteFuzzyGreaterThan(pc(0), _max_r))
     131          16 :       return MooseMeshDivision::INVALID_DIVISION_INDEX;
     132             :   }
     133             : 
     134        1264 :   const auto not_found = MooseMeshDivision::INVALID_DIVISION_INDEX;
     135        1264 :   auto ir = not_found;
     136        1264 :   const Real width = _max_r - _min_r;
     137             : 
     138             :   // Look inside the grid and on the left / back / bottom
     139        3664 :   for (const auto jr : make_range(_n_radial + 1))
     140             :   {
     141        3664 :     const auto border_r = _min_r + width * jr / _n_radial;
     142        3664 :     if (jr > 0 && jr < _n_radial && MooseUtils::absoluteFuzzyEqual(border_r, pc(0)))
     143           0 :       mooseWarning(
     144           0 :           "Querying the division index for a point of a boundary between two regions radially: " +
     145           0 :           Moose::stringify(pt));
     146        3664 :     if (border_r >= pc(0))
     147             :     {
     148        1264 :       ir = jr > 0 ? jr - 1 : 0;
     149        1264 :       break;
     150             :     }
     151             :   }
     152             : 
     153             :   // Look beyond the extent of the radial grid
     154        1264 :   if (MooseUtils::absoluteFuzzyGreaterEqual(pc(0), _max_r))
     155           0 :     ir = _n_radial - 1;
     156             : 
     157             :   // Handle edge case on widths
     158        1264 :   if (ir == not_found && MooseUtils::absoluteFuzzyEqual(width, 0))
     159           0 :     ir = 0;
     160             : 
     161             :   mooseAssert(ir != not_found, "We should have found a mesh division bin radially");
     162        1264 :   return offset + ir;
     163             : }