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