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