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 12574 : PolygonMeshGeneratorBase::validParams()
21 : {
22 12574 : InputParameters params = MeshGenerator::validParams();
23 12574 : params.addClassDescription(
24 : "A base class that contains common members for Reactor module mesh generators.");
25 :
26 12574 : return params;
27 0 : }
28 :
29 6397 : PolygonMeshGeneratorBase::PolygonMeshGeneratorBase(const InputParameters & parameters)
30 6397 : : MeshGenerator(parameters)
31 : {
32 6397 : }
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 270 : 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 270 : const Real virtual_pitch = 2.0 * primary_side_length * cos(azimuthal_angle / 360.0 * M_PI);
70 270 : const Real virtual_side_number = 360.0 / azimuthal_angle;
71 270 : 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 270 : generate_side_specific_boundaries);
101 270 : MeshTools::Modification::rotate(*mesh, rotation_angle, 0, 0);
102 270 : return mesh;
103 0 : }
104 :
105 : std::unique_ptr<ReplicatedMesh>
106 12157 : 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 12157 : quad_elem_type);
167 : }
168 :
169 : std::unique_ptr<ReplicatedMesh>
170 21861 : 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 21861 : const unsigned short order = quad_elem_type == QUAD_ELEM_TYPE::QUAD4 ? 1 : 2;
203 23040 : 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 21861 : std::vector<unsigned int> mod_ring_layers(ring_layers);
214 : std::for_each(
215 30197 : 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 21861 : std::vector<Real> mod_ring_radial_biases(ring_radial_biases);
219 21861 : std::for_each(mod_ring_radial_biases.begin(),
220 : mod_ring_radial_biases.end(),
221 30197 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
222 : // ducts_layers is similar to ring_layers
223 21861 : std::vector<unsigned int> mod_ducts_layers(ducts_layers);
224 : std::for_each(
225 2226 : 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 21861 : std::vector<Real> mod_duct_radial_biases(duct_radial_biases);
228 21861 : std::for_each(mod_duct_radial_biases.begin(),
229 : mod_duct_radial_biases.end(),
230 2226 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
231 : // Azimuthal mesh density is also doubled for order = 2
232 21861 : const unsigned int mod_num_sectors_per_side = num_sectors_per_side * order;
233 21861 : const unsigned int mod_background_intervals = background_intervals * order;
234 : // background_radial_bias is similar to ring_radial_biases
235 21861 : 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 21861 : modifiedMultiBdryLayerParamsCreator(ring_inner_boundary_layer_params, order);
239 : const auto mod_ring_outer_boundary_layer_params =
240 21861 : modifiedMultiBdryLayerParamsCreator(ring_outer_boundary_layer_params, order);
241 : const auto mod_duct_inner_boundary_layer_params =
242 21861 : modifiedMultiBdryLayerParamsCreator(duct_inner_boundary_layer_params, order);
243 : const auto mod_duct_outer_boundary_layer_params =
244 21861 : modifiedMultiBdryLayerParamsCreator(duct_outer_boundary_layer_params, order);
245 :
246 : const auto mod_background_inner_boundary_layer_params =
247 21861 : modifiedSingleBdryLayerParamsCreator(background_inner_boundary_layer_params, order);
248 : const auto mod_background_outer_boundary_layer_params =
249 21861 : 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 21861 : std::vector<Real> mod_ducts_center_dist(ducts_center_dist);
254 21861 : std::vector<Real> mod_ring_radii(ring_radii);
255 21861 : bool has_rings(ring_radii.size());
256 21861 : bool has_ducts(ducts_center_dist.size());
257 : bool has_background(background_intervals);
258 21861 : 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 21861 : biasTermsCalculator(background_radial_bias, background_intervals);
264 : const auto inner_background_bias_terms =
265 21861 : biasTermsCalculator(background_inner_boundary_layer_params.bias,
266 21861 : background_inner_boundary_layer_params.intervals);
267 : const auto outer_background_bias_terms =
268 21861 : biasTermsCalculator(background_outer_boundary_layer_params.bias,
269 21861 : 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 21861 : ring_outer_boundary_layer_params);
274 : auto duct_bias_terms = biasTermsCalculator(duct_radial_biases,
275 : ducts_layers,
276 : duct_inner_boundary_layer_params,
277 21861 : duct_outer_boundary_layer_params);
278 : // Equivalent "mod_" parts
279 : const auto mod_main_background_bias_terms =
280 21861 : biasTermsCalculator(mod_background_radial_bias, mod_background_intervals);
281 : const auto mod_inner_background_bias_terms =
282 21861 : biasTermsCalculator(mod_background_inner_boundary_layer_params.bias,
283 21861 : mod_background_inner_boundary_layer_params.intervals);
284 : const auto mod_outer_background_bias_terms =
285 21861 : biasTermsCalculator(mod_background_outer_boundary_layer_params.bias,
286 21861 : 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 21861 : 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 21861 : mod_duct_outer_boundary_layer_params);
295 :
296 : std::vector<unsigned int> total_ring_layers;
297 52058 : for (unsigned int i = 0; i < ring_layers.size(); i++)
298 30197 : 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 21861 : if (background_inner_boundary_layer_params.intervals)
302 : {
303 135 : total_ring_layers.push_back(background_inner_boundary_layer_params.intervals);
304 135 : rings_bias_terms.push_back(inner_background_bias_terms);
305 270 : ring_radii.push_back((ring_radii.empty() ? 0.0 : ring_radii.back()) +
306 135 : background_inner_boundary_layer_params.width);
307 : has_rings = true;
308 : }
309 : std::vector<unsigned int> mod_total_ring_layers;
310 52058 : for (unsigned int i = 0; i < mod_ring_layers.size(); i++)
311 30197 : mod_total_ring_layers.push_back(mod_ring_layers[i] +
312 30197 : mod_ring_inner_boundary_layer_params.intervals[i] +
313 : mod_ring_outer_boundary_layer_params.intervals[i]);
314 :
315 21861 : if (mod_background_inner_boundary_layer_params.intervals)
316 : {
317 135 : mod_total_ring_layers.push_back(mod_background_inner_boundary_layer_params.intervals);
318 135 : mod_rings_bias_terms.push_back(mod_inner_background_bias_terms);
319 270 : mod_ring_radii.push_back((mod_ring_radii.empty() ? 0.0 : mod_ring_radii.back()) +
320 135 : 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 21861 : if (background_outer_boundary_layer_params.intervals)
326 : {
327 135 : total_ducts_layers.push_back(background_outer_boundary_layer_params.intervals);
328 135 : duct_bias_terms.insert(duct_bias_terms.begin(), outer_background_bias_terms);
329 135 : ducts_center_dist.insert(ducts_center_dist.begin(),
330 : (ducts_center_dist.empty()
331 135 : ? pitch / 2.0 / std::cos(M_PI / virtual_side_number)
332 135 : : ducts_center_dist.front()) -
333 135 : background_outer_boundary_layer_params.width);
334 : has_ducts = true;
335 : }
336 24087 : for (unsigned int i = 0; i < ducts_layers.size(); i++)
337 2226 : 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 21861 : if (mod_background_outer_boundary_layer_params.intervals)
342 : {
343 135 : mod_total_ducts_layers.push_back(mod_background_outer_boundary_layer_params.intervals);
344 135 : mod_duct_bias_terms.insert(mod_duct_bias_terms.begin(), mod_outer_background_bias_terms);
345 135 : mod_ducts_center_dist.insert(mod_ducts_center_dist.begin(),
346 : (mod_ducts_center_dist.empty()
347 135 : ? pitch / 2.0 / std::cos(M_PI / virtual_side_number)
348 135 : : mod_ducts_center_dist.front()) -
349 135 : mod_background_outer_boundary_layer_params.width);
350 : // has_ducts should be modified before in the none "mod_" part
351 : }
352 24087 : for (unsigned int i = 0; i < mod_ducts_layers.size(); i++)
353 2226 : mod_total_ducts_layers.push_back(mod_ducts_layers[i] +
354 2226 : 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 21861 : : ((azimuthal_tangent.size() - 1) / order);
360 : unsigned int mod_angle_number =
361 21861 : azimuthal_tangent.size() == 0 ? mod_num_sectors_per_side : (azimuthal_tangent.size() - 1);
362 :
363 : // Geometries
364 21861 : const Real corner_to_corner =
365 21861 : pitch / std::cos(M_PI / virtual_side_number); // distance of bin center to cell corner
366 21861 : const Real corner_p[2][2] = {
367 21861 : {0.0, 0.5 * corner_to_corner},
368 21861 : {0.5 * corner_to_corner * pitch_scale_factor * std::sin(2.0 * M_PI / virtual_side_number),
369 21861 : 0.5 * corner_to_corner * pitch_scale_factor * std::cos(2.0 * M_PI / virtual_side_number)}};
370 21861 : const unsigned int div_num = angle_number / 2 + 1;
371 21861 : 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 21861 : std::vector<std::vector<Node *>> nodes(mod_div_num, std::vector<Node *>(mod_div_num));
375 21861 : if (quad_center_elements)
376 : {
377 : Real ring_radii_0;
378 :
379 2863 : if (has_rings)
380 1843 : ring_radii_0 = ring_radii.front() * mod_rings_bias_terms.front()[order - 1];
381 1020 : else if (has_ducts)
382 108 : ring_radii_0 = mod_ducts_center_dist.front() * std::cos(M_PI / virtual_side_number) *
383 108 : mod_main_background_bias_terms[order - 1];
384 : else
385 912 : 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 2863 : ring_radii_0 *=
390 2863 : center_quad_factor == 0.0 ? (((Real)div_num - 1.0) / (Real)div_num) : center_quad_factor;
391 :
392 2863 : centerNodes(*mesh, virtual_side_number, mod_div_num, ring_radii_0, nodes);
393 : }
394 : else // pin-cell center
395 18998 : mesh->add_point(Point(0.0, 0.0, 0.0));
396 :
397 : // create nodes for the ring regions
398 21861 : if (has_rings)
399 18783 : 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 21861 : 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 12427 : if (has_rings)
418 9349 : background_in = ring_radii.back();
419 : else
420 : background_in = 0;
421 :
422 12427 : if (has_ducts)
423 : {
424 1812 : 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 12427 : background_corner_radial_interval_length =
437 12427 : (background_out - background_in) / mod_background_intervals;
438 :
439 12427 : node_id_background_meta = mesh->n_nodes();
440 :
441 : // create nodes for background region
442 12427 : 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 21861 : if (has_ducts)
456 1812 : 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 21861 : if (ring_layers.empty())
478 3123 : is_central_region_independent = mod_background_inner_boundary_layer_params.intervals +
479 3123 : mod_background_intervals +
480 : mod_background_outer_boundary_layer_params.intervals ==
481 : 1;
482 : else
483 18738 : is_central_region_independent = mod_ring_layers[0] +
484 18738 : 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 21861 : if (quad_center_elements)
492 5726 : cenQuadElemDef(*mesh,
493 : div_num,
494 : block_id_shift,
495 2863 : create_outward_interface_boundaries && is_central_region_independent,
496 : boundary_id_shift,
497 : nodes,
498 2863 : (!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 18998 : cenTriElemDef(
507 : *mesh,
508 : num_sectors_per_side,
509 : azimuthal_tangent,
510 : block_id_shift,
511 18998 : create_outward_interface_boundaries && is_central_region_independent,
512 : boundary_id_shift,
513 18998 : ((!has_rings) && (!has_ducts) && (background_intervals == 1)) ||
514 9434 : ((!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 21861 : if (has_rings) // define the rings in each subdomain
529 : {
530 18783 : subdomain_rings = total_ring_layers;
531 18783 : subdomain_rings.front() -= 1; // remove the inner TRI mesh subdomain
532 18783 : if (background_inner_boundary_layer_params.intervals)
533 : {
534 135 : subdomain_rings.back() =
535 135 : background_inner_boundary_layer_params.intervals + background_intervals +
536 135 : background_outer_boundary_layer_params.intervals; // add the background region
537 135 : if (ring_radii.size() == 1)
538 45 : subdomain_rings.back() -= 1; // remove the inner TRI mesh subdomain
539 : }
540 18648 : else if (has_background)
541 9214 : subdomain_rings.push_back(background_inner_boundary_layer_params.intervals +
542 9214 : background_intervals +
543 9214 : background_outer_boundary_layer_params.intervals);
544 : }
545 : else
546 : {
547 : subdomain_rings.push_back(
548 3078 : background_inner_boundary_layer_params.intervals + background_intervals +
549 3078 : background_outer_boundary_layer_params.intervals); // add the background region
550 3078 : subdomain_rings[0] -= 1; // remove the inner TRI mesh subdomain
551 : }
552 :
553 21861 : if (has_ducts)
554 4038 : for (unsigned int i = (background_outer_boundary_layer_params.intervals > 0);
555 4038 : i < total_ducts_layers.size();
556 : i++)
557 2226 : subdomain_rings.push_back(total_ducts_layers[i]);
558 :
559 24724 : quadElemDef(*mesh,
560 : num_sectors_per_side,
561 : subdomain_rings,
562 : side_index,
563 : azimuthal_tangent,
564 : block_id_shift,
565 2863 : 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 21861 : if (tri_elem_type == TRI_ELEM_TYPE::TRI6 || quad_elem_type == QUAD_ELEM_TYPE::QUAD8)
572 733 : mesh->remove_orphaned_nodes();
573 21861 : return mesh;
574 21861 : }
575 :
576 : void
577 2863 : 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 2863 : std::make_pair(0.0, ring_radii_0 * std::cos(M_PI / virtual_side_number));
586 : const std::pair<Real, Real> p_top =
587 2863 : std::make_pair(p_bottom.second * std::sin(2.0 * M_PI / virtual_side_number),
588 2863 : p_bottom.second * std::cos(2.0 * M_PI / virtual_side_number));
589 : const std::pair<Real, Real> p_diag =
590 2863 : 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 9204 : for (unsigned int i = 0; i < div_num; i++)
628 : {
629 : unsigned int id_x = 0;
630 : unsigned int id_y = i;
631 21408 : for (unsigned int j = 0; j < 2 * i + 1; j++)
632 : {
633 15067 : std::pair<Real, Real> p1 = std::make_pair(
634 15067 : (p_origin.first * (div_num - 1 - id_x) + p_top.first * id_x) / (div_num - 1),
635 15067 : (p_origin.second * (div_num - 1 - id_x) + p_top.second * id_x) / (div_num - 1));
636 15067 : std::pair<Real, Real> p2 = std::make_pair(
637 15067 : (p_bottom.first * (div_num - 1 - id_x) + p_diag.first * id_x) / (div_num - 1),
638 15067 : (p_bottom.second * (div_num - 1 - id_x) + p_diag.second * id_x) / (div_num - 1));
639 15067 : std::pair<Real, Real> p3 = std::make_pair(
640 15067 : (p_origin.first * (div_num - 1 - id_y) + p_bottom.first * id_y) / (div_num - 1),
641 15067 : (p_origin.second * (div_num - 1 - id_y) + p_bottom.second * id_y) / (div_num - 1));
642 15067 : std::pair<Real, Real> p4 = std::make_pair(
643 15067 : (p_top.first * (div_num - 1 - id_y) + p_diag.first * id_y) / (div_num - 1),
644 15067 : (p_top.second * (div_num - 1 - id_y) + p_diag.second * id_y) / (div_num - 1));
645 15067 : std::pair<Real, Real> pc = fourPointIntercept(p1, p2, p3, p4);
646 15067 : nodes[id_x][id_y] = mesh.add_point(Point(pc.first, pc.second, 0.0));
647 15067 : if (j < i)
648 4363 : id_x++;
649 15067 : if (j >= i)
650 10704 : id_y--;
651 : }
652 : }
653 2863 : }
654 :
655 : void
656 18783 : 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 18783 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
667 :
668 : // Add nodes in pins regions
669 49115 : 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 30332 : l == 0 ? ring_radii[l] / ring_layers[l]
674 11549 : : (ring_radii[l] - ring_radii[l - 1]) / ring_layers[l];
675 :
676 : // add rings in each pin subdomain
677 135102 : for (unsigned int k = 0; k < ring_layers[l]; k++)
678 : {
679 : const Real bin_radial_distance =
680 104770 : 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 15248 : : (ring_radii[l - 1] +
684 15248 : biased_terms[l][k] * ring_layers[l] * pin_radius_interval_length);
685 104770 : const Real pin_corner_p_x = corner_p[0][0] * bin_radial_distance / (0.5 * corner_to_corner);
686 104770 : 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 104770 : mesh.add_point(Point(pin_corner_p_x, pin_corner_p_y, 0.0));
691 :
692 260920 : for (unsigned int j = 1; j <= angle_number; j++)
693 : {
694 : const Real cell_boundary_p_x =
695 156150 : corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
696 156150 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
697 1728 : : (azimuthal_tangent[j] / 2.0));
698 : const Real cell_boundary_p_y =
699 156150 : 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 156150 : cell_boundary_p_x * bin_radial_distance /
706 156150 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
707 : const Real pin_azimuthal_p_y =
708 156150 : cell_boundary_p_y * bin_radial_distance /
709 156150 : 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 156150 : mesh.add_point(Point(pin_azimuthal_p_x, pin_azimuthal_p_y, 0.0));
714 : }
715 : }
716 : }
717 18783 : }
718 :
719 : void
720 12427 : 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 12427 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
733 30269 : for (unsigned int k = 0; k < (background_intervals); k++)
734 : {
735 : const Real background_corner_p_x =
736 17842 : background_corner_distance / (0.5 * corner_to_corner) * corner_p[0][0] *
737 17842 : (background_in +
738 17842 : biased_terms[k] * background_intervals * background_corner_radial_interval_length) /
739 17842 : background_corner_distance;
740 : const Real background_corner_p_y =
741 17842 : 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 17842 : 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 17842 : mesh.add_point(Point(background_corner_p_x, background_corner_p_y, 0.0));
749 :
750 80775 : for (unsigned int j = 1; j <= angle_number; j++)
751 : {
752 : const Real cell_boundary_p_x =
753 62933 : background_corner_distance / (0.5 * corner_to_corner) *
754 62933 : (corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
755 62933 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
756 3887 : : (azimuthal_tangent[j] / 2.0)));
757 : const Real cell_boundary_p_y =
758 62933 : background_corner_distance / (0.5 * corner_to_corner) *
759 62933 : (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 62933 : cell_boundary_p_x * background_in /
766 62933 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
767 : const Real pin_boundary_p_y =
768 62933 : cell_boundary_p_y * background_in /
769 62933 : 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 62933 : std::sqrt(Utility::pow<2>(cell_boundary_p_x - pin_boundary_p_x) +
774 62933 : Utility::pow<2>(cell_boundary_p_y - pin_boundary_p_y)) /
775 62933 : background_intervals;
776 : const Real background_azimuthal_p_x =
777 62933 : cell_boundary_p_x *
778 62933 : (background_in + biased_terms[k] * background_intervals * background_radial_interval) /
779 62933 : std::sqrt(Utility::pow<2>(cell_boundary_p_x) + Utility::pow<2>(cell_boundary_p_y));
780 : const Real background_azimuthal_p_y =
781 62933 : cell_boundary_p_y *
782 : (background_in + biased_terms[k] * background_intervals * background_radial_interval) /
783 62933 : 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 62933 : mesh.add_point(Point(background_azimuthal_p_x, background_azimuthal_p_y, 0.0));
787 : }
788 : }
789 12427 : }
790 :
791 : void
792 1812 : 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 1812 : azimuthal_tangent.size() == 0 ? num_sectors_per_side : (azimuthal_tangent.size() - 1);
803 : // Add nodes in ducts regions
804 : (*ducts_center_dist)
805 1812 : .push_back(0.5 * corner_to_corner); // add hex boundary as the last element in this vector
806 1812 : std::vector<Real> duct_radius_interval_length(ducts_layers.size());
807 :
808 : Real bin_radial_distance;
809 4173 : for (unsigned int l = 0; l < ducts_layers.size(); l++)
810 : {
811 2361 : duct_radius_interval_length[l] =
812 2361 : ((*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 7746 : for (unsigned int k = 0; k < ducts_layers[l]; k++)
817 : {
818 5385 : bin_radial_distance = ((*ducts_center_dist)[l] +
819 5385 : biased_terms[l][k] * ducts_layers[l] * duct_radius_interval_length[l]);
820 5385 : const Real pin_corner_p_x = corner_p[0][0] * bin_radial_distance / (0.5 * corner_to_corner);
821 5385 : 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 5385 : mesh.add_point(Point(pin_corner_p_x, pin_corner_p_y, 0.0));
826 :
827 31149 : for (unsigned int j = 1; j <= angle_number; j++)
828 : {
829 : const Real cell_boundary_p_x =
830 25764 : corner_p[0][0] + (corner_p[1][0] - corner_p[0][0]) *
831 25764 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
832 0 : : (azimuthal_tangent[j] / 2.0));
833 : const Real cell_boundary_p_y =
834 25764 : corner_p[0][1] + (corner_p[1][1] - corner_p[0][1]) *
835 : (azimuthal_tangent.size() == 0 ? ((Real)j / (Real)angle_number)
836 25764 : : (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 25764 : const Real pin_azimuthal_p_x =
840 25764 : cell_boundary_p_x * bin_radial_distance / (0.5 * corner_to_corner);
841 25764 : const Real pin_azimuthal_p_y =
842 25764 : 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 25764 : mesh.add_point(Point(pin_azimuthal_p_x, pin_azimuthal_p_y, 0.0));
847 : }
848 : }
849 : }
850 1812 : }
851 :
852 : void
853 2863 : 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 6161 : for (unsigned int i = 0; i < div_num - 1; i++)
869 : {
870 : unsigned int id_x = 0;
871 : unsigned int id_y = i;
872 7466 : for (unsigned int j = 0; j < 2 * i + 1; j++)
873 : {
874 4168 : std::unique_ptr<Elem> new_elem;
875 4168 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
876 : {
877 3808 : new_elem = std::make_unique<Quad4>();
878 3808 : new_elem->set_node(0, nodes[id_x][id_y]);
879 3808 : new_elem->set_node(3, nodes[id_x][id_y + 1]);
880 3808 : new_elem->set_node(2, nodes[id_x + 1][id_y + 1]);
881 3808 : new_elem->set_node(1, nodes[id_x + 1][id_y]);
882 3808 : new_elem->subdomain_id() = 1 + block_id_shift;
883 : }
884 : else // QUAD8/QUAD9
885 : {
886 360 : new_elem = std::make_unique<Quad8>();
887 360 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
888 : {
889 360 : new_elem = std::make_unique<Quad9>();
890 360 : new_elem->set_node(8, nodes[id_x * 2 + 1][id_y * 2 + 1]);
891 : }
892 360 : new_elem->set_node(0, nodes[id_x * 2][id_y * 2]);
893 360 : new_elem->set_node(3, nodes[id_x * 2][id_y * 2 + 2]);
894 360 : new_elem->set_node(2, nodes[id_x * 2 + 2][id_y * 2 + 2]);
895 360 : new_elem->set_node(1, nodes[id_x * 2 + 2][id_y * 2]);
896 360 : new_elem->set_node(4, nodes[id_x * 2 + 1][id_y * 2]);
897 360 : new_elem->set_node(5, nodes[id_x * 2 + 2][id_y * 2 + 1]);
898 360 : new_elem->set_node(6, nodes[id_x * 2 + 1][id_y * 2 + 2]);
899 360 : new_elem->set_node(7, nodes[id_x * 2][id_y * 2 + 1]);
900 360 : new_elem->subdomain_id() = 1 + block_id_shift;
901 : }
902 4168 : Elem * elem_Quad = mesh.add_elem(std::move(new_elem));
903 :
904 4168 : if (id_x == 0)
905 3298 : boundary_info.add_side(elem_Quad, 3, SLICE_BEGIN);
906 4168 : if (id_y == 0)
907 3298 : boundary_info.add_side(elem_Quad, 0, SLICE_END);
908 4168 : if (j < i)
909 435 : id_x++;
910 4168 : if (j >= i)
911 3733 : id_y--;
912 4168 : }
913 : }
914 : // This loop defines the outermost layer quad elements of the central region
915 9459 : for (unsigned int i = (div_num - 1) * (div_num - 1); i < div_num * div_num - 1; i++)
916 : {
917 6596 : std::unique_ptr<Elem> new_elem;
918 6596 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
919 : {
920 6236 : new_elem = std::make_unique<Quad4>();
921 6236 : new_elem->set_node(0, mesh.node_ptr(i));
922 6236 : new_elem->set_node(3, mesh.node_ptr(i + 2 * div_num - 1));
923 6236 : new_elem->set_node(2, mesh.node_ptr(i + 2 * div_num));
924 6236 : new_elem->set_node(1, mesh.node_ptr(i + 1));
925 : }
926 : else // QUAD8/QUAD9
927 : {
928 360 : new_elem = std::make_unique<Quad8>();
929 360 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
930 : {
931 360 : new_elem = std::make_unique<Quad9>();
932 360 : new_elem->set_node(8,
933 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
934 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1 +
935 : ((div_num - 1) * 4 + 1)));
936 : }
937 360 : new_elem->set_node(0,
938 360 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
939 : (i - (div_num - 1) * (div_num - 1)) * 2));
940 360 : new_elem->set_node(3,
941 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
942 360 : (i - (div_num - 1) * (div_num - 1)) * 2 +
943 : ((div_num - 1) * 4 + 1) * 2));
944 360 : new_elem->set_node(2,
945 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
946 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2 +
947 : ((div_num - 1) * 4 + 1) * 2));
948 360 : new_elem->set_node(1,
949 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
950 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2));
951 360 : new_elem->set_node(4,
952 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
953 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1));
954 360 : new_elem->set_node(5,
955 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
956 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 2 +
957 : ((div_num - 1) * 4 + 1)));
958 360 : new_elem->set_node(6,
959 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
960 360 : (i - (div_num - 1) * (div_num - 1)) * 2 + 1 +
961 : ((div_num - 1) * 4 + 1) * 2));
962 360 : new_elem->set_node(7,
963 : mesh.node_ptr((div_num - 1) * (div_num - 1) * 4 +
964 360 : (i - (div_num - 1) * (div_num - 1)) * 2 +
965 : ((div_num - 1) * 4 + 1)));
966 : }
967 :
968 6596 : Elem * elem_Quad = mesh.add_elem(std::move(new_elem));
969 6596 : elem_Quad->subdomain_id() = 1 + block_id_shift;
970 6596 : if (create_outward_interface_boundaries)
971 1224 : boundary_info.add_side(elem_Quad, 2, 1 + boundary_id_shift);
972 6596 : if (i == (div_num - 1) * (div_num - 1))
973 2863 : boundary_info.add_side(elem_Quad, 3, SLICE_BEGIN);
974 6596 : if (i == div_num * div_num - 2)
975 2863 : boundary_info.add_side(elem_Quad, 1, SLICE_END);
976 6596 : if (assign_external_boundary)
977 : {
978 360 : boundary_info.add_side(elem_Quad, 2, OUTER_SIDESET_ID);
979 360 : if (generate_side_specific_boundaries)
980 360 : boundary_info.add_side(
981 : elem_Quad,
982 : 2,
983 540 : (i < div_num * (div_num - 1) ? OUTER_SIDESET_ID : OUTER_SIDESET_ID_ALT) + side_index);
984 : }
985 6596 : }
986 2863 : }
987 :
988 : void
989 18998 : 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 18998 : 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 18998 : : ((azimuthal_tangent.size() - 1) / order);
1004 :
1005 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1006 55617 : for (unsigned int i = 1; i <= angle_number; i++)
1007 : {
1008 36619 : std::unique_ptr<Elem> new_elem;
1009 36619 : if (tri_elem_type == TRI_ELEM_TYPE::TRI3)
1010 : {
1011 33654 : new_elem = std::make_unique<Tri3>();
1012 33654 : new_elem->set_node(0, mesh.node_ptr(0));
1013 33654 : new_elem->set_node(2, mesh.node_ptr(i));
1014 33654 : new_elem->set_node(1, mesh.node_ptr(i + 1));
1015 : }
1016 : else // TRI6/TRI7
1017 : {
1018 2965 : new_elem = std::make_unique<Tri6>();
1019 2965 : if (tri_elem_type == TRI_ELEM_TYPE::TRI7)
1020 : {
1021 1257 : new_elem = std::make_unique<Tri7>();
1022 1257 : new_elem->set_node(6, mesh.node_ptr(i * 2));
1023 : }
1024 2965 : new_elem->set_node(0, mesh.node_ptr(0));
1025 2965 : new_elem->set_node(2, mesh.node_ptr(i * 2 + angle_number * order));
1026 2965 : new_elem->set_node(1, mesh.node_ptr((i + 1) * 2 + angle_number * order));
1027 2965 : new_elem->set_node(3, mesh.node_ptr(i * 2 + 1));
1028 2965 : new_elem->set_node(5, mesh.node_ptr(i * 2 - 1));
1029 2965 : new_elem->set_node(4, mesh.node_ptr(i * 2 + 1 + angle_number * order));
1030 : }
1031 :
1032 36619 : Elem * elem = mesh.add_elem(std::move(new_elem));
1033 36619 : if (create_outward_interface_boundaries)
1034 10277 : boundary_info.add_side(elem, 1, 1 + boundary_id_shift);
1035 36619 : elem->subdomain_id() = 1 + block_id_shift;
1036 36619 : if (i == 1)
1037 18998 : boundary_info.add_side(elem, 2, SLICE_BEGIN);
1038 36619 : if (i == angle_number)
1039 18998 : boundary_info.add_side(elem, 0, SLICE_END);
1040 36619 : if (assign_external_boundary)
1041 : {
1042 4158 : boundary_info.add_side(elem, 1, OUTER_SIDESET_ID);
1043 4158 : if (generate_side_specific_boundaries)
1044 1305 : boundary_info.add_side(elem,
1045 : 1,
1046 2430 : (i <= angle_number / 2 ? OUTER_SIDESET_ID : OUTER_SIDESET_ID_ALT) +
1047 : side_index);
1048 : }
1049 36619 : }
1050 18998 : }
1051 :
1052 : void
1053 21861 : 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 21861 : 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 21861 : : ((azimuthal_tangent.size() - 1) / order);
1070 :
1071 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1072 : unsigned int j = 0;
1073 66711 : for (unsigned int k = 0; k < (subdomain_rings.size()); k++)
1074 : {
1075 146603 : for (unsigned int m = 0; m < subdomain_rings[k]; m++)
1076 : {
1077 260593 : for (unsigned int i = 1; i <= angle_number; i++)
1078 : {
1079 158840 : std::unique_ptr<Elem> new_elem;
1080 158840 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD4)
1081 : {
1082 147901 : new_elem = std::make_unique<Quad4>();
1083 147901 : new_elem->set_node(0, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * j));
1084 147901 : new_elem->set_node(1, mesh.node_ptr(nodeid_shift + i + 1 + (angle_number + 1) * j));
1085 147901 : new_elem->set_node(2, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * (j + 1) + 1));
1086 147901 : new_elem->set_node(3, mesh.node_ptr(nodeid_shift + i + (angle_number + 1) * (j + 1)));
1087 : }
1088 : else // QUAD8/QUAD9
1089 : {
1090 10939 : new_elem = std::make_unique<Quad8>();
1091 10939 : if (quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
1092 : {
1093 9411 : new_elem = std::make_unique<Quad9>();
1094 9411 : new_elem->set_node(
1095 9411 : 8, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 2)));
1096 : }
1097 10939 : new_elem->set_node(
1098 : 0,
1099 10939 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 1)));
1100 10939 : new_elem->set_node(
1101 10939 : 1, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 1)));
1102 10939 : new_elem->set_node(
1103 10939 : 2, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 3)));
1104 10939 : new_elem->set_node(
1105 : 3,
1106 10939 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 3)));
1107 10939 : new_elem->set_node(
1108 : 4, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 1)));
1109 10939 : new_elem->set_node(
1110 10939 : 5, mesh.node_ptr(nodeid_shift + i * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 2)));
1111 10939 : new_elem->set_node(
1112 : 6, mesh.node_ptr(nodeid_shift + i * 2 + (angle_number * 2 + 1) * (j * 2 + 3)));
1113 10939 : new_elem->set_node(
1114 : 7,
1115 10939 : mesh.node_ptr(nodeid_shift + (i - 1) * 2 + 1 + (angle_number * 2 + 1) * (j * 2 + 2)));
1116 : }
1117 158840 : Elem * elem = mesh.add_elem(std::move(new_elem));
1118 158840 : if (i == 1)
1119 101753 : boundary_info.add_side(elem, 3, SLICE_BEGIN);
1120 158840 : if (i == angle_number)
1121 101753 : boundary_info.add_side(elem, 1, SLICE_END);
1122 :
1123 158840 : if (subdomain_rings[0] == 0)
1124 23024 : elem->subdomain_id() = k + 1 + block_id_shift;
1125 : else
1126 135816 : elem->subdomain_id() = k + 2 + block_id_shift;
1127 :
1128 158840 : if (m == 0 && create_inward_interface_boundaries && k > 0)
1129 720 : boundary_info.add_side(elem, 0, k * 2 + boundary_id_shift);
1130 158840 : if (m == (subdomain_rings[k] - 1))
1131 : {
1132 74106 : if (k == (subdomain_rings.size() - 1))
1133 : {
1134 38697 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID);
1135 38697 : if (generate_side_specific_boundaries)
1136 : {
1137 15545 : if (i <= angle_number / 2)
1138 3528 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID + side_index);
1139 : else
1140 12017 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID_ALT + side_index);
1141 : }
1142 : }
1143 35409 : else if (create_outward_interface_boundaries)
1144 22522 : boundary_info.add_side(elem, 2, k * 2 + 1 + boundary_id_shift);
1145 : }
1146 158840 : }
1147 101753 : j++;
1148 : }
1149 : }
1150 21861 : }
1151 :
1152 : std::unique_ptr<ReplicatedMesh>
1153 27762 : 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 27762 : auto mesh = buildReplicatedMesh(2);
1164 : std::pair<Real, Real> positions_p;
1165 :
1166 : // generate node positions
1167 99042 : for (unsigned int i = 0; i <= peripheral_invervals; i++)
1168 : {
1169 221984 : for (unsigned int j = 0; j <= num_sectors_per_side / 2; j++)
1170 : {
1171 150704 : positions_p = pointInterpolate(positions_inner[0].first,
1172 150704 : positions_inner[0].second,
1173 150704 : d_positions_outer[0].first,
1174 150704 : d_positions_outer[0].second,
1175 150704 : positions_inner[1].first,
1176 150704 : positions_inner[1].second,
1177 150704 : d_positions_outer[1].first,
1178 150704 : d_positions_outer[1].second,
1179 : i,
1180 : j,
1181 : num_sectors_per_side,
1182 : peripheral_invervals);
1183 150704 : mesh->add_point(Point(positions_p.first, positions_p.second, 0.0));
1184 : }
1185 150704 : for (unsigned int j = 1; j <= num_sectors_per_side / 2; j++)
1186 : {
1187 79424 : positions_p = pointInterpolate(positions_inner[1].first,
1188 79424 : positions_inner[1].second,
1189 79424 : d_positions_outer[1].first,
1190 79424 : d_positions_outer[1].second,
1191 79424 : positions_inner[2].first,
1192 79424 : positions_inner[2].second,
1193 79424 : d_positions_outer[2].first,
1194 79424 : d_positions_outer[2].second,
1195 : i,
1196 : j,
1197 : num_sectors_per_side,
1198 : peripheral_invervals);
1199 79424 : 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 71280 : for (unsigned int i = 0; i < peripheral_invervals; i++)
1207 : {
1208 141398 : for (unsigned int j = 0; j < num_sectors_per_side; j++)
1209 : {
1210 97880 : std::unique_ptr<Elem> new_elem;
1211 :
1212 97880 : new_elem = std::make_unique<Quad4>();
1213 97880 : new_elem->set_node(0, mesh->node_ptr(j + (num_sectors_per_side + 1) * (i)));
1214 97880 : new_elem->set_node(1, mesh->node_ptr(j + 1 + (num_sectors_per_side + 1) * (i)));
1215 97880 : new_elem->set_node(2, mesh->node_ptr(j + 1 + (num_sectors_per_side + 1) * (i + 1)));
1216 97880 : new_elem->set_node(3, mesh->node_ptr(j + (num_sectors_per_side + 1) * (i + 1)));
1217 :
1218 97880 : Elem * elem = mesh->add_elem(std::move(new_elem));
1219 :
1220 : // add subdoamin and boundary IDs
1221 97880 : elem->subdomain_id() = PERIPHERAL_ID_SHIFT + id_shift;
1222 97880 : if (i == 0)
1223 : {
1224 60968 : boundary_info.add_side(elem, 0, OUTER_SIDESET_ID);
1225 60968 : if (create_inward_interface_boundaries)
1226 0 : boundary_info.add_side(elem, 0, SLICE_ALT + id_shift * 2);
1227 : }
1228 97880 : if (i == peripheral_invervals - 1)
1229 : {
1230 60968 : boundary_info.add_side(elem, 2, OUTER_SIDESET_ID);
1231 60968 : if (create_outward_interface_boundaries)
1232 14976 : boundary_info.add_side(elem, 2, SLICE_ALT + id_shift * 2 + 1);
1233 : }
1234 97880 : if (j == 0)
1235 43518 : boundary_info.add_side(elem, 3, OUTER_SIDESET_ID);
1236 97880 : if (j == num_sectors_per_side - 1)
1237 43518 : boundary_info.add_side(elem, 1, OUTER_SIDESET_ID);
1238 97880 : }
1239 : }
1240 :
1241 : // convert element to second order if needed
1242 27762 : if (quad_elem_type != QUAD_ELEM_TYPE::QUAD4)
1243 : {
1244 : // full_ordered 2nd order element --> QUAD9, otherwise QUAD8
1245 1494 : const bool full_ordered = (quad_elem_type == QUAD_ELEM_TYPE::QUAD9);
1246 1494 : mesh->all_second_order(full_ordered);
1247 : }
1248 :
1249 27762 : return mesh;
1250 0 : }
1251 :
1252 : void
1253 36 : PolygonMeshGeneratorBase::adjustPeripheralQuadraticElements(
1254 : MeshBase & out_mesh, const QUAD_ELEM_TYPE boundary_quad_elem_type) const
1255 : {
1256 36 : 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 12132 : for (auto side_item : side_list)
1262 : {
1263 : boundary_id_type boundary_id = std::get<2>(side_item);
1264 12096 : dof_id_type elem_id = std::get<0>(side_item);
1265 :
1266 12096 : if (boundary_id == OUTER_SIDESET_ID)
1267 1872 : elem_set.insert(elem_id);
1268 : }
1269 :
1270 : // adjust nodes for outer boundary elements
1271 1908 : for (const auto elem_id : elem_set)
1272 : {
1273 1872 : 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 1872 : 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 1872 : out_mesh.add_point(pt_7, elem->node_ptr(7)->id());
1282 :
1283 : // adjust central node when using QUAD9
1284 1872 : if (boundary_quad_elem_type == QUAD_ELEM_TYPE::QUAD9)
1285 : {
1286 936 : Point pt_8 = elem->true_centroid();
1287 936 : out_mesh.add_point(pt_8, elem->node_ptr(8)->id());
1288 : }
1289 : }
1290 36 : }
1291 :
1292 : std::pair<Real, Real>
1293 230128 : 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 230128 : auto position_px_inner =
1307 230128 : (pi_1_x * (num_sectors_per_side / 2.0 - j) + pi_2_x * j) / (num_sectors_per_side / 2.0);
1308 230128 : auto position_py_inner =
1309 230128 : (pi_1_y * (num_sectors_per_side / 2.0 - j) + pi_2_y * j) / (num_sectors_per_side / 2.0);
1310 230128 : auto position_px_outer =
1311 230128 : (d_po_1_x * (num_sectors_per_side / 2.0 - j) + d_po_2_x * j) / (num_sectors_per_side / 2.0);
1312 230128 : auto position_py_outer =
1313 230128 : (d_po_1_y * (num_sectors_per_side / 2.0 - j) + d_po_2_y * j) / (num_sectors_per_side / 2.0);
1314 230128 : auto position_px = position_px_inner + position_px_outer * i / peripheral_intervals;
1315 230128 : auto position_py = position_py_inner + position_py_outer * i / peripheral_intervals;
1316 230128 : return std::make_pair(position_px, position_py);
1317 : }
1318 :
1319 : void
1320 506432 : PolygonMeshGeneratorBase::nodeCoordRotate(Real & x, Real & y, const Real theta) const
1321 : {
1322 506432 : const Real x_tmp = x;
1323 506432 : const Real y_tmp = y;
1324 506432 : x = x_tmp * std::cos(theta * M_PI / 180.0) - y_tmp * std::sin(theta * M_PI / 180.0);
1325 506432 : y = x_tmp * std::sin(theta * M_PI / 180.0) + y_tmp * std::cos(theta * M_PI / 180.0);
1326 506432 : }
1327 :
1328 : void
1329 2364 : 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 309144 : 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 151026 : Real & x = (*node_ptr)(0);
1342 : Real & y = (*node_ptr)(1);
1343 151026 : if (mesh_type == CORNER_MESH)
1344 : {
1345 97614 : nodeCoordRotate(x, y, orientation);
1346 97614 : if (x >= 0.0 && y > y_max_0)
1347 5598 : y = y - y_max_0 + y_max_n;
1348 92016 : else if (x >= 0.0 && y >= y_min)
1349 14507 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1350 77509 : else if (y > -x / std::tan(unit_angle / 360.0 * M_PI) + tols && y > y_max_0)
1351 : {
1352 342 : x /= y;
1353 342 : y = y - y_max_0 + y_max_n;
1354 342 : x *= y;
1355 : }
1356 77167 : else if (y > -x / std::tan(unit_angle / 360.0 * M_PI) + tols && y >= y_min)
1357 : {
1358 2679 : x /= y;
1359 2679 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1360 2679 : x *= y;
1361 : }
1362 97614 : nodeCoordRotate(x, y, -orientation);
1363 :
1364 97614 : nodeCoordRotate(x, y, orientation - unit_angle);
1365 97614 : if (x <= 0 && y > y_max_0)
1366 5682 : y = y - y_max_0 + y_max_n;
1367 91932 : else if (x <= 0 && y >= y_min)
1368 13426 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1369 78506 : else if (y >= x / std::tan(unit_angle / 360.0 * M_PI) - tols && y > y_max_0)
1370 : {
1371 2722 : x /= y;
1372 2722 : y = y - y_max_0 + y_max_n;
1373 2722 : x *= y;
1374 : }
1375 75784 : else if (y >= x / std::tan(unit_angle / 360.0 * M_PI) - tols && y >= y_min)
1376 : {
1377 8418 : x /= y;
1378 8418 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1379 8418 : x *= y;
1380 : }
1381 97614 : nodeCoordRotate(x, y, unit_angle - orientation);
1382 : }
1383 : else
1384 : {
1385 53412 : nodeCoordRotate(x, y, orientation);
1386 53412 : if (y > y_max_0)
1387 12420 : y = y - y_max_0 + y_max_n;
1388 40992 : else if (y >= y_min)
1389 18963 : y = (y - y_min) / (y_max_0 - y_min) * (y_max_n - y_min) + y_min;
1390 53412 : nodeCoordRotate(x, y, -orientation);
1391 : }
1392 2364 : }
1393 2364 : }
1394 :
1395 : std::pair<Real, Real>
1396 15103 : 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 15103 : const Real x1 = p1.first;
1402 15103 : const Real y1 = p1.second;
1403 15103 : const Real x2 = p2.first;
1404 15103 : const Real y2 = p2.second;
1405 15103 : const Real x3 = p3.first;
1406 15103 : const Real y3 = p3.second;
1407 15103 : const Real x4 = p4.first;
1408 15103 : const Real y4 = p4.second;
1409 :
1410 15103 : Real x = -((x1 - x2) * (y3 * x4 - x3 * y4) - (x3 - x4) * (y1 * x2 - x1 * y2)) /
1411 15103 : ((y1 - y2) * (x3 - x4) - (y3 - y4) * (x1 - x2));
1412 15103 : 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 15103 : return std::make_pair(x, y);
1416 : }
1417 :
1418 : std::vector<Real>
1419 1636 : 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 1636 : mesh.get_boundary_info().build_side_list();
1433 1636 : mesh.get_boundary_info().build_node_list_from_side_list();
1434 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
1435 1636 : 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 1636 : const Real mid_azi = lower_azi <= upper_azi ? (lower_azi + upper_azi) / 2.0
1444 42 : : (lower_azi + upper_azi + 360.0) / 2.0;
1445 306886 : for (unsigned int i = 0; i < node_list.size(); ++i)
1446 305250 : if (std::get<1>(node_list[i]) == bid)
1447 : {
1448 61419 : bd_x_list.push_back((mesh.node_ref(std::get<0>(node_list[i])))(0));
1449 61419 : bd_y_list.push_back((mesh.node_ref(std::get<0>(node_list[i])))(1));
1450 61419 : bd_p_list.push_back((mesh.node_ref(std::get<0>(node_list[i]))));
1451 : }
1452 :
1453 1636 : if (calculate_origin)
1454 : {
1455 1636 : const Point origin_pt = MooseMeshUtils::meshCentroidCalculator(mesh);
1456 1636 : origin_x = origin_pt(0);
1457 1636 : 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 63055 : for (unsigned int i = 0; i < bd_x_list.size(); ++i)
1468 : {
1469 61419 : tmp_azi = atan2(bd_y_list[i] - origin_y, bd_x_list[i] - origin_x) * 180.0 / M_PI;
1470 61419 : if ((lower_azi <= upper_azi && (tmp_azi >= lower_azi - tol && tmp_azi <= upper_azi + tol)) ||
1471 2232 : (lower_azi > upper_azi && (tmp_azi >= lower_azi - tol || tmp_azi <= upper_azi + tol)))
1472 : {
1473 : azi_point_pairs.push_back(
1474 18968 : std::make_pair(return_type == ANGLE_DEGREE
1475 18968 : ? (tmp_azi - mid_azi)
1476 3027 : : (1.0 + std::cos(M_PI / num_sides) / std::sin(M_PI / num_sides) *
1477 3027 : std::tan((tmp_azi - mid_azi) / 180.0 * M_PI)),
1478 : bd_p_list[i]));
1479 : }
1480 : }
1481 1636 : 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 20604 : it != end;
1487 : it++)
1488 : {
1489 18968 : azimuthal_output.push_back(std::move(it->first));
1490 18968 : boundary_points.push_back(std::move(it->second));
1491 : }
1492 :
1493 1636 : return azimuthal_output;
1494 1636 : }
1495 :
1496 : std::vector<Real>
1497 800 : 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 1600 : tol);
1520 800 : }
1521 :
1522 : std::vector<std::vector<Real>>
1523 87444 : 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 152290 : for (unsigned int i = 0; i < radial_biases.size(); i++)
1531 129692 : 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 87444 : return bias_terms_vec;
1542 0 : }
1543 :
1544 : std::vector<Real>
1545 196339 : 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 198364 : 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 4050 : ? ((Real)(i + 1) * inner_boundary_layer_params.fraction /
1568 0 : (Real)inner_boundary_layer_params.intervals)
1569 2025 : : ((1.0 - std::pow(inner_boundary_layer_params.bias, (Real)(i + 1))) /
1570 2025 : (1.0 - std::pow(inner_boundary_layer_params.bias,
1571 : (Real)(inner_boundary_layer_params.intervals))) *
1572 : inner_boundary_layer_params.fraction));
1573 443581 : for (unsigned int i = 0; i < intervals; i++)
1574 247242 : biased_terms.push_back(inner_boundary_layer_params.fraction +
1575 : (MooseUtils::absoluteFuzzyEqual(radial_bias, 1.0)
1576 247242 : ? ((Real)(i + 1) *
1577 244477 : (1.0 - inner_boundary_layer_params.fraction -
1578 : outer_boundary_layer_params.fraction) /
1579 244477 : (Real)intervals)
1580 2765 : : ((1.0 - std::pow(radial_bias, (Real)(i + 1))) /
1581 2765 : (1.0 - std::pow(radial_bias, (Real)(intervals))) *
1582 2765 : (1.0 - inner_boundary_layer_params.fraction -
1583 : outer_boundary_layer_params.fraction))));
1584 199039 : for (unsigned int i = 0; i < outer_boundary_layer_params.intervals; i++)
1585 : biased_terms.push_back(
1586 5400 : 1.0 - outer_boundary_layer_params.fraction +
1587 : (MooseUtils::absoluteFuzzyEqual(outer_boundary_layer_params.bias, 1.0)
1588 2700 : ? ((Real)(i + 1) * outer_boundary_layer_params.fraction /
1589 0 : (Real)outer_boundary_layer_params.intervals)
1590 2700 : : ((1.0 - std::pow(outer_boundary_layer_params.bias, (Real)(i + 1))) /
1591 2700 : (1.0 - std::pow(outer_boundary_layer_params.bias,
1592 : (Real)(outer_boundary_layer_params.intervals))) *
1593 : outer_boundary_layer_params.fraction)));
1594 196339 : return biased_terms;
1595 0 : }
1596 :
1597 : void
1598 5298 : PolygonMeshGeneratorBase::addRingAndSectorIDParams(InputParameters & params)
1599 : {
1600 10596 : 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 10596 : 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 10596 : MooseEnum ring_id_option("block_wise ring_wise", "block_wise");
1607 10596 : params.addParam<MooseEnum>(
1608 : "ring_id_assign_type", ring_id_option, "Type of ring ID assignment: block_wise or ring_wise");
1609 10596 : params.addParamNamesToGroup("sector_id_name ring_id_name ring_id_assign_type", "Ring/Sector IDs");
1610 5298 : }
1611 :
1612 : void
1613 72 : 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 72 : const auto extra_id_index = mesh.add_elem_integer(id_name);
1619 : // vector to store sector ids for each element
1620 72 : auto elem_it = mesh.elements_begin();
1621 : unsigned int id = 1;
1622 : // starting element id of the current sector
1623 387 : for (unsigned int is = 0; is < num_sides; ++is)
1624 : {
1625 : // number of elements in the current sector
1626 : unsigned int nelem_sector =
1627 315 : mesh.n_elem() * num_sectors_per_side[is] /
1628 315 : (accumulate(num_sectors_per_side.begin(), num_sectors_per_side.end(), 0));
1629 : // assign sector ids to mesh
1630 15075 : for (unsigned i = 0; i < nelem_sector; ++i, ++elem_it)
1631 7380 : (*elem_it)->set_extra_integer(extra_id_index, id);
1632 : // update sector id
1633 315 : ++id;
1634 : }
1635 72 : }
1636 :
1637 : void
1638 72 : 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 72 : const auto extra_id_index = mesh.add_elem_integer(id_name);
1649 72 : auto elem_it = mesh.elements_begin();
1650 387 : for (unsigned int is = 0; is < num_sides; ++is)
1651 : {
1652 : // number of elements in the current sector
1653 315 : unsigned int nelem = mesh.n_elem() * num_sectors_per_side[is] /
1654 315 : (accumulate(num_sectors_per_side.begin(), num_sectors_per_side.end(), 0));
1655 315 : if (!ring_wise_id)
1656 : {
1657 999 : for (unsigned int ir : index_range(ring_intervals))
1658 : {
1659 : // number of elements in the current ring and sector
1660 720 : 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 720 : 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 10008 : for (unsigned i = 0; i < nelem_annular_ring; ++i, ++elem_it)
1668 4644 : (*elem_it)->set_extra_integer(extra_id_index, ir + 1);
1669 : // update number of elements in background region of current side.
1670 720 : nelem -= nelem_annular_ring;
1671 : }
1672 : }
1673 : else
1674 : {
1675 : unsigned int ir = 0;
1676 144 : for (unsigned int ir0 : index_range(ring_intervals))
1677 : {
1678 288 : for (unsigned int ir1 = 0; ir1 < ring_intervals[ir0]; ++ir1)
1679 : {
1680 : // number of elements in the current ring and sector
1681 180 : 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 180 : 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 1620 : for (unsigned i = 0; i < nelem_annular_ring; ++i, ++elem_it)
1688 720 : (*elem_it)->set_extra_integer(extra_id_index, ir + 1);
1689 : // update ring id
1690 180 : ++ir;
1691 : // update number of elements in background region of current side.
1692 180 : nelem -= nelem_annular_ring;
1693 : }
1694 : }
1695 : }
1696 : // assign ring id of 0 to the background region
1697 5499 : for (unsigned i = 0; i < nelem; ++i, ++elem_it)
1698 2592 : (*elem_it)->set_extra_integer(extra_id_index, 0);
1699 : }
1700 72 : }
1701 :
1702 : void
1703 207 : 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 531 : for (const auto id : boundary_ids)
1711 : {
1712 :
1713 324 : const boundary_id_type old_id = (!reverse) ? id : id + id_shift;
1714 324 : const boundary_id_type new_id = (!reverse) ? id + id_shift : id;
1715 : auto it = existing_boundary_ids.find(old_id);
1716 324 : if (it != existing_boundary_ids.end())
1717 324 : MooseMesh::changeBoundaryId(mesh, old_id, new_id, true);
1718 : }
1719 207 : }
1720 :
1721 : std::set<boundary_id_type>
1722 1592 : 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 1592 : if (use_interface_boundary_id_shift)
1734 : {
1735 72 : for (const auto i : make_range(pattern.size()))
1736 198 : for (const auto j : make_range(pattern[i].size()))
1737 : {
1738 144 : const auto & ids = input_interface_boundary_ids[pattern[i][j]];
1739 351 : for (const auto & id : ids)
1740 : {
1741 207 : const boundary_id_type new_id = id + interface_boundary_id_shift_pattern[i][j];
1742 : auto it = boundary_ids.find(new_id);
1743 207 : if (it != boundary_ids.end())
1744 207 : interface_boundary_ids.insert(new_id);
1745 : }
1746 : }
1747 : }
1748 : else
1749 : {
1750 4452 : for (const auto & ids : input_interface_boundary_ids)
1751 4720 : for (const auto & id : ids)
1752 : {
1753 : auto it = boundary_ids.find(id);
1754 1842 : if (it != boundary_ids.end())
1755 1779 : interface_boundary_ids.insert(id);
1756 : }
1757 : }
1758 : // add unshifted interface boundary ids for the duct & background regions
1759 1592 : if (create_interface_boundary_id)
1760 1902 : for (const auto i : make_range(num_extra_layers))
1761 : {
1762 938 : boundary_id_type id = SLICE_ALT + i * 2 + 1;
1763 : auto it = boundary_ids.find(id);
1764 938 : if (it != boundary_ids.end())
1765 315 : interface_boundary_ids.insert(id);
1766 938 : id = SLICE_ALT + i * 2;
1767 : it = boundary_ids.find(id);
1768 938 : if (it != boundary_ids.end())
1769 0 : interface_boundary_ids.insert(id);
1770 : }
1771 1592 : return interface_boundary_ids;
1772 : }
1773 :
1774 : PolygonMeshGeneratorBase::multiBdryLayerParams
1775 87444 : PolygonMeshGeneratorBase::modifiedMultiBdryLayerParamsCreator(
1776 : const multiBdryLayerParams & original_multi_bdry_layer_params, const unsigned int order) const
1777 : {
1778 87444 : 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 64846 : [&order](unsigned int & n) { n *= order; });
1782 87444 : std::for_each(mod_multi_bdry_layer_params.biases.begin(),
1783 : mod_multi_bdry_layer_params.biases.end(),
1784 64846 : [&order](Real & n) { n = std::pow(n, 1.0 / order); });
1785 87444 : return mod_multi_bdry_layer_params;
1786 : }
1787 :
1788 : PolygonMeshGeneratorBase::singleBdryLayerParams
1789 43722 : PolygonMeshGeneratorBase::modifiedSingleBdryLayerParamsCreator(
1790 : const singleBdryLayerParams & original_single_bdry_layer_params, const unsigned int order) const
1791 : {
1792 43722 : singleBdryLayerParams mod_single_bdry_layer_params(original_single_bdry_layer_params);
1793 43722 : mod_single_bdry_layer_params.intervals *= order;
1794 43722 : mod_single_bdry_layer_params.bias = std::pow(mod_single_bdry_layer_params.bias, 1.0 / order);
1795 43722 : 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 : }
|
