SCMDetailedTriAssemblyMeshGenerator.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 "SCMDetailedTriAssemblyMeshGenerator.h"
#include "TriSubChannelMesh.h"
#include <array>
#include <cmath>
#include <memory>
#include "libmesh/cell_prism6.h"
#include "libmesh/unstructured_mesh.h"

registerMooseObject("SubChannelApp", SCMDetailedTriAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           SCMDetailedTriSubChannelMeshGenerator,
                           "06/30/2027 24:00",
                           SCMDetailedTriAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           DetailedTriSubChannelMeshGenerator,
                           "06/30/2027 24:00",
                           SCMDetailedTriAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           SCMDetailedTriPinMeshGenerator,
                           "06/30/2027 24:00",
                           SCMDetailedTriAssemblyMeshGenerator);
registerMooseObjectRenamed("SubChannelApp",
                           DetailedTriPinMeshGenerator,
                           "06/30/2027 24:00",
                           SCMDetailedTriAssemblyMeshGenerator);

InputParameters
SCMDetailedTriAssemblyMeshGenerator::validParams()
{
  InputParameters params = MeshGenerator::validParams();
  params.addClassDescription(
      "Creates a detailed mesh of subchannels and pins in a triangular 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.addRequiredParam<unsigned int>(
      "nrings",
      "Number of fuel-pin rings per assembly, counting the center pin as the first ring [-]");
  params.addRequiredParam<Real>("flat_to_flat",
                                "Flat to flat distance for the hexagonal assembly [m]");
  params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
  params.addRangeCheckedParam<unsigned int>("num_sectors",
                                            16,
                                            "num_sectors>=4",
                                            "Number of azimuthal sectors used to discretize each "
                                            "circular pin cross section.");
  params.addParam<unsigned int>("block_id", 0, "Subchannel block id.");
  params.deprecateParam("block_id", "subchannel_block_id", "07/01/2027");
  params.addParam<unsigned int>("pin_block_id", 1, "Fuel pin block id.");
  params.addParam<bool>("verbose_flag", false, "Flag to print out the mesh coordinates.");
  return params;
}

SCMDetailedTriAssemblyMeshGenerator::SCMDetailedTriAssemblyMeshGenerator(
    const InputParameters & parameters)
  : MeshGenerator(parameters),
    _unheated_length_entry(getParam<Real>("unheated_length_entry")),
    _heated_length(getParam<Real>("heated_length")),
    _unheated_length_exit(getParam<Real>("unheated_length_exit")),
    _pitch(getParam<Real>("pitch")),
    _pin_diameter(getParam<Real>("pin_diameter")),
    _n_rings(getParam<unsigned int>("nrings")),
    _flat_to_flat(getParam<Real>("flat_to_flat")),
    _num_sectors(getParam<unsigned int>("num_sectors")),
    _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
    _pin_block_id(getParam<unsigned int>("pin_block_id")),
    _n_cells(getParam<unsigned int>("n_cells")),
    _nrods(0),
    _n_channels(0),
    _verbose(getParam<bool>("verbose_flag")),
    _elem_id(0)
{
  const Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;

  if (_n_rings < 2)
    paramError("nrings",
               "'nrings' must be at least 2. In this mesh generator, the center pin counts as "
               "the first ring, so a 7-pin bundle uses nrings = 2.");

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

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

  Real dz = L / _n_cells;
  for (unsigned int i = 0; i < _n_cells + 1; i++)
    _z_grid.push_back(dz * i);

  // x coordinate for the first position
  Real x0 = 0.0;
  // y coordinate for the first position
  Real y0 = 0.0;
  // x coordinate for the second position
  Real x1 = 0.0;
  // y coordinate for the second position dummy variable
  Real y1 = 0.0;
  // dummy variable
  Real a1 = 0.0;
  // dummy variable
  Real a2 = 0.0;
  // average x coordinate
  Real avg_coor_x = 0.0;
  // average y coordinate
  Real avg_coor_y = 0.0;
  // distance between two points
  Real dist = 0.0;
  // distance between two points
  Real dist0 = 0.0;
  // the indicator used while setting _gap_to_chan_map array
  std::vector<std::pair<unsigned int, unsigned int>> gap_fill;
  TriSubChannelMesh::pinPositions(_pin_position, _n_rings, _pitch, Point(0, 0));
  _nrods = _pin_position.size();
  // assign the pins to the corresponding rings
  unsigned int k = 0; // initialize the fuel Pin counter index
  _pins_in_rings.resize(_n_rings);
  _pins_in_rings[0].push_back(k++);
  for (unsigned int i = 1; i < _n_rings; i++)
    for (unsigned int j = 0; j < i * 6; j++)
      _pins_in_rings[i].push_back(k++);
  //  Given the number of pins and number of fuel Pin rings, the number of subchannels can be
  //  computed as follows:
  unsigned int chancount = 0;
  // Summing internal channels
  for (unsigned int j = 0; j < _n_rings - 1; j++)
    chancount += j * 6;
  // Adding external channels to the total count
  _n_channels = chancount + _nrods - 1 + (_n_rings - 1) * 6 + 6;

  // Utils for building the mesh
  _chan_to_pin_map.resize(_n_channels);
  _subch_type.resize(_n_channels);
  _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);
    }
  }

  // create the subchannels
  k = 0; // initialize the subchannel counter index
  for (unsigned int i = 1; i < _n_rings; i++)
  {
    // find the closest Pin at back ring
    for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
    {
      if (j == _pins_in_rings[i].size() - 1)
      {
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
        avg_coor_x =
            0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
        avg_coor_y =
            0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
      }
      else
      {
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
        avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
                            _pin_position[_pins_in_rings[i][j + 1]](0));
        avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
                            _pin_position[_pins_in_rings[i][j + 1]](1));
      }
      dist0 = 1.0e+5;
      _chan_to_pin_map[k].push_back(_pins_in_rings[i - 1][0]);
      for (unsigned int l = 0; l < _pins_in_rings[i - 1].size(); l++)
      {
        dist = std::sqrt(pow(_pin_position[_pins_in_rings[i - 1][l]](0) - avg_coor_x, 2) +
                         pow(_pin_position[_pins_in_rings[i - 1][l]](1) - avg_coor_y, 2));

        if (dist < dist0)
        {
          _chan_to_pin_map[k][2] = _pins_in_rings[i - 1][l];
          dist0 = dist;
        }
      }
      _subch_type[k] = EChannelType::CENTER;
      _orientation_map.insert(std::make_pair(k, 0.0));
      k = k + 1;
    }

    // find the closest Pin at front ring
    for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
    {
      if (j == _pins_in_rings[i].size() - 1)
      {
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
        avg_coor_x =
            0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
        avg_coor_y =
            0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
      }
      else
      {
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
        _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
        avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
                            _pin_position[_pins_in_rings[i][j + 1]](0));
        avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
                            _pin_position[_pins_in_rings[i][j + 1]](1));
      }
      // if the outermost ring, set the edge subchannels first... then the corner subchannels
      if (i == _n_rings - 1)
      {
        // add  edges
        _subch_type[k] = EChannelType::EDGE; // an edge subchannel is created
        k = k + 1;
        if (j % i == 0)
        {
          // corner subchannel
          _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
          _subch_type[k] = EChannelType::CORNER;
          k = k + 1;
        }
        // if not the outer most ring
      }
      else
      {
        dist0 = 1.0e+5;
        _chan_to_pin_map[k].push_back(_pins_in_rings[i + 1][0]);
        for (unsigned int l = 0; l < _pins_in_rings[i + 1].size(); l++)
        {
          dist = std::sqrt(pow(_pin_position[_pins_in_rings[i + 1][l]](0) - avg_coor_x, 2) +
                           pow(_pin_position[_pins_in_rings[i + 1][l]](1) - avg_coor_y, 2));
          if (dist < dist0)
          {
            _chan_to_pin_map[k][2] = _pins_in_rings[i + 1][l];
            dist0 = dist;
          }
        }
        _subch_type[k] = EChannelType::CENTER;
        _orientation_map.insert(std::make_pair(k, libMesh::pi));
        k = k + 1;
      }
    }
  }

  for (auto & pin : _chan_to_pin_map)
    pin.shrink_to_fit();

  // set the subchannel positions
  Real _duct_to_pin_gap =
      0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
  for (unsigned int i = 0; i < _n_channels; i++)
  {
    if (_subch_type[i] == EChannelType::CENTER)
    {
      _subchannel_position[i][0] =
          (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0) +
           _pin_position[_chan_to_pin_map[i][2]](0)) /
          3.0;
      _subchannel_position[i][1] =
          (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1) +
           _pin_position[_chan_to_pin_map[i][2]](1)) /
          3.0;
    }
    else if (_subch_type[i] == EChannelType::EDGE)
    {
      for (unsigned int j = 0; j < _n_channels; j++)
      {
        if (_subch_type[j] == EChannelType::CENTER &&
            ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
              _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1]) ||
             (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
              _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
        {
          x0 = _pin_position[_chan_to_pin_map[j][2]](0);
          y0 = _pin_position[_chan_to_pin_map[j][2]](1);
        }
        else if (_subch_type[j] == EChannelType::CENTER &&
                 ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
                   _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
                  (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
                   _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
        {
          x0 = _pin_position[_chan_to_pin_map[j][1]](0);
          y0 = _pin_position[_chan_to_pin_map[j][1]](1);
        }
        else if (_subch_type[j] == EChannelType::CENTER &&
                 ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
                   _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
                  (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
                   _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1])))
        {
          x0 = _pin_position[_chan_to_pin_map[j][0]](0);
          y0 = _pin_position[_chan_to_pin_map[j][0]](1);
        }
        x1 = 0.5 *
             (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0));
        y1 = 0.5 *
             (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1));
        a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
        a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
        _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
        _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
      } // j
    }
    else if (_subch_type[i] == EChannelType::CORNER)
    {
      x0 = _pin_position[0](0);
      y0 = _pin_position[0](1);
      x1 = _pin_position[_chan_to_pin_map[i][0]](0);
      y1 = _pin_position[_chan_to_pin_map[i][0]](1);
      a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
      a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
      _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
      _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
    }
  }
}

void
SCMDetailedTriAssemblyMeshGenerator::generatePin(std::unique_ptr<MeshBase> & mesh_base,
                                                 const Point & center)
{
  const Real dalpha = 360. / _num_sectors;
  const Real radius = _pin_diameter / 2.;

  // Add a center node and radial boundary nodes on each axial level so each pin is discretized into
  // triangular prism sectors.
  std::vector<std::vector<Node *>> nodes;
  nodes.resize(_n_cells + 1);
  for (unsigned int k = 0; k < _n_cells + 1; k++)
  {
    const Real elev = _z_grid[k];
    nodes[k].push_back(mesh_base->add_point(Point(center(0), center(1), elev)));
    Real alpha = 0.;
    for (unsigned int i = 0; i < _num_sectors; i++, alpha += dalpha)
    {
      const Real dx = radius * std::cos(alpha * M_PI / 180.);
      const Real dy = radius * std::sin(alpha * M_PI / 180.);
      nodes[k].push_back(mesh_base->add_point(Point(center(0) + dx, center(1) + dy, elev)));
    }
  }

  // Add the pin volume elements, linking matching radial sectors between adjacent axial levels.
  for (unsigned int k = 0; k < _n_cells; k++)
    for (unsigned int i = 0; i < _num_sectors; i++)
    {
      Elem * elem = mesh_base->add_elem(std::make_unique<Prism6>());
      elem->subdomain_id() = _pin_block_id;
      elem->set_id(_elem_id++);
      const unsigned int ctr_idx = 0;
      const unsigned int idx1 = (i % _num_sectors) + 1;
      const unsigned int idx2 = ((i + 1) % _num_sectors) + 1;
      elem->set_node(0, nodes[k][ctr_idx]);
      elem->set_node(1, nodes[k][idx1]);
      elem->set_node(2, nodes[k][idx2]);
      elem->set_node(3, nodes[k + 1][ctr_idx]);
      elem->set_node(4, nodes[k + 1][idx1]);
      elem->set_node(5, nodes[k + 1][idx2]);
    }
}

std::unique_ptr<MeshBase>
SCMDetailedTriAssemblyMeshGenerator::generate()
{
  auto mesh_base = buildMeshBaseObject();
  BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
  mesh_base->set_spatial_dimension(3);
  // Define the resolution (the number of points used to represent a circle).
  // This must be divisible by 4.
  const unsigned int theta_res_triangle = 24; // TODO: parameterize
  const unsigned int theta_res_square = 24;   // TODO: parameterize
  // Compute the number of points needed to represent one sixth and quadrant of a circle.
  const unsigned int points_per_sixth = theta_res_triangle / 6 + 1;
  const unsigned int points_per_quadrant = theta_res_square / 4 + 1;

  // Compute the points needed to represent one axial cross-flow of a subchannel.
  // For the center subchannel (sc) there is one center point plus the points from 3 side
  // circles.
  const unsigned int points_per_center = points_per_sixth * 3 + 1;
  // For the corner sc there is one center point plus the points from 1 side circle plus 3
  // corners
  const unsigned int points_per_corner = points_per_sixth * 1 + 1 + 3;
  // For the side sc there is one center point plus the points from 2 intersecting circles plus 2
  // corners
  const unsigned int points_per_side = points_per_quadrant * 2 + 1 + 2;

  if (_verbose)
  {
    _console << "Points per center: " << points_per_center << std::endl;
    _console << "Points per side: " << points_per_side << std::endl;
    _console << "Points per corner: " << points_per_corner << std::endl;
  }

  // Compute the number of elements (Prism6) which combined base creates the sub-channel
  // cross-section
  const unsigned int elems_per_center = theta_res_triangle * 3 / 6 + 3; // TODO: check
  const unsigned int elems_per_corner = theta_res_triangle / 6 + 4;
  const unsigned int elems_per_side = 2 * theta_res_square / 4 + 4;
  if (_verbose)
  {
    _console << "Elems per center: " << elems_per_center << std::endl;
    _console << "Elems per side: " << elems_per_side << std::endl;
    _console << "Elems per corner: " << elems_per_corner << std::endl;
  }

  // specify number and type of sub-channel
  unsigned int n_corner = 6;
  unsigned int n_side = (_n_rings - 1) * 6;
  unsigned int n_center = _n_channels - n_side - n_corner;
  if (_verbose)
  {
    _console << "Centers: " << n_center << std::endl;
    _console << "Sides: " << n_side << std::endl;
    _console << "Corners: " << n_corner << std::endl;
  }

  // Compute the total number of points and elements.
  const unsigned int points_per_level =
      n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
  const unsigned int elems_per_level =
      n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
  if (_verbose)
  {
    _console << "Points per level: " << points_per_level << std::endl;
    _console << "Elements per level: " << elems_per_level << std::endl;
  }
  const unsigned int n_points = points_per_level * (_n_cells + 1);
  const unsigned int n_elems = elems_per_level * _n_cells;
  const unsigned int n_pins = _pin_position.size();
  const unsigned int pin_points = n_pins > 0 ? (_n_cells + 1) * (_num_sectors + 1) * n_pins : 0;
  const unsigned int pin_elems = n_pins > 0 ? _n_cells * _num_sectors * n_pins : 0;
  if (_verbose)
  {
    _console << "Number of points: " << n_points << std::endl;
    _console << "Number of elements: " << n_elems << std::endl;
  }
  mesh_base->reserve_nodes(n_points + pin_points);
  mesh_base->reserve_elem(n_elems + pin_elems);
  // Build an array of points arranged in a circle on the xy-plane. (last and first node overlap)
  // We build for both the square discretization in the edges and the triangular discretization
  // within the mesh
  const Real radius = _pin_diameter / 2.0;
  std::array<Point, theta_res_square + 1> circle_points_square;
  {
    Real theta = 0;
    for (unsigned int i = 0; i < theta_res_square + 1; i++)
    {
      circle_points_square[i](0) = radius * std::cos(theta);
      circle_points_square[i](1) = radius * std::sin(theta);
      theta += 2.0 * libMesh::pi / theta_res_square;
    }
  }
  std::array<Point, theta_res_triangle + 1> circle_points_triangle;
  {
    Real theta = 0;
    for (unsigned int i = 0; i < theta_res_triangle + 1; i++)
    {
      circle_points_triangle[i](0) = radius * std::cos(theta);
      circle_points_triangle[i](1) = radius * std::sin(theta);
      theta += 2.0 * libMesh::pi / theta_res_triangle;
    }
  }
  // Define "quadrant center" reference points.  These will be the centers of
  // the 3 circles that represent the fuel pins.  These centers are
  // offset a little bit so that in the final mesh, there is a tiny gap between
  // neighboring subchannel cells.  That allows us to easily map a solution to
  // this detailed mesh with a nearest-neighbor search.
  const Real shrink_factor = 0.99999;
  // Quadrants are used only for the side and corner subchannels
  Real _duct_to_pin_gap =
      0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
  std::array<Point, 2> quadrant_centers_sides;
  quadrant_centers_sides[0] = Point(
      -_pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
  quadrant_centers_sides[1] = Point(
      _pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
  std::array<Point, 1> quadrant_centers_corner;
  quadrant_centers_corner[0] =
      Point(-(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::sin(libMesh::pi / 6) * shrink_factor,
            -(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::cos(libMesh::pi / 6) * shrink_factor,
            0);
  // Triangles are used for all center subchannels
  std::array<Point, 3> triangle_centers;
  triangle_centers[0] = Point(0, _pitch * std::cos(libMesh::pi / 6) * 2 / 3 * shrink_factor, 0);
  triangle_centers[1] = Point(-_pitch * 0.5 * shrink_factor,
                              -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor,
                              0);
  triangle_centers[2] = Point(
      _pitch * 0.5 * shrink_factor, -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor, 0);

  const unsigned int m_sixth = theta_res_triangle / 6;
  const unsigned int m_quarter = theta_res_square / 4;
  // Build an array of points that represent a cross section of a center subchannel
  // cell.  The points are ordered in this fashion:
  //    3   1
  //      2
  //      0
  // 4 5         8 9
  //  6         7
  std::array<Point, points_per_center> center_points{};
  center_points[0] = Point(0, 0, 0);
  {
    unsigned int start;
    for (unsigned int i = 0; i < 3; i++)
    {
      if (i == 0)
        start = 5 * (m_sixth);
      if (i == 1)
        start = 1 * (m_sixth);
      if (i == 2)
        start = 3 * (m_sixth);
      for (unsigned int ii = 0; ii < points_per_sixth; ii++)
      {
        auto c_pt = circle_points_triangle[start - ii];
        center_points[i * points_per_sixth + ii + 1] = triangle_centers[i] + c_pt;
      }
    }
  }

  // Build an array of points that represent a cross section of a top left corner subchannel
  // cell. The points are ordered in this fashion (x direction towards point 4, y direction towards
  // point 6):
  //    5       6
  //
  //      0
  // 4        2 3
  //        1
  std::array<Point, points_per_corner> corner_points{};
  corner_points[0] = Point(0, 0, 0);
  {
    for (unsigned int ii = 0; ii < points_per_sixth; ii++)
    {
      auto c_pt = circle_points_triangle[1 * m_quarter - ii];
      corner_points[ii + 1] = quadrant_centers_corner[0] + c_pt;
    }
    Real side_short = (_duct_to_pin_gap + _pin_diameter) * 0.5;
    Real side_long = (2.0 * _duct_to_pin_gap + _pin_diameter) * 0.5;
    Real side_length = std::sqrt(std::pow(side_short, 2) + std::pow(side_long, 2) -
                                 2 * side_short * side_long * std::cos(libMesh::pi / 6));
    Real angle =
        libMesh::pi - libMesh::pi / 3 -
        std::acos((-std::pow(side_long, 2) + std::pow(side_short, 2) + std::pow(side_length, 2)) /
                  (2 * side_short * side_length));
    corner_points[points_per_sixth + 1] = Point(side_length * std::cos(angle) * shrink_factor,
                                                -side_length * std::sin(angle) * shrink_factor,
                                                0);
    corner_points[points_per_sixth + 2] =
        Point(side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor *
                  std::tan(libMesh::pi / 6),
              side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
              0);
    corner_points[points_per_sixth + 3] =
        Point(-side_length * std::cos(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
              side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
              0);
  }

  // Build an array of points that represent a cross-section of a top side subchannel
  // cell.  The points are ordered in this fashion:
  // 8            7
  //
  //       0
  // 1 2        5 6
  //    3     4
  std::array<Point, points_per_side> side_points{};
  side_points[0] = Point(0, 0, 0);
  {
    for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
    {
      auto c_pt = circle_points_square[m_quarter - ii];
      side_points[ii + 1] = quadrant_centers_sides[0] + c_pt;
    }
    for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
    {
      auto c_pt = circle_points_square[2 * m_quarter - ii];
      side_points[points_per_quadrant + ii + 1] = quadrant_centers_sides[1] + c_pt;
    }
    side_points[2 * points_per_quadrant + 1] =
        Point(_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
    side_points[2 * points_per_quadrant + 2] =
        Point(-_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
  }

  int point_counter = 0;
  unsigned int node_id = 0;
  for (unsigned int i = 0; i < _n_channels; i++)
  {
    //    Real offset_x = _flat_to_flat / 2.0;
    //    Real offset_y = _flat_to_flat / 2.0;

    if (getSubchannelType(i) == EChannelType::CENTER)
    {
      for (auto z : _z_grid)
      {
        // Get height
        Point z0{0, 0, z};

        // Get suchannel position and assign to point
        auto loc_position = getSubchannelPosition(i);
        Point p0{loc_position[0], loc_position[1], 0};

        // Determine orientation of current subchannel
        auto subchannel_pins = getSubChannelPins(i);
        Point subchannel_side =
            getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
        Point base_center_orientation = {0, -1};

        // Get rotation angle for current subchannel
        Real dot_prod = 0;
        for (unsigned int lp = 0; lp < 2; lp++)
          dot_prod += base_center_orientation(lp) * subchannel_side(lp);
        auto theta =
            std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
        if (subchannel_side(0) < 0)
          theta = 2.0 * libMesh::pi - theta;

        //        Real distance_side = subchannel_side.norm();
        //        Real distance_top = getPinPosition(subchannel_pins[2]).norm();
        //        if (distance_top > distance_side)
        //                  theta += libMesh::pi * 0.0;

        theta += _orientation_map[i];

        theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;

        if (_verbose)
        {
          if (z == 0)
          {
            _console << "Subchannel Position: " << p0 << std::endl;
            auto pins = getSubChannelPins(i);
            for (auto r : pins)
              _console << r << " ";
            _console << std::endl;
            _console << "Theta: " << theta / libMesh::pi * 180. << std::endl;
          }
        }

        // Assigning points for center channels
        for (unsigned int i = 0; i < points_per_center; i++)
        {
          auto new_point = rotatePoint(center_points[i], theta) + p0 + z0;
          if (_verbose)
          {
            if (z == 0)
              _console << i << " - " << new_point << std::endl;
          }
          mesh_base->add_point(new_point, node_id++);
          point_counter += 1;
        }
      }
    }
    else if (getSubchannelType(i) == EChannelType::EDGE)
    {
      for (auto z : _z_grid)
      {
        // Get height
        Point z0{0, 0, z};

        // Get suchannel position and assign to point
        auto loc_position = getSubchannelPosition(i);
        Point p0{loc_position[0], loc_position[1], 0};

        // Determine orientation of current subchannel
        auto subchannel_pins = getSubChannelPins(i);
        Point subchannel_side =
            getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
        Point base_center_orientation = {0, 1};

        // Get rotation angle for current subchannel
        Real dot_prod = 0;
        for (unsigned int lp = 0; lp < 2; lp++)
          dot_prod += base_center_orientation(lp) * subchannel_side(lp);
        auto theta =
            std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
        if (subchannel_side(0) > 0)
          theta = 2. * libMesh::pi - theta;
        theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;

        if (_verbose)
        {
          if (z == 0)
          {
            _console << "Subchannel Position: " << p0 << std::endl;
            auto pins = getSubChannelPins(i);
            for (auto r : pins)
              _console << r << " ";
            _console << std::endl;
            _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
          }
        }

        // Assigning points for center channels
        for (unsigned int i = 0; i < points_per_side; i++)
        {
          auto new_point = rotatePoint(side_points[i], theta) + p0 + z0;
          if (_verbose)
          {
            if (z == 0)
              _console << i << " - " << new_point << std::endl;
          }
          mesh_base->add_point(new_point, node_id++);
          point_counter += 1;
        }
      }
    }
    else // getSubchannelType(i) == EChannelType::CORNER
    {
      for (auto z : _z_grid)
      {
        // Get height
        Point z0{0, 0, z};

        // Get suchannel position and assign to point
        auto loc_position = getSubchannelPosition(i);
        Point p0{loc_position[0], loc_position[1], 0};

        // Determine orientation of current subchannel
        auto subchannel_pins = getSubChannelPins(i);
        Point subchannel_side = getPinPosition(subchannel_pins[0]);
        Point base_center_orientation = {1, 1};

        // Get rotation angle for current subchannel
        Real dot_prod = 0;
        for (unsigned int lp = 0; lp < 2; lp++)
          dot_prod += base_center_orientation(lp) * subchannel_side(lp);
        auto theta =
            std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
        if (subchannel_side(0) > 0)
          theta = 2. * libMesh::pi - theta;
        theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;

        if (_verbose)
        {
          if (z == 0)
          {
            _console << "Subchannel Position: " << p0 << std::endl;
            auto pins = getSubChannelPins(i);
            for (auto r : pins)
              _console << r << " ";
            _console << std::endl;
            _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
          }
        }

        // Assigning points for center channels
        for (unsigned int i = 0; i < points_per_corner; i++)
        {
          auto new_point = rotatePoint(corner_points[i], theta) + p0 + z0;
          if (_verbose)
          {
            if (z == 0)
              _console << i << " - " << new_point << std::endl;
          }
          mesh_base->add_point(new_point, node_id++);
          point_counter += 1;
        }
      }
    }
  } // i
  if (_verbose)
    _console << "Point counter: " << point_counter << std::endl;

  int element_counter = 0;
  unsigned int elem_id = 0;
  unsigned int number_of_corner = 0;
  unsigned int number_of_side = 0;
  unsigned int number_of_center = 0;
  unsigned int elems_per_channel = 0;
  unsigned int points_per_channel = 0;
  for (unsigned int i = 0; i < _n_channels; i++)
  {
    auto subch_type = getSubchannelType(i);
    if (subch_type == EChannelType::CORNER)
    {
      number_of_corner++;
      elems_per_channel = elems_per_corner;
      points_per_channel = points_per_corner;
      if (_verbose)
        _console << "Corner" << std::endl;
    }
    else if (subch_type == EChannelType::EDGE)
    {
      number_of_side++;
      elems_per_channel = elems_per_side;
      points_per_channel = points_per_side;
      if (_verbose)
        _console << "Edge" << std::endl;
    }
    else if (subch_type == EChannelType::CENTER)
    {
      number_of_center++;
      elems_per_channel = elems_per_center;
      points_per_channel = points_per_center;
      if (_verbose)
        _console << "Center" << std::endl;
    }
    for (unsigned int iz = 0; iz < _n_cells; iz++)
    {
      unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
                                    number_of_side * points_per_side * (_n_cells + 1) +
                                    number_of_center * points_per_center * (_n_cells + 1) -
                                    points_per_channel * (_n_cells + 1);
      // index of the central node at base of cell
      unsigned int indx1 = iz * points_per_channel + elapsed_points;
      // index of the central node at top of cell
      unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;

      for (unsigned int i = 0; i < elems_per_channel; i++)
      {
        Elem * elem = mesh_base->add_elem(std::make_unique<Prism6>());
        elem->subdomain_id() = _subchannel_block_id;
        elem->set_id(elem_id++);

        if (_verbose)
          _console << "Node 0: " << *mesh_base->node_ptr(indx1) << std::endl;

        elem->set_node(0, mesh_base->node_ptr(indx1));
        elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
        if (i != elems_per_channel - 1)
          elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
        else
          elem->set_node(2, mesh_base->node_ptr(indx1 + 1));

        elem->set_node(3, mesh_base->node_ptr(indx2));
        elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
        if (i != elems_per_channel - 1)
          elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
        else
          elem->set_node(5, mesh_base->node_ptr(indx2 + 1));

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

        element_counter += 1;
      }
    }
  }
  if (_verbose)
    _console << "Element counter: " << element_counter << std::endl;
  boundary_info.sideset_name(0) = "inlet";
  boundary_info.sideset_name(1) = "outlet";
  mesh_base->subdomain_name(_subchannel_block_id) = "subchannel";
  if (n_pins > 0)
  {
    _elem_id = mesh_base->n_elem();
    for (auto & ctr : _pin_position)
      generatePin(mesh_base, ctr);
  }
  if (n_pins > 0)
    mesh_base->subdomain_name(_pin_block_id) = "fuel_pins";
  if (_verbose)
    _console << "Mesh assembly done" << std::endl;
  mesh_base->prepare_for_use();

  return mesh_base;
}

Point
SCMDetailedTriAssemblyMeshGenerator::rotatePoint(Point b, Real theta)
{
  std::vector<std::vector<Real>> A(3, std::vector<Real>(3));

  A[0] = {std::cos(theta), -std::sin(theta), 0.0};
  A[1] = {std::sin(theta), std::cos(theta), 0.0};
  A[2] = {0.0, 0.0, 1.0};

  Point rotated_vector = Point(0.0, 0.0, 0.0);
  for (unsigned int i = 0; i < 3; i++)
    for (unsigned int j = 0; j < 3; j++)
      rotated_vector(i) += A[i][j] * b(j);

  return rotated_vector;
}