diff --git a/tests/DiscreteFunctionDPkForTests.hpp b/tests/DiscreteFunctionDPkForTests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b51861d811b2d8a35c91375bc9e8dd368e2b281
--- /dev/null
+++ b/tests/DiscreteFunctionDPkForTests.hpp
@@ -0,0 +1,365 @@
+#ifndef DISCRETE_FUNCTION_DPK_FOR_TESTS_HPP
+#define DISCRETE_FUNCTION_DPK_FOR_TESTS_HPP
+
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <geometry/PrismTransformation.hpp>
+#include <geometry/PyramidTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <type_traits>
+
+namespace test_only
+{
+
+template <MeshConcept MeshType, typename DataType>
+DiscreteFunctionP0<std::remove_const_t<DataType>>
+exact_projection(const MeshType& mesh,
+ size_t degree,
+ std::function<DataType(const TinyVector<MeshType::Dimension>&)> exact_function)
+{
+ DiscreteFunctionP0<std::remove_const_t<DataType>> P0_function{mesh.meshVariant()};
+
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto xr = mesh.xr();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ auto cell_type = mesh.connectivity().cellType();
+
+ auto sum = [&exact_function, &Vj](const CellId cell_id, const auto& T,
+ const auto& qf) -> std::remove_const_t<DataType> {
+ std::remove_const_t<DataType> integral =
+ (qf.weight(0) * T.jacobianDeterminant(qf.point(0))) * exact_function(T(qf.point(0)));
+ for (size_t i_quadrarture = 1; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ integral += (qf.weight(i_quadrarture) * T.jacobianDeterminant(qf.point(i_quadrarture))) *
+ exact_function(T(qf.point(i_quadrarture)));
+ }
+ return 1. / Vj[cell_id] * integral;
+ };
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ if constexpr (MeshType::Dimension == 1) {
+ LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{degree + 1});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{degree + 2});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{degree + 2});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 3) {
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree + 3});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]]};
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 3});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
+ auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 3});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{degree + 3});
+ P0_function[cell_id] = sum(cell_id, T, qf);
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+
+ return P0_function;
+}
+
+template <MeshConcept MeshType, typename DataType, size_t NbComponents>
+DiscreteFunctionP0Vector<std::remove_const_t<DataType>>
+exact_projection(
+ const MeshType& mesh,
+ size_t degree,
+ const std::array<std::function<DataType(const TinyVector<MeshType::Dimension>&)>, NbComponents>& vector_exact)
+{
+ DiscreteFunctionP0Vector<std::remove_const_t<DataType>> P0_function_vector{mesh.meshVariant(), vector_exact.size()};
+
+ for (size_t i_component = 0; i_component < vector_exact.size(); ++i_component) {
+ auto exact_function = vector_exact[i_component];
+
+ DiscreteFunctionP0 P0_function = exact_projection(mesh, degree, vector_exact[i_component]);
+
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { P0_function_vector[cell_id][i_component] = P0_function[cell_id]; });
+ }
+
+ return P0_function_vector;
+}
+
+template <typename DataType>
+PUGS_INLINE double
+get_max_error(const DataType& x, const DataType& y)
+{
+ if constexpr (is_tiny_matrix_v<DataType>) {
+ return frobeniusNorm(x - y);
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ return l2Norm(x - y);
+ } else {
+ static_assert(std::is_arithmetic_v<DataType>, "expecting arithmetic type");
+ return std::abs(x - y);
+ }
+}
+
+template <MeshConcept MeshType, typename DataType>
+double
+max_reconstruction_error(const MeshType& mesh,
+ DiscreteFunctionDPk<MeshType::Dimension, const DataType> dpk_f,
+ std::function<DataType(const TinyVector<MeshType::Dimension>&)> exact)
+{
+ auto xr = mesh.xr();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ auto cell_type = mesh.connectivity().cellType();
+
+ double max_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ if constexpr (MeshType::Dimension == 1) {
+ Assert(cell_type[cell_id] == CellType::Line);
+ LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 3) {
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]]};
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
+ auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ }
+ }
+ return max_error;
+}
+
+template <MeshConcept MeshType, typename DataType, size_t NbComponents>
+double
+max_reconstruction_error(
+ const MeshType& mesh,
+ DiscreteFunctionDPkVector<MeshType::Dimension, const DataType> dpk_v,
+ const std::array<std::function<DataType(const TinyVector<MeshType::Dimension>&)>, NbComponents>& vector_exact)
+{
+ auto xr = mesh.xr();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ double max_error = 0;
+ auto cell_type = mesh.connectivity().cellType();
+
+ REQUIRE(NbComponents == dpk_v.numberOfComponents());
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ if constexpr (MeshType::Dimension == 1) {
+ Assert(cell_type[cell_id] == CellType::Line);
+ LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error = std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 3) {
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]]};
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
+ auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+ xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ }
+ }
+ return max_error;
+}
+
+} // namespace test_only
+
+#endif // DISCRETE_FUNCTION_DPK_FOR_TESTS_HPP
diff --git a/tests/test_PolynomialReconstruction_degree_1.cpp b/tests/test_PolynomialReconstruction_degree_1.cpp
index 4e8b8ca57d4266acb26481e91c68c83d6139c951..fb2c7b2ec4a56239e37314704b2cc30122182e8c 100644
--- a/tests/test_PolynomialReconstruction_degree_1.cpp
+++ b/tests/test_PolynomialReconstruction_degree_1.cpp
@@ -2,279 +2,30 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
-#include <Kokkos_Core.hpp>
-
#include <utils/PugsAssert.hpp>
-#include <utils/Types.hpp>
-
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/SmallVector.hpp>
-#include <analysis/GaussQuadratureDescriptor.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
-#include <geometry/CubeTransformation.hpp>
-#include <geometry/LineTransformation.hpp>
-#include <geometry/PrismTransformation.hpp>
-#include <geometry/PyramidTransformation.hpp>
-#include <geometry/SquareTransformation.hpp>
-#include <geometry/TetrahedronTransformation.hpp>
-#include <geometry/TriangleTransformation.hpp>
+
#include <mesh/Mesh.hpp>
-#include <mesh/MeshDataManager.hpp>
#include <mesh/NamedBoundaryDescriptor.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/PolynomialReconstruction.hpp>
+#include <DiscreteFunctionDPkForTests.hpp>
#include <MeshDataBaseForTests.hpp>
-#include <type_traits>
-
// clazy:excludeall=non-pod-global-static
-namespace test_only
-{
-
-template <typename DataType>
-PUGS_INLINE double
-get_max_error(const DataType& x, const DataType& y)
-{
- if constexpr (is_tiny_matrix_v<DataType>) {
- return frobeniusNorm(x - y);
- } else if constexpr (is_tiny_vector_v<DataType>) {
- return l2Norm(x - y);
- } else {
- static_assert(std::is_arithmetic_v<DataType>, "expecting arithmetic type");
- return std::abs(x - y);
- }
-}
-
-template <MeshConcept MeshType, typename DataType>
-double
-max_reconstruction_error(const MeshType& mesh,
- DiscreteFunctionDPk<MeshType::Dimension, const DataType> dpk_f,
- std::function<DataType(const TinyVector<MeshType::Dimension>&)> exact)
-{
- auto xr = mesh.xr();
- auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- double max_error = 0;
- auto cell_type = mesh.connectivity().cellType();
-
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto cell_nodes = cell_to_node_matrix[cell_id];
- if constexpr (MeshType::Dimension == 1) {
- Assert(cell_type[cell_id] == CellType::Line);
- LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
- auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- } else if constexpr (MeshType::Dimension == 2) {
- switch (cell_type[cell_id]) {
- case CellType::Triangle: {
- TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
- auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- case CellType::Quadrangle: {
- SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
- xr[cell_nodes[3]]};
- auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- default: {
- throw UnexpectedError("unexpected cell type");
- }
- }
- } else if constexpr (MeshType::Dimension == 3) {
- switch (cell_type[cell_id]) {
- case CellType::Tetrahedron: {
- TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
- auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- case CellType::Pyramid: {
- PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
- xr[cell_nodes[4]]};
- auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- case CellType::Prism: {
- PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
- xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
- auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- case CellType::Hexahedron: {
- CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
- xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
- auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
- }
- break;
- }
- default: {
- throw UnexpectedError("unexpected cell type");
- }
- }
- }
- }
- return max_error;
-}
-
-template <MeshConcept MeshType, typename DataType, size_t NbComponents>
-double
-max_reconstruction_error(
- const MeshType& mesh,
- DiscreteFunctionDPkVector<MeshType::Dimension, const DataType> dpk_v,
- const std::array<std::function<DataType(const TinyVector<MeshType::Dimension>&)>, NbComponents>& vector_exact)
-{
- auto xr = mesh.xr();
- auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- double max_error = 0;
- auto cell_type = mesh.connectivity().cellType();
-
- REQUIRE(NbComponents == dpk_v.numberOfComponents());
-
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto cell_nodes = cell_to_node_matrix[cell_id];
- if constexpr (MeshType::Dimension == 1) {
- Assert(cell_type[cell_id] == CellType::Line);
- LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
- auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error = std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- } else if constexpr (MeshType::Dimension == 2) {
- switch (cell_type[cell_id]) {
- case CellType::Triangle: {
- TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
- auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- case CellType::Quadrangle: {
- SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
- xr[cell_nodes[3]]};
- auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("unexpected cell type");
- }
- }
- } else if constexpr (MeshType::Dimension == 3) {
- switch (cell_type[cell_id]) {
- case CellType::Tetrahedron: {
- TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
- auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- case CellType::Pyramid: {
- PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
- xr[cell_nodes[4]]};
- auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- case CellType::Prism: {
- PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
- xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
- auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- case CellType::Hexahedron: {
- CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
- xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
- auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- auto x = T(qf.point(i_quadrarture));
- for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
- max_error =
- std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("unexpected cell type");
- }
- }
- }
- }
- return max_error;
-}
-
-} // namespace test_only
-
TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
+ constexpr size_t degree = 1;
+
SECTION("without symmetries")
{
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1}};
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
for (auto descriptor : descriptor_list) {
SECTION(name(descriptor.integrationMethodType()))
@@ -292,12 +43,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
@@ -324,12 +71,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
+1.4 - 0.6 * x[0], //
-0.2 + 3.1 * x[0]};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R3> uh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
@@ -358,12 +101,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
-4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R3x3> Ah{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_exact(xj[cell_id]); });
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
@@ -387,16 +126,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
[](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
[](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- for (size_t i = 0; i < vector_exact.size(); ++i) {
- Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
- }
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
@@ -418,18 +148,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
std::array<std::function<R3(const R1&)>, 2> vector_exact =
- {[](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]}; },
- [](const R1& x) -> R3 { return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]}; }};
+ {[](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+ },
+ [](const R1& x) -> R3 {
+ return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
+ }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
- Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
@@ -473,27 +199,10 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
[](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
[](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
-
- DiscreteFunctionP0<R3> uh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
-
- DiscreteFunctionP0<R3x3> Ah{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_exact(xj[cell_id]); });
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- for (size_t i = 0; i < vector_exact.size(); ++i) {
- Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
- }
- });
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions =
PolynomialReconstruction{descriptor}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
@@ -541,12 +250,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R_exact = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<double> fh{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
@@ -573,12 +278,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
+1.4 - 0.6 * x[0] + 1.3 * x[1], //
-0.2 + 3.1 * x[0] - 1.1 * x[1]};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R3> uh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
@@ -604,12 +305,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
+1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R2x2> Ah{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_exact(xj[cell_id]); });
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
@@ -634,16 +331,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
[](const R2& x) -> double { return +1.4 - 0.6 * x[0] - 2.1 * x[1]; },
[](const R2& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1]; }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- for (size_t i = 0; i < vector_exact.size(); ++i) {
- Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
- }
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
@@ -668,19 +356,15 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R_exact = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -698,12 +382,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
+1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], //
-0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R3> uh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
@@ -730,12 +410,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
//
+2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R2x2> Ah{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_exact(xj[cell_id]); });
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
descriptor.setRowWeighting(false);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
@@ -761,16 +437,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
[](const R3& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2]; },
[](const R3& x) -> double { return -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]; }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- for (size_t i = 0; i < vector_exact.size(); ++i) {
- Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
- }
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
descriptor.setPreconditioning(false);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
@@ -807,7 +474,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
- PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
@@ -844,7 +511,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
- PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
@@ -881,7 +548,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
- PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
@@ -918,13 +585,13 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
SECTION("1D")
{
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
@@ -940,12 +607,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R1> fh{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R1_exact(xj[cell_id]); });
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
@@ -967,15 +630,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
= {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; },
[](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<R1> Vh{p_mesh, 2};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
- Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
@@ -993,21 +648,21 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
SECTION("2D")
{
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
"XMAX"),
std::make_shared<NamedBoundaryDescriptor>(
"YMAX")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
"XMAX"),
std::make_shared<NamedBoundaryDescriptor>(
"YMAX")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
@@ -1026,18 +681,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R2> fh{p_mesh};
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact));
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_exact(xj[cell_id]); });
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
- auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
-
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
- double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R2_exact));
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
@@ -1047,18 +698,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
std::array<std::function<R2(const R2&)>, 2> vector_exact //
- = {[](const R2& x) -> R2 { return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)}; },
- [](const R2& x) -> R2 { return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)}; }};
-
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ = {[](const R2& x) -> R2 {
+ return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)};
+ },
+ [](const R2& x) -> R2 {
+ return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)};
+ }};
- DiscreteFunctionP0Vector<R2> Vh{p_mesh, 2};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
- Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
- });
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
@@ -1074,7 +721,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
SECTION("3D")
{
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
@@ -1083,7 +730,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
"YMAX"),
std::make_shared<NamedBoundaryDescriptor>(
"ZMAX")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
@@ -1092,7 +739,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
"YMAX"),
std::make_shared<NamedBoundaryDescriptor>(
"ZMAX")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
@@ -1113,18 +760,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)}; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R3> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R3_exact(xj[cell_id]); });
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
- auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
- double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R3_exact));
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
@@ -1134,18 +777,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh;
std::array<std::function<R3(const R3&)>, 2> vector_exact //
- = {[](const R3& x) -> R3 { return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)}; },
- [](const R3& x) -> R3 { return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)}; }};
-
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
- Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
- });
+ = {[](const R3& x) -> R3 {
+ return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)};
+ },
+ [](const R3& x) -> R3 {
+ return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)};
+ }};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);