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