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 "SCMDetailedQuadInterWrapperMeshGenerator.h"
11 : #include <array>
12 : #include <cmath>
13 : #include "libmesh/cell_prism6.h"
14 : #include "libmesh/unstructured_mesh.h"
15 :
16 : registerMooseObject("SubChannelApp", SCMDetailedQuadInterWrapperMeshGenerator);
17 : registerMooseObjectRenamed("SubChannelApp",
18 : DetailedQuadInterWrapperMeshGenerator,
19 : "06/30/2025 24:00",
20 : SCMDetailedQuadInterWrapperMeshGenerator);
21 :
22 : InputParameters
23 48 : SCMDetailedQuadInterWrapperMeshGenerator::validParams()
24 : {
25 48 : InputParameters params = MeshGenerator::validParams();
26 48 : params.addClassDescription(
27 : "Creates a detailed mesh of the inter-wrapper cells around square lattice subassemblies");
28 96 : params.addRequiredParam<Real>("assembly_pitch", "Assembly Pitch [m]");
29 96 : params.addRequiredParam<Real>("assembly_side_x",
30 : "Outer side lengths of assembly in x [m] - including duct");
31 96 : params.addRequiredParam<Real>("assembly_side_y",
32 : "Outer side lengths of assembly in y [m] - including duct");
33 96 : params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
34 96 : params.addRequiredParam<Real>("heated_length", "Heated length [m]");
35 96 : params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
36 96 : params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
37 96 : params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
38 96 : params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");
39 96 : params.addRequiredParam<Real>("side_bypass", "Half of assembly_pitch between assemblies [m]");
40 96 : params.addParam<unsigned int>("block_id", 0, "Block ID used for the mesh subdomain.");
41 48 : return params;
42 0 : }
43 :
44 24 : SCMDetailedQuadInterWrapperMeshGenerator::SCMDetailedQuadInterWrapperMeshGenerator(
45 24 : const InputParameters & parameters)
46 : : MeshGenerator(parameters),
47 24 : _unheated_length_entry(getParam<Real>("unheated_length_entry")),
48 48 : _heated_length(getParam<Real>("heated_length")),
49 48 : _unheated_length_exit(getParam<Real>("unheated_length_exit")),
50 48 : _assembly_pitch(getParam<Real>("assembly_pitch")),
51 48 : _assembly_side_x(getParam<Real>("assembly_side_x")),
52 48 : _assembly_side_y(getParam<Real>("assembly_side_y")),
53 48 : _n_cells(getParam<unsigned int>("n_cells")),
54 48 : _nx(getParam<unsigned int>("nx") + 1),
55 48 : _ny(getParam<unsigned int>("ny") + 1),
56 24 : _n_channels((_nx + 1) * (_ny + 1)),
57 48 : _side_bypass_length(getParam<Real>("side_bypass")),
58 72 : _block_id(getParam<unsigned int>("block_id"))
59 : {
60 24 : Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
61 24 : Real dz = L / _n_cells;
62 648 : for (unsigned int i = 0; i < _n_cells + 1; i++)
63 624 : _z_grid.push_back(dz * i);
64 :
65 24 : _subch_type.resize(_n_channels);
66 240 : for (unsigned int iy = 0; iy < _ny; iy++)
67 : {
68 2880 : for (unsigned int ix = 0; ix < _nx; ix++)
69 : {
70 2664 : unsigned int i_ch = _nx * iy + ix;
71 2640 : bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
72 5280 : (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
73 2664 : bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);
74 :
75 2664 : if (is_corner)
76 96 : _subch_type[i_ch] = EChannelType::CORNER;
77 2568 : else if (is_edge)
78 672 : _subch_type[i_ch] = EChannelType::EDGE;
79 : else
80 1896 : _subch_type[i_ch] = EChannelType::CENTER;
81 : }
82 : }
83 24 : }
84 :
85 : std::unique_ptr<MeshBase>
86 24 : SCMDetailedQuadInterWrapperMeshGenerator::generate()
87 : {
88 24 : auto mesh_base = buildMeshBaseObject();
89 : BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
90 24 : mesh_base->set_spatial_dimension(3);
91 : // Define the resolution (the number of points used to represent a circle).
92 : // This must be divisible by 4.
93 : const unsigned int theta_res = 16; // TODO: parameterize
94 : // Compute the number of points needed to represent one quarter of a circle.
95 : const unsigned int points_per_quad = theta_res / 4 + 1;
96 :
97 : // Compute the points needed to represent one axial cross-flow of a subchannel.
98 : // For the center subchannel (sc) there is one center point plus the points from 4 intersecting
99 : // circles.
100 : const unsigned int points_per_center = points_per_quad * 4 + 1;
101 : // For the corner sc there is one center point plus the points from 1 intersecting circle plus 3
102 : // corners
103 : const unsigned int points_per_corner = points_per_quad * 1 + 1 + 3;
104 : // For the side sc there is one center point plus the points from 2 intersecting circles plus 2
105 : // corners
106 : const unsigned int points_per_side = points_per_quad * 2 + 1 + 2;
107 :
108 : // Compute the number of elements (Prism6) which combined base creates the sub-channel
109 : // cross-section
110 : const unsigned int elems_per_center = theta_res + 4;
111 : const unsigned int elems_per_corner = theta_res / 4 + 4;
112 : const unsigned int elems_per_side = theta_res / 2 + 4;
113 :
114 : // specify number and type of sub-channel
115 : unsigned int n_center, n_side, n_corner;
116 :
117 : n_corner = 4;
118 24 : n_side = 2 * (_ny - 2) + 2 * (_nx - 2);
119 24 : n_center = _n_channels - n_side - n_corner;
120 :
121 : // Compute the total number of points and elements.
122 24 : const unsigned int points_per_level =
123 24 : n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
124 24 : const unsigned int elems_per_level =
125 24 : n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
126 24 : const unsigned int n_points = points_per_level * (_n_cells + 1);
127 24 : const unsigned int n_elems = elems_per_level * _n_cells;
128 24 : mesh_base->reserve_nodes(n_points);
129 24 : mesh_base->reserve_elem(n_elems);
130 : // Build an array of points arranged in a square on the xy-plane. (last and first node overlap)
131 : // Mapping points in a circle to a square
132 : // See ttp://arxiv.org/abs/1509.06344 for proof
133 : const Real radius = 1.0;
134 : std::array<Point, theta_res + 1> circle_points;
135 : {
136 : Real theta = 0;
137 : Real u, v;
138 432 : for (unsigned int i = 0; i < theta_res + 1; i++)
139 : {
140 408 : u = radius * std::cos(theta);
141 408 : v = radius * std::sin(theta);
142 408 : Real val_check_u = std::abs(u) - std::sqrt(2.0) / 2.0;
143 408 : Real val_check_v = std::abs(v) - std::sqrt(2.0) / 2.0;
144 408 : if (val_check_u < 1e-5 && val_check_v < 1e-5)
145 : {
146 96 : circle_points[i](0) = u * 2.0 / std::sqrt(2);
147 96 : circle_points[i](1) = v * 2.0 / std::sqrt(2);
148 : }
149 : else
150 : {
151 312 : circle_points[i](0) =
152 312 : 0.5 * std::sqrt(2. + std::pow(u, 2) - std::pow(v, 2) + 2. * u * std::sqrt(2)) -
153 312 : 0.5 * std::sqrt(2. + std::pow(u, 2) - std::pow(v, 2) - 2. * u * std::sqrt(2));
154 312 : circle_points[i](1) =
155 312 : 0.5 * std::sqrt(2. - std::pow(u, 2) + std::pow(v, 2) + 2. * v * std::sqrt(2)) -
156 312 : 0.5 * std::sqrt(2. - std::pow(u, 2) + std::pow(v, 2) - 2. * v * std::sqrt(2));
157 : }
158 408 : circle_points[i](0) *= _assembly_side_x / 2.0;
159 408 : circle_points[i](1) *= _assembly_side_y / 2.0;
160 408 : theta += 2 * M_PI / theta_res;
161 : }
162 : }
163 : // Define "quadrant center" reference points. These will be the centers of
164 : // the 4 circles that represent the fuel pins. These centers are
165 : // offset a little bit so that in the final mesh, there is a tiny assembly_pitch between
166 : // neighboring subchannel cells. That allows us to easily map a solution to
167 : // this detailed mesh with a nearest-neighbor search.
168 : const Real shrink_factor = 0.99999;
169 : std::array<Point, 4> quadrant_centers;
170 24 : quadrant_centers[0] =
171 24 : Point(_assembly_pitch * 0.5 * shrink_factor, _assembly_pitch * 0.5 * shrink_factor, 0);
172 24 : quadrant_centers[1] =
173 : Point(-_assembly_pitch * 0.5 * shrink_factor, _assembly_pitch * 0.5 * shrink_factor, 0);
174 24 : quadrant_centers[2] =
175 : Point(-_assembly_pitch * 0.5 * shrink_factor, -_assembly_pitch * 0.5 * shrink_factor, 0);
176 24 : quadrant_centers[3] =
177 : Point(_assembly_pitch * 0.5 * shrink_factor, -_assembly_pitch * 0.5 * shrink_factor, 0);
178 :
179 : const unsigned int m = theta_res / 4;
180 : // Build an array of points that represent a cross section of a center subchannel
181 : // cell. The points are ordered in this fashion:
182 : // 4 3
183 : // 6 5 2 1
184 : // 0
185 : // 7 8 * *
186 : // 9 *
187 : std::array<Point, points_per_center> center_points;
188 : {
189 : unsigned int start;
190 120 : for (unsigned int i = 0; i < 4; i++)
191 : {
192 96 : if (i == 0)
193 : start = 3 * m;
194 72 : if (i == 1)
195 : start = 4 * m;
196 72 : if (i == 2)
197 : start = 1 * m;
198 72 : if (i == 3)
199 : start = 2 * m;
200 576 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
201 : {
202 480 : auto c_pt = circle_points[start - ii];
203 480 : center_points[i * points_per_quad + ii + 1] = quadrant_centers[i] + c_pt;
204 : }
205 : }
206 : }
207 :
208 : // Build an array of points that represent a cross section of a top left corner subchannel
209 : // cell. The points are ordered in this fashion:
210 : // 5 4
211 : //
212 : // 0
213 : // 2 3
214 : // 6 1
215 : std::array<Point, points_per_corner> tl_corner_points;
216 : {
217 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
218 : {
219 120 : auto c_pt = circle_points[2 * m - ii];
220 120 : tl_corner_points[ii + 1] = quadrant_centers[3] + c_pt;
221 : }
222 24 : tl_corner_points[points_per_quad + 1] =
223 24 : Point(_assembly_pitch * 0.5 * shrink_factor, _side_bypass_length, 0);
224 24 : tl_corner_points[points_per_quad + 2] = Point(-_side_bypass_length, +_side_bypass_length, 0);
225 24 : tl_corner_points[points_per_quad + 3] =
226 : Point(-_side_bypass_length, -_assembly_pitch * 0.5 * shrink_factor, 0);
227 : }
228 :
229 : // Build an array of points that represent a cross section of a top right corner subchannel
230 : // cell. The points are ordered in this fashion:
231 : // 6 5
232 : //
233 : // 0
234 : // 1 2
235 : // 3 4
236 : std::array<Point, points_per_corner> tr_corner_points;
237 : {
238 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
239 : {
240 120 : auto c_pt = circle_points[m - ii];
241 120 : tr_corner_points[ii + 1] = quadrant_centers[2] + c_pt;
242 : }
243 24 : tr_corner_points[points_per_quad + 1] =
244 : Point(_side_bypass_length, -_assembly_pitch * 0.5 * shrink_factor, 0);
245 24 : tr_corner_points[points_per_quad + 2] = Point(_side_bypass_length, _side_bypass_length, 0);
246 24 : tr_corner_points[points_per_quad + 3] =
247 : Point(-_assembly_pitch * 0.5 * shrink_factor, _side_bypass_length, 0);
248 : }
249 :
250 : // Build an array of points that represent a cross section of a bottom left corner subchannel
251 : // cell. The points are ordered in this fashion:
252 : // 4 3
253 : // 2 1
254 : // 0
255 : //
256 : // 5 6
257 : std::array<Point, points_per_corner> bl_corner_points;
258 : {
259 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
260 : {
261 120 : auto c_pt = circle_points[3 * m - ii];
262 120 : bl_corner_points[ii + 1] = quadrant_centers[0] + c_pt;
263 : }
264 24 : bl_corner_points[points_per_quad + 1] =
265 : Point(-_side_bypass_length, _assembly_pitch * 0.5 * shrink_factor, 0);
266 24 : bl_corner_points[points_per_quad + 2] = Point(-_side_bypass_length, -_side_bypass_length, 0);
267 24 : bl_corner_points[points_per_quad + 3] =
268 : Point(_assembly_pitch * 0.5 * shrink_factor, -_side_bypass_length, 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 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
281 : {
282 120 : auto c_pt = circle_points[4 * m - ii];
283 120 : br_corner_points[ii + 1] = quadrant_centers[1] + c_pt;
284 : }
285 24 : br_corner_points[points_per_quad + 1] =
286 : Point(-_assembly_pitch * 0.5 * shrink_factor, -_side_bypass_length, 0);
287 24 : br_corner_points[points_per_quad + 2] = Point(_side_bypass_length, -_side_bypass_length, 0);
288 24 : br_corner_points[points_per_quad + 3] =
289 : Point(_side_bypass_length, _assembly_pitch * 0.5 * shrink_factor, 0);
290 : }
291 :
292 : // Build an array of points that represent a cross section of a top side subchannel
293 : // cell. The points are ordered in this fashion:
294 : // 8 7
295 : //
296 : // 0
297 : // 1 2 5 6
298 : // 3 4
299 : std::array<Point, points_per_side> top_points;
300 : {
301 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
302 : {
303 120 : auto c_pt = circle_points[m - ii];
304 120 : top_points[ii + 1] = quadrant_centers[2] + c_pt;
305 : }
306 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
307 : {
308 120 : auto c_pt = circle_points[2 * m - ii];
309 120 : top_points[points_per_quad + ii + 1] = quadrant_centers[3] + c_pt;
310 : }
311 24 : top_points[2 * points_per_quad + 1] =
312 : Point(_assembly_pitch * 0.5 * shrink_factor, _side_bypass_length, 0);
313 24 : top_points[2 * points_per_quad + 2] =
314 : Point(-_assembly_pitch * 0.5 * shrink_factor, _side_bypass_length, 0);
315 : }
316 :
317 : // Build an array of points that represent a cross section of a left side subchannel
318 : // cell. The points are ordered in this fashion:
319 : // 7 6
320 : // 5 4
321 : // 0
322 : // 2 3
323 : // 8 1
324 : std::array<Point, points_per_side> left_points;
325 : {
326 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
327 : {
328 120 : auto c_pt = circle_points[2 * m - ii];
329 120 : left_points[ii + 1] = quadrant_centers[3] + c_pt;
330 : }
331 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
332 : {
333 120 : auto c_pt = circle_points[3 * m - ii];
334 120 : left_points[points_per_quad + ii + 1] = quadrant_centers[0] + c_pt;
335 : }
336 24 : left_points[2 * points_per_quad + 1] =
337 : Point(-_side_bypass_length, _assembly_pitch * 0.5 * shrink_factor, 0);
338 24 : left_points[2 * points_per_quad + 2] =
339 : Point(-_side_bypass_length, -_assembly_pitch * 0.5 * shrink_factor, 0);
340 : }
341 :
342 : // Build an array of points that represent a cross section of a bottom side subchannel
343 : // cell. The points are ordered in this fashion:
344 : // 4 3
345 : // 6 5 2 1
346 : // 0
347 : //
348 : // 7 8
349 : std::array<Point, points_per_side> bottom_points;
350 : {
351 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
352 : {
353 120 : auto c_pt = circle_points[3 * m - ii];
354 120 : bottom_points[ii + 1] = quadrant_centers[0] + c_pt;
355 : }
356 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
357 : {
358 120 : auto c_pt = circle_points[4 * m - ii];
359 120 : bottom_points[points_per_quad + ii + 1] = quadrant_centers[1] + c_pt;
360 : }
361 24 : bottom_points[2 * points_per_quad + 1] =
362 : Point(-_assembly_pitch * 0.5 * shrink_factor, -_side_bypass_length, 0);
363 24 : bottom_points[2 * points_per_quad + 2] =
364 : Point(_assembly_pitch * 0.5 * shrink_factor, -_side_bypass_length, 0);
365 : }
366 :
367 : // Build an array of points that represent a cross section of a right side subchannel
368 : // cell. The points are ordered in this fashion:
369 : // 1 8
370 : // 3 2
371 : // 0
372 : // 4 5
373 : // 6 7
374 : std::array<Point, points_per_side> right_points;
375 : {
376 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
377 : {
378 120 : auto c_pt = circle_points[4 * m - ii];
379 120 : right_points[ii + 1] = quadrant_centers[1] + c_pt;
380 : }
381 144 : for (unsigned int ii = 0; ii < points_per_quad; ii++)
382 : {
383 120 : auto c_pt = circle_points[m - ii];
384 120 : right_points[points_per_quad + ii + 1] = quadrant_centers[2] + c_pt;
385 : }
386 24 : right_points[2 * points_per_quad + 1] =
387 : Point(_side_bypass_length, -_assembly_pitch * 0.5 * shrink_factor, 0);
388 24 : right_points[2 * points_per_quad + 2] =
389 : Point(_side_bypass_length, _assembly_pitch * 0.5 * shrink_factor, 0);
390 : }
391 :
392 : // Add the points to the mesh.
393 :
394 : unsigned int node_id = 0;
395 24 : Real offset_x = (_nx - 1) * _assembly_pitch / 2.0;
396 24 : Real offset_y = (_ny - 1) * _assembly_pitch / 2.0;
397 240 : for (unsigned int iy = 0; iy < _ny; iy++)
398 : {
399 216 : Point y0 = {0, _assembly_pitch * iy - offset_y, 0};
400 2880 : for (unsigned int ix = 0; ix < _nx; ix++)
401 : {
402 2664 : Point x0 = {_assembly_pitch * ix - offset_x, 0, 0};
403 2664 : if (ix == 0 && iy == 0) // Bottom Left corner
404 : {
405 648 : for (auto z : _z_grid)
406 : {
407 : Point z0{0, 0, z};
408 6240 : for (unsigned int i = 0; i < points_per_corner; i++)
409 5616 : mesh_base->add_point(bl_corner_points[i] + x0 + y0 + z0, node_id++);
410 : }
411 : }
412 2640 : else if (ix == _nx - 1 && iy == 0) // Bottom right corner
413 : {
414 648 : for (auto z : _z_grid)
415 : {
416 : Point z0{0, 0, z};
417 6240 : for (unsigned int i = 0; i < points_per_corner; i++)
418 5616 : mesh_base->add_point(br_corner_points[i] + x0 + y0 + z0, node_id++);
419 : }
420 : }
421 2616 : else if (ix == 0 && iy == _ny - 1) // top Left corner
422 : {
423 648 : for (auto z : _z_grid)
424 : {
425 : Point z0{0, 0, z};
426 6240 : for (unsigned int i = 0; i < points_per_corner; i++)
427 5616 : mesh_base->add_point(tl_corner_points[i] + x0 + y0 + z0, node_id++);
428 : }
429 : }
430 2592 : else if (ix == _nx - 1 && iy == _ny - 1) // top right corner
431 : {
432 648 : for (auto z : _z_grid)
433 : {
434 : Point z0{0, 0, z};
435 6240 : for (unsigned int i = 0; i < points_per_corner; i++)
436 5616 : mesh_base->add_point(tr_corner_points[i] + x0 + y0 + z0, node_id++);
437 : }
438 : }
439 2568 : else if (ix == 0 && (iy != _ny - 1 || iy != 0)) // left side
440 : {
441 7056 : for (auto z : _z_grid)
442 : {
443 : Point z0{0, 0, z};
444 96432 : for (unsigned int i = 0; i < points_per_side; i++)
445 89544 : mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
446 : }
447 : }
448 2400 : else if (ix == _nx - 1 && (iy != _ny - 1 || iy != 0)) // right side
449 : {
450 7056 : for (auto z : _z_grid)
451 : {
452 : Point z0{0, 0, z};
453 96432 : for (unsigned int i = 0; i < points_per_side; i++)
454 89544 : mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
455 : }
456 : }
457 2232 : else if (iy == 0 && (ix != _nx - 1 || ix != 0)) // bottom side
458 : {
459 7056 : for (auto z : _z_grid)
460 : {
461 : Point z0{0, 0, z};
462 96432 : for (unsigned int i = 0; i < points_per_side; i++)
463 89544 : mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
464 : }
465 : }
466 2064 : else if (iy == _ny - 1 && (ix != _nx - 1 || ix != 0)) // top side
467 : {
468 7056 : for (auto z : _z_grid)
469 : {
470 : Point z0{0, 0, z};
471 96432 : for (unsigned int i = 0; i < points_per_side; i++)
472 89544 : mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
473 : }
474 : }
475 : else // center
476 : {
477 93312 : for (auto z : _z_grid)
478 : {
479 : Point z0{0, 0, z};
480 2011152 : for (unsigned int i = 0; i < points_per_center; i++)
481 1919736 : mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
482 : }
483 : }
484 : }
485 : }
486 :
487 : // Add elements to the mesh. The elements are 6-node prisms. The
488 : // bases of these prisms form a triangulated representation of an cross-section
489 : // of a center subchannel.
490 : unsigned int elem_id = 0;
491 : unsigned int number_of_corner = 0;
492 : unsigned int number_of_side = 0;
493 : unsigned int number_of_center = 0;
494 : unsigned int elems_per_channel = 0;
495 : unsigned int points_per_channel = 0;
496 240 : for (unsigned int iy = 0; iy < _ny; iy++)
497 : {
498 2880 : for (unsigned int ix = 0; ix < _nx; ix++)
499 : {
500 2664 : unsigned int i_ch = _nx * iy + ix;
501 : auto subch_type = getSubchannelType(i_ch);
502 2664 : if (subch_type == EChannelType::CORNER)
503 : {
504 96 : number_of_corner++;
505 : elems_per_channel = elems_per_corner;
506 : points_per_channel = points_per_corner;
507 : }
508 2568 : else if (subch_type == EChannelType::EDGE)
509 : {
510 672 : number_of_side++;
511 : elems_per_channel = elems_per_side;
512 : points_per_channel = points_per_side;
513 : }
514 : else
515 : {
516 1896 : number_of_center++;
517 : elems_per_channel = elems_per_center;
518 : points_per_channel = points_per_center;
519 : }
520 121464 : for (unsigned int iz = 0; iz < _n_cells; iz++)
521 : {
522 118800 : unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
523 118800 : number_of_side * points_per_side * (_n_cells + 1) +
524 : number_of_center * points_per_center * (_n_cells + 1) -
525 118800 : points_per_channel * (_n_cells + 1);
526 : // index of the central node at base of cell
527 118800 : unsigned int indx1 = iz * points_per_channel + elapsed_points;
528 : // index of the central node at top of cell
529 118800 : unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
530 :
531 2250960 : for (unsigned int i = 0; i < elems_per_channel; i++)
532 : {
533 2132160 : Elem * elem = new Prism6;
534 2132160 : elem->subdomain_id() = _block_id;
535 2132160 : elem->set_id(elem_id++);
536 2132160 : elem = mesh_base->add_elem(elem);
537 :
538 2132160 : elem->set_node(0, mesh_base->node_ptr(indx1));
539 2132160 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
540 2132160 : if (i != elems_per_channel - 1)
541 2013360 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
542 : else
543 118800 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
544 :
545 2132160 : elem->set_node(3, mesh_base->node_ptr(indx2));
546 2132160 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
547 2132160 : if (i != elems_per_channel - 1)
548 2013360 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
549 : else
550 118800 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
551 :
552 2132160 : if (iz == 0)
553 46752 : boundary_info.add_side(elem, 0, 0);
554 2132160 : if (iz == _n_cells - 1)
555 46752 : boundary_info.add_side(elem, 4, 1);
556 : }
557 : }
558 : }
559 : }
560 24 : boundary_info.sideset_name(0) = "inlet";
561 24 : boundary_info.sideset_name(1) = "outlet";
562 24 : mesh_base->subdomain_name(_block_id) = name();
563 24 : mesh_base->prepare_for_use();
564 :
565 24 : return mesh_base;
566 0 : }
|
