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