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 "PeripheralRingMeshGenerator.h"
11 :
12 : #include "MooseMeshUtils.h"
13 : #include "MooseUtils.h"
14 : #include "LinearInterpolation.h"
15 : #include "FillBetweenPointVectorsTools.h"
16 : #include "PolygonalMeshGenerationUtils.h"
17 : #include "libmesh/int_range.h"
18 :
19 : // C++ includes
20 : #include <cmath> // for atan2
21 :
22 : registerMooseObject("ReactorApp", PeripheralRingMeshGenerator);
23 :
24 : InputParameters
25 240 : PeripheralRingMeshGenerator::validParams()
26 : {
27 240 : InputParameters params = PolygonMeshGeneratorBase::validParams();
28 480 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh to be modified.");
29 480 : params.addParam<bool>(
30 : "force_input_centroid_as_center",
31 480 : false,
32 : "Whether to enforce use of the centroid position of the input mesh as the "
33 : "center of the peripheral ring by translating the input mesh to the origin.");
34 :
35 720 : params.addRangeCheckedParam<unsigned int>(
36 : "peripheral_layer_num",
37 480 : 3,
38 : "peripheral_layer_num>0",
39 : "The radial layers of the peripheral ring to be added.");
40 720 : params.addRangeCheckedParam<Real>(
41 : "peripheral_radial_bias",
42 480 : 1.0,
43 : "peripheral_radial_bias>0",
44 : "Value used to create biasing in radial meshing for peripheral ring region.");
45 720 : params.addRangeCheckedParam<Real>(
46 : "peripheral_inner_boundary_layer_width",
47 480 : 0.0,
48 : "peripheral_inner_boundary_layer_width>=0",
49 : "Width of peripheral ring region that is assigned to be the inner boundary layer.");
50 720 : params.addRangeCheckedParam<unsigned int>(
51 : "peripheral_inner_boundary_layer_intervals",
52 480 : 1,
53 : "peripheral_inner_boundary_layer_intervals>0",
54 : "Number of radial intervals of the peripheral ring inner boundary layer");
55 720 : params.addRangeCheckedParam<Real>(
56 : "peripheral_inner_boundary_layer_bias",
57 480 : 1.0,
58 : "peripheral_inner_boundary_layer_bias>0",
59 : "Growth factor used for mesh biasing of the peripheral ring inner boundary layer.");
60 720 : params.addRangeCheckedParam<Real>(
61 : "peripheral_outer_boundary_layer_width",
62 480 : 0.0,
63 : "peripheral_outer_boundary_layer_width>=0",
64 : "Width of peripheral ring region that is assigned to be the outer boundary layer.");
65 720 : params.addRangeCheckedParam<unsigned int>(
66 : "peripheral_outer_boundary_layer_intervals",
67 480 : 1,
68 : "peripheral_outer_boundary_layer_intervals>0",
69 : "Number of radial intervals of the peripheral ring outer boundary layer");
70 720 : params.addRangeCheckedParam<Real>(
71 : "peripheral_outer_boundary_layer_bias",
72 480 : 1.0,
73 : "peripheral_outer_boundary_layer_bias>0",
74 : "Growth factor used for mesh biasing of the peripheral ring outer boundary layer.");
75 480 : params.addRequiredRangeCheckedParam<Real>("peripheral_ring_radius",
76 : "peripheral_ring_radius>0",
77 : "Radius of the peripheral ring to be added.");
78 480 : params.addParam<bool>(
79 : "preserve_volumes",
80 480 : true,
81 : "Whether the volume of the peripheral region is preserved by fixing the radius.");
82 480 : params.addRequiredParam<BoundaryName>("input_mesh_external_boundary",
83 : "The external boundary of the input mesh.");
84 480 : params.addRequiredParam<subdomain_id_type>(
85 : "peripheral_ring_block_id", "The block id assigned to the created peripheral layer.");
86 480 : params.addParam<SubdomainName>("peripheral_ring_block_name",
87 : "The block name assigned to the created peripheral layer.");
88 480 : params.addRangeCheckedParam<boundary_id_type>("external_boundary_id",
89 : "external_boundary_id>0",
90 : "Optional customized external boundary id.");
91 480 : params.addParam<BoundaryName>(
92 : "external_boundary_name", "", "Optional customized external boundary name.");
93 480 : params.addParamNamesToGroup(
94 : "peripheral_radial_bias peripheral_inner_boundary_layer_width "
95 : "peripheral_inner_boundary_layer_intervals peripheral_inner_boundary_layer_bias "
96 : "peripheral_outer_boundary_layer_width peripheral_outer_boundary_layer_intervals "
97 : "peripheral_outer_boundary_layer_bias",
98 : "Mesh Boundary Layer and Biasing Options");
99 240 : params.addClassDescription("This PeripheralRingMeshGenerator object adds a circular peripheral "
100 : "region to the input mesh.");
101 :
102 240 : return params;
103 0 : }
104 :
105 115 : PeripheralRingMeshGenerator::PeripheralRingMeshGenerator(const InputParameters & parameters)
106 : : PolygonMeshGeneratorBase(parameters),
107 0 : _input_name(getParam<MeshGeneratorName>("input")),
108 230 : _force_input_centroid_as_center(getParam<bool>("force_input_centroid_as_center")),
109 230 : _peripheral_layer_num(getParam<unsigned int>("peripheral_layer_num")),
110 230 : _peripheral_radial_bias(getParam<Real>("peripheral_radial_bias")),
111 345 : _peripheral_inner_boundary_layer_params(
112 115 : {getParam<Real>("peripheral_inner_boundary_layer_width"),
113 : 0.0,
114 230 : getParam<Real>("peripheral_inner_boundary_layer_width") > 0.0
115 242 : ? getParam<unsigned int>("peripheral_inner_boundary_layer_intervals")
116 : : 0,
117 230 : getParam<Real>("peripheral_inner_boundary_layer_bias")}),
118 345 : _peripheral_outer_boundary_layer_params(
119 115 : {getParam<Real>("peripheral_outer_boundary_layer_width"),
120 : 0.0,
121 230 : getParam<Real>("peripheral_outer_boundary_layer_width") > 0.0
122 242 : ? getParam<unsigned int>("peripheral_outer_boundary_layer_intervals")
123 : : 0,
124 230 : getParam<Real>("peripheral_outer_boundary_layer_bias")}),
125 230 : _peripheral_ring_radius(getParam<Real>("peripheral_ring_radius")),
126 230 : _preserve_volumes(getParam<bool>("preserve_volumes")),
127 115 : _input_mesh_external_boundary(getParam<BoundaryName>("input_mesh_external_boundary")),
128 230 : _peripheral_ring_block_id(getParam<subdomain_id_type>("peripheral_ring_block_id")),
129 319 : _peripheral_ring_block_name(isParamValid("peripheral_ring_block_name")
130 115 : ? getParam<SubdomainName>("peripheral_ring_block_name")
131 : : (SubdomainName) ""),
132 230 : _external_boundary_id(isParamValid("external_boundary_id")
133 115 : ? getParam<boundary_id_type>("external_boundary_id")
134 : : 0),
135 115 : _external_boundary_name(getParam<BoundaryName>("external_boundary_name")),
136 230 : _input(getMeshByName(_input_name))
137 : {
138 115 : declareMeshProperty<bool>("hexagon_peripheral_trimmability", false);
139 115 : declareMeshProperty<bool>("hexagon_center_trimmability", false);
140 115 : declareMeshProperty<bool>("square_peripheral_trimmability", false);
141 115 : declareMeshProperty<bool>("square_center_trimmability", false);
142 115 : }
143 :
144 : std::unique_ptr<MeshBase>
145 115 : PeripheralRingMeshGenerator::generate()
146 : {
147 230 : if (hasMeshProperty<bool>("hexagon_center_trimmability", _input_name))
148 18 : setMeshProperty("hexagon_center_trimmability",
149 : getMeshProperty<bool>("hexagon_center_trimmability", _input_name));
150 230 : if (hasMeshProperty<bool>("square_center_trimmability", _input_name))
151 18 : setMeshProperty("square_center_trimmability",
152 : getMeshProperty<bool>("square_center_trimmability", _input_name));
153 :
154 : // Need ReplicatedMesh for stitching
155 115 : auto input_mesh = dynamic_cast<ReplicatedMesh *>(_input.get());
156 115 : if (!input_mesh)
157 0 : paramError("input", "Input is not a replicated mesh, which is required.");
158 115 : if (*(input_mesh->elem_dimensions().begin()) != 2 ||
159 113 : *(input_mesh->elem_dimensions().rbegin()) != 2)
160 2 : paramError("input", "Only 2D meshes are supported.");
161 :
162 113 : _input_mesh_external_bid =
163 113 : MooseMeshUtils::getBoundaryID(_input_mesh_external_boundary, *input_mesh);
164 113 : if (!MooseMeshUtils::hasBoundaryName(*input_mesh, _input_mesh_external_boundary))
165 2 : paramError("input_mesh_external_boundary",
166 : "External boundary does not exist in the input mesh");
167 : // We check the element types of input mesh's external boundary here.
168 : // We support both linear and quadratic sides (i.e., EDGE2 and EDGE3), but we cannot support mixed
169 : // sides
170 111 : auto side_list_tmp = input_mesh->get_boundary_info().build_side_list();
171 : bool linear_side = false;
172 : bool quadratic_side = false;
173 295311 : for (const auto & side : side_list_tmp)
174 : {
175 295200 : if (std::get<2>(side) == _input_mesh_external_bid)
176 : {
177 18368 : const auto etype = input_mesh->elem_ptr(std::get<0>(side))->side_type(std::get<1>(side));
178 18368 : if (etype == EDGE2)
179 : linear_side = true;
180 168 : else if (etype == EDGE3)
181 : quadratic_side = true;
182 : else
183 0 : paramError("input", "Input contains elements that are not supported by this generator.");
184 : }
185 : }
186 111 : if (linear_side && quadratic_side)
187 2 : paramError("input", "Input contains mixed element types on the external boundary.");
188 109 : const unsigned int order = linear_side ? 1 : 2;
189 : if (order == 2)
190 : {
191 5 : _peripheral_outer_boundary_layer_params.bias =
192 5 : std::sqrt(_peripheral_outer_boundary_layer_params.bias);
193 5 : _peripheral_inner_boundary_layer_params.bias =
194 5 : std::sqrt(_peripheral_inner_boundary_layer_params.bias);
195 5 : _peripheral_outer_boundary_layer_params.intervals *= 2;
196 5 : _peripheral_inner_boundary_layer_params.intervals *= 2;
197 5 : _peripheral_radial_bias = std::sqrt(_peripheral_radial_bias);
198 : }
199 :
200 : // Move the centroid of the mesh to (0, 0, 0) if input centroid is enforced to be the ring center.
201 109 : const Point input_mesh_centroid = MooseMeshUtils::meshCentroidCalculator(*input_mesh);
202 109 : if (_force_input_centroid_as_center)
203 5 : MeshTools::Modification::translate(
204 : *input_mesh, -input_mesh_centroid(0), -input_mesh_centroid(1), -input_mesh_centroid(2));
205 :
206 : // Use CoM of the input mesh as its origin for azimuthal calculation
207 : const Point origin_pt =
208 109 : _force_input_centroid_as_center ? Point(0.0, 0.0, 0.0) : input_mesh_centroid;
209 : // Vessel for containing maximum radius of the boundary nodes
210 : Real max_input_mesh_node_radius;
211 :
212 : // Calculate biasing terms
213 : // For 2nd order mesh, as we need the mid points, the bias factor is square rooted.
214 : const auto main_peripheral_bias_terms =
215 109 : biasTermsCalculator(_peripheral_radial_bias, _peripheral_layer_num * order);
216 : const auto inner_peripheral_bias_terms =
217 : biasTermsCalculator(_peripheral_inner_boundary_layer_params.bias,
218 109 : _peripheral_inner_boundary_layer_params.intervals);
219 : // It is easier to create outer boundary layer inversely (inwards). Thus, 1.0 / bias is used here.
220 : // However, the input parameter definition is not affected.
221 : const auto outer_peripheral_bias_terms =
222 109 : biasTermsCalculator(1.0 / _peripheral_outer_boundary_layer_params.bias,
223 109 : _peripheral_outer_boundary_layer_params.intervals);
224 :
225 109 : const unsigned int total_peripheral_layer_num =
226 109 : _peripheral_inner_boundary_layer_params.intervals + _peripheral_layer_num * order +
227 109 : _peripheral_outer_boundary_layer_params.intervals;
228 :
229 : // Collect the external boundary nodes of the input mesh using the utility function
230 : std::vector<dof_id_type> boundary_ordered_node_list;
231 : try
232 : {
233 109 : FillBetweenPointVectorsTools::isBoundarySimpleClosedLoop(*input_mesh,
234 : max_input_mesh_node_radius,
235 : boundary_ordered_node_list,
236 : origin_pt,
237 109 : _input_mesh_external_bid);
238 : }
239 6 : catch (MooseException & e)
240 : {
241 6 : if (((std::string)e.what()).compare("The node list provided has more than one segments.") == 0)
242 2 : paramError("input_mesh_external_boundary",
243 : "This mesh generator does not work for the provided external boundary as it has "
244 : "more than one segments.");
245 : else
246 4 : paramError("input_mesh_external_boundary", e.what());
247 0 : }
248 :
249 103 : if (max_input_mesh_node_radius >= _peripheral_ring_radius)
250 2 : paramError(
251 : "peripheral_ring_radius",
252 : "The peripheral ring to be generated must be large enough to cover the entire input mesh.");
253 101 : if (!FillBetweenPointVectorsTools::isExternalBoundary(*input_mesh, _input_mesh_external_bid))
254 2 : paramError("input_mesh_external_boundary",
255 : "The boundary provided is not an external boundary.");
256 :
257 : std::vector<Point> input_ext_bd_pts;
258 : std::vector<Point> output_ext_bd_pts;
259 : std::vector<std::tuple<Real, Point, Point>> azi_points;
260 : std::vector<Real> azi_array;
261 : Real tmp_azi(0.0);
262 99 : auto mesh = buildReplicatedMesh(2);
263 :
264 : // Node counter for the external boundary
265 : unsigned int input_ext_node_num = 0;
266 :
267 : // As the node list collected before is a simple closed loop, the last node is the same as the
268 : // first node. We remove it here.
269 : boundary_ordered_node_list.pop_back();
270 : // Loop over all the boundary nodes
271 : // For quadratic sides, the middle points are included automatically
272 12071 : for (const auto i : index_range(boundary_ordered_node_list))
273 : {
274 11972 : input_ext_node_num++;
275 : // Define nodes on the inner and outer boundaries of the peripheral region.
276 11972 : input_ext_bd_pts.push_back(input_mesh->point(boundary_ordered_node_list[i]));
277 : tmp_azi =
278 11972 : atan2(input_ext_bd_pts.back()(1) - origin_pt(1), input_ext_bd_pts.back()(0) - origin_pt(0));
279 11972 : output_ext_bd_pts.push_back(Point(_peripheral_ring_radius * std::cos(tmp_azi),
280 11972 : _peripheral_ring_radius * std::sin(tmp_azi),
281 : origin_pt(2)));
282 : // a vector of tuples using azimuthal angle as the key to facilitate sorting
283 : azi_points.push_back(
284 11972 : std::make_tuple(tmp_azi, input_ext_bd_pts.back(), output_ext_bd_pts.back()));
285 : // a simple vector of azimuthal angles for radius correction purposes (in degree)
286 11972 : azi_array.push_back(tmp_azi / M_PI * 180.0);
287 : }
288 : // Only for quadratic sides, if the index of the first node after sorting is even, it is a vertex
289 : // node; otherwise, it is a midpoint node.
290 : const bool is_first_node_vertex =
291 104 : order == 1 ||
292 : !(std::distance(azi_array.begin(), std::min_element(azi_array.begin(), azi_array.end())) % 2);
293 99 : std::sort(azi_points.begin(), azi_points.end());
294 99 : std::sort(azi_array.begin(), azi_array.end());
295 :
296 : // Angles defined by three neighboring nodes on input mesh's external boundary
297 : std::vector<Real> input_bdry_angles;
298 : // Normal directions of input boundary nodes
299 : std::vector<Point> input_bdry_nd;
300 12071 : for (const auto i : index_range(azi_points))
301 : {
302 : Point p1, p2;
303 : const Point pn = Point(0.0, 0.0, 1.0);
304 11972 : if (i == 0)
305 : {
306 99 : p1 = std::get<1>(azi_points[i + 1]) - std::get<1>(azi_points[i]);
307 99 : p2 = std::get<1>(azi_points.back()) - std::get<1>(azi_points[i]);
308 : }
309 11873 : else if (i == azi_points.size() - 1)
310 : {
311 99 : p1 = std::get<1>(azi_points.front()) - std::get<1>(azi_points.back());
312 99 : p2 = std::get<1>(azi_points[i - 1]) - std::get<1>(azi_points.back());
313 : }
314 : else
315 : {
316 11774 : p1 = std::get<1>(azi_points[i + 1]) - std::get<1>(azi_points[i]);
317 11774 : p2 = std::get<1>(azi_points[i - 1]) - std::get<1>(azi_points[i]);
318 : }
319 : // Use cross point to get perpendicular direction
320 11972 : const Point p1n = (p1.cross(pn)).unit();
321 11972 : const Point p2n = -(p2.cross(pn)).unit();
322 11972 : Real tmp = p1 * p2 / p1.norm() / p2.norm();
323 : // Make sure acos() gets valid input
324 11972 : tmp = tmp > 1.0 ? 1.0 : (tmp < -1.0 ? -1.0 : tmp);
325 11972 : input_bdry_angles.push_back(acos(tmp) / 2.0);
326 23944 : input_bdry_nd.push_back((p1n + p2n).unit());
327 : }
328 :
329 : // 2D vector containing all the node positions of the peripheral region
330 99 : std::vector<std::vector<Point>> points_array(total_peripheral_layer_num + 1,
331 99 : std::vector<Point>(input_ext_node_num));
332 : // 2D vector containing all the node ids of the peripheral region
333 : std::vector<std::vector<dof_id_type>> node_id_array(total_peripheral_layer_num + 1,
334 99 : std::vector<dof_id_type>(input_ext_node_num));
335 : // Reference outmost layer of inner boundary layer
336 : std::vector<Point> ref_inner_bdry_surf;
337 : // First loop
338 12071 : for (const auto i : make_range(input_ext_node_num))
339 : {
340 : // Inner boundary nodes of the peripheral region
341 11972 : points_array[0][i] = std::get<1>(azi_points[i]);
342 : // Define outer layer of inner boundary layer
343 11972 : if (_peripheral_inner_boundary_layer_params.intervals)
344 : {
345 : // Outside point of the inner boundary layer
346 : const Point ref_inner_boundary_shift =
347 2504 : (_peripheral_inner_boundary_layer_params.width / sin(input_bdry_angles[i])) *
348 : input_bdry_nd[i];
349 2504 : ref_inner_bdry_surf.push_back(points_array[0][i] + ref_inner_boundary_shift);
350 : }
351 : }
352 :
353 99 : if (_peripheral_inner_boundary_layer_params.intervals)
354 12 : innerBdryLayerNodesDefiner(input_ext_node_num,
355 : input_bdry_angles,
356 : ref_inner_bdry_surf,
357 : inner_peripheral_bias_terms,
358 : azi_array,
359 : origin_pt,
360 : points_array);
361 99 : const Real correction_factor = _preserve_volumes
362 99 : ? PolygonalMeshGenerationUtils::radiusCorrectionFactor(
363 : azi_array, true, order, is_first_node_vertex)
364 : : 1.0;
365 : // Loop to handle outer boundary layer and main region
366 11483 : for (const auto i : make_range(input_ext_node_num))
367 : {
368 : // Outer boundary nodes of the peripheral region
369 11386 : points_array[total_peripheral_layer_num][i] = std::get<2>(azi_points[i]) * correction_factor;
370 : // Outer boundary layer points
371 11386 : if (_peripheral_outer_boundary_layer_params.intervals)
372 : {
373 : // Inner point of the outer boundary layer
374 : const Point outer_boundary_shift =
375 1918 : -Point(std::cos(std::get<0>(azi_points[i])), std::sin(std::get<0>(azi_points[i])), 0.0) *
376 : _peripheral_outer_boundary_layer_params.width;
377 1918 : for (const auto j :
378 10116 : make_range(uint(1), _peripheral_outer_boundary_layer_params.intervals + 1))
379 8198 : points_array[total_peripheral_layer_num - j][i] =
380 : points_array[total_peripheral_layer_num][i] +
381 8198 : outer_boundary_shift * outer_peripheral_bias_terms[j - 1];
382 : }
383 : // Use interpolation to get main region
384 11386 : if (MooseUtils::absoluteFuzzyGreaterEqual(
385 11386 : (points_array[_peripheral_inner_boundary_layer_params.intervals][i] - origin_pt).norm(),
386 11386 : (points_array[_peripheral_inner_boundary_layer_params.intervals +
387 22772 : _peripheral_layer_num * order][i] -
388 : origin_pt)
389 11386 : .norm()))
390 : {
391 2 : paramError("peripheral_inner_boundary_layer_width",
392 : "The summation of peripheral_inner_boundary_layer_width and "
393 : "peripheral_outer_boundary_layer_width must be smaller than the thickness of "
394 : "peripheral ring region.");
395 : }
396 :
397 28820 : for (const auto j : make_range(uint(1), _peripheral_layer_num * order))
398 17436 : points_array[j + _peripheral_inner_boundary_layer_params.intervals][i] =
399 : points_array[_peripheral_inner_boundary_layer_params.intervals][i] *
400 17436 : (1.0 - main_peripheral_bias_terms[j - 1]) +
401 : points_array[_peripheral_inner_boundary_layer_params.intervals +
402 : _peripheral_layer_num * order][i] *
403 17436 : main_peripheral_bias_terms[j - 1];
404 : }
405 97 : unsigned int num_total_nodes = (total_peripheral_layer_num + 1) * input_ext_node_num;
406 97 : std::vector<Node *> nodes(num_total_nodes); // reserve nodes pointers
407 : dof_id_type node_id = 0;
408 :
409 : // Add nodes to the new mesh
410 510 : for (const auto i : make_range(total_peripheral_layer_num + 1))
411 56793 : for (const auto j : make_range(input_ext_node_num))
412 : {
413 56380 : nodes[node_id] = mesh->add_point(points_array[i][j], node_id);
414 56380 : node_id_array[i][j] = node_id;
415 56380 : node_id++;
416 : }
417 : // Add elements to the new mesh
418 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
419 97 : const unsigned int index_shift = is_first_node_vertex ? 0 : 1;
420 363 : for (const auto i : make_range(total_peripheral_layer_num / order))
421 40498 : for (const auto j : make_range(input_ext_node_num / order))
422 : {
423 40232 : std::unique_ptr<Elem> new_elem;
424 40232 : new_elem = std::make_unique<Quad4>();
425 40232 : if (order == 2)
426 : {
427 1600 : new_elem = std::make_unique<Quad9>();
428 1600 : new_elem->set_node(4, nodes[node_id_array[i * order + 1][j * order + index_shift]]);
429 1600 : new_elem->set_node(5,
430 1600 : nodes[node_id_array[(i + 1) * order][(j * order + 1 + index_shift) %
431 1600 : input_ext_node_num]]);
432 1600 : new_elem->set_node(6,
433 1600 : nodes[node_id_array[i * order + 1][((j + 1) * order + index_shift) %
434 1600 : input_ext_node_num]]);
435 1600 : new_elem->set_node(
436 1600 : 7, nodes[node_id_array[i * order][(j * order + 1 + index_shift) % input_ext_node_num]]);
437 1600 : new_elem->set_node(8,
438 : nodes[node_id_array[i * order + 1][(j * order + 1 + index_shift) %
439 1600 : input_ext_node_num]]);
440 : }
441 40232 : new_elem->set_node(0, nodes[node_id_array[i * order][j * order + index_shift]]);
442 40232 : new_elem->set_node(1, nodes[node_id_array[(i + 1) * order][j * order + index_shift]]);
443 40232 : new_elem->set_node(2,
444 40232 : nodes[node_id_array[(i + 1) * order][((j + 1) * order + index_shift) %
445 40232 : input_ext_node_num]]);
446 40232 : new_elem->set_node(
447 40232 : 3, nodes[node_id_array[i * order][((j + 1) * order + index_shift) % input_ext_node_num]]);
448 40232 : new_elem->subdomain_id() = _peripheral_ring_block_id;
449 :
450 40232 : Elem * added_elem = mesh->add_elem(std::move(new_elem));
451 :
452 40232 : if (i == 0)
453 11188 : boundary_info.add_side(added_elem, 3, OUTER_SIDESET_ID_ALT);
454 40232 : if (i == total_peripheral_layer_num / order - 1)
455 11188 : boundary_info.add_side(added_elem, 1, OUTER_SIDESET_ID);
456 40232 : }
457 :
458 : // This would make sure that the boundary OUTER_SIDESET_ID is deleted after stitching.
459 97 : if (_input_mesh_external_bid != OUTER_SIDESET_ID)
460 36 : MooseMesh::changeBoundaryId(*input_mesh, _input_mesh_external_bid, OUTER_SIDESET_ID, false);
461 97 : mesh->prepare_for_use();
462 : // Use input_mesh here to retain the subdomain name map
463 97 : input_mesh->stitch_meshes(*mesh, OUTER_SIDESET_ID, OUTER_SIDESET_ID_ALT, TOLERANCE, true, false);
464 :
465 : // Assign subdomain name to the new block if applicable
466 194 : if (isParamValid("peripheral_ring_block_name"))
467 71 : input_mesh->subdomain_name(_peripheral_ring_block_id) = _peripheral_ring_block_name;
468 : // Assign customized external boundary id
469 97 : if (_external_boundary_id > 0)
470 0 : MooseMesh::changeBoundaryId(*input_mesh, OUTER_SIDESET_ID, _external_boundary_id, false);
471 : // Assign customized external boundary name
472 97 : if (!_external_boundary_name.empty())
473 : {
474 36 : input_mesh->get_boundary_info().sideset_name(
475 36 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
476 36 : _external_boundary_name;
477 36 : input_mesh->get_boundary_info().nodeset_name(
478 36 : _external_boundary_id > 0 ? _external_boundary_id : (boundary_id_type)OUTER_SIDESET_ID) =
479 : _external_boundary_name;
480 : }
481 :
482 97 : _input->set_isnt_prepared();
483 194 : return dynamic_pointer_cast<MeshBase>(_input);
484 97 : }
485 :
486 : void
487 12 : PeripheralRingMeshGenerator::innerBdryLayerNodesDefiner(
488 : const unsigned int input_ext_node_num,
489 : const std::vector<Real> input_bdry_angles,
490 : const std::vector<Point> ref_inner_bdry_surf,
491 : const std::vector<Real> inner_peripheral_bias_terms,
492 : const std::vector<Real> azi_array,
493 : const Point origin_pt,
494 : std::vector<std::vector<Point>> & points_array) const
495 : {
496 : // Check if any azimuthal angles are messed after inner boundary layer addition
497 12 : std::vector<bool> delete_mark(input_ext_node_num, false);
498 : std::vector<Real> interp_azi_data, interp_x_data, interp_y_data;
499 12 : std::unique_ptr<LinearInterpolation> linterp_x, linterp_y;
500 : // Mark the to-be-deleted elements
501 2516 : for (const auto i : make_range(input_ext_node_num))
502 : {
503 : // For a zig-zag external boundary, when we add a conformal layer onto it by translating the
504 : // nodes in the surface normal direction, it is possible that the azimuthal order is flipped.
505 : // As shown below, o's are the original boundary nodes, and *'s are the nodes after
506 : // translation by the boundary layer thickness.
507 : // * *
508 : // | | | /
509 : // o o-------/--*
510 : // | | / |
511 : // | outside --> | / |
512 : // | | / |
513 : // o--------o-- o-------o--
514 : // To mitigate this flipping issue, we check the node flipping here using the cross product of
515 : // neighboring node-to-origin vectors. Flipped nodes are marked and excluded during the
516 : // follow-up interpolation.
517 2504 : if (!MooseUtils::absoluteFuzzyEqual(input_bdry_angles[i], M_PI / 2.0, 1.0e-4))
518 : {
519 : unsigned int index_shift = 1;
520 : bool minus_shift_flag = true;
521 : bool plus_shift_flag = true;
522 : // Need to check both directions for multiple shifts
523 : // Half of the external boundary nodes are used as an upper limit to avert infinite loops
524 332 : while (index_shift < input_ext_node_num / 2 && (minus_shift_flag || plus_shift_flag))
525 : {
526 166 : if (((ref_inner_bdry_surf[MathUtils::euclideanMod(((int)i - (int)index_shift),
527 166 : (int)input_ext_node_num)] -
528 : origin_pt)
529 166 : .cross(ref_inner_bdry_surf[i] - origin_pt))(2) <= 0.0 &&
530 : minus_shift_flag)
531 : delete_mark[MathUtils::euclideanMod(((int)i - (int)index_shift),
532 : (int)input_ext_node_num)] = true;
533 : else
534 : minus_shift_flag = false;
535 166 : if (((ref_inner_bdry_surf[(i + index_shift) % input_ext_node_num] - origin_pt)
536 166 : .cross(ref_inner_bdry_surf[i] - origin_pt))(2) >= 0.0 &&
537 : plus_shift_flag)
538 : delete_mark[(i + index_shift) % input_ext_node_num] = true;
539 : else
540 : plus_shift_flag = false;
541 166 : index_shift++;
542 : }
543 : }
544 : }
545 : // Create vectors for interpolation
546 : // Due to the flip issue, linear interpolation is used here to mark the location of the boundary
547 : // layer's outer boundary.
548 2516 : for (const auto i : make_range(input_ext_node_num))
549 : {
550 2504 : if (!delete_mark[i])
551 : {
552 2504 : interp_azi_data.push_back(atan2(ref_inner_bdry_surf[i](1) - origin_pt(1),
553 2504 : ref_inner_bdry_surf[i](0) - origin_pt(0)));
554 2504 : interp_x_data.push_back(ref_inner_bdry_surf[i](0));
555 2504 : interp_y_data.push_back(ref_inner_bdry_surf[i](1));
556 2504 : if (interp_azi_data.size() > 1)
557 2492 : if (interp_azi_data.back() < interp_azi_data[interp_azi_data.size() - 2])
558 0 : interp_azi_data.back() += 2.0 * M_PI;
559 : }
560 : }
561 12 : const Real interp_x0 = interp_x_data.front();
562 12 : const Real interp_xt = interp_x_data.back();
563 12 : const Real interp_y0 = interp_y_data.front();
564 12 : const Real interp_yt = interp_y_data.back();
565 12 : const Real incept_azi0 = interp_azi_data.front();
566 12 : const Real incept_azit = interp_azi_data.back();
567 :
568 12 : if (MooseUtils::absoluteFuzzyGreaterThan(incept_azi0, -M_PI))
569 : {
570 5 : interp_azi_data.insert(interp_azi_data.begin(), incept_azit - 2.0 * M_PI);
571 5 : interp_x_data.insert(interp_x_data.begin(), interp_xt);
572 5 : interp_y_data.insert(interp_y_data.begin(), interp_yt);
573 : }
574 : else
575 7 : interp_azi_data.front() = -M_PI;
576 12 : if (MooseUtils::absoluteFuzzyLessThan(incept_azit, M_PI))
577 : {
578 7 : interp_azi_data.push_back(incept_azi0 + 2 * M_PI);
579 7 : interp_x_data.push_back(interp_x0);
580 7 : interp_y_data.push_back(interp_y0);
581 : }
582 : else
583 5 : interp_azi_data.back() = M_PI;
584 :
585 : // Establish interpolation
586 24 : linterp_x = std::make_unique<LinearInterpolation>(interp_azi_data, interp_x_data);
587 24 : linterp_y = std::make_unique<LinearInterpolation>(interp_azi_data, interp_y_data);
588 : // Loop to handle inner boundary layer
589 2516 : for (const auto i : make_range(input_ext_node_num))
590 : {
591 : // Outside point of the inner boundary layer
592 : // Using interpolation, the azimuthal angles do not need to be changed.
593 2504 : const Point inner_boundary_shift = Point(linterp_x->sample(azi_array[i] / 180.0 * M_PI),
594 2504 : linterp_y->sample(azi_array[i] / 180.0 * M_PI),
595 : origin_pt(2)) -
596 2504 : points_array[0][i];
597 11288 : for (const auto j : make_range(uint(1), _peripheral_inner_boundary_layer_params.intervals + 1))
598 8784 : points_array[j][i] =
599 8784 : points_array[0][i] + inner_boundary_shift * inner_peripheral_bias_terms[j - 1];
600 : }
601 24 : }
|
