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 "CartesianMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 :
13 : // libMesh includes
14 : #include "libmesh/mesh_generation.h"
15 : #include "libmesh/unstructured_mesh.h"
16 : #include "libmesh/replicated_mesh.h"
17 : #include "libmesh/point.h"
18 : #include "libmesh/elem.h"
19 : #include "libmesh/node.h"
20 :
21 : registerMooseObject("MooseApp", CartesianMeshGenerator);
22 :
23 : InputParameters
24 9195 : CartesianMeshGenerator::validParams()
25 : {
26 9195 : InputParameters params = MeshGenerator::validParams();
27 36780 : MooseEnum dims("1=1 2 3");
28 36780 : params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
29 :
30 36780 : params.addRequiredParam<std::vector<Real>>("dx", "Intervals in the X direction");
31 36780 : params.addParam<std::vector<unsigned int>>(
32 : "ix", "Number of grids in all intervals in the X direction (default to all one)");
33 36780 : params.addParam<std::vector<Real>>(
34 : "max_sizes_in_x",
35 : "Maximum element size in the X direction. This can be specified once for all 'dx' intervals "
36 : "or specified for each interval. 'ix' is then determined using the ceiling of dx/size");
37 36780 : params.addParam<std::vector<Real>>(
38 : "dy", "Intervals in the Y direction (required when dim>1 otherwise ignored)");
39 36780 : params.addParam<std::vector<unsigned int>>(
40 : "iy", "Number of grids in all intervals in the Y direction (default to all one)");
41 36780 : params.addParam<std::vector<Real>>(
42 : "max_sizes_in_y",
43 : "Maximum element size in the Y direction. This can be specified once for all 'dy' intervals "
44 : "or specified for each interval. 'iy' is then determined using the ceiling of dy/size");
45 36780 : params.addParam<std::vector<Real>>(
46 : "dz", "Intervals in the Z direction (required when dim>2 otherwise ignored)");
47 36780 : params.addParam<std::vector<unsigned int>>(
48 : "iz", "Number of grids in all intervals in the Z direction (default to all one)");
49 36780 : params.addParam<std::vector<Real>>(
50 : "max_sizes_in_z",
51 : "Maximum element size in the Z direction. This can be specified once for all 'dz' intervals "
52 : "or specified for each interval. 'iz' is then determined using the ceiling of dz/size");
53 36780 : params.addParam<std::vector<unsigned int>>("subdomain_id", "Block IDs (default to all zero)");
54 9195 : params.addClassDescription("This CartesianMeshGenerator creates a non-uniform Cartesian mesh.");
55 18390 : return params;
56 9195 : }
57 :
58 3061 : CartesianMeshGenerator::CartesianMeshGenerator(const InputParameters & parameters)
59 : : MeshGenerator(parameters),
60 3061 : _dim(getParam<MooseEnum>("dim")),
61 12244 : _dx(getParam<std::vector<Real>>("dx"))
62 : {
63 : // get all other parameters if provided and check their sizes
64 9183 : if (isParamValid("ix"))
65 : {
66 4382 : _ix = getParam<std::vector<unsigned int>>("ix");
67 2191 : if (_ix.size() != _dx.size())
68 0 : mooseError("ix must be the same size as dx");
69 7392 : for (unsigned int i = 0; i < _ix.size(); ++i)
70 5201 : if (_ix[i] == 0)
71 0 : mooseError("ix cannot be zero");
72 6573 : if (isParamValid("max_sizes_in_x"))
73 0 : paramError("max_sizes_in_x", "Cannot be provided if 'ix' is specified");
74 : }
75 2610 : else if (isParamValid("max_sizes_in_x"))
76 : {
77 8 : _ix.resize(_dx.size());
78 16 : auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_x");
79 8 : if (max_sizes.size() == 1)
80 0 : max_sizes.resize(_dx.size(), max_sizes[0]);
81 8 : if (max_sizes.size() != _dx.size())
82 0 : paramError("max_sizes_in_x", "max_sizes_in_x must be size 1 or the same size as dx");
83 32 : for (const auto i : index_range(_dx))
84 24 : if (max_sizes[i] <= 0)
85 0 : paramError("max_sizes_in_x", "max_sizes_in_x must be strictly positive.");
86 32 : for (const auto i : index_range(_dx))
87 24 : _ix[i] = std::ceil(_dx[i] / max_sizes[i]);
88 8 : }
89 : else
90 862 : _ix = std::vector<unsigned int>(_dx.size(), 1);
91 :
92 9604 : for (unsigned int i = 0; i < _dx.size(); ++i)
93 6543 : if (_dx[i] <= 0)
94 0 : mooseError("dx must be greater than zero. dx: ", Moose::stringify(_dx));
95 :
96 9183 : if (isParamValid("dy"))
97 : {
98 4050 : _dy = getParam<std::vector<Real>>("dy");
99 6204 : for (unsigned int i = 0; i < _dy.size(); ++i)
100 4179 : if (_dy[i] <= 0)
101 0 : mooseError("dy must be greater than zero. dy: ", Moose::stringify(_dy));
102 : }
103 :
104 9183 : if (isParamValid("iy"))
105 : {
106 3746 : _iy = getParam<std::vector<unsigned int>>("iy");
107 1873 : if (_iy.size() != _dy.size())
108 0 : mooseError("iy must be the same size as dy");
109 5576 : for (unsigned int i = 0; i < _iy.size(); ++i)
110 3703 : if (_iy[i] == 0)
111 0 : mooseError("iy cannot be zero");
112 5619 : if (isParamValid("max_sizes_in_y"))
113 0 : paramError("max_sizes_in_y", "Cannot be provided if 'iy' is specified");
114 : }
115 3564 : else if (isParamValid("max_sizes_in_y"))
116 : {
117 8 : _iy.resize(_dy.size());
118 16 : auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_y");
119 8 : if (max_sizes.size() == 1)
120 0 : max_sizes.resize(_dy.size(), max_sizes[0]);
121 8 : if (max_sizes.size() != _dy.size())
122 0 : paramError("max_sizes_in_y", "max_sizes_in_y must be size 1 or the same size as dy");
123 24 : for (const auto i : index_range(_dy))
124 16 : if (max_sizes[i] <= 0)
125 0 : paramError("max_sizes_in_y", "max_sizes_in_y must be strictly positive.");
126 24 : for (const auto i : index_range(_dy))
127 16 : _iy[i] = std::ceil(_dy[i] / max_sizes[i]);
128 8 : }
129 : else
130 1180 : _iy = std::vector<unsigned int>(_dy.size(), 1);
131 :
132 9183 : if (isParamValid("dz"))
133 : {
134 936 : _dz = getParam<std::vector<Real>>("dz");
135 1902 : for (unsigned int i = 0; i < _dz.size(); ++i)
136 1434 : if (_dz[i] <= 0)
137 0 : mooseError("dz must be greater than zero dz: ", Moose::stringify(_dz));
138 : }
139 :
140 9183 : if (isParamValid("iz"))
141 : {
142 814 : _iz = getParam<std::vector<unsigned int>>("iz");
143 407 : if (_iz.size() != _dz.size())
144 0 : mooseError("iz must be the same size as dz");
145 1597 : for (unsigned int i = 0; i < _iz.size(); ++i)
146 1190 : if (_iz[i] == 0)
147 0 : mooseError("iz cannot be zero");
148 1221 : if (isParamValid("max_sizes_in_z"))
149 0 : paramError("max_sizes_in_z", "Cannot be provided if 'iz' is specified");
150 : }
151 7962 : else if (isParamValid("max_sizes_in_z"))
152 : {
153 8 : _iz.resize(_dz.size());
154 16 : auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_z");
155 8 : if (max_sizes.size() == 1)
156 0 : max_sizes.resize(_dz.size(), max_sizes[0]);
157 8 : if (max_sizes.size() != _dz.size())
158 0 : paramError("max_sizes_in_z", "max_sizes_in_z must be size 1 or the same size as dz");
159 40 : for (const auto i : index_range(_dz))
160 32 : if (max_sizes[i] <= 0)
161 0 : paramError("max_sizes_in_z", "max_sizes_in_z must be strictly positive.");
162 40 : for (const auto i : index_range(_dz))
163 32 : _iz[i] = std::ceil(_dz[i] / max_sizes[i]);
164 8 : }
165 : else
166 2646 : _iz = std::vector<unsigned int>(_dz.size(), 1);
167 :
168 9183 : if (isParamValid("subdomain_id"))
169 : {
170 3476 : _subdomain_id = getParam<std::vector<unsigned int>>("subdomain_id");
171 5780 : if (isParamValid("dz") && isParamValid("dy"))
172 : {
173 283 : if (_subdomain_id.size() != _dx.size() * _dy.size() * _dz.size())
174 0 : mooseError(
175 0 : "subdomain_id must be the size of the product of sizes of dx, dy and dz. Sizes are '" +
176 0 : std::to_string(_subdomain_id.size()) + "' and '" +
177 0 : std::to_string(_dx.size() * _dy.size() * _dz.size()) + "' respectively");
178 : }
179 4365 : else if (isParamValid("dy"))
180 : {
181 1359 : if (_subdomain_id.size() != _dx.size() * _dy.size())
182 0 : mooseError(
183 0 : "subdomain_id must be the size of the product of sizes of dx and dy. Sizes are '" +
184 0 : std::to_string(_subdomain_id.size()) + "' and '" +
185 0 : std::to_string(_dx.size() * _dy.size()) + "' respectively");
186 : }
187 : else
188 : {
189 96 : if (_subdomain_id.size() != _dx.size())
190 0 : mooseError("subdomain_id must be the size of the product of sizes of dx. Sizes are '" +
191 0 : std::to_string(_subdomain_id.size()) + "' and '" + std::to_string(_dx.size()) +
192 : "' respectively");
193 : }
194 : }
195 : else
196 : {
197 3969 : if (isParamValid("dz"))
198 185 : _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size() * _dz.size(), 0);
199 3414 : else if (isParamValid("dy"))
200 198 : _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size(), 0);
201 : else
202 940 : _subdomain_id = std::vector<unsigned int>(_dx.size(), 0);
203 : }
204 :
205 : // do dimension checks and expand block IDs for all sub-grids
206 3061 : switch (_dim)
207 : {
208 1036 : case 1:
209 : {
210 1036 : _nx = 0;
211 2380 : for (unsigned int i = 0; i < _dx.size(); ++i)
212 1344 : _nx += _ix[i];
213 1036 : _ny = 1;
214 1036 : _nz = 1;
215 :
216 1036 : std::vector<unsigned int> new_id;
217 2380 : for (unsigned int i = 0; i < _dx.size(); ++i)
218 5803 : for (unsigned int ii = 0; ii < _ix[i]; ++ii)
219 4459 : new_id.push_back(_subdomain_id[i]);
220 :
221 1036 : _subdomain_id = new_id;
222 :
223 3108 : if (isParamValid("dy"))
224 0 : mooseWarning("dy provided for 1D");
225 3108 : if (isParamValid("iy"))
226 0 : mooseWarning("iy provided for 1D");
227 3108 : if (isParamValid("dz"))
228 0 : mooseWarning("dz provided for 1D");
229 3108 : if (isParamValid("iz"))
230 0 : mooseWarning("iz provided for 1D");
231 3108 : if (isParamValid("max_sizes_in_y"))
232 0 : paramWarning("max_sizes_in_y", "Should not be provided for 1D");
233 3108 : if (isParamValid("max_sizes_in_z"))
234 0 : paramWarning("max_sizes_in_z", "Should not be provided for 1D");
235 1036 : break;
236 1036 : }
237 1557 : case 2:
238 : {
239 1557 : _nx = 0;
240 5471 : for (unsigned int i = 0; i < _dx.size(); ++i)
241 3914 : _nx += _ix[i];
242 1557 : _ny = 0;
243 4780 : for (unsigned int i = 0; i < _dy.size(); ++i)
244 3223 : _ny += _iy[i];
245 1557 : _nz = 1;
246 :
247 1557 : std::vector<unsigned int> new_id;
248 4780 : for (unsigned int j = 0; j < _dy.size(); ++j)
249 11773 : for (unsigned int jj = 0; jj < _iy[j]; ++jj)
250 32710 : for (unsigned int i = 0; i < _dx.size(); ++i)
251 92607 : for (unsigned int ii = 0; ii < _ix[i]; ++ii)
252 68447 : new_id.push_back(_subdomain_id[j * _dx.size() + i]);
253 :
254 1557 : _subdomain_id = new_id;
255 :
256 4671 : if (!isParamValid("dy"))
257 0 : mooseError("dy is not provided for 2D");
258 4671 : if (isParamValid("dz"))
259 0 : mooseWarning("dz provided for 2D");
260 4671 : if (isParamValid("iz"))
261 0 : mooseWarning("iz provided for 2D");
262 4671 : if (isParamValid("max_sizes_in_z"))
263 0 : paramWarning("max_sizes_in_z", "Should not be provided for 2D");
264 1557 : break;
265 1557 : }
266 468 : case 3:
267 : {
268 468 : _nx = 0;
269 1753 : for (unsigned int i = 0; i < _dx.size(); ++i)
270 1285 : _nx += _ix[i];
271 468 : _ny = 0;
272 1424 : for (unsigned int i = 0; i < _dy.size(); ++i)
273 956 : _ny += _iy[i];
274 468 : _nz = 0;
275 1902 : for (unsigned int i = 0; i < _dz.size(); ++i)
276 1434 : _nz += _iz[i];
277 :
278 468 : std::vector<unsigned int> new_id;
279 1902 : for (unsigned int k = 0; k < _dz.size(); ++k)
280 3392 : for (unsigned int kk = 0; kk < _iz[k]; ++kk)
281 5958 : for (unsigned int j = 0; j < _dy.size(); ++j)
282 12870 : for (unsigned int jj = 0; jj < _iy[j]; ++jj)
283 34344 : for (unsigned int i = 0; i < _dx.size(); ++i)
284 59606 : for (unsigned int ii = 0; ii < _ix[i]; ++ii)
285 34132 : new_id.push_back(_subdomain_id[k * _dx.size() * _dy.size() + j * _dx.size() + i]);
286 :
287 468 : _subdomain_id = new_id;
288 :
289 1404 : if (!isParamValid("dy"))
290 0 : mooseError("dy is not provided for 3D");
291 1404 : if (!isParamValid("dz"))
292 0 : mooseError("dz is not provided for 3D");
293 468 : break;
294 468 : }
295 : }
296 3061 : }
297 :
298 : std::unique_ptr<MeshBase>
299 2884 : CartesianMeshGenerator::generate()
300 : {
301 2884 : auto mesh = buildMeshBaseObject();
302 :
303 : // switching on MooseEnum to generate the reference mesh
304 : // Note: element type is fixed
305 2884 : switch (_dim)
306 : {
307 1003 : case 1:
308 1003 : MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh), _nx, 0, _nx, EDGE2);
309 1003 : break;
310 1439 : case 2:
311 1439 : MeshTools::Generation::build_square(
312 1439 : static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, 0, _nx, 0, _ny, QUAD4);
313 1439 : break;
314 442 : case 3:
315 442 : MeshTools::Generation::build_cube(
316 442 : static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, _nz, 0, _nx, 0, _ny, 0, _nz, HEX8);
317 442 : break;
318 : }
319 :
320 : // assign block IDs
321 5768 : MeshBase::element_iterator el = mesh->active_elements_begin();
322 5768 : MeshBase::element_iterator el_end = mesh->active_elements_end();
323 99975 : for (; el != el_end; ++el)
324 : {
325 97091 : const Point p = (*el)->vertex_average();
326 97091 : unsigned int ix = std::floor(p(0));
327 97091 : unsigned int iy = std::floor(p(1));
328 97091 : unsigned int iz = std::floor(p(2));
329 97091 : unsigned int i = iz * _nx * _ny + iy * _nx + ix;
330 97091 : (*el)->subdomain_id() = _subdomain_id[i];
331 : }
332 :
333 : // adjust node coordinates
334 2884 : switch (_dim)
335 : {
336 1003 : case 1:
337 : {
338 : Real base;
339 :
340 1003 : std::vector<Real> mapx;
341 : // Note: the starting coordinate is zero
342 1003 : base = 0;
343 1003 : mapx.push_back(base);
344 2303 : for (unsigned int i = 0; i < _dx.size(); ++i)
345 : {
346 5561 : for (unsigned int j = 1; j <= _ix[i]; ++j)
347 4261 : mapx.push_back(base + _dx[i] / _ix[i] * j);
348 1300 : base += _dx[i];
349 : }
350 :
351 1003 : MeshBase::node_iterator node = mesh->active_nodes_begin();
352 1003 : MeshBase::node_iterator node_end = mesh->active_nodes_end();
353 6063 : for (; node != node_end; ++node)
354 : {
355 5060 : unsigned int i = (*(*node))(0) + 0.5;
356 5060 : (*(*node))(0) = mapx.at(i);
357 : }
358 1003 : break;
359 1003 : }
360 1439 : case 2:
361 : {
362 : Real base;
363 :
364 1439 : std::vector<Real> mapx;
365 1439 : base = 0;
366 1439 : mapx.push_back(base);
367 5078 : for (unsigned int i = 0; i < _dx.size(); ++i)
368 : {
369 11861 : for (unsigned int j = 1; j <= _ix[i]; ++j)
370 8222 : mapx.push_back(base + _dx[i] / _ix[i] * j);
371 3639 : base += _dx[i];
372 : }
373 :
374 1439 : std::vector<Real> mapy;
375 1439 : base = 0;
376 1439 : mapy.push_back(base);
377 4418 : for (unsigned int i = 0; i < _dy.size(); ++i)
378 : {
379 10864 : for (unsigned int j = 1; j <= _iy[i]; ++j)
380 7885 : mapy.push_back(base + _dy[i] / _iy[i] * j);
381 2979 : base += _dy[i];
382 : }
383 :
384 1439 : MeshBase::node_iterator node = mesh->active_nodes_begin();
385 1439 : MeshBase::node_iterator node_end = mesh->active_nodes_end();
386 80090 : for (; node != node_end; ++node)
387 : {
388 78651 : unsigned int i = (*(*node))(0) + 0.5;
389 78651 : (*(*node))(0) = mapx.at(i);
390 78651 : unsigned int j = (*(*node))(1) + 0.5;
391 78651 : (*(*node))(1) = mapy.at(j);
392 : }
393 1439 : break;
394 1439 : }
395 442 : case 3:
396 : {
397 : Real base;
398 :
399 442 : std::vector<Real> mapx;
400 442 : base = 0;
401 442 : mapx.push_back(base);
402 1667 : for (unsigned int i = 0; i < _dx.size(); ++i)
403 : {
404 2889 : for (unsigned int j = 1; j <= _ix[i]; ++j)
405 1664 : mapx.push_back(base + _dx[i] / _ix[i] * j);
406 1225 : base += _dx[i];
407 : }
408 :
409 442 : std::vector<Real> mapy;
410 442 : base = 0;
411 442 : mapy.push_back(base);
412 1355 : for (unsigned int i = 0; i < _dy.size(); ++i)
413 : {
414 2889 : for (unsigned int j = 1; j <= _iy[i]; ++j)
415 1976 : mapy.push_back(base + _dy[i] / _iy[i] * j);
416 913 : base += _dy[i];
417 : }
418 :
419 442 : std::vector<Real> mapz;
420 442 : base = 0;
421 442 : mapz.push_back(base);
422 1799 : for (unsigned int i = 0; i < _dz.size(); ++i)
423 : {
424 3210 : for (unsigned int j = 1; j <= _iz[i]; ++j)
425 1853 : mapz.push_back(base + _dz[i] / _iz[i] * j);
426 1357 : base += _dz[i];
427 : }
428 :
429 442 : MeshBase::node_iterator node = mesh->active_nodes_begin();
430 442 : MeshBase::node_iterator node_end = mesh->active_nodes_end();
431 60702 : for (; node != node_end; ++node)
432 : {
433 60260 : unsigned int i = (*(*node))(0) + 0.5;
434 60260 : (*(*node))(0) = mapx.at(i);
435 60260 : unsigned int j = (*(*node))(1) + 0.5;
436 60260 : (*(*node))(1) = mapy.at(j);
437 60260 : unsigned int k = (*(*node))(2) + 0.5;
438 60260 : (*(*node))(2) = mapz.at(k);
439 : }
440 442 : break;
441 442 : }
442 : }
443 :
444 2884 : mesh->prepare_for_use();
445 5768 : return dynamic_pointer_cast<MeshBase>(mesh);
446 2884 : }
|
