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 "HexagonalLatticeUtils.h"
11 : #include "MooseUtils.h"
12 : #include "GeometryUtils.h"
13 : #include "RotationMatrix.h"
14 :
15 : const Real HexagonalLatticeUtils::COS60 = 0.5;
16 : const Real HexagonalLatticeUtils::SIN60 = std::sqrt(3.0) / 2.0;
17 : const unsigned int HexagonalLatticeUtils::NUM_SIDES = 6;
18 :
19 201 : HexagonalLatticeUtils::HexagonalLatticeUtils(const Real bundle_inner_flat_to_flat,
20 : const Real pin_pitch,
21 : const Real pin_diameter,
22 : const Real wire_diameter,
23 : const Real wire_pitch,
24 : const unsigned int n_rings,
25 : const unsigned int axis,
26 201 : const Real rotation_around_axis)
27 201 : : _bundle_pitch(bundle_inner_flat_to_flat),
28 201 : _pin_pitch(pin_pitch),
29 201 : _pin_diameter(pin_diameter),
30 201 : _wire_diameter(wire_diameter),
31 201 : _wire_pitch(wire_pitch),
32 201 : _n_rings(n_rings),
33 201 : _axis(axis),
34 201 : _rotation_around_axis(rotation_around_axis),
35 201 : _rotation_matrix(
36 201 : RotationMatrix::rotVec2DToX(Point(std::cos(_rotation_around_axis * libMesh::pi / 180),
37 201 : -std::sin(_rotation_around_axis * libMesh::pi / 180),
38 : 0.))),
39 201 : _bundle_side_length(hexagonSide(_bundle_pitch)),
40 201 : _pin_area(M_PI * _pin_diameter * _pin_diameter / 4.0),
41 201 : _pin_circumference(M_PI * _pin_diameter),
42 201 : _wire_area(M_PI * _wire_diameter * _wire_diameter / 4.0),
43 201 : _wire_circumference(M_PI * _wire_diameter),
44 201 : _pin_surface_area_per_pitch(M_PI * _pin_diameter * _wire_pitch),
45 201 : _pin_volume_per_pitch(_pin_area * _wire_pitch),
46 201 : _P_over_D(_pin_pitch / _pin_diameter),
47 201 : _L_over_D(_wire_pitch / _pin_diameter)
48 : {
49 201 : auto idx = geom_utils::projectedIndices(_axis);
50 201 : _ix = idx.first;
51 201 : _iy = idx.second;
52 :
53 : // object is not tested and probably won't work if axis != 2
54 201 : if (_axis != 2)
55 0 : mooseError("The HexagonalLatticeUtils is currently limited to 'axis = 2'!");
56 :
57 201 : if (_pin_pitch < _pin_diameter)
58 0 : mooseError("Pin pitch of " + std::to_string(_pin_pitch) + " cannot fit pins of diameter " +
59 0 : std::to_string(_pin_diameter) + "!");
60 :
61 201 : if (_pin_pitch - _pin_diameter < _wire_diameter)
62 0 : mooseError("Wire diameter of " + std::to_string(_wire_diameter) +
63 0 : " cannot fit between pin-pin space of " +
64 0 : std::to_string(_pin_pitch - _pin_diameter) + "!");
65 :
66 201 : _unit_translation_x = {0.0, -SIN60, -SIN60, 0.0, SIN60, SIN60};
67 201 : _unit_translation_y = {1.0, COS60, -COS60, -1.0, -COS60, COS60};
68 :
69 : // Use rotation matrix on unit translation vectors
70 201 : if (rotation_around_axis != 0)
71 161 : for (const auto i : make_range(6))
72 : {
73 : const Point rotated =
74 138 : _rotation_matrix * Point(_unit_translation_x[i], _unit_translation_y[i], 0);
75 138 : _unit_translation_x[i] = rotated(0);
76 138 : _unit_translation_y[i] = rotated(1);
77 : }
78 :
79 : // compute number of each pin and channel and channel type (interior, edge channel)
80 201 : computePinAndChannelTypes();
81 :
82 201 : _corner_edge_length = (_bundle_side_length - _n_edge_channels / NUM_SIDES * _pin_pitch) / 2.0;
83 :
84 201 : computePinBundleSpacing();
85 201 : computeWireVolumeAndAreaPerPitch();
86 201 : computeFlowVolumes();
87 201 : computeWettedAreas();
88 201 : computeHydraulicDiameters();
89 201 : computePinAndDuctCoordinates();
90 201 : computeChannelPinIndices();
91 201 : computeGapIndices();
92 :
93 201 : if (_pin_bundle_spacing < _wire_diameter)
94 0 : mooseError("Specified bundle pitch " + std::to_string(_bundle_pitch) +
95 : " is too small to fit the given pins and wire wrap!");
96 201 : }
97 :
98 : unsigned int
99 6849863215 : HexagonalLatticeUtils::pins(const unsigned int n) const
100 : {
101 6849863215 : if (n <= 0)
102 : return 0;
103 6849863214 : else if (n == 1)
104 : return 1;
105 : else
106 6849808980 : return NUM_SIDES * (n - 1);
107 : }
108 :
109 : unsigned int
110 53916 : HexagonalLatticeUtils::totalPins(const unsigned int n) const
111 : {
112 : unsigned int total = 0;
113 1277569005 : for (unsigned int i = 1; i <= n; ++i)
114 1277515089 : total += pins(i);
115 :
116 53916 : return total;
117 : }
118 :
119 : unsigned int
120 5 : HexagonalLatticeUtils::rings(const unsigned int n) const
121 : {
122 : auto remaining = n;
123 : unsigned int i = 0;
124 :
125 5572346283 : while (remaining > 0)
126 : {
127 5572346278 : i += 1;
128 5572346278 : remaining -= pins(i);
129 : }
130 :
131 5 : if (n != totalPins(i))
132 2 : mooseError("Number of pins " + std::to_string(n) +
133 : " not evenly divisible in a hexagonal lattice!");
134 :
135 4 : return i;
136 : }
137 :
138 : Real
139 28620 : HexagonalLatticeUtils::hexagonSide(const Real pitch) const
140 : {
141 28620 : return pitch / std::sqrt(3.0);
142 : }
143 :
144 : Real
145 201 : HexagonalLatticeUtils::hexagonArea(const Real pitch) const
146 : {
147 201 : Real side = hexagonSide(pitch);
148 201 : return 3 * SIN60 * side * side;
149 : }
150 :
151 : Real
152 201 : HexagonalLatticeUtils::hexagonVolume(const Real pitch) const
153 : {
154 201 : return hexagonArea(pitch) * _wire_pitch;
155 : }
156 :
157 : Real
158 0 : HexagonalLatticeUtils::hexagonPitch(const Real volume) const
159 : {
160 0 : Real area = volume / _wire_pitch;
161 0 : return std::sqrt(area / SIN60);
162 : }
163 :
164 : Real
165 201 : HexagonalLatticeUtils::triangleArea(const Real side) const
166 : {
167 201 : return 0.5 * triangleHeight(side) * side;
168 : }
169 :
170 : Real
171 402 : HexagonalLatticeUtils::triangleHeight(const Real side) const
172 : {
173 402 : return SIN60 * side;
174 : }
175 :
176 : Real
177 0 : HexagonalLatticeUtils::triangleSide(const Real height) const
178 : {
179 0 : return height / SIN60;
180 : }
181 :
182 : Real
183 201 : HexagonalLatticeUtils::triangleVolume(const Real side) const
184 : {
185 201 : return triangleArea(side) * _wire_pitch;
186 : }
187 :
188 : Real
189 755 : HexagonalLatticeUtils::pinRadius() const
190 : {
191 755 : return _pin_diameter / 2.0;
192 : }
193 :
194 : channel_type::ChannelTypeEnum
195 136 : HexagonalLatticeUtils::channelType(const Point & p) const
196 : {
197 136 : const Real distance_to_wall = minDuctWallDistance(p);
198 :
199 136 : if (distance_to_wall > _pin_bundle_spacing + pinRadius())
200 : return channel_type::interior;
201 :
202 : // the point is in a corner channel if the right triangle formed by the point,
203 : // the corner, and the duct wall has a wall length less than the wall length
204 : // of a corner channel
205 58 : const Real distance_to_corner = minDuctCornerDistance(p);
206 58 : Real l = std::sqrt(distance_to_corner * distance_to_corner - distance_to_wall * distance_to_wall);
207 :
208 58 : if (l < _corner_edge_length)
209 : return channel_type::corner;
210 : else
211 38 : return channel_type::edge;
212 : }
213 :
214 : Real
215 0 : HexagonalLatticeUtils::channelSpecificSurfaceArea(
216 : const channel_type::ChannelTypeEnum & channel) const
217 : {
218 0 : switch (channel)
219 : {
220 0 : case channel_type::interior:
221 0 : return _interior_wetted_area / _interior_volume;
222 0 : case channel_type::edge:
223 0 : return _edge_wetted_area / _edge_volume;
224 0 : case channel_type::corner:
225 0 : return _corner_wetted_area / _corner_volume;
226 0 : default:
227 0 : mooseError("Unhandled ChannelTypeEnum in HexagonalLattice!");
228 : }
229 : }
230 :
231 : Real
232 0 : HexagonalLatticeUtils::channelHydraulicDiameter(const channel_type::ChannelTypeEnum & channel) const
233 : {
234 0 : switch (channel)
235 : {
236 0 : case channel_type::interior:
237 0 : return _interior_Dh;
238 0 : case channel_type::edge:
239 0 : return _edge_Dh;
240 0 : case channel_type::corner:
241 0 : return _corner_Dh;
242 0 : default:
243 0 : mooseError("Unhandled ChannelTypeEnum in HexagonalLattice!");
244 : }
245 : }
246 :
247 : Real
248 136 : HexagonalLatticeUtils::minDuctWallDistance(const Point & p) const
249 : {
250 136 : Real distance = std::numeric_limits<Real>::max();
251 952 : for (const auto i : make_range(NUM_SIDES))
252 : {
253 816 : Real a = _duct_coeffs[i][0];
254 816 : Real b = _duct_coeffs[i][1];
255 816 : Real c = _duct_coeffs[i][2];
256 :
257 816 : Real d = std::abs(a * p(_ix) + b * p(_iy) + c) / std::sqrt(a * a + b * b);
258 816 : distance = std::min(d, distance);
259 : }
260 :
261 136 : return distance;
262 : }
263 :
264 : Real
265 58 : HexagonalLatticeUtils::minDuctCornerDistance(const Point & p) const
266 : {
267 58 : return geom_utils::minDistanceToPoints(p, _duct_corners, _axis);
268 : }
269 :
270 : void
271 201 : HexagonalLatticeUtils::computePinAndDuctCoordinates()
272 : {
273 201 : Real corner_shiftx[] = {COS60, -COS60, -1, -COS60, COS60, 1};
274 201 : Real corner_shifty[] = {SIN60, SIN60, 0, -SIN60, -SIN60, 0};
275 :
276 201 : Real edge_shiftx[] = {-1, -COS60, COS60, 1, COS60, -COS60};
277 201 : Real edge_shifty[] = {0, -SIN60, -SIN60, 0, SIN60, SIN60};
278 :
279 : // Use rotation matrix on shift vectors
280 201 : if (_rotation_around_axis != 0)
281 161 : for (const auto i : make_range(6))
282 : {
283 138 : Point rotated = _rotation_matrix * Point(corner_shiftx[i], corner_shifty[i], 0);
284 138 : corner_shiftx[i] = rotated(0);
285 138 : corner_shifty[i] = rotated(1);
286 138 : rotated = _rotation_matrix * Point(edge_shiftx[i], edge_shifty[i], 0);
287 138 : edge_shiftx[i] = rotated(0);
288 138 : edge_shifty[i] = rotated(1);
289 : }
290 :
291 : // compute coordinates of the pin centers relative to the bundle's center
292 : Point center(0, 0, 0);
293 201 : _pin_centers.push_back(center);
294 :
295 472 : for (unsigned int i = 2; i <= _n_rings; ++i)
296 : {
297 271 : auto n_total_in_ring = pins(i);
298 271 : auto increment = i - 1;
299 :
300 : unsigned int d = 0;
301 :
302 2527 : for (const auto j : make_range(n_total_in_ring))
303 : {
304 2256 : unsigned int side = j / increment;
305 :
306 2256 : if (d == increment)
307 : d = 0;
308 :
309 2256 : Point center = geom_utils::projectPoint(corner_shiftx[side] * _pin_pitch * (i - 1),
310 2256 : corner_shifty[side] * _pin_pitch * (i - 1),
311 2256 : _axis);
312 :
313 : // additional shift for the edge sides
314 2256 : if (d != 0)
315 : {
316 630 : center(_ix) += edge_shiftx[side] * _pin_pitch * d;
317 630 : center(_iy) += edge_shifty[side] * _pin_pitch * d;
318 : }
319 :
320 2256 : _pin_centers.push_back(center);
321 :
322 2256 : d += 1;
323 : }
324 : }
325 :
326 : // compute corners of the hexagonal prisms that enclose each pin
327 201 : _pin_centered_corner_coordinates.resize(_n_pins);
328 201 : auto side = hexagonSide(_pin_pitch);
329 2658 : for (const auto pin : make_range(_n_pins))
330 : {
331 17199 : for (const auto i : make_range(NUM_SIDES))
332 : {
333 14742 : Point translation = geom_utils::projectPoint(
334 14742 : _unit_translation_x[i] * side, _unit_translation_y[i] * side, _axis);
335 :
336 14742 : _pin_centered_corner_coordinates[pin].push_back(translation + _pin_centers[pin]);
337 : }
338 : }
339 :
340 : // compute coordinates of duct corners relative to the bundle's center
341 201 : Real l = _bundle_side_length;
342 1407 : for (const auto i : make_range(NUM_SIDES))
343 : {
344 1206 : Point corner = geom_utils::projectPoint(corner_shiftx[i] * l, corner_shifty[i] * l, _axis);
345 1206 : _duct_corners.push_back(corner);
346 : }
347 :
348 : // compute the equations (a*x + b*y + c) defining each duct wall
349 1407 : for (const auto i : make_range(NUM_SIDES))
350 : {
351 : auto c = i;
352 1206 : unsigned int n = i == 5 ? 0 : i + 1;
353 1206 : Real slope = (_duct_corners[n](_iy) - _duct_corners[c](_iy)) /
354 1206 : (_duct_corners[n](_ix) - _duct_corners[c](_ix));
355 1206 : std::vector<Real> coeffs = {-slope, 1.0, slope * _duct_corners[c](_ix) - _duct_corners[c](_iy)};
356 1206 : _duct_coeffs.push_back(coeffs);
357 1206 : }
358 201 : }
359 :
360 : void
361 201 : HexagonalLatticeUtils::computeWettedAreas()
362 : {
363 201 : Real rod_area_per_pitch = _pin_surface_area_per_pitch + _wire_surface_area_per_pitch;
364 :
365 201 : _interior_wetted_area = 0.5 * rod_area_per_pitch;
366 201 : _edge_wetted_area = 0.5 * rod_area_per_pitch + _pin_pitch * _wire_pitch;
367 :
368 201 : Real h = (_bundle_side_length - _n_edge_channels / NUM_SIDES * _pin_pitch) / 2.0;
369 201 : _corner_wetted_area = rod_area_per_pitch / 6.0 + 2 * h * _wire_pitch;
370 201 : _wetted_area = _n_pins * rod_area_per_pitch + 6.0 * _bundle_side_length * _wire_pitch;
371 201 : }
372 :
373 : void
374 201 : HexagonalLatticeUtils::computeHydraulicDiameters()
375 : {
376 : // Because the hydraulic diameter should be different for a bundle with an
377 : // infinite wire pitch vs. a finite wire pitch (because finite = more wire spooled in),
378 : // we compute hydraulic diameters based on volumes and areas.
379 201 : _interior_Dh = 4 * _interior_flow_volume / _interior_wetted_area;
380 201 : _edge_Dh = 4 * _edge_flow_volume / _edge_wetted_area;
381 201 : _corner_Dh = 4 * _corner_flow_volume / _corner_wetted_area;
382 201 : _Dh = 4 * _flow_volume / _wetted_area;
383 201 : }
384 :
385 : void
386 201 : HexagonalLatticeUtils::computeFlowVolumes()
387 : {
388 : // It is assumed that in the interior and edge channels, that a wire is present half
389 : // of the time (because those channels contain half a pin), while in the corner channels
390 : // that a wire is present one sixth of the time (because those channels contain one sixth
391 : // of a pin).
392 201 : Real rod_volume_per_pitch = _pin_volume_per_pitch + _wire_volume_per_pitch;
393 :
394 201 : Real l = _pin_bundle_spacing + pinRadius();
395 :
396 201 : _interior_volume = triangleVolume(_pin_pitch);
397 201 : _edge_volume = _pin_pitch * l * _wire_pitch;
398 201 : _corner_volume = l * _corner_edge_length * _wire_pitch;
399 :
400 201 : _interior_flow_volume = _interior_volume - 0.5 * rod_volume_per_pitch;
401 201 : _edge_flow_volume = _edge_volume - 0.5 * rod_volume_per_pitch;
402 201 : _corner_flow_volume = _corner_volume - rod_volume_per_pitch / 6.0;
403 :
404 201 : _flow_volume = hexagonVolume(_bundle_pitch) - _n_pins * rod_volume_per_pitch;
405 201 : }
406 :
407 : unsigned int
408 542 : HexagonalLatticeUtils::interiorChannels(const unsigned int ring)
409 : {
410 542 : return (ring * 2 - 1) * NUM_SIDES;
411 : }
412 :
413 : void
414 201 : HexagonalLatticeUtils::computePinAndChannelTypes()
415 : {
416 201 : _n_pins = 0;
417 673 : for (unsigned int i = _n_rings; i > 0; --i)
418 472 : _n_pins += pins(i);
419 :
420 : // the one-ring case
421 201 : _n_interior_pins = _n_pins;
422 201 : _n_edge_pins = 0;
423 201 : _n_corner_pins = 0;
424 :
425 201 : _n_corner_channels = NUM_SIDES;
426 201 : _n_edge_channels = 0;
427 201 : _n_interior_channels = 0;
428 :
429 201 : if (_n_rings > 1)
430 : {
431 191 : _n_corner_pins = NUM_SIDES;
432 191 : _n_edge_pins = (_n_rings - 2) * NUM_SIDES;
433 191 : _n_interior_pins -= _n_corner_pins + _n_edge_pins;
434 :
435 191 : _n_edge_channels = _n_edge_pins + NUM_SIDES;
436 :
437 462 : for (unsigned int i = 1; i < _n_rings; ++i)
438 271 : _n_interior_channels += interiorChannels(i);
439 : }
440 :
441 201 : _n_channels = _n_interior_channels + _n_edge_channels + _n_corner_channels;
442 201 : }
443 :
444 : void
445 201 : HexagonalLatticeUtils::computePinBundleSpacing()
446 : {
447 : // hexagonal flat-to-flat that *just* encloses the pins (NOT including the wire)
448 201 : Real pin_h = 2.0 * (_n_rings - 1) * triangleHeight(_pin_pitch) + _pin_diameter;
449 :
450 201 : _pin_bundle_spacing = (_bundle_pitch - pin_h) / 2.0;
451 201 : }
452 :
453 : void
454 201 : HexagonalLatticeUtils::computeWireVolumeAndAreaPerPitch()
455 : {
456 : unsigned int n_wire_coils = 1;
457 : Real wire_length_per_coil =
458 201 : std::sqrt(_wire_pitch * _wire_pitch + std::pow(M_PI * (_pin_diameter + _wire_diameter), 2.0));
459 : Real wire_length = wire_length_per_coil * n_wire_coils;
460 201 : _wire_volume_per_pitch = _wire_area * wire_length;
461 201 : _wire_surface_area_per_pitch = _wire_circumference * wire_length;
462 201 : }
463 :
464 : void
465 201 : HexagonalLatticeUtils::computeChannelPinIndices()
466 : {
467 201 : _interior_channel_pin_indices.resize(_n_interior_channels);
468 201 : _edge_channel_pin_indices.resize(_n_edge_channels);
469 201 : _corner_channel_pin_indices.resize(_n_corner_channels);
470 :
471 : // 3 pins per interior channel
472 3087 : for (auto & i : _interior_channel_pin_indices)
473 2886 : i.resize(3);
474 :
475 : // 2 pins per edge channel
476 1827 : for (auto & i : _edge_channel_pin_indices)
477 1626 : i.resize(2);
478 :
479 : // 1 pin per corner channel
480 1407 : for (auto & i : _corner_channel_pin_indices)
481 1206 : i.resize(1);
482 :
483 : // for each ring of pins, and for each sector, get the "first" pin index in that ring
484 : std::vector<std::vector<unsigned int>> starts;
485 201 : starts.resize(_n_rings);
486 673 : for (const auto i : make_range(_n_rings))
487 : {
488 472 : starts[i].resize(NUM_SIDES);
489 472 : starts[i][0] = (i == 0) ? 0 : totalPins(i);
490 :
491 2832 : for (unsigned int j = 1; j < NUM_SIDES; ++j)
492 2360 : starts[i][j] = starts[i][j - 1] + i;
493 : }
494 :
495 : unsigned int c = 0;
496 472 : for (const auto i : make_range(_n_rings - 1))
497 : {
498 271 : auto channels = interiorChannels(i + 1);
499 271 : unsigned int channels_per_sector = channels / NUM_SIDES;
500 :
501 1897 : for (const auto j : make_range(NUM_SIDES))
502 : {
503 1626 : auto prev_inner = starts[i][j];
504 1626 : auto prev_outer = starts[i + 1][j];
505 :
506 4512 : for (const auto k : make_range(channels_per_sector))
507 : {
508 2886 : bool downward = k % 2 == 0;
509 :
510 2886 : if (downward)
511 : {
512 2256 : _interior_channel_pin_indices[c][0] = prev_inner;
513 2256 : _interior_channel_pin_indices[c][1] = prev_outer;
514 2256 : _interior_channel_pin_indices[c][2] = ++prev_outer;
515 : }
516 : else
517 : {
518 630 : _interior_channel_pin_indices[c][0] = prev_outer;
519 630 : _interior_channel_pin_indices[c][1] = ++prev_inner;
520 630 : _interior_channel_pin_indices[c][2] = prev_inner - 1;
521 : }
522 :
523 2886 : if (j == 5) // last sector
524 : {
525 481 : if (k == channels_per_sector - 1)
526 : {
527 271 : _interior_channel_pin_indices[c][2] = starts[i + 1][0];
528 271 : _interior_channel_pin_indices[c][0] = starts[i][0];
529 : }
530 :
531 481 : if (k == channels_per_sector - 2 && i > 0)
532 80 : _interior_channel_pin_indices[c][1] = starts[i][0];
533 : }
534 :
535 2886 : c++;
536 : }
537 : }
538 : }
539 :
540 1827 : for (const auto i : make_range(_n_edge_channels))
541 : {
542 1626 : _edge_channel_pin_indices[i][0] = starts[_n_rings - 1][0] + i;
543 1626 : _edge_channel_pin_indices[i][1] = _edge_channel_pin_indices[i][0] + 1;
544 :
545 1626 : if (i == _n_edge_channels - 1)
546 191 : _edge_channel_pin_indices[i][1] = _edge_channel_pin_indices[0][0];
547 : }
548 :
549 1407 : for (const auto i : make_range(_n_corner_channels))
550 1206 : _corner_channel_pin_indices[i][0] =
551 1206 : totalPins(_n_rings - 1) + i * (_n_edge_channels / NUM_SIDES);
552 201 : }
553 :
554 : const std::vector<Point>
555 652 : HexagonalLatticeUtils::interiorChannelCornerCoordinates(
556 : const unsigned int interior_channel_id) const
557 : {
558 : std::vector<Point> corners;
559 652 : auto pin_indices = _interior_channel_pin_indices[interior_channel_id];
560 2608 : for (const auto & pin : pin_indices)
561 1956 : corners.push_back(_pin_centers[pin]);
562 :
563 652 : return corners;
564 652 : }
565 :
566 : const std::vector<Point>
567 166 : HexagonalLatticeUtils::edgeChannelCornerCoordinates(const unsigned int edge_channel_id) const
568 : {
569 : std::vector<Point> corners;
570 :
571 166 : auto pin_indices = _edge_channel_pin_indices[edge_channel_id];
572 :
573 166 : const Point & pin1 = _pin_centers[pin_indices[0]];
574 166 : const Point & pin2 = _pin_centers[pin_indices[1]];
575 :
576 166 : Real d = pinBundleSpacing() + pinRadius();
577 :
578 166 : corners.push_back(pin1);
579 166 : corners.push_back(pin2);
580 :
581 166 : unsigned int sector = edge_channel_id / (_n_rings - 1);
582 166 : corners.push_back(pin2 +
583 166 : Point(d * _unit_translation_x[sector], d * _unit_translation_y[sector], 0.0));
584 166 : corners.push_back(pin1 +
585 166 : Point(d * _unit_translation_x[sector], d * _unit_translation_y[sector], 0.0));
586 :
587 166 : return corners;
588 166 : }
589 :
590 : const std::vector<Point>
591 50 : HexagonalLatticeUtils::cornerChannelCornerCoordinates(const unsigned int corner_channel_id) const
592 : {
593 : std::vector<Point> corners;
594 :
595 50 : auto pin_indices = _corner_channel_pin_indices[corner_channel_id];
596 50 : const Point & pin = _pin_centers[pin_indices[0]];
597 50 : corners.push_back(pin);
598 :
599 50 : Real d = pinBundleSpacing() + pinRadius();
600 :
601 50 : unsigned int side1 = corner_channel_id == 0 ? NUM_SIDES - 1 : corner_channel_id - 1;
602 : unsigned int side2 = corner_channel_id;
603 :
604 50 : corners.push_back(pin +
605 50 : Point(d * _unit_translation_x[side1], d * _unit_translation_y[side1], 0.0));
606 50 : corners.push_back(_duct_corners[corner_channel_id]);
607 50 : corners.push_back(pin +
608 50 : Point(d * _unit_translation_x[side2], d * _unit_translation_y[side2], 0.0));
609 :
610 50 : return corners;
611 50 : }
612 :
613 : Point
614 0 : HexagonalLatticeUtils::channelCentroid(const std::vector<Point> & pins) const
615 : {
616 0 : int n_pins = pins.size();
617 : Point centroid(0.0, 0.0, 0.0);
618 0 : for (const auto & p : pins)
619 : centroid += p / n_pins;
620 :
621 0 : return centroid;
622 : }
623 :
624 : unsigned int
625 28017 : HexagonalLatticeUtils::pinIndex(const Point & point) const
626 : {
627 28017 : auto side = hexagonSide(_pin_pitch);
628 :
629 440065 : for (const auto i : make_range(_n_pins))
630 : {
631 431532 : const auto & center = _pin_centers[i];
632 431532 : Real dx = center(_ix) - point(_ix);
633 431532 : Real dy = center(_iy) - point(_iy);
634 431532 : Real distance_from_pin = std::sqrt(dx * dx + dy * dy);
635 :
636 : // if we're outside the circumference of the hexagon, we're certain not
637 : // within the hexagon for this pin
638 431532 : if (distance_from_pin > side)
639 410859 : continue;
640 :
641 20673 : auto corners = _pin_centered_corner_coordinates[i];
642 20673 : if (geom_utils::pointInPolygon(point, corners, _axis))
643 : return i;
644 20673 : }
645 :
646 8533 : return _n_pins;
647 : }
648 :
649 : unsigned int
650 0 : HexagonalLatticeUtils::closestPinIndex(const Point & point) const
651 : {
652 : // If within the lattice, you must consider all pins. If outside, outer ring suffices
653 : unsigned int start_index = 0;
654 0 : if (!insideLattice(point))
655 0 : start_index = firstPinInRing(_n_rings);
656 :
657 : // Classic minimum search. If more performance is required for a large lattice we may consider
658 : // using a KD-Tree instead or checking which ring the point is part of before examining each pin
659 : auto min_distance = std::numeric_limits<Real>::max();
660 : unsigned int index_min = 0;
661 0 : for (unsigned int i = start_index; i < _n_pins; ++i)
662 : {
663 0 : const auto & center = _pin_centers[i];
664 0 : Real dx = center(_ix) - point(_ix);
665 0 : Real dy = center(_iy) - point(_iy);
666 0 : Real distance_from_pin = std::sqrt(dx * dx + dy * dy);
667 :
668 0 : if (distance_from_pin < min_distance)
669 : {
670 : min_distance = distance_from_pin;
671 : index_min = i;
672 : }
673 : }
674 :
675 0 : return index_min;
676 : }
677 :
678 : unsigned int
679 94 : HexagonalLatticeUtils::channelIndex(const Point & point) const
680 : {
681 94 : auto channel = channelType(point);
682 :
683 94 : switch (channel)
684 : {
685 54 : case channel_type::interior:
686 : {
687 652 : for (const auto i : make_range(_n_interior_channels))
688 : {
689 652 : auto corners = interiorChannelCornerCoordinates(i);
690 652 : if (geom_utils::pointInPolygon(point, corners, _axis))
691 : return i;
692 652 : }
693 : break;
694 : }
695 26 : case channel_type::edge:
696 : {
697 166 : for (const auto i : make_range(_n_edge_channels))
698 : {
699 166 : auto corners = edgeChannelCornerCoordinates(i);
700 166 : if (geom_utils::pointInPolygon(point, corners, _axis))
701 26 : return i + _n_interior_channels;
702 166 : }
703 : break;
704 : }
705 14 : case channel_type::corner:
706 : {
707 50 : for (const auto i : make_range(_n_corner_channels))
708 : {
709 50 : auto corners = cornerChannelCornerCoordinates(i);
710 50 : if (geom_utils::pointInPolygon(point, corners, _axis))
711 14 : return i + _n_interior_channels + _n_edge_channels;
712 50 : }
713 : break;
714 : }
715 0 : default:
716 0 : mooseError("Unhandled ChannelTypeEnum!");
717 : }
718 :
719 0 : mooseError(
720 0 : "Point (" + std::to_string(point(0)) + ", " + std::to_string(point(1)) + ", " +
721 0 : std::to_string(point(2)) +
722 : ") is not in any channel! This can sometimes happen "
723 : "due to:\n\n a) Points in the mesh actually being outside the domain specified with the "
724 : "HexagonalLatticeUtils.\n b) Small floating point errors - we recommend using a CONSTANT "
725 : "MONOMIAL variable with all related objects.\nYou can also try slightly decreasing the pin "
726 : "diameter and/or "
727 : "increasing the bundle pitch.");
728 : }
729 :
730 : bool
731 0 : HexagonalLatticeUtils::insideLattice(const Point & point) const
732 : {
733 0 : std::vector<Point> lattice_corners(6);
734 0 : auto side = hexagonSide(_bundle_pitch);
735 0 : for (const auto i : make_range(NUM_SIDES))
736 : {
737 0 : Point translation = geom_utils::projectPoint(
738 0 : _unit_translation_x[i] * side, _unit_translation_y[i] * side, _axis);
739 0 : lattice_corners[i] = translation;
740 : }
741 0 : return geom_utils::pointInPolygon(point, lattice_corners, _axis);
742 0 : }
743 :
744 : void
745 893 : HexagonalLatticeUtils::get2DInputPatternIndex(const unsigned int pin_index,
746 : unsigned int & row_index,
747 : unsigned int & within_row_index) const
748 : {
749 :
750 : // First compute the ring on which the pin is
751 893 : const auto ring_index = ringIndex(pin_index);
752 893 : row_index = _n_rings - ring_index;
753 893 : within_row_index = _n_rings - 1;
754 :
755 893 : const auto n_pin_sectors = pins(ring_index) / 6;
756 : // Loop around the center until we reach the pin_index
757 : // We start on the diagonal at the top
758 :
759 3572 : for (const auto local_i : make_range(pin_index - firstPinInRing(ring_index)))
760 : {
761 2679 : if (local_i < n_pin_sectors)
762 874 : within_row_index--;
763 1805 : else if (local_i < 3 * n_pin_sectors)
764 1235 : row_index++;
765 570 : else if (local_i < 4 * n_pin_sectors)
766 361 : within_row_index++;
767 209 : else if (local_i < 5 * n_pin_sectors)
768 : {
769 190 : row_index--;
770 190 : within_row_index++;
771 : }
772 : else
773 : {
774 19 : row_index--;
775 19 : within_row_index--;
776 : }
777 : }
778 :
779 : // Underflow (row_index is invalid_uint) will also fail this assert
780 : mooseAssert(row_index < 2 * _n_rings - 1, "Inverse indexing failed");
781 : mooseAssert(within_row_index < 2 * _n_rings - 1, "Inverse indexing failed");
782 893 : }
783 :
784 : void
785 201 : HexagonalLatticeUtils::computeGapIndices()
786 : {
787 : std::set<std::pair<int, int>> indices;
788 3087 : for (const auto & pins : _interior_channel_pin_indices)
789 : {
790 : std::pair<int, int> gap0 = {std::min(pins[0], pins[1]), std::max(pins[0], pins[1])};
791 : std::pair<int, int> gap1 = {std::min(pins[0], pins[2]), std::max(pins[0], pins[2])};
792 : std::pair<int, int> gap2 = {std::min(pins[2], pins[1]), std::max(pins[2], pins[1])};
793 :
794 2886 : indices.insert(gap0);
795 2886 : indices.insert(gap1);
796 2886 : indices.insert(gap2);
797 : }
798 :
799 5343 : for (const auto & it : indices)
800 5142 : _gap_indices.push_back({it.first, it.second});
801 :
802 201 : _n_interior_gaps = _gap_indices.size();
803 :
804 : // add the gaps along the peripheral channels; -1 indicates side 1, -2 indicates side 2,
805 : // and so on
806 201 : int n_edge_gaps = _n_rings;
807 201 : int pin = totalPins(_n_rings - 1);
808 1407 : for (const auto i : make_range(NUM_SIDES))
809 : {
810 4038 : for (const auto j : make_range(_n_rings))
811 2832 : _gap_indices.push_back({pin + j, -(i + 1)});
812 :
813 1206 : pin += n_edge_gaps - 1;
814 : }
815 :
816 : // fix the last gap to use the first pin
817 201 : _gap_indices.back() = {totalPins(_n_rings - 1), -int(NUM_SIDES)};
818 201 : _n_gaps = _gap_indices.size();
819 :
820 : // for each channel, determine which gaps are on that channel to find the local-to-global
821 : // indexing
822 3087 : for (const auto & pins : _interior_channel_pin_indices)
823 : {
824 : std::pair<int, int> gap0 = {std::min(pins[0], pins[1]), std::max(pins[0], pins[1])};
825 : std::pair<int, int> gap1 = {std::min(pins[1], pins[2]), std::max(pins[1], pins[2])};
826 : std::pair<int, int> gap2 = {std::min(pins[2], pins[0]), std::max(pins[2], pins[0])};
827 :
828 2886 : int loc_gap0 = globalGapIndex(gap0);
829 2886 : int loc_gap1 = globalGapIndex(gap1);
830 2886 : int loc_gap2 = globalGapIndex(gap2);
831 5772 : _local_to_global_gaps.push_back({loc_gap0, loc_gap1, loc_gap2});
832 : }
833 :
834 201 : int gap = _gap_indices.size() - _n_rings * NUM_SIDES;
835 1827 : for (const auto i : make_range(_n_edge_channels))
836 : {
837 1626 : const auto & pins = _edge_channel_pin_indices[i];
838 : std::pair<int, int> gap0 = {std::min(pins[0], pins[1]), std::max(pins[0], pins[1])};
839 1626 : int loc_gap0 = globalGapIndex(gap0);
840 3252 : _local_to_global_gaps.push_back({loc_gap0, gap + 1, gap});
841 :
842 1626 : if ((i + 1) % (_n_rings - 1) == 0)
843 1146 : gap += 2;
844 : else
845 : gap += 1;
846 : }
847 :
848 201 : int n_interior_gaps = _gap_indices.size() - _n_rings * NUM_SIDES - 1;
849 201 : n_edge_gaps = _n_rings * NUM_SIDES;
850 402 : _local_to_global_gaps.push_back({n_interior_gaps + n_edge_gaps, n_interior_gaps + 1});
851 201 : gap = n_interior_gaps + _n_rings;
852 1206 : for (unsigned int i = 1; i < _n_corner_channels; ++i)
853 : {
854 2010 : _local_to_global_gaps.push_back({gap, gap + 1});
855 1005 : gap += _n_rings;
856 : }
857 :
858 201 : _gap_points.resize(_n_gaps);
859 :
860 : // For each gap, get two points on the gap
861 5343 : for (const auto i : make_range(_n_interior_gaps))
862 : {
863 5142 : const auto & pins = _gap_indices[i];
864 5142 : Point pt1(_pin_centers[pins.first]);
865 5142 : Point pt2(_pin_centers[pins.second]);
866 5142 : _gap_centers.push_back(0.5 * (pt2 + pt1));
867 :
868 5142 : _gap_points[i] = {pt1, pt2};
869 :
870 : // for the last gap in the ring, we need to swap the ordering of pins
871 5142 : if (lastGapInRing(i))
872 : {
873 : Point tmp = pt1;
874 : pt1 = pt2;
875 : pt2 = tmp;
876 : }
877 :
878 10284 : _gap_unit_normals.push_back(geom_utils::projectedUnitNormal(pt1, pt2, _axis));
879 : }
880 :
881 201 : Real d = _pin_bundle_spacing + pinRadius();
882 3033 : for (unsigned int i = _n_interior_gaps; i < _n_gaps; ++i)
883 : {
884 2832 : const auto & pins = _gap_indices[i];
885 2832 : int side = std::abs(pins.second) - 1;
886 :
887 2832 : const auto & pt1 = _pin_centers[pins.first];
888 : const Point pt2 =
889 2832 : pt1 + Point(d * _unit_translation_x[side], d * _unit_translation_y[side], 0.0);
890 2832 : _gap_centers.push_back(0.5 * (pt2 + pt1));
891 :
892 2832 : _gap_points[i] = {pt1, pt2};
893 :
894 5664 : _gap_unit_normals.push_back(geom_utils::projectedUnitNormal(pt1, pt2, _axis));
895 : }
896 201 : }
897 :
898 : unsigned int
899 10284 : HexagonalLatticeUtils::globalGapIndex(const std::pair<int, int> & local_gap) const
900 : {
901 280146 : for (const auto i : index_range(_gap_indices))
902 : {
903 280146 : const auto gap = _gap_indices[i];
904 280146 : if (gap.first == local_gap.first && gap.second == local_gap.second)
905 10284 : return i;
906 : }
907 :
908 0 : mooseError("Failed to find local gap in global gap array!");
909 : }
910 :
911 : Real
912 31 : HexagonalLatticeUtils::distanceFromGap(const Point & pt, const unsigned int gap_index) const
913 : {
914 31 : auto p = _gap_points[gap_index];
915 62 : return geom_utils::projectedDistanceFromLine(pt, p[0], p[1], _axis);
916 31 : }
917 :
918 : unsigned int
919 11 : HexagonalLatticeUtils::gapIndex(const Point & point) const
920 : {
921 11 : const auto & channel_index = channelIndex(point);
922 11 : const auto & gap_indices = _local_to_global_gaps[channel_index];
923 :
924 : Real distance = std::numeric_limits<Real>::max();
925 : unsigned int index = 0;
926 42 : for (const auto i : index_range(gap_indices))
927 : {
928 31 : Real distance_from_gap = distanceFromGap(point, gap_indices[i]);
929 :
930 31 : if (distance_from_gap < distance)
931 : {
932 : distance = distance_from_gap;
933 19 : index = gap_indices[i];
934 : }
935 : }
936 :
937 11 : return index;
938 : }
939 :
940 : void
941 0 : HexagonalLatticeUtils::gapIndexAndDistance(const Point & point,
942 : unsigned int & index,
943 : Real & distance) const
944 : {
945 0 : const auto & channel_index = channelIndex(point);
946 0 : const auto & gap_indices = _local_to_global_gaps[channel_index];
947 :
948 0 : distance = std::numeric_limits<Real>::max();
949 0 : for (const auto i : index_range(gap_indices))
950 : {
951 0 : Real distance_from_gap = distanceFromGap(point, gap_indices[i]);
952 :
953 0 : if (distance_from_gap < distance)
954 : {
955 0 : distance = distance_from_gap;
956 0 : index = gap_indices[i];
957 : }
958 : }
959 0 : }
960 :
961 : unsigned int
962 11788 : HexagonalLatticeUtils::firstPinInRing(const unsigned int ring) const
963 : {
964 : mooseAssert(ring > 0, "Ring indexing starts at 1");
965 11788 : if (ring == 1)
966 : return 0;
967 : else
968 11692 : return totalPins(ring - 1);
969 : }
970 :
971 : unsigned int
972 12814 : HexagonalLatticeUtils::lastPinInRing(const unsigned int ring) const
973 : {
974 : mooseAssert(ring > 0, "Ring indexing starts at 1");
975 12814 : if (ring == 1)
976 : return 0;
977 : else
978 11920 : return totalPins(ring) - 1;
979 : }
980 :
981 : unsigned int
982 893 : HexagonalLatticeUtils::ringIndex(const unsigned int pin) const
983 : {
984 1919 : for (const auto i : make_range(_n_rings))
985 : {
986 1919 : if (pin <= lastPinInRing(i + 1))
987 893 : return i + 1;
988 : }
989 0 : mooseError("Should not reach here. Pin index: ", pin, " for ", _n_rings, " rings.");
990 : }
991 :
992 : bool
993 5232 : HexagonalLatticeUtils::lastGapInRing(const unsigned int gap_index) const
994 : {
995 5232 : if (gap_index >= _n_interior_gaps)
996 : return false;
997 :
998 5196 : const auto & pins = _gap_indices[gap_index];
999 :
1000 15814 : for (unsigned int i = 2; i <= _n_rings; ++i)
1001 : {
1002 : bool one_is_first_pin = false;
1003 : bool one_is_last_pin = false;
1004 :
1005 10892 : int first_pin = firstPinInRing(i);
1006 10892 : int last_pin = lastPinInRing(i);
1007 :
1008 10892 : if (pins.first == first_pin || pins.second == first_pin)
1009 : one_is_first_pin = true;
1010 :
1011 10892 : if (pins.first == last_pin || pins.second == last_pin)
1012 : one_is_last_pin = true;
1013 :
1014 10892 : if (one_is_first_pin && one_is_last_pin)
1015 : return true;
1016 : }
1017 :
1018 : return false;
1019 : }
|
