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