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