diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp
index c9f7def71f9d342cca580debf09aa8e77316c556..4576c0b64904ecb489caeb456b10271be2a076cc 100644
--- a/src/geometry/LineTransformation.hpp
+++ b/src/geometry/LineTransformation.hpp
@@ -24,12 +24,20 @@ class LineTransformation<1>
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
double
jacobianDeterminant() const
{
return m_jacobian;
}
+ PUGS_INLINE
+ double
+ jacobianDeterminant(const TinyVector<1>&) const
+ {
+ return m_jacobian;
+ }
+
PUGS_INLINE
LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
{
diff --git a/src/geometry/TetrahedronTransformation.hpp b/src/geometry/TetrahedronTransformation.hpp
index 64244266957a350e8be8f78c2eab3d5fead629e1..a81a344bfbabd0999e45b1bcab0b52519f9db4ed 100644
--- a/src/geometry/TetrahedronTransformation.hpp
+++ b/src/geometry/TetrahedronTransformation.hpp
@@ -24,12 +24,20 @@ class TetrahedronTransformation
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
double
jacobianDeterminant() const
{
return m_jacobian_determinant;
}
+ PUGS_INLINE
+ double
+ jacobianDeterminant(const TinyVector<3>&) const
+ {
+ return m_jacobian_determinant;
+ }
+
PUGS_INLINE
TetrahedronTransformation(const TinyVector<Dimension>& a,
const TinyVector<Dimension>& b,
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index 9210e0bda55bb19dfdfd877e509801fa87f0a4f0..df8444d894a7ed56c0fe3825ffbf1cc14a742531 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -26,12 +26,20 @@ class TriangleTransformation<2>
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
double
jacobianDeterminant() const
{
return m_jacobian_determinant;
}
+ PUGS_INLINE
+ double
+ jacobianDeterminant(const TinyVector<2>&) const
+ {
+ return m_jacobian_determinant;
+ }
+
PUGS_INLINE
TriangleTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b, const TinyVector<Dimension>& c)
{
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 20fc3bf61dd6a59ef82ff7bdb25dc5e32c7eef24..de7501d44b5b0aef6e50f830460e9be110456331 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -447,8 +447,7 @@ StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connec
}
Array<uint32_t> column_indices{column_indices_vector.size()};
- parallel_for(
- column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
+ parallel_for(column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
ConnectivityMatrix primal_stencil{row_map, column_indices};
@@ -743,8 +742,7 @@ StencilBuilder::_build_for_different_source_and_target(
}
Array<uint32_t> column_indices{column_indices_vector.size()};
- parallel_for(
- column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
+ parallel_for(column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
ConnectivityMatrix primal_stencil{row_map, column_indices};
@@ -836,7 +834,7 @@ StencilBuilder::buildC2C(const IConnectivity& connectivity,
if ((parallel::size() > 1) and
(stencil_descriptor.numberOfLayers() > GlobalVariableManager::instance().getNumberOfGhostLayers())) {
std::ostringstream error_msg;
- error_msg << "Stencil builder requires" << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
+ error_msg << "Stencil builder requires " << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
<< rang::fg::reset << " layers while parallel number of ghost layer is "
<< GlobalVariableManager::instance().getNumberOfGhostLayers() << ".\n";
error_msg << "Increase the number of ghost layers (using the '--number-of-ghost-layers' option).";
@@ -870,7 +868,7 @@ StencilBuilder::buildN2C(const IConnectivity& connectivity,
if ((parallel::size() > 1) and
(stencil_descriptor.numberOfLayers() > GlobalVariableManager::instance().getNumberOfGhostLayers())) {
std::ostringstream error_msg;
- error_msg << "Stencil builder requires" << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
+ error_msg << "Stencil builder requires " << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
<< rang::fg::reset << " layers while parallel number of ghost layer is "
<< GlobalVariableManager::instance().getNumberOfGhostLayers() << ".\n";
error_msg << "Increase the number of ghost layers (using the '--number-of-ghost-layers' option).";
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d7ba065480be4aac1949b813df03fb34d94d3f56..e3739abe7afd3a625b6510d9be792bdb3b3a76dd 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -280,6 +280,7 @@ add_executable (mpi_unit_tests
test_ParallelChecker_read.cpp
test_Partitioner.cpp
test_PolynomialReconstruction_degree_1.cpp
+ test_PolynomialReconstruction_degree_2.cpp
test_PolynomialReconstructionDescriptor.cpp
test_RandomEngine.cpp
test_StencilBuilder_cell2cell.cpp
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/MeshDataBaseForTests.cpp b/tests/MeshDataBaseForTests.cpp
index 1f2b485edcb7f7b3fadf35d845b59f23c1a03274..1615673d0ba28d7328ab5faef7830a774d5c75f4 100644
--- a/tests/MeshDataBaseForTests.cpp
+++ b/tests/MeshDataBaseForTests.cpp
@@ -1,11 +1,15 @@
#include <MeshDataBaseForTests.hpp>
+
#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/MeshVariant.hpp>
+#include <utils/GlobalVariableManager.hpp>
#include <utils/Messenger.hpp>
#include <utils/PugsAssert.hpp>
+#include <NbGhostLayersTester.hpp>
+
#include <filesystem>
#include <fstream>
@@ -13,17 +17,22 @@ const MeshDataBaseForTests* MeshDataBaseForTests::m_instance = nullptr;
MeshDataBaseForTests::MeshDataBaseForTests()
{
- m_cartesian_1d_mesh = CartesianMeshBuilder{TinyVector<1>{-1}, TinyVector<1>{3}, TinyVector<1, size_t>{23}}.mesh();
+ for (size_t nb_ghost_layers = 1; nb_ghost_layers <= m_max_nb_ghost_layers; ++nb_ghost_layers) {
+ NbGhostLayersTester t{nb_ghost_layers};
- m_cartesian_2d_mesh =
- CartesianMeshBuilder{TinyVector<2>{0, -1}, TinyVector<2>{3, 2}, TinyVector<2, size_t>{6, 7}}.mesh();
+ m_cartesian_1d_mesh.push_back(
+ CartesianMeshBuilder{TinyVector<1>{-1}, TinyVector<1>{3}, TinyVector<1, size_t>{23}}.mesh());
- m_cartesian_3d_mesh =
- CartesianMeshBuilder{TinyVector<3>{0, 1, 0}, TinyVector<3>{2, -1, 3}, TinyVector<3, size_t>{6, 7, 4}}.mesh();
+ m_cartesian_2d_mesh.push_back(
+ CartesianMeshBuilder{TinyVector<2>{0, -1}, TinyVector<2>{3, 2}, TinyVector<2, size_t>{6, 7}}.mesh());
- m_unordered_1d_mesh = _buildUnordered1dMesh();
- m_hybrid_2d_mesh = _buildHybrid2dMesh();
- m_hybrid_3d_mesh = _buildHybrid3dMesh();
+ m_cartesian_3d_mesh.push_back(
+ CartesianMeshBuilder{TinyVector<3>{0, 1, 0}, TinyVector<3>{2, -1, 3}, TinyVector<3, size_t>{6, 7, 4}}.mesh());
+
+ m_unordered_1d_mesh.push_back(_buildUnordered1dMesh());
+ m_hybrid_2d_mesh.push_back(_buildHybrid2dMesh());
+ m_hybrid_3d_mesh.push_back(_buildHybrid3dMesh());
+ }
}
const MeshDataBaseForTests&
@@ -50,37 +59,49 @@ MeshDataBaseForTests::destroy()
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::cartesian1DMesh() const
{
- return m_cartesian_1d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_cartesian_1d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::cartesian2DMesh() const
{
- return m_cartesian_2d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_cartesian_2d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::cartesian3DMesh() const
{
- return m_cartesian_3d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_cartesian_3d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::unordered1DMesh() const
{
- return m_unordered_1d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_unordered_1d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::hybrid2DMesh() const
{
- return m_hybrid_2d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_hybrid_2d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
MeshDataBaseForTests::hybrid3DMesh() const
{
- return m_hybrid_3d_mesh;
+ const size_t nb_ghost_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+ Assert((nb_ghost_layers >= 1) and (nb_ghost_layers <= m_max_nb_ghost_layers));
+ return m_hybrid_3d_mesh[nb_ghost_layers - 1];
}
std::shared_ptr<const MeshVariant>
diff --git a/tests/MeshDataBaseForTests.hpp b/tests/MeshDataBaseForTests.hpp
index 2b092d544271fbd1b1d0488ed73c9c5a401d548b..f0ab29fa9a927c01f16f94925fba7eb750ffd591 100644
--- a/tests/MeshDataBaseForTests.hpp
+++ b/tests/MeshDataBaseForTests.hpp
@@ -4,6 +4,7 @@
#include <array>
#include <memory>
#include <string>
+#include <vector>
class MeshVariant;
@@ -35,15 +36,17 @@ class MeshDataBaseForTests
private:
explicit MeshDataBaseForTests();
+ constexpr static size_t m_max_nb_ghost_layers = 3;
+
static const MeshDataBaseForTests* m_instance;
- std::shared_ptr<const MeshVariant> m_cartesian_1d_mesh;
- std::shared_ptr<const MeshVariant> m_cartesian_2d_mesh;
- std::shared_ptr<const MeshVariant> m_cartesian_3d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_cartesian_1d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_cartesian_2d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_cartesian_3d_mesh;
- std::shared_ptr<const MeshVariant> m_unordered_1d_mesh;
- std::shared_ptr<const MeshVariant> m_hybrid_2d_mesh;
- std::shared_ptr<const MeshVariant> m_hybrid_3d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_unordered_1d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_hybrid_2d_mesh;
+ std::vector<std::shared_ptr<const MeshVariant>> m_hybrid_3d_mesh;
std::shared_ptr<const MeshVariant> _buildUnordered1dMesh();
std::shared_ptr<const MeshVariant> _buildHybrid2dMesh();
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);
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dce627cabab012e398a271b88aafc5131e4338d5
--- /dev/null
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -0,0 +1,716 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Mesh.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 <NbGhostLayersTester.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
+{
+ constexpr size_t degree = 2;
+
+ constexpr size_t nb_ghost_layers = 2;
+ NbGhostLayersTester t{nb_ghost_layers};
+
+ SECTION("without symmetries")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0] - 1.4 * x[0] * x[0]; };
+
+ 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<1, 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));
+ }
+ }
+ }
+
+ SECTION("R^3 data")
+ {
+ using R3 = TinyVector<3>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto R3_exact = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] - x[0] * x[0], //
+ +1.4 - 0.6 * x[0] + 2 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] + 1.4 * x[0] * x[0]};
+ };
+
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
+
+ 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-14));
+ }
+ }
+ }
+
+ SECTION("R^3x3 data")
+ {
+ using R3x3 = TinyMatrix<3, 3>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto R3x3_exact = [](const R1& x) -> R3x3 {
+ return R3x3{+2.3 + 1.7 * x[0] - 2.3 * x[0] * x[0], //
+ -1.7 + 2.1 * x[0] + 1.2 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] - 2.0 * x[0] * x[0], //
+ //
+ +2.4 - 2.3 * x[0] + 1.1 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] - 0.7 * x[0] * x[0], //
+ -3.2 - 3.6 * x[0] + 0.1 * x[0] * x[0], //
+ //
+ -4.1 + 3.1 * x[0] - 0.2 * x[0] * x[0], //
+ +0.8 + 2.9 * x[0] + 4.1 * x[0] * x[0], //
+ -1.6 + 2.3 * x[0] - 1.7 * x[0] * x[0]};
+ };
+
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+
+ SECTION("R vector data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+ std::array<std::function<double(const R1&)>, 3> vector_exact =
+ {[](const R1& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[0] * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0] + 2.1 * x[0] * x[0]; },
+ [](const R1& x) -> double { return +1.4 - 0.6 * x[0] - 1.3 * x[0] * x[0]; }};
+
+ 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>>();
+
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+
+ SECTION("R3 vector data")
+ {
+ using R3 = TinyVector<3>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ 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] + 0.8 * x[0] * x[0], //
+ -1.7 + 2.1 * x[0] - 0.7 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] + 1.9 * x[0] * 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]}; }};
+
+ 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 R3>>();
+
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+
+ SECTION("list of various types")
+ {
+ using R3x3 = TinyMatrix<3>;
+ using R3 = TinyVector<3>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+
+ auto R3_exact = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] + 2.3 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] - 1.5 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] + 2.9 * x[0] * x[0]};
+ };
+
+ auto R3x3_exact = [](const R1& x) -> R3x3 {
+ return R3x3{
+ +2.3 + 1.7 * x[0] - 0.3 * x[0] * x[0],
+ -1.7 + 2.1 * x[0] + 1.4 * x[0] * x[0],
+ +1.4 - 0.6 * x[0] - 2.3 * x[0] * x[0],
+ //
+ +2.4 - 2.3 * x[0] + 1.8 * x[0] * x[0],
+ -0.2 + 3.1 * x[0] - 1.7 * x[0] * x[0],
+ -3.2 - 3.6 * x[0] + 0.7 * x[0] * x[0],
+ //
+ -4.1 + 3.1 * x[0] - 1.9 * x[0] * x[0],
+ +0.8 + 2.9 * x[0] + 2.2 * x[0] * x[0],
+ -1.6 + 2.3 * x[0] - 1.3 * x[0] * x[0],
+ };
+ };
+
+ std::array<std::function<double(const R1&)>, 3> vector_exact =
+ {[](const R1& x) -> double { return +2.3 + 1.7 * x[0] + 2.2 * x[0] * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0] - 1.9 * x[0] * x[0]; },
+ [](const R1& x) -> double { return +1.4 - 0.6 * x[0] + 3.1 * x[0] * x[0]; }};
+
+ 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,
+ std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
+ DiscreteFunctionVariant(Vh));
+
+ {
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, 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));
+ }
+
+ {
+ auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
+ 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-14));
+ }
+
+ {
+ auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ }
+
+ {
+ auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ auto R_exact = [](const R2& x) {
+ return 2.3 + 1.7 * x[0] - 1.3 * x[1] //
+ + 1.2 * x[0] * x[0] + 1.3 * x[0] * x[1] - 3.2 * x[1] * x[1];
+ };
+
+ 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<2, 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-12));
+ }
+ }
+ }
+
+ SECTION("R^3 data")
+ {
+ using R3 = TinyVector<3>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ auto R3_exact = [](const R2& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] //
+ - 2.1 * x[0] * x[0] - 2.3 * x[0] * x[1] - 3.2 * x[1] * x[1], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1] //
+ + 2.3 * x[0] * x[0] - 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] //
+ - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] - 1.1 * x[1] * x[1]};
+ };
+
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
+
+ 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));
+ }
+ }
+ }
+
+ SECTION("R^2x2 data")
+ {
+ using R2x2 = TinyMatrix<2, 2>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ auto R2x2_exact = [](const R2& x) -> R2x2 {
+ return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] //
+ - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
+ -1.7 + 2.1 * x[0] - 2.2 * x[1] //
+ - 1.2 * x[0] * x[0] + 2.1 * x[0] * x[1] - 1.3 * x[1] * x[1], //
+ +1.4 - 0.6 * x[0] - 2.1 * x[1] //
+ - 1.1 * x[0] * x[0] - 2.3 * x[0] * x[1] + 2.1 * x[1] * x[1],
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] //
+ + 2.7 * x[0] * x[0] + 2.1 * x[0] * x[1] - 2.7 * x[1] * x[1]};
+ };
+
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+
+ SECTION("vector data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ std::array<std::function<double(const R2&)>, 4> vector_exact =
+ {[](const R2& x) -> double {
+ return +2.3 + 1.7 * x[0] + 1.2 * x[1] + 1.2 * x[0] * x[0] - 2.1 * x[0] * x[1] + 3.1 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 0.7 * x[0] * x[0] + 2.2 * x[0] * x[1] - 1.6 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return +1.4 - 0.6 * x[0] - 2.1 * x[1] + 2.3 * x[0] * x[0] + 2.3 * x[0] * x[1] - 2.9 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return +2.4 - 2.3 * x[0] + 1.3 * x[1] - 2.7 * x[0] * x[0] - 1.2 * x[0] * x[1] - 0.7 * x[1] * x[1];
+ }};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+ 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] //
+ + 1.7 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
+ };
+
+ 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-12));
+ }
+ }
+ }
+
+ SECTION("R^3 data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+
+ auto R3_exact = [](const R3& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2] //
+ + 1.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2] //
+ - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2] //
+ + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2] //
+ - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2] //
+ - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2]};
+ };
+
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+
+ 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));
+ }
+ }
+ }
+
+ SECTION("R^2x2 data")
+ {
+ using R2x2 = TinyMatrix<2, 2>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+
+ auto R2x2_exact = [](const R3& x) -> R2x2 {
+ return R2x2{
+ +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2] //
+ - 1.7 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 1.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
+ -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2] //
+ + 3.7 * x[0] * x[0] + 1.3 * x[1] * x[1] + 1.6 * x[2] * x[2] //
+ - 2.1 * x[0] * x[1] - 1.5 * x[0] * x[2] - 1.7 * x[1] * x[2], //
+ //
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2] //
+ - 2.1 * x[0] * x[0] + 1.7 * x[1] * x[1] + 1.8 * x[2] * x[2] //
+ - 1.4 * x[0] * x[1] + 1.3 * x[0] * x[2] - 2.9 * x[1] * x[2], //
+ -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2] //
+ + 1.6 * x[0] * x[0] + 2.1 * x[1] * x[1] - 2.1 * x[2] * x[2] //
+ - 1.1 * x[0] * x[1] - 1.3 * x[0] * x[2] + 1.6 * x[1] * x[2], //
+ };
+ };
+
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+
+ descriptor.setRowWeighting(false);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+
+ SECTION("vector data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+#warning continue here
+ std::array<std::function<double(const R3&)>, 4> vector_exact =
+ {[](const R3& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2]; },
+ [](const R3& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2]; },
+ [](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]; }};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ descriptor.setPreconditioning(false);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("with symmetries")
+ {
+ // SECTION("1D")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
+
+ // using R1 = TinyVector<1>;
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R^1 data")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ // auto& mesh = *p_mesh;
+
+ // auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
+
+ // DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ // auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("R1 vector data")
+ // {
+ // for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ // SECTION(named_mesh.name())
+ // {
+ // auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ // auto& mesh = *p_mesh;
+
+ // std::array<std::function<R1(const R1&)>, 2> vector_exact //
+ // = {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; },
+ // [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};
+
+ // 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 R1>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+
+ // SECTION("2D")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX")}}};
+
+ // using R2 = TinyVector<2>;
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R^2 data")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ // auto& mesh = *p_mesh;
+
+ // auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
+
+ // 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));
+ // }
+
+ // SECTION("vector of R2")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ // 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)};
+ // }};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+
+ // SECTION("3D")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "ZMAX")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "ZMAX")}}};
+
+ // using R3 = TinyVector<3>;
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R^3 data")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ // 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)};
+ // };
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+
+ // 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));
+ // }
+
+ // SECTION("vector of R3")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ // 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)};
+ // }};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+ }
+}