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 "AdvancedExtruderGenerator.h"
11 : #include "MooseUtils.h"
12 : #include "MooseMeshUtils.h"
13 : #include "MathUtils.h"
14 : #include "GeometryUtils.h"
15 :
16 : #include "libmesh/boundary_info.h"
17 : #include "libmesh/function_base.h"
18 : #include "libmesh/cell_prism6.h"
19 : #include "libmesh/cell_prism18.h"
20 : #include "libmesh/cell_prism21.h"
21 : #include "libmesh/cell_hex8.h"
22 : #include "libmesh/cell_hex20.h"
23 : #include "libmesh/cell_hex27.h"
24 : #include "libmesh/cell_c0polyhedron.h"
25 : #include "libmesh/edge_edge2.h"
26 : #include "libmesh/edge_edge3.h"
27 : #include "libmesh/edge_edge4.h"
28 : #include "libmesh/face_quad4.h"
29 : #include "libmesh/face_quad8.h"
30 : #include "libmesh/face_quad9.h"
31 : #include "libmesh/face_tri3.h"
32 : #include "libmesh/face_tri6.h"
33 : #include "libmesh/face_tri7.h"
34 : #include "libmesh/face_c0polygon.h"
35 : #include "libmesh/libmesh_logging.h"
36 : #include "libmesh/mesh_communication.h"
37 : #include "libmesh/mesh_modification.h"
38 : #include "libmesh/mesh_serializer.h"
39 : #include "libmesh/mesh_tools.h"
40 : #include "libmesh/parallel.h"
41 : #include "libmesh/remote_elem.h"
42 : #include "libmesh/string_to_enum.h"
43 : #include "libmesh/unstructured_mesh.h"
44 : #include "libmesh/point.h"
45 :
46 : #include <numeric>
47 : #include <cmath>
48 :
49 : registerMooseObject("MooseApp", AdvancedExtruderGenerator);
50 : registerMooseObjectRenamed("MooseApp",
51 : FancyExtruderGenerator,
52 : "02/18/2023 24:00",
53 : AdvancedExtruderGenerator);
54 :
55 : InputParameters
56 6919 : AdvancedExtruderGenerator::validParams()
57 : {
58 6919 : InputParameters params = MeshGenerator::validParams();
59 :
60 13838 : params.addClassDescription(
61 : "Extrudes a 1D mesh into 2D, or a 2D mesh into 3D, and supports a variable height for each "
62 : "elevation, variable number of layers within each elevation, variable growth factors of "
63 : "axial element sizes within each elevation and remap subdomain_ids, boundary_ids and element "
64 : "extra integers within each elevation as well as interface boundaries between neighboring "
65 : "elevation layers, as well as following a 1D curve and modifying the radial (normal to "
66 : "the extrusion axis) extent of the geometry.");
67 :
68 27676 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh to extrude");
69 :
70 : // User input of extrusion heights / axial discretization
71 27676 : params.addParam<std::vector<Real>>("heights", {}, "The height of each elevation");
72 41514 : params.addRangeCheckedParam<std::vector<Real>>(
73 : "biases", "biases>0.0", "The axial growth factor used for mesh biasing for each elevation.");
74 27676 : params.addParam<std::vector<unsigned int>>(
75 : "num_layers",
76 : {},
77 : "The number of layers for each elevation - must be num_elevations in length!");
78 :
79 : // Swaps on every height
80 27676 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
81 : "subdomain_swaps",
82 : {},
83 : "For each row, every two entries are interpreted as a pair of "
84 : "'from' and 'to' to remap the subdomains for that elevation");
85 27676 : params.addParam<std::vector<std::vector<boundary_id_type>>>(
86 : "boundary_swaps",
87 : {},
88 : "For each row, every two entries are interpreted as a pair of "
89 : "'from' and 'to' to remap the boundaries for that elevation");
90 27676 : params.addParam<std::vector<std::string>>(
91 : "elem_integer_names_to_swap",
92 : {},
93 : "Array of element extra integer names that need to be swapped during extrusion.");
94 27676 : params.addParam<std::vector<std::vector<std::vector<dof_id_type>>>>(
95 : "elem_integers_swaps",
96 : {},
97 : "For each row, every two entries are interpreted as a pair of 'from' and 'to' to remap the "
98 : "element extra integer for that elevation. If multiple element extra integers need to be "
99 : "swapped, the enties are stacked based on the order provided in "
100 : "'elem_integer_names_to_swap' to form the third dimension.");
101 :
102 : // Direction parameter
103 27676 : params.addParam<Point>("direction",
104 : "A vector that points in the direction to extrude (note, this will be "
105 : "normalized internally - so don't worry about it here)");
106 27676 : params.addParam<MeshGeneratorName>(
107 : "extrusion_curve",
108 : "Name of the mesh generator providing the line mesh curve to be extruded along. The "
109 : "extrusion path follows the node order in the line mesh");
110 27676 : params.addParam<RealVectorValue>(
111 : "start_extrusion_direction",
112 : "A vector that points in the starting direction for extruding along a curve. This vector "
113 : "should be the tangent vector at the FIRST node of the extrusion curve. Vector will be "
114 : "normalized in code, so don't worry about it here.");
115 27676 : params.addParam<RealVectorValue>(
116 : "end_extrusion_direction",
117 : "A vector that points in the ending direction for extruding along a curve. This vector "
118 : "should be the tangent vector at the LAST node of the extrusion curve. Vector will be "
119 : "normalized in code, so don't worry about it here.");
120 :
121 : // Boundaries and interfaces
122 27676 : params.addParam<BoundaryName>(
123 : "top_boundary",
124 : "The boundary name to set on the top boundary. If omitted an ID will be generated.");
125 27676 : params.addParam<BoundaryName>(
126 : "bottom_boundary",
127 : "The boundary name to set on the bottom boundary. If omitted an ID will be generated.");
128 27676 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
129 : "upward_boundary_source_blocks", "Block ids used to generate upward interface boundaries.");
130 27676 : params.addParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids",
131 : "Upward interface boundary ids.");
132 27676 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
133 : "downward_boundary_source_blocks",
134 : "Block ids used to generate downward interface boundaries.");
135 27676 : params.addParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids",
136 : "Downward interface boundary ids.");
137 :
138 : // Radial transformations
139 20757 : params.addParam<Real>("twist_pitch",
140 13838 : 0,
141 : "Pitch for helicoidal extrusion around an axis going through the origin "
142 : "following the direction vector");
143 20757 : params.addParam<Real>("end_radial_extent",
144 13838 : 0,
145 : "If modifying the radial extent of the extruded geometry, final radial "
146 : "extent to reach at the end of the extrusion process");
147 27676 : MooseEnum radial_growth_methods("linear cubic", "linear");
148 27676 : params.addParam<MooseEnum>("radial_growth_method",
149 : radial_growth_methods,
150 : "Functional form to change radius while extruding along curve.");
151 27676 : params.addParam<Real>("start_radial_growth_rate", 1., "Starting rate of radial expansion.");
152 27676 : params.addParam<Real>("end_radial_growth_rate", 1., "Ending rate of radial expansion.");
153 :
154 27676 : params.addParamNamesToGroup(
155 : "top_boundary bottom_boundary upward_boundary_source_blocks upward_boundary_ids "
156 : "downward_boundary_source_blocks downward_boundary_ids",
157 : "Boundary Assignment");
158 27676 : params.addParamNamesToGroup(
159 : "subdomain_swaps boundary_swaps elem_integer_names_to_swap elem_integers_swaps", "ID Swap");
160 27676 : params.addParamNamesToGroup("extrusion_curve start_extrusion_direction end_extrusion_direction",
161 : "Extrusion along curve");
162 20757 : params.addParamNamesToGroup("twist_pitch end_radial_extent radial_growth_method "
163 : "start_radial_growth_rate end_radial_growth_rate",
164 : "Radial transformation");
165 :
166 13838 : return params;
167 6919 : }
168 :
169 406 : AdvancedExtruderGenerator::AdvancedExtruderGenerator(const InputParameters & parameters)
170 : : MeshGenerator(parameters),
171 406 : _input(getMesh("input")),
172 812 : _heights(getParam<std::vector<Real>>("heights")),
173 1637 : _biases(isParamValid("biases") ? getParam<std::vector<Real>>("biases")
174 799 : : std::vector<Real>(_heights.size(), 1.0)),
175 812 : _num_layers(getParam<std::vector<unsigned int>>("num_layers")),
176 812 : _subdomain_swaps(getParam<std::vector<std::vector<subdomain_id_type>>>("subdomain_swaps")),
177 812 : _boundary_swaps(getParam<std::vector<std::vector<boundary_id_type>>>("boundary_swaps")),
178 812 : _elem_integer_names_to_swap(getParam<std::vector<std::string>>("elem_integer_names_to_swap")),
179 406 : _elem_integers_swaps(
180 812 : getParam<std::vector<std::vector<std::vector<dof_id_type>>>>("elem_integers_swaps")),
181 1494 : _direction(isParamValid("direction") ? getParam<Point>("direction") : Point(0, 0, 0)),
182 812 : _extrusion_curve(getMesh("extrusion_curve", true)),
183 1218 : _start_extrusion_direction(isParamValid("start_extrusion_direction")
184 1286 : ? getParam<RealVectorValue>("start_extrusion_direction").unit()
185 744 : : Point(0, 0, 0)),
186 1218 : _end_extrusion_direction(isParamValid("end_extrusion_direction")
187 1286 : ? getParam<RealVectorValue>("end_extrusion_direction").unit()
188 744 : : Point(0, 0, 0)),
189 812 : _extrude_along_curve(isParamValid("extrusion_curve")),
190 812 : _has_top_boundary(isParamValid("top_boundary")),
191 1318 : _top_boundary(isParamValid("top_boundary") ? getParam<BoundaryName>("top_boundary") : "0"),
192 812 : _has_bottom_boundary(isParamValid("bottom_boundary")),
193 1318 : _bottom_boundary(isParamValid("bottom_boundary") ? getParam<BoundaryName>("bottom_boundary")
194 : : "0"),
195 1218 : _upward_boundary_source_blocks(
196 1218 : isParamValid("upward_boundary_source_blocks")
197 406 : ? getParam<std::vector<std::vector<subdomain_id_type>>>("upward_boundary_source_blocks")
198 393 : : std::vector<std::vector<subdomain_id_type>>(_heights.size(),
199 799 : std::vector<subdomain_id_type>())),
200 1218 : _upward_boundary_ids(
201 1218 : isParamValid("upward_boundary_ids")
202 406 : ? getParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids")
203 393 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
204 799 : std::vector<boundary_id_type>())),
205 2043 : _downward_boundary_source_blocks(isParamValid("downward_boundary_source_blocks")
206 406 : ? getParam<std::vector<std::vector<subdomain_id_type>>>(
207 : "downward_boundary_source_blocks")
208 : : std::vector<std::vector<subdomain_id_type>>(
209 799 : _heights.size(), std::vector<subdomain_id_type>())),
210 1218 : _downward_boundary_ids(
211 1218 : isParamValid("downward_boundary_ids")
212 406 : ? getParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids")
213 393 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
214 799 : std::vector<boundary_id_type>())),
215 812 : _twist_pitch(getParam<Real>("twist_pitch")),
216 812 : _end_radial_extent(getParam<Real>("end_radial_extent")),
217 812 : _radial_expansion_method(getParam<MooseEnum>("radial_growth_method")),
218 812 : _start_radial_growth_rate(getParam<Real>("start_radial_growth_rate")),
219 2030 : _end_radial_growth_rate(getParam<Real>("end_radial_growth_rate"))
220 : {
221 406 : if (_extrude_along_curve)
222 : {
223 204 : if (isParamSetByUser("heights"))
224 6 : paramError("heights", "heights cannot be set if extruding along curve!");
225 195 : if (isParamValid("biases"))
226 6 : paramError("biases", "biases cannot be set if extruding along curve!");
227 186 : if (isParamSetByUser("num_layers"))
228 6 : paramError("num_layers", "num_layers cannot be set if extruding along curve!");
229 177 : if (isParamValid("direction"))
230 6 : paramError("direction", "direction cannot be set if extruding along curve!");
231 : }
232 : else
233 : {
234 338 : if (!_direction.norm())
235 0 : paramError("direction", "Must have some length!");
236 338 : _direction /= _direction.norm();
237 : }
238 :
239 : unsigned int num_elevations;
240 394 : if (_extrude_along_curve)
241 56 : num_elevations = 1;
242 : else
243 : {
244 338 : num_elevations = _heights.size();
245 338 : if (_num_layers.size() != num_elevations)
246 0 : paramError(
247 0 : "heights", "The length of 'heights' and 'num_layers' must be the same in ", name());
248 : }
249 :
250 394 : if (_subdomain_swaps.size() && (_subdomain_swaps.size() != num_elevations))
251 : {
252 0 : if (_extrude_along_curve)
253 0 : paramError("subdomain_swaps",
254 0 : "If specified, 'subdomain_swaps' (" + std::to_string(_subdomain_swaps.size()) +
255 : ") must be equal to 1 when extruding along a curve.");
256 : else
257 : {
258 0 : paramError("subdomain_swaps",
259 0 : "If specified, 'subdomain_swaps' (" + std::to_string(_subdomain_swaps.size()) +
260 0 : ") must be the same length as 'heights' (" + std::to_string(num_elevations) +
261 : ") in ",
262 0 : name());
263 : }
264 : }
265 :
266 : try
267 : {
268 394 : MooseMeshUtils::idSwapParametersProcessor(
269 394 : name(), "subdomain_swaps", _subdomain_swaps, _subdomain_swap_pairs);
270 : }
271 0 : catch (const MooseException & e)
272 : {
273 0 : paramError("subdomain_swaps", e.what());
274 0 : }
275 :
276 394 : if (_boundary_swaps.size() && (_boundary_swaps.size() != num_elevations))
277 : {
278 0 : if (_extrude_along_curve)
279 : {
280 0 : paramError("boundary_swaps",
281 0 : "If specified, 'boundary_swaps' (" + std::to_string(_boundary_swaps.size()) +
282 0 : ") must be the same length as 'heights' (" + std::to_string(num_elevations) +
283 : ") in ",
284 0 : name());
285 : }
286 : else
287 : {
288 0 : paramError("boundary_swaps",
289 0 : "If specified, 'boundary_swaps' (" + std::to_string(_boundary_swaps.size()) +
290 : ") must be equal to 1 when extruding along a curve.");
291 : }
292 : }
293 :
294 : try
295 : {
296 394 : MooseMeshUtils::idSwapParametersProcessor(
297 394 : name(), "boundary_swaps", _boundary_swaps, _boundary_swap_pairs);
298 : }
299 0 : catch (const MooseException & e)
300 : {
301 0 : paramError("boundary_swaps", e.what());
302 0 : }
303 :
304 406 : if (_elem_integers_swaps.size() &&
305 12 : _elem_integers_swaps.size() != _elem_integer_names_to_swap.size())
306 0 : paramError("elem_integers_swaps",
307 : "If specified, 'elem_integers_swaps' must have the same length as the length of "
308 : "'elem_integer_names_to_swap'.");
309 :
310 418 : for (const auto & unit_elem_integers_swaps : _elem_integers_swaps)
311 24 : if (unit_elem_integers_swaps.size() != num_elevations)
312 0 : paramError("elem_integers_swaps",
313 : "If specified, each element of 'elem_integers_swaps' must have the same length as "
314 : "the length of 'heights'.");
315 :
316 : try
317 : {
318 394 : MooseMeshUtils::extraElemIntegerSwapParametersProcessor(name(),
319 : num_elevations,
320 394 : _elem_integer_names_to_swap.size(),
321 : _elem_integers_swaps,
322 394 : _elem_integers_swap_pairs);
323 : }
324 0 : catch (const MooseException & e)
325 : {
326 0 : paramError("elem_integers_swaps", e.what());
327 0 : }
328 :
329 394 : if (!_extrude_along_curve)
330 : {
331 338 : bool has_negative_entry = false;
332 338 : bool has_positive_entry = false;
333 960 : for (const auto & h : _heights)
334 : {
335 622 : if (h > 0.0)
336 619 : has_positive_entry = true;
337 : else
338 3 : has_negative_entry = true;
339 : }
340 :
341 338 : if (has_negative_entry && has_positive_entry)
342 6 : paramError("heights", "Cannot have both positive and negative heights!");
343 335 : if (_biases.size() != _heights.size())
344 0 : paramError("biases", "Size of this parameter, if provided, must be the same as heights.");
345 : }
346 :
347 1061 : if ((_upward_boundary_source_blocks.size() || _upward_boundary_ids.size()) &&
348 335 : (_upward_boundary_source_blocks.size() != _upward_boundary_ids.size() ||
349 335 : _upward_boundary_ids.size() != num_elevations))
350 : {
351 0 : if (_extrude_along_curve)
352 0 : paramError(
353 : "upward_boundary_ids",
354 0 : "This parameter must have the same length (" +
355 0 : std::to_string(_upward_boundary_ids.size()) + ") as upward_boundary_source_blocks (" +
356 0 : std::to_string(_upward_boundary_source_blocks.size()) + ") and num_elevations (" +
357 0 : std::to_string(num_elevations) +
358 : "). Note that the number of heights is set to 1 when extruding along a curve.");
359 : else
360 : {
361 0 : paramError("upward_boundary_ids",
362 0 : "This parameter must have the same length (" +
363 0 : std::to_string(_upward_boundary_ids.size()) +
364 0 : ") as upward_boundary_source_blocks (" +
365 0 : std::to_string(_upward_boundary_source_blocks.size()) + ") and heights (" +
366 0 : std::to_string(num_elevations) + ")");
367 : }
368 : }
369 :
370 1004 : for (unsigned int i = 0; i < _upward_boundary_source_blocks.size(); i++)
371 613 : if (_upward_boundary_source_blocks[i].size() != _upward_boundary_ids[i].size())
372 0 : paramError("upward_boundary_ids",
373 : "Every element of this parameter must have the same length as the corresponding "
374 : "element of upward_boundary_source_blocks.");
375 :
376 1061 : if ((_downward_boundary_source_blocks.size() || _downward_boundary_ids.size()) &&
377 335 : (_downward_boundary_source_blocks.size() != _downward_boundary_ids.size() ||
378 335 : _downward_boundary_ids.size() != num_elevations))
379 : {
380 0 : if (_extrude_along_curve)
381 0 : paramError(
382 : "downward_boundary_ids",
383 0 : "This parameter must have the same length (" +
384 0 : std::to_string(_downward_boundary_ids.size()) +
385 0 : ") as downward_boundary_source_blocks (" +
386 0 : std::to_string(_downward_boundary_source_blocks.size()) + ") and (" +
387 0 : std::to_string(num_elevations) +
388 : "). Note that the number of heights is set to 1 when extruding along a curve.");
389 : else
390 : {
391 0 : paramError("downward_boundary_ids",
392 0 : "This parameter must have the same length (" +
393 0 : std::to_string(_downward_boundary_ids.size()) +
394 0 : ") as downward_boundary_source_blocks (" +
395 0 : std::to_string(_downward_boundary_source_blocks.size()) + ") and heights (" +
396 0 : std::to_string(num_elevations) + ")");
397 : }
398 : }
399 1004 : for (unsigned int i = 0; i < _downward_boundary_source_blocks.size(); i++)
400 613 : if (_downward_boundary_source_blocks[i].size() != _downward_boundary_ids[i].size())
401 0 : paramError("downward_boundary_ids",
402 : "Every element of this parameter must have the same length as the corresponding "
403 : "element of downward_boundary_source_blocks.");
404 391 : }
405 :
406 : std::unique_ptr<MeshBase>
407 383 : AdvancedExtruderGenerator::generate()
408 : {
409 : // Note: bulk of this code originally from libmesh mesh_modification.C
410 : // Original copyright: Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
411 : // Original license is LGPL so it can be used here.
412 :
413 383 : auto mesh = buildMeshBaseObject(_input->mesh_dimension() + 1);
414 383 : mesh->set_mesh_dimension(_input->mesh_dimension() + 1);
415 :
416 : // Check if the element integer names are existent in the input mesh.
417 405 : for (const auto i : make_range(_elem_integer_names_to_swap.size()))
418 22 : if (_input->has_elem_integer(_elem_integer_names_to_swap[i]))
419 44 : _elem_integer_indices_to_swap.push_back(
420 22 : _input->get_elem_integer_index(_elem_integer_names_to_swap[i]));
421 : else
422 0 : paramError("elem_integer_names_to_swap",
423 : "Element ",
424 : i + 1,
425 : " of 'elem_integer_names_to_swap' in is not a valid extra element integer of the "
426 : "input mesh.");
427 :
428 : // prepare for transferring extra element integers from original mesh to the extruded mesh.
429 383 : const unsigned int num_extra_elem_integers = _input->n_elem_integers();
430 383 : std::vector<std::string> id_names;
431 :
432 405 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
433 : {
434 22 : id_names.push_back(_input->get_elem_integer_name(i));
435 22 : if (!mesh->has_elem_integer(id_names[i]))
436 22 : mesh->add_elem_integer(id_names[i]);
437 : }
438 :
439 : // retrieve subdomain/sideset/nodeset name maps
440 383 : if (!_input->preparation().has_cached_elem_data)
441 281 : _input->cache_elem_data();
442 383 : const auto & input_subdomain_map = _input->get_subdomain_name_map();
443 383 : const auto & input_sideset_map = _input->get_boundary_info().get_sideset_name_map();
444 383 : const auto & input_nodeset_map = _input->get_boundary_info().get_nodeset_name_map();
445 :
446 : // Check that the swaps source blocks are present in the mesh
447 668 : for (const auto & swap : _subdomain_swaps)
448 1242 : for (const auto i : index_range(swap))
449 957 : if (i % 2 == 0 && !MooseMeshUtils::hasSubdomainID(*_input, swap[i]))
450 6 : paramError("subdomain_swaps", "The block '", swap[i], "' was not found within the mesh");
451 :
452 : // Check that the swaps source boundaries are present in the mesh
453 398 : for (const auto & swap : _boundary_swaps)
454 117 : for (const auto i : index_range(swap))
455 99 : if (i % 2 == 0 && !MooseMeshUtils::hasBoundaryID(*_input, swap[i]))
456 6 : paramError("boundary_swaps", "The boundary '", swap[i], "' was not found within the mesh");
457 :
458 : // Check that the source blocks for layer top/bottom boundaries exist in the mesh
459 941 : for (const auto & layer_vec : _upward_boundary_source_blocks)
460 627 : for (const auto bid : layer_vec)
461 63 : if (!MooseMeshUtils::hasSubdomainID(*_input, bid))
462 6 : paramError(
463 : "upward_boundary_source_blocks", "The block '", bid, "' was not found within the mesh");
464 929 : for (const auto & layer_vec : _downward_boundary_source_blocks)
465 618 : for (const auto bid : layer_vec)
466 63 : if (!MooseMeshUtils::hasSubdomainID(*_input, bid))
467 6 : paramError("downward_boundary_source_blocks",
468 : "The block '",
469 : bid,
470 : "' was not found within the mesh");
471 :
472 : // Move the meshes as requested by the mesh generator system
473 371 : std::unique_ptr<MeshBase> input = std::move(_input);
474 371 : std::unique_ptr<MeshBase> extrusion_curve;
475 371 : if (_extrude_along_curve)
476 56 : extrusion_curve = std::move(_extrusion_curve);
477 :
478 : // We must serialize the curve to know how to extrude across all ranks
479 371 : std::unique_ptr<libMesh::MeshSerializer> serializer;
480 371 : if (_extrude_along_curve)
481 56 : serializer = std::make_unique<libMesh::MeshSerializer>(*extrusion_curve);
482 :
483 : // If we're using a distributed mesh... then make sure we don't have any remote elements hanging
484 : // around
485 371 : if (!input->is_serial())
486 : {
487 27 : input->delete_remote_elements();
488 : // This should be a no-op, and yet, it's needed for two THM tests
489 27 : mesh->delete_remote_elements();
490 : }
491 :
492 371 : if (input->n_nodes() != input->max_node_id())
493 0 : input->renumber_nodes_and_elements();
494 :
495 371 : if (input->n_nodes() != input->max_node_id())
496 0 : mooseError(
497 : "You must allow renumbering, because the extruded mesh should be contiguously numbered. "
498 : "Alternatively, you can use a separate mesh generator (MeshRepairGenerator with the "
499 : "renumber_contiguously parameter for example) to renumber the nodes contiguously.");
500 :
501 : unsigned int total_num_layers;
502 : unsigned int total_num_elevations;
503 371 : if (!_extrude_along_curve)
504 : {
505 315 : total_num_layers = std::accumulate(_num_layers.begin(), _num_layers.end(), 0);
506 315 : total_num_elevations = _heights.size();
507 : }
508 : else
509 : {
510 56 : total_num_layers = extrusion_curve->n_elem();
511 56 : total_num_elevations = 1;
512 : }
513 :
514 371 : dof_id_type orig_elem = input->n_elem();
515 371 : dof_id_type orig_nodes = input->n_nodes();
516 :
517 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
518 371 : unique_id_type orig_unique_ids = input->parallel_max_unique_id();
519 371 : bool has_poly_midnodes = false;
520 : #endif
521 :
522 371 : bool has_polygons = false;
523 371 : unsigned int order = 1;
524 :
525 371 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
526 371 : const BoundaryInfo & input_boundary_info = input->get_boundary_info();
527 :
528 : // Determine boundary IDs for the new user provided boundary names
529 371 : std::vector<BoundaryName> new_boundary_names;
530 371 : if (_has_bottom_boundary)
531 232 : new_boundary_names.push_back(_bottom_boundary);
532 371 : if (_has_top_boundary)
533 232 : new_boundary_names.push_back(_top_boundary);
534 : std::vector<boundary_id_type> new_boundary_ids =
535 371 : MooseMeshUtils::getBoundaryIDs(*input, new_boundary_names, true);
536 : const auto user_bottom_boundary_id =
537 371 : _has_bottom_boundary ? new_boundary_ids.front() : libMesh::BoundaryInfo::invalid_id;
538 : const auto user_top_boundary_id =
539 371 : _has_top_boundary ? new_boundary_ids.back() : libMesh::BoundaryInfo::invalid_id;
540 :
541 : // We know a priori how many elements we'll need
542 371 : mesh->reserve_elem(total_num_layers * orig_elem);
543 :
544 : // Look for higher order elements which introduce an extra layer
545 742 : std::set<ElemType> higher_orders = {EDGE3, EDGE4, TRI6, TRI7, QUAD8, QUAD9};
546 371 : bool extruding_quad_eights = false;
547 371 : std::vector<ElemType> types;
548 371 : MeshTools::elem_types(*input, types);
549 767 : for (const auto elem_type : types)
550 : {
551 396 : if (higher_orders.count(elem_type))
552 8 : order = 2;
553 396 : if (elem_type == QUAD8)
554 8 : extruding_quad_eights = true;
555 : }
556 371 : mesh->comm().max(order);
557 371 : mesh->comm().max(extruding_quad_eights);
558 :
559 : // Reserve for the max number possibly needed
560 371 : mesh->reserve_nodes((order * total_num_layers + 1) * orig_nodes);
561 :
562 : // Container to catch the boundary IDs handed back by the BoundaryInfo object
563 371 : std::vector<boundary_id_type> ids_to_copy;
564 :
565 : // We need to compute the distance to the axis
566 371 : Real start_radial_extent = 0;
567 371 : Point reference_point;
568 371 : if (_extrude_along_curve)
569 56 : reference_point = *(extrusion_curve->node_ptr(0));
570 315 : else if (!MooseUtils::absoluteFuzzyEqual(_twist_pitch, 0.))
571 : // Backwards compatibility: we were twisting along an axis going through the origin
572 10 : reference_point = Point(0, 0, 0);
573 : else
574 305 : reference_point = MooseMeshUtils::meshCentroidCalculator(*input);
575 :
576 371 : if (_end_radial_extent)
577 : {
578 : // the center is at the curve if we are extruding along a curve, and the centroid of the mesh
579 : // otherwise
580 40 : RealVectorValue reference_direction;
581 40 : if (_extrude_along_curve)
582 30 : reference_direction =
583 30 : _start_extrusion_direction.norm_sq() > 0
584 30 : ? _start_extrusion_direction
585 : : RealVectorValue(
586 30 : (*(extrusion_curve->node_ptr(1)) - *(extrusion_curve->node_ptr(0))).unit());
587 : else
588 10 : reference_direction = _direction;
589 :
590 : // now we measure the initial radial extent
591 : start_radial_extent =
592 40 : MooseMeshUtils::computeMaxDistanceToAxis(*input, reference_point, reference_direction);
593 : }
594 :
595 : // Compute the total extrusion distance along the axis
596 : Real total_extrusion_distance_at_axis;
597 371 : if (_extrude_along_curve)
598 : {
599 56 : if (!extrusion_curve->is_prepared())
600 56 : extrusion_curve->prepare_for_use();
601 56 : total_extrusion_distance_at_axis = MeshTools::volume(*extrusion_curve);
602 : }
603 : else
604 315 : total_extrusion_distance_at_axis = std::accumulate(_heights.begin(), _heights.end(), 0);
605 :
606 : // Create translated layers of nodes in the direction of extrusion
607 21174 : for (const auto & node : input->node_ptr_range())
608 : {
609 20803 : unsigned int current_node_layer = 0;
610 20803 : Point orig_node_to_previous;
611 20803 : Point orig_node_to_current;
612 20803 : Real sum_step_sizes = 0.;
613 20803 : Real sum_step_sizes_at_axis = 0.;
614 : // Initial distance from the node to the extrusion axis
615 20803 : Real start_node_radius = (*node - reference_point).norm();
616 :
617 : // e is the elevation layer ordering
618 53474 : for (const auto e : make_range(total_num_elevations))
619 : {
620 32671 : unsigned int num_layers = libMesh::invalid_uint;
621 32671 : Real height = std::numeric_limits<Real>::max(), bias = std::numeric_limits<Real>::max();
622 32671 : if (_extrude_along_curve)
623 : {
624 440 : num_layers = extrusion_curve->n_elem();
625 : }
626 : else
627 : {
628 32231 : num_layers = _num_layers[e];
629 32231 : height = _heights[e];
630 32231 : bias = _biases[e];
631 : }
632 :
633 : // In first layer we also add the base nodes, hence the "plus one"
634 32671 : unsigned int num_heights_at_elevation = order * num_layers + (e == 0 ? 1 : 0);
635 : // Keep track of the plane normal to the extrusion
636 32671 : RealVectorValue prev_intersecting_plane_normal_vec = _start_extrusion_direction;
637 :
638 : // k is the element layer ordering within each elevation layer
639 126305 : for (const auto k : make_range(num_heights_at_elevation))
640 : {
641 : // Compute 'orig_node_to_current', the update vector
642 : // For the first layer we don't need to move
643 93634 : if (e == 0 && k == 0)
644 20803 : orig_node_to_current.zero();
645 : else
646 : {
647 : // Shift the previous position by a certain fraction of 'height' along the extrusion
648 : // direction to get the new position.
649 :
650 : // Compute orig_node_to_current before any transformation
651 72831 : Real step_size = 0;
652 72831 : Real step_size_at_axis = 0;
653 72831 : if (_extrude_along_curve)
654 : {
655 : // current point in extrusion curve
656 6288 : const Node * P_current = extrusion_curve->node_ptr(k);
657 : // previous point in extrusion curve
658 6288 : const Node * P_prev = extrusion_curve->node_ptr(k - 1);
659 :
660 : // Quantities for the previous position of the extruded node
661 6288 : const auto old_node = orig_node_to_previous + *node;
662 6288 : RealVectorValue b_vec = old_node - *P_prev;
663 6288 : Real node_distance_to_curve = b_vec.norm();
664 :
665 : // Node does not lie exactly on the extrusion curve we are following
666 6288 : if (node_distance_to_curve > libMesh::TOLERANCE)
667 : {
668 : // normal vector to the plane to extrude the point into
669 5408 : RealVectorValue intersecting_plane_normal_vec;
670 :
671 : // Select the extrusion direction based on the curve or the user parameters
672 5408 : if (k == 1)
673 : {
674 384 : const auto P_next = extrusion_curve->node_ptr(k + 1);
675 384 : intersecting_plane_normal_vec =
676 768 : 0.5 * (_start_extrusion_direction + (*P_next - *P_current).unit());
677 : }
678 5024 : else if (k < order * num_layers - 1)
679 : {
680 4256 : const auto P_next = extrusion_curve->node_ptr(k + 1);
681 : // this approximates the derivative at the spline point
682 4256 : intersecting_plane_normal_vec = *P_next - *P_prev;
683 : }
684 : else
685 : {
686 768 : intersecting_plane_normal_vec = (_end_extrusion_direction.norm_sq() > 0)
687 768 : ? _end_extrusion_direction
688 768 : : RealVectorValue(*P_current - *P_prev);
689 : }
690 5408 : intersecting_plane_normal_vec /= intersecting_plane_normal_vec.norm();
691 :
692 5408 : Point new_node_point;
693 : // If the extrusion direction and the previous plane normal are aligned,
694 : // we can't define a rotation axis. We simply translate the points
695 10816 : if (MooseUtils::absoluteFuzzyEqual(
696 0 : prev_intersecting_plane_normal_vec.cross(intersecting_plane_normal_vec)
697 5408 : .norm_sq(),
698 10816 : 0))
699 384 : new_node_point = old_node + *P_current - *P_prev;
700 : // We use a rotation from the previous extrusion plane normal to the current one
701 : // We have tried in the past:
702 : // - assuming (P_prev, P_current, old_node, current_node) are coplanar
703 : // - assuming (P_prev, P_current) and (old_node, current_node) are parallel
704 : // Both result in a slight deformation of the shape during extrusion.
705 : else
706 : {
707 5024 : auto axis = prev_intersecting_plane_normal_vec.cross(intersecting_plane_normal_vec);
708 5024 : const auto sin_th = axis.norm();
709 5024 : axis /= sin_th;
710 5024 : Real cos_th = prev_intersecting_plane_normal_vec * intersecting_plane_normal_vec;
711 : // Rodrigues formula for rotation
712 5024 : const auto new_v = cos_th * b_vec + axis.cross(b_vec) * sin_th +
713 10048 : axis * (axis * b_vec) * (1. - cos_th);
714 :
715 : mooseAssert(MooseUtils::absoluteFuzzyEqual(new_v.norm(), b_vec.norm()),
716 : "Radial extent be conserved");
717 5024 : new_node_point = *P_current + new_v;
718 : }
719 :
720 5408 : orig_node_to_current = new_node_point - *node;
721 :
722 5408 : step_size = (orig_node_to_current - orig_node_to_previous).norm();
723 5408 : step_size_at_axis = (*P_current - *P_prev).norm();
724 5408 : prev_intersecting_plane_normal_vec = intersecting_plane_normal_vec;
725 : }
726 : // Point is on the axis of the line
727 : else
728 : {
729 880 : _direction = *P_current - *P_prev;
730 880 : _direction /= _direction.norm();
731 : mooseAssert(std::abs(_direction.norm() - 1.0) < libMesh::TOLERANCE,
732 : "Norm of direction vector is not 1!");
733 :
734 : // Calculate step size.
735 : // Note: orig_node_to_previous + *node is the vector description of the
736 : // previously-created node
737 880 : step_size = ((*P_current - (orig_node_to_previous + *node)) * _direction);
738 880 : step_size_at_axis = step_size;
739 880 : orig_node_to_current = orig_node_to_previous + _direction * step_size;
740 : }
741 : }
742 : // Extruding in a fixed direction (not along a curve)
743 : else
744 : {
745 : // Divide by the order to avoid applying the bias to the nodes within a higher order
746 : // element
747 66543 : auto layer_index = (k - (e == 0 ? 1 : 0)) / order + 1;
748 133086 : step_size = MooseUtils::absoluteFuzzyEqual(bias, 1.0)
749 66543 : ? height / (Real)num_layers / (Real)order
750 2310 : : height * std::pow(bias, (Real)(layer_index - 1)) * (1.0 - bias) /
751 2310 : (1.0 - std::pow(bias, (Real)(num_layers))) / (Real)order;
752 66543 : step_size_at_axis = step_size;
753 66543 : orig_node_to_current =
754 66543 : orig_node_to_previous +
755 133086 : _direction * step_size; // update distance from starting node to new node
756 : }
757 :
758 : // Keep track of the cumulative step size as an extrusion axis coordinate
759 72831 : sum_step_sizes += step_size;
760 72831 : sum_step_sizes_at_axis += step_size_at_axis;
761 :
762 : // Handle radial expansion. No need to perform expansion at the centerline
763 72831 : if (_end_radial_extent && start_node_radius > libMesh::TOLERANCE)
764 : {
765 : // How far along we are in the extrusion, measured at the axis (=curve when extruding
766 : // along a curve)
767 6312 : Real tm1 =
768 6312 : (sum_step_sizes_at_axis - step_size_at_axis) / total_extrusion_distance_at_axis;
769 6312 : Real t = sum_step_sizes_at_axis / total_extrusion_distance_at_axis;
770 :
771 : // Direction of radial expansion.
772 6312 : RealVectorValue node_to_extrusion_axis;
773 6312 : if (_extrude_along_curve)
774 3600 : node_to_extrusion_axis =
775 7200 : *node + orig_node_to_current - *(extrusion_curve->node_ptr(k));
776 : // The radial expansion has to be performed in the rotated frame
777 2712 : else if (!MooseUtils::absoluteFuzzyEqual(_twist_pitch, 0.))
778 0 : node_to_extrusion_axis =
779 0 : *node + orig_node_to_current - (reference_point + _direction * sum_step_sizes);
780 : // when extruding in a straight line with no twisting, we can use the original radial
781 : // direction
782 : else
783 2712 : node_to_extrusion_axis = *node - reference_point;
784 :
785 : // Fraction of the total starting radius the node lives in
786 6312 : const auto radius_scaling = start_node_radius / start_radial_extent;
787 : // Calculate weighting for expansion
788 6312 : const auto radial_ratio_m1 = AdvancedExtruderGenerator::radialExpansionRatio(tm1);
789 6312 : const auto radial_ratio = AdvancedExtruderGenerator::radialExpansionRatio(t);
790 :
791 : // Compute change in step
792 6312 : orig_node_to_current += (start_radial_extent * (radial_ratio_m1 - radial_ratio) +
793 0 : (radial_ratio - radial_ratio_m1) * _end_radial_extent) *
794 6312 : radius_scaling * node_to_extrusion_axis.unit();
795 : }
796 :
797 : // Handle helicoidal extrusion
798 72831 : if (!MooseUtils::absoluteFuzzyEqual(_twist_pitch, 0.))
799 : {
800 : // twist 1 should be 'normal' to the extruded shape
801 2772 : RealVectorValue twist1 = _direction.cross(*node);
802 : // This happens for any node on the helicoidal extrusion axis
803 2772 : if (!MooseUtils::absoluteFuzzyEqual(twist1.norm(), .0))
804 2718 : twist1 /= twist1.norm();
805 2772 : const RealVectorValue twist2 = twist1.cross(_direction);
806 :
807 : auto twist =
808 0 : (cos(2. * libMesh::pi * current_node_layer * step_size / _twist_pitch) -
809 0 : cos(2. * libMesh::pi * (current_node_layer - 1) * step_size / _twist_pitch)) *
810 2772 : twist2 +
811 0 : (sin(2. * libMesh::pi * current_node_layer * step_size / _twist_pitch) -
812 0 : sin(2. * libMesh::pi * (current_node_layer - 1) * step_size / _twist_pitch)) *
813 5544 : twist1;
814 :
815 : // If we normalize twist, we must multiply by 2 * sin(libMesh::pi * step_size /
816 : // _twist_pitch) if (!MooseUtils::absoluteFuzzyEqual(twist.norm(), .0))
817 : // twist /= twist.norm();
818 :
819 : // Get a point on the extrusion axis (around which we currently twist the geometry)
820 : // at the local elevation height to be able to compute the distance to the axis
821 2772 : Point extrusion_axis_at_elevation;
822 2772 : Point prev_extrusion_axis_at_elevation;
823 2772 : if (_extrude_along_curve)
824 : {
825 0 : extrusion_axis_at_elevation = *extrusion_curve->node_ptr(k);
826 0 : prev_extrusion_axis_at_elevation = *extrusion_curve->node_ptr(k - 1);
827 : }
828 : else
829 : {
830 : // We can't use old and current step updates because they have been "twisted"
831 2772 : extrusion_axis_at_elevation = reference_point + _direction * sum_step_sizes;
832 2772 : prev_extrusion_axis_at_elevation =
833 5544 : reference_point + _direction * (sum_step_sizes - step_size);
834 : }
835 :
836 : // Scale with how far the node is from the extrusion axis, before twisting
837 2772 : twist *= (*node + orig_node_to_current - extrusion_axis_at_elevation).norm();
838 :
839 : // No need to twist or expand the point on the axis
840 2772 : if (!MooseUtils::absoluteFuzzyEqual(twist1.norm(), .0))
841 2718 : orig_node_to_current += twist;
842 : }
843 : }
844 :
845 : // Add the new node to the mesh
846 187268 : Node * new_node = mesh->add_point(*node + orig_node_to_current,
847 93634 : node->id() + (current_node_layer * orig_nodes),
848 93634 : node->processor_id());
849 :
850 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
851 : // Let's give the base of the extruded mesh the same
852 : // unique_ids as the source mesh, in case anyone finds that
853 : // a useful map to preserve.
854 : const unique_id_type uid = (current_node_layer == 0)
855 93634 : ? node->unique_id()
856 72831 : : orig_unique_ids +
857 72831 : (current_node_layer - 1) * (orig_elem + orig_nodes) +
858 72831 : node->id();
859 93634 : new_node->set_unique_id(uid);
860 : #endif
861 :
862 : // Add the new node to the extruded boundaries
863 93634 : input_boundary_info.boundary_ids(node, ids_to_copy);
864 93634 : if (_boundary_swap_pairs.empty())
865 93634 : boundary_info.add_node(new_node, ids_to_copy);
866 : else
867 0 : for (const auto & id_to_copy : ids_to_copy)
868 : {
869 0 : boundary_info.add_node(new_node,
870 0 : _boundary_swap_pairs[e].count(id_to_copy)
871 0 : ? _boundary_swap_pairs[e][id_to_copy]
872 0 : : id_to_copy);
873 : }
874 :
875 93634 : orig_node_to_previous = orig_node_to_current;
876 93634 : current_node_layer++;
877 : }
878 : }
879 371 : }
880 :
881 371 : const auto & side_ids = input_boundary_info.get_side_boundary_ids();
882 :
883 : boundary_id_type next_side_id =
884 371 : side_ids.empty() ? 0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
885 :
886 : // side_ids may not include ids from remote elements, in which case
887 : // some processors may have underestimated the next_side_id; let's
888 : // fix that.
889 371 : input->comm().max(next_side_id);
890 :
891 : // Map to keep track of polygon sides of polygonal prisms
892 371 : std::map<std::array<unsigned int, 3>, std::shared_ptr<libMesh::Polygon>> poly_extruded_sides;
893 :
894 : // Build the extruded elements
895 17651 : for (const auto & elem : input->element_ptr_range())
896 : {
897 17280 : const ElemType etype = elem->type();
898 :
899 : // build_extrusion currently only works on coarse meshes
900 : libmesh_assert(!elem->parent());
901 :
902 17280 : unsigned int current_layer = 0;
903 :
904 43908 : for (unsigned int e = 0; e != total_num_elevations; e++)
905 : {
906 26628 : auto num_layers = !_extrude_along_curve ? _num_layers[e] : extrusion_curve->n_nodes() - 1;
907 :
908 86090 : for (unsigned int k = 0; k != num_layers; ++k)
909 : {
910 59462 : std::unique_ptr<Elem> new_elem;
911 59462 : bool is_flipped(false);
912 59462 : switch (etype)
913 : {
914 4906 : case EDGE2:
915 : {
916 4906 : new_elem = std::make_unique<Quad4>();
917 9812 : new_elem->set_node(
918 4906 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
919 9812 : new_elem->set_node(
920 4906 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
921 9812 : new_elem->set_node(
922 4906 : 2, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
923 9812 : new_elem->set_node(
924 4906 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
925 :
926 4906 : if (elem->neighbor_ptr(0) == remote_elem)
927 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
928 4906 : if (elem->neighbor_ptr(1) == remote_elem)
929 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
930 :
931 4906 : break;
932 : }
933 0 : case EDGE3:
934 : {
935 0 : new_elem = std::make_unique<Quad9>();
936 0 : new_elem->set_node(
937 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
938 0 : new_elem->set_node(
939 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
940 0 : new_elem->set_node(
941 : 2,
942 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
943 0 : new_elem->set_node(
944 : 3,
945 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
946 0 : new_elem->set_node(
947 0 : 4, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
948 0 : new_elem->set_node(
949 : 5,
950 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
951 0 : new_elem->set_node(
952 : 6,
953 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
954 0 : new_elem->set_node(
955 : 7,
956 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
957 0 : new_elem->set_node(
958 : 8,
959 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
960 :
961 0 : if (elem->neighbor_ptr(0) == remote_elem)
962 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
963 0 : if (elem->neighbor_ptr(1) == remote_elem)
964 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
965 :
966 0 : break;
967 : }
968 936 : case TRI3:
969 : {
970 936 : new_elem = std::make_unique<Prism6>();
971 1872 : new_elem->set_node(
972 936 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
973 1872 : new_elem->set_node(
974 936 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
975 1872 : new_elem->set_node(
976 936 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
977 1872 : new_elem->set_node(
978 936 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
979 1872 : new_elem->set_node(
980 936 : 4, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
981 1872 : new_elem->set_node(
982 936 : 5, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
983 :
984 936 : if (elem->neighbor_ptr(0) == remote_elem)
985 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
986 936 : if (elem->neighbor_ptr(1) == remote_elem)
987 12 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
988 936 : if (elem->neighbor_ptr(2) == remote_elem)
989 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
990 :
991 936 : if (new_elem->volume() < 0.0)
992 : {
993 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
994 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
995 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
996 0 : is_flipped = true;
997 : }
998 :
999 936 : break;
1000 : }
1001 0 : case TRI6:
1002 : {
1003 0 : new_elem = std::make_unique<Prism18>();
1004 0 : new_elem->set_node(
1005 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
1006 0 : new_elem->set_node(
1007 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
1008 0 : new_elem->set_node(
1009 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
1010 0 : new_elem->set_node(
1011 : 3,
1012 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
1013 0 : new_elem->set_node(
1014 : 4,
1015 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
1016 0 : new_elem->set_node(
1017 : 5,
1018 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
1019 0 : new_elem->set_node(
1020 0 : 6, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
1021 0 : new_elem->set_node(
1022 0 : 7, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
1023 0 : new_elem->set_node(
1024 0 : 8, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
1025 0 : new_elem->set_node(
1026 : 9,
1027 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
1028 0 : new_elem->set_node(
1029 : 10,
1030 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
1031 0 : new_elem->set_node(
1032 : 11,
1033 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
1034 0 : new_elem->set_node(
1035 : 12,
1036 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
1037 0 : new_elem->set_node(
1038 : 13,
1039 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
1040 0 : new_elem->set_node(
1041 : 14,
1042 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
1043 0 : new_elem->set_node(
1044 : 15,
1045 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
1046 0 : new_elem->set_node(
1047 : 16,
1048 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
1049 0 : new_elem->set_node(
1050 : 17,
1051 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
1052 :
1053 0 : if (elem->neighbor_ptr(0) == remote_elem)
1054 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
1055 0 : if (elem->neighbor_ptr(1) == remote_elem)
1056 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
1057 0 : if (elem->neighbor_ptr(2) == remote_elem)
1058 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
1059 :
1060 0 : if (new_elem->volume() < 0.0)
1061 : {
1062 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
1063 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
1064 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
1065 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 12);
1066 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 7, 13);
1067 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 14);
1068 0 : is_flipped = true;
1069 : }
1070 :
1071 0 : break;
1072 : }
1073 0 : case TRI7:
1074 : {
1075 0 : new_elem = std::make_unique<Prism21>();
1076 0 : new_elem->set_node(
1077 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
1078 0 : new_elem->set_node(
1079 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
1080 0 : new_elem->set_node(
1081 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
1082 0 : new_elem->set_node(
1083 : 3,
1084 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
1085 0 : new_elem->set_node(
1086 : 4,
1087 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
1088 0 : new_elem->set_node(
1089 : 5,
1090 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
1091 0 : new_elem->set_node(
1092 0 : 6, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
1093 0 : new_elem->set_node(
1094 0 : 7, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
1095 0 : new_elem->set_node(
1096 0 : 8, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
1097 0 : new_elem->set_node(
1098 : 9,
1099 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
1100 0 : new_elem->set_node(
1101 : 10,
1102 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
1103 0 : new_elem->set_node(
1104 : 11,
1105 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
1106 0 : new_elem->set_node(
1107 : 12,
1108 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
1109 0 : new_elem->set_node(
1110 : 13,
1111 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
1112 0 : new_elem->set_node(
1113 : 14,
1114 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
1115 0 : new_elem->set_node(
1116 : 15,
1117 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
1118 0 : new_elem->set_node(
1119 : 16,
1120 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
1121 0 : new_elem->set_node(
1122 : 17,
1123 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
1124 0 : new_elem->set_node(
1125 0 : 18, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
1126 0 : new_elem->set_node(
1127 : 19,
1128 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
1129 0 : new_elem->set_node(
1130 : 20,
1131 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 1) * orig_nodes)));
1132 :
1133 0 : if (elem->neighbor_ptr(0) == remote_elem)
1134 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
1135 0 : if (elem->neighbor_ptr(1) == remote_elem)
1136 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
1137 0 : if (elem->neighbor_ptr(2) == remote_elem)
1138 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
1139 :
1140 0 : if (new_elem->volume() < 0.0)
1141 : {
1142 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
1143 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
1144 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
1145 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 12);
1146 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 7, 13);
1147 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 14);
1148 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 18, 19);
1149 0 : is_flipped = true;
1150 : }
1151 :
1152 0 : break;
1153 : }
1154 53500 : case QUAD4:
1155 : {
1156 53500 : new_elem = std::make_unique<Hex8>();
1157 107000 : new_elem->set_node(
1158 53500 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
1159 107000 : new_elem->set_node(
1160 53500 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
1161 107000 : new_elem->set_node(
1162 53500 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
1163 107000 : new_elem->set_node(
1164 53500 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (current_layer * orig_nodes)));
1165 107000 : new_elem->set_node(
1166 53500 : 4, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
1167 107000 : new_elem->set_node(
1168 53500 : 5, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
1169 107000 : new_elem->set_node(
1170 53500 : 6, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
1171 107000 : new_elem->set_node(
1172 53500 : 7, mesh->node_ptr(elem->node_ptr(3)->id() + ((current_layer + 1) * orig_nodes)));
1173 :
1174 53500 : if (elem->neighbor_ptr(0) == remote_elem)
1175 94 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
1176 53500 : if (elem->neighbor_ptr(1) == remote_elem)
1177 342 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
1178 53500 : if (elem->neighbor_ptr(2) == remote_elem)
1179 128 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
1180 53500 : if (elem->neighbor_ptr(3) == remote_elem)
1181 336 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
1182 :
1183 53500 : if (new_elem->volume() < 0.0)
1184 : {
1185 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
1186 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
1187 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
1188 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
1189 2024 : is_flipped = true;
1190 : }
1191 :
1192 53500 : break;
1193 : }
1194 16 : case QUAD8:
1195 : {
1196 16 : new_elem = std::make_unique<Hex20>();
1197 32 : new_elem->set_node(
1198 16 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
1199 32 : new_elem->set_node(
1200 16 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
1201 32 : new_elem->set_node(
1202 16 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
1203 32 : new_elem->set_node(
1204 16 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
1205 32 : new_elem->set_node(
1206 : 4,
1207 16 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
1208 32 : new_elem->set_node(
1209 : 5,
1210 16 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
1211 32 : new_elem->set_node(
1212 : 6,
1213 16 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
1214 32 : new_elem->set_node(
1215 : 7,
1216 16 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
1217 32 : new_elem->set_node(
1218 16 : 8, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
1219 32 : new_elem->set_node(
1220 16 : 9, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
1221 32 : new_elem->set_node(
1222 16 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
1223 32 : new_elem->set_node(
1224 16 : 11, mesh->node_ptr(elem->node_ptr(7)->id() + (2 * current_layer * orig_nodes)));
1225 32 : new_elem->set_node(
1226 : 12,
1227 16 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
1228 32 : new_elem->set_node(
1229 : 13,
1230 16 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
1231 32 : new_elem->set_node(
1232 : 14,
1233 16 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
1234 32 : new_elem->set_node(
1235 : 15,
1236 16 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
1237 32 : new_elem->set_node(
1238 : 16,
1239 16 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
1240 32 : new_elem->set_node(
1241 : 17,
1242 16 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
1243 32 : new_elem->set_node(
1244 : 18,
1245 16 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
1246 32 : new_elem->set_node(
1247 : 19,
1248 16 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 2) * orig_nodes)));
1249 :
1250 16 : if (elem->neighbor_ptr(0) == remote_elem)
1251 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
1252 16 : if (elem->neighbor_ptr(1) == remote_elem)
1253 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
1254 16 : if (elem->neighbor_ptr(2) == remote_elem)
1255 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
1256 16 : if (elem->neighbor_ptr(3) == remote_elem)
1257 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
1258 :
1259 16 : if (new_elem->volume() < 0.0)
1260 : {
1261 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
1262 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
1263 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
1264 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
1265 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
1266 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
1267 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
1268 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
1269 8 : is_flipped = true;
1270 : }
1271 :
1272 16 : break;
1273 : }
1274 0 : case QUAD9:
1275 : {
1276 0 : new_elem = std::make_unique<Hex27>();
1277 0 : new_elem->set_node(
1278 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
1279 0 : new_elem->set_node(
1280 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
1281 0 : new_elem->set_node(
1282 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
1283 0 : new_elem->set_node(
1284 0 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
1285 0 : new_elem->set_node(
1286 : 4,
1287 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
1288 0 : new_elem->set_node(
1289 : 5,
1290 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
1291 0 : new_elem->set_node(
1292 : 6,
1293 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
1294 0 : new_elem->set_node(
1295 : 7,
1296 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
1297 0 : new_elem->set_node(
1298 0 : 8, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
1299 0 : new_elem->set_node(
1300 0 : 9, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
1301 0 : new_elem->set_node(
1302 0 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
1303 0 : new_elem->set_node(
1304 0 : 11, mesh->node_ptr(elem->node_ptr(7)->id() + (2 * current_layer * orig_nodes)));
1305 0 : new_elem->set_node(
1306 : 12,
1307 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
1308 0 : new_elem->set_node(
1309 : 13,
1310 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
1311 0 : new_elem->set_node(
1312 : 14,
1313 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
1314 0 : new_elem->set_node(
1315 : 15,
1316 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
1317 0 : new_elem->set_node(
1318 : 16,
1319 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
1320 0 : new_elem->set_node(
1321 : 17,
1322 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
1323 0 : new_elem->set_node(
1324 : 18,
1325 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
1326 0 : new_elem->set_node(
1327 : 19,
1328 0 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 2) * orig_nodes)));
1329 0 : new_elem->set_node(
1330 0 : 20, mesh->node_ptr(elem->node_ptr(8)->id() + (2 * current_layer * orig_nodes)));
1331 0 : new_elem->set_node(
1332 : 21,
1333 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
1334 0 : new_elem->set_node(
1335 : 22,
1336 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
1337 0 : new_elem->set_node(
1338 : 23,
1339 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 1) * orig_nodes)));
1340 0 : new_elem->set_node(
1341 : 24,
1342 0 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 1) * orig_nodes)));
1343 0 : new_elem->set_node(
1344 : 25,
1345 0 : mesh->node_ptr(elem->node_ptr(8)->id() + ((2 * current_layer + 2) * orig_nodes)));
1346 0 : new_elem->set_node(
1347 : 26,
1348 0 : mesh->node_ptr(elem->node_ptr(8)->id() + ((2 * current_layer + 1) * orig_nodes)));
1349 :
1350 0 : if (elem->neighbor_ptr(0) == remote_elem)
1351 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
1352 0 : if (elem->neighbor_ptr(1) == remote_elem)
1353 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
1354 0 : if (elem->neighbor_ptr(2) == remote_elem)
1355 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
1356 0 : if (elem->neighbor_ptr(3) == remote_elem)
1357 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
1358 :
1359 0 : if (new_elem->volume() < 0.0)
1360 : {
1361 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
1362 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
1363 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
1364 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
1365 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
1366 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
1367 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
1368 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
1369 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 20, 25);
1370 0 : is_flipped = true;
1371 : }
1372 :
1373 0 : break;
1374 : }
1375 104 : case libMesh::C0POLYGON:
1376 : {
1377 104 : has_polygons = true;
1378 104 : const auto num_sides = elem->n_sides();
1379 104 : std::vector<std::shared_ptr<libMesh::Polygon>> sides;
1380 104 : sides.reserve(2 + num_sides);
1381 104 : if (2 * elem->n_nodes() > libMesh::C0Polygon::max_n_nodes)
1382 0 : mooseError("Too many nodes in polygons to extrude it. Max number of the prism "
1383 0 : "polyhedral nodes after extrusion: " +
1384 0 : std::to_string(libMesh::C0Polygon::max_n_nodes));
1385 : // Make a copy of the original element to use as a side
1386 104 : auto new_ptr = std::make_shared<libMesh::C0Polygon>(num_sides);
1387 552 : for (const auto node_i : make_range(elem->n_nodes()))
1388 : {
1389 : // This one will be oriented outwards, it should be OK though, polyhedron code does
1390 : // not mind
1391 896 : new_ptr->set_node(
1392 : node_i,
1393 448 : mesh->node_ptr(elem->node_ptr(node_i)->id() + (current_layer * orig_nodes)));
1394 : }
1395 104 : sides.push_back(new_ptr);
1396 : // Form the next horizontal side
1397 104 : auto translated_side = std::make_shared<libMesh::C0Polygon>(num_sides);
1398 552 : for (const auto node_i : make_range(elem->n_nodes()))
1399 896 : translated_side->set_node(node_i,
1400 448 : mesh->node_ptr(elem->node_ptr(node_i)->id() +
1401 448 : ((current_layer + 1) * orig_nodes)));
1402 104 : sides.push_back(translated_side);
1403 :
1404 : // Form the vertical sides
1405 552 : for (const auto side_i : make_range(num_sides))
1406 : {
1407 : // If the side already exists, use that
1408 : std::array<unsigned int, 3> side_key = {
1409 : current_layer,
1410 : static_cast<unsigned int>(
1411 448 : std::min(elem->node_ptr(side_i)->id(),
1412 448 : elem->node_ptr((side_i + 1) % num_sides)->id())),
1413 : static_cast<unsigned int>(
1414 448 : std::max(elem->node_ptr(side_i)->id(),
1415 896 : elem->node_ptr((side_i + 1) % num_sides)->id()))};
1416 448 : if (poly_extruded_sides.count(side_key))
1417 : {
1418 160 : sides.push_back(poly_extruded_sides[side_key]);
1419 160 : continue;
1420 : }
1421 :
1422 : // They are all quads, but constructor expects polygons
1423 288 : auto vert_side = std::make_shared<libMesh::C0Polygon>(4);
1424 576 : vert_side->set_node(
1425 288 : 0, mesh->node_ptr(elem->node_ptr(side_i)->id() + (current_layer * orig_nodes)));
1426 576 : vert_side->set_node(1,
1427 288 : mesh->node_ptr(elem->node_ptr((side_i + 1) % num_sides)->id() +
1428 288 : (current_layer * orig_nodes)));
1429 576 : vert_side->set_node(2,
1430 288 : mesh->node_ptr(elem->node_ptr((side_i + 1) % num_sides)->id() +
1431 288 : ((current_layer + 1) * orig_nodes)));
1432 576 : vert_side->set_node(3,
1433 288 : mesh->node_ptr(elem->node_ptr(side_i)->id() +
1434 288 : ((current_layer + 1) * orig_nodes)));
1435 288 : sides.push_back(vert_side);
1436 :
1437 288 : poly_extruded_sides.insert(std::make_pair(side_key, vert_side));
1438 288 : }
1439 : mooseAssert(sides.size() == 2 + num_sides, "Unexpected size of side vector");
1440 :
1441 : // Create the element from the sides, let libMesh figure out the orientation
1442 104 : std::unique_ptr<libMesh::Node> mid_elem_node;
1443 104 : new_elem = std::make_unique<libMesh::C0Polyhedron>(sides, mid_elem_node);
1444 104 : if (mid_elem_node)
1445 : {
1446 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1447 : // Number it at the end for convenience
1448 : // use the element ID to be able to set in parallel
1449 64 : unsigned int total_new_node_layers = total_num_layers * order;
1450 64 : unsigned int last_uid = orig_unique_ids + (total_new_node_layers - 1) * orig_elem +
1451 64 : total_new_node_layers * orig_nodes + elem->unique_id();
1452 64 : mid_elem_node->set_unique_id(last_uid);
1453 64 : has_poly_midnodes = true;
1454 : #endif
1455 64 : mesh->add_node(std::move(mid_elem_node));
1456 : }
1457 :
1458 104 : break;
1459 104 : }
1460 0 : default:
1461 0 : mooseError("Extrusion is not implemented for element type " + Moose::stringify(etype));
1462 : }
1463 :
1464 59462 : new_elem->set_id(elem->id() + (current_layer * orig_elem));
1465 59462 : new_elem->processor_id() = elem->processor_id();
1466 :
1467 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1468 : // Let's give the base of the extruded mesh the same
1469 : // unique_ids as the source mesh, in case anyone finds that
1470 : // a useful map to preserve.
1471 : const unique_id_type uid = (current_layer == 0)
1472 59462 : ? elem->unique_id()
1473 42182 : : orig_unique_ids +
1474 42182 : (current_layer - 1) * (orig_elem + orig_nodes) +
1475 42182 : orig_nodes + elem->id();
1476 :
1477 59462 : new_elem->set_unique_id(uid);
1478 : #endif
1479 :
1480 : // maintain the subdomain_id
1481 59462 : new_elem->subdomain_id() = elem->subdomain_id();
1482 :
1483 : // define upward boundaries
1484 59462 : if (k == num_layers - 1)
1485 : {
1486 : // Identify the side index of the new element that is part of the upward boundary
1487 : const unsigned short top_id =
1488 26628 : new_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1489 : // Assign sideset id to the side if the element belongs to a specified
1490 : // upward_boundary_source_blocks
1491 26628 : if (_upward_boundary_source_blocks.size() || _upward_boundary_ids.size())
1492 31604 : for (unsigned int i = 0; i < _upward_boundary_source_blocks[e].size(); i++)
1493 5328 : if (new_elem->subdomain_id() == _upward_boundary_source_blocks[e][i])
1494 3552 : boundary_info.add_side(
1495 3552 : new_elem.get(), is_flipped ? 0 : top_id, _upward_boundary_ids[e][i]);
1496 : }
1497 : // define downward boundaries
1498 59462 : if (k == 0)
1499 : {
1500 : const unsigned short top_id =
1501 26628 : new_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1502 26628 : if (_downward_boundary_source_blocks.size() || _downward_boundary_ids.size())
1503 31604 : for (unsigned int i = 0; i < _downward_boundary_source_blocks[e].size(); i++)
1504 5328 : if (new_elem->subdomain_id() == _downward_boundary_source_blocks[e][i])
1505 3552 : boundary_info.add_side(
1506 3552 : new_elem.get(), is_flipped ? top_id : 0, _downward_boundary_ids[e][i]);
1507 : }
1508 :
1509 : // perform subdomain swaps
1510 59462 : if (_subdomain_swap_pairs.size())
1511 : {
1512 15492 : auto & elevation_swap_pairs = _subdomain_swap_pairs[e];
1513 :
1514 15492 : auto new_id_it = elevation_swap_pairs.find(elem->subdomain_id());
1515 :
1516 15492 : if (new_id_it != elevation_swap_pairs.end())
1517 15492 : new_elem->subdomain_id() = new_id_it->second;
1518 : }
1519 :
1520 59462 : Elem * added_elem = mesh->add_elem(std::move(new_elem));
1521 :
1522 : // maintain extra integers
1523 72134 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
1524 12672 : added_elem->set_extra_integer(i, elem->get_extra_integer(i));
1525 :
1526 59462 : if (_elem_integers_swap_pairs.size())
1527 : {
1528 19008 : for (unsigned int i = 0; i < _elem_integer_indices_to_swap.size(); i++)
1529 : {
1530 12672 : auto & elevation_extra_swap_pairs = _elem_integers_swap_pairs[i * _heights.size() + e];
1531 :
1532 12672 : auto new_extra_id_it = elevation_extra_swap_pairs.find(
1533 12672 : elem->get_extra_integer(_elem_integer_indices_to_swap[i]));
1534 :
1535 12672 : if (new_extra_id_it != elevation_extra_swap_pairs.end())
1536 7392 : added_elem->set_extra_integer(_elem_integer_indices_to_swap[i],
1537 7392 : new_extra_id_it->second);
1538 : }
1539 : }
1540 :
1541 : // Copy any old boundary ids on all sides
1542 286594 : for (auto s : elem->side_index_range())
1543 : {
1544 227132 : input_boundary_info.boundary_ids(elem, s, ids_to_copy);
1545 :
1546 227132 : if (added_elem->dim() == 3)
1547 : {
1548 : // For 2D->3D extrusion, we give the boundary IDs
1549 : // for side s on the old element to side s+1 on the
1550 : // new element. This is just a happy coincidence as
1551 : // far as I can tell...
1552 217320 : if (_boundary_swap_pairs.empty())
1553 217320 : boundary_info.add_side(added_elem, cast_int<unsigned short>(s + 1), ids_to_copy);
1554 : else
1555 0 : for (const auto & id_to_copy : ids_to_copy)
1556 0 : boundary_info.add_side(added_elem,
1557 0 : cast_int<unsigned short>(s + 1),
1558 0 : _boundary_swap_pairs[e].count(id_to_copy)
1559 0 : ? _boundary_swap_pairs[e][id_to_copy]
1560 0 : : id_to_copy);
1561 : }
1562 : else
1563 : {
1564 : // For 1D->2D extrusion, the boundary IDs map as:
1565 : // Old elem -> New elem
1566 : // 0 -> 3
1567 : // 1 -> 1
1568 : libmesh_assert_less(s, 2);
1569 9812 : const unsigned short sidemap[2] = {3, 1};
1570 9812 : if (_boundary_swap_pairs.empty())
1571 9812 : boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
1572 : else
1573 0 : for (const auto & id_to_copy : ids_to_copy)
1574 0 : boundary_info.add_side(added_elem,
1575 0 : sidemap[s],
1576 0 : _boundary_swap_pairs[e].count(id_to_copy)
1577 0 : ? _boundary_swap_pairs[e][id_to_copy]
1578 0 : : id_to_copy);
1579 : }
1580 : }
1581 :
1582 : // Give new boundary ids to bottom and top
1583 59462 : if (current_layer == 0)
1584 : {
1585 : const unsigned short top_id =
1586 17280 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1587 17280 : if (_has_bottom_boundary)
1588 : {
1589 : mooseAssert(user_bottom_boundary_id != libMesh::BoundaryInfo::invalid_id,
1590 : "We should have retrieved a proper boundary ID");
1591 9678 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, user_bottom_boundary_id);
1592 : }
1593 : else
1594 7602 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, next_side_id);
1595 : }
1596 :
1597 59462 : if (current_layer == total_num_layers - 1)
1598 : {
1599 : // For 2D->3D extrusion, the "top" ID is 1+the original
1600 : // element's number of sides. For 1D->2D extrusion, the
1601 : // "top" ID is side 2.
1602 : const unsigned short top_id =
1603 17280 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1604 :
1605 17280 : if (_has_top_boundary)
1606 : {
1607 : mooseAssert(user_top_boundary_id != libMesh::BoundaryInfo::invalid_id,
1608 : "We should have retrieved a proper boundary ID");
1609 9678 : boundary_info.add_side(added_elem, is_flipped ? 0 : top_id, user_top_boundary_id);
1610 : }
1611 : else
1612 7602 : boundary_info.add_side(
1613 7602 : added_elem, is_flipped ? 0 : top_id, cast_int<boundary_id_type>(next_side_id + 1));
1614 : }
1615 :
1616 59462 : current_layer++;
1617 59462 : }
1618 : }
1619 371 : }
1620 :
1621 371 : if (has_polygons && !input->is_serial())
1622 0 : mooseError("Distributed meshes are not supported when extruding polygons at this time.");
1623 :
1624 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1625 : // Update the value of next_unique_id based on newly created nodes and elements
1626 : // Note: Number of element layers is one less than number of node layers
1627 371 : unsigned int total_new_node_layers = total_num_layers * order;
1628 371 : unsigned int new_unique_ids = orig_unique_ids + (total_new_node_layers - 1) * orig_elem +
1629 : total_new_node_layers * orig_nodes;
1630 : // Maximum case for the unique ids: all poly elements have a midnode
1631 371 : if (has_poly_midnodes)
1632 8 : new_unique_ids += orig_elem;
1633 371 : mesh->set_next_unique_id(new_unique_ids);
1634 : #endif
1635 :
1636 : // Copy all the subdomain/sideset/nodeset name maps to the extruded mesh
1637 371 : if (!input_subdomain_map.empty())
1638 42 : mesh->set_subdomain_name_map().insert(input_subdomain_map.begin(), input_subdomain_map.end());
1639 371 : if (!input_sideset_map.empty())
1640 279 : mesh->get_boundary_info().set_sideset_name_map().insert(input_sideset_map.begin(),
1641 : input_sideset_map.end());
1642 371 : if (!input_nodeset_map.empty())
1643 160 : mesh->get_boundary_info().set_nodeset_name_map().insert(input_nodeset_map.begin(),
1644 : input_nodeset_map.end());
1645 :
1646 371 : if (_has_bottom_boundary)
1647 232 : boundary_info.sideset_name(new_boundary_ids.front()) = new_boundary_names.front();
1648 371 : if (_has_top_boundary)
1649 232 : boundary_info.sideset_name(new_boundary_ids.back()) = new_boundary_names.back();
1650 :
1651 371 : mesh->unset_is_prepared();
1652 : // Creating the layered meshes creates a lot of leftover nodes, notably in the boundary_info,
1653 : // which will crash both paraview and trigger exodiff. Best to be safe.
1654 371 : if (extruding_quad_eights)
1655 8 : mesh->prepare_for_use();
1656 :
1657 742 : return mesh;
1658 371 : }
1659 :
1660 : Real
1661 12624 : AdvancedExtruderGenerator::radialExpansionRatio(const Real t) const
1662 : {
1663 : // NOTE: All functions added to this method must obey the following: f(0)=0, f(1)=1 for all t in
1664 : // [0,1].
1665 12624 : switch (_radial_expansion_method)
1666 : {
1667 7824 : case 0:
1668 7824 : return t;
1669 : // Only linear and cubic implemented
1670 4800 : default:
1671 4800 : return (-2. + _start_radial_growth_rate + _end_radial_growth_rate) * MathUtils::pow(t, 3) +
1672 9600 : (3. - (2. * _start_radial_growth_rate + _end_radial_growth_rate)) *
1673 4800 : MathUtils::pow(t, 2) +
1674 4800 : _start_radial_growth_rate * t;
1675 : }
1676 : }
|
