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 "ExamplePatchMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 :
13 : // libMesh includes
14 : #include "libmesh/mesh_generation.h"
15 : #include "libmesh/unstructured_mesh.h"
16 : #include "libmesh/point.h"
17 : #include "libmesh/elem.h"
18 : #include "libmesh/node.h"
19 : #include "libmesh/boundary_info.h"
20 : #include "libmesh/face_quad4.h"
21 : #include "libmesh/face_quad8.h"
22 : #include "libmesh/cell_hex8.h"
23 : #include "libmesh/cell_hex20.h"
24 : #include "libmesh/utility.h"
25 :
26 : registerMooseObject("MooseApp", ExamplePatchMeshGenerator);
27 : registerMooseObjectRenamed("MooseApp",
28 : PatchMeshGenerator,
29 : "06/30/2025 24:00",
30 : ExamplePatchMeshGenerator);
31 :
32 : InputParameters
33 28618 : ExamplePatchMeshGenerator::validParams()
34 : {
35 28618 : InputParameters params = MeshGenerator::validParams();
36 28618 : MooseEnum dims("2=2 3", "2");
37 28618 : params.addParam<MooseEnum>("dim",
38 : dims,
39 : "The dimension of the mesh to be generated. Patch meshes are only "
40 : "valid in 2 or 3 dimensions.");
41 :
42 28618 : MooseEnum elem_types("QUAD4 QUAD8 HEX8 HEX20", "QUAD4");
43 28618 : params.addParam<MooseEnum>("elem_type",
44 : elem_types,
45 : "The type of element from libMesh to "
46 : "generate (default: linear element for "
47 : "requested dimension)");
48 28618 : params.addParam<Real>("x_length", 0.24, "Length of the domain in the x direction.");
49 28618 : params.addParam<Real>("y_length", 0.12, "Length of the domain in the y direction.");
50 28618 : params.addParam<Real>("z_length", 0.0, "Length of the domain in the z direction.");
51 28618 : params.addParam<Real>("x_offset", 0.0, "Offset of the Cartesian origin in the x direction.");
52 28618 : params.addParam<Real>("y_offset", 0.0, "Offset of the Cartesian origin in the y direction.");
53 28618 : params.addParam<Real>("z_offset", 0.0, "Offset of the Cartesian origin in the z direction.");
54 28618 : params.addClassDescription("Creates 2D or 3D patch meshes.");
55 57236 : return params;
56 28618 : }
57 :
58 44 : ExamplePatchMeshGenerator::ExamplePatchMeshGenerator(const InputParameters & parameters)
59 : : MeshGenerator(parameters),
60 44 : _dim(getParam<MooseEnum>("dim")),
61 44 : _elem_type(getParam<MooseEnum>("elem_type")),
62 44 : _xlength(getParam<Real>("x_length")),
63 44 : _ylength(getParam<Real>("y_length")),
64 44 : _zlength(getParam<Real>("z_length")),
65 44 : _xoffset(getParam<Real>("x_offset")),
66 44 : _yoffset(getParam<Real>("y_offset")),
67 88 : _zoffset(getParam<Real>("z_offset"))
68 : {
69 44 : if (_xlength <= 0.0)
70 0 : paramError("x_length", "Must be greater than zero");
71 :
72 44 : if (_ylength <= 0.0)
73 0 : paramError("y_length", "Must be greater than zero");
74 :
75 44 : if (_dim == 2 && (_elem_type == "HEX8" || _elem_type == "HEX20"))
76 0 : paramError("elem_type", "Must be QUAD4 or QUAD8 for 2-dimensional meshes.");
77 :
78 44 : if (_dim == 2 && _zlength != 0.0)
79 0 : paramError("z_length", "Must be zero for 2-dimensional meshes");
80 :
81 44 : if (_dim == 3 && !parameters.isParamSetByUser("elem_type"))
82 0 : _elem_type = "HEX8";
83 :
84 44 : if (_dim == 3 && (_elem_type == "QUAD4" || _elem_type == "QUAD8"))
85 0 : paramError("elem_type", "Must be HEX8 or HEX20 for 3-dimensional meshes.");
86 :
87 : // If domain dimensions are not set by user for 3-dimensions a unit cube is used
88 : // as default as per MacNeal and Harder (1985)
89 44 : if (_dim == 3)
90 : {
91 22 : if (!parameters.isParamSetByUser("x_length"))
92 22 : _xlength = 1.0;
93 22 : if (!parameters.isParamSetByUser("y_length"))
94 22 : _ylength = 1.0;
95 22 : if (!parameters.isParamSetByUser("z_length"))
96 22 : _zlength = 1.0;
97 : }
98 :
99 44 : if (_dim == 3 && _zlength <= 0.0)
100 0 : paramError("z_length", "Must be greater than zero for 3-dimensional meshes");
101 44 : }
102 :
103 : std::unique_ptr<MeshBase>
104 44 : ExamplePatchMeshGenerator::generate()
105 : {
106 44 : auto mesh = buildMeshBaseObject();
107 44 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
108 :
109 44 : unsigned num_nodes = 0;
110 44 : if (_elem_type == "QUAD4")
111 11 : num_nodes = 8;
112 33 : else if (_elem_type == "QUAD8")
113 11 : num_nodes = 20;
114 22 : else if (_elem_type == "HEX8")
115 11 : num_nodes = 16;
116 11 : else if (_elem_type == "HEX20")
117 11 : num_nodes = 48;
118 :
119 44 : std::vector<Node *> nodes(num_nodes);
120 :
121 44 : const Real xmax = _xoffset + _xlength;
122 44 : const Real ymax = _yoffset + _ylength;
123 44 : const Real zmax = _zoffset + _zlength;
124 :
125 44 : if (_dim == 2)
126 : {
127 22 : const std::vector<Real> x_interior_fractions{1.0 / 6.0, 3.0 / 4.0, 2.0 / 3.0, 1.0 / 3.0};
128 22 : const std::vector<Real> y_interior_fractions{1.0 / 6.0, 1.0 / 4.0, 2.0 / 3.0, 2.0 / 3.0};
129 :
130 : // Exterior node positions
131 22 : std::vector<Point> node_positions{{_xoffset, _yoffset, _zoffset},
132 22 : {xmax, _yoffset, _zoffset},
133 22 : {xmax, ymax, _zoffset},
134 22 : {_xoffset, ymax, _zoffset}};
135 : // Interior node positions
136 110 : for (unsigned i = 0; i < x_interior_fractions.size(); ++i)
137 176 : node_positions.push_back({_xoffset + x_interior_fractions[i] * _xlength,
138 88 : _yoffset + y_interior_fractions[i] * _ylength,
139 88 : _zoffset});
140 :
141 22 : if (_elem_type == "QUAD8")
142 : {
143 : // Exterior midside node positions
144 55 : for (unsigned i = 0; i < 4; ++i)
145 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4]));
146 :
147 : // Exterior to interior midside node positions
148 55 : for (unsigned i = 0; i < 4; ++i)
149 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
150 :
151 : // Interior midside node positions
152 55 : for (unsigned i = 4; i < 8; ++i)
153 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 4]));
154 : }
155 :
156 330 : for (unsigned node_id = 0; node_id < num_nodes; ++node_id)
157 308 : nodes[node_id] = mesh->add_point(node_positions[node_id], node_id);
158 :
159 22 : if (_elem_type == "QUAD4")
160 11 : makeQuad4Elems(*mesh, nodes);
161 : else
162 11 : makeQuad8Elems(*mesh, nodes);
163 :
164 22 : boundary_info.sideset_name(1) = "bottom";
165 22 : boundary_info.sideset_name(2) = "right";
166 22 : boundary_info.sideset_name(3) = "top";
167 22 : boundary_info.sideset_name(4) = "left";
168 :
169 22 : boundary_info.nodeset_name(100) = "bottom_left";
170 22 : boundary_info.nodeset_name(101) = "bottom_right";
171 22 : boundary_info.nodeset_name(102) = "top_right";
172 22 : boundary_info.nodeset_name(103) = "top_left";
173 22 : }
174 : else
175 : {
176 : const std::vector<Real> x_interior_fractions{
177 22 : 0.249, 0.826, 0.850, 0.273, 0.320, 0.677, 0.788, 0.165};
178 : const std::vector<Real> y_interior_fractions{
179 22 : 0.342, 0.288, 0.649, 0.750, 0.186, 0.305, 0.693, 0.745};
180 : const std::vector<Real> z_interior_fractions{
181 22 : 0.192, 0.288, 0.263, 0.230, 0.643, 0.683, 0.644, 0.702};
182 :
183 : // Exterior node positions
184 22 : std::vector<Point> node_positions{{_xoffset, _yoffset, _zoffset},
185 22 : {xmax, _yoffset, _zoffset},
186 22 : {xmax, ymax, _zoffset},
187 22 : {_xoffset, ymax, _zoffset},
188 22 : {_xoffset, _yoffset, zmax},
189 22 : {xmax, _yoffset, zmax},
190 : {xmax, ymax, zmax},
191 22 : {_xoffset, ymax, zmax}};
192 :
193 : // Interior node positions
194 198 : for (unsigned i = 0; i < x_interior_fractions.size(); ++i)
195 352 : node_positions.push_back({_xoffset + x_interior_fractions[i] * _xlength,
196 176 : _yoffset + y_interior_fractions[i] * _ylength,
197 176 : _zoffset + z_interior_fractions[i] * _zlength});
198 :
199 22 : if (_elem_type == "HEX20")
200 : {
201 : // Midside Nodes
202 :
203 : // Four on back face
204 55 : for (unsigned i = 0; i < 4; ++i)
205 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4]));
206 :
207 : // Four on front face
208 55 : for (unsigned i = 4; i < 8; ++i)
209 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 4]));
210 :
211 : // Four on remaining exterior edges
212 55 : for (unsigned i = 0; i < 4; ++i)
213 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
214 :
215 : // Four on interior hex back face
216 55 : for (unsigned i = 8; i < 12; ++i)
217 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 8]));
218 :
219 : // Four on interior hex front face
220 55 : for (unsigned i = 12; i < 16; ++i)
221 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 12]));
222 :
223 : // Four on remaining interior hex edges
224 55 : for (unsigned i = 8; i < 12; ++i)
225 44 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
226 :
227 : // Eight on remaining interior edges
228 99 : for (unsigned i = 0; i < 8; ++i)
229 88 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 8]));
230 : }
231 :
232 726 : for (unsigned node_id = 0; node_id < num_nodes; ++node_id)
233 704 : nodes[node_id] = mesh->add_point(node_positions[node_id], node_id);
234 :
235 22 : if (_elem_type == "HEX8")
236 11 : makeHex8Elems(*mesh, nodes);
237 : else
238 11 : makeHex20Elems(*mesh, nodes);
239 :
240 22 : boundary_info.sideset_name(1) = "front";
241 22 : boundary_info.sideset_name(2) = "bottom";
242 22 : boundary_info.sideset_name(3) = "left";
243 22 : boundary_info.sideset_name(4) = "right";
244 22 : boundary_info.sideset_name(5) = "top";
245 22 : boundary_info.sideset_name(6) = "back";
246 :
247 22 : boundary_info.nodeset_name(100) = "bottom_back_left";
248 22 : boundary_info.nodeset_name(101) = "bottom_back_right";
249 22 : boundary_info.nodeset_name(102) = "top_back_right";
250 22 : boundary_info.nodeset_name(103) = "top_back_left";
251 22 : boundary_info.nodeset_name(104) = "bottom_front_left";
252 22 : boundary_info.nodeset_name(105) = "bottom_front_right";
253 22 : boundary_info.nodeset_name(106) = "top_front_right";
254 22 : boundary_info.nodeset_name(107) = "top_front_left";
255 22 : }
256 :
257 44 : mesh->prepare_for_use();
258 :
259 88 : return dynamic_pointer_cast<MeshBase>(mesh);
260 44 : }
261 :
262 : void
263 11 : ExamplePatchMeshGenerator::makeQuad4Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
264 : {
265 11 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
266 : std::vector<std::vector<int>> element_connectivity{
267 66 : {0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}};
268 66 : for (unsigned i = 0; i < 5; ++i)
269 : {
270 55 : Elem * elem = mesh.add_elem(new Quad4);
271 275 : for (unsigned j = 0; j < 4; ++j)
272 220 : elem->set_node(j, nodes[element_connectivity[i][j]]);
273 :
274 55 : elem->subdomain_id() = i + 1;
275 :
276 55 : if (i < 4)
277 : {
278 44 : boundary_info.add_node(i, i + 100);
279 44 : boundary_info.add_side(elem, 0, i + 1);
280 : }
281 : }
282 33 : }
283 :
284 : void
285 11 : ExamplePatchMeshGenerator::makeQuad8Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
286 : {
287 11 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
288 : std::vector<std::vector<int>> element_connectivity{{0, 1, 5, 4, 8, 13, 16, 12},
289 : {1, 2, 6, 5, 9, 14, 17, 13},
290 : {2, 3, 7, 6, 10, 15, 18, 14},
291 : {3, 0, 4, 7, 11, 12, 19, 15},
292 66 : {4, 5, 6, 7, 16, 17, 18, 19}};
293 66 : for (unsigned i = 0; i < 5; ++i)
294 : {
295 55 : Elem * elem = mesh.add_elem(new Quad8);
296 495 : for (unsigned j = 0; j < 8; ++j)
297 440 : elem->set_node(j, nodes[element_connectivity[i][j]]);
298 :
299 55 : elem->subdomain_id() = i + 1;
300 :
301 55 : if (i < 4)
302 : {
303 44 : boundary_info.add_node(i, i + 100);
304 44 : boundary_info.add_side(elem, 0, i + 1);
305 : }
306 : }
307 33 : }
308 :
309 : void
310 11 : ExamplePatchMeshGenerator::makeHex8Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
311 : {
312 11 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
313 : std::vector<std::vector<int>> element_connectivity{{12, 13, 14, 15, 4, 5, 6, 7},
314 : {8, 9, 13, 12, 0, 1, 5, 4},
315 : {8, 12, 15, 11, 0, 4, 7, 3},
316 : {9, 10, 14, 13, 1, 2, 6, 5},
317 : {10, 11, 15, 14, 2, 3, 7, 6},
318 : {8, 11, 10, 9, 0, 3, 2, 1},
319 88 : {8, 9, 10, 11, 12, 13, 14, 15}};
320 :
321 88 : for (unsigned i = 0; i < 7; ++i)
322 : {
323 77 : Elem * elem = mesh.add_elem(new Hex8);
324 693 : for (unsigned j = 0; j < 8; ++j)
325 616 : elem->set_node(j, nodes[element_connectivity[i][j]]);
326 :
327 77 : elem->subdomain_id() = i + 1;
328 77 : if (i < 6)
329 66 : boundary_info.add_side(elem, 5, i + 1);
330 : }
331 99 : for (unsigned i = 0; i < 8; ++i)
332 88 : boundary_info.add_node(i, i + 100);
333 33 : }
334 :
335 : void
336 11 : ExamplePatchMeshGenerator::makeHex20Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
337 : {
338 11 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
339 : std::vector<std::vector<int>> element_connectivity{
340 : {12, 13, 14, 15, 4, 5, 6, 7, 32, 33, 34, 35, 44, 45, 46, 47, 20, 21, 22, 23},
341 : {8, 9, 13, 12, 0, 1, 5, 4, 28, 37, 32, 36, 40, 41, 45, 44, 16, 25, 20, 24},
342 : {8, 12, 15, 11, 0, 4, 7, 3, 36, 35, 39, 31, 40, 44, 47, 43, 24, 23, 27, 19},
343 : {9, 10, 14, 13, 1, 2, 6, 5, 29, 38, 33, 37, 41, 42, 46, 45, 17, 26, 21, 25},
344 : {10, 11, 15, 14, 2, 3, 7, 6, 30, 39, 34, 38, 42, 43, 47, 46, 18, 27, 22, 26},
345 : {8, 11, 10, 9, 0, 3, 2, 1, 31, 30, 29, 28, 40, 43, 42, 41, 19, 18, 17, 16},
346 88 : {8, 9, 10, 11, 12, 13, 14, 15, 28, 29, 30, 31, 36, 37, 38, 39, 32, 33, 34, 35}};
347 :
348 88 : for (unsigned i = 0; i < 7; ++i)
349 : {
350 77 : Elem * elem = mesh.add_elem(new Hex20);
351 1617 : for (unsigned j = 0; j < 20; ++j)
352 1540 : elem->set_node(j, nodes[element_connectivity[i][j]]);
353 :
354 77 : elem->subdomain_id() = i + 1;
355 77 : if (i < 6)
356 66 : boundary_info.add_side(elem, 5, i + 1);
357 : }
358 99 : for (unsigned i = 0; i < 8; ++i)
359 88 : boundary_info.add_node(i, i + 100);
360 33 : }
|
