diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index bf1c51121d943b37e854e295b5c472140239fc78..3753143d38ff1ea8c2c6b32472039e30cb985109 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -7,10 +7,36 @@ void
test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
uint64_t degree)
{
- PolynomialReconstructionDescriptor descriptor{PolynomialReconstructionDescriptor::IntegrationMethod::element, degree};
- {
- const size_t nb_loops = 1;
- std::cout << "** variable list at once (" << nb_loops << " times)\n";
+ std::vector descriptor_list =
+ {PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::cell_center, degree},
+ PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::element, degree},
+ PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::boundary, degree}};
+
+ [[maybe_unused]] auto x =
+ PolynomialReconstruction{
+ PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::cell_center, degree}}
+ .build(discrete_function_variant_list);
+
+ for (auto&& descriptor : descriptor_list) {
+ std::string method_name;
+ switch (descriptor.integrationMethod()) {
+ case PolynomialReconstructionDescriptor::IntegrationMethod::element: {
+ method_name = "element";
+ break;
+ }
+ case PolynomialReconstructionDescriptor::IntegrationMethod::boundary: {
+ method_name = "boundary";
+ break;
+ }
+ case PolynomialReconstructionDescriptor::IntegrationMethod::cell_center: {
+ method_name = "cell_center";
+ break;
+ }
+ }
+
+ const size_t nb_loops = 50;
+ std::cout << "** variable list at once (" << nb_loops << " times) using " << rang::fgB::yellow << method_name
+ << rang::fg::reset << "\n";
PolynomialReconstruction reconstruction{descriptor};
Timer t;
diff --git a/tests/test_QuadraticPolynomialReconstruction.cpp b/tests/test_QuadraticPolynomialReconstruction.cpp
index 69c685f544f774205effbcf304dca43ff7d84658..6b992e960c3a47d0ddf050b30ac893492edd509b 100644
--- a/tests/test_QuadraticPolynomialReconstruction.cpp
+++ b/tests/test_QuadraticPolynomialReconstruction.cpp
@@ -584,16 +584,326 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
return L1_error;
};
+ std::map<std::string, PolynomialReconstructionDescriptor> name_method_map =
+ {{"element",
+ PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::element, 2}},
+ {"boundary",
+ PolynomialReconstructionDescriptor{PolynomialReconstructionDescriptor::IntegrationMethod::boundary, 2}}};
+
+ for (auto [method_name, method_descriptor] : name_method_map) {
+ SECTION(method_name)
+ {
+ 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{method_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{method_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{method_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("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ 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::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) * T.jacobianDeterminant() * exact(x);
+ }
+ 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) * T.jacobianDeterminant(x_hat) * exact(x);
+ }
+ 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) * T.jacobianDeterminant(x_hat) * exact(x);
+ }
+ 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) * 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::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);
+ DataType diff = qf.weight(i_q) * T.jacobianDeterminant() * (exact(x) - dpkj(x));
+ error += sum_abs(diff);
+ }
+ 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);
+ DataType diff = qf.weight(i_q) * T.jacobianDeterminant(x_hat) * (exact(x) - dpkj(x));
+ error += sum_abs(diff);
+ }
+ 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);
+ DataType diff = qf.weight(i_q) * T.jacobianDeterminant(x_hat) * (exact(x) - dpkj(x));
+ error += sum_abs(diff);
+ }
+ 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);
+ 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()) {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
- auto& mesh = *p_mesh;
+ 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];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double xz = x * z;
+ const double y2 = y * y;
+ const double yz = y * z;
+ const double z2 = z * z;
- 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];
+ return 2.3 + 1.7 * x + 0.2 * y + 1.4 * z - 3.2 * x2 + 1.3 * xy - 1.6 * xz - 1.4 * y2 + 2 * yz - 1.8 * z2;
};
auto xr = mesh.xr();
@@ -614,7 +924,7 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, 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));
@@ -624,24 +934,29 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
SECTION("R^3 data")
{
- using R3 = TinyVector<3>;
-
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh;
- auto u_exact = [](const R2& X) -> R3 {
+ auto u_exact = [](const R3& X) -> R3 {
const double x = X[0];
const double y = X[1];
+ const double z = X[2];
const double x2 = x * x;
const double xy = x * y;
+ const double xz = x * z;
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};
+ const double yz = y * z;
+ const double z2 = z * z;
+
+ return R3{+2.3 + 1.7 * x - 2.2 * y + 1.8 * z //
+ + 0.3 * x2 - 1.3 * xy + 2.7 * xz + 0.7 * y2 + 2.0 * yz - 1.2 * z2, //
+ +1.4 - 0.6 * x + 1.3 * y - 3.7 * z //
+ + 3.1 * x2 - 2.4 * xy - 1.5 * xz - 1.9 * y2 + 0.2 * yz + 0.9 * z2, //
+ -0.2 + 3.1 * x - 1.1 * y + 1.9 * z //
+ - 1.3 * x2 + 2.8 * xy - 2.1 * xz + 3.2 * y2 - 2.1 * yz + 1.6 * z2};
};
auto xr = mesh.xr();
@@ -662,8 +977,7 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
- auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
-
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, 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));
}
@@ -674,23 +988,32 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
{
using R2x2 = TinyMatrix<2, 2>;
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh;
- auto A_exact = [](const R2& X) -> R2x2 {
+ auto A_exact = [](const R3& X) -> R2x2 {
const double x = X[0];
const double y = X[1];
+ const double z = X[2];
const double x2 = x * x;
const double xy = x * y;
+ const double xz = x * z;
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};
+ const double yz = y * z;
+ const double z2 = z * z;
+
+ return R2x2{+2.3 + 1.7 * x - 2.2 * y + 1.8 * z //
+ + 0.3 * x2 - 1.3 * xy + 2.7 * xz + 0.7 * y2 + 2.0 * yz - 1.2 * z2, //
+ +1.4 - 0.6 * x + 1.3 * y - 3.7 * z //
+ + 3.1 * x2 - 2.4 * xy - 1.5 * xz - 1.9 * y2 + 0.2 * yz + 0.9 * z2, //
+ //
+ -0.2 + 3.1 * x - 1.1 * y + 1.9 * z //
+ - 1.3 * x2 + 2.8 * xy - 2.1 * xz + 3.2 * y2 - 2.1 * yz + 1.6 * z2, //
+ +0.9 - 1.3 * x + 2.1 * y + 3.1 * z //
+ + 1.1 * x2 + 1.8 * xy + 1.2 * xz - 2.3 * y2 + 3.3 * yz - 1.2 * z2};
};
auto xr = mesh.xr();
@@ -709,490 +1032,13 @@ TEST_CASE("QuadraticPolynomialReconstruction", "[scheme]")
Ah[cell_id] = 1. / Vj[cell_id] * value;
});
+ descriptor.setRowWeighting(false);
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));
- }
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
}
}
}
-
- // 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));
- // }
- // }
- // }
- // }
}
}
}