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 "CartesianGridDivision.h"
11 : #include "MooseMesh.h"
12 : #include "Positions.h"
13 : #include "FEProblemBase.h"
14 :
15 : #include "libmesh/elem.h"
16 :
17 : registerMooseObject("MooseApp", CartesianGridDivision);
18 :
19 : InputParameters
20 4171 : CartesianGridDivision::validParams()
21 : {
22 4171 : InputParameters params = MeshDivision::validParams();
23 8342 : params.addClassDescription("Divide the mesh along a Cartesian grid. Numbering increases from "
24 : "bottom to top and from left to right and from back to front. The "
25 : "inner ordering is X, then Y, then Z");
26 16684 : params.addParam<Point>("bottom_left", "Bottom-back-left corner of the grid");
27 16684 : params.addParam<Point>("top_right", "Top-front-right corner of the grid");
28 16684 : params.addParam<Point>("center", "Center of the Cartesian grid");
29 16684 : params.addParam<PositionsName>("center_positions",
30 : "Positions of the centers of divided Cartesian grids");
31 16684 : params.addParam<Point>("widths", "Widths in the X, Y and Z directions");
32 25026 : params.addRequiredRangeCheckedParam<unsigned int>("nx", "nx>0", "Number of divisions in X");
33 25026 : params.addRequiredRangeCheckedParam<unsigned int>("ny", "ny>0", "Number of divisions in Y");
34 25026 : params.addRequiredRangeCheckedParam<unsigned int>("nz", "nz>0", "Number of divisions in Z");
35 8342 : params.addParam<bool>(
36 : "assign_domain_outside_grid_to_border",
37 8342 : false,
38 : "Whether to map the domain outside the grid back to the border of the grid");
39 :
40 4171 : return params;
41 0 : }
42 :
43 597 : CartesianGridDivision::CartesianGridDivision(const InputParameters & parameters)
44 : : MeshDivision(parameters),
45 1733 : _bottom_left(isParamValid("bottom_left") ? getParam<Point>("bottom_left") : Point(0, 0, 0)),
46 2330 : _top_right(isParamValid("top_right") ? getParam<Point>("top_right") : Point(0, 0, 0)),
47 1220 : _center(isParamValid("center") ? &getParam<Point>("center") : nullptr),
48 597 : _center_positions(
49 1194 : isParamValid("center_positions")
50 629 : ? &_fe_problem->getPositionsObject(getParam<PositionsName>("center_positions"))
51 : : nullptr),
52 1252 : _widths(isParamValid("widths") ? getParam<Point>("widths") : Point(_top_right - _bottom_left)),
53 1194 : _nx(getParam<unsigned int>("nx")),
54 1194 : _ny(getParam<unsigned int>("ny")),
55 1194 : _nz(getParam<unsigned int>("nz")),
56 1791 : _outside_grid_counts_as_border(getParam<bool>("assign_domain_outside_grid_to_border"))
57 : {
58 597 : CartesianGridDivision::initialize();
59 :
60 : // Check non-overlapping inputs for the dimensions of the grid
61 698 : if ((_center || _center_positions) && (isParamValid("bottom_left") || isParamValid("top_right")))
62 0 : mooseError("Either the center or the edges of the grids must be specified");
63 2970 : if ((isParamValid("top_right") + isParamValid("bottom_left") == 1) && !isParamValid("widths"))
64 0 : paramError("bottom_left",
65 : "Both bottom_left and top_right must be passed to be able to determine the width");
66 :
67 : // Pre-determine the bounds if we can
68 594 : if (!_center_positions && _center)
69 : {
70 13 : _bottom_left = *_center - _widths / 2;
71 13 : _top_right = *_center + _widths / 2;
72 : }
73 :
74 : // Check widths
75 594 : if (_widths(0) < 0)
76 0 : paramError("top_right",
77 : "Top-front-right corner must be right (X axis) of bottom-left-back corner");
78 594 : if (_widths(1) < 0)
79 0 : paramError("top_right",
80 : "Top-front-right corner must be in front (Y axis) of bottom-left-back corner");
81 594 : if (_widths(2) < 0)
82 0 : paramError("top_right",
83 : "Top-front-right corner must be on top (Z axis) of bottom-left-back corner");
84 594 : if ((_nx > 1) && MooseUtils::absoluteFuzzyEqual(_widths(0), 0))
85 0 : paramError("nx", "Subdivision number must be 1 if width is 0 in X direction");
86 594 : if ((_ny > 1) && MooseUtils::absoluteFuzzyEqual(_widths(1), 0))
87 0 : paramError("ny", "Subdivision number must be 1 if width is 0 in Y direction");
88 594 : if ((_nz > 1) && MooseUtils::absoluteFuzzyEqual(_widths(2), 0))
89 0 : paramError("nz", "Subdivision number must be 1 if width is 0 in Z direction");
90 594 : }
91 :
92 : void
93 597 : CartesianGridDivision::initialize()
94 : {
95 597 : if (!_center_positions)
96 581 : setNumDivisions(_nx * _ny * _nz);
97 : else
98 16 : setNumDivisions(_center_positions->getNumPositions() * _nx * _ny * _nz);
99 :
100 : // Check that the grid is well-defined
101 597 : if (_center_positions)
102 : {
103 16 : Real min_dist = _widths.norm();
104 16 : Real min_center_dist = _center_positions->getMinDistanceBetweenPositions();
105 : // Note that if the positions are not co-planar, the distance reported would be bigger but there
106 : // could still be an overlap. Looking at min_center_dist is not enough
107 16 : if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
108 3 : mooseWarning(
109 : "Cartesian grids centered on the positions are too close to each other (min distance: ",
110 : min_center_dist,
111 : "), closer than the extent of each grid. Mesh division is ill-defined");
112 : }
113 :
114 : // Examine mesh bounding box to see if the mesh may not be fully contained
115 594 : _mesh_fully_indexed = true;
116 594 : if (!_outside_grid_counts_as_border)
117 658 : for (auto i : make_range(LIBMESH_DIM))
118 : {
119 641 : if (_center_positions)
120 : {
121 13 : _mesh_fully_indexed = false;
122 13 : break;
123 : }
124 723 : else if (_bottom_left(i) > _mesh.getInflatedProcessorBoundingBox(0).first(i) ||
125 723 : _top_right(i) < _mesh.getInflatedProcessorBoundingBox(0).second(i))
126 : {
127 564 : _mesh_fully_indexed = false;
128 564 : break;
129 : }
130 : }
131 594 : }
132 :
133 : unsigned int
134 85283 : CartesianGridDivision::divisionIndex(const Elem & elem) const
135 : {
136 85283 : return divisionIndex(elem.vertex_average());
137 : }
138 :
139 : unsigned int
140 189432 : CartesianGridDivision::divisionIndex(const Point & pt) const
141 : {
142 189432 : unsigned int offset = 0;
143 : // Determine the local grid bounds
144 189432 : Point bottom_left, top_right, p;
145 189432 : if (_center_positions)
146 : {
147 : // If dividing using positions, find the closest position and
148 : // look at the relative position of the point compared to that position
149 640 : const bool initial = _fe_problem->getCurrentExecuteOnFlag() == EXEC_INITIAL;
150 640 : const auto nearest_grid_center_index = _center_positions->getNearestPositionIndex(pt, initial);
151 640 : offset = nearest_grid_center_index * _nx * _ny * _nz;
152 : const auto nearest_grid_center =
153 640 : _center_positions->getPosition(nearest_grid_center_index, initial);
154 640 : bottom_left = -_widths / 2;
155 640 : top_right = _widths / 2;
156 640 : p = pt - nearest_grid_center;
157 : }
158 : else
159 : {
160 188792 : bottom_left = _bottom_left;
161 188792 : top_right = _top_right;
162 188792 : p = pt;
163 : }
164 :
165 189432 : if (!_outside_grid_counts_as_border)
166 : {
167 321215 : if (MooseUtils::absoluteFuzzyLessThan(p(0), bottom_left(0)) ||
168 131783 : MooseUtils::absoluteFuzzyGreaterThan(p(0), top_right(0)))
169 72221 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
170 201272 : if (MooseUtils::absoluteFuzzyLessThan(p(1), bottom_left(1)) ||
171 84061 : MooseUtils::absoluteFuzzyGreaterThan(p(1), top_right(1)))
172 40979 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
173 152176 : if (MooseUtils::absoluteFuzzyLessThan(p(2), bottom_left(2)) ||
174 75944 : MooseUtils::absoluteFuzzyGreaterThan(p(2), top_right(2)))
175 864 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
176 : }
177 :
178 75368 : const auto not_found = MooseMeshDivision::INVALID_DIVISION_INDEX;
179 75368 : auto ix = not_found, iy = not_found, iz = not_found;
180 :
181 : // Look inside the grid and on the left / back / bottom
182 283120 : for (const auto jx : make_range(_nx + 1))
183 : {
184 283120 : const auto border_x = bottom_left(0) + _widths(0) * jx / _nx;
185 283120 : if (jx > 0 && jx < _nx && MooseUtils::absoluteFuzzyEqual(border_x, p(0)))
186 0 : mooseWarning(
187 0 : "Querying the division index for a point of a boundary between two regions in X: " +
188 0 : Moose::stringify(p));
189 283120 : if (border_x >= p(0))
190 : {
191 75368 : ix = (jx > 0) ? jx - 1 : 0;
192 75368 : break;
193 : }
194 : }
195 282774 : for (const auto jy : make_range(_ny + 1))
196 : {
197 282774 : const auto border_y = bottom_left(1) + _widths(1) * jy / _ny;
198 282774 : if (jy > 0 && jy < _ny && MooseUtils::absoluteFuzzyEqual(border_y, p(1)))
199 0 : mooseWarning(
200 0 : "Querying the division index for a point of a boundary between two regions in Y: " +
201 0 : Moose::stringify(p));
202 282774 : if (border_y >= p(1))
203 : {
204 75368 : iy = (jy > 0) ? jy - 1 : 0;
205 75368 : break;
206 : }
207 : }
208 75656 : for (const auto jz : make_range(_nz + 1))
209 : {
210 75656 : const auto border_z = bottom_left(2) + _widths(2) * jz / _nz;
211 75656 : if (jz > 0 && jz < _nz && MooseUtils::absoluteFuzzyEqual(border_z, p(2)))
212 0 : mooseWarning(
213 0 : "Querying the division index for a point of a boundary between two regions in Z: " +
214 0 : Moose::stringify(p));
215 75656 : if (border_z >= p(2))
216 : {
217 75368 : iz = (jz > 0) ? jz - 1 : 0;
218 75368 : break;
219 : }
220 : }
221 :
222 : // Look on the right / front / top of the grid
223 75368 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(0), top_right(0)))
224 0 : ix = _nx - 1;
225 75368 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(1), top_right(1)))
226 24 : iy = _ny - 1;
227 75368 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), top_right(2)))
228 75080 : iz = _nz - 1;
229 :
230 : // Handle edge case on widths
231 75368 : if (ix == not_found && MooseUtils::absoluteFuzzyEqual(_widths(0), 0))
232 0 : ix = 0;
233 75368 : if (iy == not_found && MooseUtils::absoluteFuzzyEqual(_widths(1), 0))
234 0 : iy = 0;
235 75368 : if (iz == not_found && MooseUtils::absoluteFuzzyEqual(_widths(2), 0))
236 0 : iz = 0;
237 : mooseAssert(ix != not_found, "We should have found a mesh division bin in X");
238 : mooseAssert(iy != not_found, "We should have found a mesh division bin in Y");
239 : mooseAssert(iz != not_found, "We should have found a mesh division bin in Z");
240 :
241 75368 : return offset + ix + _nx * iy + iz * _nx * _ny;
242 : }
|
