diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index fc9b88d0e76aba7b030712be95a0ad14e0855943..74fe08b7d412154d49eaae631387a277e4945ed2 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -128,8 +128,7 @@ PolynomialReconstruction::_build(
const size_t degree = 1;
const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto cell_is_owned = mesh.connectivity().cellIsOwned();
const size_t max_stencil_size = [&]() {
@@ -182,20 +181,20 @@ PolynomialReconstruction::_build(
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) {
- auto discrete_function_variant = discrete_function_variant_list[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];
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+ size_t column_begin = 0;
+ 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;
@@ -224,244 +223,12 @@ PolynomialReconstruction::_build(
discrete_function_variant->discreteFunction());
}
- 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];
- }
- }
-
- 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) {
- // const_discrete_function_dpk_variant_list.push_back(discrete_function_dpk_variant_p);
- // }
-
- return discrete_function_dpk_variant_list;
-}
-
-std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
-PolynomialReconstruction::build(
- const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
-{
- if (not hasSameMesh(discrete_function_variant_list)) {
- throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
- }
-
- auto mesh_v = getCommonMesh(discrete_function_variant_list);
-
- return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
- mesh_v->variant());
-}
-
-/**************************************/
-
-template <MeshConcept MeshType>
-[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
-PolynomialReconstruction::_build2(
- const std::shared_ptr<const MeshType>& p_mesh,
- const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
-{
- const MeshType& mesh = *p_mesh;
-
- using Rd = TinyVector<MeshType::Dimension>;
-
- const size_t number_of_columns = [&]() {
- size_t n = 0;
- for (auto discrete_function_variant_p : discrete_function_variant_list) {
- n += std::visit(
- [](auto&& discrete_function) -> size_t {
- using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
- if constexpr (std::is_same_v<data_type, double>) {
- return 1;
- } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
- return data_type::Dimension;
- } else {
- throw UnexpectedError("unexpected data type " + demangle<data_type>());
- }
- } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- return discrete_function.size();
- } else {
- throw UnexpectedError("unexpected discrete function type");
- }
- },
- discrete_function_variant_p->discreteFunction());
- }
- return n;
- }();
-
- const size_t degree = 1;
- const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
-
- Array<const double*> value_begining(discrete_function_variant_list.size());
- Array<size_t> value_size(discrete_function_variant_list.size());
-
- for (size_t i = 0; i < discrete_function_variant_list.size(); ++i) {
- DiscreteFunctionVariant discrete_function_v = *(discrete_function_variant_list[i]);
- std::visit(
- [=](auto&& discrete_function) {
- using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
- if constexpr (std::is_same_v<data_type, double>) {
- value_begining[i] = &(discrete_function[CellId{0}]);
- value_size[i] = 1;
- } else if constexpr (is_tiny_vector_v<data_type>) {
- value_begining[i] = &(discrete_function[CellId{0}][0]);
- value_size[i] = data_type::Dimension;
- } else if constexpr (is_tiny_matrix_v<data_type>) {
- value_begining[i] = &(discrete_function[CellId{0}](0, 0));
- value_size[i] = data_type::Dimension;
- } else {
- throw UnexpectedError("unexpected data type " + demangle<data_type>());
- }
- } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- throw NotImplementedError("Discrete Function P0 Vector");
- } else {
- throw UnexpectedError("unexpected discrete function type");
- }
- },
- discrete_function_v.discreteFunction());
- }
-
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- auto cell_is_owned = mesh.connectivity().cellIsOwned();
-
- const size_t max_stencil_size = [&]() {
- size_t max_size = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto stencil_cell_list = stencil[cell_id];
- if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
- max_size = stencil_cell_list.size();
- }
- }
- return max_size;
- }();
-
- Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
- Kokkos::Experimental::UniqueTokenScope::Global>
- tokens;
-
- std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
- 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::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
- mutable_discrete_function_dpk_variant_list.push_back(
- DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, degree));
- } else {
- throw NotImplementedError("discrete function type");
- }
- },
- discrete_function_variant->discreteFunction());
- }
-
- SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
-
- for (size_t i = 0; i < A_pool.size(); ++i) {
- A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
- B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
- X_pool[i] = SmallMatrix<double>(MeshType::Dimension, number_of_columns);
- }
-
- 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 i_column = 0;
-
- for (size_t i_value = 0; i_value < value_begining.size(); ++i_value) {
- const size_t i_value_size = value_size[i_value];
-
- const size_t cell_j_value_index = cell_j_id * i_value_size;
-
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
- const size_t cell_i_value_index = cell_i_id * i_value_size;
- size_t j_index = cell_j_value_index;
- size_t i_index = cell_i_value_index;
-
- for (size_t shift = 0; shift < i_value_size; ++shift, ++i_index, ++j_index) {
- B(i, i_column++) = value_begining[i_value][i_index] - value_begining[i_value][j_index];
- }
- }
- }
+ // for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ // const CellId cell_i_id = stencil_cell_list[i];
+ // for (size_t j = 0; j < number_of_columns; ++j) {
+ // B(i, j) = values[cell_i_id][j] - values[cell_j_id][j];
+ // }
+ // }
ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
@@ -475,76 +242,210 @@ PolynomialReconstruction::_build2(
}
}
- // SmallMatrix<double> X = X_pool[t];
-
- TinyMatrix<MeshType::Dimension, 1> X;
+ const SmallMatrix<double>& X = X_pool[t];
Givens::solveCollectionInPlace(A, X, B);
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) {
- auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+ for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
+ ++i_dpk_variant) {
+ const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
+
+ const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_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>;
+ [&](auto&& dpk_function, auto&& p0_function) {
+ using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
+ using P0FunctionT = std::decay_t<decltype(p0_function)>;
+ using DataType = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
+ using P0DataType = std::remove_const_t<std::decay_t<typename P0FunctionT::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_same_v<DataType, P0DataType>) {
+ if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
+ auto dpk_j = dpk_function.coefficients(cell_j_id);
+ dpk_j[0] = p0_function[cell_j_id];
- 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);
+ 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 += 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;
+ } 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("unexpected data type");
}
- column_begin += DataType::Dimension;
+ } else {
+ throw NotImplementedError("non scalar P0 functions");
}
-
} else {
- throw NotImplementedError("discrete function type");
+ throw UnexpectedError("incompatible data types");
}
},
- discrete_function_variant->discreteFunction());
+ dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
}
tokens.release(t);
}
});
+ 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) {
- // const_discrete_function_dpk_variant_list.push_back(discrete_function_dpk_variant_p);
+ // discrete_function_dpk_variant_list.push_back(
+ // std::make_shared<const DiscreteFunctionDPkVariant>(discrete_function_dpk_variant_p));
// }
return discrete_function_dpk_variant_list;
}
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
-PolynomialReconstruction::build2(
+PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
if (not hasSameMesh(discrete_function_variant_list)) {
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 0612e605b55f9a580e8e3bd13e6f46eba18c7dde..b55e84fee171e21b762e4934054f6cb2b04b82e7 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -120,11 +120,6 @@ class PolynomialReconstruction
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>> _build2(
- const std::shared_ptr<const MeshType>& p_mesh,
- const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
-
public:
template <MeshConcept MeshType, typename DataType>
[[nodiscard]] PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
@@ -193,7 +188,7 @@ class PolynomialReconstruction
}
}
- TinyVector<MeshType::Dimension> x;
+ SmallVector<double> x(MeshType::Dimension);
Givens::solveInPlace(A, x, b);
auto dpk_j = dpk.coefficients(cell_j_id);
@@ -328,9 +323,6 @@ class PolynomialReconstruction
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
- [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build2(
- const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
-
PolynomialReconstruction() {}
};
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 9ec19396e28f00119412582c9e42cd2d4ed25bcb..08527f96d0dd35838b3bbeebf5cb87db96c42b94 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -31,30 +31,63 @@ test_reconstruction(const std::shared_ptr<const MeshVariant>& mesh_v,
std::cout << "finished!\n";
}
+void
+test_reconstruction_one(const std::shared_ptr<const MeshVariant>& mesh_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_v)
+{
+ std::visit(
+ [&](auto&& p_mesh) {
+ std::visit(
+ [&](auto&& discrete_function) {
+ using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+ PolynomialReconstruction reconstruction;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+ for (size_t i = 0; i < 30; ++i) {
+ auto res = reconstruction.build(*p_mesh, discrete_function);
+ }
+ }
+ },
+ discrete_function_v->discreteFunction());
+ },
+ mesh_v->variant());
+}
+
void
test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
{
- std::cout << "** variable list at once (30 times)\n";
+ std::cout << "** variable list one by one (30 times)\n";
{
- PolynomialReconstruction reconstruction;
Timer t;
- for (size_t i = 0; i < 30; ++i) {
- auto res = reconstruction.build(discrete_function_variant_list);
+ for (auto discrete_function_v : discrete_function_variant_list) {
+ std::visit(
+ [&](auto&& discrete_function) {
+ auto mesh_v = discrete_function.meshVariant();
+ std::visit(
+ [&](auto&& p_mesh) {
+ using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+ PolynomialReconstruction reconstruction;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+ for (size_t i = 0; i < 30; ++i) {
+ auto res = reconstruction.build(*p_mesh, discrete_function);
+ }
+ }
+ },
+ mesh_v->variant());
+ },
+ discrete_function_v->discreteFunction());
}
t.pause();
std::cout << "t = " << t << '\n';
}
- std::cout << "finished!\n";
-
std::cout << "** variable list at once (30 times)\n";
{
PolynomialReconstruction reconstruction;
Timer t;
for (size_t i = 0; i < 30; ++i) {
- auto res = reconstruction.build2(discrete_function_variant_list);
+ auto res = reconstruction.build(discrete_function_variant_list);
}
t.pause();
std::cout << "t = " << t << '\n';