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 "GeneratedMesh.h"
11 :
12 : #include "MooseApp.h"
13 :
14 : #include "libmesh/mesh_generation.h"
15 : #include "libmesh/string_to_enum.h"
16 : #include "libmesh/periodic_boundaries.h"
17 : #include "libmesh/periodic_boundary_base.h"
18 : #include "libmesh/unstructured_mesh.h"
19 : #include "libmesh/node.h"
20 :
21 : // C++ includes
22 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
23 : #include <array>
24 :
25 : registerMooseObject("MooseApp", GeneratedMesh);
26 :
27 : InputParameters
28 113335 : GeneratedMesh::validParams()
29 : {
30 113335 : InputParameters params = MooseMesh::validParams();
31 :
32 113335 : MooseEnum elem_types(LIST_GEOM_ELEM); // no default
33 :
34 113335 : MooseEnum dims("1=1 2 3");
35 113335 : params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
36 :
37 113335 : params.addParam<unsigned int>("nx", 1, "Number of elements in the X direction");
38 113335 : params.addParam<unsigned int>("ny", 1, "Number of elements in the Y direction");
39 113335 : params.addParam<unsigned int>("nz", 1, "Number of elements in the Z direction");
40 113335 : params.addParam<Real>("xmin", 0.0, "Lower X Coordinate of the generated mesh");
41 113335 : params.addParam<Real>("ymin", 0.0, "Lower Y Coordinate of the generated mesh");
42 113335 : params.addParam<Real>("zmin", 0.0, "Lower Z Coordinate of the generated mesh");
43 113335 : params.addParam<Real>("xmax", 1.0, "Upper X Coordinate of the generated mesh");
44 113335 : params.addParam<Real>("ymax", 1.0, "Upper Y Coordinate of the generated mesh");
45 113335 : params.addParam<Real>("zmax", 1.0, "Upper Z Coordinate of the generated mesh");
46 113335 : params.addParam<MooseEnum>("elem_type",
47 : elem_types,
48 : "The type of element from libMesh to "
49 : "generate (default: linear element for "
50 : "requested dimension)");
51 340005 : params.addParam<bool>(
52 : "gauss_lobatto_grid",
53 226670 : false,
54 : "Grade mesh into boundaries according to Gauss-Lobatto quadrature spacing.");
55 340005 : params.addRangeCheckedParam<Real>(
56 : "bias_x",
57 226670 : 1.,
58 : "bias_x>=0.5 & bias_x<=2",
59 : "The amount by which to grow (or shrink) the cells in the x-direction.");
60 340005 : params.addRangeCheckedParam<Real>(
61 : "bias_y",
62 226670 : 1.,
63 : "bias_y>=0.5 & bias_y<=2",
64 : "The amount by which to grow (or shrink) the cells in the y-direction.");
65 340005 : params.addRangeCheckedParam<Real>(
66 : "bias_z",
67 226670 : 1.,
68 : "bias_z>=0.5 & bias_z<=2",
69 : "The amount by which to grow (or shrink) the cells in the z-direction.");
70 :
71 113335 : params.addParamNamesToGroup("dim", "Required");
72 :
73 113335 : params.addClassDescription(
74 : "Create a line, square, or cube mesh with uniformly spaced or biased elements.");
75 226670 : return params;
76 113335 : }
77 :
78 33853 : GeneratedMesh::GeneratedMesh(const InputParameters & parameters)
79 : : MooseMesh(parameters),
80 33853 : _dim(getParam<MooseEnum>("dim")),
81 33853 : _nx(getParam<unsigned int>("nx")),
82 33853 : _ny(getParam<unsigned int>("ny")),
83 33853 : _nz(getParam<unsigned int>("nz")),
84 33853 : _xmin(getParam<Real>("xmin")),
85 33853 : _xmax(getParam<Real>("xmax")),
86 33853 : _ymin(getParam<Real>("ymin")),
87 33853 : _ymax(getParam<Real>("ymax")),
88 33853 : _zmin(getParam<Real>("zmin")),
89 33853 : _zmax(getParam<Real>("zmax")),
90 33853 : _gauss_lobatto_grid(getParam<bool>("gauss_lobatto_grid")),
91 33853 : _bias_x(getParam<Real>("bias_x")),
92 33853 : _bias_y(getParam<Real>("bias_y")),
93 33853 : _bias_z(getParam<Real>("bias_z")),
94 33853 : _dims_may_have_changed(false)
95 : {
96 33853 : if (_gauss_lobatto_grid && (_bias_x != 1.0 || _bias_y != 1.0 || _bias_z != 1.0))
97 0 : mooseError("Cannot apply both Gauss-Lobatto mesh grading and biasing at the same time.");
98 :
99 : // All generated meshes are regular orthogonal meshes - until they get modified ;)
100 33853 : _regular_orthogonal_mesh = true;
101 33853 : }
102 :
103 : void
104 0 : GeneratedMesh::prepared(bool state)
105 : {
106 0 : MooseMesh::prepared(state);
107 :
108 : // Fall back on scanning the mesh for coordinates instead of using input parameters for queries
109 0 : if (!state)
110 0 : _dims_may_have_changed = true;
111 0 : }
112 :
113 : Real
114 39776 : GeneratedMesh::getMinInDimension(unsigned int component) const
115 : {
116 39776 : if (_dims_may_have_changed)
117 0 : return MooseMesh::getMinInDimension(component);
118 :
119 39776 : switch (component)
120 : {
121 824 : case 0:
122 824 : return _xmin;
123 19031 : case 1:
124 19031 : return _dim > 1 ? _ymin : 0;
125 19921 : case 2:
126 19921 : return _dim > 2 ? _zmin : 0;
127 0 : default:
128 0 : mooseError("Invalid component");
129 : }
130 : }
131 :
132 : Real
133 39716 : GeneratedMesh::getMaxInDimension(unsigned int component) const
134 : {
135 39716 : if (_dims_may_have_changed)
136 0 : return MooseMesh::getMaxInDimension(component);
137 :
138 39716 : switch (component)
139 : {
140 804 : case 0:
141 804 : return _xmax;
142 19011 : case 1:
143 19011 : return _dim > 1 ? _ymax : 0;
144 19901 : case 2:
145 19901 : return _dim > 2 ? _zmax : 0;
146 0 : default:
147 0 : mooseError("Invalid component");
148 : }
149 : }
150 :
151 : std::unique_ptr<MooseMesh>
152 1793 : GeneratedMesh::safeClone() const
153 : {
154 1793 : return _app.getFactory().copyConstruct(*this);
155 : }
156 :
157 : void
158 32017 : GeneratedMesh::buildMesh()
159 : {
160 32017 : MooseEnum elem_type_enum = getParam<MooseEnum>("elem_type");
161 :
162 32017 : if (!isParamValid("elem_type"))
163 : {
164 : // Switching on MooseEnum
165 26647 : switch (_dim)
166 : {
167 4869 : case 1:
168 4869 : elem_type_enum = "EDGE2";
169 4869 : break;
170 20675 : case 2:
171 20675 : elem_type_enum = "QUAD4";
172 20675 : break;
173 1103 : case 3:
174 1103 : elem_type_enum = "HEX8";
175 1103 : break;
176 : }
177 : }
178 :
179 32017 : ElemType elem_type = Utility::string_to_enum<ElemType>(elem_type_enum);
180 :
181 : // Switching on MooseEnum
182 32017 : switch (_dim)
183 : {
184 : // The build_XYZ mesh generation functions take an
185 : // UnstructuredMesh& as the first argument, hence the dynamic_cast.
186 5204 : case 1:
187 5204 : MeshTools::Generation::build_line(dynamic_cast<UnstructuredMesh &>(getMesh()),
188 : _nx,
189 : _xmin,
190 : _xmax,
191 : elem_type,
192 5204 : _gauss_lobatto_grid);
193 5204 : break;
194 25371 : case 2:
195 25371 : MeshTools::Generation::build_square(dynamic_cast<UnstructuredMesh &>(getMesh()),
196 : _nx,
197 : _ny,
198 : _xmin,
199 : _xmax,
200 : _ymin,
201 : _ymax,
202 : elem_type,
203 25371 : _gauss_lobatto_grid);
204 25371 : break;
205 1442 : case 3:
206 1442 : MeshTools::Generation::build_cube(dynamic_cast<UnstructuredMesh &>(getMesh()),
207 : _nx,
208 : _ny,
209 : _nz,
210 : _xmin,
211 : _xmax,
212 : _ymin,
213 : _ymax,
214 : _zmin,
215 : _zmax,
216 : elem_type,
217 1442 : _gauss_lobatto_grid);
218 1442 : break;
219 : }
220 :
221 : // Apply the bias if any exists
222 32017 : if (_bias_x != 1.0 || _bias_y != 1.0 || _bias_z != 1.0)
223 : {
224 : // Reference to the libmesh mesh
225 20 : MeshBase & mesh = getMesh();
226 :
227 : // Biases
228 : std::array<Real, LIBMESH_DIM> bias = {
229 20 : {_bias_x, _dim > 1 ? _bias_y : 1.0, _dim > 2 ? _bias_z : 1.0}};
230 :
231 : // "width" of the mesh in each direction
232 : std::array<Real, LIBMESH_DIM> width = {
233 20 : {_xmax - _xmin, _dim > 1 ? _ymax - _ymin : 0, _dim > 2 ? _zmax - _zmin : 0}};
234 :
235 : // Min mesh extent in each direction.
236 : std::array<Real, LIBMESH_DIM> mins = {
237 20 : {getMinInDimension(0), getMinInDimension(1), getMinInDimension(2)}};
238 :
239 : // Number of elements in each direction.
240 20 : std::array<unsigned int, LIBMESH_DIM> nelem = {{_nx, _dim > 1 ? _ny : 1, _dim > 2 ? _nz : 1}};
241 :
242 : // We will need the biases raised to integer powers in each
243 : // direction, so let's pre-compute those...
244 20 : std::array<std::vector<Real>, LIBMESH_DIM> pows;
245 80 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
246 : {
247 60 : pows[dir].resize(nelem[dir] + 1);
248 420 : for (unsigned int i = 0; i < pows[dir].size(); ++i)
249 360 : pows[dir][i] = std::pow(bias[dir], static_cast<int>(i));
250 : }
251 :
252 : // Loop over the nodes and move them to the desired location
253 14889 : for (auto & node_ptr : mesh.node_ptr_range())
254 : {
255 14869 : Node & node = *node_ptr;
256 :
257 59476 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
258 : {
259 44607 : if (width[dir] != 0. && bias[dir] != 1.)
260 : {
261 : // Compute the scaled "index" of the current point. This
262 : // will either be close to a whole integer or a whole
263 : // integer+0.5 for quadratic nodes.
264 29738 : Real float_index = (node(dir) - mins[dir]) * nelem[dir] / width[dir];
265 :
266 29738 : Real integer_part = 0;
267 29738 : Real fractional_part = std::modf(float_index, &integer_part);
268 :
269 : // Figure out where to move the node...
270 29738 : if (std::abs(fractional_part) < TOLERANCE || std::abs(fractional_part - 1.0) < TOLERANCE)
271 : {
272 : // If the fractional_part ~ 0.0 or 1.0, this is a vertex node, so
273 : // we don't need an average.
274 : //
275 : // Compute the "index" we are at in the current direction. We
276 : // round to the nearest integral value to get this instead
277 : // of using "integer_part", since that could be off by a
278 : // lot (e.g. we want 3.9999 to map to 4.0 instead of 3.0).
279 18145 : int index = round(float_index);
280 :
281 : mooseAssert(index >= static_cast<int>(0) && index < static_cast<int>(pows[dir].size()),
282 : "Scaled \"index\" out of range");
283 :
284 : // Move node to biased location.
285 18145 : node(dir) =
286 18145 : mins[dir] + width[dir] * (1. - pows[dir][index]) / (1. - pows[dir][nelem[dir]]);
287 : }
288 11593 : else if (std::abs(fractional_part - 0.5) < TOLERANCE)
289 : {
290 : // If the fractional_part ~ 0.5, this is a midedge/face
291 : // (i.e. quadratic) node. We don't move those with the same
292 : // bias as the vertices, instead we put them midway between
293 : // their respective vertices.
294 : //
295 : // Also, since the fractional part is nearly 0.5, we know that
296 : // the integer_part will be the index of the vertex to the
297 : // left, and integer_part+1 will be the index of the
298 : // vertex to the right.
299 23186 : node(dir) = mins[dir] +
300 11593 : width[dir] *
301 11593 : (1. - 0.5 * (pows[dir][integer_part] + pows[dir][integer_part + 1])) /
302 11593 : (1. - pows[dir][nelem[dir]]);
303 : }
304 : else
305 : {
306 : // We don't yet handle anything higher order than quadratic...
307 0 : mooseError("Unable to bias node at node(", dir, ")=", node(dir));
308 : }
309 : }
310 : }
311 20 : }
312 20 : }
313 32017 : }
|
