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 "RevolveGenerator.h"
11 : #include "PolygonalMeshGenerationUtils.h"
12 :
13 : #include "libmesh/cell_prism6.h"
14 : #include "libmesh/cell_prism15.h"
15 : #include "libmesh/cell_prism18.h"
16 : #include "libmesh/cell_prism21.h"
17 : #include "libmesh/cell_pyramid5.h"
18 : #include "libmesh/cell_pyramid13.h"
19 : #include "libmesh/cell_pyramid14.h"
20 : #include "libmesh/cell_pyramid18.h"
21 : #include "libmesh/cell_tet4.h"
22 : #include "libmesh/cell_tet10.h"
23 : #include "libmesh/cell_tet14.h"
24 : #include "libmesh/cell_hex8.h"
25 : #include "libmesh/cell_hex20.h"
26 : #include "libmesh/cell_hex27.h"
27 : #include "libmesh/face_tri3.h"
28 : #include "libmesh/face_tri7.h"
29 : #include "libmesh/face_quad4.h"
30 : #include "libmesh/face_quad9.h"
31 : #include "libmesh/point.h"
32 : #include "libmesh/mesh_tools.h"
33 :
34 : using namespace libMesh;
35 :
36 : // C++ includes
37 : #include <cmath>
38 :
39 : registerMooseObject("ReactorApp", RevolveGenerator);
40 :
41 : InputParameters
42 430 : RevolveGenerator::validParams()
43 : {
44 430 : InputParameters params = PolygonMeshGeneratorBase::validParams();
45 430 : params.addClassDescription("This RevolveGenerator object is designed to revolve a 1D mesh into "
46 : "2D, or a 2D mesh into 3D based on an axis.");
47 :
48 860 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh to revolve");
49 :
50 860 : params.addRequiredParam<Point>("axis_point", "A point on the axis of revolution");
51 :
52 860 : params.addRequiredParam<Point>("axis_direction", "The direction of the axis of revolution");
53 :
54 860 : params.addRangeCheckedParam<std::vector<Real>>(
55 : "revolving_angles",
56 : "revolving_angles<=360.0 & revolving_angles>0.0",
57 : "The angles delineating each azimuthal section of revolution around the axis in degrees");
58 :
59 860 : params.addParam<std::vector<std::vector<subdomain_id_type>>>(
60 : "subdomain_swaps",
61 : {},
62 : "For each row, every two entries are interpreted as a pair of "
63 : "'from' and 'to' to remap the subdomains for that azimuthal section");
64 :
65 860 : params.addParam<std::vector<std::vector<boundary_id_type>>>(
66 : "boundary_swaps",
67 : {},
68 : "For each row, every two entries are interpreted as a pair of "
69 : "'from' and 'to' to remap the boundaries for that elevation");
70 :
71 860 : params.addParam<std::vector<std::string>>(
72 : "elem_integer_names_to_swap",
73 : {},
74 : "Array of element extra integer names that need to be swapped during revolving.");
75 :
76 860 : params.addParam<std::vector<std::vector<std::vector<dof_id_type>>>>(
77 : "elem_integers_swaps",
78 : {},
79 : "For each row, every two entries are interpreted as a pair of 'from' and 'to' to remap the "
80 : "element extra integer for that elevation. If multiple element extra integers need to be "
81 : "swapped, the enties are stacked based on the order provided in "
82 : "'elem_integer_names_to_swap' to form the third dimension.");
83 :
84 860 : params.addParam<boundary_id_type>(
85 : "start_boundary",
86 : "The boundary ID to set on the starting boundary for a partial revolution.");
87 :
88 860 : params.addParam<boundary_id_type>(
89 : "end_boundary", "The boundary ID to set on the ending boundary for partial revolving.");
90 :
91 860 : params.addParam<bool>(
92 860 : "clockwise", true, "Revolve clockwise around the axis or not (i.e., counterclockwise)");
93 :
94 860 : params.addRequiredParam<std::vector<unsigned int>>(
95 : "nums_azimuthal_intervals",
96 : "List of the numbers of azimuthal interval discretization for each azimuthal section");
97 :
98 860 : params.addParam<bool>("preserve_volumes",
99 860 : false,
100 : "Whether the volume of the revolved mesh is preserving the circular area "
101 : "by modifying (expanding) the radius to account for polygonization.");
102 :
103 860 : params.addParamNamesToGroup("start_boundary end_boundary", "Boundary Assignment");
104 860 : params.addParamNamesToGroup(
105 : "subdomain_swaps boundary_swaps elem_integer_names_to_swap elem_integers_swaps", "ID Swap");
106 :
107 430 : return params;
108 0 : }
109 :
110 215 : RevolveGenerator::RevolveGenerator(const InputParameters & parameters)
111 : : PolygonMeshGeneratorBase(parameters),
112 215 : _input(getMesh("input")),
113 430 : _axis_point(getParam<Point>("axis_point")),
114 430 : _axis_direction(getParam<Point>("axis_direction")),
115 524 : _revolving_angles(isParamValid("revolving_angles")
116 215 : ? getParam<std::vector<Real>>("revolving_angles")
117 215 : : std::vector<Real>(1, 360.0)),
118 430 : _subdomain_swaps(getParam<std::vector<std::vector<subdomain_id_type>>>("subdomain_swaps")),
119 430 : _boundary_swaps(getParam<std::vector<std::vector<boundary_id_type>>>("boundary_swaps")),
120 430 : _elem_integer_names_to_swap(getParam<std::vector<std::string>>("elem_integer_names_to_swap")),
121 215 : _elem_integers_swaps(
122 215 : getParam<std::vector<std::vector<std::vector<dof_id_type>>>>("elem_integers_swaps")),
123 430 : _clockwise(getParam<bool>("clockwise")),
124 430 : _nums_azimuthal_intervals(getParam<std::vector<unsigned int>>("nums_azimuthal_intervals")),
125 430 : _preserve_volumes(getParam<bool>("preserve_volumes")),
126 430 : _has_start_boundary(isParamValid("start_boundary")),
127 448 : _start_boundary(isParamValid("start_boundary") ? getParam<boundary_id_type>("start_boundary")
128 : : 0),
129 430 : _has_end_boundary(isParamValid("end_boundary")),
130 448 : _end_boundary(isParamValid("end_boundary") ? getParam<boundary_id_type>("end_boundary") : 0),
131 430 : _radius_correction_factor(1.0)
132 : {
133 215 : if (_revolving_angles.size() != _nums_azimuthal_intervals.size())
134 0 : paramError("nums_azimuthal_intervals",
135 : "The number of azimuthal intervals should be the same as the number of revolving "
136 : "angles.");
137 215 : if (_subdomain_swaps.size() && (_subdomain_swaps.size() != _nums_azimuthal_intervals.size()))
138 0 : paramError(
139 : "subdomain_swaps",
140 : "If specified, 'subdomain_swaps' must be the same length as 'nums_azimuthal_intervals'.");
141 :
142 215 : if (_boundary_swaps.size() && (_boundary_swaps.size() != _nums_azimuthal_intervals.size()))
143 0 : paramError(
144 : "boundary_swaps",
145 : "If specified, 'boundary_swaps' must be the same length as 'nums_azimuthal_intervals'.");
146 :
147 233 : for (const auto & unit_elem_integers_swaps : _elem_integers_swaps)
148 18 : if (unit_elem_integers_swaps.size() != _nums_azimuthal_intervals.size())
149 0 : paramError("elem_integers_swaps",
150 : "If specified, each element of 'elem_integers_swaps' must have the same length as "
151 : "the length of 'nums_azimuthal_intervals'.");
152 :
153 215 : if (_elem_integers_swaps.size() &&
154 9 : _elem_integers_swaps.size() != _elem_integer_names_to_swap.size())
155 0 : paramError("elem_integers_swaps",
156 : "If specified, 'elem_integers_swaps' must have the same length as the length of "
157 : "'elem_integer_names_to_swap'.");
158 :
159 215 : _full_circle_revolving =
160 : MooseUtils::absoluteFuzzyEqual(
161 : std::accumulate(_revolving_angles.begin(), _revolving_angles.end(), 0), 360.0)
162 : ? true
163 : : false;
164 215 : if (MooseUtils::absoluteFuzzyGreaterThan(
165 : std::accumulate(_revolving_angles.begin(), _revolving_angles.end(), 0), 360.0))
166 0 : paramError("revolving_angles",
167 : "The sum of revolving angles should be less than or equal to 360.");
168 :
169 215 : if ((_has_start_boundary && _full_circle_revolving) ||
170 206 : (_has_end_boundary && _full_circle_revolving))
171 0 : paramError("full_circle_revolving",
172 : "starting or ending boundaries can only be assigned for partial revolving.");
173 :
174 : try
175 : {
176 215 : MooseMeshUtils::idSwapParametersProcessor(
177 215 : name(), "subdomain_swaps", _subdomain_swaps, _subdomain_swap_pairs);
178 : }
179 0 : catch (const MooseException & e)
180 : {
181 0 : paramError("subdomain_swaps", e.what());
182 0 : }
183 :
184 : try
185 : {
186 215 : MooseMeshUtils::idSwapParametersProcessor(
187 215 : name(), "boundary_swaps", _boundary_swaps, _boundary_swap_pairs);
188 : }
189 0 : catch (const MooseException & e)
190 : {
191 0 : paramError("boundary_swaps", e.what());
192 0 : }
193 :
194 : try
195 : {
196 215 : MooseMeshUtils::extraElemIntegerSwapParametersProcessor(name(),
197 215 : _nums_azimuthal_intervals.size(),
198 215 : _elem_integer_names_to_swap.size(),
199 : _elem_integers_swaps,
200 215 : _elem_integers_swap_pairs);
201 : }
202 0 : catch (const MooseException & e)
203 : {
204 0 : paramError("elem_integers_swaps", e.what());
205 0 : }
206 215 : }
207 :
208 : std::unique_ptr<MeshBase>
209 215 : RevolveGenerator::generate()
210 : {
211 : // Note: Inspired by AdvancedExtruderGenerator::generate()
212 :
213 215 : auto mesh = buildMeshBaseObject();
214 :
215 : // Only works for 1D and 2D input meshes
216 215 : if (_input->mesh_dimension() > 2)
217 2 : paramError("input", "This mesh generator only works for 1D and 2D input meshes.");
218 :
219 213 : mesh->set_mesh_dimension(_input->mesh_dimension() + 1);
220 :
221 : // Check if the element integer names are existent in the input mesh.
222 231 : for (const auto i : index_range(_elem_integer_names_to_swap))
223 18 : if (_input->has_elem_integer(_elem_integer_names_to_swap[i]))
224 18 : _elem_integer_indices_to_swap.push_back(
225 18 : _input->get_elem_integer_index(_elem_integer_names_to_swap[i]));
226 : else
227 0 : paramError("elem_integer_names_to_swap",
228 : "Element ",
229 : i + 1,
230 : " of 'elem_integer_names_to_swap' is not a valid extra element integer of the "
231 : "input mesh.");
232 :
233 : // prepare for transferring extra element integers from original mesh to the revolved mesh.
234 213 : const unsigned int num_extra_elem_integers = _input->n_elem_integers();
235 : std::vector<std::string> id_names;
236 :
237 231 : for (const auto i : make_range(num_extra_elem_integers))
238 : {
239 18 : id_names.push_back(_input->get_elem_integer_name(i));
240 18 : if (!mesh->has_elem_integer(id_names[i]))
241 36 : mesh->add_elem_integer(id_names[i]);
242 : }
243 :
244 : // retrieve subdomain/sideset/nodeset name maps
245 213 : const auto & input_subdomain_map = _input->get_subdomain_name_map();
246 : const auto & input_sideset_map = _input->get_boundary_info().get_sideset_name_map();
247 : const auto & input_nodeset_map = _input->get_boundary_info().get_nodeset_name_map();
248 :
249 : std::unique_ptr<MeshBase> input = std::move(_input);
250 :
251 : // If we're using a distributed mesh... then make sure we don't have any remote elements hanging
252 : // around
253 213 : if (!input->is_serial())
254 36 : mesh->delete_remote_elements();
255 :
256 : // check that subdomain swap sources exist in the mesh
257 : std::set<subdomain_id_type> blocks;
258 213 : input->subdomain_ids(blocks, true);
259 251 : for (const auto & swap_map : _subdomain_swap_pairs)
260 76 : for (const auto & [bid, tbid] : swap_map)
261 : {
262 : libmesh_ignore(tbid);
263 : if (blocks.count(bid) == 0)
264 2 : paramError("subdomain_swaps",
265 2 : "Source subdomain " + std::to_string(bid) + " was not found in the mesh");
266 : }
267 :
268 : // Subdomain IDs for on-axis elements must be new
269 : std::set<subdomain_id_type> subdomain_ids_set;
270 211 : input->subdomain_ids(subdomain_ids_set);
271 211 : const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin();
272 : const subdomain_id_type tri_to_pyramid_subdomain_id_shift =
273 211 : std::max((int)max_subdomain_id, 1) + 1;
274 : const subdomain_id_type tri_to_tet_subdomain_id_shift =
275 211 : std::max((int)max_subdomain_id, 1) * 2 + 1;
276 211 : const subdomain_id_type quad_to_prism_subdomain_id_shift = std::max((int)max_subdomain_id, 1) + 1;
277 : const subdomain_id_type quad_to_pyramid_subdomain_id_shift =
278 211 : std::max((int)max_subdomain_id, 1) * 2 + 1;
279 : const subdomain_id_type quad_to_hi_pyramid_subdomain_id_shift =
280 211 : std::max((int)max_subdomain_id, 1) * 3 + 1;
281 211 : const subdomain_id_type edge_to_tri_subdomain_id_shift = std::max((int)max_subdomain_id, 1) + 1;
282 :
283 : // Get the centroid of the input mesh
284 211 : const auto input_centroid = MooseMeshUtils::meshCentroidCalculator(*input);
285 211 : const Point axis_centroid_cross = (input_centroid - _axis_point).cross(_axis_direction);
286 :
287 211 : if (MooseUtils::absoluteFuzzyEqual(axis_centroid_cross.norm(), 0.0))
288 2 : mooseError("The input mesh is either across the axis or overlapped with the axis!");
289 :
290 : Real inner_product_1d(0.0);
291 : bool inner_product_1d_initialized(false);
292 : // record ids of nodes on the axis
293 : std::vector<dof_id_type> node_ids_on_axis;
294 7086 : for (const auto & node : input->node_ptr_range())
295 : {
296 3342 : const Point axis_node_cross = (*node - _axis_point).cross(_axis_direction);
297 : // if the cross product is zero, then the node is on the axis
298 3342 : if (!MooseUtils::absoluteFuzzyEqual(axis_node_cross.norm(), 0.0))
299 : {
300 3066 : if (MooseUtils::absoluteFuzzyLessThan(axis_node_cross * axis_centroid_cross, 0.0))
301 4 : mooseError("The input mesh is across the axis.");
302 3062 : else if (MooseUtils::absoluteFuzzyLessThan(axis_node_cross * axis_centroid_cross,
303 6124 : axis_centroid_cross.norm() *
304 3062 : axis_node_cross.norm()))
305 2 : mooseError("The input mesh is not in the same plane with the rotation axis.");
306 : }
307 : else
308 276 : node_ids_on_axis.push_back(node->id());
309 :
310 : // Only for 1D input mesh, we need to check if the axis is perpendicular to the input mesh
311 3336 : if (input->mesh_dimension() == 1)
312 : {
313 106 : const Real temp_inner_product = (*node - _axis_point) * _axis_direction.unit();
314 106 : if (inner_product_1d_initialized)
315 : {
316 74 : if (!MooseUtils::absoluteFuzzyEqual(temp_inner_product, inner_product_1d))
317 2 : mooseError("The 1D input mesh is not perpendicular to the rotation axis.");
318 : }
319 : else
320 : {
321 : inner_product_1d_initialized = true;
322 : inner_product_1d = temp_inner_product;
323 : }
324 : }
325 201 : }
326 :
327 : // If there are any on-axis nodes, we need to check if there are any QUAD8 elements with one
328 : // vertex on the axis. If so, we need to replace it with a QUAD9 element.
329 201 : if (!node_ids_on_axis.empty())
330 : {
331 : // Sort the vector for using set_intersection
332 120 : std::sort(node_ids_on_axis.begin(), node_ids_on_axis.end());
333 : // For QUAD8 elements with one vertex on the axis, we need to replace it with a QUAD9 element
334 : std::set<subdomain_id_type> converted_quad8_subdomain_ids;
335 1184 : for (const auto & elem : input->element_ptr_range())
336 : {
337 472 : if (elem->type() == QUAD8)
338 : {
339 : std::vector<dof_id_type> elem_vertex_node_ids;
340 360 : for (unsigned int i = 0; i < 4; i++)
341 : {
342 288 : elem_vertex_node_ids.push_back(elem->node_id(i));
343 : }
344 72 : std::sort(elem_vertex_node_ids.begin(), elem_vertex_node_ids.end());
345 : std::vector<dof_id_type> common_node_ids;
346 72 : std::set_intersection(node_ids_on_axis.begin(),
347 : node_ids_on_axis.end(),
348 : elem_vertex_node_ids.begin(),
349 : elem_vertex_node_ids.end(),
350 : std::back_inserter(common_node_ids));
351 : // Temporarily shift the subdomain ID to mark the element
352 72 : if (common_node_ids.size() == 1)
353 : {
354 : // we borrow quad_to_hi_pyramid_subdomain_id_shift here
355 9 : elem->subdomain_id() += quad_to_hi_pyramid_subdomain_id_shift;
356 9 : converted_quad8_subdomain_ids.emplace(elem->subdomain_id());
357 : }
358 72 : }
359 120 : }
360 : // Convert the recorded subdomains
361 120 : input->all_second_order_range(
362 240 : input->active_subdomain_set_elements_ptr_range(converted_quad8_subdomain_ids));
363 : // Restore the subdomain ID; we do not worry about repeated subdomain IDs because those QUAD9
364 : // will become PYRAMID and PRISM elements with new shifts
365 369 : for (auto elem : input->active_subdomain_set_elements_ptr_range(converted_quad8_subdomain_ids))
366 129 : elem->subdomain_id() -= quad_to_hi_pyramid_subdomain_id_shift;
367 : }
368 :
369 : // We should only record this info after QUAD8->QUAD9 conversion
370 201 : dof_id_type orig_elem = input->n_elem();
371 201 : dof_id_type orig_nodes = input->n_nodes();
372 :
373 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
374 : // Add the number of original elements as revolving may create two elements per layer for one
375 : // original element
376 201 : unique_id_type orig_unique_ids = input->parallel_max_unique_id() + orig_elem;
377 : #endif
378 :
379 : // get rotation vectors
380 201 : const auto rotation_vectors = rotationVectors(_axis_point, _axis_direction, input_centroid);
381 :
382 201 : unsigned int order = 1;
383 :
384 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
385 : const BoundaryInfo & input_boundary_info = input->get_boundary_info();
386 :
387 : const unsigned int total_num_azimuthal_intervals =
388 201 : std::accumulate(_nums_azimuthal_intervals.begin(), _nums_azimuthal_intervals.end(), 0);
389 : // We know a priori how many elements we'll need
390 : // In the worst case, all quad elements will become two elements per layer
391 201 : mesh->reserve_elem(total_num_azimuthal_intervals * orig_elem * 2);
392 : const dof_id_type elem_id_shift = total_num_azimuthal_intervals * orig_elem;
393 :
394 : // Look for higher order elements which introduce an extra layer
395 201 : std::set<ElemType> higher_orders = {EDGE3, TRI6, TRI7, QUAD8, QUAD9};
396 : std::vector<ElemType> types;
397 201 : MeshTools::elem_types(*input, types);
398 438 : for (const auto elem_type : types)
399 : if (higher_orders.count(elem_type))
400 78 : order = 2;
401 201 : mesh->comm().max(order);
402 :
403 : // Collect azimuthal angles and use them to calculate the correction factor if applicable
404 : std::vector<Real> azi_array;
405 429 : for (const auto & i : index_range(_revolving_angles))
406 : {
407 : const Real section_start_angle =
408 228 : azi_array.empty() ? 0.0 : (azi_array.back() + _unit_angles.back());
409 228 : _unit_angles.push_back(_revolving_angles[i] / _nums_azimuthal_intervals[i] / order);
410 1956 : for (unsigned int j = 0; j < _nums_azimuthal_intervals[i] * order; j++)
411 1728 : azi_array.push_back(section_start_angle + _unit_angles.back() * (Real)j);
412 : }
413 201 : if (_preserve_volumes)
414 : {
415 36 : _radius_correction_factor = PolygonalMeshGenerationUtils::radiusCorrectionFactor(
416 18 : azi_array, _full_circle_revolving, order);
417 :
418 : // In the meanwhile, modify the input mesh for radius correction if applicable
419 288 : for (const auto & node : input->node_ptr_range())
420 144 : nodeModification(*node);
421 : }
422 :
423 201 : mesh->reserve_nodes(
424 201 : (order * total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
425 : orig_nodes);
426 :
427 : // Container to catch the boundary IDs handed back by the BoundaryInfo object
428 : std::vector<boundary_id_type> ids_to_copy;
429 :
430 : Point old_distance;
431 : Point current_distance;
432 201 : if (!_clockwise)
433 : std::transform(_unit_angles.begin(),
434 : _unit_angles.end(),
435 : _unit_angles.begin(),
436 9 : [](auto & c) { return c * (-1.0) * M_PI / 180.0; });
437 : else
438 : std::transform(_unit_angles.begin(),
439 : _unit_angles.end(),
440 : _unit_angles.begin(),
441 219 : [](auto & c) { return c * M_PI / 180.0; });
442 : std::vector<dof_id_type> nodes_on_axis;
443 :
444 7084 : for (const auto & node : input->node_ptr_range())
445 : {
446 : // Calculate the radius and corresponding center point on the rotation axis
447 : // If the radius is 0, then the node is on the axis
448 3341 : const auto radius_and_center = getRotationCenterAndRadius(*node, _axis_point, _axis_direction);
449 : const bool isOnAxis = MooseUtils::absoluteFuzzyEqual(radius_and_center.first, 0.0);
450 3341 : if (isOnAxis)
451 : {
452 274 : nodes_on_axis.push_back(node->id());
453 : }
454 :
455 : unsigned int current_node_layer = 0;
456 :
457 : old_distance.zero();
458 :
459 3341 : const unsigned int num_rotations = _revolving_angles.size();
460 7933 : for (unsigned int e = 0; e < num_rotations; e++)
461 : {
462 4592 : auto num_layers = _nums_azimuthal_intervals[e];
463 :
464 4592 : auto angle = _unit_angles[e];
465 :
466 : const auto base_angle =
467 4592 : std::accumulate(_revolving_angles.begin(), _revolving_angles.begin() + e, 0.0) / 180.0 *
468 4592 : M_PI;
469 :
470 4592 : for (unsigned int k = 0;
471 71332 : k < order * num_layers + (e == 0 ? 1 : 0) -
472 35666 : (e == num_rotations - 1 ? (unsigned int)_full_circle_revolving : 0);
473 : ++k)
474 : {
475 : bool is_node_created(false);
476 31074 : if (!isOnAxis)
477 : {
478 : // For the first layer we don't need to move
479 28200 : if (e == 0 && k == 0)
480 : current_distance.zero();
481 : else
482 : {
483 25133 : auto layer_index = (k - (e == 0 ? 1 : 0)) + 1;
484 :
485 : // Calculate the rotation angle in XY Plane
486 : const Point vector_xy =
487 : Point(-2.0 * radius_and_center.first *
488 : std::sin((base_angle + angle * (Real)layer_index) / 2.0) *
489 : std::sin((base_angle + angle * (Real)layer_index) / 2.0),
490 25133 : 2.0 * radius_and_center.first *
491 25133 : std::sin((base_angle + angle * (Real)layer_index) / 2.0) *
492 25133 : std::cos((base_angle + angle * (Real)layer_index) / 2.0),
493 25133 : 0.0);
494 25133 : current_distance = Point(rotation_vectors[0] * vector_xy,
495 : rotation_vectors[1] * vector_xy,
496 : rotation_vectors[2] * vector_xy);
497 : }
498 :
499 : is_node_created = true;
500 : }
501 2874 : else if (e == 0 && k == 0)
502 : {
503 : // On-axis nodes are only added once
504 : current_distance.zero();
505 : is_node_created = true;
506 : }
507 :
508 : if (is_node_created)
509 : {
510 0 : Node * new_node = mesh->add_point(*node + current_distance,
511 28474 : node->id() + (current_node_layer * orig_nodes),
512 28474 : node->processor_id());
513 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
514 : // Let's give the base nodes of the revolved mesh the same
515 : // unique_ids as the source mesh, in case anyone finds that
516 : // a useful map to preserve.
517 : const unique_id_type uid =
518 : (current_node_layer == 0)
519 28474 : ? node->unique_id()
520 25133 : : (orig_unique_ids + (current_node_layer - 1) * (orig_nodes + orig_elem * 2) +
521 25133 : node->id());
522 : new_node->set_unique_id(uid);
523 : #endif
524 :
525 28474 : input_boundary_info.boundary_ids(node, ids_to_copy);
526 28474 : if (_boundary_swap_pairs.empty())
527 25549 : boundary_info.add_node(new_node, ids_to_copy);
528 : else
529 2925 : for (const auto & id_to_copy : ids_to_copy)
530 0 : boundary_info.add_node(new_node,
531 : _boundary_swap_pairs[e].count(id_to_copy)
532 0 : ? _boundary_swap_pairs[e][id_to_copy]
533 0 : : id_to_copy);
534 : }
535 :
536 31074 : current_node_layer++;
537 : }
538 : }
539 201 : }
540 :
541 4190 : for (const auto & elem : input->element_ptr_range())
542 : {
543 1894 : const ElemType etype = elem->type();
544 :
545 : // revolving currently only works on coarse meshes
546 : mooseAssert(!elem->parent(), "RevolveGenerator only works on coarse meshes.");
547 :
548 : unsigned int current_layer = 0;
549 :
550 1894 : const unsigned int num_rotations = _revolving_angles.size();
551 :
552 4868 : for (unsigned int e = 0; e < num_rotations; e++)
553 : {
554 2974 : auto num_layers = _nums_azimuthal_intervals[e];
555 :
556 15634 : for (unsigned int k = 0; k < num_layers; ++k)
557 : {
558 12660 : std::unique_ptr<Elem> new_elem;
559 12660 : std::unique_ptr<Elem> new_elem_1;
560 12660 : bool is_flipped(false);
561 : // In some cases, two elements per layer are generated by revolving one element. So we
562 : // reserve an additional flag for the potential second element.
563 12660 : bool is_flipped_additional(false);
564 12660 : dof_id_type axis_node_case(-1);
565 : std::vector<std::pair<dof_id_type, dof_id_type>> side_pairs;
566 12660 : switch (etype)
567 : {
568 144 : case EDGE2:
569 : {
570 : // Possible scenarios:
571 : // 1. None of the nodes are on the axis
572 : // Then a quad4 element is created
573 : // 2. One of the nodes is on the axis
574 : // Then a tri3 element is created
575 144 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
576 144 : if (nodes_cates.first.empty())
577 : {
578 54 : createQUADfromEDGE(QUAD4,
579 : elem,
580 : mesh,
581 : new_elem,
582 : current_layer,
583 : orig_nodes,
584 : total_num_azimuthal_intervals,
585 : side_pairs,
586 : is_flipped);
587 : }
588 : else
589 : {
590 90 : createTRIfromEDGE(nodes_cates,
591 : TRI3,
592 : elem,
593 : mesh,
594 : new_elem,
595 : current_layer,
596 : orig_nodes,
597 : total_num_azimuthal_intervals,
598 : side_pairs,
599 : axis_node_case,
600 : is_flipped);
601 : }
602 : break;
603 : }
604 144 : case EDGE3:
605 : {
606 : // Possible scenarios:
607 : // 1. None of the nodes are on the axis
608 : // Then a QUAD9 element is created
609 : // 2. One of the nodes is on the axis
610 : // Then a TRI7 element is created
611 144 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
612 144 : if (nodes_cates.first.empty())
613 : {
614 54 : createQUADfromEDGE(QUAD9,
615 : elem,
616 : mesh,
617 : new_elem,
618 : current_layer,
619 : orig_nodes,
620 : total_num_azimuthal_intervals,
621 : side_pairs,
622 : is_flipped);
623 : }
624 : else
625 : {
626 90 : createTRIfromEDGE(nodes_cates,
627 : TRI7,
628 : elem,
629 : mesh,
630 : new_elem,
631 : current_layer,
632 : orig_nodes,
633 : total_num_azimuthal_intervals,
634 : side_pairs,
635 : axis_node_case,
636 : is_flipped);
637 : }
638 : break;
639 : }
640 3564 : case TRI3:
641 : {
642 : // Possible scenarios:
643 : // 1. None of the nodes are on the axis
644 : // Then a prism6 element is created
645 : // 2. One of the nodes is on the axis
646 : // Then a pyramid5 element is created
647 : // 3. Two of the nodes are on the axis
648 : // Then a tet4 element is created
649 3564 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
650 3564 : if (nodes_cates.first.empty())
651 : {
652 3294 : createPRISMfromTRI(PRISM6,
653 : elem,
654 : mesh,
655 : new_elem,
656 : current_layer,
657 : orig_nodes,
658 : total_num_azimuthal_intervals,
659 : side_pairs,
660 : is_flipped);
661 : }
662 270 : else if (nodes_cates.first.size() == 1)
663 : {
664 216 : createPYRAMIDfromTRI(nodes_cates,
665 : PYRAMID5,
666 : elem,
667 : mesh,
668 : new_elem,
669 : current_layer,
670 : orig_nodes,
671 : total_num_azimuthal_intervals,
672 : side_pairs,
673 : axis_node_case,
674 : is_flipped);
675 : }
676 54 : else if (nodes_cates.first.size() == 2)
677 : {
678 54 : createTETfromTRI(nodes_cates,
679 : TET4,
680 : elem,
681 : mesh,
682 : new_elem,
683 : current_layer,
684 : orig_nodes,
685 : total_num_azimuthal_intervals,
686 : side_pairs,
687 : axis_node_case,
688 : is_flipped);
689 : }
690 : else
691 0 : mooseError("A degenerate TRI3 elements overlapped with the rotation axis cannot be "
692 : "revolved.");
693 :
694 : break;
695 : }
696 408 : case TRI6:
697 : {
698 : // Possible scenarios:
699 : // 1. None of the nodes are on the axis
700 : // Then a prism18 element is created
701 : // 2. One of the nodes is on the axis
702 : // Then a pyramid13 element is created
703 : // 3. Three of the nodes are on the axis
704 : // Then a tet10 element is created
705 : // NOTE: We do not support two nodes on the axis for tri6 elements
706 408 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
707 408 : if (nodes_cates.first.empty())
708 : {
709 204 : createPRISMfromTRI(PRISM18,
710 : elem,
711 : mesh,
712 : new_elem,
713 : current_layer,
714 : orig_nodes,
715 : total_num_azimuthal_intervals,
716 : side_pairs,
717 : is_flipped);
718 : }
719 204 : else if (nodes_cates.first.size() == 1)
720 : {
721 108 : createPYRAMIDfromTRI(nodes_cates,
722 : PYRAMID13,
723 : elem,
724 : mesh,
725 : new_elem,
726 : current_layer,
727 : orig_nodes,
728 : total_num_azimuthal_intervals,
729 : side_pairs,
730 : axis_node_case,
731 : is_flipped);
732 : }
733 96 : else if (nodes_cates.first.size() == 3)
734 : {
735 96 : createTETfromTRI(nodes_cates,
736 : TET10,
737 : elem,
738 : mesh,
739 : new_elem,
740 : current_layer,
741 : orig_nodes,
742 : total_num_azimuthal_intervals,
743 : side_pairs,
744 : axis_node_case,
745 : is_flipped);
746 : }
747 : else
748 0 : mooseError(
749 : "You either have a degenerate TRI6 element, or the mid-point of the "
750 : "on-axis edge is not colinear with the two vertices, which is not supported.");
751 : break;
752 : }
753 408 : case TRI7:
754 : {
755 : // Possible scenarios:
756 : // 1. None of the nodes are on the axis
757 : // Then a prism21 element is created
758 : // 2. One of the nodes is on the axis
759 : // Then a pyramid18 element is created
760 : // 3. Three of the nodes are on the axis
761 : // Then a tet14 element is created
762 : // NOTE: We do not support two nodes on the axis for tri7 elements
763 408 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
764 408 : if (nodes_cates.first.empty())
765 : {
766 204 : createPRISMfromTRI(PRISM21,
767 : elem,
768 : mesh,
769 : new_elem,
770 : current_layer,
771 : orig_nodes,
772 : total_num_azimuthal_intervals,
773 : side_pairs,
774 : is_flipped);
775 : }
776 204 : else if (nodes_cates.first.size() == 1)
777 : {
778 108 : createPYRAMIDfromTRI(nodes_cates,
779 : PYRAMID18,
780 : elem,
781 : mesh,
782 : new_elem,
783 : current_layer,
784 : orig_nodes,
785 : total_num_azimuthal_intervals,
786 : side_pairs,
787 : axis_node_case,
788 : is_flipped);
789 : }
790 96 : else if (nodes_cates.first.size() == 3)
791 : {
792 96 : createTETfromTRI(nodes_cates,
793 : TET14,
794 : elem,
795 : mesh,
796 : new_elem,
797 : current_layer,
798 : orig_nodes,
799 : total_num_azimuthal_intervals,
800 : side_pairs,
801 : axis_node_case,
802 : is_flipped);
803 : }
804 : else
805 0 : mooseError("You either have a degenerate TRI6 element, or the mid-point of the "
806 : "on-axis edge of the TRI6 element is not colinear with the two vertices, "
807 : "which is not supported.");
808 : break;
809 : }
810 7128 : case QUAD4:
811 : {
812 : // Possible scenarios:
813 : // 1. None of the nodes are on the axis
814 : // Then a hex8 element is created
815 : // 2. One of the nodes is on the axis
816 : // Then a pyramid5 element and a prism6 element are created
817 : // 3. Two of the nodes are on the axis
818 : // Then a prism6 is created
819 7128 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
820 7128 : if (nodes_cates.first.empty())
821 : {
822 6966 : createHEXfromQUAD(HEX8,
823 : elem,
824 : mesh,
825 : new_elem,
826 : current_layer,
827 : orig_nodes,
828 : total_num_azimuthal_intervals,
829 : side_pairs,
830 : is_flipped);
831 : }
832 162 : else if (nodes_cates.first.size() == 1)
833 : {
834 108 : createPYRAMIDPRISMfromQUAD(nodes_cates,
835 : PYRAMID5,
836 : PRISM6,
837 : elem,
838 : mesh,
839 : new_elem,
840 : new_elem_1,
841 : current_layer,
842 : orig_nodes,
843 : total_num_azimuthal_intervals,
844 : side_pairs,
845 : axis_node_case,
846 : is_flipped,
847 : is_flipped_additional);
848 : }
849 54 : else if (nodes_cates.first.size() == 2)
850 : {
851 54 : createPRISMfromQUAD(nodes_cates,
852 : PRISM6,
853 : elem,
854 : mesh,
855 : new_elem,
856 : current_layer,
857 : orig_nodes,
858 : total_num_azimuthal_intervals,
859 : side_pairs,
860 : axis_node_case,
861 : is_flipped);
862 : }
863 :
864 : else
865 0 : mooseError("Degenerate QUAD4 element with 3 or more aligned nodes cannot be "
866 : "azimuthally revolved");
867 :
868 : break;
869 : }
870 378 : case QUAD8:
871 : {
872 : // Possible scenarios:
873 : // 1. None of the nodes are on the axis
874 : // Then a hex20 element is created
875 : // 2. One of the nodes is on the axis
876 : // In that case, it is already converted to a QUAD9 element before,
877 : // SO we do not need to worry about this case
878 : // 3. Three of the nodes are on the axis
879 : // Then a prism15 is created
880 : // NOTE: We do not support two nodes on the axis for quad8 elements
881 378 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
882 378 : if (nodes_cates.first.empty())
883 : {
884 270 : createHEXfromQUAD(HEX20,
885 : elem,
886 : mesh,
887 : new_elem,
888 : current_layer,
889 : orig_nodes,
890 : total_num_azimuthal_intervals,
891 : side_pairs,
892 : is_flipped);
893 : }
894 108 : else if (nodes_cates.first.size() == 3)
895 : {
896 108 : createPRISMfromQUAD(nodes_cates,
897 : PRISM15,
898 : elem,
899 : mesh,
900 : new_elem,
901 : current_layer,
902 : orig_nodes,
903 : total_num_azimuthal_intervals,
904 : side_pairs,
905 : axis_node_case,
906 : is_flipped);
907 : }
908 : else
909 0 : mooseError("You either have a degenerate QUAD8 element, or the mid-point of the "
910 : "on-axis edge of the QUAD8 element is not colinear with the two vertices, "
911 : "which is not supported.");
912 :
913 : break;
914 : }
915 486 : case QUAD9:
916 : {
917 : // Possible scenarios:
918 : // 1. None of the nodes are on the axis
919 : // Then a hex27 element is created
920 : // 2. One of the nodes is on the axis
921 : // Then a pyramid14 element and a prism18 element are created
922 : // 3. Two of the nodes are on the axis
923 : // Then a prism18 is created
924 : // (we do not create prism20/21 here just to make prism18 the only possible prism
925 : // elements for simplicity)
926 486 : const auto nodes_cates = onAxisNodesIdentifier(*elem, nodes_on_axis);
927 486 : if (nodes_cates.first.empty())
928 : {
929 270 : createHEXfromQUAD(HEX27,
930 : elem,
931 : mesh,
932 : new_elem,
933 : current_layer,
934 : orig_nodes,
935 : total_num_azimuthal_intervals,
936 : side_pairs,
937 : is_flipped);
938 : }
939 216 : else if (nodes_cates.first.size() == 1)
940 : {
941 108 : createPYRAMIDPRISMfromQUAD(nodes_cates,
942 : PYRAMID14,
943 : PRISM18,
944 : elem,
945 : mesh,
946 : new_elem,
947 : new_elem_1,
948 : current_layer,
949 : orig_nodes,
950 : total_num_azimuthal_intervals,
951 : side_pairs,
952 : axis_node_case,
953 : is_flipped,
954 : is_flipped_additional);
955 : }
956 108 : else if (nodes_cates.first.size() == 3)
957 : {
958 108 : createPRISMfromQUAD(nodes_cates,
959 : PRISM18,
960 : elem,
961 : mesh,
962 : new_elem,
963 : current_layer,
964 : orig_nodes,
965 : total_num_azimuthal_intervals,
966 : side_pairs,
967 : axis_node_case,
968 : is_flipped);
969 : }
970 : else
971 0 : mooseError("You either have a degenerate QUAD9 element, or the mid-point of the "
972 : "on-axis edge of the QUAD9 element is not colinear with the two vertices, "
973 : "which is not supported.");
974 : break;
975 : }
976 0 : default:
977 0 : mooseError("The input mesh contains unsupported element type(s).");
978 : }
979 12660 : new_elem->set_id(elem->id() + (current_layer * orig_elem));
980 12660 : new_elem->processor_id() = elem->processor_id();
981 12660 : if (new_elem_1)
982 : {
983 216 : new_elem_1->set_id(elem->id() + (current_layer * orig_elem) + elem_id_shift);
984 216 : new_elem_1->processor_id() = elem->processor_id();
985 : }
986 :
987 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
988 : // Let's give the base elements of the revolved mesh the same
989 : // unique_ids as the source mesh, in case anyone finds that
990 : // a useful map to preserve.
991 : const unique_id_type uid =
992 : (current_layer == 0)
993 12660 : ? elem->unique_id()
994 10766 : : (orig_unique_ids + (current_layer - 1) * (orig_nodes + orig_elem * 2) +
995 : orig_nodes + elem->id());
996 :
997 : new_elem->set_unique_id(uid);
998 :
999 : // Special case for extra elements
1000 12660 : if (new_elem_1)
1001 : {
1002 : const unique_id_type uid_1 =
1003 : (current_layer == 0)
1004 216 : ? (elem->id() + orig_unique_ids - orig_elem)
1005 180 : : (orig_unique_ids + (current_layer - 1) * (orig_nodes + orig_elem * 2) +
1006 180 : orig_nodes + orig_elem + elem->id());
1007 :
1008 : new_elem_1->set_unique_id(uid_1);
1009 : }
1010 : #endif
1011 :
1012 : // maintain the subdomain_id
1013 12660 : switch (etype)
1014 : {
1015 : case EDGE2:
1016 144 : switch (new_elem->type())
1017 : {
1018 54 : case QUAD4:
1019 54 : new_elem->subdomain_id() = elem->subdomain_id();
1020 54 : break;
1021 90 : case TRI3:
1022 90 : new_elem->subdomain_id() = edge_to_tri_subdomain_id_shift + elem->subdomain_id();
1023 90 : break;
1024 : default:
1025 : mooseAssert(false,
1026 : "impossible element type generated by revolving an EDGE2 element");
1027 : }
1028 : break;
1029 : case EDGE3:
1030 144 : switch (new_elem->type())
1031 : {
1032 54 : case QUAD9:
1033 54 : new_elem->subdomain_id() = elem->subdomain_id();
1034 54 : break;
1035 90 : case TRI7:
1036 90 : new_elem->subdomain_id() = edge_to_tri_subdomain_id_shift + elem->subdomain_id();
1037 90 : break;
1038 : default:
1039 : mooseAssert(false,
1040 : "impossible element type generated by revolving an EDGE3 element");
1041 : }
1042 : break;
1043 : case TRI3:
1044 3564 : switch (new_elem->type())
1045 : {
1046 3294 : case PRISM6:
1047 3294 : new_elem->subdomain_id() = elem->subdomain_id();
1048 3294 : break;
1049 216 : case PYRAMID5:
1050 216 : new_elem->subdomain_id() = tri_to_pyramid_subdomain_id_shift + elem->subdomain_id();
1051 216 : break;
1052 54 : case TET4:
1053 54 : new_elem->subdomain_id() = tri_to_tet_subdomain_id_shift + elem->subdomain_id();
1054 54 : break;
1055 : default:
1056 : mooseAssert(false, "impossible element type generated by revolving a TRI3 element");
1057 : }
1058 : break;
1059 : case TRI6:
1060 408 : switch (new_elem->type())
1061 : {
1062 204 : case PRISM18:
1063 204 : new_elem->subdomain_id() = elem->subdomain_id();
1064 204 : break;
1065 108 : case PYRAMID13:
1066 108 : new_elem->subdomain_id() = tri_to_pyramid_subdomain_id_shift + elem->subdomain_id();
1067 108 : break;
1068 96 : case TET10:
1069 96 : new_elem->subdomain_id() = tri_to_tet_subdomain_id_shift + elem->subdomain_id();
1070 96 : break;
1071 : default:
1072 : mooseAssert(false, "impossible element type generated by revolving a TRI6 element");
1073 : }
1074 : break;
1075 : case TRI7:
1076 408 : switch (new_elem->type())
1077 : {
1078 204 : case PRISM21:
1079 204 : new_elem->subdomain_id() = elem->subdomain_id();
1080 204 : break;
1081 108 : case PYRAMID18:
1082 108 : new_elem->subdomain_id() = tri_to_pyramid_subdomain_id_shift + elem->subdomain_id();
1083 108 : break;
1084 96 : case TET14:
1085 96 : new_elem->subdomain_id() = tri_to_tet_subdomain_id_shift + elem->subdomain_id();
1086 96 : break;
1087 : default:
1088 : mooseAssert(false, "impossible element type generated by revolving a TRI7 element");
1089 : }
1090 : break;
1091 : case QUAD4:
1092 7128 : switch (new_elem->type())
1093 : {
1094 6966 : case HEX8:
1095 6966 : new_elem->subdomain_id() = elem->subdomain_id();
1096 6966 : break;
1097 54 : case PRISM6:
1098 54 : new_elem->subdomain_id() = quad_to_prism_subdomain_id_shift + elem->subdomain_id();
1099 54 : break;
1100 108 : case PYRAMID5:
1101 108 : new_elem->subdomain_id() =
1102 108 : quad_to_pyramid_subdomain_id_shift + elem->subdomain_id();
1103 108 : new_elem_1->subdomain_id() =
1104 108 : quad_to_prism_subdomain_id_shift + elem->subdomain_id();
1105 108 : break;
1106 : default:
1107 : mooseAssert(false,
1108 : "impossible element type generated by revolving a QUAD4 element");
1109 : }
1110 : break;
1111 : case QUAD8:
1112 378 : switch (new_elem->type())
1113 : {
1114 270 : case HEX20:
1115 270 : new_elem->subdomain_id() = elem->subdomain_id();
1116 270 : break;
1117 108 : case PRISM15:
1118 108 : new_elem->subdomain_id() = quad_to_prism_subdomain_id_shift + elem->subdomain_id();
1119 108 : break;
1120 : default:
1121 : mooseAssert(false,
1122 : "impossible element type generated by revolving a QUAD8 element");
1123 : }
1124 : break;
1125 : case QUAD9:
1126 486 : switch (new_elem->type())
1127 : {
1128 270 : case HEX27:
1129 270 : new_elem->subdomain_id() = elem->subdomain_id();
1130 270 : break;
1131 108 : case PRISM18:
1132 108 : new_elem->subdomain_id() =
1133 108 : quad_to_hi_pyramid_subdomain_id_shift + elem->subdomain_id();
1134 108 : break;
1135 108 : case PYRAMID14:
1136 108 : new_elem->subdomain_id() =
1137 108 : quad_to_pyramid_subdomain_id_shift + elem->subdomain_id();
1138 108 : new_elem_1->subdomain_id() =
1139 108 : quad_to_hi_pyramid_subdomain_id_shift + elem->subdomain_id();
1140 108 : break;
1141 : default:
1142 : mooseAssert(false,
1143 : "impossible element type generated by revolving a QUAD9 element");
1144 : }
1145 : break;
1146 12660 : default:
1147 : mooseAssert(false,
1148 : "The input mesh contains unsupported element type(s), which should have "
1149 : "been checked in prior steps in this code.");
1150 : }
1151 :
1152 12660 : if (_subdomain_swap_pairs.size())
1153 : {
1154 : auto & revolving_swap_pairs = _subdomain_swap_pairs[e];
1155 :
1156 2592 : auto new_id_it = revolving_swap_pairs.find(elem->subdomain_id());
1157 :
1158 2592 : if (new_id_it != revolving_swap_pairs.end())
1159 : {
1160 1296 : new_elem->subdomain_id() =
1161 1296 : new_elem->subdomain_id() - elem->subdomain_id() + new_id_it->second;
1162 1296 : if (new_elem_1)
1163 0 : new_elem_1->subdomain_id() =
1164 0 : new_elem_1->subdomain_id() - elem->subdomain_id() + new_id_it->second;
1165 : }
1166 : }
1167 :
1168 12660 : Elem * added_elem = mesh->add_elem(std::move(new_elem));
1169 : Elem * added_elem_1 = NULL;
1170 :
1171 12660 : if (new_elem_1)
1172 216 : added_elem_1 = mesh->add_elem(std::move(new_elem_1));
1173 :
1174 : // maintain extra integers
1175 23028 : for (unsigned int i = 0; i < num_extra_elem_integers; i++)
1176 : {
1177 10368 : added_elem->set_extra_integer(i, elem->get_extra_integer(i));
1178 10368 : if (added_elem_1)
1179 0 : added_elem_1->set_extra_integer(i, elem->get_extra_integer(i));
1180 : }
1181 :
1182 12660 : if (_elem_integers_swap_pairs.size())
1183 : {
1184 15552 : for (unsigned int i = 0; i < _elem_integer_indices_to_swap.size(); i++)
1185 : {
1186 : auto & elevation_extra_swap_pairs =
1187 10368 : _elem_integers_swap_pairs[i * _nums_azimuthal_intervals.size() + e];
1188 :
1189 : auto new_extra_id_it = elevation_extra_swap_pairs.find(
1190 10368 : elem->get_extra_integer(_elem_integer_indices_to_swap[i]));
1191 :
1192 10368 : if (new_extra_id_it != elevation_extra_swap_pairs.end())
1193 : {
1194 6912 : added_elem->set_extra_integer(_elem_integer_indices_to_swap[i],
1195 : new_extra_id_it->second);
1196 6912 : if (added_elem_1)
1197 0 : added_elem_1->set_extra_integer(_elem_integer_indices_to_swap[i],
1198 : new_extra_id_it->second);
1199 : }
1200 : }
1201 : }
1202 :
1203 : // Copy any old boundary ids on all sides
1204 58344 : for (auto s : elem->side_index_range())
1205 : {
1206 45684 : input_boundary_info.boundary_ids(elem, s, ids_to_copy);
1207 : std::vector<boundary_id_type> ids_to_copy_swapped;
1208 45684 : if (_boundary_swap_pairs.empty())
1209 36612 : ids_to_copy_swapped = ids_to_copy;
1210 : else
1211 11664 : for (const auto & id_to_copy : ids_to_copy)
1212 2592 : ids_to_copy_swapped.push_back(_boundary_swap_pairs[e].count(id_to_copy)
1213 : ? _boundary_swap_pairs[e][id_to_copy]
1214 : : id_to_copy);
1215 :
1216 45684 : switch (etype)
1217 : {
1218 288 : case EDGE2:
1219 288 : switch (added_elem->type())
1220 : {
1221 108 : case QUAD4:
1222 108 : boundary_info.add_side(
1223 108 : added_elem, cast_int<unsigned short>(s == 0 ? 3 : 1), ids_to_copy_swapped);
1224 : break;
1225 180 : case TRI3:
1226 180 : if (s != axis_node_case)
1227 90 : boundary_info.add_side(
1228 : added_elem, cast_int<unsigned short>(s), ids_to_copy_swapped);
1229 : break;
1230 : default:
1231 : mooseAssert(false,
1232 : "impossible element type generated by revolving an EDGE2 element");
1233 : }
1234 : break;
1235 288 : case EDGE3:
1236 288 : switch (added_elem->type())
1237 : {
1238 108 : case QUAD9:
1239 108 : boundary_info.add_side(
1240 108 : added_elem, cast_int<unsigned short>(s == 0 ? 3 : 1), ids_to_copy_swapped);
1241 : break;
1242 180 : case TRI7:
1243 180 : if (s != axis_node_case)
1244 90 : boundary_info.add_side(
1245 : added_elem, cast_int<unsigned short>(s), ids_to_copy_swapped);
1246 : break;
1247 : default:
1248 : mooseAssert(false,
1249 : "impossible element type generated by revolving an EDGE3 element");
1250 : }
1251 : break;
1252 10692 : case TRI3:
1253 10692 : switch (added_elem->type())
1254 : {
1255 9882 : case PRISM6:
1256 9882 : boundary_info.add_side(
1257 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1258 : break;
1259 648 : case PYRAMID5:
1260 648 : if ((s + 3 - axis_node_case) % 3 == 0)
1261 216 : boundary_info.add_side(
1262 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1263 432 : else if ((s + 3 - axis_node_case) % 3 == 1)
1264 216 : boundary_info.add_side(
1265 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1266 : else
1267 216 : boundary_info.add_side(
1268 : added_elem, cast_int<unsigned short>(1), ids_to_copy_swapped);
1269 : break;
1270 162 : case TET4:
1271 162 : if ((s + 3 - axis_node_case) % 3 == 0)
1272 54 : boundary_info.add_side(
1273 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1274 108 : else if ((s + 3 - axis_node_case) % 3 == 2)
1275 54 : boundary_info.add_side(
1276 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1277 : break;
1278 : default:
1279 : mooseAssert(false,
1280 : "impossible element type generated by revolving a TRI3 element");
1281 : }
1282 : break;
1283 1224 : case TRI6:
1284 1224 : switch (added_elem->type())
1285 : {
1286 612 : case PRISM18:
1287 612 : boundary_info.add_side(
1288 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1289 : break;
1290 324 : case PYRAMID13:
1291 324 : if ((s + 3 - axis_node_case) % 3 == 0)
1292 108 : boundary_info.add_side(
1293 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1294 216 : else if ((s + 3 - axis_node_case) % 3 == 1)
1295 108 : boundary_info.add_side(
1296 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1297 : else
1298 108 : boundary_info.add_side(
1299 : added_elem, cast_int<unsigned short>(1), ids_to_copy_swapped);
1300 : break;
1301 288 : case TET10:
1302 288 : if ((s + 3 - axis_node_case) % 3 == 0)
1303 96 : boundary_info.add_side(
1304 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1305 192 : else if ((s + 3 - axis_node_case) % 3 == 2)
1306 96 : boundary_info.add_side(
1307 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1308 : break;
1309 : default:
1310 : mooseAssert(false,
1311 : "impossible element type generated by revolving a TRI6 element");
1312 : }
1313 : break;
1314 1224 : case TRI7:
1315 1224 : switch (added_elem->type())
1316 : {
1317 612 : case PRISM21:
1318 612 : boundary_info.add_side(
1319 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1320 : break;
1321 324 : case PYRAMID18:
1322 324 : if ((s + 3 - axis_node_case) % 3 == 0)
1323 108 : boundary_info.add_side(
1324 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1325 216 : else if ((s + 3 - axis_node_case) % 3 == 1)
1326 108 : boundary_info.add_side(
1327 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1328 : else
1329 108 : boundary_info.add_side(
1330 : added_elem, cast_int<unsigned short>(1), ids_to_copy_swapped);
1331 : break;
1332 288 : case TET14:
1333 288 : if ((s + 3 - axis_node_case) % 3 == 0)
1334 96 : boundary_info.add_side(
1335 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1336 192 : else if ((s + 3 - axis_node_case) % 3 == 2)
1337 96 : boundary_info.add_side(
1338 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1339 : break;
1340 : default:
1341 : mooseAssert(false,
1342 : "impossible element type generated by revolving a TRI7 element");
1343 : }
1344 : break;
1345 28512 : case QUAD4:
1346 28512 : switch (added_elem->type())
1347 : {
1348 27864 : case HEX8:
1349 27864 : boundary_info.add_side(
1350 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1351 : break;
1352 216 : case PRISM6:
1353 216 : if ((s + 4 - axis_node_case) % 4 == 1)
1354 54 : boundary_info.add_side(
1355 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1356 162 : else if ((s + 4 - axis_node_case) % 4 == 2)
1357 54 : boundary_info.add_side(
1358 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1359 108 : else if ((s + 4 - axis_node_case) % 4 == 3)
1360 54 : boundary_info.add_side(
1361 : added_elem, cast_int<unsigned short>(0), ids_to_copy_swapped);
1362 : break;
1363 432 : case PYRAMID5:
1364 432 : if ((s + 4 - axis_node_case) % 4 == 3)
1365 108 : boundary_info.add_side(
1366 : added_elem, cast_int<unsigned short>(1), ids_to_copy_swapped);
1367 324 : else if ((s + 4 - axis_node_case) % 4 == 0)
1368 108 : boundary_info.add_side(
1369 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1370 216 : else if ((s + 4 - axis_node_case) % 4 == 1)
1371 108 : boundary_info.add_side(
1372 : added_elem_1, cast_int<unsigned short>(3), ids_to_copy_swapped);
1373 : else
1374 108 : boundary_info.add_side(
1375 : added_elem_1, cast_int<unsigned short>(2), ids_to_copy_swapped);
1376 : break;
1377 : default:
1378 : mooseAssert(false,
1379 : "impossible element type generated by revolving a QUAD4 element");
1380 : }
1381 : break;
1382 1512 : case QUAD8:
1383 1512 : switch (added_elem->type())
1384 : {
1385 1080 : case HEX20:
1386 1080 : boundary_info.add_side(
1387 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1388 : break;
1389 432 : case PRISM15:
1390 432 : if ((s + 4 - axis_node_case) % 4 == 1)
1391 108 : boundary_info.add_side(
1392 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1393 324 : else if ((s + 4 - axis_node_case) % 4 == 2)
1394 108 : boundary_info.add_side(
1395 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1396 216 : else if ((s + 4 - axis_node_case) % 4 == 3)
1397 108 : boundary_info.add_side(
1398 : added_elem, cast_int<unsigned short>(0), ids_to_copy_swapped);
1399 : break;
1400 : default:
1401 : mooseAssert(false,
1402 : "impossible element type generated by revolving a QUAD8 element");
1403 : }
1404 : break;
1405 1944 : case QUAD9:
1406 1944 : switch (added_elem->type())
1407 : {
1408 1080 : case HEX27:
1409 1080 : boundary_info.add_side(
1410 : added_elem, cast_int<unsigned short>(s + 1), ids_to_copy_swapped);
1411 : break;
1412 432 : case PRISM18:
1413 432 : if ((s + 4 - axis_node_case) % 4 == 1)
1414 108 : boundary_info.add_side(
1415 : added_elem, cast_int<unsigned short>(4), ids_to_copy_swapped);
1416 324 : else if ((s + 4 - axis_node_case) % 4 == 2)
1417 108 : boundary_info.add_side(
1418 : added_elem, cast_int<unsigned short>(2), ids_to_copy_swapped);
1419 216 : else if ((s + 4 - axis_node_case) % 4 == 3)
1420 108 : boundary_info.add_side(
1421 : added_elem, cast_int<unsigned short>(0), ids_to_copy_swapped);
1422 : break;
1423 432 : case PYRAMID14:
1424 432 : if ((s + 4 - axis_node_case) % 4 == 3)
1425 108 : boundary_info.add_side(
1426 : added_elem, cast_int<unsigned short>(1), ids_to_copy_swapped);
1427 324 : else if ((s + 4 - axis_node_case) % 4 == 0)
1428 108 : boundary_info.add_side(
1429 : added_elem, cast_int<unsigned short>(3), ids_to_copy_swapped);
1430 216 : else if ((s + 4 - axis_node_case) % 4 == 1)
1431 108 : boundary_info.add_side(
1432 : added_elem_1, cast_int<unsigned short>(3), ids_to_copy_swapped);
1433 : else
1434 108 : boundary_info.add_side(
1435 : added_elem_1, cast_int<unsigned short>(2), ids_to_copy_swapped);
1436 : break;
1437 : default:
1438 : mooseAssert(false,
1439 : "impossible element type generated by revolving a QUAD9 element");
1440 : }
1441 : break;
1442 45684 : default:
1443 : mooseAssert(false,
1444 : "The input mesh contains unsupported element type(s), which should have "
1445 : "been checked in prior steps in this code.");
1446 : }
1447 45684 : }
1448 :
1449 12660 : if (current_layer == 0 && _has_start_boundary)
1450 : {
1451 54 : boundary_info.add_side(
1452 54 : added_elem, is_flipped ? side_pairs[0].second : side_pairs[0].first, _start_boundary);
1453 54 : if (side_pairs.size() > 1)
1454 0 : boundary_info.add_side(added_elem_1,
1455 0 : is_flipped_additional ? side_pairs[1].second
1456 0 : : side_pairs[1].first,
1457 0 : _start_boundary);
1458 : }
1459 :
1460 12660 : if (current_layer == num_layers - 1 && _has_end_boundary)
1461 : {
1462 54 : boundary_info.add_side(
1463 54 : added_elem, is_flipped ? side_pairs[0].first : side_pairs[0].second, _end_boundary);
1464 54 : if (side_pairs.size() > 1)
1465 0 : boundary_info.add_side(added_elem_1,
1466 0 : is_flipped_additional ? side_pairs[1].first
1467 0 : : side_pairs[1].second,
1468 0 : _end_boundary);
1469 : }
1470 12660 : current_layer++;
1471 12660 : }
1472 : }
1473 201 : }
1474 :
1475 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1476 : // Update the value of next_unique_id based on newly created nodes and elements
1477 : // Note: the calculation here is quite conservative to ensure uniqueness
1478 201 : unsigned int total_new_node_layers = total_num_azimuthal_intervals * order;
1479 201 : unsigned int new_unique_ids = orig_unique_ids + (total_new_node_layers - 1) * orig_elem * 2 +
1480 : total_new_node_layers * orig_nodes;
1481 201 : mesh->set_next_unique_id(new_unique_ids);
1482 : #endif
1483 :
1484 : // Copy all the subdomain/sideset/nodeset name maps to the revolved mesh
1485 201 : if (!input_subdomain_map.empty())
1486 108 : mesh->set_subdomain_name_map().insert(input_subdomain_map.begin(), input_subdomain_map.end());
1487 201 : if (!input_sideset_map.empty())
1488 201 : mesh->get_boundary_info().set_sideset_name_map().insert(input_sideset_map.begin(),
1489 : input_sideset_map.end());
1490 201 : if (!input_nodeset_map.empty())
1491 84 : mesh->get_boundary_info().set_nodeset_name_map().insert(input_nodeset_map.begin(),
1492 : input_nodeset_map.end());
1493 :
1494 201 : mesh->remove_orphaned_nodes();
1495 201 : mesh->renumber_nodes_and_elements();
1496 : mesh->set_isnt_prepared();
1497 :
1498 201 : return mesh;
1499 603 : }
1500 :
1501 : std::pair<Real, Point>
1502 3341 : RevolveGenerator::getRotationCenterAndRadius(const Point & p_ext,
1503 : const Point & p_axis,
1504 : const Point & dir_axis) const
1505 : {
1506 : // First use point product to get the distance between the axis point and the projection of the
1507 : // external point on the axis
1508 3341 : const Real dist = (p_ext - p_axis) * dir_axis.unit();
1509 3341 : const Point center_pt = p_axis + dist * dir_axis.unit();
1510 : // Then get the radius
1511 3341 : const Real radius = (p_ext - center_pt).norm();
1512 3341 : return std::make_pair(radius, center_pt);
1513 : }
1514 :
1515 : std::vector<Point>
1516 201 : RevolveGenerator::rotationVectors(const Point & p_axis,
1517 : const Point & dir_axis,
1518 : const Point & p_input) const
1519 : {
1520 : // To make the rotation mathematically simple, we perform rotation in a coordination system
1521 : // (x',y',z') defined by rotation axis and the mesh to be rotated.
1522 : // z' is the rotation axis, which is trivial dir_axis.unit()
1523 201 : const Point z_prime = dir_axis.unit();
1524 : // the x' and z' should form the plane that accommodates input mesh
1525 201 : const Point x_prime = ((p_input - p_axis) - ((p_input - p_axis) * z_prime) * z_prime).unit();
1526 201 : const Point y_prime = z_prime.cross(x_prime);
1527 : // Then we transform things back to the original coordination system (x,y,z), which is trivial
1528 : // (1,0,0), (0,1,0), (0,0,1)
1529 : return {{x_prime(0), y_prime(0), z_prime(0)},
1530 : {x_prime(1), y_prime(1), z_prime(1)},
1531 201 : {x_prime(2), y_prime(2), z_prime(2)}};
1532 : }
1533 :
1534 : std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>>
1535 12660 : RevolveGenerator::onAxisNodesIdentifier(const Elem & elem,
1536 : const std::vector<dof_id_type> & nodes_on_axis) const
1537 : {
1538 : std::vector<dof_id_type> nodes_on_axis_in_elem;
1539 : std::vector<dof_id_type> nodes_not_on_axis_in_elem;
1540 65286 : for (unsigned int i = 0; i < elem.n_nodes(); i++)
1541 : {
1542 52626 : const auto node_id = elem.node_id(i);
1543 52626 : if (std::find(nodes_on_axis.begin(), nodes_on_axis.end(), node_id) != nodes_on_axis.end())
1544 : {
1545 2268 : nodes_on_axis_in_elem.push_back(i);
1546 : }
1547 : else
1548 : {
1549 50358 : nodes_not_on_axis_in_elem.push_back(i);
1550 : }
1551 : }
1552 12660 : return std::make_pair(nodes_on_axis_in_elem, nodes_not_on_axis_in_elem);
1553 12660 : }
1554 :
1555 : void
1556 126 : RevolveGenerator::nodeModification(Node & node)
1557 : {
1558 : const Point axis_component =
1559 126 : ((node - _axis_point) * _axis_direction.unit()) * _axis_direction.unit();
1560 : const Point rad_component = ((node - _axis_point) - axis_component) * _radius_correction_factor;
1561 : node = _axis_point + axis_component + rad_component;
1562 126 : }
1563 :
1564 : void
1565 108 : RevolveGenerator::createQUADfromEDGE(const ElemType quad_elem_type,
1566 : const Elem * elem,
1567 : const std::unique_ptr<MeshBase> & mesh,
1568 : std::unique_ptr<Elem> & new_elem,
1569 : const int current_layer,
1570 : const unsigned int orig_nodes,
1571 : const unsigned int total_num_azimuthal_intervals,
1572 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
1573 : bool & is_flipped) const
1574 : {
1575 108 : if (quad_elem_type != QUAD4 && quad_elem_type != QUAD9)
1576 0 : mooseError("Unsupported element type", quad_elem_type);
1577 :
1578 108 : side_pairs.push_back(std::make_pair(0, 2));
1579 108 : const unsigned int order = quad_elem_type == QUAD4 ? 1 : 2;
1580 :
1581 108 : new_elem = std::make_unique<Quad4>();
1582 108 : if (quad_elem_type == QUAD9)
1583 : {
1584 54 : new_elem = std::make_unique<Quad9>();
1585 54 : new_elem->set_node(4,
1586 54 : mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * 2 * orig_nodes)));
1587 54 : new_elem->set_node(
1588 54 : 5, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1589 54 : new_elem->set_node(
1590 : 6,
1591 54 : mesh->node_ptr(elem->node_ptr(2)->id() +
1592 54 : ((current_layer + 1) %
1593 54 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1594 54 : 2 * orig_nodes)));
1595 54 : new_elem->set_node(
1596 54 : 7, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1597 54 : new_elem->set_node(
1598 54 : 8, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1599 : }
1600 :
1601 108 : new_elem->set_node(
1602 108 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * order * orig_nodes)));
1603 108 : new_elem->set_node(
1604 108 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * order * orig_nodes)));
1605 108 : new_elem->set_node(
1606 : 3,
1607 108 : mesh->node_ptr(elem->node_ptr(0)->id() +
1608 108 : ((current_layer + 1) %
1609 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1610 108 : order * orig_nodes)));
1611 108 : new_elem->set_node(
1612 : 2,
1613 108 : mesh->node_ptr(elem->node_ptr(1)->id() +
1614 108 : ((current_layer + 1) %
1615 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1616 108 : order * orig_nodes)));
1617 :
1618 108 : if (new_elem->volume() < 0.0)
1619 : {
1620 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
1621 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 2);
1622 0 : if (quad_elem_type == QUAD9)
1623 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 4, 6);
1624 0 : is_flipped = true;
1625 : }
1626 108 : }
1627 :
1628 : void
1629 180 : RevolveGenerator::createTRIfromEDGE(
1630 : const std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>> & nodes_cates,
1631 : const ElemType tri_elem_type,
1632 : const Elem * elem,
1633 : const std::unique_ptr<MeshBase> & mesh,
1634 : std::unique_ptr<Elem> & new_elem,
1635 : const int current_layer,
1636 : const unsigned int orig_nodes,
1637 : const unsigned int total_num_azimuthal_intervals,
1638 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
1639 : dof_id_type & axis_node_case,
1640 : bool & is_flipped) const
1641 : {
1642 180 : if (tri_elem_type != TRI3 && tri_elem_type != TRI7)
1643 0 : mooseError("Unsupported element type", tri_elem_type);
1644 :
1645 180 : side_pairs.push_back(std::make_pair(0, 2));
1646 180 : const unsigned int order = tri_elem_type == TRI3 ? 1 : 2;
1647 180 : axis_node_case = nodes_cates.first.front();
1648 :
1649 180 : new_elem = std::make_unique<Tri3>();
1650 180 : if (tri_elem_type == TRI7)
1651 : {
1652 90 : new_elem = std::make_unique<Tri7>();
1653 90 : new_elem->set_node(3,
1654 90 : mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * 2 * orig_nodes)));
1655 90 : new_elem->set_node(4,
1656 90 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 2)->id() +
1657 90 : ((current_layer * 2 + 1) * orig_nodes)));
1658 90 : new_elem->set_node(
1659 : 5,
1660 90 : mesh->node_ptr(elem->node_ptr(2)->id() +
1661 90 : ((current_layer + 1) %
1662 90 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1663 90 : 2 * orig_nodes)));
1664 90 : new_elem->set_node(
1665 90 : 6, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1666 : }
1667 :
1668 180 : new_elem->set_node(0, mesh->node_ptr(elem->node_ptr(axis_node_case)->id()));
1669 180 : new_elem->set_node(1,
1670 180 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 2)->id() +
1671 180 : (current_layer * order * orig_nodes)));
1672 180 : new_elem->set_node(
1673 : 2,
1674 180 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 2)->id() +
1675 180 : ((current_layer + 1) %
1676 180 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1677 180 : order * orig_nodes)));
1678 :
1679 180 : if (new_elem->volume() < 0.0)
1680 : {
1681 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 2);
1682 0 : if (tri_elem_type == TRI7)
1683 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 5);
1684 0 : is_flipped = true;
1685 : }
1686 180 : }
1687 :
1688 : void
1689 3702 : RevolveGenerator::createPRISMfromTRI(const ElemType prism_elem_type,
1690 : const Elem * elem,
1691 : const std::unique_ptr<MeshBase> & mesh,
1692 : std::unique_ptr<Elem> & new_elem,
1693 : const int current_layer,
1694 : const unsigned int orig_nodes,
1695 : const unsigned int total_num_azimuthal_intervals,
1696 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
1697 : bool & is_flipped) const
1698 : {
1699 3702 : if (prism_elem_type != PRISM6 && prism_elem_type != PRISM18 && prism_elem_type != PRISM21)
1700 0 : mooseError("unsupported situation");
1701 :
1702 3702 : side_pairs.push_back(std::make_pair(0, 4));
1703 3702 : const unsigned int order = prism_elem_type == PRISM6 ? 1 : 2;
1704 :
1705 3702 : new_elem = std::make_unique<Prism6>();
1706 :
1707 3702 : if (order == 2)
1708 : {
1709 408 : new_elem = std::make_unique<Prism18>();
1710 408 : if (prism_elem_type == PRISM21)
1711 : {
1712 408 : new_elem = std::make_unique<Prism21>();
1713 204 : new_elem->set_node(
1714 204 : 18, mesh->node_ptr(elem->node_ptr(6)->id() + (current_layer * 2 * orig_nodes)));
1715 204 : new_elem->set_node(
1716 : 19,
1717 204 : mesh->node_ptr(elem->node_ptr(6)->id() + ((current_layer + 1) %
1718 204 : (total_num_azimuthal_intervals + 1 -
1719 204 : (unsigned int)_full_circle_revolving) *
1720 204 : 2 * orig_nodes)));
1721 204 : new_elem->set_node(
1722 204 : 20, mesh->node_ptr(elem->node_ptr(6)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1723 : }
1724 408 : new_elem->set_node(6,
1725 408 : mesh->node_ptr(elem->node_ptr(3)->id() + (current_layer * 2 * orig_nodes)));
1726 408 : new_elem->set_node(7,
1727 408 : mesh->node_ptr(elem->node_ptr(4)->id() + (current_layer * 2 * orig_nodes)));
1728 408 : new_elem->set_node(8,
1729 408 : mesh->node_ptr(elem->node_ptr(5)->id() + (current_layer * 2 * orig_nodes)));
1730 408 : new_elem->set_node(
1731 408 : 9, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1732 408 : new_elem->set_node(
1733 408 : 10, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1734 408 : new_elem->set_node(
1735 408 : 11, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1736 408 : new_elem->set_node(
1737 : 12,
1738 408 : mesh->node_ptr(elem->node_ptr(3)->id() +
1739 408 : ((current_layer + 1) %
1740 408 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1741 408 : 2 * orig_nodes)));
1742 408 : new_elem->set_node(
1743 : 13,
1744 408 : mesh->node_ptr(elem->node_ptr(4)->id() +
1745 408 : ((current_layer + 1) %
1746 408 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1747 408 : 2 * orig_nodes)));
1748 408 : new_elem->set_node(
1749 : 14,
1750 408 : mesh->node_ptr(elem->node_ptr(5)->id() +
1751 408 : ((current_layer + 1) %
1752 408 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1753 408 : 2 * orig_nodes)));
1754 408 : new_elem->set_node(
1755 408 : 15, mesh->node_ptr(elem->node_ptr(3)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1756 408 : new_elem->set_node(
1757 408 : 16, mesh->node_ptr(elem->node_ptr(4)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1758 408 : new_elem->set_node(
1759 408 : 17, mesh->node_ptr(elem->node_ptr(5)->id() + ((current_layer * 2 + 1) * orig_nodes)));
1760 : }
1761 3702 : new_elem->set_node(
1762 3702 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * order * orig_nodes)));
1763 3702 : new_elem->set_node(
1764 3702 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * order * orig_nodes)));
1765 3702 : new_elem->set_node(
1766 3702 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * order * orig_nodes)));
1767 3702 : new_elem->set_node(
1768 : 3,
1769 3702 : mesh->node_ptr(elem->node_ptr(0)->id() +
1770 3702 : ((current_layer + 1) %
1771 3702 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1772 3702 : order * orig_nodes)));
1773 3702 : new_elem->set_node(
1774 : 4,
1775 3702 : mesh->node_ptr(elem->node_ptr(1)->id() +
1776 3702 : ((current_layer + 1) %
1777 3702 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1778 3702 : order * orig_nodes)));
1779 3702 : new_elem->set_node(
1780 : 5,
1781 3702 : mesh->node_ptr(elem->node_ptr(2)->id() +
1782 3702 : ((current_layer + 1) %
1783 3702 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1784 3702 : order * orig_nodes)));
1785 :
1786 3702 : if (new_elem->volume() < 0.0)
1787 : {
1788 1110 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
1789 1110 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 4);
1790 1110 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 5);
1791 1110 : if (prism_elem_type != PRISM6)
1792 : {
1793 408 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 12);
1794 408 : MooseMeshUtils::swapNodesInElem(*new_elem, 7, 13);
1795 408 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 14);
1796 408 : if (prism_elem_type == PRISM21)
1797 : {
1798 204 : MooseMeshUtils::swapNodesInElem(*new_elem, 18, 19);
1799 : }
1800 : }
1801 1110 : is_flipped = true;
1802 : }
1803 3702 : }
1804 :
1805 : void
1806 432 : RevolveGenerator::createPYRAMIDfromTRI(
1807 : const std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>> & nodes_cates,
1808 : const ElemType pyramid_elem_type,
1809 : const Elem * elem,
1810 : const std::unique_ptr<MeshBase> & mesh,
1811 : std::unique_ptr<Elem> & new_elem,
1812 : const int current_layer,
1813 : const unsigned int orig_nodes,
1814 : const unsigned int total_num_azimuthal_intervals,
1815 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
1816 : dof_id_type & axis_node_case,
1817 : bool & is_flipped) const
1818 : {
1819 432 : if (pyramid_elem_type != PYRAMID5 && pyramid_elem_type != PYRAMID13 &&
1820 432 : pyramid_elem_type != PYRAMID18)
1821 0 : mooseError("unsupported situation");
1822 :
1823 432 : side_pairs.push_back(std::make_pair(0, 2));
1824 432 : const unsigned int order = pyramid_elem_type == PYRAMID5 ? 1 : 2;
1825 432 : axis_node_case = nodes_cates.first.front();
1826 :
1827 432 : new_elem = std::make_unique<Pyramid5>();
1828 :
1829 432 : if (order == 2)
1830 : {
1831 216 : new_elem = std::make_unique<Pyramid13>();
1832 216 : if (pyramid_elem_type == PYRAMID18)
1833 : {
1834 216 : new_elem = std::make_unique<Pyramid18>();
1835 108 : new_elem->set_node(13,
1836 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3 + 3)->id() +
1837 108 : ((current_layer * 2 + 1) * orig_nodes)));
1838 108 : new_elem->set_node(15,
1839 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3 + 3)->id() +
1840 : ((current_layer * 2 + 1) * orig_nodes)));
1841 108 : new_elem->set_node(17,
1842 108 : mesh->node_ptr(elem->node_ptr(axis_node_case + 3)->id() +
1843 : ((current_layer * 2 + 1) * orig_nodes)));
1844 108 : new_elem->set_node(
1845 108 : 14, mesh->node_ptr(elem->node_ptr(6)->id() + (current_layer * 2 * orig_nodes)));
1846 108 : new_elem->set_node(
1847 : 16,
1848 108 : mesh->node_ptr(elem->node_ptr(6)->id() + ((current_layer + 1) %
1849 108 : (total_num_azimuthal_intervals + 1 -
1850 108 : (unsigned int)_full_circle_revolving) *
1851 108 : 2 * orig_nodes)));
1852 : }
1853 216 : new_elem->set_node(6,
1854 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3)->id() +
1855 216 : ((current_layer * 2 + 1) * orig_nodes)));
1856 216 : new_elem->set_node(8,
1857 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3)->id() +
1858 : ((current_layer * 2 + 1) * orig_nodes)));
1859 216 : new_elem->set_node(5,
1860 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3 + 3)->id() +
1861 216 : (current_layer * 2 * orig_nodes)));
1862 216 : new_elem->set_node(
1863 : 7,
1864 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3 + 3)->id() +
1865 216 : ((current_layer + 1) %
1866 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1867 216 : 2 * orig_nodes)));
1868 216 : new_elem->set_node(9,
1869 216 : mesh->node_ptr(elem->node_ptr(axis_node_case + 3)->id() +
1870 : (current_layer * 2 * orig_nodes)));
1871 216 : new_elem->set_node(10,
1872 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3 + 3)->id() +
1873 : (current_layer * 2 * orig_nodes)));
1874 216 : new_elem->set_node(
1875 : 12,
1876 216 : mesh->node_ptr(elem->node_ptr(axis_node_case + 3)->id() +
1877 216 : ((current_layer + 1) %
1878 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1879 216 : 2 * orig_nodes)));
1880 216 : new_elem->set_node(
1881 : 11,
1882 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3 + 3)->id() +
1883 216 : ((current_layer + 1) %
1884 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1885 216 : 2 * orig_nodes)));
1886 : }
1887 432 : new_elem->set_node(4, mesh->node_ptr(elem->node_ptr(axis_node_case)->id()));
1888 432 : new_elem->set_node(0,
1889 432 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3)->id() +
1890 432 : (current_layer * order * orig_nodes)));
1891 432 : new_elem->set_node(1,
1892 432 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3)->id() +
1893 : (current_layer * order * orig_nodes)));
1894 432 : new_elem->set_node(
1895 : 2,
1896 432 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 3)->id() +
1897 432 : ((current_layer + 1) %
1898 432 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1899 432 : order * orig_nodes)));
1900 432 : new_elem->set_node(
1901 : 3,
1902 432 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 3)->id() +
1903 432 : ((current_layer + 1) %
1904 432 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
1905 432 : order * orig_nodes)));
1906 :
1907 432 : if (new_elem->volume() < 0.0)
1908 : {
1909 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 2);
1910 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
1911 0 : if (order == 2)
1912 : {
1913 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 5, 7);
1914 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 11);
1915 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 12);
1916 0 : if (pyramid_elem_type == PYRAMID18)
1917 : {
1918 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 14, 16);
1919 : }
1920 : }
1921 0 : is_flipped = true;
1922 : }
1923 432 : }
1924 :
1925 : void
1926 246 : RevolveGenerator::createTETfromTRI(
1927 : const std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>> & nodes_cates,
1928 : const ElemType tet_elem_type,
1929 : const Elem * elem,
1930 : const std::unique_ptr<MeshBase> & mesh,
1931 : std::unique_ptr<Elem> & new_elem,
1932 : const int current_layer,
1933 : const unsigned int orig_nodes,
1934 : const unsigned int total_num_azimuthal_intervals,
1935 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
1936 : dof_id_type & axis_node_case,
1937 : bool & is_flipped) const
1938 : {
1939 246 : if (tet_elem_type != TET4 && tet_elem_type != TET10 && tet_elem_type != TET14)
1940 0 : mooseError("unsupported situation");
1941 :
1942 246 : side_pairs.push_back(std::make_pair(0, 1));
1943 246 : const unsigned int order = tet_elem_type == TET4 ? 1 : 2;
1944 : if (order == 2)
1945 : {
1946 : // Sanity check to filter unsupported cases
1947 192 : if (nodes_cates.first[0] > 2 || nodes_cates.first[1] > 2 || nodes_cates.first[2] < 3)
1948 0 : mooseError("unsupported situation 2");
1949 : }
1950 246 : axis_node_case = nodes_cates.second.front();
1951 :
1952 246 : new_elem = std::make_unique<Tet4>();
1953 246 : if (order == 2)
1954 : {
1955 192 : const bool node_order = nodes_cates.first[1] - nodes_cates.first[0] == 1;
1956 192 : new_elem = std::make_unique<Tet10>();
1957 192 : if (tet_elem_type == TET14)
1958 : {
1959 192 : new_elem = std::make_unique<Tet14>();
1960 96 : new_elem->set_node(
1961 : 12,
1962 96 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.first[1] + 3)
1963 96 : : (nodes_cates.second.front() + 3))
1964 : ->id() +
1965 96 : ((current_layer * 2 + 1) * orig_nodes)));
1966 96 : new_elem->set_node(13,
1967 96 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.second.front() + 3)
1968 96 : : (nodes_cates.first[0] + 3))
1969 : ->id() +
1970 : ((current_layer * 2 + 1) * orig_nodes)));
1971 96 : new_elem->set_node(
1972 96 : 10, mesh->node_ptr(elem->node_ptr(6)->id() + (current_layer * 2 * orig_nodes)));
1973 96 : new_elem->set_node(
1974 : 11,
1975 96 : mesh->node_ptr(elem->node_ptr(6)->id() + ((current_layer + 1) %
1976 96 : (total_num_azimuthal_intervals + 1 -
1977 96 : (unsigned int)_full_circle_revolving) *
1978 96 : 2 * orig_nodes)));
1979 : }
1980 192 : new_elem->set_node(4,
1981 192 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.first[0] + 3)
1982 192 : : (nodes_cates.first[1] + 3))
1983 : ->id()));
1984 192 : new_elem->set_node(5,
1985 192 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.first[1] + 3)
1986 192 : : (nodes_cates.second.front() + 3))
1987 : ->id() +
1988 192 : (current_layer * 2 * orig_nodes)));
1989 192 : new_elem->set_node(6,
1990 192 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.second.front() + 3)
1991 192 : : (nodes_cates.first[0] + 3))
1992 : ->id() +
1993 : (current_layer * 2 * orig_nodes)));
1994 192 : new_elem->set_node(
1995 : 8,
1996 192 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.first[1] + 3)
1997 192 : : (nodes_cates.second.front() + 3))
1998 : ->id() +
1999 192 : ((current_layer + 1) %
2000 192 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2001 192 : 2 * orig_nodes)));
2002 192 : new_elem->set_node(
2003 : 7,
2004 192 : mesh->node_ptr(elem->node_ptr(node_order ? (nodes_cates.second.front() + 3)
2005 192 : : (nodes_cates.first[0] + 3))
2006 : ->id() +
2007 192 : ((current_layer + 1) %
2008 192 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2009 192 : 2 * orig_nodes)));
2010 192 : new_elem->set_node(9,
2011 192 : mesh->node_ptr(elem->node_ptr(nodes_cates.second.front())->id() +
2012 192 : ((current_layer * 2 + 1) * orig_nodes)));
2013 : }
2014 246 : new_elem->set_node(0, mesh->node_ptr(elem->node_ptr(nodes_cates.first[0])->id()));
2015 246 : new_elem->set_node(1, mesh->node_ptr(elem->node_ptr(nodes_cates.first[1])->id()));
2016 246 : new_elem->set_node(2,
2017 246 : mesh->node_ptr(elem->node_ptr(nodes_cates.second.front())->id() +
2018 246 : (current_layer * order * orig_nodes)));
2019 246 : new_elem->set_node(
2020 : 3,
2021 246 : mesh->node_ptr(elem->node_ptr(nodes_cates.second.front())->id() +
2022 246 : ((current_layer + 1) %
2023 246 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2024 246 : order * orig_nodes)));
2025 :
2026 246 : if (new_elem->volume() < 0.0)
2027 : {
2028 54 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 3);
2029 54 : if (order == 2)
2030 : {
2031 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 5, 8);
2032 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 7);
2033 0 : if (tet_elem_type == TET14)
2034 : {
2035 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 11);
2036 : }
2037 : }
2038 54 : is_flipped = true;
2039 : }
2040 246 : }
2041 :
2042 : void
2043 7506 : RevolveGenerator::createHEXfromQUAD(const ElemType hex_elem_type,
2044 : const Elem * elem,
2045 : const std::unique_ptr<MeshBase> & mesh,
2046 : std::unique_ptr<Elem> & new_elem,
2047 : const int current_layer,
2048 : const unsigned int orig_nodes,
2049 : const unsigned int total_num_azimuthal_intervals,
2050 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
2051 : bool & is_flipped) const
2052 : {
2053 7506 : if (hex_elem_type != HEX8 && hex_elem_type != HEX20 && hex_elem_type != HEX27)
2054 0 : mooseError("unsupported situation");
2055 :
2056 7506 : side_pairs.push_back(std::make_pair(0, 5));
2057 7506 : const unsigned int order = hex_elem_type == HEX8 ? 1 : 2;
2058 :
2059 7506 : new_elem = std::make_unique<Hex8>();
2060 7506 : if (order == 2)
2061 : {
2062 540 : new_elem = std::make_unique<Hex20>();
2063 540 : if (hex_elem_type == HEX27)
2064 : {
2065 540 : new_elem = std::make_unique<Hex27>();
2066 270 : new_elem->set_node(
2067 270 : 20, mesh->node_ptr(elem->node_ptr(8)->id() + (current_layer * 2 * orig_nodes)));
2068 270 : new_elem->set_node(
2069 : 25,
2070 270 : mesh->node_ptr(elem->node_ptr(8)->id() + ((current_layer + 1) %
2071 270 : (total_num_azimuthal_intervals + 1 -
2072 270 : (unsigned int)_full_circle_revolving) *
2073 270 : 2 * orig_nodes)));
2074 270 : new_elem->set_node(
2075 270 : 26, mesh->node_ptr(elem->node_ptr(8)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2076 270 : new_elem->set_node(
2077 270 : 21, mesh->node_ptr(elem->node_ptr(4)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2078 270 : new_elem->set_node(
2079 270 : 22, mesh->node_ptr(elem->node_ptr(5)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2080 270 : new_elem->set_node(
2081 270 : 23, mesh->node_ptr(elem->node_ptr(6)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2082 270 : new_elem->set_node(
2083 270 : 24, mesh->node_ptr(elem->node_ptr(7)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2084 : }
2085 540 : new_elem->set_node(8,
2086 540 : mesh->node_ptr(elem->node_ptr(4)->id() + (current_layer * 2 * orig_nodes)));
2087 540 : new_elem->set_node(9,
2088 540 : mesh->node_ptr(elem->node_ptr(5)->id() + (current_layer * 2 * orig_nodes)));
2089 540 : new_elem->set_node(10,
2090 540 : mesh->node_ptr(elem->node_ptr(6)->id() + (current_layer * 2 * orig_nodes)));
2091 540 : new_elem->set_node(11,
2092 540 : mesh->node_ptr(elem->node_ptr(7)->id() + (current_layer * 2 * orig_nodes)));
2093 540 : new_elem->set_node(
2094 : 16,
2095 540 : mesh->node_ptr(elem->node_ptr(4)->id() +
2096 540 : ((current_layer + 1) %
2097 540 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2098 540 : 2 * orig_nodes)));
2099 540 : new_elem->set_node(
2100 : 17,
2101 540 : mesh->node_ptr(elem->node_ptr(5)->id() +
2102 540 : ((current_layer + 1) %
2103 540 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2104 540 : 2 * orig_nodes)));
2105 540 : new_elem->set_node(
2106 : 18,
2107 540 : mesh->node_ptr(elem->node_ptr(6)->id() +
2108 540 : ((current_layer + 1) %
2109 540 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2110 540 : 2 * orig_nodes)));
2111 540 : new_elem->set_node(
2112 : 19,
2113 540 : mesh->node_ptr(elem->node_ptr(7)->id() +
2114 540 : ((current_layer + 1) %
2115 540 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2116 540 : 2 * orig_nodes)));
2117 540 : new_elem->set_node(
2118 540 : 12, mesh->node_ptr(elem->node_ptr(0)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2119 540 : new_elem->set_node(
2120 540 : 13, mesh->node_ptr(elem->node_ptr(1)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2121 540 : new_elem->set_node(
2122 540 : 14, mesh->node_ptr(elem->node_ptr(2)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2123 540 : new_elem->set_node(
2124 540 : 15, mesh->node_ptr(elem->node_ptr(3)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2125 : }
2126 7506 : new_elem->set_node(
2127 7506 : 0, mesh->node_ptr(elem->node_ptr(0)->id() + (current_layer * order * orig_nodes)));
2128 7506 : new_elem->set_node(
2129 7506 : 1, mesh->node_ptr(elem->node_ptr(1)->id() + (current_layer * order * orig_nodes)));
2130 7506 : new_elem->set_node(
2131 7506 : 2, mesh->node_ptr(elem->node_ptr(2)->id() + (current_layer * order * orig_nodes)));
2132 7506 : new_elem->set_node(
2133 7506 : 3, mesh->node_ptr(elem->node_ptr(3)->id() + (current_layer * order * orig_nodes)));
2134 7506 : new_elem->set_node(
2135 : 4,
2136 7506 : mesh->node_ptr(elem->node_ptr(0)->id() +
2137 7506 : ((current_layer + 1) %
2138 7506 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2139 7506 : order * orig_nodes)));
2140 7506 : new_elem->set_node(
2141 : 5,
2142 7506 : mesh->node_ptr(elem->node_ptr(1)->id() +
2143 7506 : ((current_layer + 1) %
2144 7506 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2145 7506 : order * orig_nodes)));
2146 7506 : new_elem->set_node(
2147 : 6,
2148 7506 : mesh->node_ptr(elem->node_ptr(2)->id() +
2149 7506 : ((current_layer + 1) %
2150 7506 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2151 7506 : order * orig_nodes)));
2152 7506 : new_elem->set_node(
2153 : 7,
2154 7506 : mesh->node_ptr(elem->node_ptr(3)->id() +
2155 7506 : ((current_layer + 1) %
2156 7506 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2157 7506 : order * orig_nodes)));
2158 :
2159 7506 : if (new_elem->volume() < 0.0)
2160 : {
2161 2214 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 4);
2162 2214 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 5);
2163 2214 : MooseMeshUtils::swapNodesInElem(*new_elem, 2, 6);
2164 2214 : MooseMeshUtils::swapNodesInElem(*new_elem, 3, 7);
2165 2214 : if (order == 2)
2166 : {
2167 540 : MooseMeshUtils::swapNodesInElem(*new_elem, 8, 16);
2168 540 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 17);
2169 540 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 18);
2170 540 : MooseMeshUtils::swapNodesInElem(*new_elem, 11, 19);
2171 540 : if (hex_elem_type == HEX27)
2172 : {
2173 270 : MooseMeshUtils::swapNodesInElem(*new_elem, 20, 25);
2174 : }
2175 : }
2176 2214 : is_flipped = true;
2177 : }
2178 7506 : }
2179 :
2180 : void
2181 270 : RevolveGenerator::createPRISMfromQUAD(
2182 : const std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>> & nodes_cates,
2183 : const ElemType prism_elem_type,
2184 : const Elem * elem,
2185 : const std::unique_ptr<MeshBase> & mesh,
2186 : std::unique_ptr<Elem> & new_elem,
2187 : const int current_layer,
2188 : const unsigned int orig_nodes,
2189 : const unsigned int total_num_azimuthal_intervals,
2190 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
2191 : dof_id_type & axis_node_case,
2192 : bool & is_flipped) const
2193 : {
2194 270 : if (prism_elem_type != PRISM6 && prism_elem_type != PRISM15 && prism_elem_type != PRISM18)
2195 0 : mooseError("unsupported situation");
2196 :
2197 270 : side_pairs.push_back(std::make_pair(1, 3));
2198 270 : const unsigned int order = prism_elem_type == PRISM6 ? 1 : 2;
2199 : if (order == 2)
2200 : {
2201 : // Sanity check to filter unsupported cases
2202 216 : if (nodes_cates.first[0] > 3 || nodes_cates.first[1] > 3 || nodes_cates.first[2] < 4)
2203 0 : mooseError("unsupported situation 2");
2204 : }
2205 : // Can only be 0-1, 1-2, 2-3, 3-0, we only consider vetices here.
2206 : // nodes_cates are natually sorted
2207 270 : const dof_id_type min_on_axis = nodes_cates.first[0];
2208 270 : const dof_id_type max_on_axis = nodes_cates.first[1];
2209 270 : axis_node_case = max_on_axis - min_on_axis == 1 ? min_on_axis : max_on_axis;
2210 :
2211 270 : new_elem = std::make_unique<Prism6>();
2212 270 : if (order == 2)
2213 : {
2214 216 : new_elem = std::make_unique<Prism15>();
2215 216 : if (prism_elem_type == PRISM18)
2216 : {
2217 216 : new_elem = std::make_unique<Prism18>();
2218 108 : new_elem->set_node(
2219 108 : 15, mesh->node_ptr(elem->node_ptr(8)->id() + (current_layer * 2 * orig_nodes)));
2220 108 : new_elem->set_node(16,
2221 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2222 108 : ((current_layer * 2 + 1) * orig_nodes)));
2223 108 : new_elem->set_node(
2224 : 17,
2225 108 : mesh->node_ptr(elem->node_ptr(8)->id() + ((current_layer + 1) %
2226 108 : (total_num_azimuthal_intervals + 1 -
2227 108 : (unsigned int)_full_circle_revolving) *
2228 108 : 2 * orig_nodes)));
2229 : }
2230 216 : new_elem->set_node(9, mesh->node_ptr(elem->node_ptr(axis_node_case + 4)->id()));
2231 216 : new_elem->set_node(10,
2232 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2233 216 : (current_layer * 2 * orig_nodes)));
2234 216 : new_elem->set_node(12,
2235 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4 + 4)->id() +
2236 : (current_layer * 2 * orig_nodes)));
2237 216 : new_elem->set_node(6,
2238 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4 + 4)->id() +
2239 : (current_layer * 2 * orig_nodes)));
2240 216 : new_elem->set_node(
2241 : 14,
2242 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4 + 4)->id() +
2243 216 : ((current_layer + 1) %
2244 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2245 216 : 2 * orig_nodes)));
2246 216 : new_elem->set_node(
2247 : 8,
2248 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4 + 4)->id() +
2249 216 : ((current_layer + 1) %
2250 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2251 216 : 2 * orig_nodes)));
2252 216 : new_elem->set_node(
2253 : 11,
2254 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2255 216 : ((current_layer + 1) %
2256 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2257 216 : 2 * orig_nodes)));
2258 216 : new_elem->set_node(7,
2259 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2260 216 : ((current_layer * 2 + 1) * orig_nodes)));
2261 216 : new_elem->set_node(13,
2262 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2263 : ((current_layer * 2 + 1) * orig_nodes)));
2264 : }
2265 270 : new_elem->set_node(0, mesh->node_ptr(elem->node_ptr(axis_node_case)->id()));
2266 270 : new_elem->set_node(3, mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id()));
2267 270 : new_elem->set_node(4,
2268 270 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2269 270 : (current_layer * order * orig_nodes)));
2270 270 : new_elem->set_node(1,
2271 270 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2272 : (current_layer * order * orig_nodes)));
2273 270 : new_elem->set_node(
2274 : 5,
2275 270 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2276 270 : ((current_layer + 1) %
2277 270 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2278 270 : order * orig_nodes)));
2279 270 : new_elem->set_node(
2280 : 2,
2281 270 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2282 270 : ((current_layer + 1) %
2283 270 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2284 270 : order * orig_nodes)));
2285 :
2286 270 : if (new_elem->volume() < 0.0)
2287 : {
2288 270 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 2);
2289 270 : MooseMeshUtils::swapNodesInElem(*new_elem, 4, 5);
2290 270 : if (order == 2)
2291 : {
2292 216 : MooseMeshUtils::swapNodesInElem(*new_elem, 12, 14);
2293 216 : MooseMeshUtils::swapNodesInElem(*new_elem, 6, 8);
2294 216 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 11);
2295 216 : if (prism_elem_type == PRISM18)
2296 : {
2297 108 : MooseMeshUtils::swapNodesInElem(*new_elem, 15, 17);
2298 : }
2299 : }
2300 270 : is_flipped = true;
2301 : }
2302 270 : }
2303 :
2304 : void
2305 216 : RevolveGenerator::createPYRAMIDPRISMfromQUAD(
2306 : const std::pair<std::vector<dof_id_type>, std::vector<dof_id_type>> & nodes_cates,
2307 : const ElemType pyramid_elem_type,
2308 : const ElemType prism_elem_type,
2309 : const Elem * elem,
2310 : const std::unique_ptr<MeshBase> & mesh,
2311 : std::unique_ptr<Elem> & new_elem,
2312 : std::unique_ptr<Elem> & new_elem_1,
2313 : const int current_layer,
2314 : const unsigned int orig_nodes,
2315 : const unsigned int total_num_azimuthal_intervals,
2316 : std::vector<std::pair<dof_id_type, dof_id_type>> & side_pairs,
2317 : dof_id_type & axis_node_case,
2318 : bool & is_flipped,
2319 : bool & is_flipped_additional) const
2320 : {
2321 216 : if (!(pyramid_elem_type == PYRAMID5 && prism_elem_type == PRISM6) &&
2322 108 : !(pyramid_elem_type == PYRAMID14 && prism_elem_type == PRISM18))
2323 0 : mooseError("unsupported situation");
2324 : const unsigned int order = pyramid_elem_type == PYRAMID5 ? 1 : 2;
2325 :
2326 216 : side_pairs.push_back(std::make_pair(0, 2));
2327 216 : axis_node_case = nodes_cates.first.front();
2328 216 : new_elem = std::make_unique<Pyramid5>();
2329 216 : if (pyramid_elem_type == PYRAMID14)
2330 : {
2331 216 : new_elem = std::make_unique<Pyramid14>();
2332 108 : new_elem->set_node(9,
2333 108 : mesh->node_ptr(elem->node_ptr(axis_node_case + 4)->id() +
2334 108 : (current_layer * 2 * orig_nodes)));
2335 108 : new_elem->set_node(10,
2336 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4 + 4)->id() +
2337 : (current_layer * 2 * orig_nodes)));
2338 108 : new_elem->set_node(5,
2339 108 : mesh->node_ptr(elem->node_ptr(8)->id() + (current_layer * 2 * orig_nodes)));
2340 108 : new_elem->set_node(8,
2341 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2342 108 : ((current_layer * 2 + 1) * orig_nodes)));
2343 108 : new_elem->set_node(6,
2344 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2345 : ((current_layer * 2 + 1) * orig_nodes)));
2346 108 : new_elem->set_node(
2347 108 : 13, mesh->node_ptr(elem->node_ptr(8)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2348 108 : new_elem->set_node(
2349 : 7,
2350 108 : mesh->node_ptr(elem->node_ptr(8)->id() +
2351 108 : ((current_layer + 1) %
2352 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2353 108 : 2 * orig_nodes)));
2354 108 : new_elem->set_node(
2355 : 11,
2356 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4 + 4)->id() +
2357 108 : ((current_layer + 1) %
2358 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2359 108 : 2 * orig_nodes)));
2360 108 : new_elem->set_node(
2361 : 12,
2362 108 : mesh->node_ptr(elem->node_ptr(axis_node_case + 4)->id() +
2363 108 : ((current_layer + 1) %
2364 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2365 108 : 2 * orig_nodes)));
2366 : }
2367 216 : new_elem->set_node(4, mesh->node_ptr(elem->node_ptr(axis_node_case)->id()));
2368 216 : new_elem->set_node(0,
2369 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2370 216 : (current_layer * order * orig_nodes)));
2371 216 : new_elem->set_node(1,
2372 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2373 : (current_layer * order * orig_nodes)));
2374 216 : new_elem->set_node(
2375 : 2,
2376 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2377 216 : ((current_layer + 1) %
2378 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2379 216 : order * orig_nodes)));
2380 216 : new_elem->set_node(
2381 : 3,
2382 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2383 216 : ((current_layer + 1) %
2384 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2385 216 : order * orig_nodes)));
2386 :
2387 216 : if (new_elem->volume() < 0.0)
2388 : {
2389 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 1, 2);
2390 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 0, 3);
2391 0 : if (pyramid_elem_type == PYRAMID14)
2392 : {
2393 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 5, 7);
2394 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 10, 11);
2395 0 : MooseMeshUtils::swapNodesInElem(*new_elem, 9, 12);
2396 : }
2397 0 : is_flipped = true;
2398 : }
2399 :
2400 216 : side_pairs.push_back(std::make_pair(0, 4));
2401 216 : new_elem_1 = std::make_unique<Prism6>();
2402 216 : if (prism_elem_type == PRISM18)
2403 : {
2404 216 : new_elem_1 = std::make_unique<Prism18>();
2405 108 : new_elem_1->set_node(
2406 108 : 6, mesh->node_ptr(elem->node_ptr(8)->id() + (current_layer * 2 * orig_nodes)));
2407 108 : new_elem_1->set_node(8,
2408 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4 + 4)->id() +
2409 : (current_layer * 2 * orig_nodes)));
2410 108 : new_elem_1->set_node(7,
2411 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2412 : (current_layer * 2 * orig_nodes)));
2413 108 : new_elem_1->set_node(9,
2414 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2415 108 : ((current_layer * 2 + 1) * orig_nodes)));
2416 108 : new_elem_1->set_node(10,
2417 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2418 : ((current_layer * 2 + 1) * orig_nodes)));
2419 108 : new_elem_1->set_node(11,
2420 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2421 : ((current_layer * 2 + 1) * orig_nodes)));
2422 108 : new_elem_1->set_node(
2423 108 : 15, mesh->node_ptr(elem->node_ptr(8)->id() + ((current_layer * 2 + 1) * orig_nodes)));
2424 108 : new_elem_1->set_node(17,
2425 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4 + 4)->id() +
2426 : ((current_layer * 2 + 1) * orig_nodes)));
2427 108 : new_elem_1->set_node(16,
2428 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2429 : ((current_layer * 2 + 1) * orig_nodes)));
2430 108 : new_elem_1->set_node(
2431 : 12,
2432 108 : mesh->node_ptr(elem->node_ptr(8)->id() +
2433 108 : ((current_layer + 1) %
2434 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2435 108 : 2 * orig_nodes)));
2436 108 : new_elem_1->set_node(
2437 : 13,
2438 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4 + 4)->id() +
2439 108 : ((current_layer + 1) %
2440 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2441 108 : 2 * orig_nodes)));
2442 108 : new_elem_1->set_node(
2443 : 14,
2444 108 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4 + 4)->id() +
2445 108 : ((current_layer + 1) %
2446 108 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2447 108 : 2 * orig_nodes)));
2448 : }
2449 216 : new_elem_1->set_node(0,
2450 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2451 : (current_layer * order * orig_nodes)));
2452 216 : new_elem_1->set_node(1,
2453 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2454 : (current_layer * order * orig_nodes)));
2455 216 : new_elem_1->set_node(2,
2456 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2457 : (current_layer * order * orig_nodes)));
2458 216 : new_elem_1->set_node(
2459 : 3,
2460 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 1) % 4)->id() +
2461 216 : ((current_layer + 1) %
2462 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2463 216 : order * orig_nodes)));
2464 216 : new_elem_1->set_node(
2465 : 4,
2466 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 3) % 4)->id() +
2467 216 : ((current_layer + 1) %
2468 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2469 216 : order * orig_nodes)));
2470 216 : new_elem_1->set_node(
2471 : 5,
2472 216 : mesh->node_ptr(elem->node_ptr((axis_node_case + 2) % 4)->id() +
2473 216 : ((current_layer + 1) %
2474 216 : (total_num_azimuthal_intervals + 1 - (unsigned int)_full_circle_revolving) *
2475 216 : order * orig_nodes)));
2476 :
2477 216 : if (new_elem_1->volume() < 0.0)
2478 : {
2479 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 0, 3);
2480 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 1, 4);
2481 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 2, 5);
2482 0 : if (prism_elem_type == PRISM18)
2483 : {
2484 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 6, 12);
2485 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 7, 13);
2486 0 : MooseMeshUtils::swapNodesInElem(*new_elem_1, 8, 14);
2487 : }
2488 0 : is_flipped_additional = true;
2489 : }
2490 216 : }
|
