CellCenterReconstructionMatrixBuilder.hpp

#ifndef CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
#define CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP

#include <algebra/ShrinkMatrixView.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/StencilArray.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>

template <MeshConcept MeshTypeT>
class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
{
 public:
  using MeshType = MeshTypeT;

  constexpr static bool handles_high_degrees = false;

 private:
  using Rd = TinyVector<MeshType::Dimension>;

  const size_t m_basis_dimension;
  const SmallArray<const Rd> m_symmetry_origin_list;
  const SmallArray<const Rd> m_symmetry_normal_list;

  const CellToCellStencilArray& m_stencil_array;
  const CellValue<const Rd> m_xj;

 public:
  template <typename MatrixType>
  void
  build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
  {
    const auto& stencil_cell_list = m_stencil_array[cell_j_id];

    const Rd& Xj = m_xj[cell_j_id];
    size_t index = 0;
    for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
      const CellId cell_i_id = stencil_cell_list[i];
      const Rd Xi_Xj         = m_xj[cell_i_id] - Xj;
      for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
        A(index, l) = Xi_Xj[l];
      }
    }
    for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
      auto& ghost_stencil  = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
      auto ghost_cell_list = ghost_stencil[cell_j_id];

      const Rd& origin = m_symmetry_origin_list[i_symmetry];
      const Rd& normal = m_symmetry_normal_list[i_symmetry];

      for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
        const CellId cell_i_id = ghost_cell_list[i];
        const Rd Xi_Xj         = symmetrize_coordinates(origin, normal, m_xj[cell_i_id]) - Xj;
        for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
          A(index, l) = Xi_Xj[l];
        }
      }
    }
  }

  CellCenterReconstructionMatrixBuilder(const MeshType& mesh,
                                        const size_t polynomial_degree,
                                        const SmallArray<const Rd>& symmetry_origin_list,
                                        const SmallArray<const Rd>& symmetry_normal_list,
                                        const CellToCellStencilArray& stencil_array)
    : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
        polynomial_degree)},
      m_symmetry_origin_list{symmetry_origin_list},
      m_symmetry_normal_list{symmetry_normal_list},
      m_stencil_array{stencil_array},
      m_xj{MeshDataManager::instance().getMeshData(mesh).xj()}
  {
    if (polynomial_degree != 1) {
      throw NormalError("cell center based reconstruction is only valid for first order");
    }
  }

  ~CellCenterReconstructionMatrixBuilder() = default;
};

#endif   // CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP