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 "SCMDetailedTriSubChannelMeshGenerator.h"
11 : #include "TriSubChannelMesh.h"
12 : #include <array>
13 : #include <cmath>
14 : #include "libmesh/cell_prism6.h"
15 : #include "libmesh/unstructured_mesh.h"
16 :
17 : registerMooseObject("SubChannelApp", SCMDetailedTriSubChannelMeshGenerator);
18 : registerMooseObjectRenamed("SubChannelApp",
19 : DetailedTriSubChannelMeshGenerator,
20 : "06/30/2025 24:00",
21 : SCMDetailedTriSubChannelMeshGenerator);
22 :
23 : InputParameters
24 108 : SCMDetailedTriSubChannelMeshGenerator::validParams()
25 : {
26 108 : InputParameters params = MeshGenerator::validParams();
27 108 : params.addClassDescription(
28 : "Creates a detailed mesh of subchannels in a triangular lattice arrangement");
29 216 : params.addRequiredParam<Real>("pitch", "Pitch [m]");
30 216 : params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
31 216 : params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
32 216 : params.addRequiredParam<Real>("heated_length", "Heated length [m]");
33 216 : params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
34 216 : params.addRequiredParam<unsigned int>("nrings", "Number of fuel Pin rings per assembly [-]");
35 216 : params.addRequiredParam<Real>("flat_to_flat",
36 : "Flat to flat distance for the hexagonal assembly [m]");
37 216 : params.addParam<unsigned int>("block_id", 0, "Block ID used for the mesh subdomain.");
38 216 : params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
39 216 : params.addParam<bool>("verbose_flag", false, "Flag to print out the mesh coordinates.");
40 108 : return params;
41 0 : }
42 :
43 54 : SCMDetailedTriSubChannelMeshGenerator::SCMDetailedTriSubChannelMeshGenerator(
44 54 : const InputParameters & parameters)
45 : : MeshGenerator(parameters),
46 54 : _unheated_length_entry(getParam<Real>("unheated_length_entry")),
47 108 : _heated_length(getParam<Real>("heated_length")),
48 108 : _unheated_length_exit(getParam<Real>("unheated_length_exit")),
49 108 : _pitch(getParam<Real>("pitch")),
50 108 : _pin_diameter(getParam<Real>("pin_diameter")),
51 108 : _n_rings(getParam<unsigned int>("nrings")),
52 108 : _flat_to_flat(getParam<Real>("flat_to_flat")),
53 108 : _block_id(getParam<unsigned int>("block_id")),
54 108 : _n_cells(getParam<unsigned int>("n_cells")),
55 162 : _verbose(getParam<bool>("verbose_flag"))
56 : {
57 54 : Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
58 54 : Real dz = L / _n_cells;
59 1908 : for (unsigned int i = 0; i < _n_cells + 1; i++)
60 1854 : _z_grid.push_back(dz * i);
61 :
62 : // x coordinate for the first position
63 : Real x0 = 0.0;
64 : // y coordinate for the first position
65 : Real y0 = 0.0;
66 : // x coordinate for the second position
67 : Real x1 = 0.0;
68 : // y coordinate for the second position dummy variable
69 : Real y1 = 0.0;
70 : // dummy variable
71 : Real a1 = 0.0;
72 : // dummy variable
73 : Real a2 = 0.0;
74 : // average x coordinate
75 : Real avg_coor_x = 0.0;
76 : // average y coordinate
77 : Real avg_coor_y = 0.0;
78 : // distance between two points
79 : Real dist = 0.0;
80 : // distance between two points
81 : Real dist0 = 0.0;
82 : // the indicator used while setting _gap_to_chan_map array
83 : std::vector<std::pair<unsigned int, unsigned int>> gap_fill;
84 54 : TriSubChannelMesh::rodPositions(_pin_position, _n_rings, _pitch, Point(0, 0));
85 54 : _nrods = _pin_position.size();
86 : // assign the pins to the corresponding rings
87 : unsigned int k = 0; // initialize the fuel Pin counter index
88 54 : _pins_in_rings.resize(_n_rings);
89 54 : _pins_in_rings[0].push_back(k++);
90 171 : for (unsigned int i = 1; i < _n_rings; i++)
91 1251 : for (unsigned int j = 0; j < i * 6; j++)
92 1134 : _pins_in_rings[i].push_back(k++);
93 : // Given the number of pins and number of fuel Pin rings, the number of subchannels can be
94 : // computed as follows:
95 : unsigned int chancount = 0.0;
96 : // Summing internal channels
97 171 : for (unsigned int j = 0; j < _n_rings - 1; j++)
98 117 : chancount += j * 6;
99 : // Adding external channels to the total count
100 54 : _n_channels = chancount + _nrods - 1 + (_n_rings - 1) * 6 + 6;
101 :
102 : // Utils for building the mesh
103 54 : _chan_to_pin_map.resize(_n_channels);
104 54 : _subch_type.resize(_n_channels);
105 54 : _subchannel_position.resize(_n_channels);
106 :
107 2646 : for (unsigned int i = 0; i < _n_channels; i++)
108 : {
109 2592 : _subchannel_position[i].reserve(3);
110 10368 : for (unsigned int j = 0; j < 3; j++)
111 : {
112 7776 : _subchannel_position.at(i).push_back(0.0);
113 : }
114 : }
115 :
116 : // create the subchannels
117 : k = 0; // initialize the subchannel counter index
118 171 : for (unsigned int i = 1; i < _n_rings; i++)
119 : {
120 : // find the closest Pin at back ring
121 1251 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
122 : {
123 1134 : if (j == _pins_in_rings[i].size() - 1)
124 : {
125 117 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
126 117 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
127 117 : avg_coor_x =
128 117 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
129 117 : avg_coor_y =
130 117 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
131 : }
132 : else
133 : {
134 1017 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
135 1017 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
136 1017 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
137 1017 : _pin_position[_pins_in_rings[i][j + 1]](0));
138 1017 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
139 1017 : _pin_position[_pins_in_rings[i][j + 1]](1));
140 : }
141 : dist0 = 1.0e+5;
142 1134 : _chan_to_pin_map[k].push_back(_pins_in_rings[i - 1][0]);
143 7290 : for (unsigned int l = 0; l < _pins_in_rings[i - 1].size(); l++)
144 : {
145 6156 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i - 1][l]](0) - avg_coor_x, 2) +
146 6156 : pow(_pin_position[_pins_in_rings[i - 1][l]](1) - avg_coor_y, 2));
147 :
148 6156 : if (dist < dist0)
149 : {
150 2736 : _chan_to_pin_map[k][2] = _pins_in_rings[i - 1][l];
151 : dist0 = dist;
152 : }
153 : }
154 1134 : _subch_type[k] = EChannelType::CENTER;
155 1134 : _orientation_map.insert(std::make_pair(k, 0.0));
156 1134 : k = k + 1;
157 : }
158 :
159 : // find the closest Pin at front ring
160 1251 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
161 : {
162 1134 : if (j == _pins_in_rings[i].size() - 1)
163 : {
164 117 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
165 117 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
166 117 : avg_coor_x =
167 117 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
168 117 : avg_coor_y =
169 117 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
170 : }
171 : else
172 : {
173 1017 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
174 1017 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
175 1017 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
176 1017 : _pin_position[_pins_in_rings[i][j + 1]](0));
177 1017 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
178 1017 : _pin_position[_pins_in_rings[i][j + 1]](1));
179 : }
180 : // if the outermost ring, set the edge subchannels first... then the corner subchannels
181 1134 : if (i == _n_rings - 1)
182 : {
183 : // add edges
184 702 : _subch_type[k] = EChannelType::EDGE; // an edge subchannel is created
185 702 : k = k + 1;
186 702 : if (j % i == 0)
187 : {
188 : // corner subchannel
189 324 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
190 324 : _subch_type[k] = EChannelType::CORNER;
191 324 : k = k + 1;
192 : }
193 : // if not the outer most ring
194 : }
195 : else
196 : {
197 : dist0 = 1.0e+5;
198 432 : _chan_to_pin_map[k].push_back(_pins_in_rings[i + 1][0]);
199 6264 : for (unsigned int l = 0; l < _pins_in_rings[i + 1].size(); l++)
200 : {
201 5832 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i + 1][l]](0) - avg_coor_x, 2) +
202 5832 : pow(_pin_position[_pins_in_rings[i + 1][l]](1) - avg_coor_y, 2));
203 5832 : if (dist < dist0)
204 : {
205 2106 : _chan_to_pin_map[k][2] = _pins_in_rings[i + 1][l];
206 : dist0 = dist;
207 : }
208 : }
209 432 : _subch_type[k] = EChannelType::CENTER;
210 432 : _orientation_map.insert(std::make_pair(k, libMesh::pi));
211 432 : k = k + 1;
212 : }
213 : }
214 : }
215 :
216 2646 : for (auto & pin : _chan_to_pin_map)
217 : pin.shrink_to_fit();
218 :
219 : // set the subchannel positions
220 54 : Real _duct_to_pin_gap =
221 54 : 0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
222 2646 : for (unsigned int i = 0; i < _n_channels; i++)
223 : {
224 2592 : if (_subch_type[i] == EChannelType::CENTER)
225 : {
226 1566 : _subchannel_position[i][0] =
227 1566 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0) +
228 1566 : _pin_position[_chan_to_pin_map[i][2]](0)) /
229 : 3.0;
230 1566 : _subchannel_position[i][1] =
231 1566 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1) +
232 1566 : _pin_position[_chan_to_pin_map[i][2]](1)) /
233 : 3.0;
234 : }
235 1026 : else if (_subch_type[i] == EChannelType::EDGE)
236 : {
237 36018 : for (unsigned int j = 0; j < _n_channels; j++)
238 : {
239 35316 : if (_subch_type[j] == EChannelType::CENTER &&
240 21708 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
241 702 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1]) ||
242 21006 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
243 702 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
244 : {
245 702 : x0 = _pin_position[_chan_to_pin_map[j][2]](0);
246 702 : y0 = _pin_position[_chan_to_pin_map[j][2]](1);
247 : }
248 34614 : else if (_subch_type[j] == EChannelType::CENTER &&
249 21006 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
250 0 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
251 21006 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
252 378 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
253 : {
254 0 : x0 = _pin_position[_chan_to_pin_map[j][1]](0);
255 0 : y0 = _pin_position[_chan_to_pin_map[j][1]](1);
256 : }
257 34614 : else if (_subch_type[j] == EChannelType::CENTER &&
258 21006 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
259 702 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
260 21006 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
261 378 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1])))
262 : {
263 0 : x0 = _pin_position[_chan_to_pin_map[j][0]](0);
264 0 : y0 = _pin_position[_chan_to_pin_map[j][0]](1);
265 : }
266 35316 : x1 = 0.5 *
267 35316 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0));
268 35316 : y1 = 0.5 *
269 35316 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1));
270 35316 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
271 35316 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
272 35316 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
273 35316 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
274 : } // j
275 : }
276 324 : else if (_subch_type[i] == EChannelType::CORNER)
277 : {
278 324 : x0 = _pin_position[0](0);
279 324 : y0 = _pin_position[0](1);
280 324 : x1 = _pin_position[_chan_to_pin_map[i][0]](0);
281 324 : y1 = _pin_position[_chan_to_pin_map[i][0]](1);
282 324 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
283 324 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
284 324 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
285 324 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
286 : }
287 :
288 : /// Special case _n_rings == 1
289 2592 : if (_n_rings == 1)
290 : {
291 0 : for (unsigned int i = 0; i < _n_channels; i++)
292 : {
293 0 : Real angle = (2 * i + 1) * libMesh::pi / 6.0;
294 0 : _subch_type[i] = EChannelType::CORNER;
295 0 : _subchannel_position[i][0] = std::cos(angle) * _flat_to_flat / 2.0;
296 0 : _subchannel_position[i][1] = std::sin(angle) * _flat_to_flat / 2.0;
297 : }
298 : }
299 : }
300 54 : }
301 :
302 : std::unique_ptr<MeshBase>
303 54 : SCMDetailedTriSubChannelMeshGenerator::generate()
304 : {
305 54 : auto mesh_base = buildMeshBaseObject();
306 : BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
307 54 : mesh_base->set_spatial_dimension(3);
308 : // Define the resolution (the number of points used to represent a circle).
309 : // This must be divisible by 4.
310 : const unsigned int theta_res_triangle = 24; // TODO: parameterize
311 : const unsigned int theta_res_square = 24; // TODO: parameterize
312 : // Compute the number of points needed to represent one sixth and quadrant of a circle.
313 : const unsigned int points_per_sixth = theta_res_triangle / 6 + 1;
314 : const unsigned int points_per_quadrant = theta_res_square / 4 + 1;
315 :
316 : // Compute the points needed to represent one axial cross-flow of a subchannel.
317 : // For the center subchannel (sc) there is one center point plus the points from 3 side
318 : // circles.
319 54 : const unsigned int points_per_center = points_per_sixth * 3 + 1;
320 : // For the corner sc there is one center point plus the points from 1 side circle plus 3
321 : // corners
322 54 : const unsigned int points_per_corner = points_per_sixth * 1 + 1 + 3;
323 : // For the side sc there is one center point plus the points from 2 intersecting circles plus 2
324 : // corners
325 54 : const unsigned int points_per_side = points_per_quadrant * 2 + 1 + 2;
326 :
327 54 : if (_verbose)
328 : {
329 9 : _console << "Points per center: " << points_per_center << std::endl;
330 9 : _console << "Points per side: " << points_per_side << std::endl;
331 9 : _console << "Points per corner: " << points_per_corner << std::endl;
332 : }
333 :
334 : // Compute the number of elements (Prism6) which combined base creates the sub-channel
335 : // cross-section
336 54 : const unsigned int elems_per_center = theta_res_triangle * 3 / 6 + 3; // TODO: check
337 54 : const unsigned int elems_per_corner = theta_res_triangle / 6 + 4;
338 54 : const unsigned int elems_per_side = 2 * theta_res_square / 4 + 4;
339 54 : if (_verbose)
340 : {
341 9 : _console << "Elems per center: " << elems_per_center << std::endl;
342 9 : _console << "Elems per side: " << elems_per_side << std::endl;
343 9 : _console << "Elems per corner: " << elems_per_corner << std::endl;
344 : }
345 :
346 : // specify number and type of sub-channel
347 : unsigned int n_center, n_side, n_corner;
348 54 : if (_n_rings == 1)
349 : {
350 0 : n_corner = 6;
351 0 : n_side = 0;
352 0 : n_center = _n_channels - n_side - n_corner;
353 : }
354 : else
355 : {
356 54 : n_corner = 6;
357 54 : n_side = (_n_rings - 1) * 6;
358 54 : n_center = _n_channels - n_side - n_corner;
359 : }
360 54 : if (_verbose)
361 : {
362 9 : _console << "Centers: " << n_center << std::endl;
363 9 : _console << "Sides: " << n_side << std::endl;
364 9 : _console << "Corners: " << n_corner << std::endl;
365 : }
366 :
367 : // Compute the total number of points and elements.
368 54 : const unsigned int points_per_level =
369 54 : n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
370 54 : const unsigned int elems_per_level =
371 54 : n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
372 54 : if (_verbose)
373 : {
374 9 : _console << "Points per level: " << points_per_level << std::endl;
375 9 : _console << "Elements per level: " << elems_per_level << std::endl;
376 : }
377 54 : const unsigned int n_points = points_per_level * (_n_cells + 1);
378 54 : const unsigned int n_elems = elems_per_level * _n_cells;
379 54 : if (_verbose)
380 : {
381 9 : _console << "Number of points: " << n_points << std::endl;
382 9 : _console << "Number of elements: " << n_elems << std::endl;
383 : }
384 54 : mesh_base->reserve_nodes(n_points);
385 54 : mesh_base->reserve_elem(n_elems);
386 : // Build an array of points arranged in a circle on the xy-plane. (last and first node overlap)
387 : // We build for both the square discretization in the edges and the triangular discretization
388 : // within the mesh
389 54 : const Real radius = _pin_diameter / 2.0;
390 : std::array<Point, theta_res_square + 1> circle_points_square;
391 : {
392 : Real theta = 0;
393 1404 : for (unsigned int i = 0; i < theta_res_square + 1; i++)
394 : {
395 1350 : circle_points_square[i](0) = radius * std::cos(theta);
396 1350 : circle_points_square[i](1) = radius * std::sin(theta);
397 1350 : theta += 2.0 * libMesh::pi / theta_res_square;
398 : }
399 : }
400 : std::array<Point, theta_res_triangle + 1> circle_points_triangle;
401 : {
402 : Real theta = 0;
403 1404 : for (unsigned int i = 0; i < theta_res_triangle + 1; i++)
404 : {
405 1350 : circle_points_triangle[i](0) = radius * std::cos(theta);
406 1350 : circle_points_triangle[i](1) = radius * std::sin(theta);
407 1350 : theta += 2.0 * libMesh::pi / theta_res_triangle;
408 : }
409 : }
410 : // Define "quadrant center" reference points. These will be the centers of
411 : // the 3 circles that represent the fuel pins. These centers are
412 : // offset a little bit so that in the final mesh, there is a tiny gap between
413 : // neighboring subchannel cells. That allows us to easily map a solution to
414 : // this detailed mesh with a nearest-neighbor search.
415 : const Real shrink_factor = 0.99999;
416 : // Quadrants are used only for the side and corner subchannels
417 54 : Real _duct_to_pin_gap =
418 54 : 0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
419 : std::array<Point, 2> quadrant_centers_sides;
420 54 : quadrant_centers_sides[0] = Point(
421 54 : -_pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
422 54 : quadrant_centers_sides[1] = Point(
423 : _pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
424 : std::array<Point, 1> quadrant_centers_corner;
425 54 : quadrant_centers_corner[0] =
426 54 : Point(-(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::sin(libMesh::pi / 6) * shrink_factor,
427 54 : -(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::cos(libMesh::pi / 6) * shrink_factor,
428 : 0);
429 : // Triangles are used for all center subchannels
430 : std::array<Point, 3> triangle_centers;
431 54 : triangle_centers[0] = Point(0, _pitch * std::cos(libMesh::pi / 6) * 2 / 3 * shrink_factor, 0);
432 54 : triangle_centers[1] = Point(-_pitch * 0.5 * shrink_factor,
433 54 : -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor,
434 : 0);
435 54 : triangle_centers[2] = Point(
436 : _pitch * 0.5 * shrink_factor, -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor, 0);
437 :
438 : const unsigned int m_sixth = theta_res_triangle / 6;
439 : const unsigned int m_quarter = theta_res_square / 4;
440 : // Build an array of points that represent a cross section of a center subchannel
441 : // cell. The points are ordered in this fashion:
442 : // 3 1
443 : // 2
444 : // 0
445 : // 4 5 8 9
446 : // 6 7
447 : std::array<Point, points_per_center> center_points;
448 : {
449 : unsigned int start;
450 216 : for (unsigned int i = 0; i < 3; i++)
451 : {
452 162 : if (i == 0)
453 : start = 5 * (m_sixth);
454 108 : if (i == 1)
455 : start = 1 * (m_sixth);
456 108 : if (i == 2)
457 : start = 3 * (m_sixth);
458 972 : for (unsigned int ii = 0; ii < points_per_sixth; ii++)
459 : {
460 810 : auto c_pt = circle_points_triangle[start - ii];
461 810 : center_points[i * points_per_sixth + ii + 1] = triangle_centers[i] + c_pt;
462 : }
463 : }
464 : }
465 :
466 : // Build an array of points that represent a cross section of a top left corner subchannel
467 : // cell. The points are ordered in this fashion (x direction towards point 4, y direction towards
468 : // point 6):
469 : // 5 6
470 : //
471 : // 0
472 : // 4 2 3
473 : // 1
474 : std::array<Point, points_per_corner> corner_points;
475 : {
476 324 : for (unsigned int ii = 0; ii < points_per_sixth; ii++)
477 : {
478 270 : auto c_pt = circle_points_triangle[1 * m_quarter - ii];
479 270 : corner_points[ii + 1] = quadrant_centers_corner[0] + c_pt;
480 : }
481 : Real side_short = (_duct_to_pin_gap + _pin_diameter) * 0.5;
482 54 : Real side_long = (2.0 * _duct_to_pin_gap + _pin_diameter) * 0.5;
483 54 : Real side_length = std::sqrt(std::pow(side_short, 2) + std::pow(side_long, 2) -
484 54 : 2 * side_short * side_long * std::cos(libMesh::pi / 6));
485 : Real angle =
486 : libMesh::pi - libMesh::pi / 3 -
487 54 : std::acos((-std::pow(side_long, 2) + std::pow(side_short, 2) + std::pow(side_length, 2)) /
488 54 : (2 * side_short * side_length));
489 54 : corner_points[points_per_sixth + 1] = Point(side_length * std::cos(angle) * shrink_factor,
490 54 : -side_length * std::sin(angle) * shrink_factor,
491 : 0);
492 54 : corner_points[points_per_sixth + 2] =
493 54 : Point(side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor *
494 : std::tan(libMesh::pi / 6),
495 54 : side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
496 : 0);
497 54 : corner_points[points_per_sixth + 3] =
498 54 : Point(-side_length * std::cos(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
499 : side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
500 : 0);
501 : }
502 :
503 : // Build an array of points that represent a cross-section of a top side subchannel
504 : // cell. The points are ordered in this fashion:
505 : // 8 7
506 : //
507 : // 0
508 : // 1 2 5 6
509 : // 3 4
510 : std::array<Point, points_per_side> side_points;
511 : {
512 432 : for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
513 : {
514 378 : auto c_pt = circle_points_square[m_quarter - ii];
515 378 : side_points[ii + 1] = quadrant_centers_sides[0] + c_pt;
516 : }
517 432 : for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
518 : {
519 378 : auto c_pt = circle_points_square[2 * m_quarter - ii];
520 378 : side_points[points_per_quadrant + ii + 1] = quadrant_centers_sides[1] + c_pt;
521 : }
522 54 : side_points[2 * points_per_quadrant + 1] =
523 54 : Point(_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
524 54 : side_points[2 * points_per_quadrant + 2] =
525 : Point(-_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
526 : }
527 :
528 54 : int point_counter = 0;
529 : unsigned int node_id = 0;
530 2646 : for (unsigned int i = 0; i < _n_channels; i++)
531 : {
532 : // Real offset_x = _flat_to_flat / 2.0;
533 : // Real offset_y = _flat_to_flat / 2.0;
534 :
535 2592 : if (getSubchannelType(i) == EChannelType::CENTER)
536 : {
537 51732 : for (auto z : _z_grid)
538 : {
539 : // Get height
540 : Point z0{0, 0, z};
541 :
542 : // Get suchannel position and assign to point
543 50166 : auto loc_position = getSubchannelPosition(i);
544 50166 : Point p0{loc_position[0], loc_position[1], 0};
545 :
546 : // Determine orientation of current subchannel
547 50166 : auto subchannel_pins = getSubChannelPins(i);
548 : Point subchannel_side =
549 50166 : getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
550 : Point base_center_orientation = {0, -1};
551 :
552 : // Get rotation angle for current subchannel
553 : Real dot_prod = 0;
554 150498 : for (unsigned int lp = 0; lp < 2; lp++)
555 100332 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
556 : auto theta =
557 50166 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
558 50166 : if (subchannel_side(0) < 0)
559 25083 : theta = 2.0 * libMesh::pi - theta;
560 :
561 : // Real distance_side = subchannel_side.norm();
562 : // Real distance_top = getPinPosition(subchannel_pins[2]).norm();
563 : // if (distance_top > distance_side)
564 : // theta += libMesh::pi * 0.0;
565 :
566 50166 : theta += _orientation_map[i];
567 :
568 50166 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
569 :
570 50166 : if (_verbose)
571 : {
572 4536 : if (z == 0)
573 : {
574 216 : _console << "Subchannel Position: " << p0 << std::endl;
575 216 : auto pins = getSubChannelPins(i);
576 864 : for (auto r : pins)
577 648 : _console << r << " ";
578 216 : _console << std::endl;
579 216 : _console << "Theta: " << theta / libMesh::pi * 180. << std::endl;
580 216 : }
581 : }
582 :
583 : // Assigning points for center channels
584 852822 : for (unsigned int i = 0; i < points_per_center; i++)
585 : {
586 802656 : auto new_point = rotatePoint(center_points[i], theta) + p0 + z0;
587 802656 : if (_verbose)
588 : {
589 72576 : if (z == 0)
590 3456 : _console << i << " - " << new_point << std::endl;
591 : }
592 802656 : mesh_base->add_point(new_point, node_id++);
593 802656 : point_counter += 1;
594 : }
595 50166 : }
596 : }
597 1026 : else if (getSubchannelType(i) == EChannelType::EDGE)
598 : {
599 24084 : for (auto z : _z_grid)
600 : {
601 : // Get height
602 : Point z0{0, 0, z};
603 :
604 : // Get suchannel position and assign to point
605 23382 : auto loc_position = getSubchannelPosition(i);
606 23382 : Point p0{loc_position[0], loc_position[1], 0};
607 :
608 : // Determine orientation of current subchannel
609 23382 : auto subchannel_pins = getSubChannelPins(i);
610 : Point subchannel_side =
611 23382 : getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
612 : Point base_center_orientation = {0, 1};
613 :
614 : // Get rotation angle for current subchannel
615 : Real dot_prod = 0;
616 70146 : for (unsigned int lp = 0; lp < 2; lp++)
617 46764 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
618 : auto theta =
619 23382 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
620 23382 : if (subchannel_side(0) > 0)
621 11691 : theta = 2. * libMesh::pi - theta;
622 23382 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
623 :
624 23382 : if (_verbose)
625 : {
626 2268 : if (z == 0)
627 : {
628 108 : _console << "Subchannel Position: " << p0 << std::endl;
629 108 : auto pins = getSubChannelPins(i);
630 324 : for (auto r : pins)
631 216 : _console << r << " ";
632 108 : _console << std::endl;
633 108 : _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
634 108 : }
635 : }
636 :
637 : // Assigning points for center channels
638 420876 : for (unsigned int i = 0; i < points_per_side; i++)
639 : {
640 397494 : auto new_point = rotatePoint(side_points[i], theta) + p0 + z0;
641 397494 : if (_verbose)
642 : {
643 38556 : if (z == 0)
644 1836 : _console << i << " - " << new_point << std::endl;
645 : }
646 397494 : mesh_base->add_point(new_point, node_id++);
647 397494 : point_counter += 1;
648 : }
649 23382 : }
650 : }
651 : else // getSubchannelType(i) == EChannelType::CORNER
652 : {
653 11448 : for (auto z : _z_grid)
654 : {
655 : // Get height
656 : Point z0{0, 0, z};
657 :
658 : // Get suchannel position and assign to point
659 11124 : auto loc_position = getSubchannelPosition(i);
660 11124 : Point p0{loc_position[0], loc_position[1], 0};
661 :
662 : // Determine orientation of current subchannel
663 11124 : auto subchannel_pins = getSubChannelPins(i);
664 11124 : Point subchannel_side = getPinPosition(subchannel_pins[0]);
665 : Point base_center_orientation = {1, 1};
666 :
667 : // Get rotation angle for current subchannel
668 : Real dot_prod = 0;
669 33372 : for (unsigned int lp = 0; lp < 2; lp++)
670 22248 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
671 : auto theta =
672 11124 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
673 11124 : if (subchannel_side(0) > 0)
674 5562 : theta = 2. * libMesh::pi - theta;
675 11124 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
676 :
677 11124 : if (_verbose)
678 : {
679 1134 : if (z == 0)
680 : {
681 54 : _console << "Subchannel Position: " << p0 << std::endl;
682 54 : auto pins = getSubChannelPins(i);
683 108 : for (auto r : pins)
684 54 : _console << r << " ";
685 54 : _console << std::endl;
686 54 : _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
687 54 : }
688 : }
689 :
690 : // Assigning points for center channels
691 111240 : for (unsigned int i = 0; i < points_per_corner; i++)
692 : {
693 100116 : auto new_point = rotatePoint(corner_points[i], theta) + p0 + z0;
694 100116 : if (_verbose)
695 : {
696 10206 : if (z == 0)
697 486 : _console << i << " - " << new_point << std::endl;
698 : }
699 100116 : mesh_base->add_point(new_point, node_id++);
700 100116 : point_counter += 1;
701 : }
702 11124 : }
703 : }
704 : } // i
705 54 : if (_verbose)
706 9 : _console << "Point counter: " << point_counter << std::endl;
707 :
708 54 : int element_counter = 0;
709 : unsigned int elem_id = 0;
710 : unsigned int number_of_corner = 0;
711 : unsigned int number_of_side = 0;
712 : unsigned int number_of_center = 0;
713 : unsigned int elems_per_channel = 0;
714 : unsigned int points_per_channel = 0;
715 2646 : for (unsigned int i = 0; i < _n_channels; i++)
716 : {
717 : auto subch_type = getSubchannelType(i);
718 2592 : if (subch_type == EChannelType::CORNER)
719 : {
720 324 : number_of_corner++;
721 : elems_per_channel = elems_per_corner;
722 : points_per_channel = points_per_corner;
723 324 : if (_verbose)
724 54 : _console << "Corner" << std::endl;
725 : }
726 2268 : else if (subch_type == EChannelType::EDGE)
727 : {
728 702 : number_of_side++;
729 : elems_per_channel = elems_per_side;
730 : points_per_channel = points_per_side;
731 702 : if (_verbose)
732 108 : _console << "Edge" << std::endl;
733 : }
734 1566 : else if (subch_type == EChannelType::CENTER)
735 : {
736 1566 : number_of_center++;
737 : elems_per_channel = elems_per_center;
738 : points_per_channel = points_per_center;
739 1566 : if (_verbose)
740 216 : _console << "Center" << std::endl;
741 : }
742 84672 : for (unsigned int iz = 0; iz < _n_cells; iz++)
743 : {
744 82080 : unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
745 82080 : number_of_side * points_per_side * (_n_cells + 1) +
746 82080 : number_of_center * points_per_center * (_n_cells + 1) -
747 82080 : points_per_channel * (_n_cells + 1);
748 : // index of the central node at base of cell
749 82080 : unsigned int indx1 = iz * points_per_channel + elapsed_points;
750 : // index of the central node at top of cell
751 82080 : unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
752 :
753 1260360 : for (unsigned int i = 0; i < elems_per_channel; i++)
754 : {
755 1178280 : Elem * elem = new Prism6;
756 1178280 : elem->subdomain_id() = _block_id;
757 1178280 : elem->set_id(elem_id++);
758 1178280 : elem = mesh_base->add_elem(elem);
759 :
760 1178280 : if (_verbose)
761 108000 : _console << "Node 0: " << *mesh_base->node_ptr(indx1) << std::endl;
762 :
763 1178280 : elem->set_node(0, mesh_base->node_ptr(indx1));
764 1178280 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
765 1178280 : if (i != elems_per_channel - 1)
766 1096200 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
767 : else
768 82080 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
769 :
770 1178280 : elem->set_node(3, mesh_base->node_ptr(indx2));
771 1178280 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
772 1178280 : if (i != elems_per_channel - 1)
773 1096200 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
774 : else
775 82080 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
776 :
777 1178280 : if (iz == 0)
778 37314 : boundary_info.add_side(elem, 0, 0);
779 1178280 : if (iz == _n_cells - 1)
780 37314 : boundary_info.add_side(elem, 4, 1);
781 :
782 1178280 : element_counter += 1;
783 : }
784 : }
785 : }
786 54 : if (_verbose)
787 9 : _console << "Element counter: " << element_counter << std::endl;
788 54 : boundary_info.sideset_name(0) = "inlet";
789 54 : boundary_info.sideset_name(1) = "outlet";
790 54 : mesh_base->subdomain_name(_block_id) = name();
791 54 : if (_verbose)
792 9 : _console << "Mesh assembly done" << std::endl;
793 54 : mesh_base->prepare_for_use();
794 :
795 54 : return mesh_base;
796 0 : }
797 :
798 : Point
799 1300266 : SCMDetailedTriSubChannelMeshGenerator::rotatePoint(Point b, Real theta)
800 : {
801 :
802 : // Building rotation matrix
803 : std::vector<std::vector<Real>> A;
804 1300266 : A.resize(3);
805 5201064 : for (std::vector<Real> a : A)
806 : {
807 3900798 : a.resize(3);
808 3900798 : }
809 :
810 1300266 : A[0] = {std::cos(theta), -std::sin(theta), 0.0};
811 1300266 : A[1] = {std::sin(theta), std::cos(theta), 0.0};
812 1300266 : A[2] = {0.0, 0.0, 1.0};
813 :
814 : // Rotating vector
815 : Point rotated_vector = Point(0.0, 0.0, 0.0);
816 5201064 : for (unsigned int i = 0; i < 3; i++)
817 : {
818 15603192 : for (unsigned int j = 0; j < 3; j++)
819 : {
820 11702394 : rotated_vector(i) += A[i][j] * b(j);
821 : }
822 : }
823 :
824 1300266 : return rotated_vector;
825 1300266 : }
826 :
827 : Point
828 0 : SCMDetailedTriSubChannelMeshGenerator::translatePoint(Point & b, Point & translation_vector)
829 : {
830 : // Translating point
831 : Point translated_vector = Point(0.0, 0.0, 0.0);
832 0 : for (unsigned int i = 0; i < 3; i++)
833 : {
834 0 : translated_vector(i) = b(i) + translation_vector(i);
835 : }
836 :
837 0 : return translated_vector;
838 : }
|
