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 "CylindricalGridDivision.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", CylindricalGridDivision);
19 :
20 : InputParameters
21 14368 : CylindricalGridDivision::validParams()
22 : {
23 14368 : InputParameters params = MeshDivision::validParams();
24 14368 : params.addClassDescription(
25 : "Divide the mesh along a cylindrical grid. The innermost numbering of divisions is the "
26 : "radial bins, then comes the azimuthal bins, then the axial bins");
27 :
28 : // Definition of the cylinder(s)
29 14368 : params.addRequiredParam<Point>("axis_direction", "Direction of the cylinder's axis");
30 14368 : params.addRequiredParam<Point>(
31 : "azimuthal_start",
32 : "Direction of the 0-azimuthal-angle vector, normal to the cylinder's axis");
33 14368 : params.addParam<Point>("center",
34 : "Point on the cylinder's axis, acting as the center of this local "
35 : "R-theta-Z coordinate based division");
36 14368 : params.addParam<PositionsName>("center_positions",
37 : "Positions of the points on the cylinders' respective axis, "
38 : "acting as the center of the local "
39 : "R-theta-Z coordinate based divisions");
40 :
41 : // Spatial bounds of the cylinder(s)
42 14368 : params.addRangeCheckedParam<Real>(
43 : "r_min", 0, "r_min>=0", "Minimum radial coordinate (for a hollow cylinder)");
44 14368 : params.addRequiredRangeCheckedParam<Real>("r_max", "r_max>0", "Maximum radial coordinate");
45 43104 : params.addParam<Real>(
46 28736 : "cylinder_axial_min", std::numeric_limits<Real>::lowest(), "Minimum axial coordinate");
47 43104 : params.addParam<Real>(
48 28736 : "cylinder_axial_max", std::numeric_limits<Real>::max(), "Maximum axial coordinate");
49 :
50 : // Number of bins
51 14368 : params.addRequiredRangeCheckedParam<unsigned int>(
52 : "n_radial", "n_radial>0", "Number of divisions in the cylinder radial direction");
53 43104 : params.addRangeCheckedParam<unsigned int>(
54 28736 : "n_azimuthal", 1, "n_azimuthal>0", "Number of divisions in the azimuthal direction");
55 43104 : params.addRangeCheckedParam<unsigned int>(
56 28736 : "n_axial", 1, "n_axial>0", "Number of divisions in the cylinder axial direction");
57 :
58 43104 : params.addParam<bool>("assign_domain_outside_grid_to_border",
59 28736 : false,
60 : "Whether to map the domain outside the grid back to the border of the grid "
61 : "(radially or axially)");
62 :
63 14368 : return params;
64 0 : }
65 :
66 56 : CylindricalGridDivision::CylindricalGridDivision(const InputParameters & parameters)
67 : : MeshDivision(parameters),
68 56 : _direction(getParam<Point>("axis_direction")),
69 56 : _center(isParamValid("center") ? &getParam<Point>("center") : nullptr),
70 56 : _center_positions(
71 112 : isParamValid("center_positions")
72 56 : ? &_fe_problem->getPositionsObject(getParam<PositionsName>("center_positions"))
73 : : nullptr),
74 56 : _azim_dir(getParam<Point>("azimuthal_start")),
75 56 : _min_r(getParam<Real>("r_min")),
76 56 : _max_r(getParam<Real>("r_max")),
77 56 : _min_z(getParam<Real>("cylinder_axial_min")),
78 56 : _max_z(getParam<Real>("cylinder_axial_max")),
79 56 : _n_radial(getParam<unsigned int>("n_radial")),
80 56 : _n_azim(getParam<unsigned int>("n_azimuthal")),
81 56 : _n_axial(getParam<unsigned int>("n_axial")),
82 112 : _outside_grid_counts_as_border(getParam<bool>("assign_domain_outside_grid_to_border"))
83 : {
84 56 : CylindricalGridDivision::initialize();
85 :
86 : // Check that we know the centers
87 52 : if (!_center && !_center_positions)
88 0 : paramError("center", "You must pass a parameter for the center of the cylindrical frame");
89 :
90 : // Check axis
91 52 : if (!MooseUtils::absoluteFuzzyEqual(_direction.norm_sq(), 1))
92 0 : paramError("axis_direction", "Axis must have a norm of 1");
93 52 : if (!MooseUtils::absoluteFuzzyEqual(_azim_dir.norm_sq(), 1))
94 0 : paramError("azimuthal_start", "Azimuthal axis must have a norm of 1");
95 :
96 : // Check non-negative size
97 52 : if (_max_r < _min_r)
98 0 : paramError("r_min", "Maximum radius must be larger than minimum radius");
99 52 : if (_max_z < _min_z)
100 0 : paramError("cylinder_axial_min", "Maximum axial extent must be larger than minimum");
101 :
102 : // Check non-zero size if subdivided
103 52 : if (_n_radial > 1 && MooseUtils::absoluteFuzzyEqual(_min_r, _max_r))
104 0 : paramError("n_radial", "Zero-thickness cylinder cannot be subdivided radially");
105 52 : if (_n_axial > 1 && MooseUtils::absoluteFuzzyEqual(_min_z, _max_z))
106 0 : paramError("n_axial", "Zero-height cylinder cannot be subdivided axially");
107 :
108 : // Check non-infinite size if subdivided
109 78 : if ((_min_z == std::numeric_limits<Real>::lowest() ||
110 104 : _max_z == std::numeric_limits<Real>::max()) &&
111 26 : _n_axial > 1)
112 0 : paramError("n_axial",
113 : "Infinite-height cylinder cannot be subdivided axially. Please specify both a "
114 : "cylinder axial minimum and maximum extent");
115 52 : }
116 :
117 : void
118 56 : CylindricalGridDivision::initialize()
119 : {
120 56 : if (!_center_positions)
121 39 : setNumDivisions(_n_radial * _n_azim * _n_axial);
122 : else
123 17 : setNumDivisions(_center_positions->getNumPositions() * _n_radial * _n_azim * _n_axial);
124 :
125 : // Check that the grid is well-defined
126 56 : if (_center_positions)
127 : {
128 17 : Real min_dist = 2 * _max_r;
129 17 : Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
130 : // Note that if the positions are not co-planar, the distance reported would be bigger but there
131 : // could still be an overlap. Looking at min_center_dist is not enough
132 17 : if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
133 4 : mooseWarning(
134 : "Cylindrical grids centered on the positions are too close to each other (min distance: ",
135 : min_center_dist,
136 : "), closer than the radial extent of each grid. Mesh division is ill-defined");
137 : }
138 :
139 : // We could alternatively check every point in the mesh but it seems expensive
140 : // A bounding box check could also suffice but the cylinders would need to be excessively large to
141 : // enclose the entire mesh
142 52 : _mesh_fully_indexed = true;
143 52 : if (!_outside_grid_counts_as_border || _center_positions)
144 52 : _mesh_fully_indexed = false;
145 52 : }
146 :
147 : unsigned int
148 14083 : CylindricalGridDivision::divisionIndex(const Elem & elem) const
149 : {
150 14083 : return divisionIndex(elem.vertex_average());
151 : }
152 :
153 : unsigned int
154 14109 : CylindricalGridDivision::divisionIndex(const Point & pt) const
155 : {
156 : // Compute coordinates of the point in the cylindrical coordinates
157 14109 : Point pc;
158 14109 : Point in_plane;
159 14109 : unsigned int offset = 0;
160 14109 : if (_center)
161 : {
162 13469 : pc(2) = (pt - *_center) * _direction;
163 13469 : in_plane = (pt - *_center) - pc(2) * _direction;
164 : }
165 : else
166 : {
167 : // If dividing using positions, find the closest position
168 640 : const bool initial = _fe_problem->getCurrentExecuteOnFlag() == EXEC_INITIAL;
169 640 : const auto nearest_center_index = _center_positions->getNearestPositionIndex(pt, initial);
170 640 : offset = nearest_center_index * _n_radial * _n_azim * _n_axial;
171 640 : const auto new_center = _center_positions->getPosition(nearest_center_index, initial);
172 640 : pc(2) = (pt - new_center) * _direction;
173 640 : in_plane = (pt - new_center) - pc(2) * _direction;
174 : }
175 :
176 14109 : pc(0) = in_plane.norm();
177 14109 : const Point pi2_dir = _azim_dir.cross(_direction);
178 14109 : pc(1) = std::atan2(in_plane * pi2_dir, in_plane * _azim_dir) + libMesh::pi;
179 :
180 14109 : if (!_outside_grid_counts_as_border)
181 : {
182 28087 : if (MooseUtils::absoluteFuzzyLessThan(pc(0), _min_r) ||
183 13978 : MooseUtils::absoluteFuzzyGreaterThan(pc(0), _max_r))
184 135 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
185 27660 : if (MooseUtils::absoluteFuzzyLessThan(pc(2), _min_z) ||
186 13686 : MooseUtils::absoluteFuzzyGreaterThan(pc(2), _max_z))
187 288 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
188 : }
189 :
190 13686 : const auto not_found = MooseMeshDivision::INVALID_DIVISION_INDEX;
191 13686 : auto ir = not_found, ia = not_found, iz = not_found;
192 13686 : const Point widths(_max_r - _min_r, 2 * libMesh::pi, _max_z - _min_z);
193 :
194 : // Look inside the grid and on the left / back / bottom
195 30674 : for (const auto jr : make_range(_n_radial + 1))
196 : {
197 30674 : const auto border_r = _min_r + widths(0) * jr / _n_radial;
198 30674 : if (jr > 0 && jr < _n_radial && MooseUtils::absoluteFuzzyEqual(border_r, pc(0)))
199 0 : mooseWarning(
200 0 : "Querying the division index for a point of a boundary between two regions radially: " +
201 0 : Moose::stringify(pt));
202 30674 : if (border_r >= pc(0))
203 : {
204 13686 : ir = jr > 0 ? jr - 1 : 0;
205 13686 : break;
206 : }
207 : }
208 47779 : for (const auto ja : make_range(_n_azim + 1))
209 : {
210 47779 : const auto border_a = widths(1) * ja / _n_azim;
211 47779 : if (ja > 0 && ja < _n_azim && MooseUtils::absoluteFuzzyEqual(border_a, pc(1)))
212 0 : mooseWarning("Querying the division index for a point of a boundary between two regions "
213 0 : "azimuthally : " +
214 0 : Moose::stringify(pt));
215 47779 : if (border_a >= pc(1))
216 : {
217 13686 : ia = ja > 0 ? ja - 1 : 0;
218 13686 : break;
219 : }
220 : }
221 30252 : for (const auto jz : make_range(_n_axial + 1))
222 : {
223 30252 : const auto border_z = _min_z + widths(2) * jz / _n_axial;
224 30252 : if (jz > 0 && jz < _n_axial && MooseUtils::absoluteFuzzyEqual(border_z, pc(2)))
225 0 : mooseWarning("Querying the division index for a point of a boundary between two axial "
226 0 : "regions along the cylinder axis: " +
227 0 : Moose::stringify(pt));
228 30252 : if (border_z >= pc(2))
229 : {
230 13686 : iz = jz > 0 ? jz - 1 : 0;
231 13686 : break;
232 : }
233 : }
234 :
235 : // Look beyond the edges of the grid
236 13686 : if (MooseUtils::absoluteFuzzyGreaterEqual(pc(0), _max_r))
237 0 : ir = _n_radial - 1;
238 13686 : if (MooseUtils::absoluteFuzzyGreaterEqual(pc(2), _max_z))
239 0 : iz = _n_axial - 1;
240 :
241 : // Handle edge case on widths
242 13686 : if (ir == not_found && MooseUtils::absoluteFuzzyEqual(widths(0), 0))
243 0 : ir = 0;
244 13686 : if (iz == not_found && MooseUtils::absoluteFuzzyEqual(widths(2), 0))
245 0 : iz = 0;
246 : mooseAssert(ir != not_found, "We should have found a mesh division bin radially");
247 : mooseAssert(ia != not_found, "We should have found a mesh division bin azimuthally");
248 : mooseAssert(iz != not_found, "We should have found a mesh division bin axially");
249 :
250 13686 : return offset + ir + _n_radial * ia + iz * _n_radial * _n_azim;
251 : }
|
