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 "GeneratedMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 :
13 : #include "libmesh/replicated_mesh.h"
14 : #include "libmesh/mesh_generation.h"
15 : #include "libmesh/mesh_modification.h"
16 : #include "libmesh/string_to_enum.h"
17 : #include "libmesh/periodic_boundaries.h"
18 : #include "libmesh/periodic_boundary_base.h"
19 : #include "libmesh/unstructured_mesh.h"
20 : #include "libmesh/elem.h"
21 :
22 : // C++ includes
23 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
24 :
25 : registerMooseObject("MooseApp", GeneratedMeshGenerator);
26 :
27 : InputParameters
28 49459 : GeneratedMeshGenerator::validParams()
29 : {
30 49459 : InputParameters params = MeshGenerator::validParams();
31 :
32 98918 : MooseEnum elem_types(LIST_GEOM_ELEM); // no default
33 :
34 197836 : MooseEnum dims("1=1 2 3");
35 197836 : params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
36 :
37 197836 : params.addParam<unsigned int>("nx", 1, "Number of elements in the X direction");
38 197836 : params.addParam<unsigned int>("ny", 1, "Number of elements in the Y direction");
39 197836 : params.addParam<unsigned int>("nz", 1, "Number of elements in the Z direction");
40 197836 : params.addParam<Real>("xmin", 0.0, "Lower X Coordinate of the generated mesh");
41 197836 : params.addParam<Real>("ymin", 0.0, "Lower Y Coordinate of the generated mesh");
42 197836 : params.addParam<Real>("zmin", 0.0, "Lower Z Coordinate of the generated mesh");
43 197836 : params.addParam<Real>("xmax", 1.0, "Upper X Coordinate of the generated mesh");
44 197836 : params.addParam<Real>("ymax", 1.0, "Upper Y Coordinate of the generated mesh");
45 197836 : params.addParam<Real>("zmax", 1.0, "Upper Z Coordinate of the generated mesh");
46 197836 : 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 197836 : params.addParam<std::vector<SubdomainID>>(
52 : "subdomain_ids",
53 : "Subdomain IDs for each element, default to all zero. If a single number is specified, that "
54 : "subdomain id is used for all elements.");
55 197836 : params.addParam<SubdomainName>("subdomain_name",
56 : "If specified, single subdomain name for all elements");
57 :
58 148377 : params.addParam<bool>(
59 : "gauss_lobatto_grid",
60 98918 : false,
61 : "Grade mesh into boundaries according to Gauss-Lobatto quadrature spacing.");
62 247295 : params.addRangeCheckedParam<Real>(
63 : "bias_x",
64 98918 : 1.,
65 : "bias_x>=0.5 & bias_x<=2",
66 : "The amount by which to grow (or shrink) the cells in the x-direction.");
67 247295 : params.addRangeCheckedParam<Real>(
68 : "bias_y",
69 98918 : 1.,
70 : "bias_y>=0.5 & bias_y<=2",
71 : "The amount by which to grow (or shrink) the cells in the y-direction.");
72 247295 : params.addRangeCheckedParam<Real>(
73 : "bias_z",
74 98918 : 1.,
75 : "bias_z>=0.5 & bias_z<=2",
76 : "The amount by which to grow (or shrink) the cells in the z-direction.");
77 :
78 197836 : params.addParam<std::string>("boundary_name_prefix",
79 : "If provided, prefix the built in boundary names with this string");
80 148377 : params.addParam<boundary_id_type>(
81 98918 : "boundary_id_offset", 0, "This offset is added to the generated boundary IDs");
82 :
83 197836 : params.addParam<std::vector<ExtraElementIDName>>("extra_element_integers",
84 : "Names of extra element integers");
85 :
86 49459 : params.addClassDescription(
87 : "Create a line, square, or cube mesh with uniformly spaced or biased elements.");
88 :
89 98918 : return params;
90 49459 : }
91 :
92 20072 : GeneratedMeshGenerator::GeneratedMeshGenerator(const InputParameters & parameters)
93 : : MeshGenerator(parameters),
94 20066 : _dim(getParam<MooseEnum>("dim")),
95 80264 : _nx(declareMeshProperty("num_elements_x", getParam<unsigned int>("nx"))),
96 80264 : _ny(declareMeshProperty("num_elements_y", getParam<unsigned int>("ny"))),
97 80264 : _nz(declareMeshProperty("num_elements_z", getParam<unsigned int>("nz"))),
98 80264 : _xmin(declareMeshProperty("xmin", getParam<Real>("xmin"))),
99 80264 : _xmax(declareMeshProperty("xmax", getParam<Real>("xmax"))),
100 80264 : _ymin(declareMeshProperty("ymin", getParam<Real>("ymin"))),
101 80264 : _ymax(declareMeshProperty("ymax", getParam<Real>("ymax"))),
102 80264 : _zmin(declareMeshProperty("zmin", getParam<Real>("zmin"))),
103 80264 : _zmax(declareMeshProperty("zmax", getParam<Real>("zmax"))),
104 40132 : _has_subdomain_ids(isParamValid("subdomain_ids")),
105 40132 : _gauss_lobatto_grid(getParam<bool>("gauss_lobatto_grid")),
106 40132 : _bias_x(getParam<Real>("bias_x")),
107 40132 : _bias_y(getParam<Real>("bias_y")),
108 40132 : _bias_z(getParam<Real>("bias_z")),
109 100330 : _boundary_name_prefix(isParamValid("boundary_name_prefix")
110 20066 : ? getParam<std::string>("boundary_name_prefix") + "_"
111 : : ""),
112 60204 : _boundary_id_offset(getParam<boundary_id_type>("boundary_id_offset"))
113 : {
114 20066 : if (_gauss_lobatto_grid && (_bias_x != 1.0 || _bias_y != 1.0 || _bias_z != 1.0))
115 0 : mooseError("Cannot apply both Gauss-Lobatto mesh grading and biasing at the same time.");
116 20066 : if (_xmax < _xmin)
117 6 : paramError("xmax", "xmax must be larger than xmin.");
118 20063 : if (_ymax < _ymin)
119 6 : paramError("ymax", "ymax must be larger than ymin.");
120 20060 : if (_zmax < _zmin)
121 6 : paramError("zmax", "zmax must be larger than zmin.");
122 20057 : }
123 :
124 : std::unique_ptr<MeshBase>
125 19072 : GeneratedMeshGenerator::generate()
126 : {
127 : // Have MOOSE construct the correct libMesh::Mesh object using Mesh block and CLI parameters.
128 19072 : auto mesh = buildMeshBaseObject(_dim);
129 :
130 57216 : if (isParamValid("extra_element_integers"))
131 : {
132 798 : for (auto & name : getParam<std::vector<ExtraElementIDName>>("extra_element_integers"))
133 222 : mesh->add_elem_integer(name);
134 : }
135 :
136 38144 : MooseEnum elem_type_enum = getParam<MooseEnum>("elem_type");
137 :
138 57216 : if (!isParamValid("elem_type"))
139 : {
140 : // Switching on MooseEnum
141 16251 : switch (_dim)
142 : {
143 4812 : case 1:
144 4812 : elem_type_enum = "EDGE2";
145 4812 : break;
146 9343 : case 2:
147 9343 : elem_type_enum = "QUAD4";
148 9343 : break;
149 2096 : case 3:
150 2096 : elem_type_enum = "HEX8";
151 2096 : break;
152 : }
153 : }
154 :
155 19072 : ElemType elem_type = Utility::string_to_enum<ElemType>(elem_type_enum);
156 :
157 : // Switching on MooseEnum
158 19072 : switch (_dim)
159 : {
160 : // The build_XYZ mesh generation functions take an
161 : // UnstructuredMesh& as the first argument, hence the static_cast.
162 4892 : case 1:
163 4892 : MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh),
164 4892 : _nx,
165 4892 : _xmin,
166 4892 : _xmax,
167 : elem_type,
168 4892 : _gauss_lobatto_grid);
169 4892 : break;
170 10902 : case 2:
171 10902 : MeshTools::Generation::build_square(static_cast<UnstructuredMesh &>(*mesh),
172 10902 : _nx,
173 10902 : _ny,
174 10902 : _xmin,
175 10902 : _xmax,
176 10902 : _ymin,
177 10902 : _ymax,
178 : elem_type,
179 10902 : _gauss_lobatto_grid);
180 10902 : break;
181 3278 : case 3:
182 3278 : MeshTools::Generation::build_cube(static_cast<UnstructuredMesh &>(*mesh),
183 3278 : _nx,
184 3278 : _ny,
185 3278 : _nz,
186 3278 : _xmin,
187 3278 : _xmax,
188 3278 : _ymin,
189 3278 : _ymax,
190 3278 : _zmin,
191 3278 : _zmax,
192 : elem_type,
193 3278 : _gauss_lobatto_grid);
194 3278 : break;
195 : }
196 :
197 19072 : if (_has_subdomain_ids)
198 : {
199 3136 : auto & bids = getParam<std::vector<SubdomainID>>("subdomain_ids");
200 1568 : if (bids.size() != _nx * _ny * _nz && bids.size() != 1)
201 0 : paramError("subdomain_ids",
202 : "Size must equal to the product of number of elements in all directions, or one.");
203 19354 : for (auto & elem : mesh->element_ptr_range())
204 : {
205 17786 : const Point p = elem->vertex_average();
206 17786 : unsigned int ix = std::floor((p(0) - _xmin) / (_xmax - _xmin) * _nx);
207 17786 : unsigned int iy = std::floor((p(1) - _ymin) / (_ymax - _ymin) * _ny);
208 17786 : unsigned int iz = std::floor((p(2) - _zmin) / (_zmax - _zmin) * _nz);
209 17786 : unsigned int i = iz * _nx * _ny + iy * _nx + ix;
210 17786 : if (bids.size() == 1)
211 6847 : elem->subdomain_id() = bids[0];
212 : else
213 10939 : elem->subdomain_id() = bids[i];
214 1568 : }
215 : }
216 :
217 57216 : if (isParamValid("subdomain_name"))
218 : {
219 750 : const auto & subdomain_name = getParam<SubdomainName>("subdomain_name");
220 1125 : if (isParamValid("subdomain_ids"))
221 : {
222 160 : const auto & bids = getParam<std::vector<SubdomainID>>("subdomain_ids");
223 80 : if (bids.size() > 1)
224 0 : paramError(
225 : "subdomain_ids",
226 : "Specifying a subdomain_name is only supported for a single entry in subdomain_ids");
227 : else
228 80 : mesh->subdomain_name(bids[0]) = subdomain_name;
229 : }
230 : else
231 295 : mesh->subdomain_name(0) = subdomain_name;
232 : }
233 :
234 : // rename and shift boundaries
235 19072 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
236 :
237 : // Copy, since we're modifying the container mid-iteration
238 19072 : const auto mesh_boundary_ids = boundary_info.get_global_boundary_ids();
239 92132 : for (auto rit = mesh_boundary_ids.rbegin(); rit != mesh_boundary_ids.rend(); ++rit)
240 : {
241 73060 : const std::string old_sideset_name = boundary_info.sideset_name(*rit);
242 73060 : const std::string old_nodeset_name = boundary_info.nodeset_name(*rit);
243 :
244 73060 : MeshTools::Modification::change_boundary_id(*mesh, *rit, *rit + _boundary_id_offset);
245 :
246 73060 : boundary_info.sideset_name(*rit + _boundary_id_offset) =
247 146120 : _boundary_name_prefix + old_sideset_name;
248 73060 : boundary_info.nodeset_name(*rit + _boundary_id_offset) =
249 146120 : _boundary_name_prefix + old_nodeset_name;
250 73060 : }
251 :
252 : // FIXME: change_boundary_id() was *supposed* to leave the mesh in
253 : // a fully prepared state, but there's a libMesh bug to fix first.
254 19072 : mesh->unset_has_boundary_id_sets();
255 :
256 : // Apply the bias if any exists
257 19072 : if (_bias_x != 1.0 || _bias_y != 1.0 || _bias_z != 1.0)
258 : {
259 156 : const auto MIN = std::numeric_limits<Real>::max();
260 :
261 : // Biases
262 : std::array<Real, LIBMESH_DIM> bias = {
263 156 : {_bias_x, _dim > 1 ? _bias_y : 1.0, _dim > 2 ? _bias_z : 1.0}};
264 :
265 : // "width" of the mesh in each direction
266 : std::array<Real, LIBMESH_DIM> width = {
267 156 : {_xmax - _xmin, _dim > 1 ? _ymax - _ymin : 0, _dim > 2 ? _zmax - _zmin : 0}};
268 :
269 : // Min mesh extent in each direction.
270 156 : std::array<Real, LIBMESH_DIM> mins = {{_xmin, _dim > 1 ? _ymin : MIN, _dim > 2 ? _zmin : MIN}};
271 :
272 : // Number of elements in each direction.
273 156 : std::array<unsigned int, LIBMESH_DIM> nelem = {{_nx, _dim > 1 ? _ny : 1, _dim > 2 ? _nz : 1}};
274 :
275 : // We will need the biases raised to integer powers in each
276 : // direction, so let's pre-compute those...
277 156 : std::array<std::vector<Real>, LIBMESH_DIM> pows;
278 624 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
279 : {
280 468 : pows[dir].resize(nelem[dir] + 1);
281 468 : pows[dir][0] = 1.0;
282 1952 : for (unsigned int i = 1; i < pows[dir].size(); ++i)
283 1484 : pows[dir][i] = pows[dir][i - 1] * bias[dir];
284 : }
285 :
286 : // Loop over the nodes and move them to the desired location
287 8492 : for (auto & node_ptr : mesh->node_ptr_range())
288 : {
289 8336 : Node & node = *node_ptr;
290 :
291 33344 : for (unsigned int dir = 0; dir < LIBMESH_DIM; ++dir)
292 : {
293 25008 : if (width[dir] != 0. && bias[dir] != 1.)
294 : {
295 : // Compute the scaled "index" of the current point. This
296 : // will either be close to a whole integer or a whole
297 : // integer+0.5 for quadratic nodes.
298 14696 : Real float_index = (node(dir) - mins[dir]) * nelem[dir] / width[dir];
299 :
300 14696 : Real integer_part = 0;
301 14696 : Real fractional_part = std::modf(float_index, &integer_part);
302 :
303 : // Figure out where to move the node...
304 14696 : if (std::abs(fractional_part) < TOLERANCE || std::abs(fractional_part - 1.0) < TOLERANCE)
305 : {
306 : // If the fractional_part ~ 0.0 or 1.0, this is a vertex node, so
307 : // we don't need an average.
308 : //
309 : // Compute the "index" we are at in the current direction. We
310 : // round to the nearest integral value to get this instead
311 : // of using "integer_part", since that could be off by a
312 : // lot (e.g. we want 3.9999 to map to 4.0 instead of 3.0).
313 14696 : int index = round(float_index);
314 :
315 : mooseAssert(index >= static_cast<int>(0) && index < static_cast<int>(pows[dir].size()),
316 : "Scaled \"index\" out of range");
317 :
318 : // Move node to biased location.
319 14696 : node(dir) =
320 14696 : mins[dir] + width[dir] * (1. - pows[dir][index]) / (1. - pows[dir][nelem[dir]]);
321 : }
322 0 : else if (std::abs(fractional_part - 0.5) < TOLERANCE)
323 : {
324 : // If the fractional_part ~ 0.5, this is a midedge/face
325 : // (i.e. quadratic) node. We don't move those with the same
326 : // bias as the vertices, instead we put them midway between
327 : // their respective vertices.
328 : //
329 : // Also, since the fractional part is nearly 0.5, we know that
330 : // the integer_part will be the index of the vertex to the
331 : // left, and integer_part+1 will be the index of the
332 : // vertex to the right.
333 0 : node(dir) = mins[dir] +
334 0 : width[dir] *
335 0 : (1. - 0.5 * (pows[dir][integer_part] + pows[dir][integer_part + 1])) /
336 0 : (1. - pows[dir][nelem[dir]]);
337 : }
338 : else
339 : {
340 : // We don't yet handle anything higher order than quadratic...
341 0 : mooseError("Unable to bias node at node(", dir, ")=", node(dir));
342 : }
343 : }
344 : }
345 156 : }
346 156 : }
347 :
348 19072 : mesh->unset_is_prepared();
349 38144 : return dynamic_pointer_cast<MeshBase>(mesh);
350 19072 : }
|
