Skip to content
Snippets Groups Projects
Commit 4926b9a0 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Instanciate relevant ReconstructionMatrixBuilderType at compile time

parent dc5cccb4
No related branches found
No related tags found
1 merge request!205High-order polynomial reconstruction
Show whitespace changes
Inline Side-by-side
Showing
with 157 additions and 99 deletions
src/scheme/PolynomialReconstruction.cpp
+ 91
− 83
View file @ 4926b9a0
...@@ -370,12 +370,13 @@ class PolynomialReconstruction::Internal ...@@ -370,12 +370,13 @@ class PolynomialReconstruction::Internal
} }
} }
template <typename ReconstructionMatrixBuilderType>
static void static void
populateDiscreteFunctionDPkByCell( populateDiscreteFunctionDPkByCell(
const CellId& cell_j_id, const CellId& cell_j_id,
const size_t& degree, const size_t& degree,
const SmallMatrix<double>& X, const SmallMatrix<double>& X,
const SmallArray<double>& mean_j_of_ejk, const ReconstructionMatrixBuilderType& reconstruction_matrix_builder,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list, const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
const std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>& const std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>&
mutable_discrete_function_dpk_variant_list) mutable_discrete_function_dpk_variant_list)
...@@ -400,11 +401,14 @@ class PolynomialReconstruction::Internal ...@@ -400,11 +401,14 @@ class PolynomialReconstruction::Internal
dpk_j[0] = p0_function[cell_j_id]; dpk_j[0] = p0_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) { if constexpr (std::is_arithmetic_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
if (degree > 1) { if (degree > 1) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i]; dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
} }
} }
}
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1]; auto& dpk_j_ip1 = dpk_j[i + 1];
...@@ -412,7 +416,9 @@ class PolynomialReconstruction::Internal ...@@ -412,7 +416,9 @@ class PolynomialReconstruction::Internal
} }
++column_begin; ++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) { } else if constexpr (is_tiny_vector_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
if (degree > 1) { if (degree > 1) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_0 = dpk_j[0]; auto& dpk_j_0 = dpk_j[0];
for (size_t k = 0; k < DataType::Dimension; ++k) { for (size_t k = 0; k < DataType::Dimension; ++k) {
...@@ -420,6 +426,7 @@ class PolynomialReconstruction::Internal ...@@ -420,6 +426,7 @@ class PolynomialReconstruction::Internal
} }
} }
} }
}
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1]; auto& dpk_j_ip1 = dpk_j[i + 1];
...@@ -429,7 +436,9 @@ class PolynomialReconstruction::Internal ...@@ -429,7 +436,9 @@ class PolynomialReconstruction::Internal
} }
column_begin += DataType::Dimension; column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) { } else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
if (degree > 1) { if (degree > 1) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_0 = dpk_j[0]; auto& dpk_j_0 = dpk_j[0];
for (size_t k = 0; k < DataType::NumberOfRows; ++k) { for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
...@@ -439,6 +448,7 @@ class PolynomialReconstruction::Internal ...@@ -439,6 +448,7 @@ class PolynomialReconstruction::Internal
} }
} }
} }
}
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1]; auto& dpk_j_ip1 = dpk_j[i + 1];
...@@ -469,11 +479,14 @@ class PolynomialReconstruction::Internal ...@@ -469,11 +479,14 @@ class PolynomialReconstruction::Internal
const size_t component_begin = l * size; const size_t component_begin = l * size;
dpk_j[component_begin] = cell_vector[l]; dpk_j[component_begin] = cell_vector[l];
if constexpr (std::is_arithmetic_v<DataType>) { if constexpr (std::is_arithmetic_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
if (degree > 1) { if (degree > 1) {
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i]; dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
} }
} }
}
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1]; auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
...@@ -481,7 +494,9 @@ class PolynomialReconstruction::Internal ...@@ -481,7 +494,9 @@ class PolynomialReconstruction::Internal
} }
++column_begin; ++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) { } else if constexpr (is_tiny_vector_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
if (degree > 1) { if (degree > 1) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_0 = dpk_j[component_begin]; auto& dpk_j_0 = dpk_j[component_begin];
for (size_t k = 0; k < DataType::Dimension; ++k) { for (size_t k = 0; k < DataType::Dimension; ++k) {
...@@ -489,6 +504,7 @@ class PolynomialReconstruction::Internal ...@@ -489,6 +504,7 @@ class PolynomialReconstruction::Internal
} }
} }
} }
}
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1]; auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
...@@ -498,12 +514,16 @@ class PolynomialReconstruction::Internal ...@@ -498,12 +514,16 @@ class PolynomialReconstruction::Internal
} }
column_begin += DataType::Dimension; column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) { } else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
if (degree > 1) { if (degree > 1) {
const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
for (size_t i = 0; i < X.numberOfRows(); ++i) { for (size_t i = 0; i < X.numberOfRows(); ++i) {
auto& dpk_j_0 = dpk_j[component_begin]; auto& dpk_j_0 = dpk_j[component_begin];
for (size_t p = 0; p < DataType::NumberOfRows; ++p) { for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) { for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
dpk_j_0(p, q) -= X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i]; dpk_j_0(p, q) -=
X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
}
} }
} }
} }
...@@ -660,12 +680,14 @@ PolynomialReconstruction::_checkDataAndSymmetriesCompatibility( ...@@ -660,12 +680,14 @@ PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
} }
} }
template <MeshConcept MeshType> template <typename ReconstructionMatrixBuilderType, MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build( PolynomialReconstruction::_build(
const std::shared_ptr<const MeshType>& p_mesh, const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{ {
static_assert(std::is_same_v<MeshType, typename ReconstructionMatrixBuilderType::MeshType>);
const MeshType& mesh = *p_mesh; const MeshType& mesh = *p_mesh;
using Rd = TinyVector<MeshType::Dimension>; using Rd = TinyVector<MeshType::Dimension>;
...@@ -748,46 +770,20 @@ PolynomialReconstruction::_build( ...@@ -748,46 +770,20 @@ PolynomialReconstruction::_build(
SmallArray<SmallVector<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency()); SmallArray<SmallVector<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency()); SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
for (size_t i = 0; i < A_pool.size(); ++i) { for (size_t i = 0; i < A_pool.size(); ++i) {
A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1); A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns); B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
G_pool[i] = SmallVector<double>(basis_dimension - 1); G_pool[i] = SmallVector<double>(basis_dimension - 1);
X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns); X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns);
mean_j_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
}
std::shared_ptr p_cell_center_reconstruction_matrix_builder =
std::make_shared<CellCenterReconstructionMatrixBuilder<MeshType>>(symmetry_origin_list, symmetry_normal_list,
stencil_array, xj);
SmallArray<std::shared_ptr<ElementIntegralReconstructionMatrixBuilder<MeshType>>>
element_integral_reconstruction_matrix_builder_pool{A_pool.size()};
for (size_t t = 0; t < element_integral_reconstruction_matrix_builder_pool.size(); ++t) {
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
element_integral_reconstruction_matrix_builder_pool[t] =
std::make_shared<ElementIntegralReconstructionMatrixBuilder<MeshType>>(*p_mesh, m_descriptor.degree(),
mean_j_of_ejk, symmetry_origin_list,
symmetry_normal_list, stencil_array);
} }
SmallArray<std::shared_ptr<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>> SmallArray<std::shared_ptr<ReconstructionMatrixBuilderType>> reconstruction_matrix_builder_pool(A_pool.size());
boundary_integral_reconstruction_matrix_builder_pool{A_pool.size()};
if constexpr (MeshType::Dimension == 2) { for (size_t t = 0; t < reconstruction_matrix_builder_pool.size(); ++t) {
for (size_t t = 0; t < boundary_integral_reconstruction_matrix_builder_pool.size(); ++t) { reconstruction_matrix_builder_pool[t] =
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t]; std::make_shared<ReconstructionMatrixBuilderType>(*p_mesh, m_descriptor.degree(), symmetry_origin_list,
boundary_integral_reconstruction_matrix_builder_pool[t] =
std::make_shared<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(*p_mesh, m_descriptor.degree(),
mean_j_of_ejk, symmetry_origin_list,
symmetry_normal_list, stencil_array); symmetry_normal_list, stencil_array);
} }
}
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) { mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
...@@ -799,25 +795,8 @@ PolynomialReconstruction::_build( ...@@ -799,25 +795,8 @@ PolynomialReconstruction::_build(
Internal<MeshType>::buildB(cell_j_id, stencil_array, discrete_function_variant_list, symmetry_normal_list, B); Internal<MeshType>::buildB(cell_j_id, stencil_array, discrete_function_variant_list, symmetry_normal_list, B);
if ((m_descriptor.degree() == 1) and ReconstructionMatrixBuilderType& reconstruction_matrix_builder = *reconstruction_matrix_builder_pool[t];
(m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) { reconstruction_matrix_builder.build(cell_j_id, A);
p_cell_center_reconstruction_matrix_builder->build(cell_j_id, A);
} else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
(m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
#warning implement boundary based reconstruction for 1d and 3d
if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
(MeshType::Dimension == 2)) {
if constexpr (MeshType::Dimension == 2) {
boundary_integral_reconstruction_matrix_builder_pool[t]->build(cell_j_id, A);
} else {
throw NotImplementedError("invalid mesh dimension");
}
} else {
element_integral_reconstruction_matrix_builder_pool[t]->build(cell_j_id, A);
}
} else {
throw UnexpectedError("invalid integration strategy");
}
if (m_descriptor.rowWeighting()) { if (m_descriptor.rowWeighting()) {
Internal<MeshType>::rowWeighting(cell_j_id, stencil_array, xj, symmetry_origin_list, symmetry_normal_list, A, Internal<MeshType>::rowWeighting(cell_j_id, stencil_array, xj, symmetry_origin_list, symmetry_normal_list, A,
...@@ -835,8 +814,9 @@ PolynomialReconstruction::_build( ...@@ -835,8 +814,9 @@ PolynomialReconstruction::_build(
Givens::solveCollectionInPlace(A, X, B); Givens::solveCollectionInPlace(A, X, B);
} }
Internal<MeshType>::populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X, Internal<MeshType>::template populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X,
mean_j_of_ejk_pool[t], discrete_function_variant_list, reconstruction_matrix_builder,
discrete_function_variant_list,
mutable_discrete_function_dpk_variant_list); mutable_discrete_function_dpk_variant_list);
tokens.release(t); tokens.release(t);
...@@ -858,6 +838,34 @@ PolynomialReconstruction::_build( ...@@ -858,6 +838,34 @@ PolynomialReconstruction::_build(
return discrete_function_dpk_variant_list; return discrete_function_dpk_variant_list;
} }
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
switch (m_descriptor.integrationMethodType()) {
case IntegrationMethodType::cell_center: {
return this->_build<CellCenterReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::boundary: {
if constexpr (MeshType::Dimension == 2) {
return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
discrete_function_variant_list);
}
[[fallthrough]];
}
case IntegrationMethodType::element: {
return this->_build<ElementIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
// LCOV_EXCL_START
default: {
throw UnexpectedError("invalid reconstruction matrix builder type");
}
// LCOV_EXCL_STOP
}
}
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::build( PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
......
src/scheme/PolynomialReconstruction.hpp
+ 5
− 0
View file @ 4926b9a0
...@@ -38,6 +38,11 @@ class PolynomialReconstruction ...@@ -38,6 +38,11 @@ class PolynomialReconstruction
const std::shared_ptr<const MeshType>& p_mesh, const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const; const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
template <typename ReconstructionMatrixBuilderType, MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
template <MeshConcept MeshType> template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build( [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
const std::shared_ptr<const MeshType>& p_mesh, const std::shared_ptr<const MeshType>& p_mesh,
......
src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+ 14
− 3
View file @ 4926b9a0
...@@ -14,9 +14,14 @@ ...@@ -14,9 +14,14 @@
#include <scheme/PolynomialReconstruction.hpp> #include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp> #include <utils/SmallArray.hpp>
template <MeshConcept MeshType> template <MeshConcept MeshTypeT>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{ {
public:
using MeshType = MeshTypeT;
constexpr static bool handles_high_degrees = true;
private: private:
using Rd = TinyVector<MeshType::Dimension>; using Rd = TinyVector<MeshType::Dimension>;
...@@ -142,6 +147,13 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder ...@@ -142,6 +147,13 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
} }
public: public:
PUGS_INLINE
SmallArray<const double>
meanjOfEjk() const
{
return m_mean_j_of_ejk;
}
template <typename MatrixType> template <typename MatrixType>
void void
build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A) build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
...@@ -188,7 +200,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder ...@@ -188,7 +200,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh, BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
const size_t polynomial_degree, const size_t polynomial_degree,
const SmallArray<double>& mean_j_of_ejk,
const SmallArray<const Rd>& symmetry_origin_list, const SmallArray<const Rd>& symmetry_origin_list,
const SmallArray<const Rd>& symmetry_normal_list, const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array) const CellToCellStencilArray& stencil_array)
...@@ -197,7 +208,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder ...@@ -197,7 +208,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)}, DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)},
m_polynomial_degree{polynomial_degree}, m_polynomial_degree{polynomial_degree},
m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{m_basis_dimension}, m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{m_basis_dimension},
m_mean_j_of_ejk{mean_j_of_ejk}, m_mean_j_of_ejk{m_basis_dimension - 1},
m_mean_i_of_ejk{m_basis_dimension - 1}, m_mean_i_of_ejk{m_basis_dimension - 1},
m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()}, m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()}, m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
......
src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
+ 18
− 7
View file @ 4926b9a0
...@@ -8,9 +8,14 @@ ...@@ -8,9 +8,14 @@
#include <scheme/PolynomialReconstruction.hpp> #include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp> #include <utils/SmallArray.hpp>
template <MeshConcept MeshType> template <MeshConcept MeshTypeT>
class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
{ {
public:
using MeshType = MeshTypeT;
constexpr static bool handles_high_degrees = false;
private: private:
using Rd = TinyVector<MeshType::Dimension>; using Rd = TinyVector<MeshType::Dimension>;
...@@ -54,16 +59,22 @@ class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder ...@@ -54,16 +59,22 @@ class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
} }
} }
CellCenterReconstructionMatrixBuilder(const SmallArray<const Rd>& symmetry_origin_list, CellCenterReconstructionMatrixBuilder(const MeshType& mesh,
const size_t polynomial_degree,
const SmallArray<const Rd>& symmetry_origin_list,
const SmallArray<const Rd>& symmetry_normal_list, const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array, const CellToCellStencilArray& stencil_array)
const CellValue<const Rd>& xj) : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
: m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(1)}, polynomial_degree)},
m_symmetry_origin_list{symmetry_origin_list}, m_symmetry_origin_list{symmetry_origin_list},
m_symmetry_normal_list{symmetry_normal_list}, m_symmetry_normal_list{symmetry_normal_list},
m_stencil_array{stencil_array}, m_stencil_array{stencil_array},
m_xj{xj} 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; ~CellCenterReconstructionMatrixBuilder() = default;
}; };
......
src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+ 16
− 6
View file @ 4926b9a0
...@@ -21,9 +21,14 @@ ...@@ -21,9 +21,14 @@
#include <scheme/PolynomialReconstruction.hpp> #include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp> #include <utils/SmallArray.hpp>
template <MeshConcept MeshType> template <MeshConcept MeshTypeT>
class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
{ {
public:
using MeshType = MeshTypeT;
constexpr static bool handles_high_degrees = true;
private: private:
using Rd = TinyVector<MeshType::Dimension>; using Rd = TinyVector<MeshType::Dimension>;
...@@ -460,6 +465,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder ...@@ -460,6 +465,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
} }
public: public:
PUGS_INLINE
SmallArray<const double>
meanjOfEjk() const
{
return m_mean_j_of_ejk;
}
template <typename MatrixType> template <typename MatrixType>
void void
build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A) build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
...@@ -502,7 +514,6 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder ...@@ -502,7 +514,6 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh, ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
const size_t polynomial_degree, const size_t polynomial_degree,
const SmallArray<double>& mean_j_of_ejk,
const SmallArray<const Rd>& symmetry_origin_list, const SmallArray<const Rd>& symmetry_origin_list,
const SmallArray<const Rd>& symmetry_normal_list, const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array) const CellToCellStencilArray& stencil_array)
...@@ -511,7 +522,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder ...@@ -511,7 +522,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
m_polynomial_degree{polynomial_degree}, m_polynomial_degree{polynomial_degree},
m_wq_detJ_ek{m_basis_dimension}, m_wq_detJ_ek{m_basis_dimension},
m_mean_j_of_ejk{mean_j_of_ejk}, m_mean_j_of_ejk{m_basis_dimension - 1},
m_mean_i_of_ejk{m_basis_dimension - 1}, m_mean_i_of_ejk{m_basis_dimension - 1},
m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()}, m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
...@@ -524,14 +535,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder ...@@ -524,14 +535,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
m_xr{mesh.xr()} m_xr{mesh.xr()}
{ {
if constexpr (MeshType::Dimension == 2) { if constexpr (MeshType::Dimension == 2) {
SmallArray<size_t> y_row_index(m_polynomial_degree + 2); SmallArray<size_t> y_row_index(m_polynomial_degree + 1);
size_t i_y = 0; size_t i_y = 0;
y_row_index[i_y++] = 0; y_row_index[i_y++] = 0;
for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) { for (ssize_t n = m_polynomial_degree + 1; n > 1; --n, ++i_y) {
y_row_index[i_y] = y_row_index[i_y - 1] + n; y_row_index[i_y] = y_row_index[i_y - 1] + n;
++i_y;
} }
m_y_row_index = y_row_index; m_y_row_index = y_row_index;
......
src/scheme/test_reconstruction.cpp
+ 13
− 0
View file @ 4926b9a0
...@@ -6,6 +6,19 @@ ...@@ -6,6 +6,19 @@
void void
test_reconstruction(std::shared_ptr<const DiscreteFunctionVariant> df, size_t degree, size_t number) test_reconstruction(std::shared_ptr<const DiscreteFunctionVariant> df, size_t degree, size_t number)
{ {
if (degree == 1) {
PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, degree};
PolynomialReconstruction rec_builder{descriptor};
Timer t;
for (size_t i = 0; i < number; ++i) {
auto rec = rec_builder.build(df);
}
std::cout << "*** Elapsed time for " << number << " cell center reconstructions: " << t.seconds() << "s ***\n";
}
{ {
PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree}; PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
PolynomialReconstruction rec_builder{descriptor}; PolynomialReconstruction rec_builder{descriptor};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment