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 15428 : CartesianGridDivision::validParams()
21 : {
22 15428 : InputParameters params = MeshDivision::validParams();
23 15428 : 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 15428 : params.addParam<Point>("bottom_left", "Bottom-back-left corner of the grid");
27 15428 : params.addParam<Point>("top_right", "Top-front-right corner of the grid");
28 15428 : params.addParam<Point>("center", "Center of the Cartesian grid");
29 15428 : params.addParam<PositionsName>("center_positions",
30 : "Positions of the centers of divided Cartesian grids");
31 15428 : params.addParam<Point>("widths", "Widths in the X, Y and Z directions");
32 15428 : params.addRequiredRangeCheckedParam<unsigned int>("nx", "nx>0", "Number of divisions in X");
33 15428 : params.addRequiredRangeCheckedParam<unsigned int>("ny", "ny>0", "Number of divisions in Y");
34 15428 : params.addRequiredRangeCheckedParam<unsigned int>("nz", "nz>0", "Number of divisions in Z");
35 46284 : params.addParam<bool>(
36 : "assign_domain_outside_grid_to_border",
37 30856 : false,
38 : "Whether to map the domain outside the grid back to the border of the grid");
39 :
40 15428 : return params;
41 0 : }
42 :
43 623 : CartesianGridDivision::CartesianGridDivision(const InputParameters & parameters)
44 : : MeshDivision(parameters),
45 623 : _bottom_left(isParamValid("bottom_left") ? getParam<Point>("bottom_left") : Point(0, 0, 0)),
46 623 : _top_right(isParamValid("top_right") ? getParam<Point>("top_right") : Point(0, 0, 0)),
47 623 : _center(isParamValid("center") ? &getParam<Point>("center") : nullptr),
48 623 : _center_positions(
49 1246 : isParamValid("center_positions")
50 623 : ? &_fe_problem->getPositionsObject(getParam<PositionsName>("center_positions"))
51 : : nullptr),
52 623 : _widths(isParamValid("widths") ? getParam<Point>("widths") : Point(_top_right - _bottom_left)),
53 623 : _nx(getParam<unsigned int>("nx")),
54 623 : _ny(getParam<unsigned int>("ny")),
55 623 : _nz(getParam<unsigned int>("nz")),
56 1246 : _outside_grid_counts_as_border(getParam<bool>("assign_domain_outside_grid_to_border"))
57 : {
58 623 : CartesianGridDivision::initialize();
59 :
60 : // Check non-overlapping inputs for the dimensions of the grid
61 619 : 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 619 : 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 619 : 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 619 : 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 619 : 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 619 : 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 619 : 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 619 : 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 619 : 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 619 : }
91 :
92 : void
93 623 : CartesianGridDivision::initialize()
94 : {
95 623 : if (!_center_positions)
96 606 : setNumDivisions(_nx * _ny * _nz);
97 : else
98 17 : setNumDivisions(_center_positions->getNumPositions() * _nx * _ny * _nz);
99 :
100 : // Check that the grid is well-defined
101 623 : if (_center_positions)
102 : {
103 17 : Real min_dist = _widths.norm();
104 17 : 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 17 : if (MooseUtils::absoluteFuzzyGreaterThan(min_dist, min_center_dist))
108 4 : 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 619 : _mesh_fully_indexed = true;
116 619 : if (!_outside_grid_counts_as_border)
117 632 : for (auto i : make_range(LIBMESH_DIM))
118 : {
119 632 : if (_center_positions)
120 : {
121 13 : _mesh_fully_indexed = false;
122 13 : break;
123 : }
124 663 : else if (_bottom_left(i) > _mesh.getInflatedProcessorBoundingBox(0).first(i) ||
125 663 : _top_right(i) < _mesh.getInflatedProcessorBoundingBox(0).second(i))
126 : {
127 606 : _mesh_fully_indexed = false;
128 606 : break;
129 : }
130 : }
131 619 : }
132 :
133 : unsigned int
134 85112 : CartesianGridDivision::divisionIndex(const Elem & elem) const
135 : {
136 85112 : return divisionIndex(elem.vertex_average());
137 : }
138 :
139 : unsigned int
140 187882 : CartesianGridDivision::divisionIndex(const Point & pt) const
141 : {
142 187882 : unsigned int offset = 0;
143 : // Determine the local grid bounds
144 187882 : Point bottom_left, top_right, p;
145 187882 : 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 187242 : bottom_left = _bottom_left;
161 187242 : top_right = _top_right;
162 187242 : p = pt;
163 : }
164 :
165 187882 : if (!_outside_grid_counts_as_border)
166 : {
167 318115 : if (MooseUtils::absoluteFuzzyLessThan(p(0), bottom_left(0)) ||
168 130233 : MooseUtils::absoluteFuzzyGreaterThan(p(0), top_right(0)))
169 72221 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
170 198172 : if (MooseUtils::absoluteFuzzyLessThan(p(1), bottom_left(1)) ||
171 82511 : MooseUtils::absoluteFuzzyGreaterThan(p(1), top_right(1)))
172 40979 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
173 149076 : if (MooseUtils::absoluteFuzzyLessThan(p(2), bottom_left(2)) ||
174 74394 : MooseUtils::absoluteFuzzyGreaterThan(p(2), top_right(2)))
175 864 : return MooseMeshDivision::INVALID_DIVISION_INDEX;
176 : }
177 :
178 73818 : const auto not_found = MooseMeshDivision::INVALID_DIVISION_INDEX;
179 73818 : auto ix = not_found, iy = not_found, iz = not_found;
180 :
181 : // Look inside the grid and on the left / back / bottom
182 279334 : for (const auto jx : make_range(_nx + 1))
183 : {
184 279334 : const auto border_x = bottom_left(0) + _widths(0) * jx / _nx;
185 279334 : 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 279334 : if (border_x >= p(0))
190 : {
191 73818 : ix = (jx > 0) ? jx - 1 : 0;
192 73818 : break;
193 : }
194 : }
195 278988 : for (const auto jy : make_range(_ny + 1))
196 : {
197 278988 : const auto border_y = bottom_left(1) + _widths(1) * jy / _ny;
198 278988 : 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 278988 : if (border_y >= p(1))
203 : {
204 73818 : iy = (jy > 0) ? jy - 1 : 0;
205 73818 : break;
206 : }
207 : }
208 74106 : for (const auto jz : make_range(_nz + 1))
209 : {
210 74106 : const auto border_z = bottom_left(2) + _widths(2) * jz / _nz;
211 74106 : 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 74106 : if (border_z >= p(2))
216 : {
217 73818 : iz = (jz > 0) ? jz - 1 : 0;
218 73818 : break;
219 : }
220 : }
221 :
222 : // Look on the right / front / top of the grid
223 73818 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(0), top_right(0)))
224 0 : ix = _nx - 1;
225 73818 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(1), top_right(1)))
226 24 : iy = _ny - 1;
227 73818 : if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), top_right(2)))
228 73530 : iz = _nz - 1;
229 :
230 : // Handle edge case on widths
231 73818 : if (ix == not_found && MooseUtils::absoluteFuzzyEqual(_widths(0), 0))
232 0 : ix = 0;
233 73818 : if (iy == not_found && MooseUtils::absoluteFuzzyEqual(_widths(1), 0))
234 0 : iy = 0;
235 73818 : 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 73818 : return offset + ix + _nx * iy + iz * _nx * _ny;
242 : }
|
