BoundaryIntegralReconstructionMatrixBuilder.hpp

#ifndef BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
#define BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP

#include <algebra/ShrinkMatrixView.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/StencilArray.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>

template <MeshConcept MeshType>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
 private:
  using Rd = TinyVector<MeshType::Dimension>;

  const MeshType& m_mesh;
  const size_t m_basis_dimension;
  const size_t m_polynomial_degree;

  const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
  SmallArray<double> m_mean_j_of_ejk;
  SmallArray<double> m_mean_i_of_ejk;

  const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
  const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
  const FaceValuePerCell<const bool> m_cell_face_is_reversed;

  const CellToCellStencilArray& m_stencil_array;

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

  const CellValue<const double> m_Vj;
  const CellValue<const Rd> m_xj;
  const NodeValue<const Rd> m_xr;

  void
  _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
                          const LineTransformation<2>& T,
                          const Rd& Xj,
                          const double Vi,
                          SmallArray<double>& mean_of_ejk)
  {
    const double velocity_perp_e1 = T.velocity()[1];

    for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
      const double wq          = quadrature.weight(i_q);
      const TinyVector<1> xi_q = quadrature.point(i_q);

      const Rd X_Xj = T(xi_q) - Xj;

      const double x_xj = X_Xj[0];
      const double y_yj = X_Xj[1];

      {
        size_t k                                   = 0;
        m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi;
        for (; k <= m_polynomial_degree; ++k) {
          m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
            x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * ((1. * k) / (k + 1));
        }

        for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
          const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
          for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
            m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
          }
        }
      }

      for (size_t k = 1; k < m_basis_dimension; ++k) {
        mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
      }
    }
  }

  void
  _computeEjkMeanByBoundary(const Rd& Xj, const CellId& cell_id, SmallArray<double>& mean_of_ejk)
  {
    const auto& quadrature =
      QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));

    mean_of_ejk.fill(0);
    auto cell_face_list = m_cell_to_face_matrix[cell_id];
    for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
      bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];

      const FaceId face_id = cell_face_list[i_face];
      auto face_node_list  = m_face_to_node_matrix[face_id];
      if (is_reversed) {
        const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
        _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk);
      } else {
        const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
        _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk);
      }
    }
  }

  void
  _computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
                                           const Rd& normal,
                                           const Rd& Xj,
                                           const CellId& cell_id,
                                           SmallArray<double>& mean_of_ejk)
  {
    const auto& quadrature =
      QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));

    mean_of_ejk.fill(0);
    auto cell_face_list = m_cell_to_face_matrix[cell_id];
    for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
      bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];

      const FaceId face_id = cell_face_list[i_face];
      auto face_node_list  = m_face_to_node_matrix[face_id];

      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);

      if (is_reversed) {
        const LineTransformation<2> T{x1, x0};
        _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk);
      } else {
        const LineTransformation<2> T{x0, x1};
        _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk);
      }
    }
  }

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

      const Rd& Xj = m_xj[cell_j_id];

      _computeEjkMeanByBoundary(Xj, cell_j_id, m_mean_j_of_ejk);

      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];

        _computeEjkMeanByBoundary(Xj, cell_i_id, m_mean_i_of_ejk);

        for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
          A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[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];

          _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, Xj, cell_i_id, m_mean_i_of_ejk);

          for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
            A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
          }
        }
      }
    } else {
      throw NotImplementedError("invalid mesh dimension");
    }
  }

  BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
                                              const size_t polynomial_degree,
                                              const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
                                              const SmallArray<double>& mean_j_of_ejk,
                                              const SmallArray<double>& mean_i_of_ejk,
                                              const SmallArray<const Rd>& symmetry_origin_list,
                                              const SmallArray<const Rd>& symmetry_normal_list,
                                              const CellToCellStencilArray& stencil_array)
    : m_mesh{mesh},
      m_basis_dimension{
        DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)},
      m_polynomial_degree{polynomial_degree},
      m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{inv_Vj_alpha_p_1_wq_X_prime_orth_ek},
      m_mean_j_of_ejk{mean_j_of_ejk},
      m_mean_i_of_ejk{mean_i_of_ejk},
      m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
      m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
      m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
      m_stencil_array{stencil_array},
      m_symmetry_origin_list{symmetry_origin_list},
      m_symmetry_normal_list{symmetry_normal_list},
      m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
      m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
      m_xr{mesh.xr()}
  {}
  BoundaryIntegralReconstructionMatrixBuilder()  = default;
  ~BoundaryIntegralReconstructionMatrixBuilder() = default;
};

#endif   // BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP