SCMQuadAssemblyMeshGenerator.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "SCMQuadAssemblyMeshGenerator.h"
#include "SubChannelMesh.h"

#include "libmesh/edge_edge2.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <memory>
#include <numeric>

registerMooseObject("SubChannelApp", SCMQuadAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           SCMQuadSubChannelMeshGenerator,
                           "06/30/2027 24:00",
                           SCMQuadAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           QuadSubChannelMeshGenerator,
                           "06/30/2027 24:00",
                           SCMQuadAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           SCMQuadPinMeshGenerator,
                           "06/30/2027 24:00",
                           SCMQuadAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           QuadPinMeshGenerator,
                           "06/30/2027 24:00",
                           SCMQuadAssemblyMeshGenerator);

InputParameters
SCMQuadAssemblyMeshGenerator::validParams()
{
  InputParameters params = MeshGenerator::validParams();
  params.addClassDescription("Creates one mesh containing both 1D subchannels and 1D pins in a "
                             "square lattice arrangement");

  params.addRequiredParam<Real>("pitch", "Pitch [m]");
  params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
  params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
  params.addRequiredParam<Real>("heated_length", "Heated length [m]");
  params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");

  params.addParam<std::vector<Real>>(
      "spacer_z", {}, "Axial location of spacers/vanes/mixing vanes [m]");
  params.addParam<std::vector<Real>>(
      "spacer_k", {}, "K-loss coefficient of spacers/vanes/mixing vanes [-]");

  params.addParam<std::vector<Real>>("z_blockage",
                                     std::vector<Real>({0.0, 0.0}),
                                     "Axial location of blockage (inlet, outlet) [m]");
  params.addParam<std::vector<unsigned int>>("index_blockage",
                                             std::vector<unsigned int>({0}),
                                             "Index of subchannels affected by blockage");
  params.addParam<std::vector<Real>>("reduction_blockage",
                                     std::vector<Real>({1.0}),
                                     "Area reduction of subchannels affected by blockage");
  params.addParam<std::vector<Real>>(
      "k_blockage", std::vector<Real>({0.0}), "Form loss coefficient of blocked subchannels");

  params.addParam<Real>("Kij", 0.5, "Lateral form loss coefficient [-]");
  params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
  params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
  params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");

  params.addRequiredParam<Real>(
      "side_gap",
      "The side gap, not to be confused with the gap between pins; this refers to the gap next "
      "to the duct or else the distance between the subchannel centroid and the duct wall. "
      "distance(edge pin center, duct wall) = pitch / 2 + side_gap [m]");

  params.addParam<unsigned int>("subchannel_block_id", 0, "Subchannel block id");
  params.addParam<unsigned int>("pin_block_id", 1, "Fuel Pin block id");

  return params;
}

SCMQuadAssemblyMeshGenerator::SCMQuadAssemblyMeshGenerator(const InputParameters & params)
  : MeshGenerator(params),
    _unheated_length_entry(getParam<Real>("unheated_length_entry")),
    _heated_length(getParam<Real>("heated_length")),
    _unheated_length_exit(getParam<Real>("unheated_length_exit")),
    _spacer_z(getParam<std::vector<Real>>("spacer_z")),
    _spacer_k(getParam<std::vector<Real>>("spacer_k")),
    _z_blockage(getParam<std::vector<Real>>("z_blockage")),
    _index_blockage(getParam<std::vector<unsigned int>>("index_blockage")),
    _reduction_blockage(getParam<std::vector<Real>>("reduction_blockage")),
    _k_blockage(getParam<std::vector<Real>>("k_blockage")),
    _kij(getParam<Real>("Kij")),
    _pitch(getParam<Real>("pitch")),
    _pin_diameter(getParam<Real>("pin_diameter")),
    _n_cells(getParam<unsigned int>("n_cells")),
    _nx(getParam<unsigned int>("nx")),
    _ny(getParam<unsigned int>("ny")),
    _n_channels(0),
    _n_gaps(0),
    _n_pins(0),
    _side_gap(getParam<Real>("side_gap")),
    _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
    _pin_block_id(getParam<unsigned int>("pin_block_id"))
{
  const Real total_length = _unheated_length_entry + _heated_length + _unheated_length_exit;

  if (_n_cells == 0)
    paramError("n_cells", "The number of axial cells must be greater than zero");

  if (total_length <= 0.0)
    mooseError("Total bundle length must be greater than zero");

  if (_nx == 0 || _ny == 0)
    mooseError("The number of subchannels must be greater than zero in each direction");

  if (_nx < 2 && _ny < 2)
    mooseError("The number of subchannels cannot be less than 2 in both directions. "
               "Smallest assembly allowed is either 2X1 or 1X2.");

  _n_channels = _nx * _ny;
  _n_gaps = (_nx - 1) * _ny + (_ny - 1) * _nx;
  _n_pins = (_nx - 1) * (_ny - 1);

  if (_spacer_z.size() != _spacer_k.size())
    mooseError("Size of vector spacer_z should equal size of spacer_k");

  for (const auto spacer_z : _spacer_z)
    if (spacer_z < 0.0 || spacer_z > total_length)
      paramError("spacer_z", "Spacer locations must be between zero and total bundle length");

  if (_z_blockage.size() != 2)
    paramError("z_blockage", "Size of vector z_blockage must be 2");

  if (_z_blockage.front() > _z_blockage.back())
    paramError("z_blockage", "z_blockage inlet location must not exceed outlet location");

  if (!_index_blockage.empty() &&
      *std::max_element(_index_blockage.begin(), _index_blockage.end()) > (_n_channels - 1))
    paramError("index_blockage", "Blocked subchannel index exceeds valid subchannel range");

  if (!_reduction_blockage.empty() &&
      *std::max_element(_reduction_blockage.begin(), _reduction_blockage.end()) > 1.0)
    paramError("reduction_blockage", "Area reduction of blocked subchannels cannot exceed 1");

  if ((_index_blockage.size() > _n_channels) || (_reduction_blockage.size() > _n_channels) ||
      (_k_blockage.size() > _n_channels))
    mooseError("Sizes of blockage-related vectors cannot exceed total number of subchannels");

  if ((_index_blockage.size() != _reduction_blockage.size()) ||
      (_index_blockage.size() != _k_blockage.size()) ||
      (_reduction_blockage.size() != _k_blockage.size()))
    mooseError("index_blockage, reduction_blockage, and k_blockage must have equal size");

  SubChannelMesh::generateZGrid(
      _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);

  initializeChannelData();
}

void
SCMQuadAssemblyMeshGenerator::initializeChannelData()
{
  // Defining the total length from 3 axial sections
  const Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;

  // Defining the position of the spacer grid in the numerical solution array
  std::vector<int> spacer_cell;
  for (const auto & elem : _spacer_z)
    spacer_cell.emplace_back(std::round(elem * _n_cells / L));

  // Defining the arrays for axial resistances
  std::vector<Real> kgrid(_n_cells + 1, 0.0);
  _k_grid.resize(_n_channels, std::vector<Real>(_n_cells + 1));

  // Summing the spacer resistance to the 1D grid resistance array
  for (unsigned int index = 0; index < spacer_cell.size(); index++)
    kgrid[spacer_cell[index]] += _spacer_k[index];

  // Creating the 2D grid resistance array
  for (unsigned int i = 0; i < _n_channels; i++)
    _k_grid[i] = kgrid;

  // Add blockage resistance to the 2D grid resistance array
  const Real dz = L / _n_cells;
  for (unsigned int i = 0; i < _n_cells + 1; i++)
    if ((dz * i >= _z_blockage.front() && dz * i <= _z_blockage.back()))
    {
      unsigned int index = 0;
      for (const auto & i_ch : _index_blockage)
      {
        _k_grid[i_ch][i] += _k_blockage[index];
        index++;
      }
    }

  // Defining the size of the maps
  _gap_to_chan_map.resize(_n_gaps);
  _gap_to_pin_map.resize(_n_gaps);
  _gapnodes.resize(_n_gaps);
  _chan_to_gap_map.resize(_n_channels);
  _chan_to_pin_map.resize(_n_channels);
  _pin_to_chan_map.resize(_n_pins);
  _sign_id_crossflow_map.resize(_n_channels);
  _gij_map.resize(_n_cells + 1);
  _subchannel_position.resize(_n_channels);

  for (unsigned int i = 0; i < _n_channels; i++)
  {
    _subchannel_position[i].reserve(3);
    for (unsigned int j = 0; j < 3; j++)
      _subchannel_position.at(i).push_back(0.0);
  }

  for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
    _gij_map[iz].reserve(_n_gaps);

  // Defining the signs for positive and negative flows
  const Real positive_flow = 1.0;
  const Real negative_flow = -1.0;

  // Defining the subchannel types
  _subch_type.resize(_n_channels);
  for (unsigned int iy = 0; iy < _ny; iy++)
    for (unsigned int ix = 0; ix < _nx; ix++)
    {
      const unsigned int i_ch = _nx * iy + ix;
      const bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
                             (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
      const bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);

      // Two channels side by side (a 1x2/2x1 grid) only occurs in verification cases; treat both
      // as CENTER. The smallest physical assembly is 2x2 subchannels around a single pin.
      if (_n_channels == 2)
        _subch_type[i_ch] = EChannelType::CENTER;
      else if (_n_channels == 4)
        _subch_type[i_ch] = EChannelType::CORNER;
      else if (is_corner)
        _subch_type[i_ch] = EChannelType::CORNER;
      else if (is_edge)
        _subch_type[i_ch] = EChannelType::EDGE;
      else
        _subch_type[i_ch] = EChannelType::CENTER;
    }

  // Index the east-west gaps.
  unsigned int i_gap = 0;
  for (unsigned int iy = 0; iy < _ny; iy++)
    for (unsigned int ix = 0; ix < _nx - 1; ix++)
    {
      const unsigned int i_ch = _nx * iy + ix;
      const unsigned int j_ch = _nx * iy + (ix + 1);
      _gap_to_chan_map[i_gap] = {i_ch, j_ch};
      _chan_to_gap_map[i_ch].push_back(i_gap);
      _chan_to_gap_map[j_ch].push_back(i_gap);
      _sign_id_crossflow_map[i_ch].push_back(positive_flow);
      _sign_id_crossflow_map[j_ch].push_back(negative_flow);

      // Make a gap size map.
      if (iy == 0 || iy == _ny - 1)
        _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
      else
        _gij_map[0].push_back(_pitch - _pin_diameter);

      ++i_gap;
    }

  // Index the north-south gaps.
  for (unsigned int iy = 0; iy < _ny - 1; iy++)
    for (unsigned int ix = 0; ix < _nx; ix++)
    {
      const unsigned int i_ch = _nx * iy + ix;
      const unsigned int j_ch = _nx * (iy + 1) + ix;
      _gap_to_chan_map[i_gap] = {i_ch, j_ch};
      _chan_to_gap_map[i_ch].push_back(i_gap);
      _chan_to_gap_map[j_ch].push_back(i_gap);
      _sign_id_crossflow_map[i_ch].push_back(positive_flow);
      _sign_id_crossflow_map[j_ch].push_back(negative_flow);

      // Make a gap size map.
      if (ix == 0 || ix == _nx - 1)
        _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
      else
        _gij_map[0].push_back(_pitch - _pin_diameter);

      ++i_gap;
    }

  for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
    _gij_map[iz] = _gij_map[0];

  // Make pin to channel map.
  for (unsigned int iy = 0; iy < _ny - 1; iy++)
    for (unsigned int ix = 0; ix < _nx - 1; ix++)
    {
      const unsigned int i_pin = (_nx - 1) * iy + ix;
      const unsigned int i_chan_1 = _nx * iy + ix;
      const unsigned int i_chan_2 = _nx * (iy + 1) + ix;
      const unsigned int i_chan_3 = _nx * (iy + 1) + (ix + 1);
      const unsigned int i_chan_4 = _nx * iy + (ix + 1);

      _pin_to_chan_map[i_pin].push_back(i_chan_1);
      _pin_to_chan_map[i_pin].push_back(i_chan_2);
      _pin_to_chan_map[i_pin].push_back(i_chan_3);
      _pin_to_chan_map[i_pin].push_back(i_chan_4);
    }

  // Set the subchannel positions so that the center of the assembly is the zero point.
  for (unsigned int iy = 0; iy < _ny; iy++)
    for (unsigned int ix = 0; ix < _nx; ix++)
    {
      const unsigned int i_ch = _nx * iy + ix;
      const Real offset_x = (_nx - 1) * _pitch / 2.0;
      const Real offset_y = (_ny - 1) * _pitch / 2.0;
      _subchannel_position[i_ch][0] = _pitch * ix - offset_x;
      _subchannel_position[i_ch][1] = _pitch * iy - offset_y;
    }

  if (_n_pins > 0)
  {
    // Make channel to pin map.
    for (unsigned int iy = 0; iy < _ny; iy++) // row
      for (unsigned int ix = 0; ix < _nx; ix++)
      {
        const unsigned int i_ch = _nx * iy + ix;

        // Corners contact 1/4 of one pin.
        if (iy == 0 && ix == 0)
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
        else if (iy == _ny - 1 && ix == 0)
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
        else if (iy == 0 && ix == _nx - 1)
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
        else if (iy == _ny - 1 && ix == _nx - 1)
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
        // Sides contact 1/4 of two pins.
        else if (iy == 0)
        {
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
        }
        else if (iy == _ny - 1)
        {
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
        }
        else if (ix == 0)
        {
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
        }
        else if (ix == _nx - 1)
        {
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
        }
        // Interior channels contact 1/4 of four pins.
        else
        {
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
          _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
        }
      }

    // Make gap to pin map.
    for (unsigned int ig = 0; ig < _n_gaps; ig++)
    {
      const auto i_ch = _gap_to_chan_map[ig].first;
      const auto j_ch = _gap_to_chan_map[ig].second;
      const auto & i_pins = _chan_to_pin_map[i_ch];
      const auto & j_pins = _chan_to_pin_map[j_ch];

      // Initialize with default values.
      _gap_to_pin_map[ig] = {10000, 10000};

      for (unsigned int i : i_pins)
        for (unsigned int j : j_pins)
          if (i == j)
          {
            if (_gap_to_pin_map[ig].first == 10000)
            {
              _gap_to_pin_map[ig].first = i;
              _gap_to_pin_map[ig].second = i;
            }
            else
              _gap_to_pin_map[ig].second = i;
          }
    }
  }

  // Reduce reserved memory in the channel-to-gap map.
  for (auto & gap : _chan_to_gap_map)
    gap.shrink_to_fit();

  // Reduce reserved memory in the channel-to-pin map.
  for (auto & pin : _chan_to_pin_map)
    pin.shrink_to_fit();

  // Reduce reserved memory in the pin-to-channel map.
  for (auto & pin : _pin_to_chan_map)
    pin.shrink_to_fit();
}

void
SCMQuadAssemblyMeshGenerator::buildSubchannelMesh(MeshBase & mesh_base,
                                                  BoundaryInfo & boundary_info)
{
  mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * _ny * _nx);
  mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * _ny * _nx);

  _nodes.resize(_nx * _ny);

  const Real offset_x = (_nx - 1) * _pitch / 2.0;
  const Real offset_y = (_ny - 1) * _pitch / 2.0;

  // Add the points in the shape of a rectilinear grid. The grid is regular on the xy-plane with a
  // spacing of `pitch` between points. The grid along z is irregular to account for pin spacers.
  // Store pointers in the _nodes array so we can keep track of which points are in which channels.
  dof_id_type node_id = mesh_base.n_nodes();
  for (unsigned int iy = 0; iy < _ny; iy++)
    for (unsigned int ix = 0; ix < _nx; ix++)
    {
      const unsigned int i_ch = _nx * iy + ix;
      _nodes[i_ch].reserve(_n_cells + 1);

      for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
        _nodes[i_ch].push_back(mesh_base.add_point(
            Point(_pitch * ix - offset_x, _pitch * iy - offset_y, _z_grid[iz]), node_id++));
    }

  // Add the elements which in this case are 2-node edges that link each subchannel's nodes
  // vertically.
  dof_id_type elem_id = mesh_base.n_elem();
  for (unsigned int iy = 0; iy < _ny; iy++)
    for (unsigned int ix = 0; ix < _nx; ix++)
      for (unsigned int iz = 0; iz < _n_cells; iz++)
      {
        Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
        elem->subdomain_id() = _subchannel_block_id;
        elem->set_id(elem_id++);

        const unsigned int i_ch = _nx * iy + ix;
        elem->set_node(0, _nodes[i_ch][iz]);
        elem->set_node(1, _nodes[i_ch][iz + 1]);

        if (iz == 0)
          boundary_info.add_side(elem, 0, 0);
        if (iz == _n_cells - 1)
          boundary_info.add_side(elem, 1, 1);
      }

  mesh_base.subdomain_name(_subchannel_block_id) = "subchannel";
}

void
SCMQuadAssemblyMeshGenerator::buildPinMesh(MeshBase & mesh_base)
{
  mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * (_ny - 1) * (_nx - 1));
  mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * (_ny - 1) * (_nx - 1));

  _pin_nodes.resize((_nx - 1) * (_ny - 1));

  // Add the points in the shape of a rectilinear grid. The grid is regular on the xy-plane with a
  // spacing of `pitch` between points. The grid along z is also regular. Store pointers in the
  // _pin_nodes array so we can keep track of which points are in which pins.
  const Real offset_x = (_nx - 2) * _pitch / 2.0;
  const Real offset_y = (_ny - 2) * _pitch / 2.0;

  dof_id_type node_id = mesh_base.n_nodes();
  for (unsigned int iy = 0; iy < _ny - 1; iy++)
    for (unsigned int ix = 0; ix < _nx - 1; ix++)
    {
      const unsigned int i_pin = (_nx - 1) * iy + ix;
      _pin_nodes[i_pin].reserve(_n_cells + 1);

      for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
        _pin_nodes[i_pin].push_back(mesh_base.add_point(
            Point(_pitch * ix - offset_x, _pitch * iy - offset_y, _z_grid[iz]), node_id++));
    }

  // Add the elements which in this case are 2-node edges that link each pin's nodes vertically.
  dof_id_type elem_id = mesh_base.n_elem();
  for (unsigned int iy = 0; iy < _ny - 1; iy++)
    for (unsigned int ix = 0; ix < _nx - 1; ix++)
      for (unsigned int iz = 0; iz < _n_cells; iz++)
      {
        Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
        elem->subdomain_id() = _pin_block_id;
        elem->set_id(elem_id++);

        const unsigned int i_pin = (_nx - 1) * iy + ix;
        elem->set_node(0, _pin_nodes[i_pin][iz]);
        elem->set_node(1, _pin_nodes[i_pin][iz + 1]);
      }

  mesh_base.subdomain_name(_pin_block_id) = "fuel_pins";
}

void
SCMQuadAssemblyMeshGenerator::transferMetadata(QuadSubChannelMesh & sch_mesh)
{
  // Move the metadata into QuadSubChannelMesh.
  sch_mesh._unheated_length_entry = _unheated_length_entry;
  sch_mesh._heated_length = _heated_length;
  sch_mesh._unheated_length_exit = _unheated_length_exit;
  sch_mesh._z_grid = _z_grid;
  sch_mesh._k_grid = _k_grid;
  sch_mesh._spacer_z = _spacer_z;
  sch_mesh._spacer_k = _spacer_k;
  sch_mesh._z_blockage = _z_blockage;
  sch_mesh._index_blockage = _index_blockage;
  sch_mesh._reduction_blockage = _reduction_blockage;
  sch_mesh._kij = _kij;
  sch_mesh._pitch = _pitch;
  sch_mesh._pin_diameter = _pin_diameter;
  sch_mesh._n_cells = _n_cells;

  sch_mesh._nx = _nx;
  sch_mesh._ny = _ny;
  sch_mesh._n_channels = _n_channels;
  sch_mesh._n_gaps = _n_gaps;
  sch_mesh._n_pins = _n_pins;
  sch_mesh._side_gap = _side_gap;

  sch_mesh._nodes = _nodes;
  sch_mesh._pin_nodes = _pin_nodes;
  sch_mesh._gapnodes = _gapnodes;
  sch_mesh._gap_to_chan_map = _gap_to_chan_map;
  sch_mesh._gap_to_pin_map = _gap_to_pin_map;
  sch_mesh._chan_to_gap_map = _chan_to_gap_map;
  sch_mesh._chan_to_pin_map = _chan_to_pin_map;
  sch_mesh._pin_to_chan_map = _pin_to_chan_map;
  sch_mesh._sign_id_crossflow_map = _sign_id_crossflow_map;
  sch_mesh._gij_map = _gij_map;
  sch_mesh._subchannel_position = _subchannel_position;
  sch_mesh._subch_type = _subch_type;

  sch_mesh._duct_mesh_exist = false;
  sch_mesh._pin_mesh_exist = (_n_pins > 0);
}

std::unique_ptr<MeshBase>
SCMQuadAssemblyMeshGenerator::generate()
{
  auto mesh_base = buildMeshBaseObject();
  BoundaryInfo & boundary_info = mesh_base->get_boundary_info();

  mesh_base->set_spatial_dimension(3);

  buildSubchannelMesh(*mesh_base, boundary_info);

  if (_n_pins > 0)
    buildPinMesh(*mesh_base);

  boundary_info.sideset_name(0) = "inlet";
  boundary_info.sideset_name(1) = "outlet";
  boundary_info.nodeset_name(0) = "inlet";
  boundary_info.nodeset_name(1) = "outlet";

  mesh_base->prepare_for_use();

  auto & sch_mesh = static_cast<QuadSubChannelMesh &>(*_mesh);
  transferMetadata(sch_mesh);
  sch_mesh.computeAssemblyHydraulicParameters();

  return mesh_base;
}