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 28954 : AdvancedExtruderGenerator::validParams()
51 : {
52 28954 : InputParameters params = MeshGenerator::validParams();
53 :
54 28954 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh to extrude");
55 :
56 28954 : 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 28954 : params.addRequiredParam<std::vector<Real>>("heights", "The height of each elevation");
64 :
65 28954 : params.addRangeCheckedParam<std::vector<Real>>(
66 : "biases", "biases>0.0", "The axial growth factor used for mesh biasing for each elevation.");
67 :
68 28954 : params.addRequiredParam<std::vector<unsigned int>>(
69 : "num_layers", "The number of layers for each elevation - must be num_elevations in length!");
70 :
71 28954 : 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 28954 : 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 28954 : 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 28954 : 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 28954 : 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 28954 : 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 28954 : 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 28954 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
110 : "upward_boundary_source_blocks", "Block ids used to generate upward interface boundaries.");
111 :
112 28954 : params.addParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids",
113 : "Upward interface boundary ids.");
114 :
115 28954 : 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 28954 : params.addParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids",
120 : "Downward interface boundary ids.");
121 28954 : 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 28954 : params.addParamNamesToGroup(
126 : "subdomain_swaps boundary_swaps elem_integer_names_to_swap elem_integers_swaps", "ID Swap");
127 86862 : params.addParam<Real>("twist_pitch",
128 57908 : 0,
129 : "Pitch for helicoidal extrusion around an axis going through the origin "
130 : "following the direction vector");
131 28954 : return params;
132 0 : }
133 :
134 214 : AdvancedExtruderGenerator::AdvancedExtruderGenerator(const InputParameters & parameters)
135 : : MeshGenerator(parameters),
136 214 : _input(getMesh("input")),
137 214 : _heights(getParam<std::vector<Real>>("heights")),
138 631 : _biases(isParamValid("biases") ? getParam<std::vector<Real>>("biases")
139 417 : : std::vector<Real>(_heights.size(), 1.0)),
140 214 : _num_layers(getParam<std::vector<unsigned int>>("num_layers")),
141 214 : _subdomain_swaps(getParam<std::vector<std::vector<subdomain_id_type>>>("subdomain_swaps")),
142 214 : _boundary_swaps(getParam<std::vector<std::vector<boundary_id_type>>>("boundary_swaps")),
143 214 : _elem_integer_names_to_swap(getParam<std::vector<std::string>>("elem_integer_names_to_swap")),
144 214 : _elem_integers_swaps(
145 214 : getParam<std::vector<std::vector<std::vector<dof_id_type>>>>("elem_integers_swaps")),
146 214 : _direction(getParam<Point>("direction")),
147 214 : _has_top_boundary(isParamValid("top_boundary")),
148 214 : _top_boundary(isParamValid("top_boundary") ? getParam<BoundaryName>("top_boundary") : "0"),
149 214 : _has_bottom_boundary(isParamValid("bottom_boundary")),
150 214 : _bottom_boundary(isParamValid("bottom_boundary") ? getParam<BoundaryName>("bottom_boundary")
151 : : "0"),
152 627 : _upward_boundary_source_blocks(
153 428 : isParamValid("upward_boundary_source_blocks")
154 214 : ? getParam<std::vector<std::vector<subdomain_id_type>>>("upward_boundary_source_blocks")
155 199 : : std::vector<std::vector<subdomain_id_type>>(_heights.size(),
156 413 : std::vector<subdomain_id_type>())),
157 627 : _upward_boundary_ids(
158 428 : isParamValid("upward_boundary_ids")
159 214 : ? getParam<std::vector<std::vector<boundary_id_type>>>("upward_boundary_ids")
160 199 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
161 413 : std::vector<boundary_id_type>())),
162 841 : _downward_boundary_source_blocks(isParamValid("downward_boundary_source_blocks")
163 214 : ? getParam<std::vector<std::vector<subdomain_id_type>>>(
164 : "downward_boundary_source_blocks")
165 : : std::vector<std::vector<subdomain_id_type>>(
166 413 : _heights.size(), std::vector<subdomain_id_type>())),
167 627 : _downward_boundary_ids(
168 428 : isParamValid("downward_boundary_ids")
169 214 : ? getParam<std::vector<std::vector<boundary_id_type>>>("downward_boundary_ids")
170 199 : : std::vector<std::vector<boundary_id_type>>(_heights.size(),
171 413 : std::vector<boundary_id_type>())),
172 856 : _twist_pitch(getParam<Real>("twist_pitch"))
173 : {
174 214 : if (!_direction.norm())
175 0 : paramError("direction", "Must have some length!");
176 :
177 : // Normalize it
178 214 : _direction /= _direction.norm();
179 :
180 214 : const auto num_elevations = _heights.size();
181 :
182 214 : 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 214 : 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 214 : MooseMeshUtils::idSwapParametersProcessor(
195 214 : 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 214 : 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 214 : MooseMeshUtils::idSwapParametersProcessor(
212 214 : 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 227 : if (_elem_integers_swaps.size() &&
220 13 : _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 240 : for (const auto & unit_elem_integers_swaps : _elem_integers_swaps)
226 26 : 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 214 : MooseMeshUtils::extraElemIntegerSwapParametersProcessor(name(),
234 : num_elevations,
235 214 : _elem_integer_names_to_swap.size(),
236 : _elem_integers_swaps,
237 214 : _elem_integers_swap_pairs);
238 : }
239 0 : catch (const MooseException & e)
240 : {
241 0 : paramError("elem_integers_swaps", e.what());
242 0 : }
243 :
244 214 : bool has_negative_entry = false;
245 214 : bool has_positive_entry = false;
246 728 : for (const auto & h : _heights)
247 : {
248 514 : if (h > 0.0)
249 510 : has_positive_entry = true;
250 : else
251 4 : has_negative_entry = true;
252 : }
253 :
254 214 : if (has_negative_entry && has_positive_entry)
255 4 : paramError("heights", "Cannot have both positive and negative heights!");
256 210 : if (_biases.size() != _heights.size())
257 0 : paramError("biases", "Size of this parameter, if provided, must be the same as heights.");
258 :
259 420 : if (_upward_boundary_source_blocks.size() != _upward_boundary_ids.size() ||
260 210 : _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 712 : for (unsigned int i = 0; i < _upward_boundary_source_blocks.size(); i++)
268 502 : 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 420 : if (_downward_boundary_source_blocks.size() != _downward_boundary_ids.size() ||
274 210 : _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 712 : for (unsigned int i = 0; i < _downward_boundary_source_blocks.size(); i++)
282 502 : 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 210 : }
287 :
288 : std::unique_ptr<MeshBase>
289 200 : 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 200 : auto mesh = buildMeshBaseObject(_input->mesh_dimension() + 1);
296 200 : mesh->set_mesh_dimension(_input->mesh_dimension() + 1);
297 :
298 : // Check if the element integer names are existent in the input mesh.
299 224 : for (unsigned int i = 0; i < _elem_integer_names_to_swap.size(); i++)
300 24 : if (_input->has_elem_integer(_elem_integer_names_to_swap[i]))
301 48 : _elem_integer_indices_to_swap.push_back(
302 24 : _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 200 : const unsigned int num_extra_elem_integers = _input->n_elem_integers();
312 200 : std::vector<std::string> id_names;
313 :
314 224 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
315 : {
316 24 : id_names.push_back(_input->get_elem_integer_name(i));
317 24 : if (!mesh->has_elem_integer(id_names[i]))
318 24 : mesh->add_elem_integer(id_names[i]);
319 : }
320 :
321 : // retrieve subdomain/sideset/nodeset name maps
322 200 : const auto & input_subdomain_map = _input->get_subdomain_name_map();
323 200 : const auto & input_sideset_map = _input->get_boundary_info().get_sideset_name_map();
324 200 : 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 488 : for (const auto & swap : _subdomain_swaps)
328 1318 : for (const auto i : index_range(swap))
329 1030 : 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 220 : 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 630 : for (const auto & layer_vec : _upward_boundary_source_blocks)
340 508 : for (const auto bid : layer_vec)
341 70 : if (!MooseMeshUtils::hasSubdomainID(*_input, bid))
342 4 : paramError(
343 : "upward_boundary_source_blocks", "The block '", bid, "' was not found within the mesh");
344 614 : for (const auto & layer_vec : _downward_boundary_source_blocks)
345 496 : for (const auto bid : layer_vec)
346 70 : 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 184 : 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 184 : if (!input->is_serial())
357 25 : mesh->delete_remote_elements();
358 :
359 184 : unsigned int total_num_layers = std::accumulate(_num_layers.begin(), _num_layers.end(), 0);
360 :
361 184 : auto total_num_elevations = _heights.size();
362 :
363 184 : dof_id_type orig_elem = input->n_elem();
364 184 : dof_id_type orig_nodes = input->n_nodes();
365 :
366 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
367 184 : unique_id_type orig_unique_ids = input->parallel_max_unique_id();
368 : #endif
369 :
370 184 : unsigned int order = 1;
371 :
372 184 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
373 184 : const BoundaryInfo & input_boundary_info = input->get_boundary_info();
374 :
375 : // Determine boundary IDs for the new user provided boundary names
376 184 : std::vector<BoundaryName> new_boundary_names;
377 184 : if (_has_bottom_boundary)
378 102 : new_boundary_names.push_back(_bottom_boundary);
379 184 : if (_has_top_boundary)
380 102 : new_boundary_names.push_back(_top_boundary);
381 : std::vector<boundary_id_type> new_boundary_ids =
382 184 : MooseMeshUtils::getBoundaryIDs(*input, new_boundary_names, true);
383 : const auto user_bottom_boundary_id =
384 184 : _has_bottom_boundary ? new_boundary_ids.front() : libMesh::BoundaryInfo::invalid_id;
385 : const auto user_top_boundary_id =
386 184 : _has_top_boundary ? new_boundary_ids.back() : libMesh::BoundaryInfo::invalid_id;
387 :
388 : // We know a priori how many elements we'll need
389 184 : mesh->reserve_elem(total_num_layers * orig_elem);
390 :
391 : // Look for higher order elements which introduce an extra layer
392 184 : std::set<ElemType> higher_orders = {EDGE3, EDGE4, TRI6, TRI7, QUAD8, QUAD9};
393 184 : bool extruding_quad_eights = false;
394 184 : std::vector<ElemType> types;
395 184 : MeshTools::elem_types(*input, types);
396 396 : for (const auto elem_type : types)
397 : {
398 212 : if (higher_orders.count(elem_type))
399 9 : order = 2;
400 212 : if (elem_type == QUAD8)
401 9 : extruding_quad_eights = true;
402 : }
403 184 : mesh->comm().max(order);
404 184 : mesh->comm().max(extruding_quad_eights);
405 :
406 : // Reserve for the max number possibly needed
407 184 : 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 184 : std::vector<boundary_id_type> ids_to_copy;
411 :
412 184 : Point old_distance;
413 184 : Point current_distance;
414 :
415 : // Create translated layers of nodes in the direction of extrusion
416 12446 : for (const auto & node : input->node_ptr_range())
417 : {
418 12262 : unsigned int current_node_layer = 0;
419 :
420 12262 : old_distance.zero();
421 :
422 : // e is the elevation layer ordering
423 36602 : for (unsigned int e = 0; e < total_num_elevations; e++)
424 : {
425 24340 : auto num_layers = _num_layers[e];
426 :
427 24340 : auto height = _heights[e];
428 :
429 24340 : auto bias = _biases[e];
430 :
431 : // k is the element layer ordering within each elevation layer
432 78916 : 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 54576 : if (e == 0 && k == 0)
436 12262 : 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 42314 : auto layer_index = (k - (e == 0 ? 1 : 0)) / order + 1;
442 :
443 42314 : const auto step_size = MooseUtils::absoluteFuzzyEqual(bias, 1.0)
444 42314 : ? height / (Real)num_layers / (Real)order
445 2555 : : height * std::pow(bias, (Real)(layer_index - 1)) *
446 2555 : (1.0 - bias) /
447 42314 : (1.0 - std::pow(bias, (Real)(num_layers))) / (Real)order;
448 :
449 42314 : current_distance = old_distance + _direction * step_size;
450 :
451 : // Handle helicoidal extrusion
452 42314 : if (!MooseUtils::absoluteFuzzyEqual(_twist_pitch, 0.))
453 : {
454 : // twist 1 should be 'normal' to the extruded shape
455 3066 : RealVectorValue twist1 = _direction.cross(*node);
456 : // This happens for any node on the helicoidal extrusion axis
457 3066 : if (!MooseUtils::absoluteFuzzyEqual(twist1.norm(), .0))
458 3006 : twist1 /= twist1.norm();
459 3066 : const RealVectorValue twist2 = twist1.cross(_direction);
460 :
461 6132 : auto twist = (cos(2. * libMesh::pi * layer_index * step_size / _twist_pitch) -
462 6132 : cos(2. * libMesh::pi * (layer_index - 1) * step_size / _twist_pitch)) *
463 3066 : twist2 +
464 6132 : (sin(2. * libMesh::pi * layer_index * step_size / _twist_pitch) -
465 6132 : sin(2. * libMesh::pi * (layer_index - 1) * step_size / _twist_pitch)) *
466 6132 : twist1;
467 3066 : twist *= std::sqrt(node->norm_sq() + libMesh::Utility::pow<2>(_direction * (*node)));
468 3066 : current_distance += twist;
469 : }
470 : }
471 :
472 109152 : Node * new_node = mesh->add_point(*node + current_distance,
473 54576 : node->id() + (current_node_layer * orig_nodes),
474 54576 : 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 54576 : ? node->unique_id()
482 42314 : : orig_unique_ids +
483 42314 : (current_node_layer - 1) * (orig_nodes + orig_elem) +
484 42314 : node->id();
485 :
486 54576 : new_node->set_unique_id(uid);
487 : #endif
488 :
489 54576 : input_boundary_info.boundary_ids(node, ids_to_copy);
490 54576 : if (_boundary_swap_pairs.empty())
491 54576 : 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 54576 : old_distance = current_distance;
502 54576 : current_node_layer++;
503 : }
504 : }
505 184 : }
506 :
507 184 : const auto & side_ids = input_boundary_info.get_side_boundary_ids();
508 :
509 : boundary_id_type next_side_id =
510 184 : 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 184 : input->comm().max(next_side_id);
516 :
517 10335 : for (const auto & elem : input->element_ptr_range())
518 : {
519 10151 : const ElemType etype = elem->type();
520 :
521 : // build_extrusion currently only works on coarse meshes
522 : libmesh_assert(!elem->parent());
523 :
524 10151 : unsigned int current_layer = 0;
525 :
526 29886 : for (unsigned int e = 0; e != total_num_elevations; e++)
527 : {
528 19735 : auto num_layers = _num_layers[e];
529 :
530 53486 : for (unsigned int k = 0; k != num_layers; ++k)
531 : {
532 33751 : std::unique_ptr<Elem> new_elem;
533 33751 : bool is_flipped(false);
534 33751 : switch (etype)
535 : {
536 11 : case EDGE2:
537 : {
538 11 : new_elem = std::make_unique<Quad4>();
539 22 : new_elem->set_node(
540 11 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
541 22 : new_elem->set_node(
542 11 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
543 22 : new_elem->set_node(
544 11 : 2, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
545 22 : new_elem->set_node(
546 11 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
547 :
548 11 : if (elem->neighbor_ptr(0) == remote_elem)
549 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
550 11 : if (elem->neighbor_ptr(1) == remote_elem)
551 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
552 :
553 11 : 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 1044 : case TRI3:
591 : {
592 1044 : new_elem = std::make_unique<Prism6>();
593 2088 : new_elem->set_node(
594 1044 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
595 2088 : new_elem->set_node(
596 1044 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
597 2088 : new_elem->set_node(
598 1044 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
599 2088 : new_elem->set_node(
600 1044 : 3, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
601 2088 : new_elem->set_node(
602 1044 : 4, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
603 2088 : new_elem->set_node(
604 1044 : 5, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
605 :
606 1044 : if (elem->neighbor_ptr(0) == remote_elem)
607 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
608 1044 : if (elem->neighbor_ptr(1) == remote_elem)
609 12 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
610 1044 : if (elem->neighbor_ptr(2) == remote_elem)
611 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
612 :
613 1044 : 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 1044 : 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 32678 : case QUAD4:
777 : {
778 32678 : new_elem = std::make_unique<Hex8>();
779 65356 : new_elem->set_node(
780 32678 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * orig_nodes)));
781 65356 : new_elem->set_node(
782 32678 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * orig_nodes)));
783 65356 : new_elem->set_node(
784 32678 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * orig_nodes)));
785 65356 : new_elem->set_node(
786 32678 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (current_layer * orig_nodes)));
787 65356 : new_elem->set_node(
788 32678 : 4, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer + 1) * orig_nodes)));
789 65356 : new_elem->set_node(
790 32678 : 5, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer + 1) * orig_nodes)));
791 65356 : new_elem->set_node(
792 32678 : 6, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer + 1) * orig_nodes)));
793 65356 : new_elem->set_node(
794 32678 : 7, mesh->node_ptr(elem->node_ptr(3)->id() + ((current_layer + 1) * orig_nodes)));
795 :
796 32678 : if (elem->neighbor_ptr(0) == remote_elem)
797 94 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
798 32678 : if (elem->neighbor_ptr(1) == remote_elem)
799 318 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
800 32678 : if (elem->neighbor_ptr(2) == remote_elem)
801 110 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
802 32678 : if (elem->neighbor_ptr(3) == remote_elem)
803 306 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
804 :
805 32678 : if (new_elem->volume() < 0.0)
806 : {
807 2241 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
808 2241 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
809 2241 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
810 2241 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
811 2241 : is_flipped = true;
812 : }
813 :
814 32678 : break;
815 : }
816 18 : case QUAD8:
817 : {
818 18 : new_elem = std::make_unique<Hex20>();
819 36 : new_elem->set_node(
820 18 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (2 * current_layer * orig_nodes)));
821 36 : new_elem->set_node(
822 18 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (2 * current_layer * orig_nodes)));
823 36 : new_elem->set_node(
824 18 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (2 * current_layer * orig_nodes)));
825 36 : new_elem->set_node(
826 18 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (2 * current_layer * orig_nodes)));
827 36 : new_elem->set_node(
828 : 4,
829 18 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 2) * orig_nodes)));
830 36 : new_elem->set_node(
831 : 5,
832 18 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 2) * orig_nodes)));
833 36 : new_elem->set_node(
834 : 6,
835 18 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 2) * orig_nodes)));
836 36 : new_elem->set_node(
837 : 7,
838 18 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 2) * orig_nodes)));
839 36 : new_elem->set_node(
840 18 : 8, mesh->node_ptr(elem->node_ptr(4)->id() + (2 * current_layer * orig_nodes)));
841 36 : new_elem->set_node(
842 18 : 9, mesh->node_ptr(elem->node_ptr(5)->id() + (2 * current_layer * orig_nodes)));
843 36 : new_elem->set_node(
844 18 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (2 * current_layer * orig_nodes)));
845 36 : new_elem->set_node(
846 18 : 11, mesh->node_ptr(elem->node_ptr(7)->id() + (2 * current_layer * orig_nodes)));
847 36 : new_elem->set_node(
848 : 12,
849 18 : mesh->node_ptr(elem->node_ptr(0)->id() + ((2 * current_layer + 1) * orig_nodes)));
850 36 : new_elem->set_node(
851 : 13,
852 18 : mesh->node_ptr(elem->node_ptr(1)->id() + ((2 * current_layer + 1) * orig_nodes)));
853 36 : new_elem->set_node(
854 : 14,
855 18 : mesh->node_ptr(elem->node_ptr(2)->id() + ((2 * current_layer + 1) * orig_nodes)));
856 36 : new_elem->set_node(
857 : 15,
858 18 : mesh->node_ptr(elem->node_ptr(3)->id() + ((2 * current_layer + 1) * orig_nodes)));
859 36 : new_elem->set_node(
860 : 16,
861 18 : mesh->node_ptr(elem->node_ptr(4)->id() + ((2 * current_layer + 2) * orig_nodes)));
862 36 : new_elem->set_node(
863 : 17,
864 18 : mesh->node_ptr(elem->node_ptr(5)->id() + ((2 * current_layer + 2) * orig_nodes)));
865 36 : new_elem->set_node(
866 : 18,
867 18 : mesh->node_ptr(elem->node_ptr(6)->id() + ((2 * current_layer + 2) * orig_nodes)));
868 36 : new_elem->set_node(
869 : 19,
870 18 : mesh->node_ptr(elem->node_ptr(7)->id() + ((2 * current_layer + 2) * orig_nodes)));
871 :
872 18 : if (elem->neighbor_ptr(0) == remote_elem)
873 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
874 18 : if (elem->neighbor_ptr(1) == remote_elem)
875 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
876 18 : if (elem->neighbor_ptr(2) == remote_elem)
877 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
878 18 : if (elem->neighbor_ptr(3) == remote_elem)
879 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
880 :
881 18 : if (new_elem->volume() < 0.0)
882 : {
883 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
884 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
885 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
886 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
887 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
888 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
889 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
890 9 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
891 9 : is_flipped = true;
892 : }
893 :
894 18 : 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 33751 : new_elem->set_id(elem->id() + (current_layer * orig_elem));
1002 33751 : 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 33751 : ? elem->unique_id()
1010 23600 : : orig_unique_ids +
1011 23600 : (current_layer - 1) * (orig_nodes + orig_elem) +
1012 23600 : orig_nodes + elem->id();
1013 :
1014 33751 : new_elem->set_unique_id(uid);
1015 : #endif
1016 :
1017 : // maintain the subdomain_id
1018 33751 : new_elem->subdomain_id() = elem->subdomain_id();
1019 :
1020 : // define upward boundaries
1021 33751 : 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 19735 : 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 25639 : for (unsigned int i = 0; i < _upward_boundary_source_blocks[e].size(); i++)
1029 5904 : if (new_elem->subdomain_id() == _upward_boundary_source_blocks[e][i])
1030 3936 : boundary_info.add_side(
1031 3936 : new_elem.get(), is_flipped ? 0 : top_id, _upward_boundary_ids[e][i]);
1032 : }
1033 : // define downward boundaries
1034 33751 : if (k == 0)
1035 : {
1036 : const unsigned short top_id =
1037 19735 : new_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1038 25639 : for (unsigned int i = 0; i < _downward_boundary_source_blocks[e].size(); i++)
1039 5904 : if (new_elem->subdomain_id() == _downward_boundary_source_blocks[e][i])
1040 3936 : boundary_info.add_side(
1041 3936 : new_elem.get(), is_flipped ? top_id : 0, _downward_boundary_ids[e][i]);
1042 : }
1043 :
1044 : // perform subdomain swaps
1045 33751 : if (_subdomain_swap_pairs.size())
1046 : {
1047 14940 : auto & elevation_swap_pairs = _subdomain_swap_pairs[e];
1048 :
1049 14940 : auto new_id_it = elevation_swap_pairs.find(elem->subdomain_id());
1050 :
1051 14940 : if (new_id_it != elevation_swap_pairs.end())
1052 14940 : new_elem->subdomain_id() = new_id_it->second;
1053 : }
1054 :
1055 33751 : Elem * added_elem = mesh->add_elem(std::move(new_elem));
1056 :
1057 : // maintain extra integers
1058 47575 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
1059 13824 : added_elem->set_extra_integer(i, elem->get_extra_integer(i));
1060 :
1061 33751 : if (_elem_integers_swap_pairs.size())
1062 : {
1063 20736 : for (unsigned int i = 0; i < _elem_integer_indices_to_swap.size(); i++)
1064 : {
1065 13824 : auto & elevation_extra_swap_pairs = _elem_integers_swap_pairs[i * _heights.size() + e];
1066 :
1067 13824 : auto new_extra_id_it = elevation_extra_swap_pairs.find(
1068 13824 : elem->get_extra_integer(_elem_integer_indices_to_swap[i]));
1069 :
1070 13824 : if (new_extra_id_it != elevation_extra_swap_pairs.end())
1071 8064 : added_elem->set_extra_integer(_elem_integer_indices_to_swap[i],
1072 8064 : new_extra_id_it->second);
1073 : }
1074 : }
1075 :
1076 : // Copy any old boundary ids on all sides
1077 167689 : for (auto s : elem->side_index_range())
1078 : {
1079 133938 : input_boundary_info.boundary_ids(elem, s, ids_to_copy);
1080 :
1081 133938 : 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 133916 : if (_boundary_swap_pairs.empty())
1088 133916 : 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 22 : const unsigned short sidemap[2] = {3, 1};
1105 22 : if (_boundary_swap_pairs.empty())
1106 22 : 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 33751 : if (current_layer == 0)
1119 : {
1120 : const unsigned short top_id =
1121 10151 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1122 10151 : 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 2854 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, user_bottom_boundary_id);
1127 : }
1128 : else
1129 7297 : boundary_info.add_side(added_elem, is_flipped ? top_id : 0, next_side_id);
1130 : }
1131 :
1132 33751 : 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 10151 : added_elem->dim() == 3 ? cast_int<unsigned short>(elem->n_sides() + 1) : 2;
1139 :
1140 10151 : 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 2854 : boundary_info.add_side(added_elem, is_flipped ? 0 : top_id, user_top_boundary_id);
1145 : }
1146 : else
1147 7297 : boundary_info.add_side(
1148 7297 : added_elem, is_flipped ? 0 : top_id, cast_int<boundary_id_type>(next_side_id + 1));
1149 : }
1150 :
1151 33751 : current_layer++;
1152 33751 : }
1153 : }
1154 184 : }
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 184 : unsigned int total_new_node_layers = total_num_layers * order;
1160 184 : unsigned int new_unique_ids = orig_unique_ids + (total_new_node_layers - 1) * orig_elem +
1161 : total_new_node_layers * orig_nodes;
1162 184 : 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 184 : if (!input_subdomain_map.empty())
1167 29 : mesh->set_subdomain_name_map().insert(input_subdomain_map.begin(), input_subdomain_map.end());
1168 184 : if (!input_sideset_map.empty())
1169 184 : mesh->get_boundary_info().set_sideset_name_map().insert(input_sideset_map.begin(),
1170 : input_sideset_map.end());
1171 184 : if (!input_nodeset_map.empty())
1172 121 : mesh->get_boundary_info().set_nodeset_name_map().insert(input_nodeset_map.begin(),
1173 : input_nodeset_map.end());
1174 :
1175 184 : if (_has_bottom_boundary)
1176 102 : boundary_info.sideset_name(new_boundary_ids.front()) = new_boundary_names.front();
1177 184 : if (_has_top_boundary)
1178 102 : boundary_info.sideset_name(new_boundary_ids.back()) = new_boundary_names.back();
1179 :
1180 184 : 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 184 : if (extruding_quad_eights)
1184 9 : mesh->prepare_for_use();
1185 :
1186 368 : return mesh;
1187 184 : }
|
