diff --git a/src/scheme/reconstruction_utils/CMakeLists.txt b/src/scheme/reconstruction_utils/CMakeLists.txt index fa1395b2dd9615050ea557e6b83f96a3de994ca8..1d4ebcb48c3d47c9c19421151257deb8a31bf9a9 100644 --- a/src/scheme/reconstruction_utils/CMakeLists.txt +++ b/src/scheme/reconstruction_utils/CMakeLists.txt @@ -2,6 +2,7 @@ add_library(PugsSchemeReconstructionUtils BoundaryIntegralReconstructionMatrixBuilder.cpp + CellCenterReconstructionMatrixBuilder.cpp ElementIntegralReconstructionMatrixBuilder.cpp ) diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp new file mode 100644 index 0000000000000000000000000000000000000000..94d1ef9a79ac1fb8f7e3cc0ece4c8eb21672bef7 --- /dev/null +++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp @@ -0,0 +1,92 @@ +#include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp> + +#include <geometry/SymmetryUtils.hpp> +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> +#include <scheme/DiscreteFunctionDPk.hpp> + +template <MeshConcept MeshTypeT> +void +PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<MeshTypeT>::build( + const CellId cell_j_id, + ShrinkMatrixView<SmallMatrix<double>>& 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]; + } + } + } +} + +template <MeshConcept MeshTypeT> +PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<MeshTypeT>::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"); + } +} + +template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<1>>::build( + const CellId, + ShrinkMatrixView<SmallMatrix<double>>&); + +template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<2>>::build( + const CellId, + ShrinkMatrixView<SmallMatrix<double>>&); + +template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<3>>::build( + const CellId, + ShrinkMatrixView<SmallMatrix<double>>&); + +template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder< + Mesh<1>>::CellCenterReconstructionMatrixBuilder(const MeshType&, + const size_t, + const SmallArray<const Rd>&, + const SmallArray<const Rd>&, + const CellToCellStencilArray&); + +template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder< + Mesh<2>>::CellCenterReconstructionMatrixBuilder(const MeshType&, + const size_t, + const SmallArray<const Rd>&, + const SmallArray<const Rd>&, + const CellToCellStencilArray&); + +template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder< + Mesh<3>>::CellCenterReconstructionMatrixBuilder(const MeshType&, + const size_t, + const SmallArray<const Rd>&, + const SmallArray<const Rd>&, + const CellToCellStencilArray&); diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp index ae3488d8125dc4690856294494c1218902427693..a8205a625dc0d62e35a33469fe2864719f8147a6 100644 --- a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp +++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp @@ -2,9 +2,9 @@ #define CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP #include <algebra/ShrinkMatrixView.hpp> -#include <mesh/Connectivity.hpp> +#include <algebra/SmallMatrix.hpp> +#include <mesh/ItemValue.hpp> #include <mesh/StencilArray.hpp> -#include <scheme/DiscreteFunctionDPk.hpp> #include <scheme/PolynomialReconstruction.hpp> #include <utils/SmallArray.hpp> @@ -27,54 +27,13 @@ class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder 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]; - } - } - } - } + void build(const CellId cell_j_id, ShrinkMatrixView<SmallMatrix<double>>& A); 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"); - } - } + const CellToCellStencilArray& stencil_array); ~CellCenterReconstructionMatrixBuilder() = default; };