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 "PolygonConcentricCircleMeshGeneratorBase.h"
11 : #include "libmesh/mesh_smoother_laplace.h"
12 : #include "MooseUtils.h"
13 : #include "PolygonalMeshGenerationUtils.h"
14 : #include "MooseMeshUtils.h"
15 :
16 : #include "libmesh/parsed_function.h"
17 : #include "libmesh/poly2tri_triangulator.h"
18 :
19 : #include <cmath>
20 :
21 : InputParameters
22 5166 : PolygonConcentricCircleMeshGeneratorBase::validParams()
23 : {
24 5166 : InputParameters params = ConcentricCircleGeneratorBase::validParams();
25 10332 : params.addRequiredRangeCheckedParam<std::vector<unsigned int>>(
26 : "num_sectors_per_side",
27 : "num_sectors_per_side>0",
28 : "Number of azimuthal sectors per polygon side (rotating counterclockwise from top right "
29 : "face).");
30 15498 : params.addRangeCheckedParam<unsigned int>(
31 : "background_intervals",
32 10332 : 1,
33 : "background_intervals>0",
34 : "Number of radial meshing intervals in background region (area "
35 : "between rings and ducts) excluding the background's boundary layers.");
36 15498 : params.addRangeCheckedParam<Real>(
37 : "background_radial_bias",
38 10332 : 1.0,
39 : "background_radial_bias>0",
40 : "Value used to create biasing in radial meshing for background region.");
41 15498 : params.addRangeCheckedParam<Real>(
42 : "background_inner_boundary_layer_width",
43 10332 : 0.0,
44 : "background_inner_boundary_layer_width>=0",
45 : "Width of background region that is assigned to be the inner boundary layer.");
46 15498 : params.addRangeCheckedParam<unsigned int>(
47 : "background_inner_boundary_layer_intervals",
48 10332 : 1,
49 : "background_inner_boundary_layer_intervals>0",
50 : "Number of radial intervals of the background inner boundary layer");
51 15498 : params.addRangeCheckedParam<Real>(
52 : "background_inner_boundary_layer_bias",
53 10332 : 1.0,
54 : "background_inner_boundary_layer_bias>0",
55 : "Growth factor used for mesh biasing of the background inner boundary layer.");
56 15498 : params.addRangeCheckedParam<Real>(
57 : "background_outer_boundary_layer_width",
58 10332 : 0.0,
59 : "background_outer_boundary_layer_width>=0",
60 : "Width of background region that is assigned to be the outer boundary layer.");
61 15498 : params.addRangeCheckedParam<unsigned int>(
62 : "background_outer_boundary_layer_intervals",
63 10332 : 1,
64 : "background_outer_boundary_layer_intervals>0",
65 : "Number of radial intervals of the background outer boundary layer");
66 15498 : params.addRangeCheckedParam<Real>(
67 : "background_outer_boundary_layer_bias",
68 10332 : 1.0,
69 : "background_outer_boundary_layer_bias>0",
70 : "Growth factor used for mesh biasing of the background outer boundary layer.");
71 10332 : params.addParam<std::vector<subdomain_id_type>>(
72 : "background_block_ids", "Optional customized block id for the background block.");
73 10332 : params.addParam<std::vector<SubdomainName>>(
74 : "background_block_names", "Optional customized block names for the background block.");
75 10332 : params.addParam<std::vector<Real>>(
76 : "duct_sizes", "Distance(s) from polygon center to duct(s) inner boundaries.");
77 10332 : MooseEnum duct_sizes_style("apothem radius", "radius");
78 10332 : params.addParam<MooseEnum>(
79 : "duct_sizes_style",
80 : duct_sizes_style,
81 : "Style in which polygon center to duct inner boundary distance is "
82 5166 : "given (apothem = center to face, radius = center to vertex). Options: " +
83 5166 : duct_sizes_style.getRawNames());
84 10332 : params.addRangeCheckedParam<std::vector<unsigned int>>(
85 : "duct_intervals",
86 : "duct_intervals>0",
87 : "Number of meshing intervals in each enclosing duct excluding duct boundary layers.");
88 10332 : params.addRangeCheckedParam<std::vector<Real>>(
89 : "duct_radial_biases",
90 : "duct_radial_biases>0",
91 : "Values used to create biasing in radial meshing for duct regions.");
92 10332 : params.addRangeCheckedParam<std::vector<Real>>(
93 : "duct_inner_boundary_layer_widths",
94 : "duct_inner_boundary_layer_widths>=0",
95 : "Widths of duct regions that are assigned to be the inner boundary layers.");
96 10332 : params.addParam<std::vector<unsigned int>>(
97 : "duct_inner_boundary_layer_intervals",
98 : "Number of radial intervals of the duct inner boundary layers");
99 10332 : params.addRangeCheckedParam<std::vector<Real>>(
100 : "duct_inner_boundary_layer_biases",
101 : "duct_inner_boundary_layer_biases>0",
102 : "Growth factors used for mesh biasing of the duct inner boundary layers.");
103 10332 : params.addRangeCheckedParam<std::vector<Real>>(
104 : "duct_outer_boundary_layer_widths",
105 : "duct_outer_boundary_layer_widths>=0",
106 : "Widths of duct regions that are assigned to be the outer boundary layers.");
107 10332 : params.addParam<std::vector<unsigned int>>(
108 : "duct_outer_boundary_layer_intervals",
109 : "Number of radial intervals of the duct outer boundary layers");
110 10332 : params.addRangeCheckedParam<std::vector<Real>>(
111 : "duct_outer_boundary_layer_biases",
112 : "duct_outer_boundary_layer_biases>0",
113 : "Growth factors used for mesh biasing of the duct outer boundary layers.");
114 10332 : params.addParam<std::vector<subdomain_id_type>>(
115 : "duct_block_ids", "Optional customized block ids for each duct geometry block.");
116 10332 : params.addParam<std::vector<SubdomainName>>(
117 : "duct_block_names", "Optional customized block names for each duct geometry block.");
118 10332 : params.addParam<bool>("uniform_mesh_on_sides",
119 10332 : false,
120 : "Whether the side elements are reorganized to have a uniform size.");
121 10332 : params.addParam<unsigned int>("smoothing_max_it",
122 10332 : 0,
123 : "Number of Laplacian smoothing iterations. This number is "
124 : "disregarded when duct_sizes is present.");
125 10332 : params.addParam<bool>(
126 : "flat_side_up",
127 10332 : false,
128 : "Whether to rotate the generated polygon mesh to ensure that one flat side faces up.");
129 :
130 10332 : params.addParam<bool>(
131 10332 : "quad_center_elements", false, "Whether the center elements are quad or triangular.");
132 10332 : params.addRangeCheckedParam<Real>(
133 : "center_quad_factor",
134 : "center_quad_factor>0¢er_quad_factor<1",
135 : "A fractional radius factor used to determine the radial positions of transition nodes in "
136 : "the center region meshed by quad elements.");
137 10332 : params.addParam<bool>("replace_inner_ring_with_delaunay_mesh",
138 10332 : false,
139 : "True to replace the inner ring mesh with a Delaunay unstructured mesh");
140 10332 : params.addParam<std::string>(
141 : "inner_ring_desired_area",
142 5166 : std::string(),
143 : "Desired area as a function of x,y; omit to use the default constant area (square of the "
144 : "smallest side length on the ring circle)");
145 :
146 10332 : params.addParamNamesToGroup(
147 : "background_block_ids background_block_names duct_block_ids duct_block_names",
148 : "Customized Subdomain/Boundary");
149 10332 : params.addParamNamesToGroup("num_sectors_per_side background_intervals duct_intervals "
150 : "uniform_mesh_on_sides",
151 : "General Mesh Density");
152 10332 : params.addParamNamesToGroup(
153 : "ring_radial_biases ring_inner_boundary_layer_biases ring_inner_boundary_layer_widths "
154 : "ring_inner_boundary_layer_intervals ring_outer_boundary_layer_biases "
155 : "ring_outer_boundary_layer_widths ring_outer_boundary_layer_intervals",
156 : "Mesh Boundary Layers and Biasing Options");
157 10332 : params.addParamNamesToGroup("quad_center_elements center_quad_factor "
158 : "replace_inner_ring_with_delaunay_mesh inner_ring_desired_area ",
159 : "Inner ring");
160 :
161 5166 : addRingAndSectorIDParams(params);
162 5166 : params.addClassDescription("This PolygonConcentricCircleMeshGeneratorBase object is a base class "
163 : "to be inherited for polygon mesh generators.");
164 :
165 5166 : return params;
166 5166 : }
167 :
168 2626 : PolygonConcentricCircleMeshGeneratorBase::PolygonConcentricCircleMeshGeneratorBase(
169 2626 : const InputParameters & parameters)
170 : : ConcentricCircleGeneratorBase(parameters),
171 7836 : _num_sides(isParamValid("num_sides")
172 7425 : ? getParam<unsigned int>("num_sides")
173 3533 : : (isParamValid("hexagon_size") ? (unsigned int)HEXAGON_NUM_SIDES
174 : : (unsigned int)SQUARE_NUM_SIDES)),
175 5224 : _duct_sizes_style(getParam<MooseEnum>("duct_sizes_style").template getEnum<PolygonSizeStyle>()),
176 8210 : _duct_sizes(isParamValid("duct_sizes") ? getParam<std::vector<Real>>("duct_sizes")
177 : : std::vector<Real>()),
178 5972 : _duct_intervals(isParamValid("duct_intervals")
179 2612 : ? getParam<std::vector<unsigned int>>("duct_intervals")
180 : : std::vector<unsigned int>()),
181 5246 : _duct_radial_biases(isParamValid("duct_radial_biases")
182 2612 : ? getParam<std::vector<Real>>("duct_radial_biases")
183 2612 : : std::vector<Real>(_duct_intervals.size(), 1.0)),
184 7960 : _duct_inner_boundary_layer_params(
185 2612 : {isParamValid("duct_inner_boundary_layer_widths")
186 2612 : ? getParam<std::vector<Real>>("duct_inner_boundary_layer_widths")
187 2588 : : std::vector<Real>(_duct_intervals.size(), 0.0),
188 : std::vector<Real>(),
189 5224 : isParamValid("duct_inner_boundary_layer_intervals")
190 2612 : ? getParam<std::vector<unsigned int>>("duct_inner_boundary_layer_intervals")
191 2592 : : std::vector<unsigned int>(_duct_intervals.size(), 0),
192 5224 : isParamValid("duct_inner_boundary_layer_biases")
193 2612 : ? getParam<std::vector<Real>>("duct_inner_boundary_layer_biases")
194 2612 : : std::vector<Real>(_duct_intervals.size(), 0.0)}),
195 7964 : _duct_outer_boundary_layer_params(
196 2612 : {isParamValid("duct_outer_boundary_layer_widths")
197 2612 : ? getParam<std::vector<Real>>("duct_outer_boundary_layer_widths")
198 2588 : : std::vector<Real>(_duct_intervals.size(), 0.0),
199 : std::vector<Real>(),
200 5224 : isParamValid("duct_outer_boundary_layer_intervals")
201 2612 : ? getParam<std::vector<unsigned int>>("duct_outer_boundary_layer_intervals")
202 2590 : : std::vector<unsigned int>(_duct_intervals.size(), 0),
203 5224 : isParamValid("duct_outer_boundary_layer_biases")
204 2612 : ? getParam<std::vector<Real>>("duct_outer_boundary_layer_biases")
205 2612 : : std::vector<Real>(_duct_intervals.size(), 0.0)}),
206 5924 : _duct_block_ids(isParamValid("duct_block_ids")
207 2612 : ? getParam<std::vector<subdomain_id_type>>("duct_block_ids")
208 : : std::vector<subdomain_id_type>()),
209 5906 : _duct_block_names(isParamValid("duct_block_names")
210 2612 : ? getParam<std::vector<SubdomainName>>("duct_block_names")
211 : : std::vector<SubdomainName>()),
212 5224 : _has_rings(isParamValid("ring_radii")),
213 5224 : _has_ducts(isParamValid("duct_sizes")),
214 2612 : _polygon_size_style(
215 2612 : isParamValid("polygon_size_style")
216 5224 : ? getParam<MooseEnum>("polygon_size_style").template getEnum<PolygonSizeStyle>()
217 3023 : : (isParamValid("hexagon_size_style")
218 3023 : ? getParam<MooseEnum>("hexagon_size_style").template getEnum<PolygonSizeStyle>()
219 2810 : : getParam<MooseEnum>("square_size_style")
220 : .template getEnum<PolygonSizeStyle>())),
221 2612 : _polygon_size(isParamValid("polygon_size")
222 7425 : ? getParam<Real>("polygon_size")
223 4157 : : (isParamValid("hexagon_size") ? getParam<Real>("hexagon_size")
224 2810 : : (getParam<Real>("square_size") /
225 : 2.0))), // square size is twice as apothem
226 5224 : _num_sectors_per_side(getParam<std::vector<unsigned int>>("num_sectors_per_side")),
227 5224 : _background_intervals(getParam<unsigned int>("background_intervals")),
228 : // background is usually a single block; however, when there are no rings, background has two
229 : // blocks.
230 5224 : _background_radial_bias(getParam<Real>("background_radial_bias")),
231 7836 : _background_inner_boundary_layer_params(
232 2612 : {getParam<Real>("background_inner_boundary_layer_width"),
233 : 0.0,
234 5224 : getParam<Real>("background_inner_boundary_layer_width") > 0.0
235 5253 : ? getParam<unsigned int>("background_inner_boundary_layer_intervals")
236 : : 0,
237 5224 : getParam<Real>("background_inner_boundary_layer_bias")}),
238 7836 : _background_outer_boundary_layer_params(
239 2612 : {getParam<Real>("background_outer_boundary_layer_width"),
240 : 0.0,
241 5224 : getParam<Real>("background_outer_boundary_layer_width") > 0.0
242 5253 : ? getParam<unsigned int>("background_outer_boundary_layer_intervals")
243 : : 0,
244 5224 : getParam<Real>("background_outer_boundary_layer_bias")}),
245 9298 : _background_block_ids(isParamValid("background_block_ids")
246 2612 : ? getParam<std::vector<subdomain_id_type>>("background_block_ids")
247 : : std::vector<subdomain_id_type>()),
248 8844 : _background_block_names(isParamValid("background_block_names")
249 2612 : ? getParam<std::vector<SubdomainName>>("background_block_names")
250 : : std::vector<SubdomainName>()),
251 5224 : _uniform_mesh_on_sides(getParam<bool>("uniform_mesh_on_sides")),
252 5224 : _quad_center_elements(getParam<bool>("quad_center_elements")),
253 5246 : _center_quad_factor(isParamValid("center_quad_factor") ? getParam<Real>("center_quad_factor")
254 : : 0.0),
255 5224 : _smoothing_max_it(getParam<unsigned int>("smoothing_max_it")),
256 6026 : _sides_to_adapt(isParamValid("sides_to_adapt")
257 2612 : ? getParam<std::vector<unsigned int>>("sides_to_adapt")
258 : : std::vector<unsigned int>()),
259 2612 : _node_id_background_meta(declareMeshProperty<dof_id_type>("node_id_background_meta", 0)),
260 5238 : _is_control_drum_meta(declareMeshProperty<bool>("is_control_drum_meta", false))
261 : {
262 7836 : declareMeshProperty<bool>("flat_side_up", getParam<bool>("flat_side_up"));
263 2612 : declareMeshProperty<unsigned int>("background_intervals_meta", 0);
264 5224 : declareMeshProperty<Real>("pattern_pitch_meta", 0.0);
265 5224 : declareMeshProperty<std::vector<Real>>("azimuthal_angle_meta", std::vector<Real>());
266 2612 : declareMeshProperty<Real>("max_radius_meta", 0.0);
267 :
268 2612 : const unsigned short tri_order = _tri_elem_type == TRI_ELEM_TYPE::TRI3 ? 1 : 2;
269 2612 : const unsigned short quad_order = _quad_elem_type == QUAD_ELEM_TYPE::QUAD4 ? 1 : 2;
270 : // 1. If the central elements are quad, then no triangular elements are generated;
271 : // 2. If the generated mesh has only one radial layer of triangular elements, then no
272 : // quad elements are generated; (For PCCMG, that layer must be background)
273 : // 3. Otherwise, both types of elements are generated.
274 2612 : _order = quad_order;
275 2612 : if (_quad_center_elements)
276 : {
277 606 : if (tri_order != quad_order)
278 18 : _tri_elem_type = quad_order == 1 ? TRI_ELEM_TYPE::TRI3 : TRI_ELEM_TYPE::TRI6;
279 : }
280 440 : else if (_ring_radii.empty() && _background_intervals == 1 &&
281 247 : _background_inner_boundary_layer_params.intervals == 0 &&
282 2253 : _background_outer_boundary_layer_params.intervals == 0 && _duct_sizes.empty())
283 : {
284 238 : _order = tri_order;
285 238 : if (tri_order != quad_order)
286 18 : _quad_elem_type = tri_order == 1 ? QUAD_ELEM_TYPE::QUAD4 : QUAD_ELEM_TYPE::QUAD9;
287 : }
288 1768 : else if (tri_order != quad_order)
289 2 : paramError("tri_element_type",
290 : "the element types of triangular and quadrilateral elements must be compatible if "
291 : "both types of elements are generated.");
292 :
293 : // This error message is only reserved for future derived classes. Neither of the current derived
294 : // classes will trigger this error.
295 2610 : if (!_sides_to_adapt.empty() && _num_sides != HEXAGON_NUM_SIDES && _num_sides != SQUARE_NUM_SIDES)
296 0 : paramError("sides_to_adapt", "If provided, the generated mesh must be a hexagon or a square.");
297 5220 : _pitch = 2.0 * (_polygon_size_style == PolygonSizeStyle::apothem
298 2610 : ? _polygon_size
299 0 : : _polygon_size * std::cos(M_PI / Real(_num_sides)));
300 5220 : declareMeshProperty<Real>("pitch_meta", _pitch);
301 2610 : if (_inward_interface_boundary_names.size() > 0 &&
302 0 : _inward_interface_boundary_names.size() != _duct_sizes.size() + _ring_radii.size())
303 0 : paramError("inward_interface_boundary_names",
304 : "If provided, the length of this parameter must be identical to the total number of "
305 : "interfaces.");
306 2610 : if (_outward_interface_boundary_names.size() > 0 &&
307 38 : _outward_interface_boundary_names.size() != _duct_sizes.size() + _ring_radii.size())
308 2 : paramError("outward_interface_boundary_names",
309 : "If provided, the length of this parameter must be identical to the total number of "
310 : "interfaces.");
311 2608 : const unsigned int num_total_background_layers =
312 2608 : _background_intervals + _background_inner_boundary_layer_params.intervals +
313 2608 : _background_outer_boundary_layer_params.intervals;
314 2608 : if ((_has_rings || num_total_background_layers == 1) && _background_block_ids.size() > 1)
315 2 : paramError(
316 : "background_block_ids",
317 : "This parameter must be either unset or have a unity length when ring_radii is "
318 : "provided or the number of background intervals (including boundary layers) is unity.");
319 2606 : if ((_has_rings || num_total_background_layers == 1) && _background_block_names.size() > 1)
320 2 : paramError(
321 : "background_block_names",
322 : "This parameter must be either unset or have a unity length when ring_radii is "
323 : "provided or the number of background intervals (including boundary layers) is unity.");
324 2604 : if (!_has_rings && num_total_background_layers > 1 && _quad_center_elements &&
325 : _background_block_ids.size() == 1)
326 18 : _background_block_ids.insert(_background_block_ids.begin(), _background_block_ids.front());
327 2604 : if ((!_has_rings && _background_intervals > 1) &&
328 262 : (!_background_block_ids.empty() && _background_block_ids.size() != 2))
329 2 : paramError("background_block_ids",
330 : "This parameter must be either unset or have a length of two when ring_radii is not "
331 : "provided and background intervals (including boundary layers) is not unity. It can "
332 : "optionally to have a unity length if `quad_center_elements` is enabled.");
333 2602 : if (!_has_rings && num_total_background_layers > 1 && _quad_center_elements &&
334 : _background_block_names.size() == 1)
335 18 : _background_block_names.insert(_background_block_names.begin(),
336 : _background_block_names.front());
337 2602 : if ((!_has_rings && _background_intervals > 1) &&
338 264 : (!_background_block_names.empty() && _background_block_names.size() != 2))
339 2 : paramError("background_block_names",
340 : "This parameter must be either unset or have a length of two when ring_radii is not "
341 : "provided and background intervals (including boundary layers) is not unity. It can "
342 : "optionally have a unity length if `quad_center_elements` is enabled.");
343 2600 : if (_num_sectors_per_side.size() != _num_sides)
344 2 : paramError("num_sectors_per_side",
345 : "This parameter must have a length that is consistent with num_sides.");
346 15899 : for (auto it = _num_sectors_per_side.begin(); it != _num_sectors_per_side.end(); ++it)
347 13303 : if (*it % 2 == 1)
348 2 : paramError("num_sectors_per_side", "This parameter must be even.");
349 2596 : declareMeshProperty("num_sectors_per_side_meta", _num_sectors_per_side);
350 :
351 7788 : if (!getParam<bool>("replace_inner_ring_with_delaunay_mesh") &&
352 7752 : isParamSetByUser("inner_ring_desired_area"))
353 0 : paramError(
354 : "inner_ring_desired_area",
355 : "This parameter should be set only when 'replace_inner_ring_with_delaunay_mesh=true'");
356 :
357 2596 : if (_has_rings)
358 : {
359 : const unsigned int num_innermost_ring_layers =
360 1953 : _ring_inner_boundary_layer_params.intervals.front() + _ring_intervals.front() +
361 1953 : _ring_outer_boundary_layer_params.intervals.front();
362 : // If conditions are met, duplicate the first element of _ring_block_ids at the start.
363 1953 : if (!_ring_block_ids.empty() && _quad_center_elements && num_innermost_ring_layers > 1 &&
364 : _ring_block_ids.size() == _ring_intervals.size())
365 9 : _ring_block_ids.insert(_ring_block_ids.begin(), _ring_block_ids.front());
366 : // check if the number of ring block ids is appropiate
367 1953 : if (!_ring_block_ids.empty() &&
368 : _ring_block_ids.size() !=
369 1482 : (_ring_intervals.size() + (unsigned int)(num_innermost_ring_layers != 1)))
370 : {
371 : // Create an ostringstream for the debug information
372 4 : std::ostringstream debug_info;
373 6 : debug_info << "quad_center_elements is : " << (_quad_center_elements ? "true" : "false")
374 : << std::endl;
375 4 : debug_info << "ring_block_ids size is : " << _ring_block_ids.size() << std::endl;
376 4 : debug_info << "ring_intervals size is : " << _ring_intervals.size() << std::endl;
377 4 : debug_info << "number of innermost ring layers is : " << num_innermost_ring_layers
378 : << std::endl;
379 :
380 : // error message
381 4 : if (!_quad_center_elements)
382 : {
383 4 : paramError(
384 : "ring_block_ids",
385 : "This parameter must have the appropriate size if it is provided: "
386 2 : "Since the number of the innermost ring layers (" +
387 2 : std::to_string(num_innermost_ring_layers) +
388 : " :first ring interval and inner/outer boundaries of the first ring) is more than "
389 : "one, and we have non-quad central elements, the size of 'ring_block_ids' must be "
390 : "equal to the size of 'ring_intervals' + 1.\n",
391 0 : debug_info.str());
392 : }
393 : else
394 : {
395 4 : paramError(
396 : "ring_block_ids",
397 : "This parameter must have the appropriate size if it is provided: "
398 2 : "Since the number of the innermost ring layers (" +
399 2 : std::to_string(num_innermost_ring_layers) +
400 : " :first ring interval and inner/outer boundaries of the first ring) is more than "
401 : "one, and we have quad central elements, the size of 'ring_block_ids' can be either"
402 : "equal to the size of 'ring_intervals' or 'ring_intervals' + 1.\n",
403 0 : debug_info.str());
404 : }
405 0 : }
406 : // If conditions are met, duplicate the first element of _ring_block_names at the start.
407 1949 : if (!_ring_block_names.empty() && _quad_center_elements && num_innermost_ring_layers > 1 &&
408 : _ring_block_names.size() == _ring_intervals.size())
409 9 : _ring_block_names.insert(_ring_block_names.begin(), _ring_block_names.front());
410 : // check if the number of ring block names is appropiate
411 1949 : if (!_ring_block_names.empty() &&
412 : _ring_block_names.size() !=
413 1455 : (_ring_intervals.size() + (unsigned int)(num_innermost_ring_layers != 1)))
414 : {
415 : // Create an ostringstream for the debug information
416 4 : std::ostringstream debug_info;
417 6 : debug_info << "quad_center_elements is : " << (_quad_center_elements ? "true" : "false")
418 : << std::endl;
419 4 : debug_info << "ring_block_names size is : " << _ring_block_names.size() << std::endl;
420 4 : debug_info << "ring_intervals size is : " << _ring_intervals.size() << std::endl;
421 4 : debug_info << "number of innermost ring layers is : " << num_innermost_ring_layers
422 : << std::endl;
423 :
424 : // error message
425 4 : if (!_quad_center_elements)
426 : {
427 4 : paramError(
428 : "ring_block_names",
429 : "This parameter must have the appropriate size if it is provided: "
430 2 : "Since the number of the innermost ring layers (" +
431 2 : std::to_string(num_innermost_ring_layers) +
432 : " :first ring interval and inner/outer boundaries of the first ring) is more than "
433 : "one, and we have non-quad central elements, the size of 'ring_block_names' must "
434 : "be equal to the size of 'ring_intervals' + 1.\n",
435 0 : debug_info.str());
436 : }
437 : else
438 : {
439 4 : paramError(
440 : "ring_block_names",
441 : "This parameter must have the appropriate size if it is provided: "
442 2 : "Since the number of the innermost ring layers (" +
443 2 : std::to_string(num_innermost_ring_layers) +
444 : " :first ring interval and inner/outer boundaries of the first ring) is more than "
445 : "one, and we have quad central elements, the size of 'ring_block_names' can be "
446 : "either equal to the size of 'ring_intervals' or 'ring_intervals' + 1.\n",
447 0 : debug_info.str());
448 : }
449 0 : }
450 4494 : for (unsigned int i = 0; i < _ring_radii.size(); i++)
451 : {
452 2549 : const Real layer_width = _ring_radii[i] - (i == 0 ? 0.0 : _ring_radii[i - 1]);
453 2549 : _ring_inner_boundary_layer_params.fractions.push_back(
454 2549 : _ring_inner_boundary_layer_params.widths[i] / layer_width);
455 2549 : _ring_outer_boundary_layer_params.fractions.push_back(
456 2549 : _ring_outer_boundary_layer_params.widths[i] / layer_width);
457 : }
458 4488 : for (unsigned int i = 0; i < _ring_inner_boundary_layer_params.fractions.size(); i++)
459 2547 : if (MooseUtils::absoluteFuzzyEqual(_ring_inner_boundary_layer_params.fractions[i], 0.0) &&
460 2527 : _ring_inner_boundary_layer_params.intervals[i] > 0)
461 2 : paramError("ring_inner_boundary_layer_intervals",
462 : "Ring inner boundary layer must have zero interval if its thickness is zero.");
463 : else if (MooseUtils::absoluteFuzzyGreaterThan(_ring_inner_boundary_layer_params.fractions[i],
464 2545 : 0.0) &&
465 20 : _ring_inner_boundary_layer_params.intervals[i] == 0)
466 2 : paramError("ring_inner_boundary_layer_intervals",
467 : "Ring inner boundary layer must have non-zero interval if its thickness is "
468 : "not zero.");
469 4472 : for (unsigned int i = 0; i < _ring_outer_boundary_layer_params.fractions.size(); i++)
470 : {
471 2537 : if (MooseUtils::absoluteFuzzyEqual(_ring_outer_boundary_layer_params.fractions[i], 0.0) &&
472 2497 : _ring_outer_boundary_layer_params.intervals[i] > 0)
473 2 : paramError("ring_outer_boundary_layer_intervals",
474 : "Ring outer boundary layer must have zero interval if its thickness is zero.");
475 : else if (MooseUtils::absoluteFuzzyGreaterThan(_ring_outer_boundary_layer_params.fractions[i],
476 2535 : 0.0) &&
477 40 : _ring_outer_boundary_layer_params.intervals[i] == 0)
478 2 : paramError("ring_outer_boundary_layer_intervals",
479 : "Ring outer boundary layer must have non-zero interval if its thickness is "
480 : "not zero.");
481 2533 : if (_ring_inner_boundary_layer_params.fractions[i] +
482 : _ring_outer_boundary_layer_params.fractions[i] >=
483 : 1.0)
484 2 : paramError("ring_inner_boundary_layer_widths",
485 : "Summation of ring_inner_boundary_layer_widths and "
486 : "ring_outer_boundary_layer_widths cannot exceeds the ring layer width.");
487 : }
488 : }
489 : // Ducts related error messages
490 2578 : if (_duct_sizes.size() != _duct_intervals.size())
491 2 : paramError("duct_sizes", "This parameter and duct_intervals must have the same length.");
492 2576 : if (_duct_sizes.size() != _duct_radial_biases.size())
493 2 : paramError("duct_sizes", "This parameter and duct_radial_biases must have the same length.");
494 2574 : if (!_duct_block_ids.empty() && _duct_block_ids.size() != _duct_intervals.size())
495 2 : paramError("duct_block_ids",
496 : "This parameter must have the same length as duct_intervals if set.");
497 2572 : if (!_duct_block_names.empty() && _duct_block_names.size() != _duct_intervals.size())
498 2 : paramError("duct_block_names",
499 : "This parameter must have the same length as duct_intervals if set.");
500 2568 : if (_duct_sizes.size() != _duct_inner_boundary_layer_params.widths.size() ||
501 2568 : _duct_sizes.size() != _duct_inner_boundary_layer_params.intervals.size() ||
502 2568 : _duct_sizes.size() != _duct_inner_boundary_layer_params.biases.size() ||
503 2568 : _duct_sizes.size() != _duct_outer_boundary_layer_params.widths.size() ||
504 5138 : _duct_sizes.size() != _duct_outer_boundary_layer_params.intervals.size() ||
505 : _duct_sizes.size() != _duct_outer_boundary_layer_params.biases.size())
506 2 : paramError("duct_sizes",
507 : "The inner and outer duct boundary layer parameters must have the same sizes as "
508 : "duct_sizes.");
509 2568 : if (_has_ducts)
510 : {
511 364 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
512 786 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
513 438 : _duct_sizes[i] /= std::cos(M_PI / Real(_num_sides));
514 468 : for (unsigned int i = 1; i < _duct_sizes.size(); i++)
515 106 : if (_duct_sizes[i] <= _duct_sizes[i - 1])
516 2 : paramError("duct_sizes", "This parameter must be strictly ascending.");
517 362 : if (_duct_sizes.front() <=
518 362 : (_has_rings ? _ring_radii.back() / std::cos(M_PI / Real(_num_sides)) : 0.0))
519 2 : paramError("duct_sizes",
520 : "This parameter must be positive and ensures no overlapping with rings.");
521 360 : if (_duct_sizes.back() >= _pitch / 2.0 / std::cos(M_PI / Real(_num_sides)))
522 2 : paramError("duct_sizes",
523 : "This parameter must ensure that ducts are smaller than the polygon size.");
524 358 : if (*std::min_element(_duct_intervals.begin(), _duct_intervals.end()) <= 0)
525 0 : paramError("duct_intervals", "Elements of this parameter must be positive.");
526 358 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
527 786 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
528 : {
529 438 : _duct_inner_boundary_layer_params.widths[i] /= std::cos(M_PI / Real(_num_sides));
530 438 : _duct_outer_boundary_layer_params.widths[i] /= std::cos(M_PI / Real(_num_sides));
531 : }
532 816 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
533 : {
534 : const Real layer_width =
535 458 : (i == _duct_sizes.size() - 1 ? _pitch / 2.0 / std::cos(M_PI / Real(_num_sides))
536 100 : : _duct_sizes[i + 1]) -
537 458 : _duct_sizes[i];
538 458 : _duct_inner_boundary_layer_params.fractions.push_back(
539 458 : _duct_inner_boundary_layer_params.widths[i] / layer_width);
540 458 : _duct_outer_boundary_layer_params.fractions.push_back(
541 458 : _duct_outer_boundary_layer_params.widths[i] / layer_width);
542 : }
543 810 : for (unsigned int i = 0; i < _duct_inner_boundary_layer_params.fractions.size(); i++)
544 456 : if (MooseUtils::absoluteFuzzyEqual(_duct_inner_boundary_layer_params.fractions[i], 0.0) &&
545 418 : _duct_inner_boundary_layer_params.intervals[i] > 0)
546 2 : paramError("duct_inner_boundary_layer_intervals",
547 : "Duct inner boundary layer must have zero interval if its thickness is zero.");
548 : else if (MooseUtils::absoluteFuzzyGreaterThan(_duct_inner_boundary_layer_params.fractions[i],
549 454 : 0.0) &&
550 38 : _duct_inner_boundary_layer_params.intervals[i] == 0)
551 2 : paramError("duct_inner_boundary_layer_intervals",
552 : "Duct inner boundary layer must have non-zero interval if its thickness is "
553 : "not zero.");
554 794 : for (unsigned int i = 0; i < _duct_outer_boundary_layer_params.fractions.size(); i++)
555 : {
556 446 : if (MooseUtils::absoluteFuzzyEqual(_duct_outer_boundary_layer_params.fractions[i], 0.0) &&
557 406 : _duct_outer_boundary_layer_params.intervals[i] > 0)
558 2 : paramError("duct_outer_boundary_layer_intervals",
559 : "Duct outer boundary layer must have zero interval if its thickness is zero.");
560 : else if (MooseUtils::absoluteFuzzyGreaterThan(_duct_outer_boundary_layer_params.fractions[i],
561 444 : 0.0) &&
562 40 : _duct_outer_boundary_layer_params.intervals[i] == 0)
563 2 : paramError("duct_outer_boundary_layer_intervals",
564 : "Duct outer boundary layer must have non-zero interval if its thickness is "
565 : "not zero.");
566 442 : if (_duct_inner_boundary_layer_params.fractions[i] +
567 : _duct_outer_boundary_layer_params.fractions[i] >=
568 : 1.0)
569 2 : paramError("duct_inner_boundary_layer_widths",
570 : "Summation of duct_inner_boundary_layer_widths and "
571 : "duct_outer_boundary_layer_widths cannot exceeds the duct layer width.");
572 : }
573 : }
574 2552 : if (MooseUtils::absoluteFuzzyGreaterThan(_background_inner_boundary_layer_params.width +
575 2552 : _background_outer_boundary_layer_params.width,
576 : 0.0))
577 : {
578 : const Real min_background_thickness =
579 29 : (_has_ducts ? (_duct_sizes.front() * std::cos(M_PI / Real(_num_sides))) : (_pitch / 2.0)) -
580 29 : (_has_rings ? _ring_radii.back() : 0.0);
581 29 : if (_background_inner_boundary_layer_params.width +
582 29 : _background_outer_boundary_layer_params.width *
583 29 : (_duct_sizes_style == PolygonSizeStyle::apothem
584 29 : ? 1.0
585 11 : : std::cos(M_PI / Real(_num_sides))) >=
586 : min_background_thickness)
587 2 : paramError("background_inner_boundary_layer_width",
588 : "The summation of background_inner_boundary_layer_width and "
589 : "background_outer_boundary_layer_width must be less than the minimum thickness of "
590 : "the background region.");
591 27 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
592 18 : _background_outer_boundary_layer_params.width /= std::cos(M_PI / Real(_num_sides));
593 : }
594 2550 : if (!_quad_center_elements && _center_quad_factor)
595 2 : paramError("center_quad_factor",
596 : "this parameter is only applicable if quad_center_elements is set true.");
597 2548 : if (_quad_center_elements)
598 602 : declareMeshProperty<subdomain_id_type>(
599 : "quad_center_block_id",
600 1419 : _has_rings ? (_ring_block_ids.empty() ? _block_id_shift + 1 : _ring_block_ids.front())
601 215 : : (_background_block_ids.empty() ? _block_id_shift + 1
602 197 : : _background_block_ids.front()));
603 : else
604 3892 : declareMeshProperty<subdomain_id_type>("quad_center_block_id",
605 : libMesh::Elem::invalid_subdomain_id);
606 :
607 : // declare metadata for internal interface boundaries
608 5096 : declareMeshProperty<bool>("interface_boundaries", false);
609 5096 : declareMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", {});
610 2548 : }
611 :
612 : std::unique_ptr<MeshBase>
613 2376 : PolygonConcentricCircleMeshGeneratorBase::generate()
614 : {
615 2376 : std::vector<std::unique_ptr<ReplicatedMesh>> input(_input_ptrs.size());
616 2760 : for (const auto i : index_range(_input_ptrs))
617 : {
618 772 : input[i] = dynamic_pointer_cast<ReplicatedMesh>(std::move(*_input_ptrs[i]));
619 386 : if (!input[i])
620 0 : mooseError("A non-replicated mesh input was supplied but replicated meshes are required.");
621 1146 : if ((_order == 1 && (*input[i]->elements_begin())->default_order() != FIRST) ||
622 412 : (_order == 2 && (*input[i]->elements_begin())->default_order() == FIRST))
623 2 : paramError("tri_element_type",
624 : "The order of the input mesh to be adapted to does not match the order of the "
625 : "mesh to be generated.");
626 : }
627 :
628 : unsigned int mesh_input_counter = 0;
629 : // azimuthal array used for radius radius_correction
630 : std::vector<Real> azimuthal_list;
631 : // loop over all sides of the polygon to collect azimuthal angles of all nodes
632 14541 : for (unsigned int mesh_index = 0; mesh_index < _num_sides; mesh_index++)
633 : {
634 : // When adaptive boundaries exist (only possible for hexagon meshes thru
635 : // `HexagonConcentricCircleAdaptiveBoundaryMeshGenerator` or square meshes thru
636 : // `CartesianConcentricCircleAdaptiveBoundaryMeshGenerator`), nodes' azimuthal angle cannot be
637 : // arithmetically obtained. Instead, `azimuthalAnglesCollector() is used to get this
638 : // information from the mesh directly.`
639 12167 : if (std::find(_sides_to_adapt.begin(), _sides_to_adapt.end(), mesh_index) !=
640 : _sides_to_adapt.end())
641 : {
642 : // The following lines only work for hexagon and square; and only a hexagon or a square needs
643 : // such functionality.
644 : Real lower_azi =
645 384 : _num_sides == 6 ? ((Real)mesh_index * 60.0 - 150.0) : ((Real)mesh_index * 90.0 - 135.0);
646 384 : Real upper_azi = _num_sides == 6 ? ((Real)((mesh_index + 1) % 6) * 60.0 - 150.0)
647 101 : : ((Real)((mesh_index + 1) % 4) * 90.0 - 135.0);
648 384 : _azimuthal_angles_array.push_back(azimuthalAnglesCollector(
649 384 : *input[mesh_input_counter], lower_azi, upper_azi, ANGLE_TANGENT, _num_sides));
650 : // loop over the _azimuthal_angles_array just collected to convert tangent to azimuthal
651 : // angles.
652 3027 : for (unsigned int i = 1; i < _azimuthal_angles_array.back().size(); i++)
653 : {
654 : azimuthal_list.push_back(
655 2643 : _num_sides == 6
656 4728 : ? ((Real)mesh_index * 60.0 - 150.0 +
657 2085 : std::atan((_azimuthal_angles_array.back()[i - 1] - 1.0) / std::sqrt(3.0)) *
658 2085 : 180.0 / M_PI)
659 558 : : ((Real)mesh_index * 90.0 - 135.0 +
660 558 : std::atan((_azimuthal_angles_array.back()[i - 1] - 1.0)) * 180.0 / M_PI));
661 : }
662 384 : mesh_input_counter++;
663 : }
664 : else
665 : {
666 11783 : _azimuthal_angles_array.push_back(std::vector<Real>());
667 43901 : for (unsigned int i = 0; i < _num_sectors_per_side[mesh_index] * _order; i++)
668 : {
669 : azimuthal_list.push_back(
670 32118 : std::atan(
671 32118 : std::tan(M_PI / _num_sides) *
672 32118 : (2.0 * (Real)i / (Real)_num_sectors_per_side[mesh_index] / (Real)_order - 1.0)) *
673 32118 : 180.0 / M_PI +
674 32118 : (Real)mesh_index * (360.0 / (Real)_num_sides) - (180.0 - 180.0 / (Real)_num_sides));
675 : }
676 : }
677 : }
678 : std::vector<Real> ring_radii_corr;
679 2374 : if (_has_rings)
680 : {
681 1779 : if (_preserve_volumes)
682 : {
683 : Real corr_factor =
684 1779 : PolygonalMeshGenerationUtils::radiusCorrectionFactor(azimuthal_list, true, _order);
685 4112 : for (unsigned int i = 0; i < _ring_radii.size(); i++)
686 2333 : ring_radii_corr.push_back(_ring_radii[i] * corr_factor);
687 : }
688 : else
689 0 : ring_radii_corr = _ring_radii;
690 1779 : if (ring_radii_corr.back() >= _pitch / 2.0)
691 2 : paramError("ring_radii",
692 : "Elements of this parameter must be smaller than polygon apothem (after volume "
693 : "preserve correction if applicable).");
694 3554 : setMeshProperty("max_radius_meta", ring_radii_corr.back());
695 : }
696 : // build the first slice of the polygon.
697 : auto mesh0 = buildSimpleSlice(ring_radii_corr,
698 2372 : _ring_intervals,
699 2372 : _ring_radial_biases,
700 2372 : _ring_inner_boundary_layer_params,
701 2372 : _ring_outer_boundary_layer_params,
702 2372 : _duct_sizes,
703 2372 : _duct_intervals,
704 2372 : _duct_radial_biases,
705 2372 : _duct_inner_boundary_layer_params,
706 2372 : _duct_outer_boundary_layer_params,
707 : _pitch,
708 : _num_sectors_per_side[0],
709 2372 : _background_intervals,
710 2372 : _background_radial_bias,
711 2372 : _background_inner_boundary_layer_params,
712 2372 : _background_outer_boundary_layer_params,
713 : _node_id_background_meta,
714 2372 : _num_sides,
715 : 1,
716 : _azimuthal_angles_array[0],
717 2372 : _block_id_shift,
718 2372 : _quad_center_elements,
719 2372 : _center_quad_factor,
720 2372 : _create_inward_interface_boundaries,
721 2372 : _create_outward_interface_boundaries,
722 2372 : _interface_boundary_id_shift,
723 2372 : _generate_side_specific_boundaries,
724 : _tri_elem_type,
725 2372 : _quad_elem_type);
726 : // This loop builds add-on slices and stitches them to the first slice
727 12157 : for (unsigned int mesh_index = 1; mesh_index < _num_sides; mesh_index++)
728 : {
729 : auto mesh_tmp = buildSimpleSlice(ring_radii_corr,
730 : _ring_intervals,
731 : _ring_radial_biases,
732 : _ring_inner_boundary_layer_params,
733 : _ring_outer_boundary_layer_params,
734 : _duct_sizes,
735 : _duct_intervals,
736 : _duct_radial_biases,
737 : _duct_inner_boundary_layer_params,
738 : _duct_outer_boundary_layer_params,
739 : _pitch,
740 : _num_sectors_per_side[mesh_index],
741 9785 : _background_intervals,
742 9785 : _background_radial_bias,
743 : _background_inner_boundary_layer_params,
744 : _background_outer_boundary_layer_params,
745 : _node_id_background_meta,
746 9785 : _num_sides,
747 : mesh_index + 1,
748 9785 : _azimuthal_angles_array[mesh_index],
749 9785 : _block_id_shift,
750 9785 : _quad_center_elements,
751 9785 : _center_quad_factor,
752 9785 : _create_inward_interface_boundaries,
753 9785 : _create_outward_interface_boundaries,
754 9785 : _interface_boundary_id_shift,
755 9785 : _generate_side_specific_boundaries,
756 : _tri_elem_type,
757 9785 : _quad_elem_type);
758 :
759 9785 : ReplicatedMesh other_mesh(*mesh_tmp);
760 9785 : MeshTools::Modification::rotate(other_mesh, 360.0 / _num_sides * mesh_index, 0, 0);
761 9785 : mesh0->prepare_for_use();
762 9785 : other_mesh.prepare_for_use();
763 9785 : mesh0->stitch_meshes(other_mesh, SLICE_BEGIN, SLICE_END, TOLERANCE, true, false);
764 9785 : other_mesh.clear();
765 9785 : }
766 :
767 : // An extra step to stich the first and last slices together
768 2372 : mesh0->stitch_surfaces(SLICE_BEGIN, SLICE_END, TOLERANCE, true, false);
769 :
770 2372 : if (!_has_rings && !_has_ducts && _background_intervals == 1)
771 256 : MooseMesh::changeBoundaryId(*mesh0, 1 + _interface_boundary_id_shift, OUTER_SIDESET_ID, false);
772 :
773 2372 : if (!_is_general_polygon)
774 : {
775 367 : setMeshProperty("azimuthal_angle_meta",
776 734 : azimuthalAnglesCollector(*mesh0, -180.0, 180.0, ANGLE_DEGREE));
777 734 : setMeshProperty("pattern_pitch_meta", _pitch);
778 : }
779 :
780 : // Move nodes on the external boundary for force uniform spacing.
781 2372 : if (_uniform_mesh_on_sides)
782 : {
783 : const Real angle_tol = 1.0e-5;
784 : std::vector<std::pair<Real, unsigned int>> node_azi_list;
785 36 : MeshTools::Modification::rotate(*mesh0, -(270.0 - 360.0 / (Real)_num_sides), 0.0, 0.0);
786 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> side_list =
787 36 : mesh0->get_boundary_info().build_side_list();
788 36 : mesh0->get_boundary_info().build_node_list_from_side_list();
789 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
790 36 : mesh0->get_boundary_info().build_node_list();
791 :
792 3276 : for (unsigned int i = 0; i < node_list.size(); ++i)
793 : {
794 3240 : if (std::get<1>(node_list[i]) == OUTER_SIDESET_ID)
795 : {
796 720 : node_azi_list.push_back(
797 720 : std::make_pair(atan2((*mesh0->node_ptr(std::get<0>(node_list[i])))(1),
798 720 : (*mesh0->node_ptr(std::get<0>(node_list[i])))(0)) *
799 720 : 180.0 / M_PI,
800 : std::get<0>(node_list[i])));
801 720 : if (std::abs(node_azi_list.back().first + 180.0) <= angle_tol)
802 0 : node_azi_list.back().first = 180.0;
803 : }
804 : }
805 36 : std::sort(node_azi_list.begin(), node_azi_list.end());
806 216 : for (unsigned int i = 0; i < _num_sides; i++)
807 : {
808 900 : for (unsigned int j = 1; j <= _num_sectors_per_side[i]; j++)
809 : {
810 720 : Real azi_corr_tmp = atan2((Real)j * 2.0 / (Real)_num_sectors_per_side[i] - 1.0,
811 720 : 1.0 / std::tan(M_PI / (Real)_num_sides));
812 720 : Real x_tmp = _pitch / 2.0;
813 720 : Real y_tmp = x_tmp * std::tan(azi_corr_tmp);
814 720 : nodeCoordRotate(
815 720 : x_tmp, y_tmp, (Real)i * 360.0 / (Real)_num_sides - (180.0 - 180.0 / (Real)_num_sides));
816 720 : Point p_tmp = Point(x_tmp, y_tmp, 0.0);
817 720 : mesh0->add_point(
818 : p_tmp,
819 720 : node_azi_list[std::accumulate(
820 720 : _num_sectors_per_side.begin(), _num_sectors_per_side.begin() + i, 0) +
821 720 : j - 1]
822 720 : .second);
823 : }
824 : }
825 36 : MeshTools::Modification::rotate(*mesh0, (270.0 - 360.0 / (Real)_num_sides), 0.0, 0.0);
826 36 : }
827 :
828 2372 : setMeshProperty("background_intervals_meta", _background_intervals);
829 :
830 2372 : if (!_has_ducts && _sides_to_adapt.empty())
831 : {
832 1677 : libMesh::LaplaceMeshSmoother lms(*mesh0);
833 1677 : lms.smooth(_smoothing_max_it);
834 : }
835 :
836 : // Set up customized Block Names and/or IDs
837 2372 : unsigned int block_it = 0;
838 2372 : unsigned ring_block_num = 0;
839 : std::vector<subdomain_id_type> block_ids_old;
840 : std::vector<subdomain_id_type> block_ids_new;
841 : std::vector<SubdomainName> block_names;
842 2372 : if (_has_rings)
843 : {
844 1777 : ringBlockIdsNamesPreparer(block_it, ring_block_num, block_ids_old, block_ids_new, block_names);
845 : }
846 : else
847 : {
848 595 : if (_background_intervals > 1)
849 : {
850 312 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
851 376 : block_ids_new.push_back(_background_block_ids.empty() ? block_ids_old.back()
852 : : _background_block_ids.front());
853 : block_names.push_back(_background_block_names.empty()
854 624 : ? (SubdomainName)std::to_string(block_ids_new.back())
855 : : _background_block_names.front());
856 312 : block_it++;
857 : }
858 : }
859 2372 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
860 4744 : block_ids_new.push_back(_background_block_ids.empty() ? block_ids_old.back()
861 : : _background_block_ids.back());
862 : block_names.push_back(_background_block_names.empty()
863 4744 : ? (SubdomainName)std::to_string(block_ids_new.back())
864 : : _background_block_names.back());
865 2372 : block_it++;
866 2372 : if (_has_ducts)
867 : {
868 766 : for (unsigned int i = 0; i < _duct_intervals.size(); i++)
869 : {
870 428 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
871 856 : block_ids_new.push_back(_duct_block_ids.empty() ? block_ids_old.back() : _duct_block_ids[i]);
872 428 : block_names.push_back(_duct_block_names.empty()
873 856 : ? (SubdomainName)std::to_string(block_ids_new.back())
874 : : _duct_block_names[i]);
875 428 : block_it++;
876 : }
877 : }
878 :
879 2372 : assignBlockIdsNames(
880 : *mesh0, block_ids_old, block_ids_new, block_names, "ConcentricCircleMeshGenerator");
881 :
882 2370 : if (_external_boundary_id > 0)
883 977 : MooseMesh::changeBoundaryId(*mesh0, OUTER_SIDESET_ID, _external_boundary_id, false);
884 2370 : if (!_external_boundary_name.empty())
885 : {
886 995 : mesh0->get_boundary_info().sideset_name(
887 995 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
888 995 : _external_boundary_name;
889 995 : mesh0->get_boundary_info().nodeset_name(
890 995 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
891 : _external_boundary_name;
892 : }
893 :
894 2370 : assignInterfaceBoundaryNames(*mesh0);
895 :
896 : // add sector ids
897 4740 : if (isParamValid("sector_id_name"))
898 90 : setSectorExtraIDs(
899 45 : *mesh0, getParam<std::string>("sector_id_name"), _num_sides, _num_sectors_per_side);
900 :
901 : // add ring ids
902 4740 : if (isParamValid("ring_id_name"))
903 162 : setRingExtraIDs(*mesh0,
904 : getParam<std::string>("ring_id_name"),
905 54 : _num_sides,
906 54 : _num_sectors_per_side,
907 : _ring_intervals,
908 108 : getParam<MooseEnum>("ring_id_assign_type") == "ring_wise",
909 54 : _quad_center_elements);
910 :
911 : // add internal side set info to metadata
912 2370 : if (_create_inward_interface_boundaries || _create_outward_interface_boundaries)
913 : {
914 3236 : setMeshProperty("interface_boundaries", true);
915 : // lists boundary ids assigned to interfaces
916 : std::set<boundary_id_type> boundary_ids = mesh0->get_boundary_info().get_boundary_ids();
917 : std::set<boundary_id_type> interface_boundary_ids;
918 1618 : if (_create_outward_interface_boundaries)
919 : {
920 1609 : const unsigned int num_boundary_ids = boundary_ids.size();
921 9843 : for (const auto i : make_range(num_boundary_ids))
922 : {
923 8234 : const unsigned int id = i * 2 + 1 + _interface_boundary_id_shift;
924 8234 : auto it = boundary_ids.find(id);
925 8234 : if (it != boundary_ids.end())
926 : {
927 2036 : boundary_ids.erase(it);
928 2036 : interface_boundary_ids.insert(id);
929 : }
930 : }
931 : }
932 1618 : if (_create_inward_interface_boundaries)
933 : {
934 9 : const unsigned int num_boundary_ids = boundary_ids.size();
935 153 : for (const auto i : make_range(num_boundary_ids))
936 : {
937 144 : const unsigned int id = i * 2 + 2 + _interface_boundary_id_shift;
938 144 : auto it = boundary_ids.find(id);
939 144 : if (it != boundary_ids.end())
940 : {
941 36 : boundary_ids.erase(it);
942 36 : interface_boundary_ids.insert(id);
943 : }
944 : }
945 : }
946 3236 : setMeshProperty("interface_boundary_ids", interface_boundary_ids);
947 : }
948 :
949 2370 : bool flat_side_up = getMeshProperty<bool>("flat_side_up", name());
950 2370 : if (flat_side_up)
951 925 : MeshTools::Modification::rotate(*mesh0, 180.0 / (Real)_num_sides, 0.0, 0.0);
952 : mesh0->set_isnt_prepared();
953 :
954 5924 : if (_has_rings && getParam<bool>("replace_inner_ring_with_delaunay_mesh"))
955 : {
956 36 : if (isParamSetByUser("quad_center_elements"))
957 0 : paramError("quad_center_elements",
958 : "Should not be set because the center elements of the inner ring will be replaced "
959 : "by a Delaunay mesh with 'replace_inner_ring_with_delaunay_mesh=true'");
960 :
961 : // remove elements of the inner ring and create an inner-ring external boundary side set
962 : std::set<Elem *> deleteable_elems;
963 4356 : for (auto & elem : mesh0->element_ptr_range())
964 2160 : if (elem->vertex_average().norm() < ring_radii_corr[0])
965 738 : deleteable_elems.insert(elem);
966 :
967 18 : const boundary_id_type boundary_id = MooseMeshUtils::getBoundaryIDs(*mesh0, {"_foo"}, true)[0];
968 : BoundaryInfo & boundary_info = mesh0->get_boundary_info();
969 738 : for (auto & elem : deleteable_elems)
970 : {
971 720 : unsigned int n_sides = elem->n_sides();
972 3240 : for (unsigned int n = 0; n != n_sides; ++n)
973 : {
974 2520 : Elem * neighbor = elem->neighbor_ptr(n);
975 2520 : if (!neighbor)
976 1080 : continue;
977 :
978 1440 : const unsigned int return_side = neighbor->which_neighbor_am_i(elem);
979 :
980 1440 : if (neighbor->neighbor_ptr(return_side) == elem)
981 : {
982 : neighbor->set_neighbor(return_side, nullptr);
983 1440 : boundary_info.add_side(neighbor, return_side, boundary_id);
984 : }
985 : }
986 :
987 720 : mesh0->delete_elem(elem);
988 : }
989 :
990 : // build the 1d mesh from the new boundary
991 18 : auto poly_mesh = MooseMeshUtils::buildBoundaryMesh(*mesh0, boundary_id);
992 : Real min_side = std::numeric_limits<Real>::max();
993 756 : for (auto & elem : poly_mesh->element_ptr_range())
994 : {
995 360 : Real l = elem->volume();
996 360 : if (l < min_side)
997 : min_side = l;
998 18 : }
999 :
1000 : // triangulate with the 1d mesh
1001 18 : libMesh::Poly2TriTriangulator poly2tri(*poly_mesh);
1002 18 : poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
1003 18 : poly2tri.set_interpolate_boundary_points(0);
1004 : poly2tri.set_refine_boundary_allowed(false);
1005 : poly2tri.set_verify_hole_boundaries(false);
1006 54 : const auto desired_area = getParam<std::string>("inner_ring_desired_area");
1007 18 : if (desired_area != "")
1008 : {
1009 9 : poly2tri.desired_area() = 0;
1010 9 : libMesh::ParsedFunction<Real> area_func{desired_area};
1011 9 : poly2tri.set_desired_area_function(&area_func);
1012 9 : }
1013 : else
1014 9 : poly2tri.desired_area() = min_side * min_side;
1015 18 : poly2tri.minimum_angle() = 0;
1016 18 : poly2tri.smooth_after_generating() = true;
1017 : // poly2tri.elem_type() is TRI3 by default
1018 18 : if (_tri_elem_type == TRI_ELEM_TYPE::TRI6)
1019 0 : poly2tri.elem_type() = libMesh::ElemType::TRI6;
1020 18 : else if (_tri_elem_type == TRI_ELEM_TYPE::TRI7)
1021 0 : poly2tri.elem_type() = libMesh::ElemType::TRI7;
1022 : // let us keep the center point
1023 18 : poly_mesh->add_point(libMesh::Point());
1024 18 : poly2tri.triangulate();
1025 : // keep the old subdomain id
1026 2628 : for (auto elem : poly_mesh->element_ptr_range())
1027 1314 : elem->subdomain_id() = block_ids_new[0];
1028 :
1029 36 : if (isParamValid("ring_id_name"))
1030 : {
1031 27 : if (_ring_intervals[0] != 1 && getParam<MooseEnum>("ring_id_assign_type") == "ring_wise")
1032 0 : paramError("replace_inner_ring_with_delaunay_mesh",
1033 : "Inner ring has multple intervals with each being assigned with a different "
1034 : "ring id, replacing inner ring with Delaunay mesh will remove this ring id "
1035 : "assign type. Either change 'ring_id_assign_type' to block_wise or set the "
1036 : "first element of 'ring_intervals' to 1 or set "
1037 : "'replace_inner_ring_with_delaunay_mesh' to false to avoid this error.");
1038 27 : auto id_name = getParam<std::string>("ring_id_name");
1039 18 : const auto extra_id_index = poly_mesh->add_elem_integer(id_name);
1040 1170 : for (auto elem : poly_mesh->element_ptr_range())
1041 585 : elem->set_extra_integer(extra_id_index, 1);
1042 : }
1043 36 : if (isParamValid("sector_id_name"))
1044 : {
1045 : // we will assign all elements in the inner ring with zero sector id
1046 27 : auto id_name = getParam<std::string>("sector_id_name");
1047 18 : const auto extra_id_index = poly_mesh->add_elem_integer(id_name);
1048 1170 : for (auto elem : poly_mesh->element_ptr_range())
1049 585 : elem->set_extra_integer(extra_id_index, 0);
1050 : }
1051 :
1052 : // stitch the triangulated mesh and the original mesh without the inner ring
1053 18 : mesh0->stitch_meshes(*poly_mesh, boundary_id, 0, TOLERANCE, true, false);
1054 18 : }
1055 :
1056 4740 : return dynamic_pointer_cast<MeshBase>(mesh0);
1057 2370 : }
|
