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 "PolygonMeshGeneratorBase.h"
11 : #include "MooseUtils.h"
12 : #include "FormattedTable.h"
13 :
14 : #include <cmath>
15 : #include <iomanip>
16 :
17 : using namespace libMesh;
18 :
19 : InputParameters
20 10370 : PolygonMeshGeneratorBase::validParams()
21 : {
22 10370 : InputParameters params = MeshGenerator::validParams();
23 10370 : params.addClassDescription(
24 : "A base class that contains common members for Reactor module mesh generators.");
25 :
26 10370 : return params;
27 0 : }
28 :
29 5295 : PolygonMeshGeneratorBase::PolygonMeshGeneratorBase(const InputParameters & parameters)
30 5295 : : MeshGenerator(parameters)
31 : {
32 5295 : }
33 :
34 : std::unique_ptr<MeshBase>
35 0 : PolygonMeshGeneratorBase::generate()
36 : {
37 0 : auto mesh = buildReplicatedMesh(2); // initiate a 2D mesh
38 0 : return dynamic_pointer_cast<MeshBase>(mesh);
39 0 : }
40 :
41 : std::unique_ptr<ReplicatedMesh>
42 210 : PolygonMeshGeneratorBase::buildGeneralSlice(
43 : std::vector<Real> ring_radii,
44 : const std::vector<unsigned int> ring_layers,
45 : const std::vector<Real> ring_radial_biases,
46 : const multiBdryLayerParams & ring_inner_boundary_layer_params,
47 : const multiBdryLayerParams & ring_outer_boundary_layer_params,
48 : std::vector<Real> ducts_center_dist,
49 : const std::vector<unsigned int> ducts_layers,
50 : const std::vector<Real> duct_radial_biases,
51 : const multiBdryLayerParams & duct_inner_boundary_layer_params,
52 : const multiBdryLayerParams & duct_outer_boundary_layer_params,
53 : const Real primary_side_length,
54 : const Real secondary_side_length,
55 : const unsigned int num_sectors_per_side,
56 : const unsigned int background_intervals,
57 : const Real background_radial_bias,
58 : const singleBdryLayerParams & background_inner_boundary_layer_params,
59 : const singleBdryLayerParams & background_outer_boundary_layer_params,
60 : dof_id_type & node_id_background_meta,
61 : const Real azimuthal_angle,
62 : const std::vector<Real> azimuthal_tangent,
63 : const unsigned int side_index,
64 : const bool quad_center_elements,
65 : const Real center_quad_factor,
66 : const Real rotation_angle,
67 : const bool generate_side_specific_boundaries)
68 : {
69 210 : const Real virtual_pitch = 2.0 * primary_side_length * cos(azimuthal_angle / 360.0 * M_PI);
70 210 : const Real virtual_side_number = 360.0 / azimuthal_angle;
71 210 : const Real pitch_scale_factor = secondary_side_length / primary_side_length;
72 :
73 : auto mesh = buildSlice(ring_radii,
74 : ring_layers,
75 : ring_radial_biases,
76 : ring_inner_boundary_layer_params,
77 : ring_outer_boundary_layer_params,
78 : ducts_center_dist,
79 : ducts_layers,
80 : duct_radial_biases,
81 : duct_inner_boundary_layer_params,
82 : duct_outer_boundary_layer_params,
83 : virtual_pitch,
84 : num_sectors_per_side,
85 : background_intervals,
86 : background_radial_bias,
87 : background_inner_boundary_layer_params,
88 : background_outer_boundary_layer_params,
89 : node_id_background_meta,
90 : virtual_side_number,
91 : side_index,
92 : azimuthal_tangent,
93 : 0,
94 : quad_center_elements,
95 : center_quad_factor,
96 : false,
97 : true,
98 : 0,
99 : pitch_scale_factor,
100 210 : generate_side_specific_boundaries);
101 210 : MeshTools::Modification::rotate(*mesh, rotation_angle, 0, 0);
102 210 : return mesh;
103 0 : }
104 :
105 : std::unique_ptr<ReplicatedMesh>
106 9819 : PolygonMeshGeneratorBase::buildSimpleSlice(
107 : std::vector<Real> ring_radii,
108 : const std::vector<unsigned int> ring_layers,
109 : const std::vector<Real> ring_radial_biases,
110 : const multiBdryLayerParams & ring_inner_boundary_layer_params,
111 : const multiBdryLayerParams & ring_outer_boundary_layer_params,
112 : std::vector<Real> ducts_center_dist,
113 : const std::vector<unsigned int> ducts_layers,
114 : const std::vector<Real> duct_radial_biases,
115 : const multiBdryLayerParams & duct_inner_boundary_layer_params,
116 : const multiBdryLayerParams & duct_outer_boundary_layer_params,
117 : const Real pitch,
118 : const unsigned int num_sectors_per_side,
119 : const unsigned int background_intervals,
120 : const Real background_radial_bias,
121 : const singleBdryLayerParams & background_inner_boundary_layer_params,
122 : const singleBdryLayerParams & background_outer_boundary_layer_params,
123 : dof_id_type & node_id_background_meta,
124 : const unsigned int side_number,
125 : const unsigned int side_index,
126 : const std::vector<Real> azimuthal_tangent,
127 : const subdomain_id_type block_id_shift,
128 : const bool quad_center_elements,
129 : const Real center_quad_factor,
130 : const bool create_inward_interface_boundaries,
131 : const bool create_outward_interface_boundaries,
132 : const boundary_id_type boundary_id_shift,
133 : const bool generate_side_specific_boundaries,
134 : const TRI_ELEM_TYPE tri_elem_type,
135 : const QUAD_ELEM_TYPE quad_elem_type)
136 : {
137 : return buildSlice(ring_radii,
138 : ring_layers,
139 : ring_radial_biases,
140 : ring_inner_boundary_layer_params,
141 : ring_outer_boundary_layer_params,
142 : ducts_center_dist,
143 : ducts_layers,
144 : duct_radial_biases,
145 : duct_inner_boundary_layer_params,
146 : duct_outer_boundary_layer_params,
147 : pitch,
148 : num_sectors_per_side,
149 : background_intervals,
150 : background_radial_bias,
151 : background_inner_boundary_layer_params,
152 : background_outer_boundary_layer_params,
153 : node_id_background_meta,
154 : side_number,
155 : side_index,
156 : azimuthal_tangent,
157 : block_id_shift,
158 : quad_center_elements,
159 : center_quad_factor,
160 : create_inward_interface_boundaries,
161 : create_outward_interface_boundaries,
162 : boundary_id_shift,
163 : 1.0,
164 : generate_side_specific_boundaries,
165 : tri_elem_type,
166 9819 : quad_elem_type);
167 : }
168 :
169 : std::unique_ptr<ReplicatedMesh>
170 18971 : PolygonMeshGeneratorBase::buildSlice(
171 : std::vector<Real> ring_radii,
172 : const std::vector<unsigned int> ring_layers,
173 : const std::vector<Real> ring_radial_biases,
174 : const multiBdryLayerParams & ring_inner_boundary_layer_params,
175 : const multiBdryLayerParams & ring_outer_boundary_layer_params,
176 : std::vector<Real> ducts_center_dist,
177 : const std::vector<unsigned int> ducts_layers,
178 : const std::vector<Real> duct_radial_biases,
179 : const multiBdryLayerParams & duct_inner_boundary_layer_params,
180 : const multiBdryLayerParams & duct_outer_boundary_layer_params,
181 : const Real pitch,
182 : const unsigned int num_sectors_per_side,
183 : const unsigned int background_intervals,
184 : const Real background_radial_bias,
185 : const singleBdryLayerParams & background_inner_boundary_layer_params,
186 : const singleBdryLayerParams & background_outer_boundary_layer_params,
187 : dof_id_type & node_id_background_meta,
188 : const Real virtual_side_number,
189 : const unsigned int side_index,
190 : const std::vector<Real> azimuthal_tangent,
191 : const subdomain_id_type block_id_shift,
192 : const bool quad_center_elements,
193 : const Real center_quad_factor,
194 : const bool create_inward_interface_boundaries,
195 : const bool create_outward_interface_boundaries,
196 : const boundary_id_type boundary_id_shift,
197 : const Real pitch_scale_factor,
198 : const bool generate_side_specific_boundaries,
199 : const TRI_ELEM_TYPE tri_elem_type,
200 : const QUAD_ELEM_TYPE quad_elem_type)
201 : {
202 18971 : const unsigned short order = quad_elem_type == QUAD_ELEM_TYPE::QUAD4 ? 1 : 2;
203 19908 : if (order != (tri_elem_type == TRI_ELEM_TYPE::TRI3 ? 1 : 2))
204 0 : mooseError("In mesh generator ",
205 : this->name(),
206 : ", an incompatible elements type combination is used when calling "
207 : "PolygonMeshGeneratorBase::buildSlice().");
208 : // In order to create quadratic elements (i.e., order = 2), we creates nodes with double mesh
209 : // density. Thus, the related parameters need to be modified accordingly. A prefix "mod_" is used
210 : // to indicate the modified parameters.
211 :
212 : // For ring_layers, modification is to double the number of layers for order = 2
213 18971 : std::vector<unsigned int> mod_ring_layers(ring_layers);
214 : std::for_each(
215 27023 : mod_ring_layers.begin(), mod_ring_layers.end(), [&order](unsigned int & n) { n *= order; });
216 : // For ring_radial_biases, modification is to take the square root of the original biases for
217 : // order = 2
218 18971 : std::vector<Real> mod_ring_radial_biases(ring_radial_biases);
219 18971 : std::for_each(mod_ring_radial_biases.begin(),
220 : mod_ring_radial_biases.end(),
221 27023 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
222 : // ducts_layers is similar to ring_layers
223 18971 : std::vector<unsigned int> mod_ducts_layers(ducts_layers);
224 : std::for_each(
225 1738 : mod_ducts_layers.begin(), mod_ducts_layers.end(), [&order](unsigned int & n) { n *= order; });
226 : // duct_radial_biases is similar to ring_radial_biases
227 18971 : std::vector<Real> mod_duct_radial_biases(duct_radial_biases);
228 18971 : std::for_each(mod_duct_radial_biases.begin(),
229 : mod_duct_radial_biases.end(),
230 1738 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
231 : // Azimuthal mesh density is also doubled for order = 2
232 18971 : const unsigned int mod_num_sectors_per_side = num_sectors_per_side * order;
233 18971 : const unsigned int mod_background_intervals = background_intervals * order;
234 : // background_radial_bias is similar to ring_radial_biases
235 18971 : const Real mod_background_radial_bias = std::pow(background_radial_bias, 1.0 / order);
236 : // Perform similar modifications for boundary layer parameters
237 : const auto mod_ring_inner_boundary_layer_params =
238 18971 : modifiedMultiBdryLayerParamsCreator(ring_inner_boundary_layer_params, order);
239 : const auto mod_ring_outer_boundary_layer_params =
240 18971 : modifiedMultiBdryLayerParamsCreator(ring_outer_boundary_layer_params, order);
241 : const auto mod_duct_inner_boundary_layer_params =
242 18971 : modifiedMultiBdryLayerParamsCreator(duct_inner_boundary_layer_params, order);
243 : const auto mod_duct_outer_boundary_layer_params =
244 18971 : modifiedMultiBdryLayerParamsCreator(duct_outer_boundary_layer_params, order);
245 :
246 : const auto mod_background_inner_boundary_layer_params =
247 18971 : modifiedSingleBdryLayerParamsCreator(background_inner_boundary_layer_params, order);
248 : const auto mod_background_outer_boundary_layer_params =
249 18971 : modifiedSingleBdryLayerParamsCreator(background_outer_boundary_layer_params, order);
250 :
251 : // The distance parameters of the rings and duct need to be modified too as they may be involved
252 : // in the boundary layer cases.
253 18971 : std::vector<Real> mod_ducts_center_dist(ducts_center_dist);
254 18971 : std::vector<Real> mod_ring_radii(ring_radii);
255 18971 : bool has_rings(ring_radii.size());
256 18971 : bool has_ducts(ducts_center_dist.size());
257 : bool has_background(background_intervals);
258 18971 : auto mesh = buildReplicatedMesh(2);
259 :
260 : // Calculate biasing terms
261 : // background region needs to be split into three parts
262 : const auto main_background_bias_terms =
263 18971 : biasTermsCalculator(background_radial_bias, background_intervals);
264 : const auto inner_background_bias_terms =
265 18971 : biasTermsCalculator(background_inner_boundary_layer_params.bias,
266 18971 : background_inner_boundary_layer_params.intervals);
267 : const auto outer_background_bias_terms =
268 18971 : biasTermsCalculator(background_outer_boundary_layer_params.bias,
269 18971 : background_outer_boundary_layer_params.intervals);
270 : auto rings_bias_terms = biasTermsCalculator(ring_radial_biases,
271 : ring_layers,
272 : ring_inner_boundary_layer_params,
273 18971 : ring_outer_boundary_layer_params);
274 : auto duct_bias_terms = biasTermsCalculator(duct_radial_biases,
275 : ducts_layers,
276 : duct_inner_boundary_layer_params,
277 18971 : duct_outer_boundary_layer_params);
278 : // Equivalent "mod_" parts
279 : const auto mod_main_background_bias_terms =
280 18971 : biasTermsCalculator(mod_background_radial_bias, mod_background_intervals);
281 : const auto mod_inner_background_bias_terms =
282 18971 : biasTermsCalculator(mod_background_inner_boundary_layer_params.bias,
283 18971 : mod_background_inner_boundary_layer_params.intervals);
284 : const auto mod_outer_background_bias_terms =
285 18971 : biasTermsCalculator(mod_background_outer_boundary_layer_params.bias,
286 18971 : mod_background_outer_boundary_layer_params.intervals);
287 : auto mod_rings_bias_terms = biasTermsCalculator(mod_ring_radial_biases,
288 : mod_ring_layers,
289 : mod_ring_inner_boundary_layer_params,
290 18971 : mod_ring_outer_boundary_layer_params);
291 : auto mod_duct_bias_terms = biasTermsCalculator(mod_duct_radial_biases,
292 : mod_ducts_layers,
293 : mod_duct_inner_boundary_layer_params,
294 18971 : mod_duct_outer_boundary_layer_params);
295 :
296 : std::vector<unsigned int> total_ring_layers;
297 45994 : for (unsigned int i = 0; i < ring_layers.size(); i++)
298 27023 : total_ring_layers.push_back(ring_layers[i] + ring_inner_boundary_layer_params.intervals[i] +
299 : ring_outer_boundary_layer_params.intervals[i]);
300 :
301 18971 : if (background_inner_boundary_layer_params.intervals)
302 : {
303 105 : total_ring_layers.push_back(background_inner_boundary_layer_params.intervals);
304 105 : rings_bias_terms.push_back(inner_background_bias_terms);
305 210 : ring_radii.push_back((ring_radii.empty() ? 0.0 : ring_radii.back()) +
306 105 : background_inner_boundary_layer_params.width);
307 : has_rings = true;
308 : }
309 : std::vector<unsigned int> mod_total_ring_layers;
310 45994 : for (unsigned int i = 0; i < mod_ring_layers.size(); i++)
311 27023 : mod_total_ring_layers.push_back(mod_ring_layers[i] +
312 27023 : mod_ring_inner_boundary_layer_params.intervals[i] +
313 : mod_ring_outer_boundary_layer_params.intervals[i]);
314 :
315 18971 : if (mod_background_inner_boundary_layer_params.intervals)
316 : {
317 105 : mod_total_ring_layers.push_back(mod_background_inner_boundary_layer_params.intervals);
318 105 : mod_rings_bias_terms.push_back(mod_inner_background_bias_terms);
319 210 : mod_ring_radii.push_back((mod_ring_radii.empty() ? 0.0 : mod_ring_radii.back()) +
320 105 : mod_background_inner_boundary_layer_params.width);
321 : // has_rings should be modified before in the none "mod_" part
322 : }
323 :
324 : std::vector<unsigned int> total_ducts_layers;
325 18971 : if (background_outer_boundary_layer_params.intervals)
326 : {
327 105 : total_ducts_layers.push_back(background_outer_boundary_layer_params.intervals);
328 105 : duct_bias_terms.insert(duct_bias_terms.begin(), outer_background_bias_terms);
329 105 : ducts_center_dist.insert(ducts_center_dist.begin(),
330 : (ducts_center_dist.empty()
331 105 : ? pitch / 2.0 / std::cos(M_PI / virtual_side_number)
332 105 : : ducts_center_dist.front()) -
333 105 : background_outer_boundary_layer_params.width);
334 : has_ducts = true;
335 : }
336 20709 : for (unsigned int i = 0; i < ducts_layers.size(); i++)
337 1738 : total_ducts_layers.push_back(ducts_layers[i] + duct_inner_boundary_layer_params.intervals[i] +
338 : duct_outer_boundary_layer_params.intervals[i]);
339 :
340 : std::vector<unsigned int> mod_total_ducts_layers;
341 18971 : if (mod_background_outer_boundary_layer_params.intervals)
342 : {
343 105 : mod_total_ducts_layers.push_back(mod_background_outer_boundary_layer_params.intervals);
344 105 : mod_duct_bias_terms.insert(mod_duct_bias_terms.begin(), mod_outer_background_bias_terms);
345 105 : mod_ducts_center_dist.insert(mod_ducts_center_dist.begin(),
346 : (mod_ducts_center_dist.empty()
347 105 : ? pitch / 2.0 / std::cos(M_PI / virtual_side_number)
348 105 : : mod_ducts_center_dist.front()) -
349 105 : mod_background_outer_boundary_layer_params.width);
350 : // has_ducts should be modified before in the none "mod_" part
351 : }
352 20709 : for (unsigned int i = 0; i < mod_ducts_layers.size(); i++)
353 1738 : mod_total_ducts_layers.push_back(mod_ducts_layers[i] +
354 1738 : mod_duct_inner_boundary_layer_params.intervals[i] +
355 : mod_duct_outer_boundary_layer_params.intervals[i]);
356 :
357 : unsigned int angle_number = azimuthal_tangent.size() == 0
358 : ? num_sectors_per_side
359 18971 : : ((azimuthal_tangent.size() - 1) / order);
360 : unsigned int mod_angle_number =
361 18971 : azimuthal_tangent.size() == 0 ? mod_num_sectors_per_side : (azimuthal_tangent.size() - 1);
362 :
363 : // Geometries
364 18971 : const Real corner_to_corner =
365 18971 : pitch / std::cos(M_PI / virtual_side_number); // distance of bin center to cell corner
366 18971 : const Real corner_p[2][2] = {
367 18971 : {0.0, 0.5 * corner_to_corner},
368 18971 : {0.5 * corner_to_corner * pitch_scale_factor * std::sin(2.0 * M_PI / virtual_side_number),
369 18971 : 0.5 * corner_to_corner * pitch_scale_factor * std::cos(2.0 * M_PI / virtual_side_number)}};
370 18971 : const unsigned int div_num = angle_number / 2 + 1;
371 18971 : const unsigned int mod_div_num = mod_angle_number / 2 + 1;
372 :
373 : // From now on, we work on the nodes, which need the "mod_" parameters
374 18971 : std::vector<std::vector<Node *>> nodes(mod_div_num, std::vector<Node *>(mod_div_num));
375 18971 : if (quad_center_elements)
376 : {
377 : Real ring_radii_0;
378 :
379 2249 : if (has_rings)
380 1449 : ring_radii_0 = ring_radii.front() * mod_rings_bias_terms.front()[order - 1];
381 800 : else if (has_ducts)
382 84 : ring_radii_0 = mod_ducts_center_dist.front() * std::cos(M_PI / virtual_side_number) *
383 84 : mod_main_background_bias_terms[order - 1];
384 : else
385 716 : ring_radii_0 = pitch / 2.0 * mod_main_background_bias_terms[order - 1];
386 : // If center_quad_factor is zero, default value (div_num - 1)/div_num is used.
387 : // We use div_num instead of mod_div_num because we are dealing wth elements here
388 : // This approach ensures that the order = 2 mesh elements are consistent with the order = 1
389 2249 : ring_radii_0 *=
390 2249 : center_quad_factor == 0.0 ? (((Real)div_num - 1.0) / (Real)div_num) : center_quad_factor;
391 :
392 2249 : centerNodes(*mesh, virtual_side_number, mod_div_num, ring_radii_0, nodes);
393 : }
394 : else // pin-cell center
395 16722 : mesh->add_point(Point(0.0, 0.0, 0.0));
396 :
397 : // create nodes for the ring regions
398 18971 : if (has_rings)
399 16477 : ringNodes(*mesh,
400 : ring_radii,
401 : mod_total_ring_layers,
402 : mod_rings_bias_terms,
403 : mod_num_sectors_per_side,
404 : corner_p,
405 : corner_to_corner,
406 : azimuthal_tangent);
407 :
408 18971 : if (has_background)
409 : {
410 : // add nodes in background region; the background region is defined as the area between the
411 : // outermost pin (if there is a pin; if no pin, the center) and the innermost hex/duct; if
412 : // _has_ducts is false, the background region is the area between the pin and enclosing hexagon
413 : Real background_corner_radial_interval_length;
414 : Real background_corner_distance;
415 : Real background_in;
416 : Real background_out; // background outer frontier
417 10029 : if (has_rings)
418 7535 : background_in = ring_radii.back();
419 : else
420 : background_in = 0;
421 :
422 10029 : if (has_ducts)
423 : {
424 1416 : background_out = mod_ducts_center_dist.front();
425 : background_corner_distance =
426 : mod_ducts_center_dist
427 : .front(); // it is the center to duct (innermost duct) corner distance
428 : }
429 : else
430 : {
431 : background_out = 0.5 * corner_to_corner;
432 : background_corner_distance =
433 : 0.5 * corner_to_corner; // it is the center to hex corner distance
434 : }
435 :
436 10029 : background_corner_radial_interval_length =
437 10029 : (background_out - background_in) / mod_background_intervals;
438 :
439 10029 : node_id_background_meta = mesh->n_nodes();
440 :
441 : // create nodes for background region
442 10029 : backgroundNodes(*mesh,
443 : mod_num_sectors_per_side,
444 : mod_background_intervals,
445 : mod_main_background_bias_terms,
446 : background_corner_distance,
447 : background_corner_radial_interval_length,
448 : corner_p,
449 : corner_to_corner,
450 : background_in,
451 : azimuthal_tangent);
452 : }
453 :
454 : // create nodes for duct regions
455 18971 : if (has_ducts)
456 1416 : ductNodes(*mesh,
457 : &mod_ducts_center_dist,
458 : mod_total_ducts_layers,
459 : mod_duct_bias_terms,
460 : mod_num_sectors_per_side,
461 : corner_p,
462 : corner_to_corner,
463 : azimuthal_tangent);
464 :
465 : // See if the central region is the only part of the innermost part
466 : // The central region of the slice is special.
467 : // Unlike the outer regions, which are layered quad elements,
468 : // the central region is either a layer of tri elements or a specially-patterned quad elements.
469 : // If there is at least one `ring` defined in the slice,
470 : // the central region must belong to the innermost (first) ring.
471 : // Otherwise the central region belongs to the `background`
472 : // In either case, if the innermost ring or background has only one radial interval,
473 : // the central region is an independent ring or background
474 : // Otherwise, the central region and one or several quad element layers together form the
475 : // innermost ring or background
476 : bool is_central_region_independent;
477 18971 : if (ring_layers.empty())
478 2529 : is_central_region_independent = mod_background_inner_boundary_layer_params.intervals +
479 2529 : mod_background_intervals +
480 : mod_background_outer_boundary_layer_params.intervals ==
481 : 1;
482 : else
483 16442 : is_central_region_independent = mod_ring_layers[0] +
484 16442 : mod_ring_inner_boundary_layer_params.intervals[0] +
485 : mod_ring_outer_boundary_layer_params.intervals[0] ==
486 : 1;
487 :
488 : // From now on, we work on the elements, which need the none "mod_" parameters
489 : // Assign elements, boundaries, and subdomains;
490 : // Add Tri3/Tri6/Tri7 or Quad4/Quad8/Quad9 mesh into innermost (central) region
491 18971 : if (quad_center_elements)
492 4498 : cenQuadElemDef(*mesh,
493 : div_num,
494 : block_id_shift,
495 2249 : create_outward_interface_boundaries && is_central_region_independent,
496 : boundary_id_shift,
497 : nodes,
498 2249 : (!has_rings) && (!has_ducts) && (background_intervals == 1),
499 : // Note here, has_ring means either there are ring regions or background inner
500 : // boundary layer; has_ducts means either there are duct regions or background
501 : // outer boundary layer. Same in cenTriElemDef()
502 : side_index,
503 : generate_side_specific_boundaries,
504 : quad_elem_type);
505 : else
506 16722 : cenTriElemDef(
507 : *mesh,
508 : num_sectors_per_side,
509 : azimuthal_tangent,
510 : block_id_shift,
511 16722 : create_outward_interface_boundaries && is_central_region_independent,
512 : boundary_id_shift,
513 16722 : ((!has_rings) && (!has_ducts) && (background_intervals == 1)) ||
514 8942 : ((!has_background) &&
515 : (std::accumulate(total_ring_layers.begin(), total_ring_layers.end(), 0) == 1)),
516 : // Only for ACCG, it is possible that the entire mesh is a single-layer ring.
517 : // cenQuadElemDef() does not need this as it does not work for ACCG.
518 : side_index,
519 : generate_side_specific_boundaries,
520 : tri_elem_type);
521 :
522 : // Add Quad4 mesh into outer circle
523 : // total number of mesh should be all the rings for pin regions + background regions;
524 : // total number of quad mesh should be total number of mesh -1 (-1 is because the inner circle for
525 : // tri/quad mesh has been added above)
526 :
527 : std::vector<unsigned int> subdomain_rings;
528 18971 : if (has_rings) // define the rings in each subdomain
529 : {
530 16477 : subdomain_rings = total_ring_layers;
531 16477 : subdomain_rings.front() -= 1; // remove the inner TRI mesh subdomain
532 16477 : if (background_inner_boundary_layer_params.intervals)
533 : {
534 105 : subdomain_rings.back() =
535 105 : background_inner_boundary_layer_params.intervals + background_intervals +
536 105 : background_outer_boundary_layer_params.intervals; // add the background region
537 105 : if (ring_radii.size() == 1)
538 35 : subdomain_rings.back() -= 1; // remove the inner TRI mesh subdomain
539 : }
540 16372 : else if (has_background)
541 7430 : subdomain_rings.push_back(background_inner_boundary_layer_params.intervals +
542 7430 : background_intervals +
543 7430 : background_outer_boundary_layer_params.intervals);
544 : }
545 : else
546 : {
547 : subdomain_rings.push_back(
548 2494 : background_inner_boundary_layer_params.intervals + background_intervals +
549 2494 : background_outer_boundary_layer_params.intervals); // add the background region
550 2494 : subdomain_rings[0] -= 1; // remove the inner TRI mesh subdomain
551 : }
552 :
553 18971 : if (has_ducts)
554 3154 : for (unsigned int i = (background_outer_boundary_layer_params.intervals > 0);
555 3154 : i < total_ducts_layers.size();
556 : i++)
557 1738 : subdomain_rings.push_back(total_ducts_layers[i]);
558 :
559 21220 : quadElemDef(*mesh,
560 : num_sectors_per_side,
561 : subdomain_rings,
562 : side_index,
563 : azimuthal_tangent,
564 : block_id_shift,
565 2249 : quad_center_elements ? (mod_div_num * mod_div_num - 1) : 0,
566 : create_inward_interface_boundaries,
567 : create_outward_interface_boundaries,
568 : boundary_id_shift,
569 : generate_side_specific_boundaries,
570 : quad_elem_type);
571 18971 : if (tri_elem_type == TRI_ELEM_TYPE::TRI6 || quad_elem_type == QUAD_ELEM_TYPE::QUAD8)
572 579 : mesh->remove_orphaned_nodes();
573 18971 : return mesh;
574 18971 : }
575 :
576 : void
577 2249 : PolygonMeshGeneratorBase::centerNodes(ReplicatedMesh & mesh,
578 : const Real virtual_side_number,
579 : const unsigned int div_num,
580 : const Real ring_radii_0,
581 : std::vector<std::vector<Node *>> & nodes) const
582 : {
583 : const std::pair<Real, Real> p_origin = std::make_pair(0.0, 0.0);
584 : const std::pair<Real, Real> p_bottom =
585 2249 : std::make_pair(0.0, ring_radii_0 * std::cos(M_PI / virtual_side_number));
586 : const std::pair<Real, Real> p_top =
587 2249 : std::make_pair(p_bottom.second * std::sin(2.0 * M_PI / virtual_side_number),
588 2249 : p_bottom.second * std::cos(2.0 * M_PI / virtual_side_number));
589 : const std::pair<Real, Real> p_diag =
590 2249 : std::make_pair(ring_radii_0 * std::sin(M_PI / virtual_side_number),
591 : ring_radii_0 * std::cos(M_PI / virtual_side_number));
592 :
593 : // The four vertices of the central quad region are defined above.
594 : // The following loops transverse all the nodes within this central quad region by moving p1 thru
595 : // p4 and calculate the four-point intercept (pc).
596 : // p_top------o-------p4--------o-----p_diag
597 : // | | | | |
598 : // | | | | |
599 : // o--------o--------o--------o--------o
600 : // | | | | |
601 : // | | | | |
602 : // p1--------o-------pc--------o-------p2
603 : // | | | | |
604 : // | | | | |
605 : // o--------o--------o--------o--------o
606 : // | | | | |
607 : // | | | | |
608 : // p_origin-----o-------p3--------o----p_bottom
609 : //
610 : // The loops are designed to transverse the nodes as shown below to facilitate elements
611 : // and sides creation.
612 : //
613 : // 25-------24-------23-------22-------21
614 : // | | | | |
615 : // | | | | |
616 : // 16-------15-------14-------13-------20
617 : // | | | | |
618 : // | | | | |
619 : // 9--------8------- 7-------12-------19
620 : // | | | | |
621 : // | | | | |
622 : // 4--------3--------6-------11-------18
623 : // | | | | |
624 : // | | | | |
625 : // 1--------2--------5-------10-------17
626 :
627 7228 : for (unsigned int i = 0; i < div_num; i++)
628 : {
629 : unsigned int id_x = 0;
630 : unsigned int id_y = i;
631 16800 : for (unsigned int j = 0; j < 2 * i + 1; j++)
632 : {
633 11821 : std::pair<Real, Real> p1 = std::make_pair(
634 11821 : (p_origin.first * (div_num - 1 - id_x) + p_top.first * id_x) / (div_num - 1),
635 11821 : (p_origin.second * (div_num - 1 - id_x) + p_top.second * id_x) / (div_num - 1));
636 11821 : std::pair<Real, Real> p2 = std::make_pair(
637 11821 : (p_bottom.first * (div_num - 1 - id_x) + p_diag.first * id_x) / (div_num - 1),
638 11821 : (p_bottom.second * (div_num - 1 - id_x) + p_diag.second * id_x) / (div_num - 1));
639 11821 : std::pair<Real, Real> p3 = std::make_pair(
640 11821 : (p_origin.first * (div_num - 1 - id_y) + p_bottom.first * id_y) / (div_num - 1),
641 11821 : (p_origin.second * (div_num - 1 - id_y) + p_bottom.second * id_y) / (div_num - 1));
642 11821 : std::pair<Real, Real> p4 = std::make_pair(
643 11821 : (p_top.first * (div_num - 1 - id_y) + p_diag.first * id_y) / (div_num - 1),
644 11821 : (p_top.second * (div_num - 1 - id_y) + p_diag.second * id_y) / (div_num - 1));
645 11821 : std::pair<Real, Real> pc = fourPointIntercept(p1, p2, p3, p4);
646 11821 : nodes[id_x][id_y] = mesh.add_point(Point(pc.first, pc.second, 0.0));
647 11821 : if (j < i)
648 3421 : id_x++;
649 11821 : if (j >= i)
650 8400 : id_y--;
651 : }
652 : }
653 2249 : }
654 :
655 : void
656 16477 : PolygonMeshGeneratorBase::ringNodes(ReplicatedMesh & mesh,
657 : const std::vector<Real> ring_radii,
658 : const std::vector<unsigned int> ring_layers,
659 : const std::vector<std::vector<Real>> biased_terms,
660 : const unsigned int num_sectors_per_side,
661 : const Real corner_p[2][2],
662 : const Real corner_to_corner,
663 : const std::vector<Real> azimuthal_tangent) const
664 : {
665 : const unsigned int angle_number =
666 16477 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
667 :
668 : // Add nodes in pins regions
669 43605 : for (unsigned int l = 0; l < ring_layers.size(); l++)
670 : {
671 : // the pin radius interval for each ring_radii/subdomain
672 : const Real pin_radius_interval_length =
673 27128 : l == 0 ? ring_radii[l] / ring_layers[l]
674 10651 : : (ring_radii[l] - ring_radii[l - 1]) / ring_layers[l];
675 :
676 : // add rings in each pin subdomain
677 126596 : for (unsigned int k = 0; k < ring_layers[l]; k++)
678 : {
679 : const Real bin_radial_distance =
680 99468 : l == 0 ? (biased_terms[l][k] * ring_layers[l] *
681 : pin_radius_interval_length) // this is from the cell/pin center to
682 : // the first circle
683 13538 : : (ring_radii[l - 1] +
684 13538 : biased_terms[l][k] * ring_layers[l] * pin_radius_interval_length);
685 99468 : const Real pin_corner_p_x = corner_p[0][0] * bin_radial_distance / (0.5 * corner_to_corner);
686 99468 : const Real pin_corner_p_y = corner_p[0][1] * bin_radial_distance / (0.5 * corner_to_corner);
687 :
688 : // pin_corner_p(s) are the points in the pin region, on the bins towards the six corners,
689 : // at different intervals
690 99468 : mesh.add_point(Point(pin_corner_p_x, pin_corner_p_y, 0.0));
691 :
692 239632 : for (unsigned int j = 1; j <= angle_number; j++)
693 : {
694 : const Real cell_boundary_p_x =
695 140164 : corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
696 140164 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
697 1440 : : (azimuthal_tangent[j] / 2.0));
698 : const Real cell_boundary_p_y =
699 140164 : corner_p[0][1] + (corner_p[1][1] - corner_p[0][1]) *
700 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
701 : : (azimuthal_tangent[j] / 2.0));
702 : // cell_boundary_p(s) are the points on the cell's six boundaries (flat sides) at
703 : // different azimuthal angles
704 : const Real pin_azimuthal_p_x =
705 140164 : cell_boundary_p_x * bin_radial_distance /
706 140164 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
707 : const Real pin_azimuthal_p_y =
708 140164 : cell_boundary_p_y * bin_radial_distance /
709 140164 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
710 :
711 : // pin_azimuthal_p are the points on the bins towards different azimuthal angles, at
712 : // different intervals; excluding the ones produced by pin_corner_p
713 140164 : mesh.add_point(Point(pin_azimuthal_p_x, pin_azimuthal_p_y, 0.0));
714 : }
715 : }
716 : }
717 16477 : }
718 :
719 : void
720 10029 : PolygonMeshGeneratorBase::backgroundNodes(ReplicatedMesh & mesh,
721 : const unsigned int num_sectors_per_side,
722 : const unsigned int background_intervals,
723 : const std::vector<Real> biased_terms,
724 : const Real background_corner_distance,
725 : const Real background_corner_radial_interval_length,
726 : const Real corner_p[2][2],
727 : const Real corner_to_corner,
728 : const Real background_in,
729 : const std::vector<Real> azimuthal_tangent) const
730 : {
731 : unsigned int angle_number =
732 10029 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
733 24407 : for (unsigned int k = 0; k < (background_intervals); k++)
734 : {
735 : const Real background_corner_p_x =
736 14378 : background_corner_distance / (0.5 * corner_to_corner) * corner_p[0][0] *
737 14378 : (background_in +
738 14378 : biased_terms[k] * background_intervals * background_corner_radial_interval_length) /
739 14378 : background_corner_distance;
740 : const Real background_corner_p_y =
741 14378 : background_corner_distance / (0.5 * corner_to_corner) * corner_p[0][1] *
742 : (background_in +
743 : biased_terms[k] * background_intervals * background_corner_radial_interval_length) /
744 14378 : background_corner_distance;
745 :
746 : // background_corner_p(s) are the points in the background region, on the bins towards the six
747 : // corners, at different intervals
748 14378 : mesh.add_point(Point(background_corner_p_x, background_corner_p_y, 0.0));
749 :
750 65215 : for (unsigned int j = 1; j <= angle_number; j++)
751 : {
752 : const Real cell_boundary_p_x =
753 50837 : background_corner_distance / (0.5 * corner_to_corner) *
754 50837 : (corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
755 50837 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
756 3299 : : (azimuthal_tangent[j] / 2.0)));
757 : const Real cell_boundary_p_y =
758 50837 : background_corner_distance / (0.5 * corner_to_corner) *
759 50837 : (corner_p[0][1] + (corner_p[1][1] - corner_p[0][1]) *
760 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
761 : : (azimuthal_tangent[j] / 2.0)));
762 : // cell_boundary_p(s) are the points on the cell's six boundaries (flat sides) at different
763 : // azimuthal angles
764 : const Real pin_boundary_p_x =
765 50837 : cell_boundary_p_x * background_in /
766 50837 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
767 : const Real pin_boundary_p_y =
768 50837 : cell_boundary_p_y * background_in /
769 50837 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
770 : // pin_boundary_p(s) are the points on pin boundary (outside ring) at different azimuthal
771 : // angles
772 : const Real background_radial_interval =
773 50837 : std::sqrt(Utility::pow<2>(cell_boundary_p_x - pin_boundary_p_x) +
774 50837 : Utility::pow<2>(cell_boundary_p_y - pin_boundary_p_y)) /
775 50837 : background_intervals;
776 : const Real background_azimuthal_p_x =
777 50837 : cell_boundary_p_x *
778 50837 : (background_in + biased_terms[k] * background_intervals * background_radial_interval) /
779 50837 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
780 : const Real background_azimuthal_p_y =
781 50837 : cell_boundary_p_y *
782 : (background_in + biased_terms[k] * background_intervals * background_radial_interval) /
783 50837 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
784 : // background_azimuthal_p are the points on the bins towards different azimuthal angles, at
785 : // different intervals; excluding the ones produced by background_corner_p
786 50837 : mesh.add_point(Point(background_azimuthal_p_x, background_azimuthal_p_y, 0.0));
787 : }
788 : }
789 10029 : }
790 :
791 : void
792 1416 : PolygonMeshGeneratorBase::ductNodes(ReplicatedMesh & mesh,
793 : std::vector<Real> * const ducts_center_dist,
794 : const std::vector<unsigned int> ducts_layers,
795 : const std::vector<std::vector<Real>> biased_terms,
796 : const unsigned int num_sectors_per_side,
797 : const Real corner_p[2][2],
798 : const Real corner_to_corner,
799 : const std::vector<Real> azimuthal_tangent) const
800 : {
801 : unsigned int angle_number =
802 1416 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
803 : // Add nodes in ducts regions
804 : (*ducts_center_dist)
805 1416 : .push_back(0.5 * corner_to_corner); // add hex boundary as the last element in this vector
806 1416 : std::vector<Real> duct_radius_interval_length(ducts_layers.size());
807 :
808 : Real bin_radial_distance;
809 3259 : for (unsigned int l = 0; l < ducts_layers.size(); l++)
810 : {
811 1843 : duct_radius_interval_length[l] =
812 1843 : ((*ducts_center_dist)[l + 1] - (*ducts_center_dist)[l]) /
813 : ducts_layers[l]; // the pin radius interval for each ring_radii/subdomain
814 :
815 : // add rings in each pin subdomain
816 6038 : for (unsigned int k = 0; k < ducts_layers[l]; k++)
817 : {
818 4195 : bin_radial_distance = ((*ducts_center_dist)[l] +
819 4195 : biased_terms[l][k] * ducts_layers[l] * duct_radius_interval_length[l]);
820 4195 : const Real pin_corner_p_x = corner_p[0][0] * bin_radial_distance / (0.5 * corner_to_corner);
821 4195 : const Real pin_corner_p_y = corner_p[0][1] * bin_radial_distance / (0.5 * corner_to_corner);
822 :
823 : // pin_corner_p(s) are the points in the pin region, on the bins towards the six corners,
824 : // at different intervals
825 4195 : mesh.add_point(Point(pin_corner_p_x, pin_corner_p_y, 0.0));
826 :
827 24247 : for (unsigned int j = 1; j <= angle_number; j++)
828 : {
829 : const Real cell_boundary_p_x =
830 20052 : corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
831 20052 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
832 0 : : (azimuthal_tangent[j] / 2.0));
833 : const Real cell_boundary_p_y =
834 20052 : corner_p[0][1] + (corner_p[1][1] - corner_p[0][1]) *
835 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
836 20052 : : (azimuthal_tangent[j] / 2.0));
837 : // cell_boundary_p(s) are the points on the cell's six boundaries (flat sides) at
838 : // different azimuthal angles
839 20052 : const Real pin_azimuthal_p_x =
840 20052 : cell_boundary_p_x * bin_radial_distance / (0.5 * corner_to_corner);
841 20052 : const Real pin_azimuthal_p_y =
842 20052 : cell_boundary_p_y * bin_radial_distance / (0.5 * corner_to_corner);
843 :
844 : // pin_azimuthal_p are the points on the bins towards different azimuthal angles, at
845 : // different intervals; excluding the ones produced by pin_corner_p
846 20052 : mesh.add_point(Point(pin_azimuthal_p_x, pin_azimuthal_p_y, 0.0));
847 : }
848 : }
849 : }
850 1416 : }
851 :
852 : void
853 2249 : PolygonMeshGeneratorBase::cenQuadElemDef(ReplicatedMesh & mesh,
854 : const unsigned int div_num,
855 : const subdomain_id_type block_id_shift,
856 : const bool create_outward_interface_boundaries,
857 : const boundary_id_type boundary_id_shift,
858 : std::vector<std::vector<Node *>> & nodes,
859 : const bool assign_external_boundary,
860 : const unsigned int side_index,
861 : const bool generate_side_specific_boundaries,
862 : const QUAD_ELEM_TYPE quad_elem_type) const
863 : {
864 :
865 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
866 :
867 : // This loop defines quad elements for the central regions except for the outermost layer
868 4839 : for (unsigned int i = 0; i < div_num - 1; i++)
869 : {
870 : unsigned int id_x = 0;
871 : unsigned int id_y = i;
872 5862 : for (unsigned int j = 0; j < 2 * i + 1; j++)
873 : {
874 3272 : std::unique_ptr<Elem> new_elem;
875 3272 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
876 : {
877 2992 : new_elem = std::make_unique<Quad4>();
878 2992 : new_elem->set_node(0, nodes[id_x][id_y]);
879 2992 : new_elem->set_node(3, nodes[id_x][id_y + 1]);
880 2992 : new_elem->set_node(2, nodes[id_x + 1][id_y + 1]);
881 2992 : new_elem->set_node(1, nodes[id_x + 1][id_y]);
882 2992 : new_elem->subdomain_id() = 1 + block_id_shift;
883 : }
884 : else // QUAD8/QUAD9
885 : {
886 280 : new_elem = std::make_unique<Quad8>();
887 280 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
888 : {
889 280 : new_elem = std::make_unique<Quad9>();
890 280 : new_elem->set_node(8, nodes[id_x * 2 + 1][id_y * 2 + 1]);
891 : }
892 280 : new_elem->set_node(0, nodes[id_x * 2][id_y * 2]);
893 280 : new_elem->set_node(3, nodes[id_x * 2][id_y * 2 + 2]);
894 280 : new_elem->set_node(2, nodes[id_x * 2 + 2][id_y * 2 + 2]);
895 280 : new_elem->set_node(1, nodes[id_x * 2 + 2][id_y * 2]);
896 280 : new_elem->set_node(4, nodes[id_x * 2 + 1][id_y * 2]);
897 280 : new_elem->set_node(5, nodes[id_x * 2 + 2][id_y * 2 + 1]);
898 280 : new_elem->set_node(6, nodes[id_x * 2 + 1][id_y * 2 + 2]);
899 280 : new_elem->set_node(7, nodes[id_x * 2][id_y * 2 + 1]);
900 280 : new_elem->subdomain_id() = 1 + block_id_shift;
901 : }
902 3272 : Elem * elem_Quad = mesh.add_elem(std::move(new_elem));
903 :
904 3272 : if (id_x == 0)
905 2590 : boundary_info.add_side(elem_Quad, 3, SLICE_BEGIN);
906 3272 : if (id_y == 0)
907 2590 : boundary_info.add_side(elem_Quad, 0, SLICE_END);
908 3272 : if (j < i)
909 341 : id_x++;
910 3272 : if (j >= i)
911 2931 : id_y--;
912 3272 : }
913 : }
914 : // This loop defines the outermost layer quad elements of the central region
915 7429 : for (unsigned int i = (div_num - 1) * (div_num - 1); i < div_num * div_num - 1; i++)
916 : {
917 5180 : std::unique_ptr<Elem> new_elem;
918 5180 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
919 : {
920 4900 : new_elem = std::make_unique<Quad4>();
921 4900 : new_elem->set_node(0, mesh.node_ptr(i));
922 4900 : new_elem->set_node(3, mesh.node_ptr(i + 2 * div_num - 1));
923 4900 : new_elem->set_node(2, mesh.node_ptr(i + 2 * div_num));
924 4900 : new_elem->set_node(1, mesh.node_ptr(i + 1));
925 : }
926 : else // QUAD8/QUAD9
927 : {
928 280 : new_elem = std::make_unique<Quad8>();
929 280 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
930 : {
931 280 : new_elem = std::make_unique<Quad9>();
932 280 : new_elem->set_node(8,
933 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
934 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1 +
935 : ((div_num - 1) * 4 + 1)));
936 : }
937 280 : new_elem->set_node(0,
938 280 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
939 : (i - (div_num - 1) * (div_num - 1)) * 2));
940 280 : new_elem->set_node(3,
941 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
942 280 : (i - (div_num - 1) * (div_num - 1)) * 2 +
943 : ((div_num - 1) * 4 + 1) * 2));
944 280 : new_elem->set_node(2,
945 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
946 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2 +
947 : ((div_num - 1) * 4 + 1) * 2));
948 280 : new_elem->set_node(1,
949 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
950 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2));
951 280 : new_elem->set_node(4,
952 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
953 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1));
954 280 : new_elem->set_node(5,
955 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
956 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2 +
957 : ((div_num - 1) * 4 + 1)));
958 280 : new_elem->set_node(6,
959 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
960 280 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1 +
961 : ((div_num - 1) * 4 + 1) * 2));
962 280 : new_elem->set_node(7,
963 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
964 280 : (i - (div_num - 1) * (div_num - 1)) * 2 +
965 : ((div_num - 1) * 4 + 1)));
966 : }
967 :
968 5180 : Elem * elem_Quad = mesh.add_elem(std::move(new_elem));
969 5180 : elem_Quad->subdomain_id() = 1 + block_id_shift;
970 5180 : if (create_outward_interface_boundaries)
971 952 : boundary_info.add_side(elem_Quad, 2, 1 + boundary_id_shift);
972 5180 : if (i == (div_num - 1) * (div_num - 1))
973 2249 : boundary_info.add_side(elem_Quad, 3, SLICE_BEGIN);
974 5180 : if (i == div_num * div_num - 2)
975 2249 : boundary_info.add_side(elem_Quad, 1, SLICE_END);
976 5180 : if (assign_external_boundary)
977 : {
978 280 : boundary_info.add_side(elem_Quad, 2, OUTER_SIDESET_ID);
979 280 : if (generate_side_specific_boundaries)
980 280 : boundary_info.add_side(
981 : elem_Quad,
982 : 2,
983 420 : (i < div_num * (div_num - 1) ? OUTER_SIDESET_ID : OUTER_SIDESET_ID_ALT) + side_index);
984 : }
985 5180 : }
986 2249 : }
987 :
988 : void
989 16722 : PolygonMeshGeneratorBase::cenTriElemDef(ReplicatedMesh & mesh,
990 : const unsigned int num_sectors_per_side,
991 : const std::vector<Real> azimuthal_tangent,
992 : const subdomain_id_type block_id_shift,
993 : const bool create_outward_interface_boundaries,
994 : const boundary_id_type boundary_id_shift,
995 : const bool assign_external_boundary,
996 : const unsigned int side_index,
997 : const bool generate_side_specific_boundaries,
998 : const TRI_ELEM_TYPE tri_elem_type) const
999 : {
1000 16722 : const unsigned short order = tri_elem_type == TRI_ELEM_TYPE::TRI3 ? 1 : 2;
1001 : unsigned int angle_number = azimuthal_tangent.size() == 0
1002 : ? num_sectors_per_side
1003 16722 : : ((azimuthal_tangent.size() - 1) / order);
1004 :
1005 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1006 47769 : for (unsigned int i = 1; i <= angle_number; i++)
1007 : {
1008 31047 : std::unique_ptr<Elem> new_elem;
1009 31047 : if (tri_elem_type == TRI_ELEM_TYPE::TRI3)
1010 : {
1011 28652 : new_elem = std::make_unique<Tri3>();
1012 28652 : new_elem->set_node(0, mesh.node_ptr(0));
1013 28652 : new_elem->set_node(2, mesh.node_ptr(i));
1014 28652 : new_elem->set_node(1, mesh.node_ptr(i + 1));
1015 : }
1016 : else // TRI6/TRI7
1017 : {
1018 2395 : new_elem = std::make_unique<Tri6>();
1019 2395 : if (tri_elem_type == TRI_ELEM_TYPE::TRI7)
1020 : {
1021 1031 : new_elem = std::make_unique<Tri7>();
1022 1031 : new_elem->set_node(6, mesh.node_ptr(i * 2));
1023 : }
1024 2395 : new_elem->set_node(0, mesh.node_ptr(0));
1025 2395 : new_elem->set_node(2, mesh.node_ptr(i * 2 + angle_number * order));
1026 2395 : new_elem->set_node(1, mesh.node_ptr((i + 1) * 2 + angle_number * order));
1027 2395 : new_elem->set_node(3, mesh.node_ptr(i * 2 + 1));
1028 2395 : new_elem->set_node(5, mesh.node_ptr(i * 2 - 1));
1029 2395 : new_elem->set_node(4, mesh.node_ptr(i * 2 + 1 + angle_number * order));
1030 : }
1031 :
1032 31047 : Elem * elem = mesh.add_elem(std::move(new_elem));
1033 31047 : if (create_outward_interface_boundaries)
1034 8317 : boundary_info.add_side(elem, 1, 1 + boundary_id_shift);
1035 31047 : elem->subdomain_id() = 1 + block_id_shift;
1036 31047 : if (i == 1)
1037 16722 : boundary_info.add_side(elem, 2, SLICE_BEGIN);
1038 31047 : if (i == angle_number)
1039 16722 : boundary_info.add_side(elem, 0, SLICE_END);
1040 31047 : if (assign_external_boundary)
1041 : {
1042 3364 : boundary_info.add_side(elem, 1, OUTER_SIDESET_ID);
1043 3364 : if (generate_side_specific_boundaries)
1044 1015 : boundary_info.add_side(elem,
1045 : 1,
1046 1890 : (i <= angle_number / 2 ? OUTER_SIDESET_ID : OUTER_SIDESET_ID_ALT) +
1047 : side_index);
1048 : }
1049 31047 : }
1050 16722 : }
1051 :
1052 : void
1053 18971 : PolygonMeshGeneratorBase::quadElemDef(ReplicatedMesh & mesh,
1054 : const unsigned int num_sectors_per_side,
1055 : const std::vector<unsigned int> subdomain_rings,
1056 : const unsigned int side_index,
1057 : const std::vector<Real> azimuthal_tangent,
1058 : const subdomain_id_type block_id_shift,
1059 : const dof_id_type nodeid_shift,
1060 : const bool create_inward_interface_boundaries,
1061 : const bool create_outward_interface_boundaries,
1062 : const boundary_id_type boundary_id_shift,
1063 : const bool generate_side_specific_boundaries,
1064 : const QUAD_ELEM_TYPE quad_elem_type) const
1065 : {
1066 18971 : const unsigned short order = quad_elem_type == QUAD_ELEM_TYPE::QUAD4 ? 1 : 2;
1067 : unsigned int angle_number = azimuthal_tangent.size() == 0
1068 : ? num_sectors_per_side
1069 18971 : : ((azimuthal_tangent.size() - 1) / order);
1070 :
1071 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1072 : unsigned int j = 0;
1073 57761 : for (unsigned int k = 0; k < (subdomain_rings.size()); k++)
1074 : {
1075 134411 : for (unsigned int m = 0; m < subdomain_rings[k]; m++)
1076 : {
1077 236631 : for (unsigned int i = 1; i <= angle_number; i++)
1078 : {
1079 141010 : std::unique_ptr<Elem> new_elem;
1080 141010 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
1081 : {
1082 132413 : new_elem = std::make_unique<Quad4>();
1083 132413 : new_elem->set_node(0, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * j));
1084 132413 : new_elem->set_node(1, mesh.node_ptr(nodeid_shift + i + 1 + (angle_number + 1) * j));
1085 132413 : new_elem->set_node(2, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * (j + 1) + 1));
1086 132413 : new_elem->set_node(3, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * (j + 1)));
1087 : }
1088 : else // QUAD8/QUAD9
1089 : {
1090 8597 : new_elem = std::make_unique<Quad8>();
1091 8597 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
1092 : {
1093 7373 : new_elem = std::make_unique<Quad9>();
1094 7373 : new_elem->set_node(
1095 7373 : 8, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 2)));
1096 : }
1097 8597 : new_elem->set_node(
1098 : 0,
1099 8597 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 1)));
1100 8597 : new_elem->set_node(
1101 8597 : 1, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 1)));
1102 8597 : new_elem->set_node(
1103 8597 : 2, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 3)));
1104 8597 : new_elem->set_node(
1105 : 3,
1106 8597 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 3)));
1107 8597 : new_elem->set_node(
1108 : 4, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 1)));
1109 8597 : new_elem->set_node(
1110 8597 : 5, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 2)));
1111 8597 : new_elem->set_node(
1112 : 6, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 3)));
1113 8597 : new_elem->set_node(
1114 : 7,
1115 8597 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 2)));
1116 : }
1117 141010 : Elem * elem = mesh.add_elem(std::move(new_elem));
1118 141010 : if (i == 1)
1119 95621 : boundary_info.add_side(elem, 3, SLICE_BEGIN);
1120 141010 : if (i == angle_number)
1121 95621 : boundary_info.add_side(elem, 1, SLICE_END);
1122 :
1123 141010 : if (subdomain_rings[0] == 0)
1124 18496 : elem->subdomain_id() = k + 1 + block_id_shift;
1125 : else
1126 122514 : elem->subdomain_id() = k + 2 + block_id_shift;
1127 :
1128 141010 : if (m == 0 && create_inward_interface_boundaries && k > 0)
1129 560 : boundary_info.add_side(elem, 0, k * 2 + boundary_id_shift);
1130 141010 : if (m == (subdomain_rings[k] - 1))
1131 : {
1132 62182 : if (k == (subdomain_rings.size() - 1))
1133 : {
1134 32583 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID);
1135 32583 : if (generate_side_specific_boundaries)
1136 : {
1137 13775 : if (i <= angle_number / 2)
1138 2784 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID + side_index);
1139 : else
1140 10991 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID_ALT + side_index);
1141 : }
1142 : }
1143 29599 : else if (create_outward_interface_boundaries)
1144 17958 : boundary_info.add_side(elem, 2, k * 2 + 1 + boundary_id_shift);
1145 : }
1146 141010 : }
1147 95621 : j++;
1148 : }
1149 : }
1150 18971 : }
1151 :
1152 : std::unique_ptr<ReplicatedMesh>
1153 22326 : PolygonMeshGeneratorBase::buildSimplePeripheral(
1154 : const unsigned int num_sectors_per_side,
1155 : const unsigned int peripheral_invervals,
1156 : const std::vector<std::pair<Real, Real>> & positions_inner,
1157 : const std::vector<std::pair<Real, Real>> & d_positions_outer,
1158 : const subdomain_id_type id_shift,
1159 : const QUAD_ELEM_TYPE quad_elem_type,
1160 : const bool create_inward_interface_boundaries,
1161 : const bool create_outward_interface_boundaries)
1162 : {
1163 22326 : auto mesh = buildReplicatedMesh(2);
1164 : std::pair<Real, Real> positions_p;
1165 :
1166 : // generate node positions
1167 79454 : for (unsigned int i = 0; i <= peripheral_invervals; i++)
1168 : {
1169 177992 : for (unsigned int j = 0; j <= num_sectors_per_side / 2; j++)
1170 : {
1171 120864 : positions_p = pointInterpolate(positions_inner[0].first,
1172 120864 : positions_inner[0].second,
1173 120864 : d_positions_outer[0].first,
1174 120864 : d_positions_outer[0].second,
1175 120864 : positions_inner[1].first,
1176 120864 : positions_inner[1].second,
1177 120864 : d_positions_outer[1].first,
1178 120864 : d_positions_outer[1].second,
1179 : i,
1180 : j,
1181 : num_sectors_per_side,
1182 : peripheral_invervals);
1183 120864 : mesh->add_point(Point(positions_p.first, positions_p.second, 0.0));
1184 : }
1185 120864 : for (unsigned int j = 1; j <= num_sectors_per_side / 2; j++)
1186 : {
1187 63736 : positions_p = pointInterpolate(positions_inner[1].first,
1188 63736 : positions_inner[1].second,
1189 63736 : d_positions_outer[1].first,
1190 63736 : d_positions_outer[1].second,
1191 63736 : positions_inner[2].first,
1192 63736 : positions_inner[2].second,
1193 63736 : d_positions_outer[2].first,
1194 63736 : d_positions_outer[2].second,
1195 : i,
1196 : j,
1197 : num_sectors_per_side,
1198 : peripheral_invervals);
1199 63736 : mesh->add_point(Point(positions_p.first, positions_p.second, 0.0));
1200 : }
1201 : }
1202 :
1203 : // element definition
1204 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
1205 :
1206 57128 : for (unsigned int i = 0; i < peripheral_invervals; i++)
1207 : {
1208 113162 : for (unsigned int j = 0; j < num_sectors_per_side; j++)
1209 : {
1210 78360 : std::unique_ptr<Elem> new_elem;
1211 :
1212 78360 : new_elem = std::make_unique<Quad4>();
1213 78360 : new_elem->set_node(0, mesh->node_ptr(j + (num_sectors_per_side + 1) * (i)));
1214 78360 : new_elem->set_node(1, mesh->node_ptr(j + 1 + (num_sectors_per_side + 1) * (i)));
1215 78360 : new_elem->set_node(2, mesh->node_ptr(j + 1 + (num_sectors_per_side + 1) * (i + 1)));
1216 78360 : new_elem->set_node(3, mesh->node_ptr(j + (num_sectors_per_side + 1) * (i + 1)));
1217 :
1218 78360 : Elem * elem = mesh->add_elem(std::move(new_elem));
1219 :
1220 : // add subdoamin and boundary IDs
1221 78360 : elem->subdomain_id() = PERIPHERAL_ID_SHIFT + id_shift;
1222 78360 : if (i == 0)
1223 : {
1224 49112 : boundary_info.add_side(elem, 0, OUTER_SIDESET_ID);
1225 49112 : if (create_inward_interface_boundaries)
1226 0 : boundary_info.add_side(elem, 0, SLICE_ALT + id_shift * 2);
1227 : }
1228 78360 : if (i == peripheral_invervals - 1)
1229 : {
1230 49112 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID);
1231 49112 : if (create_outward_interface_boundaries)
1232 11776 : boundary_info.add_side(elem, 2, SLICE_ALT + id_shift * 2 + 1);
1233 : }
1234 78360 : if (j == 0)
1235 34802 : boundary_info.add_side(elem, 3, OUTER_SIDESET_ID);
1236 78360 : if (j == num_sectors_per_side - 1)
1237 34802 : boundary_info.add_side(elem, 1, OUTER_SIDESET_ID);
1238 78360 : }
1239 : }
1240 :
1241 : // convert element to second order if needed
1242 22326 : if (quad_elem_type != QUAD_ELEM_TYPE::QUAD4)
1243 : {
1244 : // full_ordered 2nd order element --> QUAD9, otherwise QUAD8
1245 1162 : const bool full_ordered = (quad_elem_type == QUAD_ELEM_TYPE::QUAD9);
1246 1162 : mesh->all_second_order(full_ordered);
1247 : }
1248 :
1249 22326 : return mesh;
1250 0 : }
1251 :
1252 : void
1253 28 : PolygonMeshGeneratorBase::adjustPeripheralQuadraticElements(
1254 : MeshBase & out_mesh, const QUAD_ELEM_TYPE boundary_quad_elem_type) const
1255 : {
1256 28 : const auto side_list = out_mesh.get_boundary_info().build_side_list();
1257 :
1258 : // select out elements on outer boundary
1259 : // std::set used to filter duplicate elem_ids
1260 : std::set<dof_id_type> elem_set;
1261 9436 : for (auto side_item : side_list)
1262 : {
1263 : boundary_id_type boundary_id = std::get<2>(side_item);
1264 9408 : dof_id_type elem_id = std::get<0>(side_item);
1265 :
1266 9408 : if (boundary_id == OUTER_SIDESET_ID)
1267 1456 : elem_set.insert(elem_id);
1268 : }
1269 :
1270 : // adjust nodes for outer boundary elements
1271 1484 : for (const auto elem_id : elem_set)
1272 : {
1273 1456 : Elem * elem = out_mesh.elem_ptr(elem_id);
1274 :
1275 : // adjust right side mid-edge node
1276 : Point pt_5 = (elem->point(1) + elem->point(2)) / 2.0;
1277 1456 : out_mesh.add_point(pt_5, elem->node_ptr(5)->id());
1278 :
1279 : // adjust left side mid-edge node
1280 : Point pt_7 = (elem->point(0) + elem->point(3)) / 2.0;
1281 1456 : out_mesh.add_point(pt_7, elem->node_ptr(7)->id());
1282 :
1283 : // adjust central node when using QUAD9
1284 1456 : if (boundary_quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
1285 : {
1286 728 : Point pt_8 = elem->true_centroid();
1287 728 : out_mesh.add_point(pt_8, elem->node_ptr(8)->id());
1288 : }
1289 : }
1290 28 : }
1291 :
1292 : std::pair<Real, Real>
1293 184600 : PolygonMeshGeneratorBase::pointInterpolate(const Real pi_1_x,
1294 : const Real pi_1_y,
1295 : const Real d_po_1_x,
1296 : const Real d_po_1_y,
1297 : const Real pi_2_x,
1298 : const Real pi_2_y,
1299 : const Real d_po_2_x,
1300 : const Real d_po_2_y,
1301 : const unsigned int i,
1302 : const unsigned int j,
1303 : const unsigned int num_sectors_per_side,
1304 : const unsigned int peripheral_intervals) const
1305 : {
1306 184600 : auto position_px_inner =
1307 184600 : (pi_1_x * (num_sectors_per_side / 2.0 - j) + pi_2_x * j) / (num_sectors_per_side / 2.0);
1308 184600 : auto position_py_inner =
1309 184600 : (pi_1_y * (num_sectors_per_side / 2.0 - j) + pi_2_y * j) / (num_sectors_per_side / 2.0);
1310 184600 : auto position_px_outer =
1311 184600 : (d_po_1_x * (num_sectors_per_side / 2.0 - j) + d_po_2_x * j) / (num_sectors_per_side / 2.0);
1312 184600 : auto position_py_outer =
1313 184600 : (d_po_1_y * (num_sectors_per_side / 2.0 - j) + d_po_2_y * j) / (num_sectors_per_side / 2.0);
1314 184600 : auto position_px = position_px_inner + position_px_outer * i / peripheral_intervals;
1315 184600 : auto position_py = position_py_inner + position_py_outer * i / peripheral_intervals;
1316 184600 : return std::make_pair(position_px, position_py);
1317 : }
1318 :
1319 : void
1320 390606 : PolygonMeshGeneratorBase::nodeCoordRotate(Real & x, Real & y, const Real theta) const
1321 : {
1322 390606 : const Real x_tmp = x;
1323 390606 : const Real y_tmp = y;
1324 390606 : x = x_tmp * std::cos(theta * M_PI / 180.0) - y_tmp * std::sin(theta * M_PI / 180.0);
1325 390606 : y = x_tmp * std::sin(theta * M_PI / 180.0) + y_tmp * std::cos(theta * M_PI / 180.0);
1326 390606 : }
1327 :
1328 : void
1329 1818 : PolygonMeshGeneratorBase::cutOffPolyDeform(MeshBase & mesh,
1330 : const Real orientation,
1331 : const Real y_max_0,
1332 : const Real y_max_n,
1333 : const Real y_min,
1334 : const unsigned int mesh_type,
1335 : const Real unit_angle,
1336 : const Real tols) const
1337 : {
1338 239150 : for (auto & node_ptr : as_range(mesh.nodes_begin(), mesh.nodes_end()))
1339 : {
1340 : // This function can definitely be optimized in future for better efficiency.
1341 116848 : Real & x = (*node_ptr)(0);
1342 : Real & y = (*node_ptr)(1);
1343 116848 : if (mesh_type == CORNER_MESH)
1344 : {
1345 74836 : nodeCoordRotate(x, y, orientation);
1346 74836 : if (x >= 0.0 && y > y_max_0)
1347 4292 : y = y - y_max_0 + y_max_n;
1348 70544 : else if (x >= 0.0 && y >= y_min)
1349 11139 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1350 59405 : else if (y > -x / std::tan(unit_angle / 360.0 * M_PI) + tols && y > y_max_0)
1351 : {
1352 266 : x /= y;
1353 266 : y = y - y_max_0 + y_max_n;
1354 266 : x *= y;
1355 : }
1356 59139 : else if (y > -x / std::tan(unit_angle / 360.0 * M_PI) + tols && y >= y_min)
1357 : {
1358 2082 : x /= y;
1359 2082 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1360 2082 : x *= y;
1361 : }
1362 74836 : nodeCoordRotate(x, y, -orientation);
1363 :
1364 74836 : nodeCoordRotate(x, y, orientation - unit_angle);
1365 74836 : if (x <= 0 && y > y_max_0)
1366 4374 : y = y - y_max_0 + y_max_n;
1367 70462 : else if (x <= 0 && y >= y_min)
1368 10305 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1369 60157 : else if (y >= x / std::tan(unit_angle / 360.0 * M_PI) - tols && y > y_max_0)
1370 : {
1371 2094 : x /= y;
1372 2094 : y = y - y_max_0 + y_max_n;
1373 2094 : x *= y;
1374 : }
1375 58063 : else if (y >= x / std::tan(unit_angle / 360.0 * M_PI) - tols && y >= y_min)
1376 : {
1377 6452 : x /= y;
1378 6452 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1379 6452 : x *= y;
1380 : }
1381 74836 : nodeCoordRotate(x, y, unit_angle - orientation);
1382 : }
1383 : else
1384 : {
1385 42012 : nodeCoordRotate(x, y, orientation);
1386 42012 : if (y > y_max_0)
1387 9852 : y = y - y_max_0 + y_max_n;
1388 32160 : else if (y >= y_min)
1389 14881 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1390 42012 : nodeCoordRotate(x, y, -orientation);
1391 : }
1392 1818 : }
1393 1818 : }
1394 :
1395 : std::pair<Real, Real>
1396 11849 : PolygonMeshGeneratorBase::fourPointIntercept(const std::pair<Real, Real> & p1,
1397 : const std::pair<Real, Real> & p2,
1398 : const std::pair<Real, Real> & p3,
1399 : const std::pair<Real, Real> & p4) const
1400 : {
1401 11849 : const Real x1 = p1.first;
1402 11849 : const Real y1 = p1.second;
1403 11849 : const Real x2 = p2.first;
1404 11849 : const Real y2 = p2.second;
1405 11849 : const Real x3 = p3.first;
1406 11849 : const Real y3 = p3.second;
1407 11849 : const Real x4 = p4.first;
1408 11849 : const Real y4 = p4.second;
1409 :
1410 11849 : Real x = -((x1 - x2) * (y3 * x4 - x3 * y4) - (x3 - x4) * (y1 * x2 - x1 * y2)) /
1411 11849 : ((y1 - y2) * (x3 - x4) - (y3 - y4) * (x1 - x2));
1412 11849 : Real y = -((y1 - y2) * (y3 * x4 - x3 * y4) - (y3 - y4) * (y1 * x2 - x1 * y2)) /
1413 : ((y1 - y2) * (x3 - x4) - (y3 - y4) * (x1 - x2));
1414 :
1415 11849 : return std::make_pair(x, y);
1416 : }
1417 :
1418 : std::vector<Real>
1419 1320 : PolygonMeshGeneratorBase::azimuthalAnglesCollector(ReplicatedMesh & mesh,
1420 : std::vector<Point> & boundary_points,
1421 : const Real lower_azi,
1422 : const Real upper_azi,
1423 : const unsigned int return_type,
1424 : const unsigned int num_sides,
1425 : const boundary_id_type bid,
1426 : const bool calculate_origin,
1427 : const Real input_origin_x,
1428 : const Real input_origin_y,
1429 : const Real tol) const
1430 : {
1431 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> side_list =
1432 1320 : mesh.get_boundary_info().build_side_list();
1433 2640 : mesh.get_boundary_info().build_node_list_from_side_list();
1434 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
1435 1320 : mesh.get_boundary_info().build_node_list();
1436 :
1437 : std::vector<Real> bd_x_list;
1438 : std::vector<Real> bd_y_list;
1439 : std::vector<Point> bd_p_list;
1440 : Real origin_x = 0.0;
1441 : Real origin_y = 0.0;
1442 : Real tmp_azi;
1443 1320 : const Real mid_azi = lower_azi <= upper_azi ? (lower_azi + upper_azi) / 2.0
1444 34 : : (lower_azi + upper_azi + 360.0) / 2.0;
1445 242138 : for (unsigned int i = 0; i < node_list.size(); ++i)
1446 240818 : if (std::get<1>(node_list[i]) == bid)
1447 : {
1448 49099 : bd_x_list.push_back((mesh.node_ref(std::get<0>(node_list[i])))(0));
1449 49099 : bd_y_list.push_back((mesh.node_ref(std::get<0>(node_list[i])))(1));
1450 49099 : bd_p_list.push_back((mesh.node_ref(std::get<0>(node_list[i]))));
1451 : }
1452 :
1453 1320 : if (calculate_origin)
1454 : {
1455 1320 : const Point origin_pt = MooseMeshUtils::meshCentroidCalculator(mesh);
1456 1320 : origin_x = origin_pt(0);
1457 1320 : origin_y = origin_pt(1);
1458 : }
1459 : else
1460 : {
1461 : origin_x = input_origin_x;
1462 : origin_y = input_origin_y;
1463 : }
1464 :
1465 : std::vector<std::pair<Real, Point>> azi_point_pairs;
1466 :
1467 50419 : for (unsigned int i = 0; i < bd_x_list.size(); ++i)
1468 : {
1469 49099 : tmp_azi = atan2(bd_y_list[i] - origin_y, bd_x_list[i] - origin_x) * 180.0 / M_PI;
1470 49099 : if ((lower_azi <= upper_azi && (tmp_azi >= lower_azi - tol && tmp_azi <= upper_azi + tol)) ||
1471 1784 : (lower_azi > upper_azi && (tmp_azi >= lower_azi - tol || tmp_azi <= upper_azi + tol)))
1472 : {
1473 : azi_point_pairs.push_back(
1474 15396 : std::make_pair(return_type == ANGLE_DEGREE
1475 15396 : ? (tmp_azi - mid_azi)
1476 2527 : : (1.0 + std::cos(M_PI / num_sides) / std::sin(M_PI / num_sides) *
1477 2527 : std::tan((tmp_azi - mid_azi) / 180.0 * M_PI)),
1478 : bd_p_list[i]));
1479 : }
1480 : }
1481 1320 : std::sort(azi_point_pairs.begin(), azi_point_pairs.end());
1482 :
1483 : std::vector<Real> azimuthal_output;
1484 : for (auto it = std::make_move_iterator(azi_point_pairs.begin()),
1485 : end = std::make_move_iterator(azi_point_pairs.end());
1486 16716 : it != end;
1487 : it++)
1488 : {
1489 15396 : azimuthal_output.push_back(std::move(it->first));
1490 15396 : boundary_points.push_back(std::move(it->second));
1491 : }
1492 :
1493 1320 : return azimuthal_output;
1494 1320 : }
1495 :
1496 : std::vector<Real>
1497 668 : PolygonMeshGeneratorBase::azimuthalAnglesCollector(ReplicatedMesh & mesh,
1498 : const Real lower_azi,
1499 : const Real upper_azi,
1500 : const unsigned int return_type,
1501 : const unsigned int num_sides,
1502 : const boundary_id_type bid,
1503 : const bool calculate_origin,
1504 : const Real input_origin_x,
1505 : const Real input_origin_y,
1506 : const Real tol) const
1507 : {
1508 : std::vector<Point> boundary_points;
1509 : return azimuthalAnglesCollector(mesh,
1510 : boundary_points,
1511 : lower_azi,
1512 : upper_azi,
1513 : return_type,
1514 : num_sides,
1515 : bid,
1516 : calculate_origin,
1517 : input_origin_x,
1518 : input_origin_y,
1519 1336 : tol);
1520 668 : }
1521 :
1522 : std::vector<std::vector<Real>>
1523 75884 : PolygonMeshGeneratorBase::biasTermsCalculator(
1524 : const std::vector<Real> radial_biases,
1525 : const std::vector<unsigned int> intervals,
1526 : const multiBdryLayerParams inner_boundary_layer_params,
1527 : const multiBdryLayerParams outer_boundary_layer_params) const
1528 : {
1529 : std::vector<std::vector<Real>> bias_terms_vec;
1530 133406 : for (unsigned int i = 0; i < radial_biases.size(); i++)
1531 115044 : bias_terms_vec.push_back(biasTermsCalculator(radial_biases[i],
1532 : intervals[i],
1533 : {0.0,
1534 : inner_boundary_layer_params.fractions[i],
1535 : inner_boundary_layer_params.intervals[i],
1536 : inner_boundary_layer_params.biases[i]},
1537 : {0.0,
1538 : outer_boundary_layer_params.fractions[i],
1539 : outer_boundary_layer_params.intervals[i],
1540 : outer_boundary_layer_params.biases[i]}));
1541 75884 : return bias_terms_vec;
1542 0 : }
1543 :
1544 : std::vector<Real>
1545 171639 : PolygonMeshGeneratorBase::biasTermsCalculator(
1546 : const Real radial_bias,
1547 : const unsigned int intervals,
1548 : const singleBdryLayerParams inner_boundary_layer_params,
1549 : const singleBdryLayerParams outer_boundary_layer_params) const
1550 : {
1551 : // To get biased indices:
1552 : // If no bias is involved, namely bias factor = 1.0, the increment in indices is uniform.
1553 : // Thus, (i + 1) is used to get such linearly increasing indices.
1554 : // If a non-trivial bias factor q is used, the increment in the indices is geometric
1555 : // progression. So, if first (i = 0) increment is 1.0, second (i = 1) is q, third (i = 2) is
1556 : // q^2,..., last or n_interval'th is q^(n_interval - 1). Then, the summation of the first (i +
1557 : // 1) increments over the summation of all n_interval increments is the (i + 1)th index The
1558 : // summation of the first (i + 1) increments is (1.0 - q^(i + 1)) / (1 - q); The summation of
1559 : // all n_interval increments is (1.0 - q^n_interval) / (1 - q); Thus, the index is (1.0 - q^(i +
1560 : // 1)) / (1.0 - q^n_interval)
1561 : // This approach is used by inner boundary layer, main region, outer boundary layer separately.
1562 :
1563 : std::vector<Real> biased_terms;
1564 173214 : for (unsigned int i = 0; i < inner_boundary_layer_params.intervals; i++)
1565 : biased_terms.push_back(
1566 : MooseUtils::absoluteFuzzyEqual(inner_boundary_layer_params.bias, 1.0)
1567 3150 : ? ((Real)(i + 1) * inner_boundary_layer_params.fraction /
1568 0 : (Real)inner_boundary_layer_params.intervals)
1569 1575 : : ((1.0 - std::pow(inner_boundary_layer_params.bias, (Real)(i + 1))) /
1570 1575 : (1.0 - std::pow(inner_boundary_layer_params.bias,
1571 : (Real)(inner_boundary_layer_params.intervals))) *
1572 : inner_boundary_layer_params.fraction));
1573 400933 : for (unsigned int i = 0; i < intervals; i++)
1574 229294 : biased_terms.push_back(inner_boundary_layer_params.fraction +
1575 : (MooseUtils::absoluteFuzzyEqual(radial_bias, 1.0)
1576 229294 : ? ((Real)(i + 1) *
1577 227109 : (1.0 - inner_boundary_layer_params.fraction -
1578 : outer_boundary_layer_params.fraction) /
1579 227109 : (Real)intervals)
1580 2185 : : ((1.0 - std::pow(radial_bias, (Real)(i + 1))) /
1581 2185 : (1.0 - std::pow(radial_bias, (Real)(intervals))) *
1582 2185 : (1.0 - inner_boundary_layer_params.fraction -
1583 : outer_boundary_layer_params.fraction))));
1584 173739 : for (unsigned int i = 0; i < outer_boundary_layer_params.intervals; i++)
1585 : biased_terms.push_back(
1586 4200 : 1.0 - outer_boundary_layer_params.fraction +
1587 : (MooseUtils::absoluteFuzzyEqual(outer_boundary_layer_params.bias, 1.0)
1588 2100 : ? ((Real)(i + 1) * outer_boundary_layer_params.fraction /
1589 0 : (Real)outer_boundary_layer_params.intervals)
1590 2100 : : ((1.0 - std::pow(outer_boundary_layer_params.bias, (Real)(i + 1))) /
1591 2100 : (1.0 - std::pow(outer_boundary_layer_params.bias,
1592 : (Real)(outer_boundary_layer_params.intervals))) *
1593 : outer_boundary_layer_params.fraction)));
1594 171639 : return biased_terms;
1595 0 : }
1596 :
1597 : void
1598 4302 : PolygonMeshGeneratorBase::addRingAndSectorIDParams(InputParameters & params)
1599 : {
1600 8604 : params.addParam<std::string>("sector_id_name",
1601 : "Name of integer (reporting) ID for sector regions to use the "
1602 : "reporting ID for azimuthal sector regions of ring geometry block.");
1603 8604 : params.addParam<std::string>("ring_id_name",
1604 : "Name of integer (reporting) ID for ring regions to use the "
1605 : "reporting ID for annular regions of ring geometry block.");
1606 8604 : MooseEnum ring_id_option("block_wise ring_wise", "block_wise");
1607 8604 : params.addParam<MooseEnum>(
1608 : "ring_id_assign_type", ring_id_option, "Type of ring ID assignment: block_wise or ring_wise");
1609 8604 : params.addParamNamesToGroup("sector_id_name ring_id_name ring_id_assign_type", "Ring/Sector IDs");
1610 4302 : }
1611 :
1612 : void
1613 56 : PolygonMeshGeneratorBase::setSectorExtraIDs(MeshBase & mesh,
1614 : const std::string id_name,
1615 : const unsigned int num_sides,
1616 : const std::vector<unsigned int> num_sectors_per_side)
1617 : {
1618 56 : const auto extra_id_index = mesh.add_elem_integer(id_name);
1619 : // vector to store sector ids for each element
1620 56 : auto elem_it = mesh.elements_begin();
1621 : unsigned int id = 1;
1622 : // starting element id of the current sector
1623 301 : for (unsigned int is = 0; is < num_sides; ++is)
1624 : {
1625 : // number of elements in the current sector
1626 : unsigned int nelem_sector =
1627 245 : mesh.n_elem() * num_sectors_per_side[is] /
1628 245 : (accumulate(num_sectors_per_side.begin(), num_sectors_per_side.end(), 0));
1629 : // assign sector ids to mesh
1630 11725 : for (unsigned i = 0; i < nelem_sector; ++i, ++elem_it)
1631 5740 : (*elem_it)->set_extra_integer(extra_id_index, id);
1632 : // update sector id
1633 245 : ++id;
1634 : }
1635 56 : }
1636 :
1637 : void
1638 56 : PolygonMeshGeneratorBase::setRingExtraIDs(MeshBase & mesh,
1639 : const std::string id_name,
1640 : const unsigned int num_sides,
1641 : const std::vector<unsigned int> num_sectors_per_side,
1642 : const std::vector<unsigned int> ring_intervals,
1643 : const bool ring_wise_id,
1644 : const bool quad_center_elements)
1645 : {
1646 : // this function assumes that elements are ordered by rings (inner) then by sectors (outer
1647 : // ordering)
1648 56 : const auto extra_id_index = mesh.add_elem_integer(id_name);
1649 56 : auto elem_it = mesh.elements_begin();
1650 301 : for (unsigned int is = 0; is < num_sides; ++is)
1651 : {
1652 : // number of elements in the current sector
1653 245 : unsigned int nelem = mesh.n_elem() * num_sectors_per_side[is] /
1654 245 : (accumulate(num_sectors_per_side.begin(), num_sectors_per_side.end(), 0));
1655 245 : if (!ring_wise_id)
1656 : {
1657 777 : for (unsigned int ir : index_range(ring_intervals))
1658 : {
1659 : // number of elements in the current ring and sector
1660 560 : unsigned int nelem_annular_ring = num_sectors_per_side[is] * ring_intervals[ir];
1661 : // if _quad_center_elements is true, the number of elements in center ring are
1662 : // _num_sectors_per_side[is] * _num_sectors_per_side[is] / 4
1663 560 : if (quad_center_elements && ir == 0)
1664 0 : nelem_annular_ring = num_sectors_per_side[is] * (ring_intervals[ir] - 1) +
1665 0 : num_sectors_per_side[is] * num_sectors_per_side[is] / 4;
1666 : // assign ring id
1667 7784 : for (unsigned i = 0; i < nelem_annular_ring; ++i, ++elem_it)
1668 3612 : (*elem_it)->set_extra_integer(extra_id_index, ir + 1);
1669 : // update number of elements in background region of current side.
1670 560 : nelem -= nelem_annular_ring;
1671 : }
1672 : }
1673 : else
1674 : {
1675 : unsigned int ir = 0;
1676 112 : for (unsigned int ir0 : index_range(ring_intervals))
1677 : {
1678 224 : for (unsigned int ir1 = 0; ir1 < ring_intervals[ir0]; ++ir1)
1679 : {
1680 : // number of elements in the current ring and sector
1681 140 : unsigned int nelem_annular_ring = num_sectors_per_side[is];
1682 : // if _quad_center_elements is true, the number of elements in center ring are
1683 : // _num_sectors_per_side[is] * _num_sectors_per_side[is] / 4
1684 140 : if (quad_center_elements && ir == 0)
1685 0 : nelem_annular_ring = num_sectors_per_side[is] * num_sectors_per_side[is] / 4;
1686 : // assign ring id
1687 1260 : for (unsigned i = 0; i < nelem_annular_ring; ++i, ++elem_it)
1688 560 : (*elem_it)->set_extra_integer(extra_id_index, ir + 1);
1689 : // update ring id
1690 140 : ++ir;
1691 : // update number of elements in background region of current side.
1692 140 : nelem -= nelem_annular_ring;
1693 : }
1694 : }
1695 : }
1696 : // assign ring id of 0 to the background region
1697 4277 : for (unsigned i = 0; i < nelem; ++i, ++elem_it)
1698 2016 : (*elem_it)->set_extra_integer(extra_id_index, 0);
1699 : }
1700 56 : }
1701 :
1702 : void
1703 161 : PolygonMeshGeneratorBase::reassignBoundaryIDs(MeshBase & mesh,
1704 : const boundary_id_type id_shift,
1705 : const std::set<boundary_id_type> & boundary_ids,
1706 : const bool reverse)
1707 : {
1708 : const std::set<boundary_id_type> existing_boundary_ids =
1709 : mesh.get_boundary_info().get_boundary_ids();
1710 413 : for (const auto id : boundary_ids)
1711 : {
1712 :
1713 252 : const boundary_id_type old_id = (!reverse) ? id : id + id_shift;
1714 252 : const boundary_id_type new_id = (!reverse) ? id + id_shift : id;
1715 : auto it = existing_boundary_ids.find(old_id);
1716 252 : if (it != existing_boundary_ids.end())
1717 252 : MooseMesh::changeBoundaryId(mesh, old_id, new_id, true);
1718 : }
1719 161 : }
1720 :
1721 : std::set<boundary_id_type>
1722 1276 : PolygonMeshGeneratorBase::getInterfaceBoundaryIDs(
1723 : const std::vector<std::vector<unsigned int>> & pattern,
1724 : const std::vector<std::vector<boundary_id_type>> & interface_boundary_id_shift_pattern,
1725 : const std::set<boundary_id_type> & boundary_ids,
1726 : const std::vector<std::set<boundary_id_type>> & input_interface_boundary_ids,
1727 : const bool use_interface_boundary_id_shift,
1728 : const bool create_interface_boundary_id,
1729 : const unsigned int num_extra_layers) const
1730 : {
1731 : std::set<boundary_id_type> interface_boundary_ids;
1732 : // add existing interface boundary ids from input meshes
1733 1276 : if (use_interface_boundary_id_shift)
1734 : {
1735 56 : for (const auto i : make_range(pattern.size()))
1736 154 : for (const auto j : make_range(pattern[i].size()))
1737 : {
1738 112 : const auto & ids = input_interface_boundary_ids[pattern[i][j]];
1739 273 : for (const auto & id : ids)
1740 : {
1741 161 : const boundary_id_type new_id = id + interface_boundary_id_shift_pattern[i][j];
1742 : auto it = boundary_ids.find(new_id);
1743 161 : if (it != boundary_ids.end())
1744 161 : interface_boundary_ids.insert(new_id);
1745 : }
1746 : }
1747 : }
1748 : else
1749 : {
1750 3582 : for (const auto & ids : input_interface_boundary_ids)
1751 3790 : for (const auto & id : ids)
1752 : {
1753 : auto it = boundary_ids.find(id);
1754 1470 : if (it != boundary_ids.end())
1755 1421 : interface_boundary_ids.insert(id);
1756 : }
1757 : }
1758 : // add unshifted interface boundary ids for the duct & background regions
1759 1276 : if (create_interface_boundary_id)
1760 1504 : for (const auto i : make_range(num_extra_layers))
1761 : {
1762 742 : boundary_id_type id = SLICE_ALT + i * 2 + 1;
1763 : auto it = boundary_ids.find(id);
1764 742 : if (it != boundary_ids.end())
1765 249 : interface_boundary_ids.insert(id);
1766 742 : id = SLICE_ALT + i * 2;
1767 : it = boundary_ids.find(id);
1768 742 : if (it != boundary_ids.end())
1769 0 : interface_boundary_ids.insert(id);
1770 : }
1771 1276 : return interface_boundary_ids;
1772 : }
1773 :
1774 : PolygonMeshGeneratorBase::multiBdryLayerParams
1775 75884 : PolygonMeshGeneratorBase::modifiedMultiBdryLayerParamsCreator(
1776 : const multiBdryLayerParams & original_multi_bdry_layer_params, const unsigned int order) const
1777 : {
1778 75884 : multiBdryLayerParams mod_multi_bdry_layer_params(original_multi_bdry_layer_params);
1779 : std::for_each(mod_multi_bdry_layer_params.intervals.begin(),
1780 : mod_multi_bdry_layer_params.intervals.end(),
1781 57522 : [&order](unsigned int & n) { n *= order; });
1782 75884 : std::for_each(mod_multi_bdry_layer_params.biases.begin(),
1783 : mod_multi_bdry_layer_params.biases.end(),
1784 57522 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
1785 75884 : return mod_multi_bdry_layer_params;
1786 : }
1787 :
1788 : PolygonMeshGeneratorBase::singleBdryLayerParams
1789 37942 : PolygonMeshGeneratorBase::modifiedSingleBdryLayerParamsCreator(
1790 : const singleBdryLayerParams & original_single_bdry_layer_params, const unsigned int order) const
1791 : {
1792 37942 : singleBdryLayerParams mod_single_bdry_layer_params(original_single_bdry_layer_params);
1793 37942 : mod_single_bdry_layer_params.intervals *= order;
1794 37942 : mod_single_bdry_layer_params.bias = std::pow(mod_single_bdry_layer_params.bias, 1.0 / order);
1795 37942 : return mod_single_bdry_layer_params;
1796 : }
1797 :
1798 : std::string
1799 4 : PolygonMeshGeneratorBase::pitchMetaDataErrorGenerator(
1800 : const std::vector<MeshGeneratorName> & input_names,
1801 : const std::vector<Real> & metadata_vals,
1802 : const std::string & metadata_name) const
1803 : {
1804 4 : FormattedTable table;
1805 12 : for (unsigned int i = 0; i < input_names.size(); i++)
1806 : {
1807 8 : table.addRow(i);
1808 16 : table.addData<std::string>("input name", (std::string)input_names[i]);
1809 8 : table.addData<Real>(metadata_name, metadata_vals[i]);
1810 : }
1811 : table.outputTimeColumn(false);
1812 4 : std::stringstream detailed_error;
1813 4 : table.printTable(detailed_error);
1814 8 : return "\n" + detailed_error.str();
1815 4 : }
|
