diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 44313040017f24ee7d79a8dec35697b3ebae1bc7..34015395ae6d57c5e0c2850333ad90ef65211ef1 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -1,6 +1,11 @@
#include <mesh/Connectivity.hpp>
#include <mesh/StencilBuilder.hpp>
+#include <utils/Messenger.hpp>
+
+#warning remove after rework
+#include <set>
+
template <typename ConnectivityType>
Array<const uint32_t>
StencilBuilder::_getRowMap(const ConnectivityType& connectivity) const
@@ -104,26 +109,91 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
template <typename ConnectivityType>
Stencil
-StencilBuilder::_build(const ConnectivityType& connectivity) const
+StencilBuilder::_build(const ConnectivityType& connectivity, size_t number_of_layers) const
{
- Array<const uint32_t> row_map = this->_getRowMap(connectivity);
- Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
+ if (number_of_layers == 1) {
+ Array<const uint32_t> row_map = this->_getRowMap(connectivity);
+ Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
+
+ return ConnectivityMatrix{row_map, column_indices};
+ } else {
+ if (number_of_layers > 2) {
+ throw NotImplementedError("number of layers too large");
+ }
+ if (parallel::size() > 1) {
+ throw NotImplementedError("stencils with more than 1 layers are not defined in parallel");
+ }
+
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
+
+ Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+ row_map[0] = 0;
+
+ std::vector<CellId> column_indices_vector;
- return ConnectivityMatrix{row_map, column_indices};
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ if (cell_is_owned[cell_id]) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (size_t i_node_1 = 0; i_node_1 < cell_to_node_matrix[cell_id].size(); ++i_node_1) {
+ const NodeId layer_1_node_id = cell_to_node_matrix[cell_id][i_node_1];
+
+ for (size_t i_cell_1 = 0; i_cell_1 < node_to_cell_matrix[layer_1_node_id].size(); ++i_cell_1) {
+ CellId cell_1_id = node_to_cell_matrix[layer_1_node_id][i_cell_1];
+
+ for (size_t i_node_2 = 0; i_node_2 < cell_to_node_matrix[cell_1_id].size(); ++i_node_2) {
+ const NodeId layer_2_node_id = cell_to_node_matrix[cell_1_id][i_node_2];
+
+ for (size_t i_cell_2 = 0; i_cell_2 < node_to_cell_matrix[layer_2_node_id].size(); ++i_cell_2) {
+ CellId cell_2_id = node_to_cell_matrix[layer_2_node_id][i_cell_2];
+
+ if (cell_2_id != cell_id) {
+ cell_set.insert(cell_2_id);
+ }
+ }
+ }
+ }
+ }
+
+ for (auto stencil_cell_id : cell_set) {
+ column_indices_vector.push_back(stencil_cell_id);
+ }
+ row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+ }
+ }
+
+ if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
+ throw UnexpectedError("incorrect stencil size");
+ }
+
+ Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
+ column_indices.fill(std::numeric_limits<uint32_t>::max());
+
+ for (size_t i = 0; i < column_indices.size(); ++i) {
+ column_indices[i] = column_indices_vector[i];
+ }
+
+ return ConnectivityMatrix{row_map, column_indices};
+ }
}
Stencil
-StencilBuilder::build(const IConnectivity& connectivity) const
+StencilBuilder::build(const IConnectivity& connectivity, size_t number_of_layers) const
{
switch (connectivity.dimension()) {
case 1: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity));
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers);
}
case 2: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity));
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers);
}
case 3: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity));
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers);
}
default: {
throw UnexpectedError("invalid connectivity dimension");
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 1743a2d719fd920d121b604e61732d9f6450197d..53918e64ac9234c396986f2bfef3f90510f287e0 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -15,10 +15,10 @@ class StencilBuilder
const Array<const uint32_t>& row_map) const;
template <typename ConnectivityType>
- Stencil _build(const ConnectivityType& connectivity) const;
+ Stencil _build(const ConnectivityType& connectivity, size_t number_of_layers) const;
friend class StencilManager;
- Stencil build(const IConnectivity& connectivity) const;
+ Stencil build(const IConnectivity& connectivity, size_t number_of_layers) const;
public:
StencilBuilder() = default;
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
index 58190f002739ca0070ac14c35214598c4476d8d9..120a95e0dfc54a38f9ea7a4ab94bbe4607014981 100644
--- a/src/mesh/StencilManager.cpp
+++ b/src/mesh/StencilManager.cpp
@@ -18,11 +18,11 @@ StencilManager::destroy()
Assert(m_instance != nullptr, "StencilManager was not created");
// LCOV_EXCL_START
- if (m_instance->m_connectivity_id_to_stencil_map.size() > 0) {
+ if (m_instance->m_stored_stencil_map.size() > 0) {
std::stringstream error;
error << ": some connectivities are still registered\n";
- for (const auto& [id, stencil_p] : m_instance->m_connectivity_id_to_stencil_map) {
- error << " - connectivity of id " << rang::fgB::magenta << id << rang::style::reset << '\n';
+ for (const auto& [key, stencil_p] : m_instance->m_stored_stencil_map) {
+ error << " - connectivity of id " << rang::fgB::magenta << key.connectivity_id << rang::style::reset << '\n';
}
throw UnexpectedError(error.str());
// LCOV_EXCL_STOP
@@ -32,14 +32,14 @@ StencilManager::destroy()
}
const Stencil&
-StencilManager::getStencil(const IConnectivity& connectivity)
+StencilManager::getStencil(const IConnectivity& connectivity, size_t degree)
{
- if (not m_connectivity_id_to_stencil_map.contains(connectivity.id())) {
- m_connectivity_id_to_stencil_map[connectivity.id()] =
- std::make_shared<Stencil>(StencilBuilder{}.build(connectivity));
+ if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree})) {
+ m_stored_stencil_map[Key{connectivity.id(), degree}] =
+ std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree));
}
- return *m_connectivity_id_to_stencil_map.at(connectivity.id());
+ return *m_stored_stencil_map.at(Key{connectivity.id(), degree});
}
void
@@ -48,9 +48,9 @@ StencilManager::deleteConnectivity(const size_t connectivity_id)
bool has_removed = false;
do {
has_removed = false;
- for (const auto& [id, p_stencil] : m_connectivity_id_to_stencil_map) {
- if (connectivity_id == id) {
- m_connectivity_id_to_stencil_map.erase(connectivity_id);
+ for (const auto& [key, p_stencil] : m_stored_stencil_map) {
+ if (connectivity_id == key.connectivity_id) {
+ m_stored_stencil_map.erase(key);
has_removed = true;
break;
}
diff --git a/src/mesh/StencilManager.hpp b/src/mesh/StencilManager.hpp
index d1503e60f64eb44e18ed03c228d5c754afcbe248..f2912be1fca197aab6b9e010c600de3d8e507827 100644
--- a/src/mesh/StencilManager.hpp
+++ b/src/mesh/StencilManager.hpp
@@ -15,7 +15,28 @@ class StencilManager
static StencilManager* m_instance;
- std::unordered_map<size_t, std::shared_ptr<const Stencil>> m_connectivity_id_to_stencil_map;
+ struct Key
+ {
+ size_t connectivity_id;
+ size_t degree;
+
+ PUGS_INLINE bool
+ operator==(const Key& k) const
+ {
+ return (connectivity_id == k.connectivity_id) and (degree == k.degree);
+ }
+ };
+ struct HashKey
+ {
+ size_t
+ operator()(const Key& key) const
+ {
+ return (std::hash<decltype(Key::connectivity_id)>()(key.connectivity_id)) ^
+ (std::hash<decltype(Key::degree)>()(key.degree) << 1);
+ }
+ };
+
+ std::unordered_map<Key, std::shared_ptr<const Stencil>, HashKey> m_stored_stencil_map;
public:
static void create();
@@ -31,7 +52,7 @@ class StencilManager
void deleteConnectivity(const size_t connectivity_id);
- const Stencil& getStencil(const IConnectivity& i_connectivity);
+ const Stencil& getStencil(const IConnectivity& i_connectivity, size_t degree = 1);
StencilManager(const StencilManager&) = delete;
StencilManager(StencilManager&&) = delete;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 0d0c829c7ce0d7f91c635dd59ae7636717ab7e2b..270d663c3363471ecc05de12a5d1bc8661e2a19f 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -253,7 +253,7 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
}
for (size_t i_z = 1; i_z <= degree; ++i_z) {
- const size_t begin_i_z_1 = i_z * (i_z + 1) * (i_z + 2) / 6 - 1;
+ const size_t begin_i_z_1 = ((i_z - 1) * (3 * (degree + 2) * (degree + 3) + (i_z - (3 * degree + 8)) * i_z)) / 6;
for (size_t i_y = 0; i_y <= degree - i_z; ++i_y) {
const size_t begin_i_y_1 = (i_y * (2 * degree - i_y + 3)) / 2 + begin_i_z_1;
for (size_t l = 0; l <= degree - i_y - i_z; ++l) {
@@ -550,7 +550,7 @@ PolynomialReconstruction::_build(
const size_t basis_dimension =
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
- const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
+ const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity(), m_descriptor.degree());
auto xr = mesh.xr();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
@@ -604,16 +604,6 @@ PolynomialReconstruction::_build(
inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
}
-#warning remove after rework
- [[maybe_unused]] const int INTEGRATE = []() {
- auto value = std::getenv("INTEGRATE");
- if (value == nullptr) {
- return 0;
- } else {
- return std::atoi(value);
- }
- }();
-
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
@@ -680,7 +670,8 @@ PolynomialReconstruction::_build(
ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
- if ((m_descriptor.degree() == 1) and (INTEGRATE == 0)) {
+ if ((m_descriptor.degree() == 1) and
+ (m_descriptor.integrationMethod() == PolynomialReconstructionDescriptor::IntegrationMethod::cell_center)) {
const Rd& Xj = xj[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
@@ -690,8 +681,12 @@ PolynomialReconstruction::_build(
}
}
- } else if ((INTEGRATE == 2) or (INTEGRATE == 1)) {
- if ((INTEGRATE == 2) and (MeshType::Dimension == 2)) {
+ } else if ((m_descriptor.integrationMethod() ==
+ PolynomialReconstructionDescriptor::IntegrationMethod::element) or
+ (m_descriptor.integrationMethod() ==
+ PolynomialReconstructionDescriptor::IntegrationMethod::boundary)) {
+ if ((m_descriptor.integrationMethod() == PolynomialReconstructionDescriptor::IntegrationMethod::boundary) and
+ (MeshType::Dimension == 2)) {
if constexpr (MeshType::Dimension == 2) {
SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
@@ -817,12 +812,29 @@ PolynomialReconstruction::_build(
dpk_j[0] = p0_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) {
+ if (m_descriptor.degree() > 1) {
+ auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ for (size_t i = 0; i < basis_dimension - 1; ++i) {
+ dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
+ }
+
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
dpk_j_ip1 = X(i, column_begin);
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
+ if (m_descriptor.degree() > 1) {
+ auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ for (size_t i = 0; i < basis_dimension - 1; ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ }
+ }
+ }
+
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::Dimension; ++k) {
@@ -831,6 +843,19 @@ PolynomialReconstruction::_build(
}
column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) {
+ if (m_descriptor.degree() > 1) {
+ auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ for (size_t i = 0; i < basis_dimension - 1; ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ dpk_j_0(k, l) -=
+ X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
+ }
+ }
+ }
+ }
+
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
@@ -860,6 +885,13 @@ PolynomialReconstruction::_build(
const size_t component_begin = l * size;
dpk_j[component_begin] = cell_vector[l];
if constexpr (std::is_arithmetic_v<DataType>) {
+ if (m_descriptor.degree() > 1) {
+ auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ for (size_t i = 0; i < basis_dimension - 1; ++i) {
+ dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
+ }
+
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
dpk_j_ip1 = X(i, column_begin);
diff --git a/src/scheme/PolynomialReconstructionDescriptor.hpp b/src/scheme/PolynomialReconstructionDescriptor.hpp
index a7394500106b2f8a43c5e62f9ee74912b4d8fee5..b561f7080b3c47d4bd2029a0d0a03791e02c4411 100644
--- a/src/scheme/PolynomialReconstructionDescriptor.hpp
+++ b/src/scheme/PolynomialReconstructionDescriptor.hpp
@@ -7,12 +7,31 @@
class PolynomialReconstructionDescriptor
{
+ public:
+ enum class IntegrationMethod
+ {
+ cell_center, // use exact integrals for degree 1 polynomials
+ // using evaluation at mass center
+ element, // use element based quadrature to compute
+ // polynomial integrals
+ boundary // use divergence theorem to compute polynomial
+ // integrals
+ };
+
private:
+ IntegrationMethod m_integration_method;
size_t m_degree;
bool m_preconditioning = true;
bool m_row_weighting = true;
public:
+ PUGS_INLINE
+ IntegrationMethod
+ integrationMethod() const
+ {
+ return m_integration_method;
+ }
+
PUGS_INLINE
size_t
degree() const
@@ -55,7 +74,9 @@ class PolynomialReconstructionDescriptor
m_row_weighting = row_weighting;
}
- PolynomialReconstructionDescriptor(const size_t degree) : m_degree{degree} {}
+ PolynomialReconstructionDescriptor(const IntegrationMethod integration_method, const size_t degree)
+ : m_integration_method{integration_method}, m_degree{degree}
+ {}
PolynomialReconstructionDescriptor() = delete;
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 4ba54bda2dc6a586e14e16fb01afee0417bb2476..bf1c51121d943b37e854e295b5c472140239fc78 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -7,7 +7,7 @@ void
test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
uint64_t degree)
{
- PolynomialReconstructionDescriptor descriptor{degree};
+ PolynomialReconstructionDescriptor descriptor{PolynomialReconstructionDescriptor::IntegrationMethod::element, degree};
{
const size_t nb_loops = 1;
std::cout << "** variable list at once (" << nb_loops << " times)\n";
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b242c06dc5ff8c5d0882b46f962926890b51e433..e2dff14cced997684da1d214bfd20f9df19df664 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -128,6 +128,7 @@ add_executable (unit_tests
test_PugsUtils.cpp
test_PyramidGaussQuadrature.cpp
test_PyramidTransformation.cpp
+ test_QuadraticPolynomialReconstruction.cpp
test_QuadratureType.cpp
test_RefId.cpp
test_RefItemList.cpp
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index b7c7019f57e6227de20f746900d87416e0a400ff..dfa092bb109c2540892e0d921e11683a4be3dfd4 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -24,7 +24,9 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
SECTION("degree 1")
{
- PolynomialReconstructionDescriptor descriptor{1};
+ PolynomialReconstructionDescriptor::IntegrationMethod method =
+ PolynomialReconstructionDescriptor::IntegrationMethod::cell_center;
+ PolynomialReconstructionDescriptor descriptor{method, 1};
SECTION("1D")
{
diff --git a/tests/test_PolynomialReconstructionDescriptor.cpp b/tests/test_PolynomialReconstructionDescriptor.cpp
index 6b7cd2853c783449f8b8f9167c06241a47b00b38..eda18e3d6447461347cf4018f2165700bcb9fdb4 100644
--- a/tests/test_PolynomialReconstructionDescriptor.cpp
+++ b/tests/test_PolynomialReconstructionDescriptor.cpp
@@ -8,9 +8,12 @@
TEST_CASE("PolynomialReconstructionDescriptor", "[scheme]")
{
+#warning tests are not completed
SECTION("degree 1")
{
- PolynomialReconstructionDescriptor descriptor{1};
+ PolynomialReconstructionDescriptor::IntegrationMethod method =
+ PolynomialReconstructionDescriptor::IntegrationMethod::cell_center;
+ PolynomialReconstructionDescriptor descriptor{method, 1};
REQUIRE(descriptor.degree() == 1);
REQUIRE(descriptor.preconditioning() == true);
@@ -19,7 +22,10 @@ TEST_CASE("PolynomialReconstructionDescriptor", "[scheme]")
SECTION("degree 4")
{
- PolynomialReconstructionDescriptor descriptor{4};
+ PolynomialReconstructionDescriptor::IntegrationMethod method =
+ PolynomialReconstructionDescriptor::IntegrationMethod::cell_center;
+
+ PolynomialReconstructionDescriptor descriptor{method, 4};
REQUIRE(descriptor.degree() == 4);
REQUIRE(descriptor.preconditioning() == true);
@@ -28,7 +34,9 @@ TEST_CASE("PolynomialReconstructionDescriptor", "[scheme]")
SECTION("utlities")
{
- PolynomialReconstructionDescriptor descriptor{3};
+ PolynomialReconstructionDescriptor::IntegrationMethod method =
+ PolynomialReconstructionDescriptor::IntegrationMethod::cell_center;
+ PolynomialReconstructionDescriptor descriptor{method, 3};
REQUIRE(descriptor.degree() == 3);
REQUIRE(descriptor.preconditioning() == true);
diff --git a/tests/test_QuadraticPolynomialReconstruction.cpp b/tests/test_QuadraticPolynomialReconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..69c685f544f774205effbcf304dca43ff7d84658
--- /dev/null
+++ b/tests/test_QuadraticPolynomialReconstruction.cpp
@@ -0,0 +1,1198 @@
+#include <catch2/catch_approx.hpp>
+#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/GaussLegendreQuadratureDescriptor.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 <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
+{
+ SECTION("degree 2")
+ {
+ PolynomialReconstructionDescriptor::IntegrationMethod method =
+ PolynomialReconstructionDescriptor::IntegrationMethod::element;
+ PolynomialReconstructionDescriptor descriptor{method, 2};
+
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{2});
+
+ 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 f_exact = [](const R1& x) { return 2.3 + 1.7 * x[0] - 3.2 * x[0] * x[0]; };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ double value = 0;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * f_exact(T(qf.point(i_q))) * T.jacobianDeterminant();
+ }
+
+ fh[cell_id] = value / Vj[cell_id];
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_fh[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ error += std::abs(qf.weight(i_q) * (f_exact(T(qf.point(i_q))) - dpkj(T(qf.point(i_q)))) *
+ T.jacobianDeterminant());
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ 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 u_exact = [](const R1& X) -> R3 {
+ const double x = X[0];
+
+ return R3{2.3 + 1.7 * x - 3.2 * x * x, //
+ 7 + 5 * x + 3 * x * x, //
+ -4 + 3.5 * x + 2.7 * x * x};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<R3> uh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3 value = zero;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * u_exact(T(qf.point(i_q)));
+ }
+
+ uh[cell_id] = 1. / Vj[cell_id] * value;
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_fh[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ R3 diff =
+ qf.weight(i_q) * T.jacobianDeterminant() * (u_exact(T(qf.point(i_q))) - dpkj(T(qf.point(i_q))));
+ error += std::abs(diff[0]) + std::abs(diff[1]) + std::abs(diff[2]);
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ 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 A_exact = [](const R1& X) -> R3x3 {
+ const double x = X[0];
+ const double x2 = x * x;
+ return R3x3{+2.3 + 1.7 * x + 2.9 * x2, //
+ -1.7 + 2.1 * x + 1.7 * x2, //
+ +1.4 - 0.6 * x - 0.5 * x2, //
+ //
+ +2.4 - 2.3 * x + 2.3 * x2, //
+ -0.2 + 3.1 * x - 1.9 * x2, //
+ -3.2 - 3.6 * x - 0.3 * x2, //
+ //
+ -4.1 + 3.1 * x - 1.3 * x2, //
+ +0.8 + 2.9 * x + 2.1 * x2, //
+ -1.6 + 2.3 * x + 0.8 * x2};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ DiscreteFunctionP0<R3x3> Ah{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3x3 value = zero;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * A_exact(T(qf.point(i_q)));
+ }
+
+ Ah[cell_id] = 1. / Vj[cell_id] * value;
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_Ah[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ R3x3 diff =
+ qf.weight(i_q) * T.jacobianDeterminant() * (A_exact(T(qf.point(i_q))) - dpkj(T(qf.point(i_q))));
+ for (size_t i = 0; i < diff.numberOfRows(); ++i) {
+ for (size_t j = 0; j < diff.numberOfColumns(); ++j) {
+ error += std::abs(diff(i, j));
+ }
+ }
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ SECTION("R 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;
+
+ auto u_exact = [](const R1& X) -> R3 {
+ const double x = X[0];
+
+ return R3{2.3 + 1.7 * x - 3.2 * x * x, //
+ 7 + 5 * x + 3 * x * x, //
+ -4 + 3.5 * x + 2.7 * x * x};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0Vector<double> uh{p_mesh, 3};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3 value = zero;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * u_exact(T(qf.point(i_q)));
+ }
+
+ value *= 1. / Vj[cell_id];
+
+ for (size_t i = 0; i < value.dimension(); ++i) {
+ uh[cell_id][i] = value[i];
+ }
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ for (size_t l = 0; l < dpk_uh.numberOfComponents(); ++l) {
+ error += std::abs(qf.weight(i_q) * T.jacobianDeterminant() *
+ (u_exact(T(qf.point(i_q)))[l] - dpk_uh(cell_id, l)(T(qf.point(i_q)))));
+ }
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ 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 f_exact = [](const R1& x) { return 2.3 + 1.7 * x[0] - 3.2 * x[0] * x[0]; };
+
+ auto u_exact = [](const R1& X) -> R3 {
+ const double x = X[0];
+
+ return R3{2.3 + 1.7 * x - 3.2 * x * x, //
+ 7 + 5 * x + 3 * x * x, //
+ -4 + 3.5 * x + 2.7 * x * x};
+ };
+
+ auto A_exact = [](const R1& X) -> R3x3 {
+ const double x = X[0];
+ const double x2 = x * x;
+ return R3x3{+2.3 + 1.7 * x + 2.9 * x2, //
+ -1.7 + 2.1 * x + 1.7 * x2, //
+ +1.4 - 0.6 * x - 0.5 * x2, //
+ //
+ +2.4 - 2.3 * x + 2.3 * x2, //
+ -0.2 + 3.1 * x - 1.9 * x2, //
+ -3.2 - 3.6 * x - 0.3 * x2, //
+ //
+ -4.1 + 3.1 * x - 1.3 * x2, //
+ +0.8 + 2.9 * x + 2.1 * x2, //
+ -1.6 + 2.3 * x + 0.8 * x2};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0Vector<double> vh{p_mesh, 3};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3 value = zero;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * u_exact(T(qf.point(i_q)));
+ }
+
+ value *= 1. / Vj[cell_id];
+
+ for (size_t i = 0; i < value.dimension(); ++i) {
+ vh[cell_id][i] = value[i];
+ }
+ });
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ double value = 0;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * f_exact(T(qf.point(i_q))) * T.jacobianDeterminant();
+ }
+
+ fh[cell_id] = value / Vj[cell_id];
+ });
+
+ DiscreteFunctionP0<R3x3> Ah{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3x3 value = zero;
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * A_exact(T(qf.point(i_q)));
+ }
+
+ Ah[cell_id] = 1. / Vj[cell_id] * value;
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah, vh, fh);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+ auto dpk_vh = reconstructions[1]->get<DiscreteFunctionDPkVector<1, const double>>();
+ auto dpk_fh = reconstructions[2]->get<DiscreteFunctionDPk<1, const double>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ for (size_t l = 0; l < dpk_vh.numberOfComponents(); ++l) {
+ error += std::abs(qf.weight(i_q) * T.jacobianDeterminant() *
+ (u_exact(T(qf.point(i_q)))[l] - dpk_vh(cell_id, l)(T(qf.point(i_q)))));
+ }
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_fh[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ error += std::abs(qf.weight(i_q) * (f_exact(T(qf.point(i_q))) - dpkj(T(qf.point(i_q)))) *
+ T.jacobianDeterminant());
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_Ah[cell_id];
+
+ LineTransformation<1> T{xr[cell_node_list[0]], xr[cell_node_list[1]]};
+
+ double error = 0;
+
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ R3x3 diff =
+ qf.weight(i_q) * T.jacobianDeterminant() * (A_exact(T(qf.point(i_q))) - dpkj(T(qf.point(i_q))));
+ for (size_t i = 0; i < diff.numberOfRows(); ++i) {
+ for (size_t j = 0; j < diff.numberOfColumns(); ++j) {
+ error += std::abs(diff(i, j));
+ }
+ }
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ auto integrate_in_cell = [](const CellType& cell_type, const auto& cell_node_list, const auto& xr,
+ const auto& exact, auto& value) {
+ switch (cell_type) {
+ case CellType::Triangle: {
+ TriangleTransformation<2> T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]]};
+ QuadratureFormula<2> qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{2});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ value += qf.weight(i_q) * T.jacobianDeterminant() * exact(T(qf.point(i_q)));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<2> T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]};
+ QuadratureFormula<2> qf =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R2 x_hat = qf.point(i_q);
+ const R2 x = T(x_hat);
+ value += qf.weight(i_q) * T.jacobianDeterminant(x_hat) * exact(x);
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type");
+ }
+ }
+ };
+
+ auto compute_L1_error = [](const auto& cell_type, const auto& mesh, const auto& exact, const auto& dpk) {
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& xr = mesh.xr();
+
+ using DataType = typename std::decay_t<decltype(dpk)>::data_type;
+
+ auto sum_abs = [](const DataType& diff) -> double {
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ return std::abs(diff);
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ double sum = 0;
+ for (size_t i = 0; i < diff.dimension(); ++i) {
+ sum += std::abs(diff[i]);
+ }
+ return sum;
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ double sum = 0;
+ for (size_t i = 0; i < diff.numberOfRows(); ++i) {
+ for (size_t j = 0; j < diff.numberOfRows(); ++j) {
+ sum += std::abs(diff(i, j));
+ }
+ }
+ return sum;
+ } else {
+ throw UnexpectedError("unexpected value type");
+ }
+ };
+
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk[cell_id];
+
+ double error = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<2> T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]]};
+ QuadratureFormula<2> qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{2});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R2 x = T(qf.point(i_q));
+ DataType diff = qf.weight(i_q) * T.jacobianDeterminant() * (exact(x) - dpkj(x));
+ error += sum_abs(diff);
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<2> T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]};
+ QuadratureFormula<2> qf =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R2 x_hat = qf.point(i_q);
+ const R2 x = T(x_hat);
+ DataType diff = qf.weight(i_q) * T.jacobianDeterminant(x_hat) * (exact(x) - dpkj(x));
+ error += sum_abs(diff);
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type");
+ }
+ }
+
+ L1_error += error;
+ }
+ return L1_error;
+ };
+
+ 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 f_exact = [](const R2& x) {
+ return 2.3 + 1.7 * x[0] + 0.2 * x[1] - 3.2 * x[0] * x[0] + 1.3 * x[0] * x[1] - 1.4 * x[1] * x[1];
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_type = mesh.connectivity().cellType();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ double value = 0;
+ integrate_in_cell(cell_type[cell_id], cell_node_list, xr, f_exact, value);
+ fh[cell_id] = value / Vj[cell_id];
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+
+ double L1_error = compute_L1_error(cell_type, mesh, f_exact, dpk_fh);
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+
+ 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 u_exact = [](const R2& X) -> R3 {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+
+ return R3{+2.3 + 1.7 * x + 1.3 * y - 3.2 * x2 + 1.6 * xy + 0.9 * y2, //
+ +7.0 + 5.0 * x - 2.4 * y + 3.0 * x2 - 0.8 * xy + 1.1 * y2, //
+ -4.0 + 3.5 * x - 0.7 * y + 2.7 * x2 - 2.1 * xy - 1.8 * y2};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_type = mesh.connectivity().cellType();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<R3> uh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R3 value = zero;
+ integrate_in_cell(cell_type[cell_id], cell_node_list, xr, u_exact, value);
+ uh[cell_id] = 1. / Vj[cell_id] * value;
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
+
+ double L1_error = compute_L1_error(cell_type, mesh, u_exact, dpk_uh);
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+
+ 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 A_exact = [](const R2& X) -> R2x2 {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+
+ return R2x2{+2.3 + 1.7 * x + 1.2 * y + 1.1 * x2 + 0.7 * xy - 0.8 * y2, //
+ -1.7 + 2.1 * x - 2.2 * y - 1.9 * x2 + 0.2 * xy + 1.2 * y2, //
+ +1.4 - 0.6 * x - 2.1 * y + 0.3 * x2 - 3.1 * xy + 1.3 * y2, //
+ +2.4 - 2.3 * x + 1.3 * y - 0.8 * x2 + 1.2 * xy - 2.0 * y2};
+ };
+
+ auto xr = mesh.xr();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_type = mesh.connectivity().cellType();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<R2x2> Ah{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ R2x2 value = zero;
+ integrate_in_cell(cell_type[cell_id], cell_node_list, xr, A_exact, value);
+ Ah[cell_id] = 1. / Vj[cell_id] * value;
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+
+ double L1_error = compute_L1_error(cell_type, mesh, A_exact, dpk_Ah);
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+
+ // SECTION("vector data")
+ // {
+ // using R4 = TinyVector<4>;
+
+ // 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 vector_affine = [](const R2& x) -> R4 {
+ // return R4{+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();
+
+ // DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+ // parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ // auto vector = vector_affine(xj[cell_id]);
+ // for (size_t i = 0; i < vector.dimension(); ++i) {
+ // Vh[cell_id][i] = vector[i];
+ // }
+ // });
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+
+ // {
+ // double max_mean_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // auto vector = vector_affine(xj[cell_id]);
+ // for (size_t i = 0; i < vector.dimension(); ++i) {
+ // max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ // }
+
+ // {
+ // double max_x_slope_error = 0;
+ // const R4 slope{+1.7, +2.1, -0.6, -2.3};
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // for (size_t i = 0; i < slope.dimension(); ++i) {
+ // const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
+ // dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
+
+ // max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ // }
+
+ // {
+ // double max_y_slope_error = 0;
+ // const R4 slope{+1.2, -2.2, -2.1, +1.3};
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // for (size_t i = 0; i < slope.dimension(); ++i) {
+ // const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
+ // dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
+
+ // max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ // }
+ // }
+ // }
+ // }
+ }
+
+ 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 f_exact = [](const R3& X) {
+ const double& x = X[0];
+ const double& y = X[1];
+ const double& z = X[2];
+
+ return 2.3 + 1.7 * x + 0.2 * y + 1.4 * z - 3.2 * x * x + 1.3 * x * y - 1.6 * x * z - 1.4 * y * y +
+ 2 * y * z - 1.8 * z * z;
+ };
+
+ auto xr = mesh.xr();
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ auto cell_type = mesh.connectivity().cellType();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ double value = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{2});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ value += qf.weight(i_q) * f_exact(x) * T.jacobianDeterminant();
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ value += qf.weight(i_q) * f_exact(x) * T.jacobianDeterminant(x_hat);
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]]};
+ QuadratureFormula<3> qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ value += qf.weight(i_q) * f_exact(x) * T.jacobianDeterminant(x_hat);
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
+ xr[cell_node_list[6]], xr[cell_node_list[7]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ value += qf.weight(i_q) * f_exact(x) * T.jacobianDeterminant(x_hat);
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type");
+ }
+ }
+
+ fh[cell_id] = value / Vj[cell_id];
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+
+ {
+ double L1_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const auto dpkj = dpk_fh[cell_id];
+
+ double error = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{2});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ error += std::abs(qf.weight(i_q) * (f_exact(x) - dpkj(x)) * T.jacobianDeterminant());
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ error += std::abs(qf.weight(i_q) * (f_exact(x) - dpkj(x)) * T.jacobianDeterminant(x_hat));
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]]};
+ QuadratureFormula<3> qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ error += std::abs(qf.weight(i_q) * (f_exact(x) - dpkj(x)) * T.jacobianDeterminant(x_hat));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T{xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
+ xr[cell_node_list[6]], xr[cell_node_list[7]]};
+ QuadratureFormula<3> qf =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{3});
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ const R3 x_hat = qf.point(i_q);
+ const R3 x = T(x_hat);
+ error += std::abs(qf.weight(i_q) * (f_exact(x) - dpkj(x)) * T.jacobianDeterminant(x_hat));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid cell type");
+ }
+ }
+
+ L1_error += error;
+ }
+ REQUIRE(parallel::allReduceMax(L1_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ // 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_affine = [](const R3& x) -> R3 {
+ // return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
+ // +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_affine(xj[cell_id]); });
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+
+ // {
+ // double max_mean_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // max_mean_error =
+ // std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ // }
+
+ // {
+ // double max_x_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+ // dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+ // max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // {
+ // double max_y_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+ // dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+ // max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // {
+ // double max_z_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+ // dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+ // max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_z_slope_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_affine = [](const R3& x) -> R2x2 {
+ // return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+ // //
+ // +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_affine(xj[cell_id]); });
+
+ // descriptor.setRowWeighting(false);
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ // auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
+
+ // {
+ // double max_mean_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // max_mean_error =
+ // std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ // }
+
+ // {
+ // double max_x_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+ // dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+ // max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1,
+ // //
+ // -2.3,
+ // +3.1}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ // }
+
+ // {
+ // double max_y_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+ // dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+ // max_y_slope_error =
+ // std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
+ // 1.3, +0.8}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // {
+ // double max_z_slope_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+ // dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+ // max_z_slope_error =
+ // std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
+ // +1.4, -1.8}));
+ // }
+ // REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+
+ // SECTION("vector data")
+ // {
+ // using R4 = TinyVector<4>;
+
+ // 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 vector_affine = [](const R3& x) -> R4 {
+ // return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+ // //
+ // +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();
+
+ // DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+ // parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ // auto vector = vector_affine(xj[cell_id]);
+ // for (size_t i = 0; i < vector.dimension(); ++i) {
+ // Vh[cell_id][i] = vector[i];
+ // }
+ // });
+
+ // descriptor.setPreconditioning(false);
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+
+ // {
+ // double max_mean_error = 0;
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // auto vector = vector_affine(xj[cell_id]);
+ // for (size_t i = 0; i < vector.dimension(); ++i) {
+ // max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ // }
+
+ // {
+ // double max_x_slope_error = 0;
+ // const R4 slope{+1.7, 2.1, -2.3, +3.1};
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // for (size_t i = 0; i < slope.dimension(); ++i) {
+ // const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
+ // dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
+
+ // max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ // }
+
+ // {
+ // double max_y_slope_error = 0;
+ // const R4 slope{+1.2, -2.2, 1.3, +0.8};
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // for (size_t i = 0; i < slope.dimension(); ++i) {
+ // const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
+ // dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
+
+ // max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // {
+ // double max_z_slope_error = 0;
+ // const R4 slope{-1.3, -2.4, +1.4, -1.8};
+ // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ // for (size_t i = 0; i < slope.dimension(); ++i) {
+ // const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
+ // dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
+
+ // max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+ // }
+ // }
+ // REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
+ // }
+ // }
+ // }
+ // }
+ }
+ }
+}