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 4190 : PolygonConcentricCircleMeshGeneratorBase::validParams()
23 : {
24 4190 : InputParameters params = ConcentricCircleGeneratorBase::validParams();
25 8380 : 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 12570 : params.addRangeCheckedParam<unsigned int>(
31 : "background_intervals",
32 8380 : 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 12570 : params.addRangeCheckedParam<Real>(
37 : "background_radial_bias",
38 8380 : 1.0,
39 : "background_radial_bias>0",
40 : "Value used to create biasing in radial meshing for background region.");
41 12570 : params.addRangeCheckedParam<Real>(
42 : "background_inner_boundary_layer_width",
43 8380 : 0.0,
44 : "background_inner_boundary_layer_width>=0",
45 : "Width of background region that is assigned to be the inner boundary layer.");
46 12570 : params.addRangeCheckedParam<unsigned int>(
47 : "background_inner_boundary_layer_intervals",
48 8380 : 1,
49 : "background_inner_boundary_layer_intervals>0",
50 : "Number of radial intervals of the background inner boundary layer");
51 12570 : params.addRangeCheckedParam<Real>(
52 : "background_inner_boundary_layer_bias",
53 8380 : 1.0,
54 : "background_inner_boundary_layer_bias>0",
55 : "Growth factor used for mesh biasing of the background inner boundary layer.");
56 12570 : params.addRangeCheckedParam<Real>(
57 : "background_outer_boundary_layer_width",
58 8380 : 0.0,
59 : "background_outer_boundary_layer_width>=0",
60 : "Width of background region that is assigned to be the outer boundary layer.");
61 12570 : params.addRangeCheckedParam<unsigned int>(
62 : "background_outer_boundary_layer_intervals",
63 8380 : 1,
64 : "background_outer_boundary_layer_intervals>0",
65 : "Number of radial intervals of the background outer boundary layer");
66 12570 : params.addRangeCheckedParam<Real>(
67 : "background_outer_boundary_layer_bias",
68 8380 : 1.0,
69 : "background_outer_boundary_layer_bias>0",
70 : "Growth factor used for mesh biasing of the background outer boundary layer.");
71 8380 : params.addParam<std::vector<subdomain_id_type>>(
72 : "background_block_ids", "Optional customized block id for the background block.");
73 8380 : params.addParam<std::vector<SubdomainName>>(
74 : "background_block_names", "Optional customized block names for the background block.");
75 8380 : params.addParam<std::vector<Real>>(
76 : "duct_sizes", "Distance(s) from polygon center to duct(s) inner boundaries.");
77 8380 : MooseEnum duct_sizes_style("apothem radius", "radius");
78 8380 : params.addParam<MooseEnum>(
79 : "duct_sizes_style",
80 : duct_sizes_style,
81 : "Style in which polygon center to duct inner boundary distance is "
82 4190 : "given (apothem = center to face, radius = center to vertex). Options: " +
83 4190 : duct_sizes_style.getRawNames());
84 8380 : 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 8380 : 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 8380 : 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 8380 : params.addParam<std::vector<unsigned int>>(
97 : "duct_inner_boundary_layer_intervals",
98 : "Number of radial intervals of the duct inner boundary layers");
99 8380 : 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 8380 : 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 8380 : params.addParam<std::vector<unsigned int>>(
108 : "duct_outer_boundary_layer_intervals",
109 : "Number of radial intervals of the duct outer boundary layers");
110 8380 : 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 8380 : params.addParam<std::vector<subdomain_id_type>>(
115 : "duct_block_ids", "Optional customized block ids for each duct geometry block.");
116 8380 : params.addParam<std::vector<SubdomainName>>(
117 : "duct_block_names", "Optional customized block names for each duct geometry block.");
118 8380 : params.addParam<bool>("uniform_mesh_on_sides",
119 8380 : false,
120 : "Whether the side elements are reorganized to have a uniform size.");
121 8380 : params.addParam<unsigned int>("smoothing_max_it",
122 8380 : 0,
123 : "Number of Laplacian smoothing iterations. This number is "
124 : "disregarded when duct_sizes is present.");
125 8380 : params.addParam<bool>(
126 : "flat_side_up",
127 8380 : false,
128 : "Whether to rotate the generated polygon mesh to ensure that one flat side faces up.");
129 :
130 8380 : params.addParam<bool>(
131 8380 : "quad_center_elements", false, "Whether the center elements are quad or triangular.");
132 8380 : 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 8380 : params.addParam<bool>("replace_inner_ring_with_delaunay_mesh",
138 8380 : false,
139 : "True to replace the inner ring mesh with a Delaunay unstructured mesh");
140 8380 : params.addParam<std::string>(
141 : "inner_ring_desired_area",
142 4190 : 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 8380 : params.addParamNamesToGroup(
147 : "background_block_ids background_block_names duct_block_ids duct_block_names",
148 : "Customized Subdomain/Boundary");
149 8380 : params.addParamNamesToGroup("num_sectors_per_side background_intervals duct_intervals "
150 : "uniform_mesh_on_sides",
151 : "General Mesh Density");
152 8380 : 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 8380 : params.addParamNamesToGroup("quad_center_elements center_quad_factor "
158 : "replace_inner_ring_with_delaunay_mesh inner_ring_desired_area ",
159 : "Inner ring");
160 :
161 4190 : addRingAndSectorIDParams(params);
162 4190 : params.addClassDescription("This PolygonConcentricCircleMeshGeneratorBase object is a base class "
163 : "to be inherited for polygon mesh generators.");
164 :
165 4190 : return params;
166 4190 : }
167 :
168 2138 : PolygonConcentricCircleMeshGeneratorBase::PolygonConcentricCircleMeshGeneratorBase(
169 2138 : const InputParameters & parameters)
170 : : ConcentricCircleGeneratorBase(parameters),
171 6372 : _num_sides(isParamValid("num_sides")
172 6043 : ? getParam<unsigned int>("num_sides")
173 2865 : : (isParamValid("hexagon_size") ? (unsigned int)HEXAGON_NUM_SIDES
174 : : (unsigned int)SQUARE_NUM_SIDES)),
175 4248 : _duct_sizes_style(getParam<MooseEnum>("duct_sizes_style").template getEnum<PolygonSizeStyle>()),
176 6672 : _duct_sizes(isParamValid("duct_sizes") ? getParam<std::vector<Real>>("duct_sizes")
177 : : std::vector<Real>()),
178 4848 : _duct_intervals(isParamValid("duct_intervals")
179 2124 : ? getParam<std::vector<unsigned int>>("duct_intervals")
180 : : std::vector<unsigned int>()),
181 4266 : _duct_radial_biases(isParamValid("duct_radial_biases")
182 2124 : ? getParam<std::vector<Real>>("duct_radial_biases")
183 2124 : : std::vector<Real>(_duct_intervals.size(), 1.0)),
184 6472 : _duct_inner_boundary_layer_params(
185 2124 : {isParamValid("duct_inner_boundary_layer_widths")
186 2124 : ? getParam<std::vector<Real>>("duct_inner_boundary_layer_widths")
187 2104 : : std::vector<Real>(_duct_intervals.size(), 0.0),
188 : std::vector<Real>(),
189 4248 : isParamValid("duct_inner_boundary_layer_intervals")
190 2124 : ? getParam<std::vector<unsigned int>>("duct_inner_boundary_layer_intervals")
191 2108 : : std::vector<unsigned int>(_duct_intervals.size(), 0),
192 4248 : isParamValid("duct_inner_boundary_layer_biases")
193 2124 : ? getParam<std::vector<Real>>("duct_inner_boundary_layer_biases")
194 2124 : : std::vector<Real>(_duct_intervals.size(), 0.0)}),
195 6476 : _duct_outer_boundary_layer_params(
196 2124 : {isParamValid("duct_outer_boundary_layer_widths")
197 2124 : ? getParam<std::vector<Real>>("duct_outer_boundary_layer_widths")
198 2104 : : std::vector<Real>(_duct_intervals.size(), 0.0),
199 : std::vector<Real>(),
200 4248 : isParamValid("duct_outer_boundary_layer_intervals")
201 2124 : ? getParam<std::vector<unsigned int>>("duct_outer_boundary_layer_intervals")
202 2106 : : std::vector<unsigned int>(_duct_intervals.size(), 0),
203 4248 : isParamValid("duct_outer_boundary_layer_biases")
204 2124 : ? getParam<std::vector<Real>>("duct_outer_boundary_layer_biases")
205 2124 : : std::vector<Real>(_duct_intervals.size(), 0.0)}),
206 4800 : _duct_block_ids(isParamValid("duct_block_ids")
207 2124 : ? getParam<std::vector<subdomain_id_type>>("duct_block_ids")
208 : : std::vector<subdomain_id_type>()),
209 4786 : _duct_block_names(isParamValid("duct_block_names")
210 2124 : ? getParam<std::vector<SubdomainName>>("duct_block_names")
211 : : std::vector<SubdomainName>()),
212 4248 : _has_rings(isParamValid("ring_radii")),
213 4248 : _has_ducts(isParamValid("duct_sizes")),
214 2124 : _polygon_size_style(
215 2124 : isParamValid("polygon_size_style")
216 4248 : ? getParam<MooseEnum>("polygon_size_style").template getEnum<PolygonSizeStyle>()
217 2453 : : (isParamValid("hexagon_size_style")
218 2453 : ? getParam<MooseEnum>("hexagon_size_style").template getEnum<PolygonSizeStyle>()
219 2290 : : getParam<MooseEnum>("square_size_style")
220 : .template getEnum<PolygonSizeStyle>())),
221 2124 : _polygon_size(isParamValid("polygon_size")
222 6043 : ? getParam<Real>("polygon_size")
223 3357 : : (isParamValid("hexagon_size") ? getParam<Real>("hexagon_size")
224 2290 : : (getParam<Real>("square_size") /
225 : 2.0))), // square size is twice as apothem
226 4248 : _num_sectors_per_side(getParam<std::vector<unsigned int>>("num_sectors_per_side")),
227 4248 : _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 4248 : _background_radial_bias(getParam<Real>("background_radial_bias")),
231 6372 : _background_inner_boundary_layer_params(
232 2124 : {getParam<Real>("background_inner_boundary_layer_width"),
233 : 0.0,
234 4248 : getParam<Real>("background_inner_boundary_layer_width") > 0.0
235 4271 : ? getParam<unsigned int>("background_inner_boundary_layer_intervals")
236 : : 0,
237 4248 : getParam<Real>("background_inner_boundary_layer_bias")}),
238 6372 : _background_outer_boundary_layer_params(
239 2124 : {getParam<Real>("background_outer_boundary_layer_width"),
240 : 0.0,
241 4248 : getParam<Real>("background_outer_boundary_layer_width") > 0.0
242 4271 : ? getParam<unsigned int>("background_outer_boundary_layer_intervals")
243 : : 0,
244 4248 : getParam<Real>("background_outer_boundary_layer_bias")}),
245 7578 : _background_block_ids(isParamValid("background_block_ids")
246 2124 : ? getParam<std::vector<subdomain_id_type>>("background_block_ids")
247 : : std::vector<subdomain_id_type>()),
248 7208 : _background_block_names(isParamValid("background_block_names")
249 2124 : ? getParam<std::vector<SubdomainName>>("background_block_names")
250 : : std::vector<SubdomainName>()),
251 4248 : _uniform_mesh_on_sides(getParam<bool>("uniform_mesh_on_sides")),
252 4248 : _quad_center_elements(getParam<bool>("quad_center_elements")),
253 4266 : _center_quad_factor(isParamValid("center_quad_factor") ? getParam<Real>("center_quad_factor")
254 : : 0.0),
255 4248 : _smoothing_max_it(getParam<unsigned int>("smoothing_max_it")),
256 4886 : _sides_to_adapt(isParamValid("sides_to_adapt")
257 2124 : ? getParam<std::vector<unsigned int>>("sides_to_adapt")
258 : : std::vector<unsigned int>()),
259 2124 : _node_id_background_meta(declareMeshProperty<dof_id_type>("node_id_background_meta", 0)),
260 4262 : _is_control_drum_meta(declareMeshProperty<bool>("is_control_drum_meta", false))
261 : {
262 6372 : declareMeshProperty<bool>("flat_side_up", getParam<bool>("flat_side_up"));
263 2124 : declareMeshProperty<unsigned int>("background_intervals_meta", 0);
264 4248 : declareMeshProperty<Real>("pattern_pitch_meta", 0.0);
265 4248 : declareMeshProperty<std::vector<Real>>("azimuthal_angle_meta", std::vector<Real>());
266 2124 : declareMeshProperty<Real>("max_radius_meta", 0.0);
267 :
268 2124 : const unsigned short tri_order = _tri_elem_type == TRI_ELEM_TYPE::TRI3 ? 1 : 2;
269 2124 : 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 2124 : _order = quad_order;
275 2124 : if (_quad_center_elements)
276 : {
277 482 : if (tri_order != quad_order)
278 14 : _tri_elem_type = quad_order == 1 ? TRI_ELEM_TYPE::TRI3 : TRI_ELEM_TYPE::TRI6;
279 : }
280 370 : else if (_ring_radii.empty() && _background_intervals == 1 &&
281 201 : _background_inner_boundary_layer_params.intervals == 0 &&
282 1843 : _background_outer_boundary_layer_params.intervals == 0 && _duct_sizes.empty())
283 : {
284 194 : _order = tri_order;
285 194 : if (tri_order != quad_order)
286 14 : _quad_elem_type = tri_order == 1 ? QUAD_ELEM_TYPE::QUAD4 : QUAD_ELEM_TYPE::QUAD9;
287 : }
288 1448 : 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 2122 : 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 4244 : _pitch = 2.0 * (_polygon_size_style == PolygonSizeStyle::apothem
298 2122 : ? _polygon_size
299 0 : : _polygon_size * std::cos(M_PI / Real(_num_sides)));
300 4244 : declareMeshProperty<Real>("pitch_meta", _pitch);
301 2122 : 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 2122 : if (_outward_interface_boundary_names.size() > 0 &&
307 30 : _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 2120 : const unsigned int num_total_background_layers =
312 2120 : _background_intervals + _background_inner_boundary_layer_params.intervals +
313 2120 : _background_outer_boundary_layer_params.intervals;
314 2120 : 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 2118 : 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 2116 : if (!_has_rings && num_total_background_layers > 1 && _quad_center_elements &&
325 : _background_block_ids.size() == 1)
326 14 : _background_block_ids.insert(_background_block_ids.begin(), _background_block_ids.front());
327 2116 : if ((!_has_rings && _background_intervals > 1) &&
328 208 : (!_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 2114 : if (!_has_rings && num_total_background_layers > 1 && _quad_center_elements &&
334 : _background_block_names.size() == 1)
335 14 : _background_block_names.insert(_background_block_names.begin(),
336 : _background_block_names.front());
337 2114 : if ((!_has_rings && _background_intervals > 1) &&
338 210 : (!_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 2112 : 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 12905 : for (auto it = _num_sectors_per_side.begin(); it != _num_sectors_per_side.end(); ++it)
347 10797 : if (*it % 2 == 1)
348 2 : paramError("num_sectors_per_side", "This parameter must be even.");
349 2108 : declareMeshProperty("num_sectors_per_side_meta", _num_sectors_per_side);
350 :
351 6324 : if (!getParam<bool>("replace_inner_ring_with_delaunay_mesh") &&
352 6296 : 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 2108 : if (_has_rings)
358 : {
359 : const unsigned int num_innermost_ring_layers =
360 1579 : _ring_inner_boundary_layer_params.intervals.front() + _ring_intervals.front() +
361 1579 : _ring_outer_boundary_layer_params.intervals.front();
362 : // If conditions are met, duplicate the first element of _ring_block_ids at the start.
363 1579 : if (!_ring_block_ids.empty() && _quad_center_elements && num_innermost_ring_layers > 1 &&
364 : _ring_block_ids.size() == _ring_intervals.size())
365 7 : _ring_block_ids.insert(_ring_block_ids.begin(), _ring_block_ids.front());
366 : // check if the number of ring block ids is appropiate
367 1579 : if (!_ring_block_ids.empty() &&
368 : _ring_block_ids.size() !=
369 1218 : (_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 1575 : if (!_ring_block_names.empty() && _quad_center_elements && num_innermost_ring_layers > 1 &&
408 : _ring_block_names.size() == _ring_intervals.size())
409 7 : _ring_block_names.insert(_ring_block_names.begin(), _ring_block_names.front());
410 : // check if the number of ring block names is appropiate
411 1575 : if (!_ring_block_names.empty() &&
412 : _ring_block_names.size() !=
413 1197 : (_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 3610 : for (unsigned int i = 0; i < _ring_radii.size(); i++)
451 : {
452 2039 : const Real layer_width = _ring_radii[i] - (i == 0 ? 0.0 : _ring_radii[i - 1]);
453 2039 : _ring_inner_boundary_layer_params.fractions.push_back(
454 2039 : _ring_inner_boundary_layer_params.widths[i] / layer_width);
455 2039 : _ring_outer_boundary_layer_params.fractions.push_back(
456 2039 : _ring_outer_boundary_layer_params.widths[i] / layer_width);
457 : }
458 3604 : for (unsigned int i = 0; i < _ring_inner_boundary_layer_params.fractions.size(); i++)
459 2037 : if (MooseUtils::absoluteFuzzyEqual(_ring_inner_boundary_layer_params.fractions[i], 0.0) &&
460 2021 : _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 2035 : 0.0) &&
465 16 : _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 3588 : for (unsigned int i = 0; i < _ring_outer_boundary_layer_params.fractions.size(); i++)
470 : {
471 2027 : if (MooseUtils::absoluteFuzzyEqual(_ring_outer_boundary_layer_params.fractions[i], 0.0) &&
472 1995 : _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 2025 : 0.0) &&
477 32 : _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 2023 : 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 2090 : if (_duct_sizes.size() != _duct_intervals.size())
491 2 : paramError("duct_sizes", "This parameter and duct_intervals must have the same length.");
492 2088 : 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 2086 : 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 2084 : 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 2080 : if (_duct_sizes.size() != _duct_inner_boundary_layer_params.widths.size() ||
501 2080 : _duct_sizes.size() != _duct_inner_boundary_layer_params.intervals.size() ||
502 2080 : _duct_sizes.size() != _duct_inner_boundary_layer_params.biases.size() ||
503 2080 : _duct_sizes.size() != _duct_outer_boundary_layer_params.widths.size() ||
504 4162 : _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 2080 : if (_has_ducts)
510 : {
511 290 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
512 618 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
513 344 : _duct_sizes[i] /= std::cos(M_PI / Real(_num_sides));
514 374 : for (unsigned int i = 1; i < _duct_sizes.size(); i++)
515 86 : if (_duct_sizes[i] <= _duct_sizes[i - 1])
516 2 : paramError("duct_sizes", "This parameter must be strictly ascending.");
517 288 : if (_duct_sizes.front() <=
518 288 : (_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 286 : 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 284 : if (*std::min_element(_duct_intervals.begin(), _duct_intervals.end()) <= 0)
525 0 : paramError("duct_intervals", "Elements of this parameter must be positive.");
526 284 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
527 618 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
528 : {
529 344 : _duct_inner_boundary_layer_params.widths[i] /= std::cos(M_PI / Real(_num_sides));
530 344 : _duct_outer_boundary_layer_params.widths[i] /= std::cos(M_PI / Real(_num_sides));
531 : }
532 648 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
533 : {
534 : const Real layer_width =
535 364 : (i == _duct_sizes.size() - 1 ? _pitch / 2.0 / std::cos(M_PI / Real(_num_sides))
536 80 : : _duct_sizes[i + 1]) -
537 364 : _duct_sizes[i];
538 364 : _duct_inner_boundary_layer_params.fractions.push_back(
539 364 : _duct_inner_boundary_layer_params.widths[i] / layer_width);
540 364 : _duct_outer_boundary_layer_params.fractions.push_back(
541 364 : _duct_outer_boundary_layer_params.widths[i] / layer_width);
542 : }
543 642 : for (unsigned int i = 0; i < _duct_inner_boundary_layer_params.fractions.size(); i++)
544 362 : if (MooseUtils::absoluteFuzzyEqual(_duct_inner_boundary_layer_params.fractions[i], 0.0) &&
545 332 : _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 360 : 0.0) &&
550 30 : _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 626 : for (unsigned int i = 0; i < _duct_outer_boundary_layer_params.fractions.size(); i++)
555 : {
556 352 : if (MooseUtils::absoluteFuzzyEqual(_duct_outer_boundary_layer_params.fractions[i], 0.0) &&
557 320 : _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 350 : 0.0) &&
562 32 : _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 348 : 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 2064 : if (MooseUtils::absoluteFuzzyGreaterThan(_background_inner_boundary_layer_params.width +
575 2064 : _background_outer_boundary_layer_params.width,
576 : 0.0))
577 : {
578 : const Real min_background_thickness =
579 23 : (_has_ducts ? (_duct_sizes.front() * std::cos(M_PI / Real(_num_sides))) : (_pitch / 2.0)) -
580 23 : (_has_rings ? _ring_radii.back() : 0.0);
581 23 : if (_background_inner_boundary_layer_params.width +
582 23 : _background_outer_boundary_layer_params.width *
583 23 : (_duct_sizes_style == PolygonSizeStyle::apothem
584 23 : ? 1.0
585 9 : : 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 21 : if (_duct_sizes_style == PolygonSizeStyle::apothem)
592 14 : _background_outer_boundary_layer_params.width /= std::cos(M_PI / Real(_num_sides));
593 : }
594 2062 : 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 2060 : if (_quad_center_elements)
598 478 : declareMeshProperty<subdomain_id_type>(
599 : "quad_center_block_id",
600 1127 : _has_rings ? (_ring_block_ids.empty() ? _block_id_shift + 1 : _ring_block_ids.front())
601 171 : : (_background_block_ids.empty() ? _block_id_shift + 1
602 157 : : _background_block_ids.front()));
603 : else
604 3164 : declareMeshProperty<subdomain_id_type>("quad_center_block_id",
605 : libMesh::Elem::invalid_subdomain_id);
606 :
607 : // declare metadata for internal interface boundaries
608 4120 : declareMeshProperty<bool>("interface_boundaries", false);
609 4120 : declareMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", {});
610 2060 : }
611 :
612 : std::unique_ptr<MeshBase>
613 1916 : PolygonConcentricCircleMeshGeneratorBase::generate()
614 : {
615 1916 : std::vector<std::unique_ptr<ReplicatedMesh>> input(_input_ptrs.size());
616 2236 : for (const auto i : index_range(_input_ptrs))
617 : {
618 644 : input[i] = dynamic_pointer_cast<ReplicatedMesh>(std::move(*_input_ptrs[i]));
619 322 : if (!input[i])
620 0 : mooseError("A non-replicated mesh input was supplied but replicated meshes are required.");
621 954 : if ((_order == 1 && (*input[i]->elements_begin())->default_order() != FIRST) ||
622 348 : (_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 11743 : 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 9829 : 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 320 : _num_sides == 6 ? ((Real)mesh_index * 60.0 - 150.0) : ((Real)mesh_index * 90.0 - 135.0);
646 320 : Real upper_azi = _num_sides == 6 ? ((Real)((mesh_index + 1) % 6) * 60.0 - 150.0)
647 83 : : ((Real)((mesh_index + 1) % 4) * 90.0 - 135.0);
648 320 : _azimuthal_angles_array.push_back(azimuthalAnglesCollector(
649 320 : *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 2527 : for (unsigned int i = 1; i < _azimuthal_angles_array.back().size(); i++)
653 : {
654 : azimuthal_list.push_back(
655 2207 : _num_sides == 6
656 3940 : ? ((Real)mesh_index * 60.0 - 150.0 +
657 1733 : std::atan((_azimuthal_angles_array.back()[i - 1] - 1.0) / std::sqrt(3.0)) *
658 1733 : 180.0 / M_PI)
659 474 : : ((Real)mesh_index * 90.0 - 135.0 +
660 474 : std::atan((_azimuthal_angles_array.back()[i - 1] - 1.0)) * 180.0 / M_PI));
661 : }
662 320 : mesh_input_counter++;
663 : }
664 : else
665 : {
666 9509 : _azimuthal_angles_array.push_back(std::vector<Real>());
667 35447 : for (unsigned int i = 0; i < _num_sectors_per_side[mesh_index] * _order; i++)
668 : {
669 : azimuthal_list.push_back(
670 25938 : std::atan(
671 25938 : std::tan(M_PI / _num_sides) *
672 25938 : (2.0 * (Real)i / (Real)_num_sectors_per_side[mesh_index] / (Real)_order - 1.0)) *
673 25938 : 180.0 / M_PI +
674 25938 : (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 1914 : if (_has_rings)
680 : {
681 1433 : if (_preserve_volumes)
682 : {
683 : Real corr_factor =
684 1433 : PolygonalMeshGenerationUtils::radiusCorrectionFactor(azimuthal_list, true, _order);
685 3308 : for (unsigned int i = 0; i < _ring_radii.size(); i++)
686 1875 : ring_radii_corr.push_back(_ring_radii[i] * corr_factor);
687 : }
688 : else
689 0 : ring_radii_corr = _ring_radii;
690 1433 : 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 2862 : 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 1912 : _ring_intervals,
699 1912 : _ring_radial_biases,
700 1912 : _ring_inner_boundary_layer_params,
701 1912 : _ring_outer_boundary_layer_params,
702 1912 : _duct_sizes,
703 1912 : _duct_intervals,
704 1912 : _duct_radial_biases,
705 1912 : _duct_inner_boundary_layer_params,
706 1912 : _duct_outer_boundary_layer_params,
707 : _pitch,
708 : _num_sectors_per_side[0],
709 1912 : _background_intervals,
710 1912 : _background_radial_bias,
711 1912 : _background_inner_boundary_layer_params,
712 1912 : _background_outer_boundary_layer_params,
713 : _node_id_background_meta,
714 1912 : _num_sides,
715 : 1,
716 : _azimuthal_angles_array[0],
717 1912 : _block_id_shift,
718 1912 : _quad_center_elements,
719 1912 : _center_quad_factor,
720 1912 : _create_inward_interface_boundaries,
721 1912 : _create_outward_interface_boundaries,
722 1912 : _interface_boundary_id_shift,
723 1912 : _generate_side_specific_boundaries,
724 : _tri_elem_type,
725 1912 : _quad_elem_type);
726 : // This loop builds add-on slices and stitches them to the first slice
727 9819 : 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 7907 : _background_intervals,
742 7907 : _background_radial_bias,
743 : _background_inner_boundary_layer_params,
744 : _background_outer_boundary_layer_params,
745 : _node_id_background_meta,
746 7907 : _num_sides,
747 : mesh_index + 1,
748 7907 : _azimuthal_angles_array[mesh_index],
749 7907 : _block_id_shift,
750 7907 : _quad_center_elements,
751 7907 : _center_quad_factor,
752 7907 : _create_inward_interface_boundaries,
753 7907 : _create_outward_interface_boundaries,
754 7907 : _interface_boundary_id_shift,
755 7907 : _generate_side_specific_boundaries,
756 : _tri_elem_type,
757 7907 : _quad_elem_type);
758 :
759 7907 : ReplicatedMesh other_mesh(*mesh_tmp);
760 7907 : MeshTools::Modification::rotate(other_mesh, 360.0 / _num_sides * mesh_index, 0, 0);
761 7907 : mesh0->prepare_for_use();
762 7907 : other_mesh.prepare_for_use();
763 7907 : mesh0->stitch_meshes(other_mesh, SLICE_BEGIN, SLICE_END, TOLERANCE, true, false);
764 7907 : other_mesh.clear();
765 7907 : }
766 :
767 : // An extra step to stich the first and last slices together
768 1912 : mesh0->stitch_surfaces(SLICE_BEGIN, SLICE_END, TOLERANCE, true, false);
769 :
770 1912 : if (!_has_rings && !_has_ducts && _background_intervals == 1)
771 208 : MooseMesh::changeBoundaryId(*mesh0, 1 + _interface_boundary_id_shift, OUTER_SIDESET_ID, false);
772 :
773 1912 : if (!_is_general_polygon)
774 : {
775 309 : setMeshProperty("azimuthal_angle_meta",
776 618 : azimuthalAnglesCollector(*mesh0, -180.0, 180.0, ANGLE_DEGREE));
777 618 : setMeshProperty("pattern_pitch_meta", _pitch);
778 : }
779 :
780 : // Move nodes on the external boundary for force uniform spacing.
781 1912 : 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 28 : 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 28 : mesh0->get_boundary_info().build_side_list();
788 56 : mesh0->get_boundary_info().build_node_list_from_side_list();
789 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
790 28 : mesh0->get_boundary_info().build_node_list();
791 :
792 2548 : for (unsigned int i = 0; i < node_list.size(); ++i)
793 : {
794 2520 : if (std::get<1>(node_list[i]) == OUTER_SIDESET_ID)
795 : {
796 560 : node_azi_list.push_back(
797 560 : std::make_pair(atan2((*mesh0->node_ptr(std::get<0>(node_list[i])))(1),
798 560 : (*mesh0->node_ptr(std::get<0>(node_list[i])))(0)) *
799 560 : 180.0 / M_PI,
800 : std::get<0>(node_list[i])));
801 560 : if (std::abs(node_azi_list.back().first + 180.0) <= angle_tol)
802 0 : node_azi_list.back().first = 180.0;
803 : }
804 : }
805 28 : std::sort(node_azi_list.begin(), node_azi_list.end());
806 168 : for (unsigned int i = 0; i < _num_sides; i++)
807 : {
808 700 : for (unsigned int j = 1; j <= _num_sectors_per_side[i]; j++)
809 : {
810 560 : Real azi_corr_tmp = atan2((Real)j * 2.0 / (Real)_num_sectors_per_side[i] - 1.0,
811 560 : 1.0 / std::tan(M_PI / (Real)_num_sides));
812 560 : Real x_tmp = _pitch / 2.0;
813 560 : Real y_tmp = x_tmp * std::tan(azi_corr_tmp);
814 560 : nodeCoordRotate(
815 560 : x_tmp, y_tmp, (Real)i * 360.0 / (Real)_num_sides - (180.0 - 180.0 / (Real)_num_sides));
816 560 : Point p_tmp = Point(x_tmp, y_tmp, 0.0);
817 560 : mesh0->add_point(
818 : p_tmp,
819 560 : node_azi_list[std::accumulate(
820 560 : _num_sectors_per_side.begin(), _num_sectors_per_side.begin() + i, 0) +
821 560 : j - 1]
822 560 : .second);
823 : }
824 : }
825 28 : MeshTools::Modification::rotate(*mesh0, (270.0 - 360.0 / (Real)_num_sides), 0.0, 0.0);
826 28 : }
827 :
828 1912 : setMeshProperty("background_intervals_meta", _background_intervals);
829 :
830 1912 : if (!_has_ducts && _sides_to_adapt.empty())
831 : {
832 1349 : libMesh::LaplaceMeshSmoother lms(*mesh0, _smoothing_max_it);
833 1349 : lms.smooth();
834 : }
835 :
836 : // Set up customized Block Names and/or IDs
837 1912 : unsigned int block_it = 0;
838 1912 : 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 1912 : if (_has_rings)
843 : {
844 1431 : ringBlockIdsNamesPreparer(block_it, ring_block_num, block_ids_old, block_ids_new, block_names);
845 : }
846 : else
847 : {
848 481 : if (_background_intervals > 1)
849 : {
850 252 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
851 310 : 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 504 : ? (SubdomainName)std::to_string(block_ids_new.back())
855 : : _background_block_names.front());
856 252 : block_it++;
857 : }
858 : }
859 1912 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
860 3824 : 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 3824 : ? (SubdomainName)std::to_string(block_ids_new.back())
864 : : _background_block_names.back());
865 1912 : block_it++;
866 1912 : if (_has_ducts)
867 : {
868 598 : for (unsigned int i = 0; i < _duct_intervals.size(); i++)
869 : {
870 334 : block_ids_old.push_back(_block_id_shift + 1 + block_it);
871 668 : block_ids_new.push_back(_duct_block_ids.empty() ? block_ids_old.back() : _duct_block_ids[i]);
872 334 : block_names.push_back(_duct_block_names.empty()
873 668 : ? (SubdomainName)std::to_string(block_ids_new.back())
874 : : _duct_block_names[i]);
875 334 : block_it++;
876 : }
877 : }
878 :
879 1912 : assignBlockIdsNames(
880 : *mesh0, block_ids_old, block_ids_new, block_names, "ConcentricCircleMeshGenerator");
881 :
882 1910 : if (_external_boundary_id > 0)
883 781 : MooseMesh::changeBoundaryId(*mesh0, OUTER_SIDESET_ID, _external_boundary_id, false);
884 1910 : if (!_external_boundary_name.empty())
885 : {
886 795 : mesh0->get_boundary_info().sideset_name(
887 795 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
888 795 : _external_boundary_name;
889 795 : mesh0->get_boundary_info().nodeset_name(
890 795 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
891 : _external_boundary_name;
892 : }
893 :
894 1910 : assignInterfaceBoundaryNames(*mesh0);
895 :
896 : // add sector ids
897 3820 : if (isParamValid("sector_id_name"))
898 70 : setSectorExtraIDs(
899 35 : *mesh0, getParam<std::string>("sector_id_name"), _num_sides, _num_sectors_per_side);
900 :
901 : // add ring ids
902 3820 : if (isParamValid("ring_id_name"))
903 126 : setRingExtraIDs(*mesh0,
904 : getParam<std::string>("ring_id_name"),
905 42 : _num_sides,
906 42 : _num_sectors_per_side,
907 : _ring_intervals,
908 84 : getParam<MooseEnum>("ring_id_assign_type") == "ring_wise",
909 42 : _quad_center_elements);
910 :
911 : // add internal side set info to metadata
912 1910 : if (_create_inward_interface_boundaries || _create_outward_interface_boundaries)
913 : {
914 2608 : 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 1304 : if (_create_outward_interface_boundaries)
919 : {
920 1297 : const unsigned int num_boundary_ids = boundary_ids.size();
921 7845 : for (const auto i : make_range(num_boundary_ids))
922 : {
923 6548 : const unsigned int id = i * 2 + 1 + _interface_boundary_id_shift;
924 6548 : auto it = boundary_ids.find(id);
925 6548 : if (it != boundary_ids.end())
926 : {
927 1632 : boundary_ids.erase(it);
928 1632 : interface_boundary_ids.insert(id);
929 : }
930 : }
931 : }
932 1304 : if (_create_inward_interface_boundaries)
933 : {
934 7 : const unsigned int num_boundary_ids = boundary_ids.size();
935 119 : for (const auto i : make_range(num_boundary_ids))
936 : {
937 112 : const unsigned int id = i * 2 + 2 + _interface_boundary_id_shift;
938 112 : auto it = boundary_ids.find(id);
939 112 : if (it != boundary_ids.end())
940 : {
941 28 : boundary_ids.erase(it);
942 28 : interface_boundary_ids.insert(id);
943 : }
944 : }
945 : }
946 2608 : setMeshProperty("interface_boundary_ids", interface_boundary_ids);
947 : }
948 :
949 1910 : bool flat_side_up = getMeshProperty<bool>("flat_side_up", name());
950 1910 : if (flat_side_up)
951 739 : MeshTools::Modification::rotate(*mesh0, 180.0 / (Real)_num_sides, 0.0, 0.0);
952 : mesh0->unset_is_prepared();
953 :
954 4772 : if (_has_rings && getParam<bool>("replace_inner_ring_with_delaunay_mesh"))
955 : {
956 28 : 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 3388 : for (auto & elem : mesh0->element_ptr_range())
964 1680 : if (elem->vertex_average().norm() < ring_radii_corr[0])
965 574 : deleteable_elems.insert(elem);
966 :
967 14 : const boundary_id_type boundary_id = MooseMeshUtils::getBoundaryIDs(*mesh0, {"_foo"}, true)[0];
968 : BoundaryInfo & boundary_info = mesh0->get_boundary_info();
969 574 : for (auto & elem : deleteable_elems)
970 : {
971 560 : unsigned int n_sides = elem->n_sides();
972 2520 : for (unsigned int n = 0; n != n_sides; ++n)
973 : {
974 1960 : Elem * neighbor = elem->neighbor_ptr(n);
975 1960 : if (!neighbor)
976 840 : continue;
977 :
978 1120 : const unsigned int return_side = neighbor->which_neighbor_am_i(elem);
979 :
980 1120 : if (neighbor->neighbor_ptr(return_side) == elem)
981 : {
982 : neighbor->set_neighbor(return_side, nullptr);
983 1120 : boundary_info.add_side(neighbor, return_side, boundary_id);
984 : }
985 : }
986 :
987 560 : mesh0->delete_elem(elem);
988 : }
989 :
990 : // build the 1d mesh from the new boundary
991 14 : auto poly_mesh = MooseMeshUtils::buildBoundaryMesh(*mesh0, boundary_id);
992 : Real min_side = std::numeric_limits<Real>::max();
993 588 : for (auto & elem : poly_mesh->element_ptr_range())
994 : {
995 280 : Real l = elem->volume();
996 280 : if (l < min_side)
997 : min_side = l;
998 14 : }
999 :
1000 : // triangulate with the 1d mesh
1001 14 : libMesh::Poly2TriTriangulator poly2tri(*poly_mesh);
1002 14 : poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
1003 14 : poly2tri.set_interpolate_boundary_points(0);
1004 : poly2tri.set_refine_boundary_allowed(false);
1005 : poly2tri.set_verify_hole_boundaries(false);
1006 42 : const auto desired_area = getParam<std::string>("inner_ring_desired_area");
1007 14 : if (desired_area != "")
1008 : {
1009 7 : poly2tri.desired_area() = 0;
1010 7 : libMesh::ParsedFunction<Real> area_func{desired_area};
1011 7 : poly2tri.set_desired_area_function(&area_func);
1012 7 : }
1013 : else
1014 7 : poly2tri.desired_area() = min_side * min_side;
1015 14 : poly2tri.minimum_angle() = 0;
1016 14 : poly2tri.smooth_after_generating() = true;
1017 : // poly2tri.elem_type() is TRI3 by default
1018 14 : if (_tri_elem_type == TRI_ELEM_TYPE::TRI6)
1019 0 : poly2tri.elem_type() = libMesh::ElemType::TRI6;
1020 14 : 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 14 : poly_mesh->add_point(libMesh::Point());
1024 14 : poly2tri.triangulate();
1025 : // keep the old subdomain id
1026 2044 : for (auto elem : poly_mesh->element_ptr_range())
1027 1022 : elem->subdomain_id() = block_ids_new[0];
1028 :
1029 28 : if (isParamValid("ring_id_name"))
1030 : {
1031 21 : 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 21 : auto id_name = getParam<std::string>("ring_id_name");
1039 14 : const auto extra_id_index = poly_mesh->add_elem_integer(id_name);
1040 910 : for (auto elem : poly_mesh->element_ptr_range())
1041 455 : elem->set_extra_integer(extra_id_index, 1);
1042 : }
1043 28 : if (isParamValid("sector_id_name"))
1044 : {
1045 : // we will assign all elements in the inner ring with zero sector id
1046 21 : auto id_name = getParam<std::string>("sector_id_name");
1047 14 : const auto extra_id_index = poly_mesh->add_elem_integer(id_name);
1048 910 : for (auto elem : poly_mesh->element_ptr_range())
1049 455 : elem->set_extra_integer(extra_id_index, 0);
1050 : }
1051 :
1052 : // stitch the triangulated mesh and the original mesh without the inner ring
1053 14 : mesh0->stitch_meshes(*poly_mesh, boundary_id, 0, TOLERANCE, true, false);
1054 14 : }
1055 :
1056 3820 : return dynamic_pointer_cast<MeshBase>(mesh0);
1057 1910 : }
|
