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 :
14 : #include "libmesh/boundary_info.h"
15 : #include "libmesh/function_base.h"
16 : #include "libmesh/cell_prism6.h"
17 : #include "libmesh/cell_prism18.h"
18 : #include "libmesh/cell_prism21.h"
19 : #include "libmesh/cell_hex8.h"
20 : #include "libmesh/cell_hex20.h"
21 : #include "libmesh/cell_hex27.h"
22 : #include "libmesh/edge_edge2.h"
23 : #include "libmesh/edge_edge3.h"
24 : #include "libmesh/edge_edge4.h"
25 : #include "libmesh/face_quad4.h"
26 : #include "libmesh/face_quad8.h"
27 : #include "libmesh/face_quad9.h"
28 : #include "libmesh/face_tri3.h"
29 : #include "libmesh/face_tri6.h"
30 : #include "libmesh/face_tri7.h"
31 : #include "libmesh/libmesh_logging.h"
32 : #include "libmesh/mesh_communication.h"
33 : #include "libmesh/mesh_modification.h"
34 : #include "libmesh/mesh_tools.h"
35 : #include "libmesh/parallel.h"
36 : #include "libmesh/remote_elem.h"
37 : #include "libmesh/string_to_enum.h"
38 : #include "libmesh/unstructured_mesh.h"
39 : #include "libmesh/point.h"
40 :
41 : #include <numeric>
42 :
43 : registerMooseObject("MooseApp", AdvancedExtruderGenerator);
44 : registerMooseObjectRenamed("MooseApp",
45 : FancyExtruderGenerator,
46 : "02/18/2023 24:00",
47 : AdvancedExtruderGenerator);
48 :
49 : InputParameters
50 28920 : AdvancedExtruderGenerator::validParams()
51 : {
52 28920 : InputParameters params = MeshGenerator::validParams();
53 :
54 28920 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh to extrude");
55 :
56 28920 : params.addClassDescription(
57 : "Extrudes a 1D mesh into 2D, or a 2D mesh into 3D, can have a variable height for each "
58 : "elevation, variable number of layers within each elevation, variable growth factors of "
59 : "axial element sizes within each elevation and remap subdomain_ids, boundary_ids and element "
60 : "extra integers within each elevation as well as interface boundaries between neighboring "
61 : "elevation layers.");
62 :
63 28920 : params.addRequiredParam<std::vector<Real>>("heights", "The height of each elevation");
64 :
65 28920 : params.addRangeCheckedParam<std::vector<Real>>(
66 : "biases", "biases>0.0", "The axial growth factor used for mesh biasing for each elevation.");
67 :
68 28920 : params.addRequiredParam<std::vector<unsigned int>>(
69 : "num_layers", "The number of layers for each elevation - must be num_elevations in length!");
70 :
71 28920 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
72 : "subdomain_swaps",
73 : {},
74 : "For each row, every two entries are interpreted as a pair of "
75 : "'from' and 'to' to remap the subdomains for that elevation");
76 :
77 28920 : params.addParam<std::vector<std::vector<boundary_id_type>>>(
78 : "boundary_swaps",
79 : {},
80 : "For each row, every two entries are interpreted as a pair of "
81 : "'from' and 'to' to remap the boundaries for that elevation");
82 :
83 28920 : params.addParam<std::vector<std::string>>(
84 : "elem_integer_names_to_swap",
85 : {},
86 : "Array of element extra integer names that need to be swapped during extrusion.");
87 :
88 28920 : params.addParam<std::vector<std::vector<std::vector<dof_id_type>>>>(
89 : "elem_integers_swaps",
90 : {},
91 : "For each row, every two entries are interpreted as a pair of 'from' and 'to' to remap the "
92 : "element extra integer for that elevation. If multiple element extra integers need to be "
93 : "swapped, the enties are stacked based on the order provided in "
94 : "'elem_integer_names_to_swap' to form the third dimension.");
95 :
96 28920 : params.addRequiredParam<Point>(
97 : "direction",
98 : "A vector that points in the direction to extrude (note, this will be "
99 : "normalized internally - so don't worry about it here)");
100 :
101 28920 : params.addParam<BoundaryName>(
102 : "top_boundary",
103 : "The boundary name to set on the top boundary. If omitted an ID will be generated.");
104 :
105 28920 : params.addParam<BoundaryName>(
106 : "bottom_boundary",
107 : "The boundary name to set on the bottom boundary. If omitted an ID will be generated.");
108 :
109 28920 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
110 : "upward_boundary_source_blocks", "Block ids used to generate upward interface boundaries.");
111 :
112 28920 : params.addParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids",
113 : "Upward interface boundary ids.");
114 :
115 28920 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
116 : "downward_boundary_source_blocks",
117 : "Block ids used to generate downward interface boundaries.");
118 :
119 28920 : params.addParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids",
120 : "Downward interface boundary ids.");
121 28920 : params.addParamNamesToGroup(
122 : "top_boundary bottom_boundary upward_boundary_source_blocks upward_boundary_ids "
123 : "downward_boundary_source_blocks downward_boundary_ids",
124 : "Boundary Assignment");
125 28920 : params.addParamNamesToGroup(
126 : "subdomain_swaps boundary_swaps elem_integer_names_to_swap elem_integers_swaps", "ID Swap");
127 86760 : params.addParam<Real>("twist_pitch",
128 57840 : 0,
129 : "Pitch for helicoidal extrusion around an axis going through the origin "
130 : "following the direction vector");
131 28920 : return params;
132 0 : }
133 :
134 197 : AdvancedExtruderGenerator::AdvancedExtruderGenerator(const InputParameters & parameters)
135 : : MeshGenerator(parameters),
136 197 : _input(getMesh("input")),
137 197 : _heights(getParam<std::vector<Real>>("heights")),
138 581 : _biases(isParamValid("biases") ? getParam<std::vector<Real>>("biases")
139 384 : : std::vector<Real>(_heights.size(), 1.0)),
140 197 : _num_layers(getParam<std::vector<unsigned int>>("num_layers")),
141 197 : _subdomain_swaps(getParam<std::vector<std::vector<subdomain_id_type>>>("subdomain_swaps")),
142 197 : _boundary_swaps(getParam<std::vector<std::vector<boundary_id_type>>>("boundary_swaps")),
143 197 : _elem_integer_names_to_swap(getParam<std::vector<std::string>>("elem_integer_names_to_swap")),
144 197 : _elem_integers_swaps(
145 197 : getParam<std::vector<std::vector<std::vector<dof_id_type>>>>("elem_integers_swaps")),
146 197 : _direction(getParam<Point>("direction")),
147 197 : _has_top_boundary(isParamValid("top_boundary")),
148 197 : _top_boundary(isParamValid("top_boundary") ? getParam<BoundaryName>("top_boundary") : "0"),
149 197 : _has_bottom_boundary(isParamValid("bottom_boundary")),
150 197 : _bottom_boundary(isParamValid("bottom_boundary") ? getParam<BoundaryName>("bottom_boundary")
151 : : "0"),
152 577 : _upward_boundary_source_blocks(
153 394 : isParamValid("upward_boundary_source_blocks")
154 197 : ? getParam<std::vector<std::vector<subdomain_id_type>>>("upward_boundary_source_blocks")
155 183 : : std::vector<std::vector<subdomain_id_type>>(_heights.size(),
156 380 : std::vector<subdomain_id_type>())),
157 577 : _upward_boundary_ids(
158 394 : isParamValid("upward_boundary_ids")
159 197 : ? getParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids")
160 183 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
161 380 : std::vector<boundary_id_type>())),
162 774 : _downward_boundary_source_blocks(isParamValid("downward_boundary_source_blocks")
163 197 : ? getParam<std::vector<std::vector<subdomain_id_type>>>(
164 : "downward_boundary_source_blocks")
165 : : std::vector<std::vector<subdomain_id_type>>(
166 380 : _heights.size(), std::vector<subdomain_id_type>())),
167 577 : _downward_boundary_ids(
168 394 : isParamValid("downward_boundary_ids")
169 197 : ? getParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids")
170 183 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
171 380 : std::vector<boundary_id_type>())),
172 788 : _twist_pitch(getParam<Real>("twist_pitch"))
173 : {
174 197 : if (!_direction.norm())
175 0 : paramError("direction", "Must have some length!");
176 :
177 : // Normalize it
178 197 : _direction /= _direction.norm();
179 :
180 197 : const auto num_elevations = _heights.size();
181 :
182 197 : if (_num_layers.size() != num_elevations)
183 0 : paramError("heights", "The length of 'heights' and 'num_layers' must be the same in ", name());
184 :
185 197 : if (_subdomain_swaps.size() && (_subdomain_swaps.size() != num_elevations))
186 0 : paramError("subdomain_swaps",
187 0 : "If specified, 'subdomain_swaps' (" + std::to_string(_subdomain_swaps.size()) +
188 0 : ") must be the same length as 'heights' (" + std::to_string(num_elevations) +
189 : ") in ",
190 0 : name());
191 :
192 : try
193 : {
194 197 : MooseMeshUtils::idSwapParametersProcessor(
195 197 : name(), "subdomain_swaps", _subdomain_swaps, _subdomain_swap_pairs);
196 : }
197 0 : catch (const MooseException & e)
198 : {
199 0 : paramError("subdomain_swaps", e.what());
200 0 : }
201 :
202 197 : if (_boundary_swaps.size() && (_boundary_swaps.size() != num_elevations))
203 0 : paramError("boundary_swaps",
204 0 : "If specified, 'boundary_swaps' (" + std::to_string(_boundary_swaps.size()) +
205 0 : ") must be the same length as 'heights' (" + std::to_string(num_elevations) +
206 : ") in ",
207 0 : name());
208 :
209 : try
210 : {
211 197 : MooseMeshUtils::idSwapParametersProcessor(
212 197 : name(), "boundary_swaps", _boundary_swaps, _boundary_swap_pairs);
213 : }
214 0 : catch (const MooseException & e)
215 : {
216 0 : paramError("boundary_swaps", e.what());
217 0 : }
218 :
219 209 : if (_elem_integers_swaps.size() &&
220 12 : _elem_integers_swaps.size() != _elem_integer_names_to_swap.size())
221 0 : paramError("elem_integers_swaps",
222 : "If specified, 'elem_integers_swaps' must have the same length as the length of "
223 : "'elem_integer_names_to_swap'.");
224 :
225 221 : for (const auto & unit_elem_integers_swaps : _elem_integers_swaps)
226 24 : if (unit_elem_integers_swaps.size() != num_elevations)
227 0 : paramError("elem_integers_swaps",
228 : "If specified, each element of 'elem_integers_swaps' must have the same length as "
229 : "the length of 'heights'.");
230 :
231 : try
232 : {
233 197 : MooseMeshUtils::extraElemIntegerSwapParametersProcessor(name(),
234 : num_elevations,
235 197 : _elem_integer_names_to_swap.size(),
236 : _elem_integers_swaps,
237 197 : _elem_integers_swap_pairs);
238 : }
239 0 : catch (const MooseException & e)
240 : {
241 0 : paramError("elem_integers_swaps", e.what());
242 0 : }
243 :
244 197 : bool has_negative_entry = false;
245 197 : bool has_positive_entry = false;
246 672 : for (const auto & h : _heights)
247 : {
248 475 : if (h > 0.0)
249 471 : has_positive_entry = true;
250 : else
251 4 : has_negative_entry = true;
252 : }
253 :
254 197 : if (has_negative_entry && has_positive_entry)
255 4 : paramError("heights", "Cannot have both positive and negative heights!");
256 193 : if (_biases.size() != _heights.size())
257 0 : paramError("biases", "Size of this parameter, if provided, must be the same as heights.");
258 :
259 386 : if (_upward_boundary_source_blocks.size() != _upward_boundary_ids.size() ||
260 193 : _upward_boundary_ids.size() != num_elevations)
261 0 : paramError("upward_boundary_ids",
262 0 : "This parameter must have the same length (" +
263 0 : std::to_string(_upward_boundary_ids.size()) +
264 0 : ") as upward_boundary_source_blocks (" +
265 0 : std::to_string(_upward_boundary_source_blocks.size()) + ") and heights (" +
266 0 : std::to_string(num_elevations) + ")");
267 656 : for (unsigned int i = 0; i < _upward_boundary_source_blocks.size(); i++)
268 463 : if (_upward_boundary_source_blocks[i].size() != _upward_boundary_ids[i].size())
269 0 : paramError("upward_boundary_ids",
270 : "Every element of this parameter must have the same length as the corresponding "
271 : "element of upward_boundary_source_blocks.");
272 :
273 386 : if (_downward_boundary_source_blocks.size() != _downward_boundary_ids.size() ||
274 193 : _downward_boundary_ids.size() != num_elevations)
275 0 : paramError("downward_boundary_ids",
276 0 : "This parameter must have the same length (" +
277 0 : std::to_string(_downward_boundary_ids.size()) +
278 0 : ") as downward_boundary_source_blocks (" +
279 0 : std::to_string(_downward_boundary_source_blocks.size()) + ") and heights (" +
280 0 : std::to_string(num_elevations) + ")");
281 656 : for (unsigned int i = 0; i < _downward_boundary_source_blocks.size(); i++)
282 463 : if (_downward_boundary_source_blocks[i].size() != _downward_boundary_ids[i].size())
283 0 : paramError("downward_boundary_ids",
284 : "Every element of this parameter must have the same length as the corresponding "
285 : "element of downward_boundary_source_blocks.");
286 193 : }
287 :
288 : std::unique_ptr<MeshBase>
289 183 : AdvancedExtruderGenerator::generate()
290 : {
291 : // Note: bulk of this code originally from libmesh mesh_modification.C
292 : // Original copyright: Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
293 : // Original license is LGPL so it can be used here.
294 :
295 183 : auto mesh = buildMeshBaseObject(_input->mesh_dimension() + 1);
296 183 : mesh->set_mesh_dimension(_input->mesh_dimension() + 1);
297 :
298 : // Check if the element integer names are existent in the input mesh.
299 205 : for (unsigned int i = 0; i < _elem_integer_names_to_swap.size(); i++)
300 22 : if (_input->has_elem_integer(_elem_integer_names_to_swap[i]))
301 44 : _elem_integer_indices_to_swap.push_back(
302 22 : _input->get_elem_integer_index(_elem_integer_names_to_swap[i]));
303 : else
304 0 : paramError("elem_integer_names_to_swap",
305 : "Element ",
306 : i + 1,
307 : " of 'elem_integer_names_to_swap' in is not a valid extra element integer of the "
308 : "input mesh.");
309 :
310 : // prepare for transferring extra element integers from original mesh to the extruded mesh.
311 183 : const unsigned int num_extra_elem_integers = _input->n_elem_integers();
312 183 : std::vector<std::string> id_names;
313 :
314 205 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
315 : {
316 22 : id_names.push_back(_input->get_elem_integer_name(i));
317 22 : if (!mesh->has_elem_integer(id_names[i]))
318 22 : mesh->add_elem_integer(id_names[i]);
319 : }
320 :
321 : // retrieve subdomain/sideset/nodeset name maps
322 183 : const auto & input_subdomain_map = _input->get_subdomain_name_map();
323 183 : const auto & input_sideset_map = _input->get_boundary_info().get_sideset_name_map();
324 183 : const auto & input_nodeset_map = _input->get_boundary_info().get_nodeset_name_map();
325 :
326 : // Check that the swaps source blocks are present in the mesh
327 447 : for (const auto & swap : _subdomain_swaps)
328 1216 : for (const auto i : index_range(swap))
329 952 : if (i % 2 == 0 && !MooseMeshUtils::hasSubdomainID(*_input, swap[i]))
330 4 : paramError("subdomain_swaps", "The block '", swap[i], "' was not found within the mesh");
331 :
332 : // Check that the swaps source boundaries are present in the mesh
333 203 : for (const auto & swap : _boundary_swaps)
334 156 : for (const auto i : index_range(swap))
335 132 : if (i % 2 == 0 && !MooseMeshUtils::hasBoundaryID(*_input, swap[i]))
336 4 : paramError("boundary_swaps", "The boundary '", swap[i], "' was not found within the mesh");
337 :
338 : // Check that the source blocks for layer top/bottom boundaries exist in the mesh
339 574 : for (const auto & layer_vec : _upward_boundary_source_blocks)
340 463 : for (const auto bid : layer_vec)
341 64 : if (!MooseMeshUtils::hasSubdomainID(*_input, bid))
342 4 : paramError(
343 : "upward_boundary_source_blocks", "The block '", bid, "' was not found within the mesh");
344 558 : for (const auto & layer_vec : _downward_boundary_source_blocks)
345 451 : for (const auto bid : layer_vec)
346 64 : if (!MooseMeshUtils::hasSubdomainID(*_input, bid))
347 4 : paramError("downward_boundary_source_blocks",
348 : "The block '",
349 : bid,
350 : "' was not found within the mesh");
351 :
352 167 : std::unique_ptr<MeshBase> input = std::move(_input);
353 :
354 : // If we're using a distributed mesh... then make sure we don't have any remote elements hanging
355 : // around
356 167 : if (!input->is_serial())
357 25 : mesh->delete_remote_elements();
358 :
359 167 : unsigned int total_num_layers = std::accumulate(_num_layers.begin(), _num_layers.end(), 0);
360 :
361 167 : auto total_num_elevations = _heights.size();
362 :
363 167 : dof_id_type orig_elem = input->n_elem();
364 167 : dof_id_type orig_nodes = input->n_nodes();
365 :
366 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
367 167 : unique_id_type orig_unique_ids = input->parallel_max_unique_id();
368 : #endif
369 :
370 167 : unsigned int order = 1;
371 :
372 167 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
373 167 : const BoundaryInfo & input_boundary_info = input->get_boundary_info();
374 :
375 : // Determine boundary IDs for the new user provided boundary names
376 167 : std::vector<BoundaryName> new_boundary_names;
377 167 : if (_has_bottom_boundary)
378 92 : new_boundary_names.push_back(_bottom_boundary);
379 167 : if (_has_top_boundary)
380 92 : new_boundary_names.push_back(_top_boundary);
381 : std::vector<boundary_id_type> new_boundary_ids =
382 167 : MooseMeshUtils::getBoundaryIDs(*input, new_boundary_names, true);
383 : const auto user_bottom_boundary_id =
384 167 : _has_bottom_boundary ? new_boundary_ids.front() : libMesh::BoundaryInfo::invalid_id;
385 : const auto user_top_boundary_id =
386 167 : _has_top_boundary ? new_boundary_ids.back() : libMesh::BoundaryInfo::invalid_id;
387 :
388 : // We know a priori how many elements we'll need
389 167 : mesh->reserve_elem(total_num_layers * orig_elem);
390 :
391 : // Look for higher order elements which introduce an extra layer
392 167 : std::set<ElemType> higher_orders = {EDGE3, EDGE4, TRI6, TRI7, QUAD8, QUAD9};
393 167 : bool extruding_quad_eights = false;
394 167 : std::vector<ElemType> types;
395 167 : MeshTools::elem_types(*input, types);
396 359 : for (const auto elem_type : types)
397 : {
398 192 : if (higher_orders.count(elem_type))
399 8 : order = 2;
400 192 : if (elem_type == QUAD8)
401 8 : extruding_quad_eights = true;
402 : }
403 167 : mesh->comm().max(order);
404 167 : mesh->comm().max(extruding_quad_eights);
405 :
406 : // Reserve for the max number possibly needed
407 167 : mesh->reserve_nodes((order * total_num_layers + 1) * orig_nodes);
408 :
409 : // Container to catch the boundary IDs handed back by the BoundaryInfo object
410 167 : std::vector<boundary_id_type> ids_to_copy;
411 :
412 167 : Point old_distance;
413 167 : Point current_distance;
414 :
415 : // Create translated layers of nodes in the direction of extrusion
416 11276 : for (const auto & node : input->node_ptr_range())
417 : {
418 11109 : unsigned int current_node_layer = 0;
419 :
420 11109 : old_distance.zero();
421 :
422 : // e is the elevation layer ordering
423 33214 : for (unsigned int e = 0; e < total_num_elevations; e++)
424 : {
425 22105 : auto num_layers = _num_layers[e];
426 :
427 22105 : auto height = _heights[e];
428 :
429 22105 : auto bias = _biases[e];
430 :
431 : // k is the element layer ordering within each elevation layer
432 71683 : for (unsigned int k = 0; k < order * num_layers + (e == 0 ? 1 : 0); ++k)
433 : {
434 : // For the first layer we don't need to move
435 49578 : if (e == 0 && k == 0)
436 11109 : current_distance.zero();
437 : else
438 : {
439 : // Shift the previous position by a certain fraction of 'height' along the extrusion
440 : // direction to get the new position.
441 38469 : auto layer_index = (k - (e == 0 ? 1 : 0)) / order + 1;
442 :
443 38469 : const auto step_size = MooseUtils::absoluteFuzzyEqual(bias, 1.0)
444 38469 : ? height / (Real)num_layers / (Real)order
445 2310 : : height * std::pow(bias, (Real)(layer_index - 1)) *
446 2310 : (1.0 - bias) /
447 38469 : (1.0 - std::pow(bias, (Real)(num_layers))) / (Real)order;
448 :
449 38469 : current_distance = old_distance + _direction * step_size;
450 :
451 : // Handle helicoidal extrusion
452 38469 : if (!MooseUtils::absoluteFuzzyEqual(_twist_pitch, 0.))
453 : {
454 : // twist 1 should be 'normal' to the extruded shape
455 2772 : RealVectorValue twist1 = _direction.cross(*node);
456 : // This happens for any node on the helicoidal extrusion axis
457 2772 : if (!MooseUtils::absoluteFuzzyEqual(twist1.norm(), .0))
458 2718 : twist1 /= twist1.norm();
459 2772 : const RealVectorValue twist2 = twist1.cross(_direction);
460 :
461 5544 : auto twist = (cos(2. * libMesh::pi * layer_index * step_size / _twist_pitch) -
462 5544 : cos(2. * libMesh::pi * (layer_index - 1) * step_size / _twist_pitch)) *
463 2772 : twist2 +
464 5544 : (sin(2. * libMesh::pi * layer_index * step_size / _twist_pitch) -
465 5544 : sin(2. * libMesh::pi * (layer_index - 1) * step_size / _twist_pitch)) *
466 5544 : twist1;
467 2772 : twist *= std::sqrt(node->norm_sq() + libMesh::Utility::pow<2>(_direction * (*node)));
468 2772 : current_distance += twist;
469 : }
470 : }
471 :
472 99156 : Node * new_node = mesh->add_point(*node + current_distance,
473 49578 : node->id() + (current_node_layer * orig_nodes),
474 49578 : node->processor_id());
475 :
476 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
477 : // Let's give the base of the extruded mesh the same
478 : // unique_ids as the source mesh, in case anyone finds that
479 : // a useful map to preserve.
480 : const unique_id_type uid = (current_node_layer == 0)
481 49578 : ? node->unique_id()
482 38469 : : orig_unique_ids +
483 38469 : (current_node_layer - 1) * (orig_nodes + orig_elem) +
484 38469 : node->id();
485 :
486 49578 : new_node->set_unique_id(uid);
487 : #endif
488 :
489 49578 : input_boundary_info.boundary_ids(node, ids_to_copy);
490 49578 : if (_boundary_swap_pairs.empty())
491 49578 : boundary_info.add_node(new_node, ids_to_copy);
492 : else
493 0 : for (const auto & id_to_copy : ids_to_copy)
494 : {
495 0 : boundary_info.add_node(new_node,
496 0 : _boundary_swap_pairs[e].count(id_to_copy)
497 0 : ? _boundary_swap_pairs[e][id_to_copy]
498 0 : : id_to_copy);
499 : }
500 :
501 49578 : old_distance = current_distance;
502 49578 : current_node_layer++;
503 : }
504 : }
505 167 : }
506 :
507 167 : const auto & side_ids = input_boundary_info.get_side_boundary_ids();
508 :
509 : boundary_id_type next_side_id =
510 167 : side_ids.empty() ? 0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
511 :
512 : // side_ids may not include ids from remote elements, in which case
513 : // some processors may have underestimated the next_side_id; let's
514 : // fix that.
515 167 : input->comm().max(next_side_id);
516 :
517 9357 : for (const auto & elem : input->element_ptr_range())
518 : {
519 9190 : const ElemType etype = elem->type();
520 :
521 : // build_extrusion currently only works on coarse meshes
522 : libmesh_assert(!elem->parent());
523 :
524 9190 : unsigned int current_layer = 0;
525 :
526 27100 : for (unsigned int e = 0; e != total_num_elevations; e++)
527 : {
528 17910 : auto num_layers = _num_layers[e];
529 :
530 48576 : for (unsigned int k = 0; k != num_layers; ++k)
531 : {
532 30666 : std::unique_ptr<Elem> new_elem;
533 30666 : bool is_flipped(false);
534 30666 : switch (etype)
535 : {
536 10 : case EDGE2:
537 : {
538 10 : new_elem = std::make_unique<Quad4>();
539 20 : new_elem->set_node(
540 10 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
541 20 : new_elem->set_node(
542 10 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
543 20 : new_elem->set_node(
544 10 : 2, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
545 20 : new_elem->set_node(
546 10 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
547 :
548 10 : if (elem->neighbor_ptr(0) == remote_elem)
549 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
550 10 : if (elem->neighbor_ptr(1) == remote_elem)
551 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
552 :
553 10 : break;
554 : }
555 0 : case EDGE3:
556 : {
557 0 : new_elem = std::make_unique<Quad9>();
558 0 : new_elem->set_node(
559 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
560 0 : new_elem->set_node(
561 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
562 0 : new_elem->set_node(
563 : 2,
564 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
565 0 : new_elem->set_node(
566 : 3,
567 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
568 0 : new_elem->set_node(
569 0 : 4, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
570 0 : new_elem->set_node(
571 : 5,
572 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
573 0 : new_elem->set_node(
574 : 6,
575 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
576 0 : new_elem->set_node(
577 : 7,
578 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
579 0 : new_elem->set_node(
580 : 8,
581 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
582 :
583 0 : if (elem->neighbor_ptr(0) == remote_elem)
584 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
585 0 : if (elem->neighbor_ptr(1) == remote_elem)
586 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
587 :
588 0 : break;
589 : }
590 936 : case TRI3:
591 : {
592 936 : new_elem = std::make_unique<Prism6>();
593 1872 : new_elem->set_node(
594 936 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
595 1872 : new_elem->set_node(
596 936 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
597 1872 : new_elem->set_node(
598 936 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
599 1872 : new_elem->set_node(
600 936 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
601 1872 : new_elem->set_node(
602 936 : 4, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
603 1872 : new_elem->set_node(
604 936 : 5, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
605 :
606 936 : if (elem->neighbor_ptr(0) == remote_elem)
607 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
608 936 : if (elem->neighbor_ptr(1) == remote_elem)
609 12 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
610 936 : if (elem->neighbor_ptr(2) == remote_elem)
611 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
612 :
613 936 : if (new_elem->volume() < 0.0)
614 : {
615 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
616 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
617 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
618 0 : is_flipped = true;
619 : }
620 :
621 936 : break;
622 : }
623 0 : case TRI6:
624 : {
625 0 : new_elem = std::make_unique<Prism18>();
626 0 : new_elem->set_node(
627 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
628 0 : new_elem->set_node(
629 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
630 0 : new_elem->set_node(
631 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
632 0 : new_elem->set_node(
633 : 3,
634 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
635 0 : new_elem->set_node(
636 : 4,
637 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
638 0 : new_elem->set_node(
639 : 5,
640 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
641 0 : new_elem->set_node(
642 0 : 6, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
643 0 : new_elem->set_node(
644 0 : 7, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
645 0 : new_elem->set_node(
646 0 : 8, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
647 0 : new_elem->set_node(
648 : 9,
649 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
650 0 : new_elem->set_node(
651 : 10,
652 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
653 0 : new_elem->set_node(
654 : 11,
655 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
656 0 : new_elem->set_node(
657 : 12,
658 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
659 0 : new_elem->set_node(
660 : 13,
661 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
662 0 : new_elem->set_node(
663 : 14,
664 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
665 0 : new_elem->set_node(
666 : 15,
667 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
668 0 : new_elem->set_node(
669 : 16,
670 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
671 0 : new_elem->set_node(
672 : 17,
673 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
674 :
675 0 : if (elem->neighbor_ptr(0) == remote_elem)
676 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
677 0 : if (elem->neighbor_ptr(1) == remote_elem)
678 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
679 0 : if (elem->neighbor_ptr(2) == remote_elem)
680 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
681 :
682 0 : if (new_elem->volume() < 0.0)
683 : {
684 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
685 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
686 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
687 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 12);
688 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 7, 13);
689 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 14);
690 0 : is_flipped = true;
691 : }
692 :
693 0 : break;
694 : }
695 0 : case TRI7:
696 : {
697 0 : new_elem = std::make_unique<Prism21>();
698 0 : new_elem->set_node(
699 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
700 0 : new_elem->set_node(
701 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
702 0 : new_elem->set_node(
703 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
704 0 : new_elem->set_node(
705 : 3,
706 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
707 0 : new_elem->set_node(
708 : 4,
709 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
710 0 : new_elem->set_node(
711 : 5,
712 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
713 0 : new_elem->set_node(
714 0 : 6, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
715 0 : new_elem->set_node(
716 0 : 7, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
717 0 : new_elem->set_node(
718 0 : 8, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
719 0 : new_elem->set_node(
720 : 9,
721 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
722 0 : new_elem->set_node(
723 : 10,
724 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
725 0 : new_elem->set_node(
726 : 11,
727 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
728 0 : new_elem->set_node(
729 : 12,
730 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
731 0 : new_elem->set_node(
732 : 13,
733 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
734 0 : new_elem->set_node(
735 : 14,
736 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
737 0 : new_elem->set_node(
738 : 15,
739 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
740 0 : new_elem->set_node(
741 : 16,
742 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
743 0 : new_elem->set_node(
744 : 17,
745 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
746 0 : new_elem->set_node(
747 0 : 18, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
748 0 : new_elem->set_node(
749 : 19,
750 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
751 0 : new_elem->set_node(
752 : 20,
753 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 1) * orig_nodes)));
754 :
755 0 : if (elem->neighbor_ptr(0) == remote_elem)
756 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
757 0 : if (elem->neighbor_ptr(1) == remote_elem)
758 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
759 0 : if (elem->neighbor_ptr(2) == remote_elem)
760 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
761 :
762 0 : if (new_elem->volume() < 0.0)
763 : {
764 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
765 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
766 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
767 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 12);
768 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 7, 13);
769 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 14);
770 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 18, 19);
771 0 : is_flipped = true;
772 : }
773 :
774 0 : break;
775 : }
776 29704 : case QUAD4:
777 : {
778 29704 : new_elem = std::make_unique<Hex8>();
779 59408 : new_elem->set_node(
780 29704 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
781 59408 : new_elem->set_node(
782 29704 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
783 59408 : new_elem->set_node(
784 29704 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
785 59408 : new_elem->set_node(
786 29704 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (current_layer * orig_nodes)));
787 59408 : new_elem->set_node(
788 29704 : 4, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
789 59408 : new_elem->set_node(
790 29704 : 5, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
791 59408 : new_elem->set_node(
792 29704 : 6, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
793 59408 : new_elem->set_node(
794 29704 : 7, mesh->node_ptr(elem->node_ptr(3)->id() + ((current_layer + 1) * orig_nodes)));
795 :
796 29704 : if (elem->neighbor_ptr(0) == remote_elem)
797 94 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
798 29704 : if (elem->neighbor_ptr(1) == remote_elem)
799 318 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
800 29704 : if (elem->neighbor_ptr(2) == remote_elem)
801 110 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
802 29704 : if (elem->neighbor_ptr(3) == remote_elem)
803 306 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
804 :
805 29704 : if (new_elem->volume() < 0.0)
806 : {
807 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
808 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
809 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
810 2024 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
811 2024 : is_flipped = true;
812 : }
813 :
814 29704 : break;
815 : }
816 16 : case QUAD8:
817 : {
818 16 : new_elem = std::make_unique<Hex20>();
819 32 : new_elem->set_node(
820 16 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
821 32 : new_elem->set_node(
822 16 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
823 32 : new_elem->set_node(
824 16 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
825 32 : new_elem->set_node(
826 16 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
827 32 : new_elem->set_node(
828 : 4,
829 16 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
830 32 : new_elem->set_node(
831 : 5,
832 16 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
833 32 : new_elem->set_node(
834 : 6,
835 16 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
836 32 : new_elem->set_node(
837 : 7,
838 16 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
839 32 : new_elem->set_node(
840 16 : 8, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
841 32 : new_elem->set_node(
842 16 : 9, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
843 32 : new_elem->set_node(
844 16 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
845 32 : new_elem->set_node(
846 16 : 11, mesh->node_ptr(elem->node_ptr(7)->id() + (2 * current_layer * orig_nodes)));
847 32 : new_elem->set_node(
848 : 12,
849 16 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
850 32 : new_elem->set_node(
851 : 13,
852 16 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
853 32 : new_elem->set_node(
854 : 14,
855 16 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
856 32 : new_elem->set_node(
857 : 15,
858 16 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
859 32 : new_elem->set_node(
860 : 16,
861 16 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
862 32 : new_elem->set_node(
863 : 17,
864 16 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
865 32 : new_elem->set_node(
866 : 18,
867 16 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
868 32 : new_elem->set_node(
869 : 19,
870 16 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 2) * orig_nodes)));
871 :
872 16 : if (elem->neighbor_ptr(0) == remote_elem)
873 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
874 16 : if (elem->neighbor_ptr(1) == remote_elem)
875 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
876 16 : if (elem->neighbor_ptr(2) == remote_elem)
877 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
878 16 : if (elem->neighbor_ptr(3) == remote_elem)
879 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
880 :
881 16 : if (new_elem->volume() < 0.0)
882 : {
883 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
884 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
885 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
886 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
887 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
888 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
889 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
890 8 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
891 8 : is_flipped = true;
892 : }
893 :
894 16 : break;
895 : }
896 0 : case QUAD9:
897 : {
898 0 : new_elem = std::make_unique<Hex27>();
899 0 : new_elem->set_node(
900 0 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
901 0 : new_elem->set_node(
902 0 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
903 0 : new_elem->set_node(
904 0 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
905 0 : new_elem->set_node(
906 0 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
907 0 : new_elem->set_node(
908 : 4,
909 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
910 0 : new_elem->set_node(
911 : 5,
912 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
913 0 : new_elem->set_node(
914 : 6,
915 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
916 0 : new_elem->set_node(
917 : 7,
918 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
919 0 : new_elem->set_node(
920 0 : 8, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
921 0 : new_elem->set_node(
922 0 : 9, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
923 0 : new_elem->set_node(
924 0 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
925 0 : new_elem->set_node(
926 0 : 11, mesh->node_ptr(elem->node_ptr(7)->id() + (2 * current_layer * orig_nodes)));
927 0 : new_elem->set_node(
928 : 12,
929 0 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
930 0 : new_elem->set_node(
931 : 13,
932 0 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
933 0 : new_elem->set_node(
934 : 14,
935 0 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
936 0 : new_elem->set_node(
937 : 15,
938 0 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
939 0 : new_elem->set_node(
940 : 16,
941 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
942 0 : new_elem->set_node(
943 : 17,
944 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
945 0 : new_elem->set_node(
946 : 18,
947 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
948 0 : new_elem->set_node(
949 : 19,
950 0 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 2) * orig_nodes)));
951 0 : new_elem->set_node(
952 0 : 20, mesh->node_ptr(elem->node_ptr(8)->id() + (2 * current_layer * orig_nodes)));
953 0 : new_elem->set_node(
954 : 21,
955 0 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 1) * orig_nodes)));
956 0 : new_elem->set_node(
957 : 22,
958 0 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 1) * orig_nodes)));
959 0 : new_elem->set_node(
960 : 23,
961 0 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 1) * orig_nodes)));
962 0 : new_elem->set_node(
963 : 24,
964 0 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 1) * orig_nodes)));
965 0 : new_elem->set_node(
966 : 25,
967 0 : mesh->node_ptr(elem->node_ptr(8)->id() + ((2 * current_layer + 2) * orig_nodes)));
968 0 : new_elem->set_node(
969 : 26,
970 0 : mesh->node_ptr(elem->node_ptr(8)->id() + ((2 * current_layer + 1) * orig_nodes)));
971 :
972 0 : if (elem->neighbor_ptr(0) == remote_elem)
973 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
974 0 : if (elem->neighbor_ptr(1) == remote_elem)
975 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
976 0 : if (elem->neighbor_ptr(2) == remote_elem)
977 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
978 0 : if (elem->neighbor_ptr(3) == remote_elem)
979 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
980 :
981 0 : if (new_elem->volume() < 0.0)
982 : {
983 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
984 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
985 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
986 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
987 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
988 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
989 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
990 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
991 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 20, 25);
992 0 : is_flipped = true;
993 : }
994 :
995 0 : break;
996 : }
997 0 : default:
998 0 : mooseError("Extrusion is not implemented for element type " + Moose::stringify(etype));
999 : }
1000 :
1001 30666 : new_elem->set_id(elem->id() + (current_layer * orig_elem));
1002 30666 : new_elem->processor_id() = elem->processor_id();
1003 :
1004 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1005 : // Let's give the base of the extruded mesh the same
1006 : // unique_ids as the source mesh, in case anyone finds that
1007 : // a useful map to preserve.
1008 : const unique_id_type uid = (current_layer == 0)
1009 30666 : ? elem->unique_id()
1010 21476 : : orig_unique_ids +
1011 21476 : (current_layer - 1) * (orig_nodes + orig_elem) +
1012 21476 : orig_nodes + elem->id();
1013 :
1014 30666 : new_elem->set_unique_id(uid);
1015 : #endif
1016 :
1017 : // maintain the subdomain_id
1018 30666 : new_elem->subdomain_id() = elem->subdomain_id();
1019 :
1020 : // define upward boundaries
1021 30666 : if (k == num_layers - 1)
1022 : {
1023 : // Identify the side index of the new element that is part of the upward boundary
1024 : const unsigned short top_id =
1025 17910 : new_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1026 : // Assign sideset id to the side if the element belongs to a specified
1027 : // upward_boundary_source_blocks
1028 23238 : for (unsigned int i = 0; i < _upward_boundary_source_blocks[e].size(); i++)
1029 5328 : if (new_elem->subdomain_id() == _upward_boundary_source_blocks[e][i])
1030 3552 : boundary_info.add_side(
1031 3552 : new_elem.get(), is_flipped ? 0 : top_id, _upward_boundary_ids[e][i]);
1032 : }
1033 : // define downward boundaries
1034 30666 : if (k == 0)
1035 : {
1036 : const unsigned short top_id =
1037 17910 : new_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1038 23238 : for (unsigned int i = 0; i < _downward_boundary_source_blocks[e].size(); i++)
1039 5328 : if (new_elem->subdomain_id() == _downward_boundary_source_blocks[e][i])
1040 3552 : boundary_info.add_side(
1041 3552 : new_elem.get(), is_flipped ? top_id : 0, _downward_boundary_ids[e][i]);
1042 : }
1043 :
1044 : // perform subdomain swaps
1045 30666 : if (_subdomain_swap_pairs.size())
1046 : {
1047 13476 : auto & elevation_swap_pairs = _subdomain_swap_pairs[e];
1048 :
1049 13476 : auto new_id_it = elevation_swap_pairs.find(elem->subdomain_id());
1050 :
1051 13476 : if (new_id_it != elevation_swap_pairs.end())
1052 13476 : new_elem->subdomain_id() = new_id_it->second;
1053 : }
1054 :
1055 30666 : Elem * added_elem = mesh->add_elem(std::move(new_elem));
1056 :
1057 : // maintain extra integers
1058 43338 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
1059 12672 : added_elem->set_extra_integer(i, elem->get_extra_integer(i));
1060 :
1061 30666 : if (_elem_integers_swap_pairs.size())
1062 : {
1063 19008 : for (unsigned int i = 0; i < _elem_integer_indices_to_swap.size(); i++)
1064 : {
1065 12672 : auto & elevation_extra_swap_pairs = _elem_integers_swap_pairs[i * _heights.size() + e];
1066 :
1067 12672 : auto new_extra_id_it = elevation_extra_swap_pairs.find(
1068 12672 : elem->get_extra_integer(_elem_integer_indices_to_swap[i]));
1069 :
1070 12672 : if (new_extra_id_it != elevation_extra_swap_pairs.end())
1071 7392 : added_elem->set_extra_integer(_elem_integer_indices_to_swap[i],
1072 7392 : new_extra_id_it->second);
1073 : }
1074 : }
1075 :
1076 : // Copy any old boundary ids on all sides
1077 152374 : for (auto s : elem->side_index_range())
1078 : {
1079 121708 : input_boundary_info.boundary_ids(elem, s, ids_to_copy);
1080 :
1081 121708 : if (added_elem->dim() == 3)
1082 : {
1083 : // For 2D->3D extrusion, we give the boundary IDs
1084 : // for side s on the old element to side s+1 on the
1085 : // new element. This is just a happy coincidence as
1086 : // far as I can tell...
1087 121688 : if (_boundary_swap_pairs.empty())
1088 121688 : boundary_info.add_side(added_elem, cast_int<unsigned short>(s + 1), ids_to_copy);
1089 : else
1090 0 : for (const auto & id_to_copy : ids_to_copy)
1091 0 : boundary_info.add_side(added_elem,
1092 0 : cast_int<unsigned short>(s + 1),
1093 0 : _boundary_swap_pairs[e].count(id_to_copy)
1094 0 : ? _boundary_swap_pairs[e][id_to_copy]
1095 0 : : id_to_copy);
1096 : }
1097 : else
1098 : {
1099 : // For 1D->2D extrusion, the boundary IDs map as:
1100 : // Old elem -> New elem
1101 : // 0 -> 3
1102 : // 1 -> 1
1103 : libmesh_assert_less(s, 2);
1104 20 : const unsigned short sidemap[2] = {3, 1};
1105 20 : if (_boundary_swap_pairs.empty())
1106 20 : boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
1107 : else
1108 0 : for (const auto & id_to_copy : ids_to_copy)
1109 0 : boundary_info.add_side(added_elem,
1110 0 : sidemap[s],
1111 0 : _boundary_swap_pairs[e].count(id_to_copy)
1112 0 : ? _boundary_swap_pairs[e][id_to_copy]
1113 0 : : id_to_copy);
1114 : }
1115 : }
1116 :
1117 : // Give new boundary ids to bottom and top
1118 30666 : if (current_layer == 0)
1119 : {
1120 : const unsigned short top_id =
1121 9190 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1122 9190 : if (_has_bottom_boundary)
1123 : {
1124 : mooseAssert(user_bottom_boundary_id != libMesh::BoundaryInfo::invalid_id,
1125 : "We should have retrieved a proper boundary ID");
1126 2574 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, user_bottom_boundary_id);
1127 : }
1128 : else
1129 6616 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, next_side_id);
1130 : }
1131 :
1132 30666 : if (current_layer == total_num_layers - 1)
1133 : {
1134 : // For 2D->3D extrusion, the "top" ID is 1+the original
1135 : // element's number of sides. For 1D->2D extrusion, the
1136 : // "top" ID is side 2.
1137 : const unsigned short top_id =
1138 9190 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1139 :
1140 9190 : if (_has_top_boundary)
1141 : {
1142 : mooseAssert(user_top_boundary_id != libMesh::BoundaryInfo::invalid_id,
1143 : "We should have retrieved a proper boundary ID");
1144 2574 : boundary_info.add_side(added_elem, is_flipped ? 0 : top_id, user_top_boundary_id);
1145 : }
1146 : else
1147 6616 : boundary_info.add_side(
1148 6616 : added_elem, is_flipped ? 0 : top_id, cast_int<boundary_id_type>(next_side_id + 1));
1149 : }
1150 :
1151 30666 : current_layer++;
1152 30666 : }
1153 : }
1154 167 : }
1155 :
1156 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1157 : // Update the value of next_unique_id based on newly created nodes and elements
1158 : // Note: Number of element layers is one less than number of node layers
1159 167 : unsigned int total_new_node_layers = total_num_layers * order;
1160 167 : unsigned int new_unique_ids = orig_unique_ids + (total_new_node_layers - 1) * orig_elem +
1161 : total_new_node_layers * orig_nodes;
1162 167 : mesh->set_next_unique_id(new_unique_ids);
1163 : #endif
1164 :
1165 : // Copy all the subdomain/sideset/nodeset name maps to the extruded mesh
1166 167 : if (!input_subdomain_map.empty())
1167 26 : mesh->set_subdomain_name_map().insert(input_subdomain_map.begin(), input_subdomain_map.end());
1168 167 : if (!input_sideset_map.empty())
1169 167 : mesh->get_boundary_info().set_sideset_name_map().insert(input_sideset_map.begin(),
1170 : input_sideset_map.end());
1171 167 : if (!input_nodeset_map.empty())
1172 110 : mesh->get_boundary_info().set_nodeset_name_map().insert(input_nodeset_map.begin(),
1173 : input_nodeset_map.end());
1174 :
1175 167 : if (_has_bottom_boundary)
1176 92 : boundary_info.sideset_name(new_boundary_ids.front()) = new_boundary_names.front();
1177 167 : if (_has_top_boundary)
1178 92 : boundary_info.sideset_name(new_boundary_ids.back()) = new_boundary_names.back();
1179 :
1180 167 : mesh->set_isnt_prepared();
1181 : // Creating the layered meshes creates a lot of leftover nodes, notably in the boundary_info,
1182 : // which will crash both paraview and trigger exodiff. Best to be safe.
1183 167 : if (extruding_quad_eights)
1184 8 : mesh->prepare_for_use();
1185 :
1186 334 : return mesh;
1187 167 : }
|
