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 "SCMDetailedQuadSubChannelMeshGenerator.h"
11 : #include <array>
12 : #include <cmath>
13 : #include "libmesh/cell_prism6.h"
14 : #include "libmesh/unstructured_mesh.h"
15 :
16 : registerMooseObject("SubChannelApp", SCMDetailedQuadSubChannelMeshGenerator);
17 :
18 : InputParameters
19 156 : SCMDetailedQuadSubChannelMeshGenerator::validParams()
20 : {
21 156 : InputParameters params = MeshGenerator::validParams();
22 156 : params.addClassDescription(
23 : "Creates a detailed mesh of subchannels in a square lattice arrangement");
24 312 : params.addRequiredParam<Real>("pitch", "Pitch [m]");
25 312 : params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
26 312 : params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
27 312 : params.addRequiredParam<Real>("heated_length", "Heated length [m]");
28 312 : params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
29 312 : params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
30 312 : params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
31 312 : params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");
32 312 : params.addRequiredParam<Real>(
33 : "gap",
34 : "The side gap, not to be confused with the gap between pins, this refers to the gap "
35 : "next to the duct or else the distance between the subchannel centroid to the duct wall."
36 : "distance(edge pin center, duct wall) = pitch / 2 + side_gap [m]");
37 312 : params.deprecateParam("gap", "side_gap", "08/06/2026");
38 312 : params.addParam<unsigned int>("block_id", 0, "Block ID used for the mesh subdomain.");
39 156 : return params;
40 0 : }
41 :
42 78 : SCMDetailedQuadSubChannelMeshGenerator::SCMDetailedQuadSubChannelMeshGenerator(
43 78 : const InputParameters & parameters)
44 : : MeshGenerator(parameters),
45 78 : _unheated_length_entry(getParam<Real>("unheated_length_entry")),
46 156 : _heated_length(getParam<Real>("heated_length")),
47 156 : _unheated_length_exit(getParam<Real>("unheated_length_exit")),
48 156 : _pitch(getParam<Real>("pitch")),
49 156 : _pin_diameter(getParam<Real>("pin_diameter")),
50 156 : _n_cells(getParam<unsigned int>("n_cells")),
51 156 : _nx(getParam<unsigned int>("nx")),
52 156 : _ny(getParam<unsigned int>("ny")),
53 78 : _n_channels(_nx * _ny),
54 156 : _side_gap(getParam<Real>("side_gap")),
55 234 : _block_id(getParam<unsigned int>("block_id"))
56 : {
57 78 : Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
58 78 : Real dz = L / _n_cells;
59 2241 : for (unsigned int i = 0; i < _n_cells + 1; i++)
60 2163 : _z_grid.push_back(dz * i);
61 :
62 78 : _subchannel_position.resize(_n_channels);
63 1080 : for (unsigned int i = 0; i < _n_channels; i++)
64 : {
65 1002 : _subchannel_position[i].reserve(3);
66 4008 : for (unsigned int j = 0; j < 3; j++)
67 : {
68 3006 : _subchannel_position.at(i).push_back(0.0);
69 : }
70 : }
71 :
72 78 : _subch_type.resize(_n_channels);
73 314 : for (unsigned int iy = 0; iy < _ny; iy++)
74 : {
75 1238 : for (unsigned int ix = 0; ix < _nx; ix++)
76 : {
77 1002 : unsigned int i_ch = _nx * iy + ix;
78 924 : bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
79 1862 : (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
80 1002 : bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);
81 :
82 1002 : if (is_corner)
83 258 : _subch_type[i_ch] = EChannelType::CORNER;
84 744 : else if (is_edge)
85 454 : _subch_type[i_ch] = EChannelType::EDGE;
86 : else
87 290 : _subch_type[i_ch] = EChannelType::CENTER;
88 :
89 : // set the subchannel positions so that the center of the assembly is the zero point
90 1002 : Real offset_x = (_nx - 1) * _pitch / 2.0;
91 1002 : Real offset_y = (_ny - 1) * _pitch / 2.0;
92 1002 : _subchannel_position[i_ch][0] = _pitch * ix - offset_x;
93 1002 : _subchannel_position[i_ch][1] = _pitch * iy - offset_y;
94 : }
95 : }
96 78 : }
97 :
98 : std::unique_ptr<MeshBase>
99 78 : SCMDetailedQuadSubChannelMeshGenerator::generate()
100 : {
101 78 : auto mesh_base = buildMeshBaseObject();
102 : BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
103 78 : mesh_base->set_spatial_dimension(3);
104 : // Define the resolution (the number of points used to represent a circle).
105 : // This must be divisible by 4.
106 : const unsigned int theta_res = 16; // TODO: parameterize
107 : // Compute the number of points needed to represent one quarter of a circle.
108 : const unsigned int points_per_quad = theta_res / 4 + 1;
109 :
110 : // Compute the points needed to represent one axial cross-flow of a subchannel.
111 : // For the center subchannel (sc) there is one center point plus the points from 4 intersecting
112 : // circles.
113 : const unsigned int points_per_center = points_per_quad * 4 + 1;
114 : // For the corner sc there is one center point plus the points from 1 intersecting circle plus 3
115 : // corners
116 : const unsigned int points_per_corner = points_per_quad * 1 + 1 + 3;
117 : // For the side sc there is one center point plus the points from 2 intersecting circles plus 2
118 : // corners
119 : const unsigned int points_per_side = points_per_quad * 2 + 1 + 2;
120 :
121 : // Compute the number of elements (Prism6) which combined base creates the sub-channel
122 : // cross-section
123 : const unsigned int elems_per_center = theta_res + 4;
124 : const unsigned int elems_per_corner = theta_res / 4 + 4;
125 : const unsigned int elems_per_side = theta_res / 2 + 4;
126 :
127 : // specify number and type of sub-channel
128 : unsigned int n_center, n_side, n_corner;
129 78 : if (_n_channels == 2)
130 : {
131 : n_corner = 0;
132 : n_side = _n_channels;
133 : n_center = _n_channels - n_side - n_corner;
134 : }
135 65 : else if (_n_channels > 2 && (_ny == 1 || _nx == 1))
136 : {
137 : n_corner = 0;
138 : n_side = 2;
139 14 : n_center = _n_channels - n_side - n_corner;
140 : }
141 : else
142 : {
143 : n_corner = 4;
144 51 : n_side = 2 * (_ny - 2) + 2 * (_nx - 2);
145 51 : n_center = _n_channels - n_side - n_corner;
146 : }
147 :
148 : // Compute the total number of points and elements.
149 78 : const unsigned int points_per_level =
150 78 : n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
151 78 : const unsigned int elems_per_level =
152 78 : n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
153 78 : const unsigned int n_points = points_per_level * (_n_cells + 1);
154 78 : const unsigned int n_elems = elems_per_level * _n_cells;
155 78 : mesh_base->reserve_nodes(n_points);
156 78 : mesh_base->reserve_elem(n_elems);
157 : // Build an array of points arranged in a circle on the xy-plane. (last and first node overlap)
158 78 : const Real radius = _pin_diameter / 2.0;
159 : std::array<Point, theta_res + 1> circle_points;
160 : {
161 : Real theta = 0;
162 1404 : for (unsigned int i = 0; i < theta_res + 1; i++)
163 : {
164 1326 : circle_points[i](0) = radius * std::cos(theta);
165 1326 : circle_points[i](1) = radius * std::sin(theta);
166 1326 : theta += 2 * M_PI / theta_res;
167 : }
168 : }
169 : // Define "quadrant center" reference points. These will be the centers of
170 : // the 4 circles that represent the fuel pins. These centers are
171 : // offset a little bit so that in the final mesh, there is a tiny gap between
172 : // neighboring subchannel cells. That allows us to easily map a solution to
173 : // this detailed mesh with a nearest-neighbor search.
174 : const Real shrink_factor = 0.99999;
175 : std::array<Point, 4> quadrant_centers;
176 78 : quadrant_centers[0] = Point(_pitch * 0.5 * shrink_factor, _pitch * 0.5 * shrink_factor, 0);
177 78 : quadrant_centers[1] = Point(-_pitch * 0.5 * shrink_factor, _pitch * 0.5 * shrink_factor, 0);
178 78 : quadrant_centers[2] = Point(-_pitch * 0.5 * shrink_factor, -_pitch * 0.5 * shrink_factor, 0);
179 78 : quadrant_centers[3] = Point(_pitch * 0.5 * shrink_factor, -_pitch * 0.5 * shrink_factor, 0);
180 :
181 : const unsigned int m = theta_res / 4;
182 : // Build an array of points that represent a cross section of a center subchannel
183 : // cell. The points are ordered in this fashion:
184 : // 4 3
185 : // 6 5 2 1
186 : // 0
187 : // 7 8 * *
188 : // 9 *
189 : std::array<Point, points_per_center> center_points;
190 : {
191 : unsigned int start;
192 390 : for (unsigned int i = 0; i < 4; i++)
193 : {
194 312 : if (i == 0)
195 : start = 3 * m;
196 234 : if (i == 1)
197 : start = 4 * m;
198 234 : if (i == 2)
199 : start = 1 * m;
200 234 : if (i == 3)
201 : start = 2 * m;
202 1872 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
203 : {
204 1560 : auto c_pt = circle_points[start - ii];
205 1560 : center_points[i * points_per_quad + ii + 1] = quadrant_centers[i] + c_pt;
206 : }
207 : }
208 : }
209 :
210 : // Build an array of points that represent a cross section of a top left corner subchannel
211 : // cell. The points are ordered in this fashion:
212 : // 5 4
213 : //
214 : // 0
215 : // 2 3
216 : // 6 1
217 : std::array<Point, points_per_corner> tl_corner_points;
218 : {
219 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
220 : {
221 390 : auto c_pt = circle_points[2 * m - ii];
222 390 : tl_corner_points[ii + 1] = quadrant_centers[3] + c_pt;
223 : }
224 78 : tl_corner_points[points_per_quad + 1] = Point(_pitch * 0.5 * shrink_factor, _side_gap, 0);
225 78 : tl_corner_points[points_per_quad + 2] = Point(-_side_gap, _side_gap, 0);
226 78 : tl_corner_points[points_per_quad + 3] = Point(-_side_gap, -_pitch * 0.5 * shrink_factor, 0);
227 : }
228 :
229 : // Build an array of points that represent a cross section of a top right corner subchannel
230 : // cell. The points are ordered in this fashion:
231 : // 6 5
232 : //
233 : // 0
234 : // 1 2
235 : // 3 4
236 : std::array<Point, points_per_corner> tr_corner_points;
237 : {
238 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
239 : {
240 390 : auto c_pt = circle_points[m - ii];
241 390 : tr_corner_points[ii + 1] = quadrant_centers[2] + c_pt;
242 : }
243 78 : tr_corner_points[points_per_quad + 1] = Point(_side_gap, -_pitch * 0.5 * shrink_factor, 0);
244 78 : tr_corner_points[points_per_quad + 2] = Point(_side_gap, _side_gap, 0);
245 78 : tr_corner_points[points_per_quad + 3] = Point(-_pitch * 0.5 * shrink_factor, _side_gap, 0);
246 : }
247 :
248 : // Build an array of points that represent a cross section of a bottom left corner subchannel
249 : // cell. The points are ordered in this fashion:
250 : // 4 3
251 : // 2 1
252 : // 0
253 : //
254 : // 5 6
255 : std::array<Point, points_per_corner> bl_corner_points;
256 : {
257 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
258 : {
259 390 : auto c_pt = circle_points[3 * m - ii];
260 390 : bl_corner_points[ii + 1] = quadrant_centers[0] + c_pt;
261 : }
262 78 : bl_corner_points[points_per_quad + 1] = Point(-_side_gap, _pitch * 0.5 * shrink_factor, 0);
263 78 : bl_corner_points[points_per_quad + 2] = Point(-_side_gap, -_side_gap, 0);
264 78 : bl_corner_points[points_per_quad + 3] = Point(_pitch * 0.5 * shrink_factor, -_side_gap, 0);
265 : }
266 :
267 : // Build an array of points that represent a cross section of a bottom right corner subchannel
268 : // cell. The points are ordered in this fashion:
269 : // 1 6
270 : // 3 2
271 : // 0
272 : //
273 : // 4 5
274 : std::array<Point, points_per_corner> br_corner_points;
275 : {
276 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
277 : {
278 390 : auto c_pt = circle_points[4 * m - ii];
279 390 : br_corner_points[ii + 1] = quadrant_centers[1] + c_pt;
280 : }
281 78 : br_corner_points[points_per_quad + 1] = Point(-_pitch * 0.5 * shrink_factor, -_side_gap, 0);
282 78 : br_corner_points[points_per_quad + 2] = Point(_side_gap, -_side_gap, 0);
283 78 : br_corner_points[points_per_quad + 3] = Point(_side_gap, _pitch * 0.5 * shrink_factor, 0);
284 : }
285 :
286 : // Build an array of points that represent a cross section of a top side subchannel
287 : // cell. The points are ordered in this fashion:
288 : // 8 7
289 : //
290 : // 0
291 : // 1 2 5 6
292 : // 3 4
293 : std::array<Point, points_per_side> top_points;
294 : {
295 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
296 : {
297 390 : auto c_pt = circle_points[m - ii];
298 390 : top_points[ii + 1] = quadrant_centers[2] + c_pt;
299 : }
300 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
301 : {
302 390 : auto c_pt = circle_points[2 * m - ii];
303 390 : top_points[points_per_quad + ii + 1] = quadrant_centers[3] + c_pt;
304 : }
305 78 : top_points[2 * points_per_quad + 1] = Point(_pitch * 0.5 * shrink_factor, _side_gap, 0);
306 78 : top_points[2 * points_per_quad + 2] = Point(-_pitch * 0.5 * shrink_factor, _side_gap, 0);
307 : }
308 :
309 : // Build an array of points that represent a cross section of a left side subchannel
310 : // cell. The points are ordered in this fashion:
311 : // 7 6
312 : // 5 4
313 : // 0
314 : // 2 3
315 : // 8 1
316 : std::array<Point, points_per_side> left_points;
317 : {
318 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
319 : {
320 390 : auto c_pt = circle_points[2 * m - ii];
321 390 : left_points[ii + 1] = quadrant_centers[3] + c_pt;
322 : }
323 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
324 : {
325 390 : auto c_pt = circle_points[3 * m - ii];
326 390 : left_points[points_per_quad + ii + 1] = quadrant_centers[0] + c_pt;
327 : }
328 78 : left_points[2 * points_per_quad + 1] = Point(-_side_gap, _pitch * 0.5 * shrink_factor, 0);
329 78 : left_points[2 * points_per_quad + 2] = Point(-_side_gap, -_pitch * 0.5 * shrink_factor, 0);
330 : }
331 :
332 : // Build an array of points that represent a cross section of a bottom side subchannel
333 : // cell. The points are ordered in this fashion:
334 : // 4 3
335 : // 6 5 2 1
336 : // 0
337 : //
338 : // 7 8
339 : std::array<Point, points_per_side> bottom_points;
340 : {
341 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
342 : {
343 390 : auto c_pt = circle_points[3 * m - ii];
344 390 : bottom_points[ii + 1] = quadrant_centers[0] + c_pt;
345 : }
346 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
347 : {
348 390 : auto c_pt = circle_points[4 * m - ii];
349 390 : bottom_points[points_per_quad + ii + 1] = quadrant_centers[1] + c_pt;
350 : }
351 78 : bottom_points[2 * points_per_quad + 1] = Point(-_pitch * 0.5 * shrink_factor, -_side_gap, 0);
352 78 : bottom_points[2 * points_per_quad + 2] = Point(_pitch * 0.5 * shrink_factor, -_side_gap, 0);
353 : }
354 :
355 : // Build an array of points that represent a cross section of a right side subchannel
356 : // cell. The points are ordered in this fashion:
357 : // 1 8
358 : // 3 2
359 : // 0
360 : // 4 5
361 : // 6 7
362 : std::array<Point, points_per_side> right_points;
363 : {
364 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
365 : {
366 390 : auto c_pt = circle_points[4 * m - ii];
367 390 : right_points[ii + 1] = quadrant_centers[1] + c_pt;
368 : }
369 468 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
370 : {
371 390 : auto c_pt = circle_points[m - ii];
372 390 : right_points[points_per_quad + ii + 1] = quadrant_centers[2] + c_pt;
373 : }
374 78 : right_points[2 * points_per_quad + 1] = Point(_side_gap, -_pitch * 0.5 * shrink_factor, 0);
375 78 : right_points[2 * points_per_quad + 2] = Point(_side_gap, _pitch * 0.5 * shrink_factor, 0);
376 : }
377 :
378 : // Add the points to the mesh.
379 78 : if (_n_channels == 2)
380 : {
381 : unsigned int node_id = 0;
382 13 : Real offset_x = (_nx - 1) * _pitch / 2.0;
383 13 : Real offset_y = (_ny - 1) * _pitch / 2.0;
384 33 : for (unsigned int iy = 0; iy < _ny; iy++)
385 : {
386 20 : Point y0 = {0, _pitch * iy - offset_y, 0};
387 46 : for (unsigned int ix = 0; ix < _nx; ix++)
388 : {
389 26 : Point x0 = {_pitch * ix - offset_x, 0, 0};
390 1392 : for (auto z : _z_grid)
391 : {
392 : Point z0{0, 0, z};
393 1366 : if (_nx == 1 && iy == 0) // vertical orientation
394 : {
395 1078 : for (unsigned int i = 0; i < points_per_side; i++)
396 1001 : mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
397 : }
398 1366 : if (_nx == 1 && iy == 1) // vertical orientation
399 : {
400 1078 : for (unsigned int i = 0; i < points_per_side; i++)
401 1001 : mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
402 : }
403 1366 : if (_ny == 1 && ix == 0) // horizontal orientation
404 : {
405 8484 : for (unsigned int i = 0; i < points_per_side; i++)
406 7878 : mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
407 : }
408 1366 : if (_ny == 1 && ix == 1) // horizontal orientation
409 : {
410 8484 : for (unsigned int i = 0; i < points_per_side; i++)
411 7878 : mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
412 : }
413 : }
414 : }
415 : }
416 : }
417 65 : else if (_n_channels > 2 && (_ny == 1 || _nx == 1))
418 : {
419 : unsigned int node_id = 0;
420 14 : Real offset_x = (_nx - 1) * _pitch / 2.0;
421 14 : Real offset_y = (_ny - 1) * _pitch / 2.0;
422 42 : for (unsigned int iy = 0; iy < _ny; iy++)
423 : {
424 28 : Point y0 = {0, _pitch * iy - offset_y, 0};
425 70 : for (unsigned int ix = 0; ix < _nx; ix++)
426 : {
427 42 : Point x0 = {_pitch * ix - offset_x, 0, 0};
428 504 : for (auto z : _z_grid)
429 : {
430 : Point z0{0, 0, z};
431 462 : if (_nx == 1)
432 : {
433 231 : if (iy == 0)
434 : {
435 1078 : for (unsigned int i = 0; i < points_per_side; i++)
436 1001 : mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
437 : }
438 154 : else if (iy == _ny - 1)
439 : {
440 1078 : for (unsigned int i = 0; i < points_per_side; i++)
441 1001 : mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
442 : }
443 : else
444 : {
445 1694 : for (unsigned int i = 0; i < points_per_center; i++)
446 1617 : mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
447 : }
448 : }
449 231 : else if (_ny == 1)
450 : {
451 231 : if (ix == 0)
452 : {
453 1078 : for (unsigned int i = 0; i < points_per_side; i++)
454 1001 : mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
455 : }
456 154 : else if (ix == _nx - 1)
457 : {
458 1078 : for (unsigned int i = 0; i < points_per_side; i++)
459 1001 : mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
460 : }
461 : else
462 : {
463 1694 : for (unsigned int i = 0; i < points_per_center; i++)
464 1617 : mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
465 : }
466 : }
467 : }
468 : }
469 : }
470 : }
471 : else
472 : {
473 : unsigned int node_id = 0;
474 51 : Real offset_x = (_nx - 1) * _pitch / 2.0;
475 51 : Real offset_y = (_ny - 1) * _pitch / 2.0;
476 239 : for (unsigned int iy = 0; iy < _ny; iy++)
477 : {
478 188 : Point y0 = {0, _pitch * iy - offset_y, 0};
479 1122 : for (unsigned int ix = 0; ix < _nx; ix++)
480 : {
481 934 : Point x0 = {_pitch * ix - offset_x, 0, 0};
482 934 : if (ix == 0 && iy == 0) // Bottom Left corner
483 : {
484 1377 : for (auto z : _z_grid)
485 : {
486 : Point z0{0, 0, z};
487 13260 : for (unsigned int i = 0; i < points_per_corner; i++)
488 11934 : mesh_base->add_point(bl_corner_points[i] + x0 + y0 + z0, node_id++);
489 : }
490 : }
491 883 : else if (ix == _nx - 1 && iy == 0) // Bottom right corner
492 : {
493 1377 : for (auto z : _z_grid)
494 : {
495 : Point z0{0, 0, z};
496 13260 : for (unsigned int i = 0; i < points_per_corner; i++)
497 11934 : mesh_base->add_point(br_corner_points[i] + x0 + y0 + z0, node_id++);
498 : }
499 : }
500 832 : else if (ix == 0 && iy == _ny - 1) // top Left corner
501 : {
502 1377 : for (auto z : _z_grid)
503 : {
504 : Point z0{0, 0, z};
505 13260 : for (unsigned int i = 0; i < points_per_corner; i++)
506 11934 : mesh_base->add_point(tl_corner_points[i] + x0 + y0 + z0, node_id++);
507 : }
508 : }
509 781 : else if (ix == _nx - 1 && iy == _ny - 1) // top right corner
510 : {
511 1377 : for (auto z : _z_grid)
512 : {
513 : Point z0{0, 0, z};
514 13260 : for (unsigned int i = 0; i < points_per_corner; i++)
515 11934 : mesh_base->add_point(tr_corner_points[i] + x0 + y0 + z0, node_id++);
516 : }
517 : }
518 730 : else if (ix == 0 && (iy != _ny - 1 || iy != 0)) // left side
519 : {
520 2907 : for (auto z : _z_grid)
521 : {
522 : Point z0{0, 0, z};
523 39494 : for (unsigned int i = 0; i < points_per_side; i++)
524 36673 : mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
525 : }
526 : }
527 644 : else if (ix == _nx - 1 && (iy != _ny - 1 || iy != 0)) // right side
528 : {
529 2907 : for (auto z : _z_grid)
530 : {
531 : Point z0{0, 0, z};
532 39494 : for (unsigned int i = 0; i < points_per_side; i++)
533 36673 : mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
534 : }
535 : }
536 558 : else if (iy == 0 && (ix != _nx - 1 || ix != 0)) // bottom side
537 : {
538 5307 : for (auto z : _z_grid)
539 : {
540 : Point z0{0, 0, z};
541 72422 : for (unsigned int i = 0; i < points_per_side; i++)
542 67249 : mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
543 : }
544 : }
545 424 : else if (iy == _ny - 1 && (ix != _nx - 1 || ix != 0)) // top side
546 : {
547 5307 : for (auto z : _z_grid)
548 : {
549 : Point z0{0, 0, z};
550 72422 : for (unsigned int i = 0; i < points_per_side; i++)
551 67249 : mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
552 : }
553 : }
554 : else // center
555 : {
556 11619 : for (auto z : _z_grid)
557 : {
558 : Point z0{0, 0, z};
559 249238 : for (unsigned int i = 0; i < points_per_center; i++)
560 237909 : mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
561 : }
562 : }
563 : }
564 : }
565 : }
566 :
567 : // Add elements to the mesh. The elements are 6-node prisms. The
568 : // bases of these prisms form a triangulated representation of a cross-section
569 : // of a center subchannel.
570 78 : if (_n_channels == 2)
571 : {
572 : unsigned int elem_id = 0;
573 33 : for (unsigned int iy = 0; iy < _ny; iy++)
574 : {
575 46 : for (unsigned int ix = 0; ix < _nx; ix++)
576 : {
577 26 : unsigned int i_ch = _nx * iy + ix;
578 1366 : for (unsigned int iz = 0; iz < _n_cells; iz++)
579 : {
580 17420 : for (unsigned int i = 0; i < elems_per_side; i++)
581 : {
582 16080 : Elem * elem = new Prism6;
583 16080 : elem->subdomain_id() = _block_id;
584 16080 : elem->set_id(elem_id++);
585 16080 : elem = mesh_base->add_elem(elem);
586 : // index of the central node at base of cell
587 16080 : unsigned int indx1 = iz * points_per_side + points_per_side * (_n_cells + 1) * i_ch;
588 : // index of the central node at top of cell
589 16080 : unsigned int indx2 =
590 : (iz + 1) * points_per_side + points_per_side * (_n_cells + 1) * i_ch;
591 : unsigned int elems_per_channel = elems_per_side;
592 16080 : elem->set_node(0, mesh_base->node_ptr(indx1));
593 16080 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
594 16080 : if (i != elems_per_channel - 1)
595 14740 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
596 : else
597 1340 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
598 :
599 16080 : elem->set_node(3, mesh_base->node_ptr(indx2));
600 16080 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
601 16080 : if (i != elems_per_channel - 1)
602 14740 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
603 : else
604 1340 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
605 :
606 16080 : if (iz == 0)
607 312 : boundary_info.add_side(elem, 0, 0);
608 16080 : if (iz == _n_cells - 1)
609 312 : boundary_info.add_side(elem, 4, 1);
610 : }
611 : }
612 : }
613 : }
614 13 : boundary_info.sideset_name(0) = "inlet";
615 13 : boundary_info.sideset_name(1) = "outlet";
616 13 : mesh_base->subdomain_name(_block_id) = name();
617 13 : mesh_base->prepare_for_use();
618 : }
619 65 : else if (_n_channels > 2 && (_ny == 1 || _nx == 1))
620 : {
621 : unsigned int elem_id = 0;
622 : unsigned int number_of_corner = 0;
623 : unsigned int number_of_side = 0;
624 : unsigned int number_of_center = 0;
625 : unsigned int elems_per_channel = 0;
626 : unsigned int points_per_channel = 0;
627 42 : for (unsigned int iy = 0; iy < _ny; iy++)
628 : {
629 70 : for (unsigned int ix = 0; ix < _nx; ix++)
630 : {
631 42 : unsigned int i_ch = _nx * iy + ix;
632 : auto subch_type = getSubchannelType(i_ch);
633 : // note that in this case i use side geometry for corner subchannel
634 42 : if (subch_type == EChannelType::CORNER)
635 : {
636 28 : number_of_side++;
637 : elems_per_channel = elems_per_side;
638 : points_per_channel = points_per_side;
639 : }
640 : // note that in this case i use center geometry for edge subchannel
641 14 : else if (subch_type == EChannelType::EDGE)
642 : {
643 14 : number_of_center++;
644 : elems_per_channel = elems_per_center;
645 : points_per_channel = points_per_center;
646 : }
647 462 : for (unsigned int iz = 0; iz < _n_cells; iz++)
648 : {
649 420 : unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
650 420 : number_of_side * points_per_side * (_n_cells + 1) +
651 : number_of_center * points_per_center * (_n_cells + 1) -
652 420 : points_per_channel * (_n_cells + 1);
653 : // index of the central node at base of cell
654 420 : unsigned int indx1 = iz * points_per_channel + elapsed_points;
655 : // index of the central node at top of cell
656 420 : unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
657 :
658 6580 : for (unsigned int i = 0; i < elems_per_channel; i++)
659 : {
660 6160 : Elem * elem = new Prism6;
661 6160 : elem->subdomain_id() = _block_id;
662 6160 : elem->set_id(elem_id++);
663 6160 : elem = mesh_base->add_elem(elem);
664 :
665 6160 : elem->set_node(0, mesh_base->node_ptr(indx1));
666 6160 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
667 6160 : if (i != elems_per_channel - 1)
668 5740 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
669 : else
670 420 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
671 :
672 6160 : elem->set_node(3, mesh_base->node_ptr(indx2));
673 6160 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
674 6160 : if (i != elems_per_channel - 1)
675 5740 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
676 : else
677 420 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
678 :
679 6160 : if (iz == 0)
680 616 : boundary_info.add_side(elem, 0, 0);
681 6160 : if (iz == _n_cells - 1)
682 616 : boundary_info.add_side(elem, 4, 1);
683 : }
684 : }
685 : }
686 : }
687 14 : boundary_info.sideset_name(0) = "inlet";
688 14 : boundary_info.sideset_name(1) = "outlet";
689 14 : mesh_base->subdomain_name(_block_id) = name();
690 14 : mesh_base->prepare_for_use();
691 : }
692 : else
693 : {
694 : unsigned int elem_id = 0;
695 : unsigned int number_of_corner = 0;
696 : unsigned int number_of_side = 0;
697 : unsigned int number_of_center = 0;
698 : unsigned int elems_per_channel = 0;
699 : unsigned int points_per_channel = 0;
700 239 : for (unsigned int iy = 0; iy < _ny; iy++)
701 : {
702 1122 : for (unsigned int ix = 0; ix < _nx; ix++)
703 : {
704 934 : unsigned int i_ch = _nx * iy + ix;
705 : auto subch_type = getSubchannelType(i_ch);
706 934 : if (subch_type == EChannelType::CORNER)
707 : {
708 204 : number_of_corner++;
709 : elems_per_channel = elems_per_corner;
710 : points_per_channel = points_per_corner;
711 : }
712 730 : else if (subch_type == EChannelType::EDGE)
713 : {
714 440 : number_of_side++;
715 : elems_per_channel = elems_per_side;
716 : points_per_channel = points_per_side;
717 : }
718 : else
719 : {
720 290 : number_of_center++;
721 : elems_per_channel = elems_per_center;
722 : points_per_channel = points_per_center;
723 : }
724 32621 : for (unsigned int iz = 0; iz < _n_cells; iz++)
725 : {
726 31687 : unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
727 31687 : number_of_side * points_per_side * (_n_cells + 1) +
728 : number_of_center * points_per_center * (_n_cells + 1) -
729 31687 : points_per_channel * (_n_cells + 1);
730 : // index of the central node at base of cell
731 31687 : unsigned int indx1 = iz * points_per_channel + elapsed_points;
732 : // index of the central node at top of cell
733 31687 : unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
734 :
735 479843 : for (unsigned int i = 0; i < elems_per_channel; i++)
736 : {
737 448156 : Elem * elem = new Prism6;
738 448156 : elem->subdomain_id() = _block_id;
739 448156 : elem->set_id(elem_id++);
740 448156 : elem = mesh_base->add_elem(elem);
741 :
742 448156 : elem->set_node(0, mesh_base->node_ptr(indx1));
743 448156 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
744 448156 : if (i != elems_per_channel - 1)
745 416469 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
746 : else
747 31687 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
748 :
749 448156 : elem->set_node(3, mesh_base->node_ptr(indx2));
750 448156 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
751 448156 : if (i != elems_per_channel - 1)
752 416469 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
753 : else
754 31687 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
755 :
756 448156 : if (iz == 0)
757 12712 : boundary_info.add_side(elem, 0, 0);
758 448156 : if (iz == _n_cells - 1)
759 12712 : boundary_info.add_side(elem, 4, 1);
760 : }
761 : }
762 : }
763 : }
764 51 : boundary_info.sideset_name(0) = "inlet";
765 51 : boundary_info.sideset_name(1) = "outlet";
766 51 : mesh_base->subdomain_name(_block_id) = name();
767 51 : mesh_base->prepare_for_use();
768 : }
769 :
770 78 : return mesh_base;
771 0 : }
|
