diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index e29e8b8da42f6d897ed96398be3a2a0bd6c82691..b9cfb57e151e5e8f3ba3df6b33b215e3805df345 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -824,7 +824,7 @@ PolynomialReconstruction::_build(
return this->_build<CellCenterReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::boundary: {
- if constexpr (MeshType::Dimension == 2) {
+ if constexpr (MeshType::Dimension < 3) {
return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
discrete_function_variant_list);
}
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
index 0e3c01a7d4857ad39ff3be8555141be6b47c4caf..5c7c5af3e4d39a9b903becfaeaade4a8aef4a6ca 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -9,8 +9,8 @@
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
-template <>
-class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::ConnectivityData
+template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::SpecificDimensionalData
{
private:
const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
@@ -18,8 +18,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh
const FaceValuePerCell<const bool> m_cell_face_is_reversed;
public:
- PUGS_INLINE
- const auto&
+ PUGS_INLINE const auto&
cellToFaceMatrix() const noexcept
{
return m_cell_to_face_matrix;
@@ -39,13 +38,33 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh
return m_cell_face_is_reversed;
}
- ConnectivityData(const Connectivity<2>& connectivity)
+ SpecificDimensionalData(const Connectivity<2>& connectivity)
: m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
{}
- ~ConnectivityData() = default;
+ ~SpecificDimensionalData() = default;
+};
+
+template <>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::SpecificDimensionalData
+{
+ private:
+ const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_matrix;
+
+ public:
+ PUGS_INLINE
+ const auto&
+ cellToNodeMatrix() const noexcept
+ {
+ return m_cell_to_node_matrix;
+ }
+
+ SpecificDimensionalData(const Connectivity<1>& connectivity) : m_cell_to_node_matrix{connectivity.cellToNodeMatrix()}
+ {}
+
+ ~SpecificDimensionalData() = default;
};
template <MeshConcept MeshTypeT>
@@ -70,23 +89,23 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
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;
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * velocity_perp_e1;
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] * m_antiderivative_coef[k]; // ((1. * k) / (k + 1));
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((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];
+ m_tmp_array[k++] = y_yj * m_tmp_array[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];
+ mean_of_ejk[k - 1] += m_tmp_array[k];
}
}
}
@@ -103,9 +122,9 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const double inv_Vi = 1. / m_Vj[cell_id];
- const auto& face_is_reversed = m_dimensional_data->cellFaceIsReversed()[cell_id];
- const auto& cell_face_list = m_dimensional_data->cellToFaceMatrix()[cell_id];
- const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
mean_of_ejk.fill(0);
for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
@@ -135,9 +154,9 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const double inv_Vi = 1. / m_Vj[cell_id];
- const auto& face_is_reversed = m_dimensional_data->cellFaceIsReversed()[cell_id];
- const auto& cell_face_list = m_dimensional_data->cellToFaceMatrix()[cell_id];
- const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
mean_of_ejk.fill(0);
for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
@@ -157,13 +176,87 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
}
}
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::_computeEjkMeanByBoundary(
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (m_xr[cell_node_list[1]] - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (m_xr[cell_node_list[0]] - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[0]]) - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[1]]) - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
template <MeshConcept MeshTypeT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
const CellId cell_j_id,
ShrinkMatrixView<SmallMatrix<double>>& A)
{
- if constexpr (MeshType::Dimension == 2) {
+ if constexpr (MeshType::Dimension < 3) {
const auto& stencil_cell_list = m_stencil_array[cell_j_id];
const Rd& Xj = m_xj[cell_j_id];
@@ -213,11 +306,11 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
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{m_basis_dimension},
+ m_tmp_array{m_basis_dimension},
m_mean_j_of_ejk{m_basis_dimension - 1},
m_mean_i_of_ejk{m_basis_dimension - 1},
- m_dimensional_data{std::make_shared<ConnectivityData>(mesh.connectivity())},
+ m_dimensional_data_ptr{std::make_shared<SpecificDimensionalData>(mesh.connectivity())},
m_stencil_array{stencil_array},
m_symmetry_origin_list{symmetry_origin_list},
@@ -226,20 +319,30 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
m_xr{mesh.xr()}
{
- if constexpr (MeshType::Dimension == 2) {
- SmallArray<double> antiderivative_coef(m_polynomial_degree + 1);
- for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
- antiderivative_coef[k] = ((1. * k) / (k + 1));
- }
-
- m_antiderivative_coef = antiderivative_coef;
+ SmallArray<double> antiderivative_correction_coef(m_polynomial_degree + 1);
+ for (size_t k = 0; k < antiderivative_correction_coef.size(); ++k) {
+ // The antiderivative of x^k is k/(k+1) times the antiderivative of x^(k-1)
+ antiderivative_correction_coef[k] = ((1. * k) / (k + 1));
}
+
+ m_antiderivative_correction_coef = antiderivative_correction_coef;
}
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
const CellId,
ShrinkMatrixView<SmallMatrix<double>>&);
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
Mesh<2>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
const size_t,
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index cb638fa2bef710f8ffe92257c546911b158f4b42..824aa902796cb0989c500b8b5eb05c8f0b4393fc 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -26,12 +26,12 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
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;
+ const SmallArray<double> m_tmp_array;
SmallArray<double> m_mean_j_of_ejk;
SmallArray<double> m_mean_i_of_ejk;
- class ConnectivityData;
- std::shared_ptr<ConnectivityData> m_dimensional_data;
+ class SpecificDimensionalData;
+ std::shared_ptr<SpecificDimensionalData> m_dimensional_data_ptr;
const CellToCellStencilArray& m_stencil_array;
@@ -42,7 +42,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const CellValue<const Rd> m_xj;
const NodeValue<const Rd> m_xr;
- SmallArray<const double> m_antiderivative_coef;
+ SmallArray<const double> m_antiderivative_correction_coef;
template <typename ConformTransformationT>
void _computeEjkBoundaryMean(const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
index a1c5c058975c17f9a5ec7e40c41f28b9abb18a66..db13612c1d8c5efad94a6f890251e366826d014b 100644
--- a/tests/test_PolynomialReconstruction_degree_2.cpp
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -118,7 +118,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -139,7 +139,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -172,7 +172,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -219,7 +219,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -231,7 +231,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}