diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 464ce1a0f532b263e74a42d795dcb25d9ef84f12..3123063361c458385f34a2ab02b983bbfcad8517 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -370,12 +370,13 @@ class PolynomialReconstruction::Internal
}
}
+ template <typename ReconstructionMatrixBuilderType>
static void
populateDiscreteFunctionDPkByCell(
const CellId& cell_j_id,
const size_t& degree,
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<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>&
mutable_discrete_function_dpk_variant_list)
@@ -400,9 +401,12 @@ class PolynomialReconstruction::Internal
dpk_j[0] = p0_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
}
}
@@ -412,11 +416,14 @@ class PolynomialReconstruction::Internal
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- auto& dpk_j_0 = dpk_j[0];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ }
}
}
}
@@ -429,12 +436,15 @@ class PolynomialReconstruction::Internal
}
column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- auto& dpk_j_0 = dpk_j[0];
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- dpk_j_0(k, l) -= X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ dpk_j_0(k, l) -= X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
+ }
}
}
}
@@ -469,9 +479,12 @@ class PolynomialReconstruction::Internal
const size_t component_begin = l * size;
dpk_j[component_begin] = cell_vector[l];
if constexpr (std::is_arithmetic_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ if (degree > 1) {
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
}
}
@@ -481,11 +494,14 @@ class PolynomialReconstruction::Internal
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- auto& dpk_j_0 = dpk_j[component_begin];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[component_begin];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ }
}
}
}
@@ -498,12 +514,16 @@ class PolynomialReconstruction::Internal
}
column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) {
- if (degree > 1) {
- for (size_t i = 0; i < X.numberOfRows(); ++i) {
- auto& dpk_j_0 = dpk_j[component_begin];
- for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
- 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];
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[component_begin];
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ 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];
+ }
}
}
}
@@ -660,12 +680,14 @@ PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
}
}
-template <MeshConcept MeshType>
+template <typename ReconstructionMatrixBuilderType, 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
{
+ static_assert(std::is_same_v<MeshType, typename ReconstructionMatrixBuilderType::MeshType>);
+
const MeshType& mesh = *p_mesh;
using Rd = TinyVector<MeshType::Dimension>;
@@ -748,45 +770,19 @@ PolynomialReconstruction::_build(
SmallArray<SmallVector<double>> G_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) {
A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
G_pool[i] = SmallVector<double>(basis_dimension - 1);
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>>>
- boundary_integral_reconstruction_matrix_builder_pool{A_pool.size()};
-
- if constexpr (MeshType::Dimension == 2) {
- for (size_t t = 0; t < boundary_integral_reconstruction_matrix_builder_pool.size(); ++t) {
- SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ SmallArray<std::shared_ptr<ReconstructionMatrixBuilderType>> reconstruction_matrix_builder_pool(A_pool.size());
- 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);
- }
+ for (size_t t = 0; t < reconstruction_matrix_builder_pool.size(); ++t) {
+ reconstruction_matrix_builder_pool[t] =
+ std::make_shared<ReconstructionMatrixBuilderType>(*p_mesh, m_descriptor.degree(), symmetry_origin_list,
+ symmetry_normal_list, stencil_array);
}
parallel_for(
@@ -799,25 +795,8 @@ PolynomialReconstruction::_build(
Internal<MeshType>::buildB(cell_j_id, stencil_array, discrete_function_variant_list, symmetry_normal_list, B);
- if ((m_descriptor.degree() == 1) and
- (m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
- 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");
- }
+ ReconstructionMatrixBuilderType& reconstruction_matrix_builder = *reconstruction_matrix_builder_pool[t];
+ reconstruction_matrix_builder.build(cell_j_id, A);
if (m_descriptor.rowWeighting()) {
Internal<MeshType>::rowWeighting(cell_j_id, stencil_array, xj, symmetry_origin_list, symmetry_normal_list, A,
@@ -835,9 +814,10 @@ PolynomialReconstruction::_build(
Givens::solveCollectionInPlace(A, X, B);
}
- Internal<MeshType>::populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X,
- mean_j_of_ejk_pool[t], discrete_function_variant_list,
- mutable_discrete_function_dpk_variant_list);
+ Internal<MeshType>::template populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X,
+ reconstruction_matrix_builder,
+ discrete_function_variant_list,
+ mutable_discrete_function_dpk_variant_list);
tokens.release(t);
}
@@ -858,6 +838,34 @@ PolynomialReconstruction::_build(
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>>
PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index fabd3e8ee30ed5ae5f5176539560ed343d54ef93..4f76844711ada4641b69234c5ac787c472247a37 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -38,6 +38,11 @@ class PolynomialReconstruction
const std::shared_ptr<const MeshType>& p_mesh,
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>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
const std::shared_ptr<const MeshType>& p_mesh,
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index d22ecbb140709829a99c34f45f42f4efc0ed3e56..105c4640e4d534c75c0f956b39d8347846f7ac48 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -14,9 +14,14 @@
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
-template <MeshConcept MeshType>
+template <MeshConcept MeshTypeT>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = true;
+
private:
using Rd = TinyVector<MeshType::Dimension>;
@@ -142,6 +147,13 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
}
public:
+ PUGS_INLINE
+ SmallArray<const double>
+ meanjOfEjk() const
+ {
+ return m_mean_j_of_ejk;
+ }
+
template <typename MatrixType>
void
build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
@@ -188,7 +200,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
const size_t polynomial_degree,
- const SmallArray<double>& mean_j_of_ejk,
const SmallArray<const Rd>& symmetry_origin_list,
const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array)
@@ -197,7 +208,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
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_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_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
index 9f3ee95247f8bfbd870dfbaf62983b7b90a16aa4..ae3488d8125dc4690856294494c1218902427693 100644
--- a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
@@ -8,9 +8,14 @@
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
-template <MeshConcept MeshType>
+template <MeshConcept MeshTypeT>
class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = false;
+
private:
using Rd = TinyVector<MeshType::Dimension>;
@@ -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 CellToCellStencilArray& stencil_array,
- const CellValue<const Rd>& xj)
- : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(1)},
+ 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{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;
};
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 1e5b7d715270137f0c9ba09a433c38155939a7ee..fd877ccdbfa2d7de8eb32150c6b4ad310d17c4bc 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -21,9 +21,14 @@
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
-template <MeshConcept MeshType>
+template <MeshConcept MeshTypeT>
class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = true;
+
private:
using Rd = TinyVector<MeshType::Dimension>;
@@ -460,6 +465,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
}
public:
+ PUGS_INLINE
+ SmallArray<const double>
+ meanjOfEjk() const
+ {
+ return m_mean_j_of_ejk;
+ }
+
template <typename MatrixType>
void
build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
@@ -502,7 +514,6 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
const size_t polynomial_degree,
- const SmallArray<double>& mean_j_of_ejk,
const SmallArray<const Rd>& symmetry_origin_list,
const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array)
@@ -511,7 +522,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
m_polynomial_degree{polynomial_degree},
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_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
@@ -524,14 +535,13 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
m_xr{mesh.xr()}
{
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;
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;
- ++i_y;
}
m_y_row_index = y_row_index;
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index c0bfa81f18e5adc36050f39f41a8024f446d2318..1353ab46149f333ca4ba2cdf8e72f2011932919d 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -6,6 +6,19 @@
void
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};
PolynomialReconstruction rec_builder{descriptor};