From ae1cbc1060ac59874c06a058ccdc30ae28575a08 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 16 Jul 2024 19:23:20 +0200
Subject: [PATCH] Begin code clean-up
---
src/language/modules/SchemeModule.cpp | 11 -
src/scheme/PolynomialReconstruction.cpp | 454 +++++++++++++-----------
src/scheme/PolynomialReconstruction.hpp | 238 +------------
src/scheme/test_reconstruction.cpp | 51 ---
src/scheme/test_reconstruction.hpp | 3 -
5 files changed, 261 insertions(+), 496 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index f96e96f13..0f3dadaad 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -671,17 +671,6 @@ SchemeModule::SchemeModule()
));
-#warning REMOVE AFTER TESTS FINISHED
- this->_addBuiltinFunction("test_reconstruction",
- std::function(
-
- [](const std::shared_ptr<const MeshVariant> mesh_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_v) -> void {
- test_reconstruction(mesh_v, discrete_function_v);
- }
-
- ));
-
#warning REMOVE AFTER TESTS FINISHED
this->_addBuiltinFunction("test_reconstruction",
std::function(
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 69d24b9e5..f9e55f19c 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -110,6 +110,225 @@ class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
~MutableDiscreteFunctionDPkVariant() = default;
};
+template <typename ConformTransformationT>
+void
+_computeEjkMean(const QuadratureFormula<1>& quadrature,
+ const ConformTransformationT& T,
+ const TinyVector<1>& Xj,
+ const double Vi,
+ const size_t degree,
+ const size_t dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ using Rd = TinyVector<1>;
+ mean_of_ejk.fill(0);
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const Rd& xi_q = quadrature.point(i_q);
+ const double detT = T.jacobianDeterminant();
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+
+ {
+ size_t k = 0;
+ inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ for (; k <= degree; ++k) {
+ inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
+ }
+ }
+
+ for (size_t k = 1; k < dimension; ++k) {
+ mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
+ }
+ }
+}
+
+template <typename ConformTransformationT>
+void
+_computeEjkMean(const QuadratureFormula<2>& quadrature,
+ const ConformTransformationT& T,
+ const TinyVector<2>& Xj,
+ const double Vi,
+ const size_t degree,
+ const size_t dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ using Rd = TinyVector<2>;
+
+ mean_of_ejk.fill(0);
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const Rd& xi_q = quadrature.point(i_q);
+ const double detT = [&] {
+ if constexpr (std::is_same_v<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
+ return T.jacobianDeterminant();
+ } else {
+ return T.jacobianDeterminant(xi_q);
+ }
+ }();
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+
+ {
+ size_t k = 0;
+ inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ for (; k <= degree; ++k) {
+ inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
+ }
+
+ for (size_t i_y = 1; i_y <= degree; ++i_y) {
+ const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
+ for (size_t l = 0; l <= degree - i_y; ++l) {
+ inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ for (size_t k = 1; k < dimension; ++k) {
+ mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
+ }
+ }
+}
+
+template <typename ConformTransformationT>
+void
+_computeEjkMean(const QuadratureFormula<3>& quadrature,
+ const ConformTransformationT& T,
+ const TinyVector<3>& Xj,
+ const double Vi,
+ const size_t degree,
+ const size_t dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ using Rd = TinyVector<3>;
+ mean_of_ejk.fill(0);
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const Rd& xi_q = quadrature.point(i_q);
+ const double detT = [&] {
+ if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
+ return T.jacobianDeterminant();
+ } else {
+ return T.jacobianDeterminant(xi_q);
+ }
+ }();
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+ const double z_zj = X_Xj[2];
+
+ {
+ size_t k = 0;
+ inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ for (; k <= degree; ++k) {
+ inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
+ }
+
+ for (size_t i_y = 1; i_y <= degree; ++i_y) {
+ const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
+ for (size_t l = 0; l <= degree - i_y; ++l) {
+ inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= degree; ++i_z) {
+ const size_t begin_i_z_1 = i_z * (i_z + 1) * (i_z + 2) / 6 - 1;
+ for (size_t i_y = 0; i_y <= degree - i_z; ++i_y) {
+ const size_t begin_i_y_1 = (i_y * (2 * degree - i_y + 3)) / 2 + begin_i_z_1;
+ for (size_t l = 0; l <= degree - i_y - i_z; ++l) {
+ inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < dimension; ++k) {
+ mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
+ }
+ }
+}
+
+size_t
+PolynomialReconstruction::_getNumberOfColumns(
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+ size_t number_of_columns = 0;
+ for (auto discrete_function_variant_p : discrete_function_variant_list) {
+ number_of_columns += 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 {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected data type " + demangle<data_type>());
+ // LCOV_EXCL_STOP
+ }
+ } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+ return discrete_function.size();
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
+ }
+ },
+ discrete_function_variant_p->discreteFunction());
+ }
+ return number_of_columns;
+}
+
+template <MeshConcept MeshType>
+std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>
+PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
+ const std::shared_ptr<const MeshType>& p_mesh,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+ 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, m_descriptor.degree()));
+ } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+ using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
+ mutable_discrete_function_dpk_variant_list.push_back(
+ DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
+ discrete_function.size()));
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
+ }
+ },
+ discrete_function_variant->discreteFunction());
+ }
+
+ return mutable_discrete_function_dpk_variant_list;
+}
+
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
@@ -120,35 +339,7 @@ PolynomialReconstruction::_build(
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 {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected data type " + demangle<data_type>());
- // LCOV_EXCL_STOP
- }
- } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- return discrete_function.size();
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected discrete function type");
- // LCOV_EXCL_STOP
- }
- },
- discrete_function_variant_p->discreteFunction());
- }
- return n;
- }();
+ const size_t number_of_columns = this->_getNumberOfColumns(discrete_function_variant_list);
const size_t basis_dimension =
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
@@ -178,31 +369,8 @@ PolynomialReconstruction::_build(
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, m_descriptor.degree()));
- } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
- mutable_discrete_function_dpk_variant_list.push_back(
- DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
- discrete_function.size()));
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected discrete function type");
- // LCOV_EXCL_STOP
- }
- },
- discrete_function_variant->discreteFunction());
- }
+ auto mutable_discrete_function_dpk_variant_list =
+ this->_createMutableDiscreteFunctionDPKVariantList(p_mesh, discrete_function_variant_list);
SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
@@ -310,39 +478,11 @@ PolynomialReconstruction::_build(
}
}
} else {
- SmallArray<double> inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
- SmallArray<double> mean_j_of_ejk = mean_j_of_ejk_pool[t];
- SmallArray<double> mean_i_of_ejk = mean_i_of_ejk_pool[t];
+ SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
+ SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
if constexpr (MeshType::Dimension == 1) {
- auto compute_mean_ejk = [&inv_Vj_wq_detJ_ek](const auto& quadrature, const auto& T, const Rd& Xj,
- const double Vi, const size_t degree, const size_t dimension,
- SmallArray<double>& mean_of_ejk) {
- mean_of_ejk.fill(0);
-
- for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
- const double wq = quadrature.weight(i_q);
- const Rd& xi_q = quadrature.point(i_q);
- const double detT = T.jacobianDeterminant();
-
- const Rd X_Xj = T(xi_q) - Xj;
-
- const double x_xj = X_Xj[0];
-
- {
- size_t k = 0;
- inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
- for (; k <= degree; ++k) {
- inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
- }
- }
-
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
- }
- }
- };
-
const Rd& Xj = xj[cell_j_id];
if (cell_type[cell_j_id] == CellType::Line) {
@@ -353,7 +493,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
@@ -370,8 +511,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
@@ -382,48 +523,6 @@ PolynomialReconstruction::_build(
}
}
} else if constexpr (MeshType::Dimension == 2) {
- auto compute_mean_ejk = [&inv_Vj_wq_detJ_ek](const auto& quadrature, const auto& T, const Rd& Xj,
- const double Vi, const size_t degree, const size_t dimension,
- SmallArray<double>& mean_of_ejk) {
- mean_of_ejk.fill(0);
-
- for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
- const double wq = quadrature.weight(i_q);
- const Rd& xi_q = quadrature.point(i_q);
- const double detT = [&] {
- if constexpr (std::is_same_v<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
- return T.jacobianDeterminant();
- } else {
- return T.jacobianDeterminant(xi_q);
- }
- }();
-
- const Rd X_Xj = T(xi_q) - Xj;
-
- const double x_xj = X_Xj[0];
- const double y_yj = X_Xj[1];
-
- {
- size_t k = 0;
- inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
- for (; k <= degree; ++k) {
- inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
- }
-
- for (size_t i_y = 1; i_y <= degree; ++i_y) {
- const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
- for (size_t l = 0; l <= degree - i_y; ++l) {
- inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
- }
- }
- }
-
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
- }
- }
- };
-
const Rd& Xj = xj[cell_j_id];
if (cell_type[cell_j_id] == CellType::Triangle) {
@@ -434,7 +533,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else if (cell_type[cell_j_id] == CellType::Quadrangle) {
auto cell_node = cell_to_node_matrix[cell_j_id];
@@ -444,7 +544,8 @@ PolynomialReconstruction::_build(
const auto& quadrature = QuadratureManager::instance().getSquareFormula(
GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
@@ -461,8 +562,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else if (cell_type[cell_i_id] == CellType::Quadrangle) {
const SquareTransformation<2> T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]]};
@@ -470,8 +571,8 @@ PolynomialReconstruction::_build(
const auto& quadrature = QuadratureManager::instance().getSquareFormula(
GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
@@ -483,59 +584,6 @@ PolynomialReconstruction::_build(
}
} else if constexpr (MeshType::Dimension == 3) {
- auto compute_mean_ejk = [&inv_Vj_wq_detJ_ek](const auto& quadrature, const auto& T, const Rd& Xj,
- const double Vi, const size_t degree, const size_t dimension,
- SmallArray<double>& mean_of_ejk) {
- mean_of_ejk.fill(0);
-
- for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
- const double wq = quadrature.weight(i_q);
- const Rd& xi_q = quadrature.point(i_q);
- const double detT = [&] {
- if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
- return T.jacobianDeterminant();
- } else {
- return T.jacobianDeterminant(xi_q);
- }
- }();
-
- const Rd X_Xj = T(xi_q) - Xj;
-
- const double x_xj = X_Xj[0];
- const double y_yj = X_Xj[1];
- const double z_zj = X_Xj[2];
-
- {
- size_t k = 0;
- inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
- for (; k <= degree; ++k) {
- inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
- }
-
- for (size_t i_y = 1; i_y <= degree; ++i_y) {
- const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
- for (size_t l = 0; l <= degree - i_y; ++l) {
- inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
- }
- }
-
- for (size_t i_z = 1; i_z <= degree; ++i_z) {
- const size_t begin_i_z_1 = i_z * (i_z + 1) * (i_z + 2) / 6 - 1;
- for (size_t i_y = 0; i_y <= degree - i_z; ++i_y) {
- const size_t begin_i_y_1 = (i_y * (2 * degree - i_y + 3)) / 2 + begin_i_z_1;
- for (size_t l = 0; l <= degree - i_y - i_z; ++l) {
- inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
- }
- }
- }
- }
-
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
- }
- }
- };
-
const Rd& Xj = xj[cell_j_id];
if (cell_type[cell_j_id] == CellType::Tetrahedron) {
@@ -546,7 +594,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else if (cell_type[cell_j_id] == CellType::Prism) {
auto cell_node = cell_to_node_matrix[cell_j_id];
@@ -557,7 +606,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else if (cell_type[cell_j_id] == CellType::Pyramid) {
auto cell_node = cell_to_node_matrix[cell_j_id];
@@ -568,7 +618,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else if (cell_type[cell_j_id] == CellType::Hexahedron) {
auto cell_node = cell_to_node_matrix[cell_j_id];
@@ -579,7 +630,8 @@ PolynomialReconstruction::_build(
const auto& quadrature = QuadratureManager::instance().getCubeFormula(
GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_j_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
@@ -597,8 +649,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else if (cell_type[cell_i_id] == CellType::Prism) {
const PrismTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], //
@@ -607,8 +659,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else if (cell_type[cell_i_id] == CellType::Pyramid) {
const PyramidTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
@@ -617,8 +669,8 @@ PolynomialReconstruction::_build(
const auto& quadrature =
QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else if (cell_type[cell_i_id] == CellType::Hexahedron) {
const CubeTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
@@ -627,8 +679,8 @@ PolynomialReconstruction::_build(
const auto& quadrature = QuadratureManager::instance().getCubeFormula(
GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
- compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- mean_i_of_ejk);
+ _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ inv_Vj_wq_detJ_ek, mean_i_of_ejk);
} else {
throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index c0da0f260..f52314872 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -29,242 +29,20 @@ class PolynomialReconstruction
const PolynomialReconstructionDescriptor m_descriptor;
+ size_t _getNumberOfColumns(
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
+ template <MeshConcept MeshType>
+ std::vector<MutableDiscreteFunctionDPkVariant> _createMutableDiscreteFunctionDPKVariantList(
+ 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,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
public:
- template <MeshConcept MeshType, typename DataType>
- PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
- build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function) const
- {
- using Rd = TinyVector<MeshType::Dimension>;
-
- const size_t degree = 1;
- const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
-
- 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;
-
- DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
-
- if constexpr (is_polygonal_mesh_v<MeshType>) {
- if constexpr (std::is_arithmetic_v<DataType>) {
- SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallVector<double>> b_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] = SmallVector<double>(max_stencil_size);
- }
-
- 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];
-
- ShrinkVectorView b(b_pool[t], stencil_cell_list.size());
-
- const double pj = p0_function[cell_j_id];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
-
- b[i] = p0_function[cell_i_id] - pj;
- }
-
- 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];
- }
- }
-
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
- const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
- for (size_t l = 0; l < MeshType::Dimension; ++l) {
- A(i, l) *= weight;
- }
- b[i] *= weight;
- }
-
- SmallVector<double> x(MeshType::Dimension);
- Givens::solveInPlace(A, x, b);
-
- auto dpk_j = dpk.coefficients(cell_j_id);
-
- dpk_j[0] = pj;
-
- for (size_t l = 0; l < MeshType::Dimension; ++l) {
- dpk_j[1 + l] = x[l];
- }
-
- tokens.release(t);
- }
- });
- } else if constexpr (is_tiny_vector_v<DataType>) {
- SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallMatrix<double>> B_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, DataType::Dimension);
- }
-
- 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());
-
- const DataType& pj = p0_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 = p0_function[cell_i_id] - pj;
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- B(i, k) = pi_pj[k];
- }
- }
-
- 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];
- }
- }
-
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
- const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
- for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
- A(i, l) *= weight;
- }
- for (size_t l = 0; l < DataType::Dimension; ++l) {
- B(i, l) *= weight;
- }
- }
-
- TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
- Givens::solveCollectionInPlace(A, X, B);
-
- auto dpk_j = dpk.coefficients(cell_j_id);
-
- dpk_j[0] = pj;
- 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, k);
- }
- }
-
- tokens.release(t);
- }
- });
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallMatrix<double>> B_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, DataType::Dimension);
- }
-
- 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());
-
- const DataType& pj = p0_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 = p0_function[cell_i_id] - pj;
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(i, k * DataType::NumberOfColumns + l) = pi_pj(k, l);
- }
- }
- }
-
- 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];
- }
- }
-
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
- const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
- for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
- A(i, l) *= weight;
- }
- for (size_t l = 0; l < DataType::Dimension; ++l) {
- B(i, l) *= weight;
- }
- }
-
- TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
- Givens::solveCollectionInPlace(A, X, B);
-
- auto dpk_j = dpk.coefficients(cell_j_id);
- dpk_j[0] = pj;
-
- 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, k * DataType::NumberOfColumns + l);
- }
- }
- }
-
- tokens.release(t);
- }
- });
- } else {
- throw NotImplementedError("dealing with " + demangle<DataType>());
- }
- }
-
- synchronize(dpk.cellArrays());
- return dpk;
- }
-
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 42ec98cc9..3d255a547 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -3,57 +3,6 @@
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/Timer.hpp>
-void
-test_reconstruction(const std::shared_ptr<const MeshVariant>& mesh_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_v)
-{
- std::cout << "** one variable at a time (30 times)\n";
- PolynomialReconstructionDescriptor descriptor{1};
-
- std::visit(
- [&](auto&& p_mesh) {
- std::visit(
- [&](auto&& discrete_function) {
- using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- PolynomialReconstruction reconstruction{descriptor};
- if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- Timer t;
- for (size_t i = 0; i < 30; ++i) {
- auto res = reconstruction.build(*p_mesh, discrete_function);
- }
- t.pause();
- std::cout << "t = " << t << '\n';
- }
- },
- discrete_function_v->discreteFunction());
- },
- mesh_v->variant());
-
- 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)
-{
- PolynomialReconstructionDescriptor descriptor{1};
- std::visit(
- [&](auto&& p_mesh) {
- std::visit(
- [&](auto&& discrete_function) {
- using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- PolynomialReconstruction reconstruction{descriptor};
- 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)
{
diff --git a/src/scheme/test_reconstruction.hpp b/src/scheme/test_reconstruction.hpp
index f74bf0781..6620acd44 100644
--- a/src/scheme/test_reconstruction.hpp
+++ b/src/scheme/test_reconstruction.hpp
@@ -7,9 +7,6 @@
#include <vector>
-void test_reconstruction(const std::shared_ptr<const MeshVariant>& mesh_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_variant_list);
-
void test_reconstruction(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list);
--
GitLab