LCOV - code coverage report
Current view: top level - src/meshdivisions - HexagonalGridDivision.C (source / functions) Hit Total Coverage
Test: idaholab/moose reactor: #31405 (292dce) with base fef103 Lines: 75 86 87.2 %
Date: 2025-09-04 07:56:24 Functions: 5 5 100.0 %
Legend: Lines: hit not hit

          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 "HexagonalGridDivision.h"
      11             : #include "MooseMesh.h"
      12             : #include "HexagonalLatticeUtils.h"
      13             : #include "Positions.h"
      14             : 
      15             : #include "libmesh/elem.h"
      16             : 
      17             : registerMooseObject("ReactorApp", HexagonalGridDivision);
      18             : 
      19             : InputParameters
      20          62 : HexagonalGridDivision::validParams()
      21             : {
      22          62 :   InputParameters params = MeshDivision::validParams();
      23          62 :   params.addClassDescription(
      24             :       "Divide the mesh along a hexagonal grid. Numbering of pin divisions increases first "
      25             :       "counterclockwise, then expanding outwards from the inner ring, then axially. "
      26             :       "Inner-numbering is within a radial ring, outer-numbering is axial divisions");
      27             : 
      28         124 :   params.addParam<Point>("center", "Center of the hexagonal grid");
      29         124 :   params.addParam<PositionsName>("center_positions", "Centers of the hexagonal grids");
      30             : 
      31         124 :   params.addRequiredRangeCheckedParam<Real>(
      32             :       "lattice_flat_to_flat",
      33             :       "lattice_flat_to_flat>0",
      34             :       "Distance between two (inner) opposite sides of a lattice. Also known as bundle pitch or "
      35             :       "inner flat-to-flat distance");
      36         124 :   params.addRequiredRangeCheckedParam<Real>("pin_pitch", "pin_pitch>0", "Distance between pins");
      37             : 
      38         124 :   params.addRequiredParam<Real>("z_min", "Minimal axial extent of the lattice");
      39         124 :   params.addRequiredParam<Real>("z_max", "Maximum axial extent of the lattice");
      40         124 :   params.addRequiredRangeCheckedParam<unsigned int>("nr", "nr>0", "Number of hexagonal rings");
      41         124 :   params.addRequiredRangeCheckedParam<unsigned int>("nz", "nz>0", "Number of divisions in Z");
      42         124 :   params.addParam<bool>(
      43             :       "assign_domain_outside_grid_to_border",
      44         124 :       false,
      45             :       "Whether to map the domain outside the grid back to the border of the grid");
      46         124 :   params.addParam<Real>("rotation_around_axis",
      47         124 :                         0.,
      48             :                         "Rotation angle to apply to the underlying hexagonal lattice (in degrees)");
      49             : 
      50          62 :   return params;
      51           0 : }
      52             : 
      53          34 : HexagonalGridDivision::HexagonalGridDivision(const InputParameters & parameters)
      54             :   : MeshDivision(parameters),
      55          90 :     _center(isParamValid("center") ? getParam<Point>("center") : Point(0, 0, 0)),
      56          34 :     _center_positions(
      57          34 :         isParamValid("center_positions")
      58          58 :             ? &_fe_problem->getPositionsObject(getParam<PositionsName>("center_positions"))
      59             :             : nullptr),
      60          68 :     _lattice_flat_to_flat(getParam<Real>("lattice_flat_to_flat")),
      61          68 :     _pin_pitch(getParam<Real>("pin_pitch")),
      62          68 :     _z_axis_index(MooseEnum("X Y Z", "Z")),
      63          68 :     _min_z(getParam<Real>("z_min")),
      64          68 :     _max_z(getParam<Real>("z_max")),
      65          68 :     _nr(getParam<unsigned int>("nr")),
      66          68 :     _nz(getParam<unsigned int>("nz")),
      67         102 :     _outside_grid_counts_as_border(getParam<bool>("assign_domain_outside_grid_to_border"))
      68             : {
      69          34 :   HexagonalGridDivision::initialize();
      70             : 
      71         102 :   if (!isParamValid("center") && !_center_positions)
      72           0 :     paramError("center", "A center must be provided, or a Positions object for the centers");
      73          34 :   if (_pin_pitch > _lattice_flat_to_flat)
      74           0 :     mooseError("lattice_flat_to_flat", "Pin pitch should be smaller than bundle pitch");
      75          34 :   if ((_nz > 1) && MooseUtils::absoluteFuzzyEqual(_max_z, _min_z))
      76           0 :     paramError("nz", "Subdivision number must be 1 if width is 0 in Z direction");
      77          34 : }
      78             : 
      79             : void
      80          34 : HexagonalGridDivision::initialize()
      81             : {
      82             :   // We make very large pins so they cover the entire position
      83          68 :   _hex_latt = std::make_unique<HexagonalLatticeUtils>(_lattice_flat_to_flat,
      84             :                                                       _pin_pitch,
      85          34 :                                                       _pin_pitch,
      86          68 :                                                       0.,
      87          34 :                                                       1.,
      88          34 :                                                       _nr,
      89          34 :                                                       _z_axis_index,
      90             :                                                       getParam<Real>("rotation_around_axis"));
      91             : 
      92          34 :   if (!_center_positions)
      93          22 :     setNumDivisions(_hex_latt->totalPins(_nr) * _nz);
      94             :   else
      95          12 :     setNumDivisions(_center_positions->getNumPositions() * _hex_latt->totalPins(_nr) * _nz);
      96             : 
      97             :   // Check that the grid is well-defined
      98          34 :   if (_center_positions)
      99             :   {
     100          12 :     const Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
     101             :     // Note that if the positions are not aligned on a hexagonal lattice themselves,
     102             :     // this bound is not sufficiently strict. The simplest example would be non-coplanar
     103             :     // points, which can be a great distance away axially but be on the same axis
     104          12 :     if (MooseUtils::absoluteFuzzyGreaterThan(_lattice_flat_to_flat, min_center_dist))
     105           0 :       mooseWarning(
     106             :           "Hexagonal grids centered on the positions are too close to each other (min distance: ",
     107             :           min_center_dist,
     108             :           "), closer than the extent of each grid (",
     109             :           _lattice_flat_to_flat,
     110             :           "). Mesh division is ill-defined ");
     111             :   }
     112          34 : }
     113             : 
     114             : unsigned int
     115       28008 : HexagonalGridDivision::divisionIndex(const Elem & elem) const
     116             : {
     117       28008 :   return divisionIndex(elem.vertex_average());
     118             : }
     119             : 
     120             : unsigned int
     121       28008 : HexagonalGridDivision::divisionIndex(const Point & pt) const
     122             : {
     123             :   unsigned int offset = 0;
     124             : 
     125             :   // Get point in the coordinates of the lattice. This can involve a projection due to
     126             :   // the axis of the lattice, or simply a translation if there are lattices distributed
     127             :   // using positions
     128             :   Point pc;
     129       28008 :   if (_center_positions)
     130             :   {
     131             :     // If dividing using positions, find the closest position and
     132             :     // look at the relative position of the point compared to that position
     133        8568 :     const bool initial = _fe_problem->getCurrentExecuteOnFlag() == EXEC_INITIAL;
     134        8568 :     const auto nearest_grid_center_index = _center_positions->getNearestPositionIndex(pt, initial);
     135        8568 :     offset = nearest_grid_center_index * _hex_latt->totalPins(_nr) * _nz;
     136             :     const auto nearest_grid_center =
     137        8568 :         _center_positions->getPosition(nearest_grid_center_index, initial);
     138             : 
     139             :     // Project in local hexagonal grid
     140        8568 :     pc = pt - nearest_grid_center;
     141             :   }
     142             :   else
     143       19440 :     pc = pt - _center;
     144             : 
     145             :   // Get radial division index, using the channel as the pins are 0-radius
     146             :   // The logic in get pin index requires getting the point in the plane of the pin centers
     147       28008 :   auto ir = _hex_latt->pinIndex(pc);
     148             :   const auto n_pins = _hex_latt->nPins();
     149             : 
     150       28008 :   if (!_outside_grid_counts_as_border)
     151             :   {
     152       28008 :     if (ir == n_pins)
     153        8532 :       return MooseMeshDivision::INVALID_DIVISION_INDEX;
     154       19476 :     if (MooseUtils::absoluteFuzzyLessThan(pc(_z_axis_index), _min_z) ||
     155             :         MooseUtils::absoluteFuzzyGreaterThan(pc(_z_axis_index), _max_z))
     156           0 :       return MooseMeshDivision::INVALID_DIVISION_INDEX;
     157             :   }
     158             : 
     159             :   // If too far from the grid to have a valid radial index, use the closest pin
     160       19476 :   if (ir == n_pins)
     161           0 :     ir = _hex_latt->closestPinIndex(pc);
     162             : 
     163             :   // Find axial index
     164       19476 :   const auto not_found = MooseMeshDivision::INVALID_DIVISION_INDEX;
     165             :   auto iz = not_found;
     166       35001 :   for (const auto jz : make_range(_nz + 1))
     167             :   {
     168       35001 :     const auto border_z = _min_z + (_max_z - _min_z) * jz / _nz;
     169       35001 :     if (jz > 0 && jz < _nz && MooseUtils::absoluteFuzzyEqual(border_z, pc(_z_axis_index)))
     170           0 :       mooseWarning(
     171           0 :           "Querying the division index for a point of a boundary between two regions in Z: " +
     172           0 :               Moose::stringify(pt),
     173             :           ", in local hex grid frame: ",
     174           0 :           Moose::stringify(pc));
     175       35001 :     if (border_z >= pc(_z_axis_index))
     176             :     {
     177       19476 :       iz = (jz > 0) ? jz - 1 : 0;
     178             :       break;
     179             :     }
     180             :   }
     181             : 
     182             :   // Look on the top of the grid
     183       19476 :   if (MooseUtils::absoluteFuzzyGreaterEqual(pc(_z_axis_index), _max_z))
     184        7056 :     iz = _nz - 1;
     185             : 
     186             :   // Handle edge case on widths
     187       19476 :   if (iz == not_found && MooseUtils::absoluteFuzzyEqual(_max_z - _min_z, 0))
     188             :     iz = 0;
     189             :   mooseAssert(ir != not_found, "We should have found a mesh division bin radially");
     190             :   mooseAssert(iz != not_found, "We should have found a mesh division bin in Z");
     191             : 
     192       19476 :   const auto n_radial = _hex_latt->totalPins(_nr);
     193       19476 :   return offset + ir + iz * n_radial;
     194             : }

Generated by: LCOV version 1.14