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 28610 : ExamplePatchMeshGenerator::validParams()
34 : {
35 28610 : InputParameters params = MeshGenerator::validParams();
36 28610 : MooseEnum dims("2=2 3", "2");
37 28610 : 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 28610 : MooseEnum elem_types("QUAD4 QUAD8 HEX8 HEX20", "QUAD4");
43 28610 : 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 28610 : params.addParam<Real>("x_length", 0.24, "Length of the domain in the x direction.");
49 28610 : params.addParam<Real>("y_length", 0.12, "Length of the domain in the y direction.");
50 28610 : params.addParam<Real>("z_length", 0.0, "Length of the domain in the z direction.");
51 28610 : params.addParam<Real>("x_offset", 0.0, "Offset of the Cartesian origin in the x direction.");
52 28610 : params.addParam<Real>("y_offset", 0.0, "Offset of the Cartesian origin in the y direction.");
53 28610 : params.addParam<Real>("z_offset", 0.0, "Offset of the Cartesian origin in the z direction.");
54 28610 : params.addClassDescription("Creates 2D or 3D patch meshes.");
55 57220 : return params;
56 28610 : }
57 :
58 40 : ExamplePatchMeshGenerator::ExamplePatchMeshGenerator(const InputParameters & parameters)
59 : : MeshGenerator(parameters),
60 40 : _dim(getParam<MooseEnum>("dim")),
61 40 : _elem_type(getParam<MooseEnum>("elem_type")),
62 40 : _xlength(getParam<Real>("x_length")),
63 40 : _ylength(getParam<Real>("y_length")),
64 40 : _zlength(getParam<Real>("z_length")),
65 40 : _xoffset(getParam<Real>("x_offset")),
66 40 : _yoffset(getParam<Real>("y_offset")),
67 80 : _zoffset(getParam<Real>("z_offset"))
68 : {
69 40 : if (_xlength <= 0.0)
70 0 : paramError("x_length", "Must be greater than zero");
71 :
72 40 : if (_ylength <= 0.0)
73 0 : paramError("y_length", "Must be greater than zero");
74 :
75 40 : 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 40 : if (_dim == 2 && _zlength != 0.0)
79 0 : paramError("z_length", "Must be zero for 2-dimensional meshes");
80 :
81 40 : if (_dim == 3 && !parameters.isParamSetByUser("elem_type"))
82 0 : _elem_type = "HEX8";
83 :
84 40 : 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 40 : if (_dim == 3)
90 : {
91 20 : if (!parameters.isParamSetByUser("x_length"))
92 20 : _xlength = 1.0;
93 20 : if (!parameters.isParamSetByUser("y_length"))
94 20 : _ylength = 1.0;
95 20 : if (!parameters.isParamSetByUser("z_length"))
96 20 : _zlength = 1.0;
97 : }
98 :
99 40 : if (_dim == 3 && _zlength <= 0.0)
100 0 : paramError("z_length", "Must be greater than zero for 3-dimensional meshes");
101 40 : }
102 :
103 : std::unique_ptr<MeshBase>
104 40 : ExamplePatchMeshGenerator::generate()
105 : {
106 40 : auto mesh = buildMeshBaseObject();
107 40 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
108 :
109 40 : unsigned num_nodes = 0;
110 40 : if (_elem_type == "QUAD4")
111 10 : num_nodes = 8;
112 30 : else if (_elem_type == "QUAD8")
113 10 : num_nodes = 20;
114 20 : else if (_elem_type == "HEX8")
115 10 : num_nodes = 16;
116 10 : else if (_elem_type == "HEX20")
117 10 : num_nodes = 48;
118 :
119 40 : std::vector<Node *> nodes(num_nodes);
120 :
121 40 : const Real xmax = _xoffset + _xlength;
122 40 : const Real ymax = _yoffset + _ylength;
123 40 : const Real zmax = _zoffset + _zlength;
124 :
125 40 : if (_dim == 2)
126 : {
127 20 : const std::vector<Real> x_interior_fractions{1.0 / 6.0, 3.0 / 4.0, 2.0 / 3.0, 1.0 / 3.0};
128 20 : 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 20 : std::vector<Point> node_positions{{_xoffset, _yoffset, _zoffset},
132 20 : {xmax, _yoffset, _zoffset},
133 20 : {xmax, ymax, _zoffset},
134 20 : {_xoffset, ymax, _zoffset}};
135 : // Interior node positions
136 100 : for (unsigned i = 0; i < x_interior_fractions.size(); ++i)
137 160 : node_positions.push_back({_xoffset + x_interior_fractions[i] * _xlength,
138 80 : _yoffset + y_interior_fractions[i] * _ylength,
139 80 : _zoffset});
140 :
141 20 : if (_elem_type == "QUAD8")
142 : {
143 : // Exterior midside node positions
144 50 : for (unsigned i = 0; i < 4; ++i)
145 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4]));
146 :
147 : // Exterior to interior midside node positions
148 50 : for (unsigned i = 0; i < 4; ++i)
149 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
150 :
151 : // Interior midside node positions
152 50 : for (unsigned i = 4; i < 8; ++i)
153 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 4]));
154 : }
155 :
156 300 : for (unsigned node_id = 0; node_id < num_nodes; ++node_id)
157 280 : nodes[node_id] = mesh->add_point(node_positions[node_id], node_id);
158 :
159 20 : if (_elem_type == "QUAD4")
160 10 : makeQuad4Elems(*mesh, nodes);
161 : else
162 10 : makeQuad8Elems(*mesh, nodes);
163 :
164 20 : boundary_info.sideset_name(1) = "bottom";
165 20 : boundary_info.sideset_name(2) = "right";
166 20 : boundary_info.sideset_name(3) = "top";
167 20 : boundary_info.sideset_name(4) = "left";
168 :
169 20 : boundary_info.nodeset_name(100) = "bottom_left";
170 20 : boundary_info.nodeset_name(101) = "bottom_right";
171 20 : boundary_info.nodeset_name(102) = "top_right";
172 20 : boundary_info.nodeset_name(103) = "top_left";
173 20 : }
174 : else
175 : {
176 : const std::vector<Real> x_interior_fractions{
177 20 : 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 20 : 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 20 : 0.192, 0.288, 0.263, 0.230, 0.643, 0.683, 0.644, 0.702};
182 :
183 : // Exterior node positions
184 20 : std::vector<Point> node_positions{{_xoffset, _yoffset, _zoffset},
185 20 : {xmax, _yoffset, _zoffset},
186 20 : {xmax, ymax, _zoffset},
187 20 : {_xoffset, ymax, _zoffset},
188 20 : {_xoffset, _yoffset, zmax},
189 20 : {xmax, _yoffset, zmax},
190 : {xmax, ymax, zmax},
191 20 : {_xoffset, ymax, zmax}};
192 :
193 : // Interior node positions
194 180 : for (unsigned i = 0; i < x_interior_fractions.size(); ++i)
195 320 : node_positions.push_back({_xoffset + x_interior_fractions[i] * _xlength,
196 160 : _yoffset + y_interior_fractions[i] * _ylength,
197 160 : _zoffset + z_interior_fractions[i] * _zlength});
198 :
199 20 : if (_elem_type == "HEX20")
200 : {
201 : // Midside Nodes
202 :
203 : // Four on back face
204 50 : for (unsigned i = 0; i < 4; ++i)
205 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4]));
206 :
207 : // Four on front face
208 50 : for (unsigned i = 4; i < 8; ++i)
209 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 4]));
210 :
211 : // Four on remaining exterior edges
212 50 : for (unsigned i = 0; i < 4; ++i)
213 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
214 :
215 : // Four on interior hex back face
216 50 : for (unsigned i = 8; i < 12; ++i)
217 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 8]));
218 :
219 : // Four on interior hex front face
220 50 : for (unsigned i = 12; i < 16; ++i)
221 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[(i + 1) % 4 + 12]));
222 :
223 : // Four on remaining interior hex edges
224 50 : for (unsigned i = 8; i < 12; ++i)
225 40 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 4]));
226 :
227 : // Eight on remaining interior edges
228 90 : for (unsigned i = 0; i < 8; ++i)
229 80 : node_positions.push_back(0.5 * (node_positions[i] + node_positions[i + 8]));
230 : }
231 :
232 660 : for (unsigned node_id = 0; node_id < num_nodes; ++node_id)
233 640 : nodes[node_id] = mesh->add_point(node_positions[node_id], node_id);
234 :
235 20 : if (_elem_type == "HEX8")
236 10 : makeHex8Elems(*mesh, nodes);
237 : else
238 10 : makeHex20Elems(*mesh, nodes);
239 :
240 20 : boundary_info.sideset_name(1) = "front";
241 20 : boundary_info.sideset_name(2) = "bottom";
242 20 : boundary_info.sideset_name(3) = "left";
243 20 : boundary_info.sideset_name(4) = "right";
244 20 : boundary_info.sideset_name(5) = "top";
245 20 : boundary_info.sideset_name(6) = "back";
246 :
247 20 : boundary_info.nodeset_name(100) = "bottom_back_left";
248 20 : boundary_info.nodeset_name(101) = "bottom_back_right";
249 20 : boundary_info.nodeset_name(102) = "top_back_right";
250 20 : boundary_info.nodeset_name(103) = "top_back_left";
251 20 : boundary_info.nodeset_name(104) = "bottom_front_left";
252 20 : boundary_info.nodeset_name(105) = "bottom_front_right";
253 20 : boundary_info.nodeset_name(106) = "top_front_right";
254 20 : boundary_info.nodeset_name(107) = "top_front_left";
255 20 : }
256 :
257 40 : mesh->prepare_for_use();
258 :
259 80 : return dynamic_pointer_cast<MeshBase>(mesh);
260 40 : }
261 :
262 : void
263 10 : ExamplePatchMeshGenerator::makeQuad4Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
264 : {
265 10 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
266 : std::vector<std::vector<int>> element_connectivity{
267 60 : {0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}};
268 60 : for (unsigned i = 0; i < 5; ++i)
269 : {
270 50 : Elem * elem = mesh.add_elem(new Quad4);
271 250 : for (unsigned j = 0; j < 4; ++j)
272 200 : elem->set_node(j, nodes[element_connectivity[i][j]]);
273 :
274 50 : elem->subdomain_id() = i + 1;
275 :
276 50 : if (i < 4)
277 : {
278 40 : boundary_info.add_node(i, i + 100);
279 40 : boundary_info.add_side(elem, 0, i + 1);
280 : }
281 : }
282 30 : }
283 :
284 : void
285 10 : ExamplePatchMeshGenerator::makeQuad8Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
286 : {
287 10 : 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 60 : {4, 5, 6, 7, 16, 17, 18, 19}};
293 60 : for (unsigned i = 0; i < 5; ++i)
294 : {
295 50 : Elem * elem = mesh.add_elem(new Quad8);
296 450 : for (unsigned j = 0; j < 8; ++j)
297 400 : elem->set_node(j, nodes[element_connectivity[i][j]]);
298 :
299 50 : elem->subdomain_id() = i + 1;
300 :
301 50 : if (i < 4)
302 : {
303 40 : boundary_info.add_node(i, i + 100);
304 40 : boundary_info.add_side(elem, 0, i + 1);
305 : }
306 : }
307 30 : }
308 :
309 : void
310 10 : ExamplePatchMeshGenerator::makeHex8Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
311 : {
312 10 : 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 80 : {8, 9, 10, 11, 12, 13, 14, 15}};
320 :
321 80 : for (unsigned i = 0; i < 7; ++i)
322 : {
323 70 : Elem * elem = mesh.add_elem(new Hex8);
324 630 : for (unsigned j = 0; j < 8; ++j)
325 560 : elem->set_node(j, nodes[element_connectivity[i][j]]);
326 :
327 70 : elem->subdomain_id() = i + 1;
328 70 : if (i < 6)
329 60 : boundary_info.add_side(elem, 5, i + 1);
330 : }
331 90 : for (unsigned i = 0; i < 8; ++i)
332 80 : boundary_info.add_node(i, i + 100);
333 30 : }
334 :
335 : void
336 10 : ExamplePatchMeshGenerator::makeHex20Elems(MeshBase & mesh, const std::vector<Node *> & nodes)
337 : {
338 10 : 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 80 : {8, 9, 10, 11, 12, 13, 14, 15, 28, 29, 30, 31, 36, 37, 38, 39, 32, 33, 34, 35}};
347 :
348 80 : for (unsigned i = 0; i < 7; ++i)
349 : {
350 70 : Elem * elem = mesh.add_elem(new Hex20);
351 1470 : for (unsigned j = 0; j < 20; ++j)
352 1400 : elem->set_node(j, nodes[element_connectivity[i][j]]);
353 :
354 70 : elem->subdomain_id() = i + 1;
355 70 : if (i < 6)
356 60 : boundary_info.add_side(elem, 5, i + 1);
357 : }
358 90 : for (unsigned i = 0; i < 8; ++i)
359 80 : boundary_info.add_node(i, i + 100);
360 30 : }
|
