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 "TriSubChannelMesh.h"
11 : #include <cmath>
12 : #include "libmesh/node.h"
13 :
14 : registerMooseObject("SubChannelApp", TriSubChannelMesh);
15 :
16 : InputParameters
17 368 : TriSubChannelMesh::validParams()
18 : {
19 368 : InputParameters params = SubChannelMesh::validParams();
20 368 : params.addClassDescription("Creates an subchannel mesh container for a triangular "
21 : "lattice arrangement");
22 368 : return params;
23 0 : }
24 :
25 184 : TriSubChannelMesh::TriSubChannelMesh(const InputParameters & params)
26 184 : : SubChannelMesh(params), _pin_mesh_exist(false), _duct_mesh_exist(false)
27 : {
28 184 : }
29 :
30 0 : TriSubChannelMesh::TriSubChannelMesh(const TriSubChannelMesh & other_mesh)
31 : : SubChannelMesh(other_mesh),
32 0 : _n_rings(other_mesh._n_rings),
33 0 : _n_channels(other_mesh._n_channels),
34 0 : _flat_to_flat(other_mesh._flat_to_flat),
35 0 : _dwire(other_mesh._dwire),
36 0 : _hwire(other_mesh._hwire),
37 0 : _duct_to_pin_gap(other_mesh._duct_to_pin_gap),
38 0 : _nodes(other_mesh._nodes),
39 0 : _duct_nodes(other_mesh._duct_nodes),
40 : _chan_to_duct_node_map(other_mesh._chan_to_duct_node_map),
41 : _duct_node_to_chan_map(other_mesh._duct_node_to_chan_map),
42 0 : _gap_to_chan_map(other_mesh._gap_to_chan_map),
43 0 : _gap_to_pin_map(other_mesh._gap_to_pin_map),
44 0 : _chan_to_gap_map(other_mesh._chan_to_gap_map),
45 0 : _sign_id_crossflow_map(other_mesh._sign_id_crossflow_map),
46 0 : _gij_map(other_mesh._gij_map),
47 0 : _pin_position(other_mesh._pin_position),
48 0 : _pins_in_rings(other_mesh._pins_in_rings),
49 0 : _chan_to_pin_map(other_mesh._chan_to_pin_map),
50 0 : _npins(other_mesh._npins),
51 0 : _n_gaps(other_mesh._n_gaps),
52 0 : _subch_type(other_mesh._subch_type),
53 0 : _gap_type(other_mesh._gap_type),
54 0 : _gap_pairs_sf(other_mesh._gap_pairs_sf),
55 0 : _chan_pairs_sf(other_mesh._chan_pairs_sf),
56 0 : _pin_to_chan_map(other_mesh._pin_to_chan_map),
57 0 : _pin_mesh_exist(other_mesh._pin_mesh_exist),
58 0 : _duct_mesh_exist(other_mesh._duct_mesh_exist)
59 : {
60 0 : _subchannel_position = other_mesh._subchannel_position;
61 0 : }
62 :
63 : std::unique_ptr<MooseMesh>
64 0 : TriSubChannelMesh::safeClone() const
65 : {
66 0 : return _app.getFactory().copyConstruct(*this);
67 : }
68 :
69 : unsigned int
70 970632 : TriSubChannelMesh::getSubchannelIndexFromPoint(const Point & p) const
71 : {
72 : /// Function that returns the subchannel index given a point
73 970632 : return this->channelIndex(p);
74 : }
75 :
76 : unsigned int
77 1675692 : TriSubChannelMesh::channelIndex(const Point & p) const
78 : {
79 : /// Function that returns the subchannel index given a point
80 : /// Determining a channel index given a point
81 : /// Looping over all subchannels to determine the closest one to the point
82 : /// Special treatment for edge and corner subchannels since deformed elements may lead to wrong transfers
83 :
84 : // Distances to determine the closest subchannel
85 : Real distance0 = 1.0e+8; // Dummy distance to keep updating the closes distance found
86 : Real distance1; // Distance to be updated with the current subchannel distance
87 : unsigned int j = 0; // Index to keep track of the closest point to subchannel
88 :
89 : // Projecting point into hexahedral coordinated to determine if the point belongs to a center
90 : // subchannel
91 1675692 : Real distance_outer_ring = _flat_to_flat / 2 - _duct_to_pin_gap - _pin_diameter / 2;
92 1675692 : Real channel_distance = std::sqrt(std::pow(p(0), 2) + std::pow(p(1), 2));
93 1675692 : Real angle = std::abs(std::atan(p(1) / p(0)));
94 1675692 : Real projection_angle =
95 1675692 : angle - libMesh::pi / 6 - std::trunc(angle / (libMesh::pi / 3)) * (libMesh::pi / 3);
96 1675692 : channel_distance = channel_distance * std::cos(projection_angle);
97 :
98 : // Projecting point on top edge to determine if the point is a corner or edge subchannel by x
99 : // coordinate
100 : Real loc_angle = std::atan(p(1) / p(0));
101 1675692 : if (p(0) <= 0)
102 837846 : loc_angle += libMesh::pi;
103 837846 : else if (p(0) >= 0 && p(1) <= 0)
104 431982 : loc_angle += 2 * libMesh::pi;
105 1675692 : Real rem_ang = std::trunc(loc_angle / (libMesh::pi / 3)) * (libMesh::pi / 3) - libMesh::pi / 3;
106 1675692 : Real x_coord_new = (std::cos(-rem_ang) * p(0) - std::sin(-rem_ang) * p(1));
107 1675692 : Real x_lim = (_n_rings - 1) * _pitch / 2.0;
108 :
109 : // looping over all channels
110 123459300 : for (unsigned int i = 0; i < _n_channels; i++)
111 : {
112 :
113 121783608 : if (_n_rings == 1)
114 : {
115 : Real angle = std::atan(p(1) / p(0));
116 0 : if ((i * libMesh::pi / 6.0 < angle) && (angle <= (i + 1) * libMesh::pi / 6.0))
117 0 : return i;
118 : }
119 : else
120 : {
121 : // Distance from the point to subchannel
122 121783608 : distance1 = std::sqrt(std::pow((p(0) - _subchannel_position[i][0]), 2.0) +
123 121783608 : std::pow((p(1) - _subchannel_position[i][1]), 2.0));
124 :
125 : // If subchannel belongs to center ring
126 121783608 : if (channel_distance < distance_outer_ring)
127 : {
128 86486400 : if ((distance1 < distance0) && (_subch_type[i] == EChannelType::CENTER))
129 : {
130 : j = i;
131 : distance0 = distance1;
132 : } // if
133 : } // if
134 : // If subchannel belongs to outer ring
135 : else
136 : {
137 35297208 : if ((distance1 < distance0) &&
138 26501238 : (_subch_type[i] == EChannelType::EDGE || _subch_type[i] == EChannelType::CORNER))
139 : {
140 3264822 : if (((x_coord_new > x_lim) || (x_coord_new < -x_lim)) &&
141 : _subch_type[i] == EChannelType::CORNER)
142 : {
143 : j = i;
144 : distance0 = distance1;
145 : } // if
146 3039852 : else if (((x_coord_new > x_lim) || (x_coord_new > -x_lim)) &&
147 : _subch_type[i] == EChannelType::EDGE)
148 : {
149 : j = i;
150 : distance0 = distance1;
151 : }
152 : }
153 : }
154 : }
155 : }
156 : return j;
157 : }
158 :
159 : void
160 184 : TriSubChannelMesh::buildMesh()
161 : {
162 184 : }
163 :
164 : unsigned int
165 56382 : TriSubChannelMesh::getPinIndexFromPoint(const Point & p) const
166 : {
167 : /// Function that returns the pin number given a point
168 :
169 56382 : return this->pinIndex(p);
170 : }
171 :
172 : unsigned int
173 56382 : TriSubChannelMesh::pinIndex(const Point & p) const
174 : {
175 : /// Function that returns the pin number given a point
176 :
177 : // Define the current ring
178 : Real distance_rod;
179 : Real d0 = 1e5;
180 : unsigned int current_rod = 0;
181 :
182 : std::vector<Point> positions;
183 : Point center(0, 0);
184 56382 : this->rodPositions(positions, _n_rings, _pitch, center);
185 5168820 : for (unsigned int i = 0; i < _npins; i++)
186 : {
187 5112438 : Real x_dist = positions[i](0) - p(0);
188 5112438 : Real y_dist = positions[i](1) - p(1);
189 5112438 : distance_rod = std::sqrt(std::pow(x_dist, 2) + std::pow(y_dist, 2));
190 5112438 : if (distance_rod < d0)
191 : {
192 : d0 = distance_rod;
193 : current_rod = i;
194 : }
195 : }
196 :
197 56382 : return current_rod;
198 56382 : }
199 :
200 : void
201 56638 : TriSubChannelMesh::rodPositions(std::vector<Point> & positions,
202 : unsigned int nrings,
203 : Real pitch,
204 : Point center)
205 : {
206 : /// Defining parameters
207 : // distance: it is the distance to the next Pin
208 : //
209 : Real theta = 0.0;
210 : Real dtheta = 0.0;
211 : Real distance = 0.0;
212 : Real theta1 = 0.0;
213 : Real theta_corrected = 0.0;
214 : Real pi = libMesh::pi;
215 : unsigned int k = 0;
216 56638 : positions.emplace_back(0.0, 0.0);
217 319046 : for (unsigned int i = 1; i < nrings; i++)
218 : {
219 262408 : dtheta = 2.0 * pi / (i * 6);
220 : theta = 0.0;
221 5326738 : for (unsigned int j = 0; j < i * 6; j++)
222 : {
223 : k = k + 1;
224 5064330 : theta1 = fmod(theta + 1.0e-10, pi / 3.0);
225 5064330 : distance = std::sqrt((pow(i * pitch, 2) + pow(theta1 / dtheta * pitch, 2) -
226 5064330 : 2.0 * i * pitch * (theta1 / dtheta * pitch) * std::cos(pi / 3.0)));
227 5064330 : double argument = 1.0 / (i * pitch) / distance / 2.0 *
228 5064330 : (std::pow(i * pitch, 2.0) + std::pow(distance, 2.0) -
229 5064330 : std::pow(theta1 / dtheta * pitch, 2.0));
230 : // Check if the argument to std::acos() is within the valid range [-1, 1]
231 5064330 : if (argument >= -1.0 && argument <= 1.0)
232 : {
233 4975362 : theta_corrected = std::acos(argument);
234 : }
235 : else if (argument > 1.0)
236 : {
237 : theta_corrected = 0.0;
238 : }
239 : else
240 : {
241 : theta_corrected = pi;
242 : }
243 5064330 : if (theta1 < 1.0e-6)
244 : {
245 : theta_corrected = theta;
246 : }
247 : else
248 : {
249 3489882 : if (theta > pi / 3.0 && theta <= 2.0 / 3.0 * pi)
250 581647 : theta_corrected = theta_corrected + pi / 3.0;
251 2326588 : else if (theta > 2.0 / 3.0 * pi && theta <= pi)
252 581647 : theta_corrected = theta_corrected + 2.0 / 3.0 * pi;
253 1744941 : else if (theta > pi && theta <= 4.0 / 3.0 * pi)
254 581647 : theta_corrected = theta_corrected + pi;
255 1163294 : else if (theta > 4.0 / 3.0 * pi && theta <= 5.0 / 3.0 * pi)
256 581647 : theta_corrected = theta_corrected + 4.0 / 3.0 * pi;
257 581647 : else if (theta > 5.0 / 3.0 * pi && theta <= 2.0 * pi)
258 581647 : theta_corrected = theta_corrected + 5.0 / 3.0 * pi;
259 : }
260 10128660 : positions.emplace_back(center(0) + distance * std::cos(theta_corrected),
261 5064330 : center(1) + distance * std::sin(theta_corrected));
262 5064330 : theta = theta + dtheta;
263 : } // j
264 : } // i
265 56638 : }
266 :
267 : void
268 40 : TriSubChannelMesh::setChannelToDuctMaps(const std::vector<Node *> & duct_nodes)
269 : {
270 : const Real tol = 1e-10;
271 8626 : for (size_t i = 0; i < duct_nodes.size(); i++)
272 : {
273 : int min_chan = 0;
274 : Real min_dist = std::numeric_limits<double>::max();
275 8586 : Point ductpos((*duct_nodes[i])(0), (*duct_nodes[i])(1), 0);
276 1483758 : for (size_t j = 0; j < _subchannel_position.size(); j++)
277 : {
278 1475172 : Point chanpos(_subchannel_position[j][0], _subchannel_position[j][1], 0);
279 1475172 : auto dist = (chanpos - ductpos).norm();
280 1475172 : if (dist < min_dist)
281 : {
282 : min_dist = dist;
283 128067 : min_chan = j;
284 : }
285 : }
286 :
287 8586 : Node * chan_node = nullptr;
288 52326 : for (auto cn : _nodes[min_chan])
289 : {
290 52326 : if (std::abs((*cn)(2) - (*duct_nodes[i])(2)) < tol)
291 : {
292 8586 : chan_node = cn;
293 8586 : break;
294 : }
295 : }
296 :
297 8586 : if (chan_node == nullptr)
298 0 : mooseError("failed to find matching channel node for duct node");
299 :
300 8586 : _duct_node_to_chan_map[duct_nodes[i]] = chan_node;
301 8586 : _chan_to_duct_node_map[chan_node] = duct_nodes[i];
302 : }
303 :
304 40 : _duct_nodes = duct_nodes;
305 40 : }
|
