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 118416 : GeneratedMesh::validParams()
29 : {
30 118416 : InputParameters params = MooseMesh::validParams();
31 :
32 118416 : MooseEnum elem_types(LIST_GEOM_ELEM); // no default
33 :
34 118416 : MooseEnum dims("1=1 2 3");
35 118416 : params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
36 :
37 118416 : params.addParam<unsigned int>("nx", 1, "Number of elements in the X direction");
38 118416 : params.addParam<unsigned int>("ny", 1, "Number of elements in the Y direction");
39 118416 : params.addParam<unsigned int>("nz", 1, "Number of elements in the Z direction");
40 118416 : params.addParam<Real>("xmin", 0.0, "Lower X Coordinate of the generated mesh");
41 118416 : params.addParam<Real>("ymin", 0.0, "Lower Y Coordinate of the generated mesh");
42 118416 : params.addParam<Real>("zmin", 0.0, "Lower Z Coordinate of the generated mesh");
43 118416 : params.addParam<Real>("xmax", 1.0, "Upper X Coordinate of the generated mesh");
44 118416 : params.addParam<Real>("ymax", 1.0, "Upper Y Coordinate of the generated mesh");
45 118416 : params.addParam<Real>("zmax", 1.0, "Upper Z Coordinate of the generated mesh");
46 118416 : 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 355248 : params.addParam<bool>(
52 : "gauss_lobatto_grid",
53 236832 : false,
54 : "Grade mesh into boundaries according to Gauss-Lobatto quadrature spacing.");
55 355248 : params.addRangeCheckedParam<Real>(
56 : "bias_x",
57 236832 : 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 355248 : params.addRangeCheckedParam<Real>(
61 : "bias_y",
62 236832 : 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 355248 : params.addRangeCheckedParam<Real>(
66 : "bias_z",
67 236832 : 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 118416 : params.addParamNamesToGroup("dim", "Required");
72 :
73 118416 : params.addClassDescription(
74 : "Create a line, square, or cube mesh with uniformly spaced or biased elements.");
75 236832 : return params;
76 118416 : }
77 :
78 36276 : GeneratedMesh::GeneratedMesh(const InputParameters & parameters)
79 : : MooseMesh(parameters),
80 36276 : _dim(getParam<MooseEnum>("dim")),
81 36276 : _nx(getParam<unsigned int>("nx")),
82 36276 : _ny(getParam<unsigned int>("ny")),
83 36276 : _nz(getParam<unsigned int>("nz")),
84 36276 : _xmin(getParam<Real>("xmin")),
85 36276 : _xmax(getParam<Real>("xmax")),
86 36276 : _ymin(getParam<Real>("ymin")),
87 36276 : _ymax(getParam<Real>("ymax")),
88 36276 : _zmin(getParam<Real>("zmin")),
89 36276 : _zmax(getParam<Real>("zmax")),
90 36276 : _gauss_lobatto_grid(getParam<bool>("gauss_lobatto_grid")),
91 36276 : _bias_x(getParam<Real>("bias_x")),
92 36276 : _bias_y(getParam<Real>("bias_y")),
93 36276 : _bias_z(getParam<Real>("bias_z")),
94 36276 : _dims_may_have_changed(false)
95 : {
96 36276 : 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 36276 : _regular_orthogonal_mesh = true;
101 36276 : }
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 43367 : GeneratedMesh::getMinInDimension(unsigned int component) const
115 : {
116 43367 : if (_dims_may_have_changed)
117 0 : return MooseMesh::getMinInDimension(component);
118 :
119 43367 : switch (component)
120 : {
121 894 : case 0:
122 894 : return _xmin;
123 20744 : case 1:
124 20744 : return _dim > 1 ? _ymin : 0;
125 21729 : case 2:
126 21729 : return _dim > 2 ? _zmin : 0;
127 0 : default:
128 0 : mooseError("Invalid component");
129 : }
130 : }
131 :
132 : Real
133 43301 : GeneratedMesh::getMaxInDimension(unsigned int component) const
134 : {
135 43301 : if (_dims_may_have_changed)
136 0 : return MooseMesh::getMaxInDimension(component);
137 :
138 43301 : switch (component)
139 : {
140 872 : case 0:
141 872 : return _xmax;
142 20722 : case 1:
143 20722 : return _dim > 1 ? _ymax : 0;
144 21707 : case 2:
145 21707 : return _dim > 2 ? _zmax : 0;
146 0 : default:
147 0 : mooseError("Invalid component");
148 : }
149 : }
150 :
151 : std::unique_ptr<MooseMesh>
152 1943 : GeneratedMesh::safeClone() const
153 : {
154 1943 : return _app.getFactory().copyConstruct(*this);
155 : }
156 :
157 : void
158 34435 : GeneratedMesh::buildMesh()
159 : {
160 34435 : MooseEnum elem_type_enum = getParam<MooseEnum>("elem_type");
161 :
162 34435 : if (!isParamValid("elem_type"))
163 : {
164 : // Switching on MooseEnum
165 28694 : switch (_dim)
166 : {
167 5222 : case 1:
168 5222 : elem_type_enum = "EDGE2";
169 5222 : break;
170 22275 : case 2:
171 22275 : elem_type_enum = "QUAD4";
172 22275 : break;
173 1197 : case 3:
174 1197 : elem_type_enum = "HEX8";
175 1197 : break;
176 : }
177 : }
178 :
179 34435 : ElemType elem_type = Utility::string_to_enum<ElemType>(elem_type_enum);
180 :
181 : // Switching on MooseEnum
182 34435 : switch (_dim)
183 : {
184 : // The build_XYZ mesh generation functions take an
185 : // UnstructuredMesh& as the first argument, hence the dynamic_cast.
186 5577 : case 1:
187 5577 : MeshTools::Generation::build_line(dynamic_cast<UnstructuredMesh &>(getMesh()),
188 : _nx,
189 : _xmin,
190 : _xmax,
191 : elem_type,
192 5577 : _gauss_lobatto_grid);
193 5577 : break;
194 27291 : case 2:
195 27291 : MeshTools::Generation::build_square(dynamic_cast<UnstructuredMesh &>(getMesh()),
196 : _nx,
197 : _ny,
198 : _xmin,
199 : _xmax,
200 : _ymin,
201 : _ymax,
202 : elem_type,
203 27291 : _gauss_lobatto_grid);
204 27291 : break;
205 1567 : case 3:
206 1567 : 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 1567 : _gauss_lobatto_grid);
218 1567 : break;
219 : }
220 :
221 : // Apply the bias if any exists
222 34435 : if (_bias_x != 1.0 || _bias_y != 1.0 || _bias_z != 1.0)
223 : {
224 : // Reference to the libmesh mesh
225 22 : MeshBase & mesh = getMesh();
226 :
227 : // Biases
228 : std::array<Real, LIBMESH_DIM> bias = {
229 22 : {_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 22 : {_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 22 : {getMinInDimension(0), getMinInDimension(1), getMinInDimension(2)}};
238 :
239 : // Number of elements in each direction.
240 22 : 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 22 : std::array<std::vector<Real>, LIBMESH_DIM> pows;
245 88 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
246 : {
247 66 : pows[dir].resize(nelem[dir] + 1);
248 462 : for (unsigned int i = 0; i < pows[dir].size(); ++i)
249 396 : 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 16438 : for (auto & node_ptr : mesh.node_ptr_range())
254 : {
255 16416 : Node & node = *node_ptr;
256 :
257 65664 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
258 : {
259 49248 : 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 32832 : Real float_index = (node(dir) - mins[dir]) * nelem[dir] / width[dir];
265 :
266 32832 : Real integer_part = 0;
267 32832 : Real fractional_part = std::modf(float_index, &integer_part);
268 :
269 : // Figure out where to move the node...
270 32832 : 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 20029 : 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 20029 : node(dir) =
286 20029 : mins[dir] + width[dir] * (1. - pows[dir][index]) / (1. - pows[dir][nelem[dir]]);
287 : }
288 12803 : 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 25606 : node(dir) = mins[dir] +
300 12803 : width[dir] *
301 12803 : (1. - 0.5 * (pows[dir][integer_part] + pows[dir][integer_part + 1])) /
302 12803 : (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 22 : }
312 22 : }
313 34435 : }
|
