LCOV - code coverage report
Current view: top level - src/meshgenerators - SCMQuadDuctMeshGenerator.C (source / functions) Hit Total Coverage
Test: idaholab/moose subchannel: #32971 (54bef8) with base c6cf66 Lines: 102 110 92.7 %
Date: 2026-05-29 20:40:47 Functions: 8 9 88.9 %
Legend: Lines: hit not hit

          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 "SCMQuadDuctMeshGenerator.h"
      11             : #include "QuadSubChannelMesh.h"
      12             : #include <cmath>
      13             : #include "libmesh/unstructured_mesh.h"
      14             : #include "libmesh/face_quad4.h"
      15             : 
      16             : registerMooseObject("SubChannelApp", SCMQuadDuctMeshGenerator);
      17             : 
      18             : InputParameters
      19          10 : SCMQuadDuctMeshGenerator::validParams()
      20             : {
      21          10 :   InputParameters params = MeshGenerator::validParams();
      22          10 :   params.addClassDescription("Creates a mesh of 2D duct cells around a square-lattice subassembly");
      23          20 :   params.addRequiredParam<MeshGeneratorName>("input", "The corresponding subchannel mesh");
      24          20 :   params.addParam<unsigned int>("block_id", 2, "Subdomain id for the duct mesh cells");
      25          20 :   params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
      26          20 :   params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
      27          20 :   params.addRequiredParam<Real>("heated_length", "Heated length [m]");
      28          20 :   params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
      29          20 :   params.addRangeCheckedParam<Real>("pitch", "pitch > 0", "Lattice pitch (must be positive)");
      30          20 :   params.addRangeCheckedParam<unsigned int>(
      31             :       "nx",
      32             :       "nx > 1",
      33             :       "Number of channels in the x direction for the subchannel assembly. Must be more than 1 to "
      34             :       "built a duct[-]");
      35          20 :   params.addRangeCheckedParam<unsigned int>(
      36             :       "ny",
      37             :       "ny > 1",
      38             :       "Number of channels in the y direction for the subchannel assembly. Must be more than 1 to "
      39             :       "built a duct[-]");
      40          20 :   params.addRequiredParam<Real>("side_gap",
      41             :                                 "Gap between duct wall and outer pin lattice: distance(edge pin "
      42             :                                 "center, duct wall) = pitch/2 + side_gap [m]");
      43          10 :   return params;
      44           0 : }
      45             : 
      46           5 : SCMQuadDuctMeshGenerator::SCMQuadDuctMeshGenerator(const InputParameters & params)
      47             :   : MeshGenerator(params),
      48           5 :     _input(getMesh("input")),
      49          10 :     _n_cells(getParam<unsigned int>("n_cells")),
      50          10 :     _unheated_length_entry(getParam<Real>("unheated_length_entry")),
      51          10 :     _heated_length(getParam<Real>("heated_length")),
      52          10 :     _unheated_length_exit(getParam<Real>("unheated_length_exit")),
      53          10 :     _block_id(getParam<unsigned int>("block_id")),
      54          10 :     _pitch(getParam<Real>("pitch")),
      55          10 :     _nx(getParam<unsigned int>("nx")),
      56          10 :     _ny(getParam<unsigned int>("ny")),
      57          15 :     _side_gap(getParam<Real>("side_gap"))
      58             : {
      59           5 :   SubChannelMesh::generateZGrid(
      60           5 :       _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);
      61           5 : }
      62             : 
      63             : std::unique_ptr<MeshBase>
      64           5 : SCMQuadDuctMeshGenerator::generate()
      65             : {
      66           5 :   std::unique_ptr<MeshBase> mesh_base = std::move(_input);
      67           5 :   mesh_base->set_mesh_dimension(3);
      68             : 
      69             :   std::vector<Point> cross_sec;
      70           5 :   ductCrossSec(cross_sec, _nx, _ny, _pitch, _side_gap);
      71             :   std::vector<Point> points;
      72           5 :   ductPoints(points, cross_sec, _z_grid);
      73             :   std::vector<std::vector<size_t>> elem_point_indices;
      74           5 :   ductElems(elem_point_indices, _z_grid.size(), cross_sec.size());
      75             :   std::vector<Node *> duct_nodes;
      76           5 :   buildDuct(mesh_base, duct_nodes, points, elem_point_indices, _block_id);
      77           5 :   mesh_base->subdomain_name(_block_id) = name();
      78             : 
      79           5 :   mesh_base->prepare_for_use();
      80             : 
      81             :   // Mirror the Tri variant: provide mapping hooks into the subchannel mesh
      82           5 :   auto & sch_mesh = static_cast<QuadSubChannelMesh &>(*_mesh);
      83           5 :   sch_mesh.setChannelToDuctMaps(duct_nodes);
      84             : 
      85           5 :   return mesh_base;
      86           5 : }
      87             : 
      88             : size_t
      89        3570 : SCMQuadDuctMeshGenerator::ductPointIndex(unsigned int points_per_layer,
      90             :                                          unsigned int layer,
      91             :                                          unsigned int point) const
      92             : {
      93        3570 :   return layer * points_per_layer + point;
      94             : }
      95             : 
      96             : void
      97           0 : SCMQuadDuctMeshGenerator::ductCorners(std::vector<Point> & corners,
      98             :                                       Real half_x,
      99             :                                       Real half_y,
     100             :                                       const Point & center) const
     101             : {
     102           0 :   corners.resize(4);
     103           0 :   corners[0] = center + Point(-half_x, -half_y, 0);
     104           0 :   corners[1] = center + Point(half_x, -half_y, 0);
     105           0 :   corners[2] = center + Point(half_x, half_y, 0);
     106           0 :   corners[3] = center + Point(-half_x, half_y, 0);
     107           0 : }
     108             : 
     109             : void
     110           5 : SCMQuadDuctMeshGenerator::ductCrossSec(std::vector<Point> & cross_sec,
     111             :                                        unsigned int nx,
     112             :                                        unsigned int ny,
     113             :                                        Real pitch,
     114             :                                        Real side_gap) const
     115             : {
     116           5 :   cross_sec.clear();
     117             : 
     118           5 :   const Real half_x = 0.5 * (nx - 1) * pitch + side_gap;
     119           5 :   const Real half_y = 0.5 * (ny - 1) * pitch + side_gap;
     120             : 
     121             :   // Subchannel-aligned grid extents (no side_gap)
     122             :   const Real x0 = -0.5 * (nx - 1) * pitch;
     123             :   const Real y0 = -0.5 * (ny - 1) * pitch;
     124             : 
     125             :   // ---- Bottom edge (y = -half_y): nx points, left -> right ----
     126          30 :   for (unsigned int i = 0; i < nx; ++i)
     127             :   {
     128             :     Real x;
     129          25 :     if (i == 0)
     130           5 :       x = -half_x;
     131          20 :     else if (i == nx - 1)
     132           5 :       x = half_x;
     133             :     else
     134          15 :       x = x0 + i * pitch; // aligned with subchannel x grid
     135             : 
     136          25 :     cross_sec.emplace_back(x, -half_y, 0.0);
     137             :   }
     138             : 
     139             :   // ---- Right edge (x = +half_x): (ny-2) points, bottom -> top (exclude corners) ----
     140          15 :   for (unsigned int j = 1; j + 1 < ny; ++j)
     141             :   {
     142          10 :     const Real y = y0 + j * pitch; // aligned with subchannel y grid
     143          10 :     cross_sec.emplace_back(half_x, y, 0.0);
     144             :   }
     145             : 
     146             :   // ---- Top edge (y = +half_y): nx points, right -> left ----
     147          30 :   for (unsigned int i = 0; i < nx; ++i)
     148             :   {
     149          25 :     const unsigned int ii = nx - 1 - i;
     150             : 
     151             :     Real x;
     152          25 :     if (ii == 0)
     153           5 :       x = -half_x;
     154          20 :     else if (ii == nx - 1)
     155           5 :       x = half_x;
     156             :     else
     157          15 :       x = x0 + ii * pitch;
     158             : 
     159          25 :     cross_sec.emplace_back(x, half_y, 0.0);
     160             :   }
     161             : 
     162             :   // ---- Left edge (x = -half_x): (ny-2) points, top -> bottom (exclude corners) ----
     163          15 :   for (unsigned int j = 1; j + 1 < ny; ++j)
     164             :   {
     165          10 :     const unsigned int jj = ny - 1 - j;
     166          10 :     const Real y = y0 + jj * pitch; // aligned with subchannel y grid
     167          10 :     cross_sec.emplace_back(-half_x, y, 0.0);
     168             :   }
     169             : 
     170             :   // Total points: 2*nx + 2*ny - 4
     171           5 : }
     172             : 
     173             : void
     174           5 : SCMQuadDuctMeshGenerator::ductPoints(std::vector<Point> & points,
     175             :                                      const std::vector<Point> & cross_sec,
     176             :                                      const std::vector<Real> & z_layers) const
     177             : {
     178           5 :   points.resize(cross_sec.size() * z_layers.size());
     179          60 :   for (size_t i = 0; i < z_layers.size(); i++)
     180         825 :     for (size_t j = 0; j < cross_sec.size(); j++)
     181         770 :       points[ductPointIndex(cross_sec.size(), i, j)] =
     182         770 :           Point(cross_sec[j](0), cross_sec[j](1), z_layers[i]);
     183           5 : }
     184             : 
     185             : void
     186           5 : SCMQuadDuctMeshGenerator::ductElems(std::vector<std::vector<size_t>> & elem_point_indices,
     187             :                                     unsigned int n_layers,
     188             :                                     unsigned int points_per_layer) const
     189             : {
     190           5 :   elem_point_indices.clear();
     191          55 :   for (unsigned int i = 0; i < n_layers - 1; i++)
     192             :   {
     193             :     const unsigned int bottom = i;
     194          50 :     const unsigned int top = i + 1;
     195         750 :     for (unsigned int j = 0; j < points_per_layer; j++)
     196             :     {
     197             :       const unsigned int left = j;
     198         700 :       const unsigned int right = (j + 1) % points_per_layer;
     199        2100 :       elem_point_indices.push_back({ductPointIndex(points_per_layer, bottom, left),
     200         700 :                                     ductPointIndex(points_per_layer, bottom, right),
     201         700 :                                     ductPointIndex(points_per_layer, top, right),
     202         700 :                                     ductPointIndex(points_per_layer, top, left)});
     203             :     }
     204             :   }
     205           5 : }
     206             : 
     207             : void
     208           5 : SCMQuadDuctMeshGenerator::buildDuct(std::unique_ptr<MeshBase> & mesh,
     209             :                                     std::vector<Node *> & duct_nodes,
     210             :                                     const std::vector<Point> & points,
     211             :                                     const std::vector<std::vector<size_t>> & elem_point_indices,
     212             :                                     SubdomainID block) const
     213             : {
     214             :   // Create mesh nodes for all duct points and keep a local index -> Node* map
     215           5 :   duct_nodes.clear();
     216           5 :   duct_nodes.reserve(points.size());
     217         775 :   for (const auto & p : points)
     218         770 :     duct_nodes.push_back(mesh->add_point(p));
     219             : 
     220             :   // Create QUAD4 surface elements using libMesh factory style
     221         705 :   for (const auto & elem_indices : elem_point_indices)
     222             :   {
     223             :     mooseAssert(elem_indices.size() == 4,
     224             :                 "Expected 4 node indices per element when building QUAD4 elements.");
     225             : 
     226         700 :     auto elem = Elem::build(ElemType::QUAD4);
     227         700 :     elem->subdomain_id() = block;
     228             : 
     229             :     // Set the 4 nodes of the QUAD4
     230        3500 :     for (unsigned int i = 0; i < 4; ++i)
     231             :     {
     232        2800 :       const auto idx = elem_indices[i];
     233             :       mooseAssert(idx < duct_nodes.size(), "Element node index out of range.");
     234        2800 :       elem->set_node(i, duct_nodes[idx]);
     235             :     }
     236             : 
     237             :     // Hand ownership to the mesh
     238         700 :     mesh->add_elem(elem.release());
     239         700 :   }
     240           5 : }

Generated by: LCOV version 1.14