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 "PatternedHexMeshGenerator.h"
11 : #include "ReportingIDGeneratorUtils.h"
12 : #include "MooseUtils.h"
13 : #include "MooseMeshUtils.h"
14 :
15 : // C++ includes
16 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
17 : #include <fstream> // used to generate the optional control drum position file
18 :
19 : registerMooseObject("ReactorApp", PatternedHexMeshGenerator);
20 :
21 : InputParameters
22 1994 : PatternedHexMeshGenerator::validParams()
23 : {
24 1994 : InputParameters params = PolygonMeshGeneratorBase::validParams();
25 3988 : params.addRequiredParam<std::vector<MeshGeneratorName>>("inputs", "The input MeshGenerators.");
26 3988 : params.addRequiredRangeCheckedParam<std::vector<std::vector<unsigned int>>>(
27 : "pattern",
28 : "pattern>=0",
29 : "A double-indexed hexagonal-shaped array starting with the upper-left corner.");
30 3988 : MooseEnum pattern_boundary("none hexagon", "hexagon");
31 3988 : params.addParam<MooseEnum>(
32 : "pattern_boundary", pattern_boundary, "The boundary shape of the patterned mesh.");
33 3988 : params.addParam<bool>(
34 : "generate_core_metadata",
35 3988 : false,
36 : "A Boolean parameter that controls whether the core related metadata "
37 : "is generated for other MOOSE objects such as 'MultiControlDrumFunction' or not.");
38 5982 : params.addRangeCheckedParam<unsigned int>("background_intervals",
39 3988 : 3,
40 : "background_intervals>0",
41 : "Radial intervals in the assembly peripheral region.");
42 3988 : params.addRangeCheckedParam<Real>("hexagon_size",
43 : "hexagon_size>0.0",
44 : "Size of the outmost hexagon boundary to be generated; this is "
45 : "required only when pattern type is 'hexagon'.");
46 3988 : MooseEnum hexagon_size_style("apothem radius", "apothem");
47 3988 : params.addParam<MooseEnum>(
48 : "hexagon_size_style",
49 : hexagon_size_style,
50 : "Style in which the hexagon size is given (default: apothem i.e. half-pitch).");
51 3988 : params.addRangeCheckedParam<std::vector<Real>>(
52 : "duct_sizes", "duct_sizes>0.0", "Distance(s) from center to duct(s) inner boundaries.");
53 3988 : MooseEnum duct_sizes_style("apothem radius", "apothem");
54 3988 : params.addParam<MooseEnum>("duct_sizes_style",
55 : duct_sizes_style,
56 : "Style in which hexagon center to duct distance(s) is given (apothem "
57 : "= center-to-face, radius = center-to-vertex).");
58 3988 : params.addRangeCheckedParam<std::vector<unsigned int>>(
59 : "duct_intervals", "duct_intervals>0", "Number of meshing intervals in each enclosing duct.");
60 3988 : params.addParam<bool>("uniform_mesh_on_sides",
61 3988 : false,
62 : "Whether the side elements are reorganized to have a uniform size.");
63 3988 : params.addParam<bool>("generate_control_drum_positions_file",
64 3988 : false,
65 : "Whether a positions file is generated in the core mesh mode.");
66 3988 : params.addParam<bool>("assign_control_drum_id",
67 3988 : false,
68 : "Whether control drum id is assigned to the mesh as an extra integer.");
69 1994 : std::string position_file_default = "positions_meta.data";
70 3988 : params.addParam<std::string>(
71 : "position_file", position_file_default, "Data file name to store control drum positions.");
72 : // A hexagon pattern_boundary mesh can be used in "inputs" of `PatternedHexMeshGenerator` after a
73 : // 90-degree rotation.
74 3988 : params.addParam<Real>(
75 : "rotate_angle",
76 3988 : 90.0,
77 : "Rotate the entire patterned mesh by a certain degrees that is defined here.");
78 3988 : params.addParam<subdomain_id_type>(
79 : "background_block_id",
80 : "Optional customized block id for the background block in 'assembly' mode; must be provided "
81 : "along with 'duct_block_ids' if 'duct_sizes' is provided.");
82 3988 : params.addParam<SubdomainName>(
83 : "background_block_name",
84 : "Optional customized block name for the background block in 'assembly' mode; must be "
85 : "provided along with 'duct_block_names' if 'duct_sizes' is provided.");
86 3988 : params.addParam<std::vector<subdomain_id_type>>(
87 : "duct_block_ids",
88 : "Optional customized block ids for each duct geometry block in 'assembly' mode; must be "
89 : "provided along with 'background_block_id'.");
90 1994 : params.addParam<std::vector<SubdomainName>>(
91 : "duct_block_names",
92 1994 : std::vector<SubdomainName>(),
93 : "Optional customized block names for each duct geometry block in 'assembly' mode; must be "
94 : "provided along with 'background_block_name'.");
95 3988 : params.addRangeCheckedParam<boundary_id_type>("external_boundary_id",
96 : "external_boundary_id>0",
97 : "Optional customized external boundary id.");
98 3988 : params.addParam<bool>("create_inward_interface_boundaries",
99 3988 : false,
100 : "Whether the inward interface boundary sidesets are created.");
101 3988 : params.addParam<bool>("create_outward_interface_boundaries",
102 3988 : true,
103 : "Whether the outward interface boundary sidesets are created.");
104 3988 : params.addParam<BoundaryName>(
105 1994 : "external_boundary_name", BoundaryName(), "Optional customized external boundary name.");
106 3988 : params.addParam<bool>("deform_non_circular_region",
107 3988 : true,
108 : "Whether the non-circular region (outside the rings) can be deformed.");
109 3988 : params.addParam<std::vector<std::string>>("id_name", "List of extra integer ID set names");
110 3988 : params.addParam<std::vector<MeshGeneratorName>>(
111 : "exclude_id", "Name of input meshes to be excluded in ID generation");
112 3988 : std::vector<MooseEnum> option = {MooseEnum("cell pattern manual", "cell")};
113 3988 : params.addParam<std::vector<MooseEnum>>(
114 : "assign_type", option, "List of integer ID assignment types");
115 3988 : params.addParam<std::vector<std::vector<std::vector<dof_id_type>>>>(
116 : "id_pattern",
117 : "User-defined element IDs. A double-indexed array starting with the upper-left corner. When "
118 : "providing multiple patterns, each pattern should be separated using '|'");
119 3988 : params.addParam<std::vector<std::vector<boundary_id_type>>>(
120 : "interface_boundary_id_shift_pattern",
121 : "User-defined shift values for each pattern cell. A double-indexed array starting with the "
122 : "upper-left corner.");
123 3988 : MooseEnum quad_elem_type("QUAD4 QUAD8 QUAD9", "QUAD4");
124 3988 : params.addParam<MooseEnum>(
125 : "boundary_region_element_type",
126 : quad_elem_type,
127 : "Type of the quadrilateral elements to be generated in the boundary region.");
128 3988 : params.addParam<bool>(
129 : "allow_unused_inputs",
130 3988 : false,
131 : "Whether additional input assemblies can be part of inputs without being used in lattice");
132 3988 : params.addParamNamesToGroup(
133 : "background_block_id background_block_name duct_block_ids boundary_region_element_type "
134 : "duct_block_names external_boundary_id external_boundary_name "
135 : "create_inward_interface_boundaries create_outward_interface_boundaries",
136 : "Customized Subdomain/Boundary");
137 3988 : params.addParamNamesToGroup(
138 : "generate_control_drum_positions_file assign_control_drum_id position_file", "Control Drum");
139 3988 : params.addParamNamesToGroup(
140 : "background_intervals duct_intervals uniform_mesh_on_sides deform_non_circular_region",
141 : "Mesh Density");
142 3988 : params.addParamNamesToGroup("id_name exclude_id assign_type id_pattern", "Reporting ID");
143 1994 : params.addClassDescription(
144 : "This PatternedHexMeshGenerator source code assembles hexagonal meshes into a hexagonal grid "
145 : "and optionally forces the outer boundary to be hexagonal and/or adds a duct.");
146 :
147 1994 : return params;
148 7976 : }
149 :
150 1007 : PatternedHexMeshGenerator::PatternedHexMeshGenerator(const InputParameters & parameters)
151 : : PolygonMeshGeneratorBase(parameters),
152 1007 : _mesh_ptrs(getMeshes("inputs")),
153 2014 : _input_names(getParam<std::vector<MeshGeneratorName>>("inputs")),
154 2014 : _pattern(getParam<std::vector<std::vector<unsigned int>>>("pattern")),
155 2014 : _pattern_boundary(getParam<MooseEnum>("pattern_boundary")),
156 2014 : _generate_core_metadata(getParam<bool>("generate_core_metadata")),
157 2014 : _background_intervals(getParam<unsigned int>("background_intervals")),
158 2014 : _has_assembly_duct(isParamValid("duct_sizes")),
159 3322 : _duct_sizes(isParamValid("duct_sizes") ? getParam<std::vector<Real>>("duct_sizes")
160 : : std::vector<Real>()),
161 2014 : _duct_sizes_style(getParam<MooseEnum>("duct_sizes_style").template getEnum<PolygonSizeStyle>()),
162 2616 : _duct_intervals(isParamValid("duct_intervals")
163 1007 : ? getParam<std::vector<unsigned int>>("duct_intervals")
164 : : std::vector<unsigned int>()),
165 2014 : _uniform_mesh_on_sides(getParam<bool>("uniform_mesh_on_sides")),
166 2014 : _generate_control_drum_positions_file(getParam<bool>("generate_control_drum_positions_file")),
167 2014 : _assign_control_drum_id(getParam<bool>("assign_control_drum_id")),
168 2014 : _rotate_angle(getParam<Real>("rotate_angle")),
169 2500 : _duct_block_ids(isParamValid("duct_block_ids")
170 1007 : ? getParam<std::vector<subdomain_id_type>>("duct_block_ids")
171 : : std::vector<subdomain_id_type>()),
172 2014 : _duct_block_names(getParam<std::vector<SubdomainName>>("duct_block_names")),
173 2427 : _external_boundary_id(isParamValid("external_boundary_id")
174 1833 : ? getParam<boundary_id_type>("external_boundary_id")
175 : : 0),
176 1007 : _external_boundary_name(getParam<BoundaryName>("external_boundary_name")),
177 2014 : _create_inward_interface_boundaries(getParam<bool>("create_inward_interface_boundaries")),
178 2014 : _create_outward_interface_boundaries(getParam<bool>("create_outward_interface_boundaries")),
179 1007 : _hexagon_size_style(
180 1007 : getParam<MooseEnum>("hexagon_size_style").template getEnum<PolygonSizeStyle>()),
181 2014 : _deform_non_circular_region(getParam<bool>("deform_non_circular_region")),
182 2014 : _use_reporting_id(isParamValid("id_name")),
183 2014 : _use_exclude_id(isParamValid("exclude_id")),
184 2014 : _use_interface_boundary_id_shift(isParamValid("interface_boundary_id_shift_pattern")),
185 1007 : _boundary_quad_elem_type(
186 1007 : getParam<MooseEnum>("boundary_region_element_type").template getEnum<QUAD_ELEM_TYPE>()),
187 3021 : _allow_unused_inputs(getParam<bool>("allow_unused_inputs"))
188 : {
189 1007 : declareMeshProperty("pattern_pitch_meta", 0.0);
190 1007 : declareMeshProperty("input_pitch_meta", 0.0);
191 2014 : declareMeshProperty<bool>("is_control_drum_meta", false);
192 2014 : declareMeshProperty<std::vector<Point>>("control_drum_positions", std::vector<Point>());
193 2014 : declareMeshProperty<std::vector<Real>>("control_drum_angles", std::vector<Real>());
194 1007 : declareMeshProperty<std::vector<std::vector<Real>>>("control_drums_azimuthal_meta",
195 1007 : std::vector<std::vector<Real>>());
196 3021 : declareMeshProperty<std::string>("position_file_name", getParam<std::string>("position_file"));
197 1007 : declareMeshProperty<bool>("hexagon_peripheral_trimmability", !_generate_core_metadata);
198 1007 : declareMeshProperty<bool>("hexagon_center_trimmability", true);
199 1007 : declareMeshProperty<bool>("peripheral_modifier_compatible", _pattern_boundary == "hexagon");
200 :
201 1007 : const unsigned int n_pattern_layers = _pattern.size();
202 1007 : declareMeshProperty("pattern_size", n_pattern_layers);
203 1007 : if (n_pattern_layers % 2 == 0)
204 2 : paramError(
205 : "pattern",
206 : "The length (layer number) of this parameter must be odd to ensure a hexagonal shape.");
207 1005 : if (n_pattern_layers == 1)
208 2 : paramError("pattern", "The length (layer number) of this parameter must be larger than unity.");
209 : std::vector<unsigned int> pattern_max_array;
210 : std::vector<unsigned int> pattern_1d;
211 3188 : for (unsigned int i = 0; i <= n_pattern_layers / 2; i++)
212 : {
213 4374 : pattern_max_array.push_back(*std::max_element(_pattern[i].begin(), _pattern[i].end()));
214 2187 : pattern_max_array.push_back(*std::max_element(_pattern[n_pattern_layers - 1 - i].begin(),
215 2187 : _pattern[n_pattern_layers - 1 - i].end()));
216 2187 : pattern_1d.insert(pattern_1d.end(), _pattern[i].begin(), _pattern[i].end());
217 2187 : pattern_1d.insert(pattern_1d.end(),
218 : _pattern[n_pattern_layers - 1 - i].begin(),
219 2187 : _pattern[n_pattern_layers - 1 - i].end());
220 2187 : if (_pattern[i].size() != n_pattern_layers / 2 + i + 1 ||
221 : _pattern[n_pattern_layers - 1 - i].size() != n_pattern_layers / 2 + i + 1)
222 2 : paramError(
223 : "pattern",
224 : "The two-dimentional array parameter pattern must have a correct hexagonal shape.");
225 : }
226 1001 : if (*std::max_element(pattern_max_array.begin(), pattern_max_array.end()) >= _input_names.size())
227 2 : paramError("pattern",
228 : "Elements of this parameter must be smaller than the length of input_name.");
229 1998 : if (std::set<unsigned int>(pattern_1d.begin(), pattern_1d.end()).size() < _input_names.size() &&
230 38 : !_allow_unused_inputs)
231 2 : paramError("pattern",
232 : "All the meshes provided in inputs must be used in the lattice pattern. To bypass "
233 : "this requirement, set 'allow_unused_inputs = true'");
234 :
235 1994 : if (isParamValid("background_block_id"))
236 : {
237 1443 : _peripheral_block_ids.push_back(getParam<subdomain_id_type>("background_block_id"));
238 481 : _peripheral_block_ids.insert(
239 : _peripheral_block_ids.end(), _duct_block_ids.begin(), _duct_block_ids.end());
240 : }
241 516 : else if (!_duct_block_ids.empty())
242 2 : paramError("background_block_id",
243 : "This parameter and duct_block_ids must be "
244 : "provided simultaneously.");
245 1990 : if (isParamValid("background_block_name"))
246 : {
247 1410 : _peripheral_block_names.push_back(getParam<SubdomainName>("background_block_name"));
248 470 : _peripheral_block_names.insert(
249 : _peripheral_block_names.end(), _duct_block_names.begin(), _duct_block_names.end());
250 : }
251 525 : else if (!_duct_block_names.empty())
252 2 : paramError("background_block_name",
253 : "This parameter and duct_block_names must be provided simultaneously.");
254 :
255 993 : if (_pattern_boundary == "hexagon")
256 : {
257 738 : for (unsigned int i = 1; i < _duct_sizes.size(); i++)
258 69 : if (_duct_sizes[i] <= _duct_sizes[i - 1])
259 2 : paramError("duct_sizes", "This parameter must be strictly ascending.");
260 669 : if (!_peripheral_block_ids.empty() && _peripheral_block_ids.size() != _duct_sizes.size() + 1)
261 2 : paramError("duct_block_ids",
262 : "This parameter, if provided, must have a length equal to length of duct_sizes.");
263 667 : if (!_peripheral_block_names.empty() &&
264 468 : _peripheral_block_names.size() != _duct_sizes.size() + 1)
265 2 : paramError("duct_block_names",
266 : "This parameter, if provided, must have a length equal to length of duct_sizes.");
267 : }
268 : else
269 : {
270 322 : if (!_peripheral_block_ids.empty() || !_peripheral_block_names.empty())
271 2 : paramError("background_block_id",
272 : "This parameter and background_block_name must not be set when the "
273 : "pattern_boundary is none.");
274 : }
275 :
276 985 : if (_use_interface_boundary_id_shift)
277 : {
278 : // check "interface_boundary_id_shift_pattern" parameter
279 : _interface_boundary_id_shift_pattern =
280 27 : getParam<std::vector<std::vector<boundary_id_type>>>("interface_boundary_id_shift_pattern");
281 9 : if (_interface_boundary_id_shift_pattern.size() != _pattern.size())
282 0 : paramError("interface_boundary_id_shift_pattern",
283 : "This parameter, if provided, should have the same two-dimensional array shape as "
284 : "the 'pattern' parameter. First dimension does not match");
285 36 : for (const auto i : make_range(_pattern.size()))
286 27 : if (_interface_boundary_id_shift_pattern[i].size() != _pattern[i].size())
287 0 : paramError("interface_boundary_id_shift_pattern",
288 : "This parameter, if provided, should have the same two-dimensional array shape "
289 : "as the 'pattern' parameter.");
290 : }
291 : // declare metadata for internal interface boundaries
292 1970 : declareMeshProperty<bool>("interface_boundaries", false);
293 1970 : declareMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", {});
294 :
295 985 : if (_use_reporting_id)
296 : {
297 : // get reporting id name input
298 1602 : _reporting_id_names = getParam<std::vector<std::string>>("id_name");
299 534 : const unsigned int num_reporting_ids = _reporting_id_names.size();
300 : // get reporting id assign type input
301 1602 : const auto input_assign_types = getParam<std::vector<MooseEnum>>("assign_type");
302 534 : if (input_assign_types.size() != num_reporting_ids)
303 0 : paramError("assign_type", "This parameter must have a length equal to length of id_name.");
304 : // list of reporting id names using manual id patterns;
305 : std::vector<std::string> manual_ids;
306 1095 : for (const auto i : make_range(num_reporting_ids))
307 : {
308 561 : _assign_types.push_back(
309 561 : input_assign_types[i].getEnum<ReportingIDGeneratorUtils::AssignType>());
310 561 : if (_assign_types[i] == ReportingIDGeneratorUtils::AssignType::manual)
311 27 : manual_ids.push_back(_reporting_id_names[i]);
312 : }
313 : // processing "id_pattern" input parameter
314 588 : if (manual_ids.size() > 0 && !isParamValid("id_pattern"))
315 0 : paramError("id_pattern", "required when 'manual' is defined in \"assign_type\"");
316 1068 : if (isParamValid("id_pattern"))
317 : {
318 : const auto input_id_patterns =
319 54 : getParam<std::vector<std::vector<std::vector<dof_id_type>>>>("id_pattern");
320 18 : if (input_id_patterns.size() != manual_ids.size())
321 0 : paramError("id_pattern",
322 : "The number of patterns must be equal to the number of 'manual' types defined "
323 : "in \"assign_type\".");
324 45 : for (unsigned int i = 0; i < manual_ids.size(); ++i)
325 27 : _id_patterns[manual_ids[i]] = input_id_patterns[i];
326 18 : }
327 : // processing exlude id
328 534 : _exclude_ids.resize(_input_names.size());
329 534 : if (_use_exclude_id)
330 : {
331 : std::vector<MeshGeneratorName> exclude_id_name =
332 240 : getParam<std::vector<MeshGeneratorName>>("exclude_id");
333 530 : for (unsigned int i = 0; i < _input_names.size(); ++i)
334 : {
335 : _exclude_ids[i] = false;
336 700 : for (auto input_name : exclude_id_name)
337 410 : if (_input_names[i] == input_name)
338 : {
339 : _exclude_ids[i] = true;
340 : break;
341 : }
342 : }
343 120 : }
344 : else
345 1033 : for (unsigned int i = 0; i < _input_names.size(); ++i)
346 : _exclude_ids[i] = false;
347 534 : }
348 985 : }
349 :
350 : std::unique_ptr<MeshBase>
351 957 : PatternedHexMeshGenerator::generate()
352 : {
353 957 : std::vector<std::unique_ptr<ReplicatedMesh>> meshes(_input_names.size());
354 2652 : for (const auto i : index_range(_input_names))
355 : {
356 3394 : meshes[i] = dynamic_pointer_cast<ReplicatedMesh>(std::move(*_mesh_ptrs[i]));
357 1697 : if (!meshes[i])
358 0 : paramError("inputs",
359 : "Mesh '",
360 0 : _input_names[i],
361 : "' is not a replicated mesh. Only replicated meshes are supported");
362 :
363 : // throw an error message if the input mesh has a flat side up
364 3394 : if (hasMeshProperty<bool>("flat_side_up", _input_names[i]))
365 714 : if (getMeshProperty<bool>("flat_side_up", _input_names[i]))
366 2 : paramError("inputs",
367 : "Mesh '",
368 2 : _input_names[i],
369 : "' has a flat side facing up, which is not supported.");
370 : }
371 :
372 : std::vector<Real> pitch_array;
373 : std::vector<unsigned int> num_sectors_per_side_array;
374 : std::vector<unsigned int> num_sectors_per_side_array_tmp;
375 : std::vector<std::vector<Real>> control_drum_azimuthal_array;
376 : std::vector<unsigned int> background_intervals_array;
377 : std::vector<dof_id_type> node_id_background_array;
378 : std::vector<Real> max_radius_array;
379 : std::vector<bool> is_control_drum_array;
380 : Real max_radius_global(0.0);
381 : std::vector<Real> pattern_pitch_array;
382 :
383 955 : if (_pattern_boundary == "none" && _generate_core_metadata)
384 : {
385 1014 : for (MooseIndex(_input_names) i = 0; i < _input_names.size(); ++i)
386 : {
387 : // throw an error message if the input mesh does not contain the required meta data
388 1524 : if (!hasMeshProperty<Real>("pattern_pitch_meta", _input_names[i]))
389 0 : mooseError(
390 : "In PatternedHexMeshGenerator ",
391 0 : _name,
392 : ": the unit hexagonal input mesh does not contain appropriate meta data "
393 : "required for generating a core mesh. Involved input mesh: ",
394 0 : _input_names[i],
395 : "; metadata issue: 'pattern_pitch_meta' is missing. Note that "
396 : "'generate_core_metadata' is set to true, which"
397 : "means that the mesh generator is producing a core mesh by stitching the input "
398 : "assembly meshes together. Therefore,"
399 : "the input meshes must contain the metadata of assembly meshes, which can "
400 : "usually be either automatically assigned "
401 : "by using another PatternedHexMeshGenerator with 'generate_core_metadata' set as "
402 : "false and 'pattern_boundary' set as hexagon, or manually assigned by "
403 : "AddMetaDataGenerator.");
404 1524 : pattern_pitch_array.push_back(getMeshProperty<Real>("pattern_pitch_meta", _input_names[i]));
405 : // throw an error message if the input mesh contains non-sense meta data
406 762 : if (pattern_pitch_array.back() == 0.0)
407 2 : mooseError(
408 : "In PatternedHexMeshGenerator ",
409 2 : _name,
410 : ": the unit hexagonal input mesh does not contain appropriate meta data "
411 : "required for generating a core mesh. Involved input mesh: ",
412 2 : _input_names[i],
413 : "; metadata issue: 'pattern_pitch_meta' is zero. Note that "
414 : "'generate_core_metadata' is set to true, which"
415 : "means that the mesh generator is producing a core mesh by stitching the input "
416 : "assembly meshes together. Therefore,"
417 : "the input meshes must contain the metadata of assembly meshes, which can "
418 : "usually be either automatically assigned "
419 : "by using another PatternedHexMeshGenerator with 'generate_core_metadata' set as "
420 : "false and 'pattern_boundary' set as hexagon, or manually assigned by "
421 : "AddMetaDataGenerator.");
422 760 : is_control_drum_array.push_back(
423 760 : getMeshProperty<bool>("is_control_drum_meta", _input_names[i]));
424 760 : control_drum_azimuthal_array.push_back(
425 760 : getMeshProperty<bool>("is_control_drum_meta", _input_names[i])
426 1673 : ? getMeshProperty<std::vector<Real>>("azimuthal_angle_meta", _input_names[i])
427 : : std::vector<Real>());
428 : }
429 252 : if (!MooseUtils::absoluteFuzzyEqual(
430 : *std::max_element(pattern_pitch_array.begin(), pattern_pitch_array.end()),
431 : *std::min_element(pattern_pitch_array.begin(), pattern_pitch_array.end())))
432 4 : mooseError(
433 : "In PatternedHexMeshGenerator ",
434 2 : _name,
435 : ": pattern_pitch metadata values of all input mesh generators must be identical "
436 : "when pattern_boundary is 'none' and generate_core_metadata is true. Please check the "
437 : "parameters of the mesh generators that produce the input meshes."
438 : "Note that some of these mesh generators, such as "
439 : "HexagonConcentricCircleAdaptiveBoundaryMeshGenerator and FlexiblePatternGenerator,"
440 : "may have different definitions of hexagon size in their input parameters. Please refer "
441 : "to the documentation of these mesh generators.",
442 2 : pitchMetaDataErrorGenerator(_input_names, pattern_pitch_array, "pattern_pitch_meta"));
443 : else
444 : {
445 250 : _pattern_pitch = pattern_pitch_array.front();
446 500 : setMeshProperty("input_pitch_meta", _pattern_pitch);
447 : }
448 : }
449 : else
450 : {
451 701 : if (_pattern_boundary == "hexagon")
452 : {
453 1286 : if (!isParamValid("hexagon_size"))
454 2 : paramError("hexagon_size",
455 : "This parameter must be provided when pattern_boundary is hexagon.");
456 : else
457 1282 : _pattern_pitch = 2.0 * (_hexagon_size_style == PolygonSizeStyle::apothem
458 1923 : ? getParam<Real>("hexagon_size")
459 641 : : getParam<Real>("hexagon_size") * std::cos(M_PI / 6.0));
460 : }
461 1626 : for (MooseIndex(_input_names) i = 0; i < _input_names.size(); ++i)
462 : {
463 : // throw an error message if the input mesh does not contain the required meta data
464 1862 : if (!hasMeshProperty<Real>("pitch_meta", _input_names[i]))
465 2 : mooseError("In PatternedHexMeshGenerator ",
466 2 : _name,
467 : ": the unit hexagonal input mesh does not contain appropriate meta data "
468 : "required for generating an assembly. Involved input mesh: ",
469 2 : _input_names[i],
470 : "; metadata issue: 'pitch_meta' is missing");
471 929 : pitch_array.push_back(getMeshProperty<Real>("pitch_meta", _input_names[i]));
472 :
473 : num_sectors_per_side_array_tmp =
474 1858 : getMeshProperty<std::vector<unsigned int>>("num_sectors_per_side_meta", _input_names[i]);
475 929 : if (*std::max_element(num_sectors_per_side_array_tmp.begin(),
476 : num_sectors_per_side_array_tmp.end()) !=
477 929 : *std::min_element(num_sectors_per_side_array_tmp.begin(),
478 : num_sectors_per_side_array_tmp.end()))
479 0 : mooseError("In PatternedHexMeshGenerator ",
480 0 : _name,
481 : ": num_sectors_per_side metadata values of all six sides of each input mesh "
482 : "generator must be identical.");
483 929 : if (num_sectors_per_side_array_tmp.front() == 1 && _pattern_boundary == "hexagon")
484 2 : paramError(
485 : "inputs",
486 : "for each input mesh, the number of sectors on each side must be greater than unity.");
487 927 : num_sectors_per_side_array.push_back(*num_sectors_per_side_array_tmp.begin());
488 927 : background_intervals_array.push_back(
489 927 : getMeshProperty<unsigned int>("background_intervals_meta", _input_names[i]));
490 927 : node_id_background_array.push_back(
491 927 : getMeshProperty<dof_id_type>("node_id_background_meta", _input_names[i]));
492 1854 : max_radius_array.push_back(getMeshProperty<Real>("max_radius_meta", _input_names[i]));
493 : }
494 695 : max_radius_global = *max_element(max_radius_array.begin(), max_radius_array.end());
495 695 : if (!MooseUtils::absoluteFuzzyEqual(*std::max_element(pitch_array.begin(), pitch_array.end()),
496 : *std::min_element(pitch_array.begin(), pitch_array.end())))
497 0 : mooseError("In PatternedHexMeshGenerator ",
498 0 : _name,
499 : ": pitch metadata values of all input mesh generators must be identical. Please "
500 : "check the "
501 : "parameters of the mesh generators that produce the input meshes.",
502 0 : pitchMetaDataErrorGenerator(_input_names, pitch_array, "pitch_meta"));
503 1390 : setMeshProperty("input_pitch_meta", pitch_array.front());
504 695 : if (*std::max_element(num_sectors_per_side_array.begin(), num_sectors_per_side_array.end()) !=
505 695 : *std::min_element(num_sectors_per_side_array.begin(), num_sectors_per_side_array.end()))
506 0 : mooseError(
507 : "In PatternedHexMeshGenerator ",
508 0 : _name,
509 : ": num_sectors_per_side metadata values of all input mesh generators must be identical.");
510 : }
511 :
512 : std::vector<Real> extra_dist;
513 : Real extra_dist_shift(0.0);
514 : Real y_min(0.0);
515 : Real y_max_0(0.0);
516 : Real y_max_n(0.0);
517 945 : const Real extra_dist_tol = _pattern_boundary == "hexagon" ? pitch_array.front() / 10.0 : 0.0;
518 945 : const Real extra_dist_shift_0 = _pattern_boundary == "hexagon" ? pitch_array.front() / 5.0 : 0.0;
519 : std::vector<unsigned int> peripheral_duct_intervals;
520 945 : if (_pattern_boundary == "hexagon")
521 : {
522 639 : if (_has_assembly_duct)
523 : {
524 645 : for (unsigned int i = 0; i < _duct_sizes.size(); i++)
525 : {
526 354 : if (_duct_sizes_style == PolygonSizeStyle::radius)
527 0 : _duct_sizes[i] /= std::cos(M_PI / 6.0);
528 354 : extra_dist.push_back(0.5 *
529 354 : (_duct_sizes[i] * 2.0 - pitch_array.front() / std::sqrt(3.0) *
530 354 : (Real)((_pattern.size() / 2) * 3 + 2)));
531 354 : peripheral_duct_intervals.push_back(_duct_intervals[i]);
532 : }
533 291 : if (_duct_sizes.back() >= _pattern_pitch / 2.0)
534 2 : paramError("duct_sizes",
535 : "The duct sizes should not exceed the size of the hexagonal boundary.");
536 : }
537 : // calculate the distance between the larger hexagon boundary and the boundary of stitched unit
538 : // hexagons this is used to decide whether deformation is needed when cut-off happens or when
539 : // the distance is small.
540 637 : extra_dist.push_back(0.5 * (_pattern_pitch - pitch_array.front() / std::sqrt(3.0) *
541 637 : (Real)((_pattern.size() / 2) * 3 + 2)));
542 637 : peripheral_duct_intervals.insert(peripheral_duct_intervals.begin(), _background_intervals);
543 :
544 : // In same cases, when the external hexagon size is small enough, the external hexagon boundary
545 : // may either be very close to the input hexagon meshes that are near the boundary or even cut
546 : // off the these hexagons. As long as the ring regions are not cut off, the input hexagons can
547 : // be deformed to accommodate the external hexagon shape. This block sets up the range of mesh
548 : // region that needs to be deformed.
549 637 : if (extra_dist.front() <= extra_dist_tol)
550 : {
551 159 : extra_dist_shift = extra_dist_shift_0 - extra_dist.front();
552 480 : for (Real & d : extra_dist)
553 321 : d += extra_dist_shift;
554 159 : y_min = _deform_non_circular_region
555 159 : ? max_radius_global // Currently use this, ideally this should be the max of the
556 : // outer layer radii
557 20 : : (pitch_array.front() / std::sqrt(3.0));
558 159 : y_max_0 = pitch_array.front() / std::sqrt(3.0) + extra_dist.front();
559 159 : y_max_n = y_max_0 - extra_dist_shift;
560 159 : if (y_max_n <= y_min)
561 2 : mooseError("In PatternedHexMeshGenerator ",
562 2 : _name,
563 : ": the assembly is cut off so much that the internal structure that should not "
564 : "be altered is compromised.");
565 : }
566 : }
567 :
568 1882 : setMeshProperty("pattern_pitch_meta", _pattern_pitch);
569 :
570 : // create a list of interface boundary ids for each input mesh
571 : // NOTE: list of interface boundary ids is stored in mesh metadata
572 : std::vector<std::set<boundary_id_type>> input_interface_boundary_ids;
573 941 : input_interface_boundary_ids.resize(_input_names.size());
574 941 : if (_use_interface_boundary_id_shift)
575 : {
576 18 : for (const auto i : make_range(_input_names.size()))
577 : {
578 18 : if (!hasMeshProperty<bool>("interface_boundaries", _input_names[i]))
579 0 : mooseError("Metadata 'interface_boundaries' could not be found on the input mesh: ",
580 0 : _input_names[i]);
581 9 : if (!getMeshProperty<bool>("interface_boundaries", _input_names[i]))
582 0 : mooseError("Interface boundary ids were not constructed in the input mesh",
583 0 : _input_names[i]);
584 18 : if (!hasMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", _input_names[i]))
585 0 : mooseError("Metadata 'interface_boundary_ids' could not be found on the input mesh: ",
586 0 : _input_names[i]);
587 : }
588 : }
589 2620 : for (const auto i : make_range(_input_names.size()))
590 1679 : if (_use_interface_boundary_id_shift ||
591 3349 : hasMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", _input_names[i]))
592 : input_interface_boundary_ids[i] =
593 1640 : getMeshProperty<std::set<boundary_id_type>>("interface_boundary_ids", _input_names[i]);
594 :
595 : int x_mov = 0;
596 :
597 1247 : const Real input_pitch((_pattern_boundary == "hexagon" || !_generate_core_metadata)
598 1247 : ? pitch_array.front()
599 : : _pattern_pitch);
600 : std::vector<Real> control_drum_positions_x;
601 : std::vector<Real> control_drum_positions_y;
602 : std::vector<std::vector<Real>> control_drum_azimuthals;
603 :
604 941 : std::unique_ptr<ReplicatedMesh> out_mesh;
605 :
606 4122 : for (unsigned i = 0; i < _pattern.size(); i++)
607 : {
608 3183 : Real deltax = -x_mov * input_pitch / 2;
609 3183 : if (i == _pattern.size() - 1)
610 939 : x_mov--;
611 2244 : else if (_pattern[i].size() < _pattern[i + 1].size())
612 1122 : x_mov++;
613 : else
614 1122 : x_mov--;
615 3183 : Real deltay = -(Real)(i)*input_pitch * std::sin(M_PI / 3);
616 :
617 12204 : for (unsigned int j = 0; j < _pattern[i].size(); j++)
618 : {
619 9023 : const auto pattern = _pattern[i][j];
620 9023 : ReplicatedMesh & pattern_mesh = *meshes[pattern];
621 :
622 : // No pattern boundary, so no need to wrap with a peripheral mesh
623 : // Use the inputs as they stand and translate accordingly later
624 9023 : if (_pattern_boundary == "none")
625 : {
626 2460 : if (_generate_core_metadata && is_control_drum_array[pattern] == 1)
627 : {
628 153 : control_drum_positions_x.push_back(deltax + j * input_pitch);
629 153 : control_drum_positions_y.push_back(deltay);
630 153 : control_drum_azimuthals.push_back(control_drum_azimuthal_array[pattern]);
631 : }
632 :
633 2460 : if (j == 0 && i == 0)
634 : {
635 612 : out_mesh = dynamic_pointer_cast<ReplicatedMesh>(pattern_mesh.clone());
636 306 : if (_assign_control_drum_id && _generate_core_metadata)
637 306 : out_mesh->add_elem_integer(
638 : "control_drum_id",
639 : true,
640 : is_control_drum_array[pattern] ? control_drum_azimuthals.size() : 0);
641 : // Reassign interface boundary ids
642 306 : if (_use_interface_boundary_id_shift)
643 9 : reassignBoundaryIDs(*out_mesh,
644 : _interface_boundary_id_shift_pattern[i][j],
645 : input_interface_boundary_ids[pattern]);
646 5040 : continue;
647 : }
648 : }
649 : // Has a pattern boundary so we need to wrap with a peripheral mesh
650 : else // has a pattern boundary
651 : {
652 : Real rotation_angle = std::numeric_limits<Real>::max();
653 : Real orientation = std::numeric_limits<Real>::max();
654 : unsigned int mesh_type = std::numeric_limits<unsigned int>::max();
655 : bool on_periphery = true;
656 :
657 6563 : if (j == 0 && i == 0)
658 : {
659 : rotation_angle = 60.;
660 : mesh_type = CORNER_MESH;
661 : orientation = 0.;
662 : }
663 5928 : else if (j == 0 && i == _pattern.size() - 1)
664 : {
665 : rotation_angle = 180.;
666 : mesh_type = CORNER_MESH;
667 : orientation = -120.;
668 : }
669 5293 : else if (j == 0 && i == (_pattern.size() - 1) / 2)
670 : {
671 : rotation_angle = 120.;
672 : mesh_type = CORNER_MESH;
673 : orientation = -60.;
674 : }
675 4658 : else if (j == _pattern[i].size() - 1 && i == 0)
676 : {
677 : rotation_angle = 0.;
678 : mesh_type = CORNER_MESH;
679 : orientation = 60.;
680 : }
681 4023 : else if (j == _pattern[i].size() - 1 && i == _pattern.size() - 1)
682 : {
683 : rotation_angle = -120.;
684 : mesh_type = CORNER_MESH;
685 : orientation = 180.;
686 : }
687 3388 : else if (j == _pattern[i].size() - 1 && i == (_pattern.size() - 1) / 2)
688 : {
689 : rotation_angle = -60.;
690 : mesh_type = CORNER_MESH;
691 : orientation = 120.;
692 : }
693 2753 : else if (i == 0)
694 : {
695 : rotation_angle = 0.;
696 : mesh_type = BOUNDARY_MESH;
697 : orientation = 0.;
698 : }
699 2599 : else if (i == _pattern.size() - 1)
700 : {
701 : rotation_angle = 180.;
702 : mesh_type = BOUNDARY_MESH;
703 : orientation = -180.;
704 : }
705 2445 : else if (j == 0 && i < (_pattern.size() - 1) / 2)
706 : {
707 : rotation_angle = 60.;
708 : mesh_type = BOUNDARY_MESH;
709 : orientation = -60.;
710 : }
711 2291 : else if (j == 0 && i > (_pattern.size() - 1) / 2)
712 : {
713 : rotation_angle = 120.;
714 : mesh_type = BOUNDARY_MESH;
715 : orientation = -120.;
716 : }
717 2137 : else if (j == _pattern[i].size() - 1 && i < (_pattern.size() - 1) / 2)
718 : {
719 : rotation_angle = -60.;
720 : mesh_type = BOUNDARY_MESH;
721 : orientation = 60.;
722 : }
723 1983 : else if (j == _pattern[i].size() - 1 && i > (_pattern.size() - 1) / 2)
724 : {
725 : rotation_angle = -120.;
726 : mesh_type = BOUNDARY_MESH;
727 : orientation = 120;
728 : }
729 : else
730 : on_periphery = false;
731 :
732 6563 : if (on_periphery)
733 : {
734 4734 : auto tmp_peripheral_mesh = dynamic_pointer_cast<ReplicatedMesh>(pattern_mesh.clone());
735 4734 : addPeripheralMesh(*tmp_peripheral_mesh,
736 : pattern,
737 : pitch_array.front(),
738 : extra_dist,
739 : num_sectors_per_side_array,
740 : peripheral_duct_intervals,
741 : rotation_angle,
742 : mesh_type);
743 :
744 4734 : if (extra_dist_shift != 0)
745 1320 : cutOffPolyDeform(*tmp_peripheral_mesh, orientation, y_max_0, y_max_n, y_min, mesh_type);
746 :
747 : // Reassign interface boundary ids
748 4734 : if (_use_interface_boundary_id_shift)
749 0 : reassignBoundaryIDs(*tmp_peripheral_mesh,
750 : _interface_boundary_id_shift_pattern[i][j],
751 : input_interface_boundary_ids[pattern]);
752 :
753 4734 : if (i == 0 && j == 0)
754 : out_mesh = std::move(tmp_peripheral_mesh);
755 : else
756 : {
757 : // Retrieve subdomain name map from the mesh to be stitched and insert it to the main
758 : // subdomain map
759 : const auto & increment_subdomain_map = tmp_peripheral_mesh->get_subdomain_name_map();
760 4099 : out_mesh->set_subdomain_name_map().insert(increment_subdomain_map.begin(),
761 : increment_subdomain_map.end());
762 :
763 4099 : MeshTools::Modification::translate(
764 4099 : *tmp_peripheral_mesh, deltax + j * input_pitch, deltay, 0);
765 4099 : out_mesh->stitch_meshes(*tmp_peripheral_mesh,
766 : OUTER_SIDESET_ID,
767 : OUTER_SIDESET_ID,
768 : TOLERANCE,
769 : /*clear_stitched_boundary_ids=*/true,
770 : /*verbose=*/false);
771 : }
772 :
773 : continue;
774 4734 : }
775 : }
776 :
777 3983 : if (_assign_control_drum_id && _pattern_boundary == "none" && _generate_core_metadata)
778 1829 : pattern_mesh.add_elem_integer(
779 : "control_drum_id",
780 : true,
781 : is_control_drum_array[pattern] ? control_drum_azimuthals.size() : 0);
782 :
783 : // Translate to correct position
784 3983 : MeshTools::Modification::translate(pattern_mesh, deltax + j * input_pitch, deltay, 0);
785 :
786 : // Define a reference map variable for subdomain map
787 : auto & main_subdomain_map = out_mesh->set_subdomain_name_map();
788 : // Retrieve subdomain name map from the mesh to be stitched and insert it to the main
789 : // subdomain map
790 : const auto & increment_subdomain_map = pattern_mesh.get_subdomain_name_map();
791 3983 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
792 : // Check if one SubdomainName is shared by more than one subdomain ids
793 : std::set<SubdomainName> main_subdomain_map_name_list;
794 22179 : for (auto const & id_name_pair : main_subdomain_map)
795 18196 : main_subdomain_map_name_list.emplace(id_name_pair.second);
796 3983 : if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
797 2 : paramError("inputs", "The input meshes contain subdomain name maps with conflicts.");
798 : // Shift interface boundary ids
799 3981 : if (_use_interface_boundary_id_shift)
800 54 : reassignBoundaryIDs(pattern_mesh,
801 : _interface_boundary_id_shift_pattern[i][j],
802 : input_interface_boundary_ids[pattern]);
803 :
804 3981 : out_mesh->stitch_meshes(pattern_mesh,
805 : OUTER_SIDESET_ID,
806 : OUTER_SIDESET_ID,
807 : TOLERANCE,
808 : /*clear_stitched_boundary_ids=*/false,
809 : /*verbose=*/false);
810 :
811 : // Translate back now that we've stitched so that anyone else that uses this mesh has it at
812 : // the origin
813 3981 : MeshTools::Modification::translate(pattern_mesh, -(deltax + j * input_pitch), -deltay, 0);
814 : // Roll back the changes in interface boundary ids for the same reason
815 3981 : if (_use_interface_boundary_id_shift)
816 54 : reassignBoundaryIDs(pattern_mesh,
817 : _interface_boundary_id_shift_pattern[i][j],
818 : input_interface_boundary_ids[pattern],
819 : true);
820 : }
821 : }
822 :
823 : // Check if OUTER_SIDESET_ID is really the external boundary. Correct if needed.
824 939 : auto side_list = out_mesh->get_boundary_info().build_side_list();
825 454785 : for (auto & sl : side_list)
826 : {
827 453846 : if (std::get<2>(sl) == OUTER_SIDESET_ID)
828 123456 : if (out_mesh->elem_ptr(std::get<0>(sl))->neighbor_ptr(std::get<1>(sl)) != nullptr)
829 54814 : out_mesh->get_boundary_info().remove_side(
830 54814 : out_mesh->elem_ptr(std::get<0>(sl)), std::get<1>(sl), std::get<2>(sl));
831 : }
832 939 : out_mesh->get_boundary_info().clear_boundary_node_ids();
833 :
834 939 : out_mesh->get_boundary_info().build_node_list_from_side_list();
835 939 : const auto node_list = out_mesh->get_boundary_info().build_node_list();
836 :
837 : std::vector<Real> bd_x_list;
838 : std::vector<Real> bd_y_list;
839 : std::vector<std::pair<Real, dof_id_type>> node_azi_list;
840 939 : const Point origin_pt = MooseMeshUtils::meshCentroidCalculator(*out_mesh);
841 939 : const Real origin_x = origin_pt(0);
842 939 : const Real origin_y = origin_pt(1);
843 :
844 939 : MeshTools::Modification::translate(*out_mesh, -origin_x, -origin_y, 0);
845 :
846 939 : if (_uniform_mesh_on_sides)
847 : {
848 : const Real azi_tol = 1E-8;
849 32229 : for (unsigned int i = 0; i < node_list.size(); ++i)
850 : {
851 32184 : if (std::get<1>(node_list[i]) == OUTER_SIDESET_ID)
852 : {
853 : node_azi_list.push_back(
854 9288 : std::make_pair(atan2(out_mesh->node_ref(std::get<0>(node_list[i]))(1),
855 4644 : out_mesh->node_ref(std::get<0>(node_list[i]))(0)) *
856 4644 : 180.0 / M_PI,
857 : std::get<0>(node_list[i])));
858 : // correct the last node's angle value (180.0) if it becomes a negative value.
859 4644 : if (node_azi_list.back().first + 180.0 <= azi_tol)
860 24 : node_azi_list.back().first = 180;
861 : }
862 : }
863 45 : std::sort(node_azi_list.begin(), node_azi_list.end());
864 45 : const unsigned int side_intervals = node_azi_list.size() / HEXAGON_NUM_SIDES;
865 315 : for (unsigned int i = 0; i < HEXAGON_NUM_SIDES; i++)
866 : {
867 4914 : for (unsigned int j = 1; j <= side_intervals; j++)
868 : {
869 4644 : Real azi_corr_tmp = atan2((Real)j * 2.0 / (Real)side_intervals - 1.0, std::sqrt(3.0));
870 4644 : Real x_tmp = _pattern_pitch / 2.0;
871 4644 : Real y_tmp = x_tmp * std::tan(azi_corr_tmp);
872 4644 : nodeCoordRotate(x_tmp, y_tmp, (Real)i * 60.0 - 150.0);
873 4644 : Point p_tmp = Point(x_tmp, y_tmp, 0.0);
874 4644 : out_mesh->add_point(p_tmp, node_azi_list[i * side_intervals + j - 1].second);
875 : }
876 : }
877 :
878 : // if quadratic elements are used, additional nodes need to be adjusted based on the new
879 : // boundary node locations. adjust side mid-edge nodes to the midpoints of the corner
880 : // points, and if QUAD9, adjust center point to new centroid.
881 45 : if (_boundary_quad_elem_type != QUAD_ELEM_TYPE::QUAD4)
882 18 : adjustPeripheralQuadraticElements(*out_mesh, _boundary_quad_elem_type);
883 : }
884 :
885 939 : MeshTools::Modification::rotate(*out_mesh, _rotate_angle, 0.0, 0.0);
886 :
887 : // This combination of input parameters is usually used for core mesh generation by stitching
888 : // assembly meshes together.
889 939 : if (_pattern_boundary == "none" && _generate_core_metadata)
890 : {
891 : const Real azi_tol = 1E-8;
892 : std::vector<std::tuple<Real, Point, std::vector<Real>, dof_id_type>> control_drum_tmp;
893 : std::vector<dof_id_type> control_drum_id_sorted;
894 403 : for (unsigned int i = 0; i < control_drum_positions_x.size(); ++i)
895 : {
896 153 : control_drum_positions_x[i] -= origin_x;
897 153 : control_drum_positions_y[i] -= origin_y;
898 153 : nodeCoordRotate(control_drum_positions_x[i], control_drum_positions_y[i], _rotate_angle);
899 153 : Real cd_angle = atan2(control_drum_positions_y[i], control_drum_positions_x[i]);
900 :
901 3105 : for (unsigned int j = 0; j < control_drum_azimuthals[i].size(); j++)
902 : {
903 2952 : control_drum_azimuthals[i][j] += _rotate_angle;
904 2952 : control_drum_azimuthals[i][j] =
905 2952 : atan2(std::sin(control_drum_azimuthals[i][j] / 180.0 * M_PI),
906 2952 : std::cos(control_drum_azimuthals[i][j] / 180.0 * M_PI)) /
907 2952 : M_PI * 180.0; // quick way to move to -M_PI to M_PI
908 : }
909 153 : std::sort(control_drum_azimuthals[i].begin(), control_drum_azimuthals[i].end());
910 :
911 153 : if (std::abs(cd_angle) < azi_tol)
912 : cd_angle = 0;
913 153 : else if (cd_angle < 0.0)
914 81 : cd_angle += 2 * M_PI;
915 153 : control_drum_tmp.push_back(
916 153 : std::make_tuple(cd_angle,
917 153 : Point(control_drum_positions_x[i], control_drum_positions_y[i], 0.0),
918 : control_drum_azimuthals[i],
919 153 : i + 1)); // control drum index to help sort control_drum_id
920 : }
921 250 : std::sort(control_drum_tmp.begin(), control_drum_tmp.end());
922 : std::vector<Point> control_drum_positions;
923 : std::vector<Real> control_drum_angles;
924 : std::vector<std::vector<Real>> control_drums_azimuthal_meta;
925 403 : for (unsigned int i = 0; i < control_drum_tmp.size(); ++i)
926 : {
927 153 : control_drum_positions.push_back(std::get<1>(control_drum_tmp[i]));
928 153 : control_drum_angles.push_back(std::get<0>(control_drum_tmp[i]));
929 153 : control_drums_azimuthal_meta.push_back(std::get<2>(control_drum_tmp[i]));
930 153 : control_drum_id_sorted.push_back(std::get<3>(control_drum_tmp[i]));
931 : }
932 250 : setMeshProperty("control_drum_positions", control_drum_positions);
933 250 : setMeshProperty("control_drum_angles", control_drum_angles);
934 250 : setMeshProperty("control_drums_azimuthal_meta", control_drums_azimuthal_meta);
935 :
936 250 : if (_assign_control_drum_id)
937 : {
938 145 : unsigned int id = out_mesh->get_elem_integer_index("control_drum_id");
939 648378 : for (const auto & elem : out_mesh->element_ptr_range())
940 : {
941 324044 : dof_id_type unsorted_control_drum_id = elem->get_extra_integer(id);
942 324044 : if (unsorted_control_drum_id != 0)
943 : {
944 5904 : auto sorted_iter = std::find(control_drum_id_sorted.begin(),
945 : control_drum_id_sorted.end(),
946 : unsorted_control_drum_id);
947 5904 : elem->set_extra_integer(id,
948 5904 : std::distance(control_drum_id_sorted.begin(), sorted_iter) + 1);
949 : }
950 145 : }
951 : }
952 :
953 250 : if (_generate_control_drum_positions_file)
954 : {
955 33 : std::string position_file_name = getMeshProperty<std::string>("position_file_name", name());
956 33 : std::ofstream pos_file(position_file_name);
957 186 : for (unsigned int i = 0; i < control_drum_positions.size(); ++i)
958 153 : pos_file << control_drum_positions[i](0) << " " << control_drum_positions[i](1) << " 0.0\n";
959 33 : pos_file.close();
960 33 : }
961 250 : }
962 :
963 : // before return, add reporting IDs if _use_reporting_id is set true
964 939 : if (_use_reporting_id)
965 528 : addReportingIDs(*out_mesh, meshes);
966 :
967 : // Assign customized peripheral block ids and names
968 939 : if (!_peripheral_block_ids.empty())
969 422198 : for (const auto & elem : out_mesh->active_element_ptr_range())
970 462840 : for (subdomain_id_type i = PERIPHERAL_ID_SHIFT; i <= PERIPHERAL_ID_SHIFT + _duct_sizes.size();
971 : i++)
972 298512 : if (elem->subdomain_id() == i)
973 : {
974 46308 : elem->subdomain_id() = _peripheral_block_ids[i - PERIPHERAL_ID_SHIFT];
975 46308 : break;
976 463 : }
977 939 : if (!_peripheral_block_names.empty())
978 : {
979 1134 : for (unsigned i = 0; i < _peripheral_block_names.size(); i++)
980 680 : out_mesh->subdomain_name(_peripheral_block_ids.empty() ? (PERIPHERAL_ID_SHIFT + i)
981 : : _peripheral_block_ids[i]) =
982 680 : _peripheral_block_names[i];
983 : }
984 : // Assign customized outer surface boundary id and name
985 939 : if (_external_boundary_id > 0)
986 409 : MooseMesh::changeBoundaryId(*out_mesh, OUTER_SIDESET_ID, _external_boundary_id, false);
987 939 : if (!_external_boundary_name.empty())
988 : {
989 409 : out_mesh->get_boundary_info().sideset_name(
990 409 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
991 409 : _external_boundary_name;
992 409 : out_mesh->get_boundary_info().nodeset_name(
993 409 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
994 : _external_boundary_name;
995 : }
996 : // Merge the boundary name maps of all the input meshed to generate the output mesh's boundary
997 : // name maps
998 : auto & new_sideset_map = out_mesh->get_boundary_info().set_sideset_name_map();
999 : auto & new_nodeset_map = out_mesh->get_boundary_info().set_nodeset_name_map();
1000 2614 : for (unsigned int i = 0; i < meshes.size(); i++)
1001 : {
1002 : const auto input_sideset_map = meshes[i]->get_boundary_info().get_sideset_name_map();
1003 1675 : new_sideset_map.insert(input_sideset_map.begin(), input_sideset_map.end());
1004 : const auto input_nodeset_map = meshes[i]->get_boundary_info().get_nodeset_name_map();
1005 1675 : new_nodeset_map.insert(input_nodeset_map.begin(), input_nodeset_map.end());
1006 : }
1007 :
1008 : // set mesh metadata related with interface boundary ids
1009 : const std::set<boundary_id_type> boundary_ids = out_mesh->get_boundary_info().get_boundary_ids();
1010 : const std::set<boundary_id_type> interface_boundary_ids = getInterfaceBoundaryIDs(
1011 : _pattern,
1012 939 : _interface_boundary_id_shift_pattern,
1013 : boundary_ids,
1014 : input_interface_boundary_ids,
1015 939 : _use_interface_boundary_id_shift,
1016 939 : _create_inward_interface_boundaries || _create_outward_interface_boundaries,
1017 1878 : extra_dist.size());
1018 939 : if (interface_boundary_ids.size() > 0)
1019 : {
1020 564 : setMeshProperty("interface_boundaries", true);
1021 1128 : setMeshProperty("interface_boundary_ids", interface_boundary_ids);
1022 : }
1023 :
1024 : out_mesh->set_isnt_prepared();
1025 939 : auto mesh = dynamic_pointer_cast<MeshBase>(out_mesh);
1026 939 : return mesh;
1027 1878 : }
1028 :
1029 : void
1030 4734 : PatternedHexMeshGenerator::addPeripheralMesh(
1031 : ReplicatedMesh & mesh,
1032 : const unsigned int pattern,
1033 : const Real pitch,
1034 : const std::vector<Real> & extra_dist,
1035 : const std::vector<unsigned int> & num_sectors_per_side_array,
1036 : const std::vector<unsigned int> & peripheral_duct_intervals,
1037 : const Real rotation_angle,
1038 : const unsigned int mesh_type)
1039 : {
1040 : std::vector<std::pair<Real, Real>> positions_inner;
1041 : std::vector<std::pair<Real, Real>> d_positions_outer;
1042 :
1043 : std::vector<std::vector<unsigned int>> peripheral_point_index;
1044 : std::vector<std::pair<Real, Real>> sub_positions_inner;
1045 : std::vector<std::pair<Real, Real>> sub_d_positions_outer;
1046 :
1047 4734 : if (mesh_type == CORNER_MESH)
1048 : // corner mesh has three sides that need peripheral meshes.
1049 : // each element has three sub-elements, representing beginning, middle, and ending azimuthal
1050 : // points
1051 19050 : peripheral_point_index = {{0, 1, 2}, {2, 3, 4}, {4, 5, 6}};
1052 : else
1053 : // side mesh has two sides that need peripheral meshes.
1054 3696 : peripheral_point_index = {{0, 1, 2}, {2, 7, 8}};
1055 :
1056 : // extra_dist includes background and ducts.
1057 : // Loop to calculate the positions of the boundaries.
1058 12396 : for (unsigned int i = 0; i < extra_dist.size(); i++)
1059 : {
1060 7662 : positionSetup(positions_inner,
1061 : d_positions_outer,
1062 2928 : i == 0 ? 0.0 : extra_dist[i - 1],
1063 : extra_dist[i],
1064 : pitch,
1065 : i);
1066 :
1067 : // Loop for all applicable sides that need peripheral mesh (3 for corner and 2 for edge)
1068 28908 : for (unsigned int peripheral_index = 0; peripheral_index < peripheral_point_index.size();
1069 : peripheral_index++)
1070 : {
1071 : // Loop for beginning, middle and ending positions of a side
1072 84984 : for (unsigned int vector_index = 0; vector_index < 3; vector_index++)
1073 : {
1074 63738 : sub_positions_inner.push_back(
1075 63738 : positions_inner[peripheral_point_index[peripheral_index][vector_index]]);
1076 63738 : sub_d_positions_outer.push_back(
1077 63738 : d_positions_outer[peripheral_point_index[peripheral_index][vector_index]]);
1078 : }
1079 21246 : auto meshp0 = buildSimplePeripheral(num_sectors_per_side_array[pattern],
1080 : peripheral_duct_intervals[i],
1081 : sub_positions_inner,
1082 : sub_d_positions_outer,
1083 : i,
1084 : _boundary_quad_elem_type,
1085 21246 : _create_inward_interface_boundaries,
1086 21246 : (i != extra_dist.size() - 1) &&
1087 21246 : _create_outward_interface_boundaries);
1088 :
1089 : // The other_mesh must be prepared before stitching
1090 21246 : meshp0->prepare_for_use();
1091 :
1092 : // rotate the peripheral mesh to the desired side of the hexagon.
1093 21246 : MeshTools::Modification::rotate(*meshp0, rotation_angle, 0, 0);
1094 21246 : mesh.stitch_meshes(*meshp0, OUTER_SIDESET_ID, OUTER_SIDESET_ID, TOLERANCE, true, false);
1095 21246 : sub_positions_inner.resize(0);
1096 21246 : sub_d_positions_outer.resize(0);
1097 21246 : }
1098 : }
1099 14202 : }
1100 :
1101 : void
1102 7662 : PatternedHexMeshGenerator::positionSetup(std::vector<std::pair<Real, Real>> & positions_inner,
1103 : std::vector<std::pair<Real, Real>> & d_positions_outer,
1104 : const Real extra_dist_in,
1105 : const Real extra_dist_out,
1106 : const Real pitch,
1107 : const unsigned int radial_index) const
1108 : {
1109 7662 : positions_inner.resize(0);
1110 7662 : d_positions_outer.resize(0);
1111 :
1112 : // Nine sets of positions are generated here, as shown below.
1113 : // CORNER MESH Peripheral {0 1 2}, {2 3 4} and {4 5 6}
1114 : // 3 2 1 0
1115 : // \ : : :
1116 : // \ : : :
1117 : // 4. . : :
1118 : // ` . :
1119 : // 5. | |
1120 : // ` | |
1121 : // 6. | |
1122 : // ` | |
1123 : // . .
1124 : // .
1125 : //
1126 : // EDGE MESH Peripheral {0 1 2} and {2 7 8}
1127 : // 8 7 2 1 0
1128 : // : : : : :
1129 : // : : : : :
1130 : // : : : :
1131 : // : :
1132 : // | |
1133 : // | |
1134 : // | |
1135 : // | |
1136 : // . .
1137 : // .
1138 :
1139 : positions_inner.push_back(
1140 15324 : std::make_pair(-pitch / 2.0,
1141 7662 : std::sqrt(3.0) * pitch / 6.0 +
1142 7662 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 6.0 + extra_dist_in : 0.0)));
1143 15324 : positions_inner.push_back(std::make_pair(
1144 7662 : -pitch / 4.0,
1145 7662 : std::sqrt(3.0) * pitch / 4.0 +
1146 7662 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in : 0.0)));
1147 15324 : positions_inner.push_back(std::make_pair(
1148 12396 : 0.0, std::sqrt(3.0) * pitch / 3.0 + (radial_index != 0 ? extra_dist_in : 0.0)));
1149 15324 : positions_inner.push_back(std::make_pair(
1150 7662 : pitch / 4.0 + (radial_index != 0 ? pitch / 12.0 + std::sqrt(3.0) * extra_dist_in / 3.0 : 0.0),
1151 : std::sqrt(3.0) * pitch / 4.0 +
1152 : (radial_index != 0 ? pitch * std::sqrt(3.0) / 12.0 + extra_dist_in : 0.0)));
1153 15324 : positions_inner.push_back(std::make_pair(
1154 7662 : pitch / 2.0 + (radial_index != 0 ? std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1155 7662 : std::sqrt(3.0) * pitch / 6.0 + (radial_index != 0 ? extra_dist_in / 2.0 : 0.0)));
1156 15324 : positions_inner.push_back(std::make_pair(
1157 7662 : pitch / 2.0 + (radial_index != 0 ? pitch / 8.0 + std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1158 7662 : 0.0 + (radial_index != 0 ? std::sqrt(3.0) * pitch / 24.0 + extra_dist_in / 2.0 : 0.0)));
1159 15324 : positions_inner.push_back(std::make_pair(
1160 7662 : pitch / 2.0 + (radial_index != 0 ? pitch / 4.0 + std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1161 7662 : -std::sqrt(3.0) * pitch / 6.0 +
1162 7662 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in / 2.0 : 0.0)));
1163 7662 : positions_inner.push_back(std::make_pair(
1164 : pitch / 4.0,
1165 : std::sqrt(3.0) * pitch / 4.0 +
1166 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in : 0.0)));
1167 : positions_inner.push_back(
1168 7662 : std::make_pair(pitch / 2.0,
1169 : std::sqrt(3.0) * pitch / 6.0 +
1170 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 6.0 + extra_dist_in : 0.0)));
1171 :
1172 : d_positions_outer.push_back(
1173 15324 : std::make_pair(0.0,
1174 7662 : std::sqrt(3.0) * pitch / 6.0 + extra_dist_out -
1175 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 6.0 + extra_dist_in : 0.0)));
1176 15324 : d_positions_outer.push_back(std::make_pair(
1177 : 0.0,
1178 7662 : std::sqrt(3.0) * pitch / 12.0 + extra_dist_out -
1179 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in : 0.0)));
1180 : d_positions_outer.push_back(
1181 7662 : std::make_pair(0.0, extra_dist_out - (radial_index != 0 ? extra_dist_in : 0.0)));
1182 15324 : d_positions_outer.push_back(std::make_pair(
1183 7662 : pitch / 12.0 + std::sqrt(3.0) * extra_dist_out / 3.0 -
1184 : (radial_index != 0 ? pitch / 12.0 + std::sqrt(3.0) * extra_dist_in / 3.0 : 0.0),
1185 : pitch * std::sqrt(3.0) / 12.0 + extra_dist_out -
1186 : (radial_index != 0 ? pitch * std::sqrt(3.0) / 12.0 + extra_dist_in : 0.0)));
1187 : d_positions_outer.push_back(
1188 15324 : std::make_pair(std::sqrt(3.0) * extra_dist_out / 2.0 -
1189 : (radial_index != 0 ? std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1190 7662 : extra_dist_out / 2.0 - (radial_index != 0 ? extra_dist_in / 2.0 : 0.0)));
1191 15324 : d_positions_outer.push_back(std::make_pair(
1192 7662 : pitch / 8.0 + std::sqrt(3.0) * extra_dist_out / 2.0 -
1193 : (radial_index != 0 ? pitch / 8.0 + std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1194 7662 : std::sqrt(3.0) * pitch / 24.0 + extra_dist_out / 2.0 -
1195 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 24.0 + extra_dist_in / 2.0 : 0.0)));
1196 15324 : d_positions_outer.push_back(std::make_pair(
1197 7662 : pitch / 4.0 + std::sqrt(3.0) * extra_dist_out / 2.0 -
1198 : (radial_index != 0 ? pitch / 4.0 + std::sqrt(3.0) * extra_dist_in / 2.0 : 0.0),
1199 7662 : std::sqrt(3.0) * pitch / 12.0 + extra_dist_out / 2.0 -
1200 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in / 2.0 : 0.0)));
1201 7662 : d_positions_outer.push_back(std::make_pair(
1202 : 0.0,
1203 : std::sqrt(3.0) * pitch / 12.0 + extra_dist_out -
1204 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 12.0 + extra_dist_in : 0.0)));
1205 : d_positions_outer.push_back(
1206 7662 : std::make_pair(0.0,
1207 : std::sqrt(3.0) * pitch / 6.0 + extra_dist_out -
1208 : (radial_index != 0 ? std::sqrt(3.0) * pitch / 6.0 + extra_dist_in : 0.0)));
1209 7662 : }
1210 :
1211 : void
1212 528 : PatternedHexMeshGenerator::addReportingIDs(
1213 : MeshBase & mesh, const std::vector<std::unique_ptr<ReplicatedMesh>> & from_meshes) const
1214 : {
1215 528 : const unsigned int num_reporting_ids = _reporting_id_names.size();
1216 1083 : for (unsigned int i = 0; i < num_reporting_ids; ++i)
1217 : {
1218 555 : const std::string element_id_name = _reporting_id_names[i];
1219 : unsigned int extra_id_index;
1220 555 : if (!mesh.has_elem_integer(element_id_name))
1221 1110 : extra_id_index = mesh.add_elem_integer(element_id_name);
1222 : else
1223 : {
1224 0 : extra_id_index = mesh.get_elem_integer_index(element_id_name);
1225 0 : paramWarning(
1226 : "id_name", "An element integer with the name '", element_id_name, "' already exists");
1227 : }
1228 :
1229 : // assign reporting IDs to individual elements
1230 : // NOTE: background block id should be set "PERIPHERAL_ID_SHIFT" because this function is called
1231 : // before assigning the user-defined background block id
1232 : std::set<subdomain_id_type> background_block_ids =
1233 1665 : (isParamValid("background_block_id")) ? std::set<subdomain_id_type>({PERIPHERAL_ID_SHIFT})
1234 555 : : std::set<subdomain_id_type>();
1235 :
1236 : const bool using_manual_id =
1237 555 : (_assign_types[i] == ReportingIDGeneratorUtils::AssignType::manual);
1238 555 : ReportingIDGeneratorUtils::assignReportingIDs(mesh,
1239 : extra_id_index,
1240 : _assign_types[i],
1241 555 : _use_exclude_id,
1242 555 : _exclude_ids,
1243 555 : _pattern_boundary == "hexagon",
1244 : background_block_ids,
1245 : from_meshes,
1246 : _pattern,
1247 : (using_manual_id)
1248 1110 : ? _id_patterns.at(element_id_name)
1249 : : std::vector<std::vector<dof_id_type>>());
1250 : }
1251 528 : }
|
