Merge pull request #26311 from GiudGiud/PR_hex_div

Hexagonal division object

framework/doc/content/source/meshdivisions/CartesianGridDivision.md

@@ -8,6 +8,15 @@ Points that lie outside the Cartesian grid may be assigned to the grid outer bin
 using the [!param](/MeshDivisions/CartesianGridDivision/assign_domain_outside_grid_to_border)
 parameter.
 
+Using a [Positions](syntax/Positions/index.md) object as the [!param](/MeshDivisions/CartesianGridDivision/center_positions)
+parameter, multiple Cartesian grids can be created around each position computed by that object. The division index
+of a point is then:
+
+!equation
+\text{division index} = (i - 1) N_{\text{single division}} + \text{division index in Cartesian grid centered around position i}
+
+with $i$ the index in the `Positions` object of the position nearest from the point and $N_{\text{single division}}$ the number of divisions for a single Cartesian grid, based on the X/Y/Z discretization specified.
+
 !alert note
 For points lying within the standard tolerance of an internal boundary of the Cartesian grid, this object
 will output a warning. If you do not mind the indetermination on which bins they belong to but do mind

framework/doc/content/source/meshdivisions/CylindricalGridDivision.md

@@ -12,6 +12,15 @@ parameter.
 The [!param](/MeshDivisions/CylindricalGridDivision/center) should be a point on the axis of the cylinder
 that also corresponds to the reference for the minimal ([!param](/MeshDivisions/CylindricalGridDivision/cylinder_axial_min)) and maximal ([!param](/MeshDivisions/CylindricalGridDivision/cylinder_axial_max)) axial extent of the cylinder.
 
+Using a [Positions](syntax/Positions/index.md) object as the [!param](/MeshDivisions/CylindricalGridDivision/center_positions)
+parameter, multiple cylindrical grids can be created around each position computed by that object. The division index
+of a point is then:
+
+!equation
+\text{division index} = (i - 1) N_{\text{single division}} + \text{division index in cylindrical grid centered around position i}
+
+with $i$ the index in the `Positions` object of the position nearest from the point and $N_{\text{single division}}$ the number of divisions for a single cylindrical grid, based on the number of rings and axial discretization specified.
+
 !alert note
 We have not implemented restrictions in the azimuthal direction so the entire ($0$, $2 \pi$) arc will be split.
 This is a desirable extension of this object.

framework/doc/content/source/meshdivisions/SphericalGridDivision.md

@@ -9,6 +9,15 @@ Points that lie outside the spherical grid may be assigned to the closest grid o
 using the [!param](/MeshDivisions/SphericalGridDivision/assign_domain_outside_grid_to_border)
 parameter.
 
+Using a [Positions](syntax/Positions/index.md) object as the [!param](/MeshDivisions/CylindricalGridDivision/center_positions)
+parameter, multiple cylindrical grids can be created around each position computed by that object. The division index
+of a point is then:
+
+!equation
+\text{division index} = (i - 1) N_{\text{single division}} + \text{division index in cylindrical grid centered around position i}
+
+with $i$ the index in the `Positions` object of the position nearest from the point and $N_{\text{single division}}$ the number of divisions for a single spherical grid, based on the number of rings specified.
+
 !alert note
 We have not implemented binning in the azimuthal nor toroidal coordinates nor restrictions in those angular coordinates.
 This is a desirable extension of this object.

framework/include/utils/GeometryUtils.h

@@ -0,0 +1,155 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Moose.h"
+#include "libmesh/point.h"
+#include "libmesh/mesh_tools.h"
+
+namespace geom_utils
+{
+/**
+ * Check whether a point is equal to zero
+ * @param[in] pt point
+ * @return whether point is equal to the zero point
+ */
+bool isPointZero(const Point & pt);
+
+/**
+ * Get the unit vector for a point parameter
+ * @param[in] pt point
+ * @param[in] name name of the parameter
+ */
+Point unitVector(const Point & pt, const std::string & name);
+
+/**
+ * Rotate point about an axis
+ * @param[in] p point
+ * @param[in] angle angle to rotate (radians)
+ * @param[in] axis axis expressed as vector
+ * @return rotated point
+ */
+Point rotatePointAboutAxis(const Point & p, const Real angle, const Point & axis);
+
+/**
+ * Get the minimum distance from a point to another set of points, in the plane
+ * perpendicular to the specified axis
+ * @param[in] pt point
+ * @param[in] candidates set of points we will find the nearest of
+ * @param[in] axis axis perpendicular to the plane of the polygon
+ * @return minimum distance to the provided points
+ */
+Real minDistanceToPoints(const Point & pt,
+                         const std::vector<Point> & candidates,
+                         const unsigned int axis);
+
+/**
+ * Get the corner coordinates of a regular 2-D polygon, assuming a face of the polygon
+ * is parallel to the x0 axis
+ * @param[in] num_sides number of sides to polygon
+ * @param[in] radius distance from polygon center to a corner
+ * @param[in] axis axis perpendicular to the plane of the polygon
+ * @return corner coordinates
+ */
+std::vector<Point>
+polygonCorners(const unsigned int num_sides, const Real radius, const unsigned int axis);
+
+/**
+ * Get the indices of the plane perpendicular to the specified axis.
+ * For example, if the axis is the y-axis (1), then this will return
+ * (0, 2), indicating that the coordinates of a general 3-D point once
+ * projected onto the x-z plane can be obtained with the 0 and 2 indices.
+ * @param[in] axis axis perpendicular to the projection plane
+ * @return indices of coordinates on plane
+ */
+std::pair<unsigned int, unsigned int> projectedIndices(const unsigned int axis);
+
+/**
+ * Given two coordinates, construct a point in the 2-D plane perpendicular to the
+ * specified axis.
+ * @param[in] x0 first coordinate
+ * @param[in] x1 second coordinate
+ * @param[in] axis axis perpendicular to the projection plane
+ * @return point
+ */
+Point projectPoint(const Real x0, const Real x1, const unsigned int axis);
+
+/**
+ * Get the unit normal vector between two points (which are first projected onto
+ * the plane perpendicular to the 'axis'), such that the cross product of
+ * the unit normal with the line from pt1 to pt2 has a positive 'axis' component.
+ * @param[in] pt1 first point for line
+ * @param[in] pt2 second point for line
+ * @param[in] axis project points onto plane perpendicular to this axis
+ * @return unit normal
+ */
+Point projectedUnitNormal(Point pt1, Point pt2, const unsigned int axis);
+
+/**
+ * Compute the distance from a 3-D line, provided in terms of two points on the line
+ * @param[in] pt point of interest
+ * @param[in] line0 first point on line
+ * @param[in] line1 second point on line
+ * @return distance from line
+ */
+Real distanceFromLine(const Point & pt, const Point & line0, const Point & line1);
+
+/**
+ * Compute the distance from a 3-D line, provided in terms of two points on the line.
+ * Both the input point and the points on the line are projected into the 2-d plane
+ * perpendicular to the specified axis.
+ * @param[in] pt point of interest
+ * @param[in] line0 first point on line
+ * @param[in] line1 second point on line
+ * @param[in] axis axis index (0 = x, 1 = y, 2 = z) perpendicular to the projection plane
+ * @return distance from line
+ */
+Real projectedDistanceFromLine(Point pt, Point line0, Point line1, const unsigned int axis);
+
+/**
+ * If positive, point is on the positive side of the half space (and vice versa). Because
+ * a 3-D line does not have a negative or positive "side," you must provide the 'axis'
+ * perpendicular to the plane into which the point and line are first projected.
+ * @param[in] pt1 point of interest
+ * @param[in] pt2 one end point of line
+ * @param[in] pt3 other end point of line
+ * @param[in] axis axis perpendicular to plane onto which point and line are first projected
+ * @return half space of line
+ */
+Real projectedLineHalfSpace(Point pt1, Point pt2, Point pt3, const unsigned int axis);
+
+/**
+ * Whether a point is in 2-D a polygon in the plane perpendicular to the specified
+ * axis, given by corner points
+ * @param[in] point point of interest
+ * @param[in] corners corner points of polygon
+ * @param[in] axis axis perpendicular to the plane of the polygon
+ * @return whether point is inside the polygon
+ */
+bool
+pointInPolygon(const Point & point, const std::vector<Point> & corners, const unsigned int axis);
+
+/**
+ * Whether a point is on the edge of a 2-D polygon in the plane perpendicular to
+ * the specified axis, given its corner points
+ * @param[in] point point of interest
+ * @param[in] corners corner points of polygon
+ * @param[in] axis axis perpendicular to the plane of the polygon
+ * @return whether point is on edge of polygon
+ */
+bool pointOnEdge(const Point & point, const std::vector<Point> & corners, const unsigned int axis);
+
+/**
+ * Get corner points of a bounding box, with side length re-scaled
+ * @param[in] box bounding box to start from
+ * @param[in] factor by which to multiply the bounding box side
+ */
+std::vector<Point> boxCorners(const BoundingBox & box, const Real factor);
+} // end of namespace geom_utils

framework/src/meshdivisions/CartesianGridDivision.C

@@ -92,7 +92,10 @@ CartesianGridDivision::CartesianGridDivision(const InputParameters & parameters)
 void
 CartesianGridDivision::initialize()
 {
-  setNumDivisions(_nx * _ny * _nz);
+  if (!_center_positions)
+    setNumDivisions(_nx * _ny * _nz);
+  else
+    setNumDivisions(_center_positions->getNumPositions() * _nx * _ny * _nz);
 
   // Check that the grid is well-defined
   if (_center_positions)
@@ -101,7 +104,7 @@ CartesianGridDivision::initialize()
     Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
     // Note that if the positions are not co-planar, the distance reported would be bigger but there
     // could still be an overlap. Looking at min_center_dist is not enough
-    if (min_dist > min_center_dist)
+    if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
       mooseError(
           "Cartesian grids centered on the positions are too close to each other (min distance: ",
           min_center_dist,

framework/src/meshdivisions/CylindricalGridDivision.C

@@ -117,7 +117,10 @@ CylindricalGridDivision::CylindricalGridDivision(const InputParameters & paramet
 void
 CylindricalGridDivision::initialize()
 {
-  setNumDivisions(_n_radial * _n_azim * _n_axial);
+  if (!_center_positions)
+    setNumDivisions(_n_radial * _n_azim * _n_axial);
+  else
+    setNumDivisions(_center_positions->getNumPositions() * _n_radial * _n_azim * _n_axial);
 
   // Check that the grid is well-defined
   if (_center_positions)
@@ -126,7 +129,7 @@ CylindricalGridDivision::initialize()
     Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
     // Note that if the positions are not co-planar, the distance reported would be bigger but there
     // could still be an overlap. Looking at min_center_dist is not enough
-    if (min_dist > min_center_dist)
+    if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
       mooseError(
           "Cylindrical grids centered on the positions are too close to each other (min distance: ",
           min_center_dist,

framework/src/meshdivisions/SphericalGridDivision.C

@@ -74,7 +74,10 @@ SphericalGridDivision::SphericalGridDivision(const InputParameters & parameters)
 void
 SphericalGridDivision::initialize()
 {
-  setNumDivisions(_n_radial);
+  if (!_center_positions)
+    setNumDivisions(_n_radial);
+  else
+    setNumDivisions(_center_positions->getNumPositions() * _n_radial);
 
   // Check that the grid is well-defined
   if (_center_positions)
@@ -83,7 +86,7 @@ SphericalGridDivision::initialize()
     Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
     // Note that if the positions are not co-planar, the distance reported would be bigger but there
     // could still be an overlap. Looking at min_center_dist is not enough
-    if (min_dist > min_center_dist)
+    if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
       mooseError(
           "Spherical grids centered on the positions are too close to each other (min distance: ",
           min_center_dist,