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