diff --git a/src/scheme/DiscreteFunctionDPk.hpp b/src/scheme/DiscreteFunctionDPk.hpp
index 64be52ef27ad72ae9412c7eecc2ecb610e4624ca..2690646b7f88aee4e4e05be88aaeae7fb1b90275 100644
--- a/src/scheme/DiscreteFunctionDPk.hpp
+++ b/src/scheme/DiscreteFunctionDPk.hpp
@@ -135,7 +135,7 @@ class PolynomialCenteredCanonicalBasisView
if constexpr (Dimension == 1) {
const double x_xj = (x - m_xj)[0];
- DataType result = m_coefficients[m_degree];
+ std::remove_const_t<DataType> result = m_coefficients[m_degree];
for (ssize_t i_coeffiencient = m_degree - 1; i_coeffiencient >= 0; --i_coeffiencient) {
result = x_xj * result + m_coefficients[i_coeffiencient];
}
@@ -148,9 +148,9 @@ class PolynomialCenteredCanonicalBasisView
size_t i = m_coefficients.size() - 1;
- DataType result = m_coefficients[i--];
+ std::remove_const_t<DataType> result = m_coefficients[i--];
for (ssize_t i_y = m_degree - 1; i_y >= 0; --i_y) {
- DataType x_result = m_coefficients[i--];
+ std::remove_const_t<DataType> x_result = m_coefficients[i--];
for (ssize_t i_x = m_degree - i_y - 1; i_x >= 0; --i_x) {
x_result = x_xj * x_result + m_coefficients[i--];
}
@@ -166,11 +166,11 @@ class PolynomialCenteredCanonicalBasisView
size_t i = m_coefficients.size() - 1;
- DataType result = m_coefficients[i--];
+ std::remove_const_t<DataType> result = m_coefficients[i--];
for (ssize_t i_z = m_degree - 1; i_z >= 0; --i_z) {
- DataType y_result = m_coefficients[i--];
+ std::remove_const_t<DataType> y_result = m_coefficients[i--];
for (ssize_t i_y = m_degree - i_z - 1; i_y >= 0; --i_y) {
- DataType x_result = m_coefficients[i--];
+ std::remove_const_t<DataType> x_result = m_coefficients[i--];
for (ssize_t i_x = m_degree - i_z - i_y - 1; i_x >= 0; --i_x) {
x_result = x_xj * x_result + m_coefficients[i--];
}
@@ -297,7 +297,7 @@ class DiscreteFunctionDPk
m_xj{MeshDataManager::instance().getMeshData(*m_mesh_v->get<Mesh<Dimension>>()).xj()}
{}
- DiscreteFunctionDPk(const std::shared_ptr<const MeshVariant>& mesh_v, const CellValue<DataType>& cell_array)
+ DiscreteFunctionDPk(const std::shared_ptr<const MeshVariant>& mesh_v, const CellArray<DataType>& cell_array)
: m_mesh_v{mesh_v},
m_degree{BasisView::degreeFromDimension(cell_array.sizeOfArrays())},
m_by_cell_coefficients{cell_array},
@@ -312,7 +312,7 @@ class DiscreteFunctionDPk
{}
template <MeshConcept MeshType>
- DiscreteFunctionDPk(const std::shared_ptr<const MeshType>& p_mesh, const CellValue<DataType>& cell_array)
+ DiscreteFunctionDPk(const std::shared_ptr<const MeshType>& p_mesh, const CellArray<DataType>& cell_array)
: DiscreteFunctionDPk{p_mesh->meshVariant(), cell_array}
{
Assert(m_mesh_v->connectivity().id() == cell_array.connectivity_ptr()->id());
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 74fe08b7d412154d49eaae631387a277e4945ed2..ed6a9cfc69a975180834bc006ae33152ab87d768 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -305,141 +305,17 @@ PolynomialReconstruction::_build(
}
});
- return {};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
- // if (cell_is_owned[cell_j_id]) {
- // const int32_t t = tokens.acquire();
-
- // auto stencil_cell_list = stencil[cell_j_id];
- // ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
-
- // size_t column_begin = 0;
- // for (size_t i_discrete_function_variant = 0;
- // i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
- // // const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
-
- // // std::visit(
- // // [&](auto&& discrete_function) {
- // // using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- // // if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- // // using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
-
- // // const DataType& pj = discrete_function[cell_j_id];
- // // for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- // // const CellId cell_i_id = stencil_cell_list[i];
- // // const DataType& pi_pj = discrete_function[cell_i_id] - pj;
- // // if constexpr (std::is_arithmetic_v<DataType>) {
- // // B(i, column_begin) = pi_pj;
- // // } else if constexpr (is_tiny_vector_v<DataType>) {
- // // for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
- // // B(i, kB) = pi_pj[k];
- // // }
- // // } else if constexpr (is_tiny_matrix_v<DataType>) {
- // // for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- // // for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- // // B(i, column_begin + k * DataType::NumberOfColumns + l) = pi_pj(k, l);
- // // }
- // // }
- // // }
- // // }
-
- // // if constexpr (std::is_arithmetic_v<DataType>) {
- // // ++column_begin;
- // // } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
- // // column_begin += DataType::Dimension;
- // // }
- // // }
- // // },
- // // discrete_function_variant->discreteFunction());
-
- // using DataType = double;
-
- // const DataType& pj = my_discrete_function[cell_j_id];
- // for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- // const CellId cell_i_id = stencil_cell_list[i];
- // const DataType& pi_pj = my_discrete_function[cell_i_id] - pj;
- // B(i, column_begin) = pi_pj;
- // }
-
- // ++column_begin;
- // }
- // ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
-
- // const Rd& Xj = xj[cell_j_id];
-
- // for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- // const CellId cell_i_id = stencil_cell_list[i];
- // const Rd Xi_Xj = xj[cell_i_id] - Xj;
- // for (size_t l = 0; l < MeshType::Dimension; ++l) {
- // A(i, l) = Xi_Xj[l];
- // }
- // }
-
- // const SmallMatrix<double>& X = X_pool[t];
- // Givens::solveCollectionInPlace(A, X, B);
-
- // column_begin = 0;
- // for (size_t i_discrete_function_variant = 0;
- // i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
- // auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
-
- // std::visit(
- // [&](auto&& discrete_function) {
- // using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- // if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- // using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
-
- // const DataType& pj = discrete_function[cell_j_id];
- // auto discrete_function_dpk = mutable_discrete_function_dpk_variant_list[i_discrete_function_variant]
- // .get<DiscreteFunctionDPk<MeshType::Dimension, DataType>>();
- // auto dpk_j = discrete_function_dpk.coefficients(cell_j_id);
- // dpk_j[0] = pj;
-
- // if constexpr (std::is_arithmetic_v<DataType>) {
- // for (size_t i = 0; i < MeshType::Dimension; ++i) {
- // auto& dpk_j_ip1 = dpk_j[i + 1];
- // dpk_j_ip1 = X(i, column_begin);
- // }
- // ++column_begin;
- // } else if constexpr (is_tiny_vector_v<DataType>) {
- // for (size_t i = 0; i < MeshType::Dimension; ++i) {
- // auto& dpk_j_ip1 = dpk_j[i + 1];
- // for (size_t k = 0; k < DataType::Dimension; ++k) {
- // dpk_j_ip1[k] = X(i, column_begin + k);
- // }
- // }
- // column_begin += DataType::Dimension;
- // } else if constexpr (is_tiny_matrix_v<DataType>) {
- // for (size_t i = 0; i < MeshType::Dimension; ++i) {
- // auto& dpk_j_ip1 = dpk_j[i + 1];
- // for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- // for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- // dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
- // }
- // }
- // }
- // column_begin += DataType::Dimension;
- // }
-
- // } else {
- // throw NotImplementedError("discrete function type");
- // }
- // },
- // discrete_function_variant->discreteFunction());
- // }
-
- // tokens.release(t);
- // }
- // });
-
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
- // for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
- // discrete_function_dpk_variant_list.push_back(
- // std::make_shared<const DiscreteFunctionDPkVariant>(discrete_function_dpk_variant_p));
- // }
+ for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
+ std::visit(
+ [&](auto&& mutable_function_dpk) {
+ synchronize(mutable_function_dpk.cellArrays());
+ discrete_function_dpk_variant_list.push_back(
+ std::make_shared<DiscreteFunctionDPkVariant>(mutable_function_dpk));
+ },
+ discrete_function_dpk_variant_p.mutableDiscreteFunctionDPk());
+ }
return discrete_function_dpk_variant_list;
}
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 72036d09fd9995cbf8ef8a1635dfbb7fa0c17983..d90d0e84a153d67834942756ddb18cc746c5b05b 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -12,10 +12,21 @@
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <MeshDataBaseForTests.hpp>
+template <typename... DataType>
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
+build_list(DiscreteFunctionP0<DataType>... input)
+{
+ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector;
+ variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(input...));
+
+ return variant_vector;
+}
+
// clazy:excludeall=non-pod-global-static
TEST_CASE("PolynomialReconstruction", "[scheme]")
@@ -40,8 +51,9 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ auto reconstructions = PolynomialReconstruction{}.build(build_list(fh));
+
+ auto dpk = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
{
double max_mean_error = 0;