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 "AzimuthalBlockSplitGenerator.h"
11 : #include "MooseMeshUtils.h"
12 : #include "PolygonalMeshGenerationUtils.h"
13 :
14 : // C++ includes
15 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
16 :
17 : registerMooseObject("ReactorApp", AzimuthalBlockSplitGenerator);
18 :
19 : InputParameters
20 124 : AzimuthalBlockSplitGenerator::validParams()
21 : {
22 124 : InputParameters params = PolygonMeshGeneratorBase::validParams();
23 248 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh to be modified.");
24 248 : params.addRequiredParam<std::vector<SubdomainName>>(
25 : "old_blocks", "The list of blocks in the input mesh that need to be modified.");
26 248 : params.addRequiredParam<std::vector<subdomain_id_type>>(
27 : "new_block_ids", "The block IDs to be used for the new selected azimuthal angle blocks.");
28 248 : params.addParam<std::vector<SubdomainName>>(
29 : "new_block_names",
30 : {},
31 : "The optional block names to be used for the new selected azimulathal angle blocks.");
32 248 : params.addParam<bool>("preserve_volumes",
33 248 : true,
34 : "Volume of concentric circles can be preserved using this function.");
35 248 : params.addRequiredRangeCheckedParam<Real>("start_angle",
36 : "start_angle>=0.0 & start_angle<360.0",
37 : "Starting azimuthal angle of the new block.");
38 248 : params.addRequiredRangeCheckedParam<Real>("angle_range",
39 : "angle_range>0.0 & angle_range<=180.0",
40 : "Azimuthal angle range of the new block.");
41 124 : params.addClassDescription(
42 : "This AzimuthalBlockSplitGenerator object takes in a polygon/hexagon concentric circle mesh "
43 : "and renames blocks on a user-defined azimuthal segment / wedge of the mesh.");
44 :
45 124 : return params;
46 0 : }
47 :
48 63 : AzimuthalBlockSplitGenerator::AzimuthalBlockSplitGenerator(const InputParameters & parameters)
49 : : PolygonMeshGeneratorBase(parameters),
50 63 : _input_name(getParam<MeshGeneratorName>("input")),
51 126 : _new_block_ids(getParam<std::vector<subdomain_id_type>>("new_block_ids")),
52 126 : _new_block_names(getParam<std::vector<SubdomainName>>("new_block_names")),
53 126 : _preserve_volumes(getParam<bool>("preserve_volumes")),
54 126 : _start_angle(getParam<Real>("start_angle") + 90.0),
55 126 : _angle_range(getParam<Real>("angle_range")),
56 63 : _azimuthal_angle_meta(
57 126 : declareMeshProperty<std::vector<Real>>("azimuthal_angle_meta", std::vector<Real>())),
58 126 : _input(getMeshByName(_input_name))
59 : {
60 63 : declareMeshProperty<Real>("pattern_pitch_meta", 0.0);
61 126 : declareMeshProperty<bool>("is_control_drum_meta", true);
62 63 : if (!_new_block_names.empty() && _new_block_names.size() != _new_block_ids.size())
63 2 : paramError("new_block_names",
64 : "This parameter, if provided, must have the same size as new_block_ids.");
65 : _num_sectors_per_side_meta =
66 61 : getMeshProperty<std::vector<unsigned int>>("num_sectors_per_side_meta", _input_name);
67 61 : _start_angle = _start_angle >= 360.0 ? _start_angle - 360.0 : _start_angle;
68 61 : _end_angle = (_start_angle + _angle_range >= 360.0) ? (_start_angle + _angle_range - 360.0)
69 : : (_start_angle + _angle_range);
70 61 : }
71 :
72 : std::unique_ptr<MeshBase>
73 61 : AzimuthalBlockSplitGenerator::generate()
74 : {
75 61 : auto replicated_mesh_ptr = dynamic_cast<ReplicatedMesh *>(_input.get());
76 61 : if (!replicated_mesh_ptr)
77 0 : paramError("input", "Input is not a replicated mesh, which is required");
78 :
79 : ReplicatedMesh & mesh = *replicated_mesh_ptr;
80 :
81 : // Check the order of the input mesh's elements
82 : // Meanwhile, record the vertex average of each element for future comparison
83 : unsigned short order = 0;
84 : std::map<dof_id_type, Point> vertex_avgs;
85 6468 : for (const auto & elem : mesh.element_ptr_range())
86 : {
87 6350 : switch (elem->type())
88 : {
89 5338 : case TRI3:
90 : case QUAD4:
91 5338 : if (order == 2)
92 0 : paramError("input", "This mesh contains elements of different orders.");
93 : order = 1;
94 : break;
95 1010 : case TRI6:
96 : case TRI7:
97 : case QUAD8:
98 : case QUAD9:
99 1010 : if (order == 1)
100 2 : paramError("input", "This mesh contains elements of different orders.");
101 : order = 2;
102 : break;
103 2 : default:
104 2 : paramError("input", "This mesh contains elements of unsupported type.");
105 : }
106 6346 : vertex_avgs[elem->id()] = elem->vertex_average();
107 57 : }
108 :
109 : _old_block_ids =
110 171 : MooseMeshUtils::getSubdomainIDs(mesh, getParam<std::vector<SubdomainName>>("old_blocks"));
111 57 : if (_old_block_ids.size() != _new_block_ids.size())
112 4 : paramError("new_block_ids", "This parameter must have the same size as old_blocks.");
113 :
114 : // Check that the block ids/names exist in the mesh
115 : std::set<SubdomainID> mesh_blocks;
116 53 : mesh.subdomain_ids(mesh_blocks);
117 :
118 190 : for (unsigned int i = 0; i < _old_block_ids.size(); i++)
119 139 : if (_old_block_ids[i] == Moose::INVALID_BLOCK_ID || !mesh_blocks.count(_old_block_ids[i]))
120 2 : paramError("old_blocks",
121 : "This parameter contains blocks that do not exist in the input mesh. Error "
122 2 : "triggered by block: " +
123 2 : getParam<std::vector<SubdomainName>>("old_blocks")[i]);
124 :
125 51 : if (std::find(_old_block_ids.begin(),
126 : _old_block_ids.end(),
127 : getMeshProperty<subdomain_id_type>("quad_center_block_id", _input_name)) !=
128 : _old_block_ids.end())
129 2 : paramError(
130 : "old_blocks",
131 : "This parameter contains a block that involves center quad elements, azimuthal splitting "
132 : "is currently not supported in this case.");
133 :
134 49 : MeshTools::Modification::rotate(mesh, 90.0, 0.0, 0.0);
135 49 : _azimuthal_angle_meta = azimuthalAnglesCollector(mesh, -180.0, 180.0, ANGLE_DEGREE);
136 45 : std::vector<Real> azimuthal_angles_vtx;
137 : const bool is_first_value_vertex =
138 49 : order == 1 ? true : MooseUtils::absoluteFuzzyEqual(_azimuthal_angle_meta[0], -180.0);
139 : if (order == 1)
140 40 : azimuthal_angles_vtx = _azimuthal_angle_meta;
141 : else
142 513 : for (const auto & i : make_range(_azimuthal_angle_meta.size()))
143 504 : if (i % 2 != is_first_value_vertex)
144 252 : azimuthal_angles_vtx.push_back(_azimuthal_angle_meta[i]);
145 :
146 : // So that this class works for both derived classes of PolygonMeshGeneratorBase
147 98 : auto pattern_pitch_meta = std::max(getMeshProperty<Real>("pitch_meta", _input_name),
148 49 : getMeshProperty<Real>("pattern_pitch_meta", _input_name));
149 49 : setMeshProperty("pattern_pitch_meta", pattern_pitch_meta);
150 :
151 : Real radiusCorrectionFactor_original =
152 49 : _preserve_volumes ? PolygonalMeshGenerationUtils::radiusCorrectionFactor(
153 49 : _azimuthal_angle_meta, true, order, is_first_value_vertex)
154 : : 1.0;
155 :
156 49 : const Real azi_tol = 1.0E-6;
157 49 : const Real rad_tol = 1.0e-6;
158 49 : const Real side_angular_range = 360.0 / (Real)_num_sectors_per_side_meta.size();
159 : const Real side_angular_shift =
160 49 : _num_sectors_per_side_meta.size() % 2 == 0 ? 0.0 : side_angular_range / 2.0;
161 :
162 49 : _start_angle = _start_angle > 180.0 ? _start_angle - 360.0 : _start_angle;
163 49 : _end_angle = _end_angle > 180.0 ? _end_angle - 360.0 : _end_angle;
164 :
165 : Real original_start_down;
166 49 : std::vector<Real>::iterator original_start_down_it;
167 : Real original_start_up;
168 49 : std::vector<Real>::iterator original_start_up_it;
169 :
170 : // Identify the mesh azimuthal angles of the elements that need to be modified for the starting
171 : // edge of the absorber region.
172 49 : angleIdentifier(_start_angle,
173 : original_start_down,
174 : original_start_down_it,
175 : original_start_up,
176 : original_start_up_it,
177 : azimuthal_angles_vtx);
178 :
179 : Real original_end_down;
180 49 : std::vector<Real>::iterator original_end_down_it;
181 : Real original_end_up;
182 49 : std::vector<Real>::iterator original_end_up_it;
183 :
184 : // Identify the mesh azimuthal angles of the elements that need to be modified for the ending
185 : // edge of the absorber region.
186 49 : angleIdentifier(_end_angle,
187 : original_end_down,
188 : original_end_down_it,
189 : original_end_up,
190 : original_end_up_it,
191 : azimuthal_angles_vtx);
192 :
193 : Real azi_to_mod_start;
194 : Real azi_to_keep_start;
195 :
196 : // For the two mesh azimuthal angles identified, determine which need to be moved to the starting
197 : // edge position.
198 49 : angleModifier(side_angular_shift,
199 : side_angular_range,
200 : azi_tol,
201 : _start_angle,
202 : original_start_down,
203 : original_start_down_it,
204 : original_start_up,
205 : original_start_up_it,
206 : azi_to_keep_start,
207 : azi_to_mod_start);
208 :
209 : Real azi_to_mod_end;
210 : Real azi_to_keep_end;
211 :
212 : // For the two mesh azimuthal angles identified, determine which need to be moved to the ending
213 : // edge position.
214 49 : angleModifier(side_angular_shift,
215 : side_angular_range,
216 : azi_tol,
217 : _end_angle,
218 : original_end_down,
219 : original_end_down_it,
220 : original_end_up,
221 : original_end_up_it,
222 : azi_to_keep_end,
223 : azi_to_mod_end);
224 :
225 49 : if (azi_to_mod_end == azi_to_mod_start)
226 2 : paramError("angle_range",
227 : "The azimuthal intervals of the input mesh are too coarse for this parameter.");
228 :
229 92 : const auto & side_list = mesh.get_boundary_info().build_side_list();
230 :
231 : // See if the block that contains the external boundary is modified or not
232 47 : bool external_block_change(false);
233 3887 : for (unsigned int i = 0; i < side_list.size(); i++)
234 : {
235 3840 : if (std::get<2>(side_list[i]) == OUTER_SIDESET_ID)
236 : {
237 1280 : dof_id_type block_id_tmp = mesh.elem_ref(std::get<0>(side_list[i])).subdomain_id();
238 1280 : if (std::find(_old_block_ids.begin(), _old_block_ids.end(), block_id_tmp) !=
239 : _old_block_ids.end())
240 252 : external_block_change = true;
241 : }
242 : }
243 :
244 47 : mesh.get_boundary_info().build_node_list_from_side_list();
245 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
246 92 : mesh.get_boundary_info().build_node_list();
247 :
248 45 : std::vector<std::pair<Real, dof_id_type>> node_id_keep_start;
249 45 : std::vector<std::pair<Real, dof_id_type>> node_id_mod_start;
250 45 : std::vector<std::pair<Real, dof_id_type>> node_id_keep_end;
251 45 : std::vector<std::pair<Real, dof_id_type>> node_id_mod_end;
252 :
253 : // Determine which nodes are involved in the elements that are intercepted by the starting/ending
254 : // angles
255 : Real max_quad_elem_rad(0.0);
256 47 : if (mesh.elem_ref(0).n_vertices() == 4)
257 81 : for (dof_id_type i = 0;
258 90 : i < (_num_sectors_per_side_meta[0] / 2 + 1) * (_num_sectors_per_side_meta[0] / 2 + 1);
259 : i++)
260 : max_quad_elem_rad =
261 81 : max_quad_elem_rad < mesh.node_ref(i).norm() ? mesh.node_ref(i).norm() : max_quad_elem_rad;
262 16667 : for (const auto & node_ptr : as_range(mesh.nodes_begin(), mesh.nodes_end()))
263 : {
264 8263 : const Node & node = *node_ptr;
265 8263 : Real node_azi = atan2(node(1), node(0)) * 180.0 / M_PI;
266 8263 : Real node_rad = std::sqrt(node(0) * node(0) + node(1) * node(1));
267 8263 : if (node_rad > max_quad_elem_rad + rad_tol &&
268 7892 : (std::abs(node_azi - azi_to_mod_start) < azi_tol ||
269 7668 : std::abs(std::abs(node_azi - azi_to_mod_start) - 360.0) < azi_tol))
270 224 : node_id_mod_start.push_back(std::make_pair(node_rad, node.id()));
271 8263 : if (node_rad > max_quad_elem_rad + rad_tol &&
272 7892 : (std::abs(node_azi - azi_to_keep_start) < azi_tol ||
273 7668 : std::abs(std::abs(node_azi - azi_to_keep_start) - 360.0) < azi_tol))
274 224 : node_id_keep_start.push_back(std::make_pair(node_rad, node.id()));
275 8263 : if (node_rad > max_quad_elem_rad + rad_tol &&
276 7892 : (std::abs(node_azi - azi_to_mod_end) < azi_tol ||
277 7668 : std::abs(std::abs(node_azi - azi_to_mod_end) - 360.0) < azi_tol))
278 224 : node_id_mod_end.push_back(std::make_pair(node_rad, node.id()));
279 8263 : if (node_rad > max_quad_elem_rad + rad_tol &&
280 7892 : (std::abs(node_azi - azi_to_keep_end) < azi_tol ||
281 7668 : std::abs(std::abs(node_azi - azi_to_keep_end) - 360.0) < azi_tol))
282 224 : node_id_keep_end.push_back(std::make_pair(node_rad, node.id()));
283 47 : }
284 : // Sort the involved nodes using radius as the key; this facilitates the determination of circular
285 : // regions nodes
286 47 : std::sort(node_id_mod_start.begin(), node_id_mod_start.end());
287 47 : std::sort(node_id_keep_start.begin(), node_id_keep_start.end());
288 47 : std::sort(node_id_mod_end.begin(), node_id_mod_end.end());
289 47 : std::sort(node_id_keep_end.begin(), node_id_keep_end.end());
290 45 : std::vector<Real> circular_rad_list;
291 45 : std::vector<Real> non_circular_rad_list;
292 :
293 : // Modify the nodes with the azimuthal angles identified before.
294 47 : nodeModifier(mesh,
295 : node_id_mod_start,
296 : node_id_keep_start,
297 : circular_rad_list,
298 : non_circular_rad_list,
299 : node_list,
300 : _start_angle,
301 : external_block_change,
302 : rad_tol);
303 47 : nodeModifier(mesh,
304 : node_id_mod_end,
305 : node_id_keep_end,
306 : circular_rad_list,
307 : non_circular_rad_list,
308 : node_list,
309 : _end_angle,
310 : external_block_change,
311 : rad_tol);
312 :
313 : const Real max_circular_radius =
314 47 : *std::max_element(circular_rad_list.begin(), circular_rad_list.end());
315 : const Real min_non_circular_radius =
316 47 : *std::min_element(non_circular_rad_list.begin(), non_circular_rad_list.end());
317 :
318 : // Before re-correcting the radii, correct the mid-point of the elements that are altered
319 47 : if (order == 2)
320 2034 : for (const auto & elem : mesh.element_ptr_range())
321 : {
322 1008 : const Point & original_vertex_avg = vertex_avgs[elem->id()];
323 1008 : const Point & new_vertex_avg = elem->vertex_average();
324 1008 : if (MooseUtils::absoluteFuzzyGreaterThan((original_vertex_avg - new_vertex_avg).norm(), 0.0))
325 : {
326 1008 : if (elem->type() == TRI6 || elem->type() == TRI7)
327 : {
328 252 : *elem->node_ptr(3) = midPointCorrector(
329 252 : *elem->node_ptr(0), *elem->node_ptr(1), max_circular_radius, rad_tol);
330 252 : *elem->node_ptr(4) = midPointCorrector(
331 252 : *elem->node_ptr(1), *elem->node_ptr(2), max_circular_radius, rad_tol);
332 252 : *elem->node_ptr(5) = midPointCorrector(
333 252 : *elem->node_ptr(2), *elem->node_ptr(0), max_circular_radius, rad_tol);
334 252 : if (elem->type() == TRI7)
335 : *elem->node_ptr(6) = (*elem->node_ptr(0) + *elem->node_ptr(1) + *elem->node_ptr(2) +
336 0 : *elem->node_ptr(3) + *elem->node_ptr(4) + *elem->node_ptr(5)) /
337 : 6.0;
338 : }
339 756 : else if (elem->type() == QUAD8 || elem->type() == QUAD9)
340 : {
341 756 : *elem->node_ptr(4) = midPointCorrector(
342 756 : *elem->node_ptr(0), *elem->node_ptr(1), max_circular_radius, rad_tol);
343 756 : *elem->node_ptr(5) = midPointCorrector(
344 756 : *elem->node_ptr(1), *elem->node_ptr(2), max_circular_radius, rad_tol);
345 756 : *elem->node_ptr(6) = midPointCorrector(
346 756 : *elem->node_ptr(2), *elem->node_ptr(3), max_circular_radius, rad_tol);
347 756 : *elem->node_ptr(7) = midPointCorrector(
348 756 : *elem->node_ptr(3), *elem->node_ptr(0), max_circular_radius, rad_tol);
349 756 : if (elem->type() == QUAD9)
350 : *elem->node_ptr(8) = (*elem->node_ptr(0) + *elem->node_ptr(1) + *elem->node_ptr(2) +
351 : *elem->node_ptr(3) + *elem->node_ptr(4) + *elem->node_ptr(5) +
352 756 : *elem->node_ptr(6) + *elem->node_ptr(7)) /
353 : 8.0;
354 : }
355 : }
356 9 : }
357 :
358 45 : std::vector<Real> new_azimuthal_angle;
359 16667 : for (const auto & node_ptr : as_range(mesh.nodes_begin(), mesh.nodes_end()))
360 : {
361 8263 : const Node & node = *node_ptr;
362 8263 : if (MooseUtils::absoluteFuzzyEqual(
363 8263 : std::sqrt(node(0) * node(0) + node(1) * node(1)), max_circular_radius, rad_tol))
364 : {
365 1532 : Real node_azi = atan2(node(1), node(0)) * 180.0 / M_PI;
366 1532 : new_azimuthal_angle.push_back(node_azi);
367 : }
368 47 : }
369 47 : std::sort(new_azimuthal_angle.begin(), new_azimuthal_angle.end());
370 47 : _azimuthal_angle_meta = new_azimuthal_angle;
371 : const bool is_first_value_vertex_mod =
372 47 : order == 1 ? true : MooseUtils::absoluteFuzzyEqual(_azimuthal_angle_meta[0], -180.0);
373 :
374 : Real radiusCorrectionFactor_mod =
375 47 : _preserve_volumes ? PolygonalMeshGenerationUtils::radiusCorrectionFactor(
376 47 : _azimuthal_angle_meta, true, order, is_first_value_vertex_mod)
377 : : 1.0;
378 :
379 : // Re-correct Radii
380 47 : if (_preserve_volumes)
381 : {
382 47 : if (min_non_circular_radius <
383 47 : max_circular_radius / radiusCorrectionFactor_original * radiusCorrectionFactor_mod)
384 2 : mooseError("In AzimuthalBlockSplitGenerator ",
385 2 : _name,
386 : ": the circular region is overlapped with background region after correction.");
387 16209 : for (Node * node_ptr : as_range(mesh.nodes_begin(), mesh.nodes_end()))
388 : {
389 : Node & node = *node_ptr;
390 8037 : Real node_rad = std::sqrt(node(0) * node(0) + node(1) * node(1));
391 : // Any nodes with radii smaller than the threshold belong to circular regions
392 8037 : if (node_rad > rad_tol && node_rad <= max_circular_radius + rad_tol)
393 : {
394 6012 : const Real node_azi = atan2(node(1), node(0));
395 6012 : const Real node_rad_corr =
396 6012 : node_rad / radiusCorrectionFactor_original * radiusCorrectionFactor_mod;
397 6012 : node(0) = node_rad_corr * std::cos(node_azi);
398 6012 : node(1) = node_rad_corr * std::sin(node_azi);
399 : }
400 45 : }
401 : }
402 : // Assign New Block IDs
403 162 : for (unsigned int block_id_index = 0; block_id_index < _old_block_ids.size(); block_id_index++)
404 : for (auto & elem :
405 234 : as_range(mesh.active_subdomain_elements_begin(_old_block_ids[block_id_index]),
406 6948 : mesh.active_subdomain_elements_end(_old_block_ids[block_id_index])))
407 : {
408 3240 : auto p_cent = elem->true_centroid();
409 3240 : auto p_cent_azi = atan2(p_cent(1), p_cent(0)) * 180.0 / M_PI;
410 3240 : if (_start_angle < _end_angle && p_cent_azi >= _start_angle && p_cent_azi <= _end_angle)
411 702 : elem->subdomain_id() = _new_block_ids[block_id_index];
412 2538 : else if (_start_angle > _end_angle &&
413 0 : (p_cent_azi >= _start_angle || p_cent_azi <= _end_angle))
414 0 : elem->subdomain_id() = _new_block_ids[block_id_index];
415 117 : }
416 : // Assign new Block Names if applicable
417 162 : for (unsigned int i = 0; i < _new_block_names.size(); i++)
418 117 : mesh.subdomain_name(_new_block_ids[i]) = _new_block_names[i];
419 45 : MeshTools::Modification::rotate(mesh, -90.0, 0.0, 0.0);
420 1521 : for (unsigned int i = 0; i < _azimuthal_angle_meta.size(); i++)
421 1476 : _azimuthal_angle_meta[i] = (_azimuthal_angle_meta[i] - 90.0 <= -180.0)
422 1476 : ? (_azimuthal_angle_meta[i] + 270.0)
423 : : _azimuthal_angle_meta[i] - 90.0;
424 45 : std::sort(_azimuthal_angle_meta.begin(), _azimuthal_angle_meta.end());
425 :
426 90 : return std::move(_input);
427 : }
428 :
429 : void
430 98 : AzimuthalBlockSplitGenerator::angleIdentifier(const Real & terminal_angle,
431 : Real & original_down,
432 : std::vector<Real>::iterator & original_down_it,
433 : Real & original_up,
434 : std::vector<Real>::iterator & original_up_it,
435 : std::vector<Real> & azimuthal_angles_vtx)
436 : {
437 :
438 : auto term_up =
439 98 : std::lower_bound(azimuthal_angles_vtx.begin(), azimuthal_angles_vtx.end(), terminal_angle);
440 98 : if (term_up == azimuthal_angles_vtx.begin())
441 : {
442 0 : original_up = *term_up;
443 0 : original_up_it = term_up;
444 0 : original_down = -180.0;
445 0 : original_down_it = azimuthal_angles_vtx.end() - 1;
446 : }
447 98 : else if (term_up == azimuthal_angles_vtx.end())
448 : {
449 0 : original_down = azimuthal_angles_vtx.back();
450 0 : original_down_it = azimuthal_angles_vtx.end() - 1;
451 0 : original_up = 180.0;
452 0 : original_up_it = azimuthal_angles_vtx.begin();
453 : }
454 : else
455 : {
456 98 : original_down = *(term_up - 1);
457 98 : original_down_it = term_up - 1;
458 98 : original_up = *term_up;
459 98 : original_up_it = term_up;
460 : }
461 98 : }
462 :
463 : void
464 98 : AzimuthalBlockSplitGenerator::angleModifier(const Real & side_angular_shift,
465 : const Real & side_angular_range,
466 : const Real & azi_tol,
467 : const Real & terminal_angle,
468 : const Real & original_down,
469 : std::vector<Real>::iterator & original_down_it,
470 : const Real & original_up,
471 : std::vector<Real>::iterator & original_up_it,
472 : Real & azi_to_keep,
473 : Real & azi_to_mod)
474 : {
475 : // Half interval is used to help determine which of the two identified azimuthal angles needs to
476 : // be moved.
477 98 : const Real half_interval = (original_up - original_down) / 2.0;
478 : // In this case, the lower azimuthal angle matches a vertex position, while the target angle is
479 : // not overlapped with the same vertex position. Thus the lower azimuthal angle is not moved
480 : // anyway.
481 98 : if (std::abs((original_down + side_angular_shift) / side_angular_range -
482 98 : std::round((original_down + side_angular_shift) / side_angular_range)) < azi_tol &&
483 45 : std::abs(original_down - terminal_angle) > azi_tol)
484 : {
485 45 : azi_to_keep = original_down;
486 45 : azi_to_mod = original_up;
487 45 : *original_up_it = terminal_angle;
488 : }
489 : // In this case, the upper azimuthal angle matches a vertex position, while the target angle is
490 : // not overlapped with the same vertex position. Thus the upper azimuthal angle is not moved
491 : // anyway.
492 53 : else if (std::abs((original_up + side_angular_shift) / side_angular_range -
493 53 : std::round((original_up + side_angular_shift) / side_angular_range)) <
494 53 : azi_tol &&
495 51 : std::abs(original_up - terminal_angle) > azi_tol)
496 : {
497 51 : azi_to_keep = original_up;
498 51 : azi_to_mod = original_down;
499 51 : *original_down_it = terminal_angle;
500 : }
501 : // Move upper azimuthal angle as it is closer to the target angle.
502 2 : else if (terminal_angle - original_down > half_interval)
503 : {
504 2 : azi_to_keep = original_down;
505 2 : azi_to_mod = original_up;
506 2 : *original_up_it = terminal_angle;
507 : }
508 : // Move lower azimuthal angle as it is closer to the target angle.
509 : else
510 : {
511 0 : azi_to_keep = original_up;
512 0 : azi_to_mod = original_down;
513 0 : *original_down_it = terminal_angle;
514 : }
515 98 : }
516 :
517 : void
518 94 : AzimuthalBlockSplitGenerator::nodeModifier(
519 : ReplicatedMesh & mesh,
520 : const std::vector<std::pair<Real, dof_id_type>> & node_id_mod,
521 : const std::vector<std::pair<Real, dof_id_type>> & node_id_keep,
522 : std::vector<Real> & circular_rad_list,
523 : std::vector<Real> & non_circular_rad_list,
524 : const std::vector<std::tuple<dof_id_type, boundary_id_type>> & node_list,
525 : const Real & term_angle,
526 : const bool & external_block_change,
527 : const Real & rad_tol)
528 : {
529 542 : for (unsigned int i = 0; i < node_id_mod.size(); i++)
530 : {
531 : // Circular regions, radius is not altered
532 448 : if (std::abs(node_id_mod[i].first - node_id_keep[i].first) < rad_tol)
533 : {
534 336 : Node & node_mod = mesh.node_ref(node_id_mod[i].second);
535 336 : circular_rad_list.push_back(node_id_mod[i].first);
536 336 : node_mod(0) = node_id_mod[i].first * std::cos(term_angle * M_PI / 180.0);
537 336 : node_mod(1) = node_id_mod[i].first * std::sin(term_angle * M_PI / 180.0);
538 : }
539 : else
540 : {
541 112 : non_circular_rad_list.push_back(node_id_mod[i].first);
542 : // For external boundary nodes, if the external block is not modified, the boundary nodes are
543 : // not moved.
544 : // Non-circular range, use intercept instead
545 112 : if (std::find(node_list.begin(),
546 : node_list.end(),
547 206 : std::make_tuple(node_id_mod[i].second, OUTER_SIDESET_ID)) == node_list.end() ||
548 94 : external_block_change)
549 : {
550 36 : Node & node_mod = mesh.node_ref(node_id_mod[i].second);
551 36 : const Node & node_keep = mesh.node_ref(node_id_keep[i].second);
552 36 : std::pair<Real, Real> pair_tmp = fourPointIntercept(
553 36 : std::make_pair(node_mod(0), node_mod(1)),
554 36 : std::make_pair(node_keep(0), node_keep(1)),
555 72 : std::make_pair(0.0, 0.0),
556 72 : std::make_pair(2.0 * node_id_mod[i].first * std::cos(term_angle * M_PI / 180.0),
557 36 : 2.0 * node_id_mod[i].first * std::sin(term_angle * M_PI / 180.0)));
558 36 : node_mod(0) = pair_tmp.first;
559 36 : node_mod(1) = pair_tmp.second;
560 : }
561 : }
562 : }
563 94 : }
564 :
565 : Point
566 3780 : AzimuthalBlockSplitGenerator::midPointCorrector(const Point vertex_0,
567 : const Point vertex_1,
568 : const Real max_radius,
569 : const Real rad_tol)
570 : {
571 : // Check if two vertices have the same radius
572 3780 : const Real r_vertex_0 = std::sqrt(vertex_0(0) * vertex_0(0) + vertex_0(1) * vertex_0(1));
573 3780 : const Real r_vertex_1 = std::sqrt(vertex_1(0) * vertex_1(0) + vertex_1(1) * vertex_1(1));
574 : // If both vertices have the same radius and they are located in the circular region,
575 : // their midpoint should have the same radius and has the average azimuthal angle of the two
576 : // vertices.
577 3780 : if (std::abs(r_vertex_0 - r_vertex_1) < rad_tol && r_vertex_0 < max_radius + rad_tol)
578 1512 : return (vertex_0 + vertex_1).unit() * r_vertex_0;
579 : else
580 : return (vertex_0 + vertex_1) / 2.0;
581 : }
|
