https://mooseframework.inl.gov
SCMDetailedQuadSubChannelMeshGenerator.C
Go to the documentation of this file.
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 
11 #include <array>
12 #include <cmath>
13 #include "libmesh/cell_prism6.h"
14 #include "libmesh/unstructured_mesh.h"
15 
17 registerMooseObjectRenamed("SubChannelApp",
18  DetailedQuadSubChannelMeshGenerator,
19  "06/30/2025 24:00",
21 
24 {
26  params.addClassDescription(
27  "Creates a detailed mesh of subchannels in a square lattice arrangement");
28  params.addRequiredParam<Real>("pitch", "Pitch [m]");
29  params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
30  params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
31  params.addRequiredParam<Real>("heated_length", "Heated length [m]");
32  params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
33  params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
34  params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
35  params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");
36  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  params.deprecateParam("gap", "side_gap", "08/06/2026");
42  params.addParam<unsigned int>("block_id", 0, "Block ID used for the mesh subdomain.");
43  return params;
44 }
45 
47  const InputParameters & parameters)
48  : MeshGenerator(parameters),
49  _unheated_length_entry(getParam<Real>("unheated_length_entry")),
50  _heated_length(getParam<Real>("heated_length")),
51  _unheated_length_exit(getParam<Real>("unheated_length_exit")),
52  _pitch(getParam<Real>("pitch")),
53  _pin_diameter(getParam<Real>("pin_diameter")),
54  _n_cells(getParam<unsigned int>("n_cells")),
55  _nx(getParam<unsigned int>("nx")),
56  _ny(getParam<unsigned int>("ny")),
57  _n_channels(_nx * _ny),
58  _side_gap(getParam<Real>("side_gap")),
59  _block_id(getParam<unsigned int>("block_id"))
60 {
62  Real dz = L / _n_cells;
63  for (unsigned int i = 0; i < _n_cells + 1; i++)
64  _z_grid.push_back(dz * i);
65 
67  for (unsigned int i = 0; i < _n_channels; i++)
68  {
69  _subchannel_position[i].reserve(3);
70  for (unsigned int j = 0; j < 3; j++)
71  {
72  _subchannel_position.at(i).push_back(0.0);
73  }
74  }
75 
76  _subch_type.resize(_n_channels);
77  for (unsigned int iy = 0; iy < _ny; iy++)
78  {
79  for (unsigned int ix = 0; ix < _nx; ix++)
80  {
81  unsigned int i_ch = _nx * iy + ix;
82  bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
83  (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
84  bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);
85 
86  if (is_corner)
88  else if (is_edge)
90  else
92 
93  // set the subchannel positions so that the center of the assembly is the zero point
94  Real offset_x = (_nx - 1) * _pitch / 2.0;
95  Real offset_y = (_ny - 1) * _pitch / 2.0;
96  _subchannel_position[i_ch][0] = _pitch * ix - offset_x;
97  _subchannel_position[i_ch][1] = _pitch * iy - offset_y;
98  }
99  }
100 }
101 
102 std::unique_ptr<MeshBase>
104 {
105  auto mesh_base = buildMeshBaseObject();
106  BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
107  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  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  else if (_n_channels > 2 && (_ny == 1 || _nx == 1))
140  {
141  n_corner = 0;
142  n_side = 2;
143  n_center = _n_channels - n_side - n_corner;
144  }
145  else
146  {
147  n_corner = 4;
148  n_side = 2 * (_ny - 2) + 2 * (_nx - 2);
149  n_center = _n_channels - n_side - n_corner;
150  }
151 
152  // Compute the total number of points and elements.
153  const unsigned int points_per_level =
154  n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
155  const unsigned int elems_per_level =
156  n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
157  const unsigned int n_points = points_per_level * (_n_cells + 1);
158  const unsigned int n_elems = elems_per_level * _n_cells;
159  mesh_base->reserve_nodes(n_points);
160  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  const Real radius = _pin_diameter / 2.0;
163  std::array<Point, theta_res + 1> circle_points;
164  {
165  Real theta = 0;
166  for (unsigned int i = 0; i < theta_res + 1; i++)
167  {
168  circle_points[i](0) = radius * std::cos(theta);
169  circle_points[i](1) = radius * std::sin(theta);
170  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  quadrant_centers[0] = Point(_pitch * 0.5 * shrink_factor, _pitch * 0.5 * shrink_factor, 0);
181  quadrant_centers[1] = Point(-_pitch * 0.5 * shrink_factor, _pitch * 0.5 * shrink_factor, 0);
182  quadrant_centers[2] = Point(-_pitch * 0.5 * shrink_factor, -_pitch * 0.5 * shrink_factor, 0);
183  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  for (unsigned int i = 0; i < 4; i++)
197  {
198  if (i == 0)
199  start = 3 * m;
200  if (i == 1)
201  start = 4 * m;
202  if (i == 2)
203  start = 1 * m;
204  if (i == 3)
205  start = 2 * m;
206  for (unsigned int ii = 0; ii < points_per_quad; ii++)
207  {
208  auto c_pt = circle_points[start - ii];
209  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
224  {
225  auto c_pt = circle_points[2 * m - ii];
226  tl_corner_points[ii + 1] = quadrant_centers[3] + c_pt;
227  }
228  tl_corner_points[points_per_quad + 1] = Point(_pitch * 0.5 * shrink_factor, _side_gap, 0);
229  tl_corner_points[points_per_quad + 2] = Point(-_side_gap, _side_gap, 0);
230  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
243  {
244  auto c_pt = circle_points[m - ii];
245  tr_corner_points[ii + 1] = quadrant_centers[2] + c_pt;
246  }
247  tr_corner_points[points_per_quad + 1] = Point(_side_gap, -_pitch * 0.5 * shrink_factor, 0);
248  tr_corner_points[points_per_quad + 2] = Point(_side_gap, _side_gap, 0);
249  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
262  {
263  auto c_pt = circle_points[3 * m - ii];
264  bl_corner_points[ii + 1] = quadrant_centers[0] + c_pt;
265  }
266  bl_corner_points[points_per_quad + 1] = Point(-_side_gap, _pitch * 0.5 * shrink_factor, 0);
267  bl_corner_points[points_per_quad + 2] = Point(-_side_gap, -_side_gap, 0);
268  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
281  {
282  auto c_pt = circle_points[4 * m - ii];
283  br_corner_points[ii + 1] = quadrant_centers[1] + c_pt;
284  }
285  br_corner_points[points_per_quad + 1] = Point(-_pitch * 0.5 * shrink_factor, -_side_gap, 0);
286  br_corner_points[points_per_quad + 2] = Point(_side_gap, -_side_gap, 0);
287  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
300  {
301  auto c_pt = circle_points[m - ii];
302  top_points[ii + 1] = quadrant_centers[2] + c_pt;
303  }
304  for (unsigned int ii = 0; ii < points_per_quad; ii++)
305  {
306  auto c_pt = circle_points[2 * m - ii];
307  top_points[points_per_quad + ii + 1] = quadrant_centers[3] + c_pt;
308  }
309  top_points[2 * points_per_quad + 1] = Point(_pitch * 0.5 * shrink_factor, _side_gap, 0);
310  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
323  {
324  auto c_pt = circle_points[2 * m - ii];
325  left_points[ii + 1] = quadrant_centers[3] + c_pt;
326  }
327  for (unsigned int ii = 0; ii < points_per_quad; ii++)
328  {
329  auto c_pt = circle_points[3 * m - ii];
330  left_points[points_per_quad + ii + 1] = quadrant_centers[0] + c_pt;
331  }
332  left_points[2 * points_per_quad + 1] = Point(-_side_gap, _pitch * 0.5 * shrink_factor, 0);
333  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
346  {
347  auto c_pt = circle_points[3 * m - ii];
348  bottom_points[ii + 1] = quadrant_centers[0] + c_pt;
349  }
350  for (unsigned int ii = 0; ii < points_per_quad; ii++)
351  {
352  auto c_pt = circle_points[4 * m - ii];
353  bottom_points[points_per_quad + ii + 1] = quadrant_centers[1] + c_pt;
354  }
355  bottom_points[2 * points_per_quad + 1] = Point(-_pitch * 0.5 * shrink_factor, -_side_gap, 0);
356  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  for (unsigned int ii = 0; ii < points_per_quad; ii++)
369  {
370  auto c_pt = circle_points[4 * m - ii];
371  right_points[ii + 1] = quadrant_centers[1] + c_pt;
372  }
373  for (unsigned int ii = 0; ii < points_per_quad; ii++)
374  {
375  auto c_pt = circle_points[m - ii];
376  right_points[points_per_quad + ii + 1] = quadrant_centers[2] + c_pt;
377  }
378  right_points[2 * points_per_quad + 1] = Point(_side_gap, -_pitch * 0.5 * shrink_factor, 0);
379  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  if (_n_channels == 2)
384  {
385  unsigned int node_id = 0;
386  Real offset_x = (_nx - 1) * _pitch / 2.0;
387  Real offset_y = (_ny - 1) * _pitch / 2.0;
388  for (unsigned int iy = 0; iy < _ny; iy++)
389  {
390  Point y0 = {0, _pitch * iy - offset_y, 0};
391  for (unsigned int ix = 0; ix < _nx; ix++)
392  {
393  Point x0 = {_pitch * ix - offset_x, 0, 0};
394  for (auto z : _z_grid)
395  {
396  Point z0{0, 0, z};
397  if (_nx == 1 && iy == 0) // vertical orientation
398  {
399  for (unsigned int i = 0; i < points_per_side; i++)
400  mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
401  }
402  if (_nx == 1 && iy == 1) // vertical orientation
403  {
404  for (unsigned int i = 0; i < points_per_side; i++)
405  mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
406  }
407  if (_ny == 1 && ix == 0) // horizontal orientation
408  {
409  for (unsigned int i = 0; i < points_per_side; i++)
410  mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
411  }
412  if (_ny == 1 && ix == 1) // horizontal orientation
413  {
414  for (unsigned int i = 0; i < points_per_side; i++)
415  mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
416  }
417  }
418  }
419  }
420  }
421  else if (_n_channels > 2 && (_ny == 1 || _nx == 1))
422  {
423  unsigned int node_id = 0;
424  Real offset_x = (_nx - 1) * _pitch / 2.0;
425  Real offset_y = (_ny - 1) * _pitch / 2.0;
426  for (unsigned int iy = 0; iy < _ny; iy++)
427  {
428  Point y0 = {0, _pitch * iy - offset_y, 0};
429  for (unsigned int ix = 0; ix < _nx; ix++)
430  {
431  Point x0 = {_pitch * ix - offset_x, 0, 0};
432  for (auto z : _z_grid)
433  {
434  Point z0{0, 0, z};
435  if (_nx == 1)
436  {
437  if (iy == 0)
438  {
439  for (unsigned int i = 0; i < points_per_side; i++)
440  mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
441  }
442  else if (iy == _ny - 1)
443  {
444  for (unsigned int i = 0; i < points_per_side; i++)
445  mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
446  }
447  else
448  {
449  for (unsigned int i = 0; i < points_per_center; i++)
450  mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
451  }
452  }
453  else if (_ny == 1)
454  {
455  if (ix == 0)
456  {
457  for (unsigned int i = 0; i < points_per_side; i++)
458  mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
459  }
460  else if (ix == _nx - 1)
461  {
462  for (unsigned int i = 0; i < points_per_side; i++)
463  mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
464  }
465  else
466  {
467  for (unsigned int i = 0; i < points_per_center; i++)
468  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  Real offset_x = (_nx - 1) * _pitch / 2.0;
479  Real offset_y = (_ny - 1) * _pitch / 2.0;
480  for (unsigned int iy = 0; iy < _ny; iy++)
481  {
482  Point y0 = {0, _pitch * iy - offset_y, 0};
483  for (unsigned int ix = 0; ix < _nx; ix++)
484  {
485  Point x0 = {_pitch * ix - offset_x, 0, 0};
486  if (ix == 0 && iy == 0) // Bottom Left corner
487  {
488  for (auto z : _z_grid)
489  {
490  Point z0{0, 0, z};
491  for (unsigned int i = 0; i < points_per_corner; i++)
492  mesh_base->add_point(bl_corner_points[i] + x0 + y0 + z0, node_id++);
493  }
494  }
495  else if (ix == _nx - 1 && iy == 0) // Bottom right corner
496  {
497  for (auto z : _z_grid)
498  {
499  Point z0{0, 0, z};
500  for (unsigned int i = 0; i < points_per_corner; i++)
501  mesh_base->add_point(br_corner_points[i] + x0 + y0 + z0, node_id++);
502  }
503  }
504  else if (ix == 0 && iy == _ny - 1) // top Left corner
505  {
506  for (auto z : _z_grid)
507  {
508  Point z0{0, 0, z};
509  for (unsigned int i = 0; i < points_per_corner; i++)
510  mesh_base->add_point(tl_corner_points[i] + x0 + y0 + z0, node_id++);
511  }
512  }
513  else if (ix == _nx - 1 && iy == _ny - 1) // top right corner
514  {
515  for (auto z : _z_grid)
516  {
517  Point z0{0, 0, z};
518  for (unsigned int i = 0; i < points_per_corner; i++)
519  mesh_base->add_point(tr_corner_points[i] + x0 + y0 + z0, node_id++);
520  }
521  }
522  else if (ix == 0 && (iy != _ny - 1 || iy != 0)) // left side
523  {
524  for (auto z : _z_grid)
525  {
526  Point z0{0, 0, z};
527  for (unsigned int i = 0; i < points_per_side; i++)
528  mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
529  }
530  }
531  else if (ix == _nx - 1 && (iy != _ny - 1 || iy != 0)) // right side
532  {
533  for (auto z : _z_grid)
534  {
535  Point z0{0, 0, z};
536  for (unsigned int i = 0; i < points_per_side; i++)
537  mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
538  }
539  }
540  else if (iy == 0 && (ix != _nx - 1 || ix != 0)) // bottom side
541  {
542  for (auto z : _z_grid)
543  {
544  Point z0{0, 0, z};
545  for (unsigned int i = 0; i < points_per_side; i++)
546  mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
547  }
548  }
549  else if (iy == _ny - 1 && (ix != _nx - 1 || ix != 0)) // top side
550  {
551  for (auto z : _z_grid)
552  {
553  Point z0{0, 0, z};
554  for (unsigned int i = 0; i < points_per_side; i++)
555  mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
556  }
557  }
558  else // center
559  {
560  for (auto z : _z_grid)
561  {
562  Point z0{0, 0, z};
563  for (unsigned int i = 0; i < points_per_center; i++)
564  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  if (_n_channels == 2)
575  {
576  unsigned int elem_id = 0;
577  for (unsigned int iy = 0; iy < _ny; iy++)
578  {
579  for (unsigned int ix = 0; ix < _nx; ix++)
580  {
581  unsigned int i_ch = _nx * iy + ix;
582  for (unsigned int iz = 0; iz < _n_cells; iz++)
583  {
584  for (unsigned int i = 0; i < elems_per_side; i++)
585  {
586  Elem * elem = new Prism6;
587  elem->subdomain_id() = _block_id;
588  elem->set_id(elem_id++);
589  elem = mesh_base->add_elem(elem);
590  // index of the central node at base of cell
591  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  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  elem->set_node(0, mesh_base->node_ptr(indx1));
597  elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
598  if (i != elems_per_channel - 1)
599  elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
600  else
601  elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
602 
603  elem->set_node(3, mesh_base->node_ptr(indx2));
604  elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
605  if (i != elems_per_channel - 1)
606  elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
607  else
608  elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
609 
610  if (iz == 0)
611  boundary_info.add_side(elem, 0, 0);
612  if (iz == _n_cells - 1)
613  boundary_info.add_side(elem, 4, 1);
614  }
615  }
616  }
617  }
618  boundary_info.sideset_name(0) = "inlet";
619  boundary_info.sideset_name(1) = "outlet";
620  mesh_base->subdomain_name(_block_id) = name();
621  mesh_base->prepare_for_use();
622  }
623  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  for (unsigned int iy = 0; iy < _ny; iy++)
632  {
633  for (unsigned int ix = 0; ix < _nx; ix++)
634  {
635  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  if (subch_type == EChannelType::CORNER)
639  {
640  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  else if (subch_type == EChannelType::EDGE)
646  {
647  number_of_center++;
648  elems_per_channel = elems_per_center;
649  points_per_channel = points_per_center;
650  }
651  for (unsigned int iz = 0; iz < _n_cells; iz++)
652  {
653  unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
654  number_of_side * points_per_side * (_n_cells + 1) +
655  number_of_center * points_per_center * (_n_cells + 1) -
656  points_per_channel * (_n_cells + 1);
657  // index of the central node at base of cell
658  unsigned int indx1 = iz * points_per_channel + elapsed_points;
659  // index of the central node at top of cell
660  unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
661 
662  for (unsigned int i = 0; i < elems_per_channel; i++)
663  {
664  Elem * elem = new Prism6;
665  elem->subdomain_id() = _block_id;
666  elem->set_id(elem_id++);
667  elem = mesh_base->add_elem(elem);
668 
669  elem->set_node(0, mesh_base->node_ptr(indx1));
670  elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
671  if (i != elems_per_channel - 1)
672  elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
673  else
674  elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
675 
676  elem->set_node(3, mesh_base->node_ptr(indx2));
677  elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
678  if (i != elems_per_channel - 1)
679  elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
680  else
681  elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
682 
683  if (iz == 0)
684  boundary_info.add_side(elem, 0, 0);
685  if (iz == _n_cells - 1)
686  boundary_info.add_side(elem, 4, 1);
687  }
688  }
689  }
690  }
691  boundary_info.sideset_name(0) = "inlet";
692  boundary_info.sideset_name(1) = "outlet";
693  mesh_base->subdomain_name(_block_id) = name();
694  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  for (unsigned int iy = 0; iy < _ny; iy++)
705  {
706  for (unsigned int ix = 0; ix < _nx; ix++)
707  {
708  unsigned int i_ch = _nx * iy + ix;
709  auto subch_type = getSubchannelType(i_ch);
710  if (subch_type == EChannelType::CORNER)
711  {
712  number_of_corner++;
713  elems_per_channel = elems_per_corner;
714  points_per_channel = points_per_corner;
715  }
716  else if (subch_type == EChannelType::EDGE)
717  {
718  number_of_side++;
719  elems_per_channel = elems_per_side;
720  points_per_channel = points_per_side;
721  }
722  else
723  {
724  number_of_center++;
725  elems_per_channel = elems_per_center;
726  points_per_channel = points_per_center;
727  }
728  for (unsigned int iz = 0; iz < _n_cells; iz++)
729  {
730  unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
731  number_of_side * points_per_side * (_n_cells + 1) +
732  number_of_center * points_per_center * (_n_cells + 1) -
733  points_per_channel * (_n_cells + 1);
734  // index of the central node at base of cell
735  unsigned int indx1 = iz * points_per_channel + elapsed_points;
736  // index of the central node at top of cell
737  unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
738 
739  for (unsigned int i = 0; i < elems_per_channel; i++)
740  {
741  Elem * elem = new Prism6;
742  elem->subdomain_id() = _block_id;
743  elem->set_id(elem_id++);
744  elem = mesh_base->add_elem(elem);
745 
746  elem->set_node(0, mesh_base->node_ptr(indx1));
747  elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
748  if (i != elems_per_channel - 1)
749  elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
750  else
751  elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
752 
753  elem->set_node(3, mesh_base->node_ptr(indx2));
754  elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
755  if (i != elems_per_channel - 1)
756  elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
757  else
758  elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
759 
760  if (iz == 0)
761  boundary_info.add_side(elem, 0, 0);
762  if (iz == _n_cells - 1)
763  boundary_info.add_side(elem, 4, 1);
764  }
765  }
766  }
767  }
768  boundary_info.sideset_name(0) = "inlet";
769  boundary_info.sideset_name(1) = "outlet";
770  mesh_base->subdomain_name(_block_id) = name();
771  mesh_base->prepare_for_use();
772  }
773 
774  return mesh_base;
775 }
registerMooseObject("SubChannelApp", SCMDetailedQuadSubChannelMeshGenerator)
registerMooseObjectRenamed("SubChannelApp", DetailedQuadSubChannelMeshGenerator, "06/30/2025 24:00", SCMDetailedQuadSubChannelMeshGenerator)
const Real radius
const Real _pitch
Distance between the neighbor fuel pins, pitch.
const unsigned int _n_cells
Number of cells in the axial direction.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const Real _heated_length
heated length of the fuel Pin
std::vector< Real > _z_grid
axial location of nodes
unsigned int _n_channels
Total number of subchannels.
void addRequiredParam(const std::string &name, const std::string &doc_string)
const unsigned int _nx
Number of subchannels in the x direction.
const std::string & name() const
std::vector< std::vector< Real > > _subchannel_position
x,y coordinates of the subchannel centroids
void deprecateParam(const std::string &old_name, const std::string &new_name, const std::string &removal_date)
const Real _side_gap
The side gap, not to be confused with the gap between pins, this refers to the gap next to the duct o...
static InputParameters validParams()
SCMDetailedQuadSubChannelMeshGenerator(const InputParameters &parameters)
const unsigned int _ny
Number of subchannels in the y direction.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Mesh generator that builds a 3D mesh representing quadrilateral subchannels.
const Real _unheated_length_exit
unheated length of the fuel Pin at the exit of the assembly
std::vector< EChannelType > _subch_type
Subchannel type.
void addClassDescription(const std::string &doc_string)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
EChannelType getSubchannelType(unsigned int index) const
std::unique_ptr< MeshBase > buildMeshBaseObject(unsigned int dim=libMesh::invalid_uint)
const Real _unheated_length_entry
unheated length of the fuel Pin at the entry of the assembly
const unsigned int & _block_id
Subdomain ID used for the mesh block.
virtual std::unique_ptr< MeshBase > generate() override
void ErrorVector unsigned int