diff --git a/src/analysis/CubeGaussQuadrature.cpp b/src/analysis/CubeGaussQuadrature.cpp index f1aaf9fd23dfdd55cd5ac3eced298a2b83101204..8f49fccbeb62b5e4e06da65fbccc0050bb0a0ce8 100644 --- a/src/analysis/CubeGaussQuadrature.cpp +++ b/src/analysis/CubeGaussQuadrature.cpp @@ -492,10 +492,11 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } - + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(CubeGaussQuadrature::max_degree) + " on cubes"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(CubeGaussQuadrature::max_degree) + " on cubes"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/PrismGaussQuadrature.cpp b/src/analysis/PrismGaussQuadrature.cpp index 10c8fab3720ab794299522452e5528d25f38e394..623d3f87864b08fc7b65586ebf09453568dabcf1 100644 --- a/src/analysis/PrismGaussQuadrature.cpp +++ b/src/analysis/PrismGaussQuadrature.cpp @@ -796,9 +796,11 @@ PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(PrismGaussQuadrature::max_degree) + "on prisms"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(PrismGaussQuadrature::max_degree) + "on prisms"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/PyramidGaussQuadrature.cpp b/src/analysis/PyramidGaussQuadrature.cpp index 4b4d17359486f0f44bff24dfed013c790f50a36e..48eb4ec3421f1f6498420c9c851994520dc1b486 100644 --- a/src/analysis/PyramidGaussQuadrature.cpp +++ b/src/analysis/PyramidGaussQuadrature.cpp @@ -968,9 +968,11 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(PyramidGaussQuadrature::max_degree) + " on pyramids"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(PyramidGaussQuadrature::max_degree) + " on pyramids"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/QuadratureManager.hpp b/src/analysis/QuadratureManager.hpp index b0f81c543c0612181cc5728b2c7ef36fc1e25515..04655f291f8ff75e8c1207c4b6171a7b7f461a05 100644 --- a/src/analysis/QuadratureManager.hpp +++ b/src/analysis/QuadratureManager.hpp @@ -47,6 +47,10 @@ class QuadratureManager const QuadratureFormula<1>& getLineFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxLineDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxLineDegree(quadrature_descriptor.type())) + " on lines"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: case QuadratureType::GaussLegendre: { @@ -56,28 +60,40 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_line_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2]; } + // LCOV_EXCL_START default: { - throw UnexpectedError("invalid quadrature type"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on lines"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<2>& getTriangleFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxTriangleDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxTriangleDegree(quadrature_descriptor.type())) + " on triangles"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_triangle_gauss_formula_list[quadrature_descriptor.degree() - 1]; } + // LCOV_EXCL_START default: { - throw UnexpectedError(quadrature_descriptor.name() + " is not defined on triangles"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on triangles"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<2>& getSquareFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxSquareDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxSquareDegree(quadrature_descriptor.type())) + " on squares"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_square_gauss_formula_list[quadrature_descriptor.degree() / 2]; @@ -88,54 +104,78 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_square_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2]; } + // LCOV_EXCL_START default: { - throw UnexpectedError("invalid quadrature type"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on squares"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<3>& getTetrahedronFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxTetrahedronDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxTetrahedronDegree(quadrature_descriptor.type())) + " on tetrahedra"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_tetrahedron_gauss_formula_list[quadrature_descriptor.degree() - 1]; } + // LCOV_EXCL_START default: { - throw UnexpectedError(quadrature_descriptor.name() + " is not defined on tetrahedron"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on tetrahedra"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<3>& getPrismFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxPrismDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxPrismDegree(quadrature_descriptor.type())) + " on prisms"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_prism_gauss_formula_list[quadrature_descriptor.degree() - 1]; } + // LCOV_EXCL_START default: { - throw UnexpectedError(quadrature_descriptor.name() + " is not defined on prism"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on prisms"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<3>& getPyramidFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxPyramidDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxPyramidDegree(quadrature_descriptor.type())) + " on pyramids"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_pyramid_gauss_formula_list[quadrature_descriptor.degree() - 1]; } + // LCOV_EXCL_START default: { - throw UnexpectedError(quadrature_descriptor.name() + " is not defined on pyramid"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on pyramid"); } + // LCOV_EXCL_STOP } } const QuadratureFormula<3>& getCubeFormula(const IQuadratureDescriptor& quadrature_descriptor) const { + if (quadrature_descriptor.degree() > maxCubeDegree(quadrature_descriptor.type())) { + throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " + + std::to_string(maxCubeDegree(quadrature_descriptor.type())) + " on cubes"); + } switch (quadrature_descriptor.type()) { case QuadratureType::Gauss: { return m_cube_gauss_formula_list[quadrature_descriptor.degree() / 2]; @@ -147,9 +187,11 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_cube_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2]; } + // LCOV_EXCL_START default: { - throw UnexpectedError("invalid quadrature type"); + throw UnexpectedError(::name(quadrature_descriptor.type()) + " is not defined on cubes"); } + // LCOV_EXCL_STOP } } @@ -164,9 +206,11 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_line_gauss_lobatto_formula_list.size() * 2 - 1; } + // LCOV_EXCL_START default: { throw UnexpectedError("invalid quadrature type"); } + // LCOV_EXCL_STOP } } @@ -183,9 +227,11 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_square_gauss_lobatto_formula_list.size() * 2 - 1; } + // LCOV_EXCL_START default: { throw UnexpectedError("invalid quadrature type"); } + // LCOV_EXCL_STOP } } @@ -196,9 +242,11 @@ class QuadratureManager case QuadratureType::Gauss: { return m_triangle_gauss_formula_list.size(); } + // LCOV_EXCL_START default: { - throw UnexpectedError(name(type) + " is not defined on triangle"); + throw UnexpectedError(::name(type) + " is not defined on triangle"); } + // LCOV_EXCL_STOP } } @@ -215,9 +263,11 @@ class QuadratureManager case QuadratureType::GaussLobatto: { return m_cube_gauss_lobatto_formula_list.size() * 2 - 1; } + // LCOV_EXCL_START default: { throw UnexpectedError("invalid quadrature type"); } + // LCOV_EXCL_STOP } } @@ -228,9 +278,11 @@ class QuadratureManager case QuadratureType::Gauss: { return m_prism_gauss_formula_list.size(); } + // LCOV_EXCL_START default: { throw UnexpectedError(::name(type) + " is not defined on prism"); } + // LCOV_EXCL_STOP } } @@ -241,9 +293,11 @@ class QuadratureManager case QuadratureType::Gauss: { return m_pyramid_gauss_formula_list.size(); } + // LCOV_EXCL_START default: { throw UnexpectedError(::name(type) + " is not defined on pyramid"); } + // LCOV_EXCL_STOP } } @@ -254,9 +308,11 @@ class QuadratureManager case QuadratureType::Gauss: { return m_tetrahedron_gauss_formula_list.size(); } + // LCOV_EXCL_START default: { throw UnexpectedError(::name(type) + " is not defined on tetrahedron"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/QuadratureType.hpp b/src/analysis/QuadratureType.hpp index 1c67daa0daec9f4692e909194a5f2861245f0c7a..232a5fda0376852e9e9667ba1bd6bfd6c12287ac 100644 --- a/src/analysis/QuadratureType.hpp +++ b/src/analysis/QuadratureType.hpp @@ -27,9 +27,11 @@ name(QuadratureType type) case QuadratureType::GaussLobatto: { return "Gauss-Lobatto"; } + // LCOV_EXCL_START default: { throw UnexpectedError("unknown quadrature type name"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/SquareGaussQuadrature.cpp b/src/analysis/SquareGaussQuadrature.cpp index 20e27c8c02d9ebce1c2305834aeda3b6430a0d27..2597e013b1b4f7ede3b69981fb43bfd4e901aff7 100644 --- a/src/analysis/SquareGaussQuadrature.cpp +++ b/src/analysis/SquareGaussQuadrature.cpp @@ -319,9 +319,11 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(SquareGaussQuadrature::max_degree) + " on squares"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(SquareGaussQuadrature::max_degree) + " on squares"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/TensorialGaussLegendreQuadrature.cpp b/src/analysis/TensorialGaussLegendreQuadrature.cpp index 634ac9f9559ee87a6926dccf28d6c966aede0acc..a306f44135b7a4cf31ddc6e906f0352ca1bde067 100644 --- a/src/analysis/TensorialGaussLegendreQuadrature.cpp +++ b/src/analysis/TensorialGaussLegendreQuadrature.cpp @@ -229,10 +229,12 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr weight_list[11] = 0.0471753363865118271946160; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss-Legendre quadrature formulae handle degrees up to " + - std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree)); + throw UnexpectedError("Gauss-Legendre quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree)); } + // LCOV_EXCL_STOP } m_point_list = point_list; diff --git a/src/analysis/TensorialGaussLobattoQuadrature.cpp b/src/analysis/TensorialGaussLobattoQuadrature.cpp index d848943751a6c5e3487b5fcc24394c66985a271a..f53b49434335643684be80e20d5d750dde1c4e93 100644 --- a/src/analysis/TensorialGaussLobattoQuadrature.cpp +++ b/src/analysis/TensorialGaussLobattoQuadrature.cpp @@ -118,10 +118,12 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre weight_list[7] = 0.0357142857142857142857143; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss-Lobatto quadrature formulae handle degrees up to " + - std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree)); + throw UnexpectedError("Gauss-Lobatto quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree)); } + // LCOV_EXCL_STOP } m_point_list = point_list; diff --git a/src/analysis/TetrahedronGaussQuadrature.cpp b/src/analysis/TetrahedronGaussQuadrature.cpp index e5ee82d3fb149ab7c50647d98092754f146dcab1..49fd5efe685c765aedccae3ee1b2c59d82915dd3 100644 --- a/src/analysis/TetrahedronGaussQuadrature.cpp +++ b/src/analysis/TetrahedronGaussQuadrature.cpp @@ -680,9 +680,11 @@ TetrahedronGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(TetrahedronGaussQuadrature::max_degree) + " on tetrahedra"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(TetrahedronGaussQuadrature::max_degree) + " on tetrahedra"); } + // LCOV_EXCL_STOP } } diff --git a/src/analysis/TriangleGaussQuadrature.cpp b/src/analysis/TriangleGaussQuadrature.cpp index 19013da54f79268329416b5499fd80748612b4cc..627edc4ff8293a7a00b689b5602fa2f8d20c3911 100644 --- a/src/analysis/TriangleGaussQuadrature.cpp +++ b/src/analysis/TriangleGaussQuadrature.cpp @@ -482,9 +482,11 @@ TriangleGaussQuadrature::_buildPointAndWeightLists(const size_t degree) m_weight_list = weight_list; break; } + // LCOV_EXCL_START default: { - throw NormalError("Gauss quadrature formulae handle degrees up to " + - std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles"); + throw UnexpectedError("Gauss quadrature formulae handle degrees up to " + + std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles"); } + // LCOV_EXCL_STOP } } diff --git a/src/scheme/FaceIntegrator.hpp b/src/scheme/FaceIntegrator.hpp index 30ee39077ea12d898b4bdee00f12970be79e8753..657ca7a88a3a959f2856bc34ae8938ab91561785 100644 --- a/src/scheme/FaceIntegrator.hpp +++ b/src/scheme/FaceIntegrator.hpp @@ -144,9 +144,11 @@ class FaceIntegrator } break; } + // LCOV_EXCL_START default: { invalid_face = std::make_pair(true, face_id); } + // LCOV_EXCL_STOP } }); @@ -159,7 +161,9 @@ class FaceIntegrator // LCOV_EXCL_STOP } else { + // LCOV_EXCL_START throw UnexpectedError("invalid dimension for face integration"); + // LCOV_EXCL_STOP } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6bccbfae744ed002e0edcde321205f309e88105a..18d5656925cae0cb5d0cba949bebb2ce74ee4092 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -54,10 +54,12 @@ add_executable (unit_tests test_BuiltinFunctionProcessor.cpp test_CastArray.cpp test_CellIntegrator.cpp + test_CellType.cpp test_ConsoleManager.cpp test_CG.cpp - test_ContinueProcessor.cpp test_ConcatExpressionProcessor.cpp + test_ConsoleManager.cpp + test_ContinueProcessor.cpp test_CRSGraph.cpp test_CRSMatrix.cpp test_CRSMatrixDescriptor.cpp @@ -91,9 +93,13 @@ add_executable (unit_tests test_FunctionProcessor.cpp test_FunctionSymbolId.cpp test_FunctionTable.cpp + test_GaussLegendreQuadratureDescriptor.cpp + test_GaussLobattoQuadratureDescriptor.cpp + test_GaussQuadratureDescriptor.cpp test_IfProcessor.cpp test_IncDecExpressionProcessor.cpp test_INodeProcessor.cpp + test_ItemId.cpp test_ItemType.cpp test_LinearSolver.cpp test_LinearSolverOptions.cpp @@ -117,6 +123,9 @@ add_executable (unit_tests test_PugsUtils.cpp test_PyramidGaussQuadrature.cpp test_PyramidTransformation.cpp + test_QuadratureType.cpp + test_RefId.cpp + test_RefItemList.cpp test_RevisionInfo.cpp test_SmallArray.cpp test_SmallVector.cpp diff --git a/tests/test_CellIntegrator.cpp b/tests/test_CellIntegrator.cpp index 49f750afb6ebd7a28135949b2d6aca2473b24c44..ebd322d14c8ea5df4053163dcae9b60d02318a53 100644 --- a/tests/test_CellIntegrator.cpp +++ b/tests/test_CellIntegrator.cpp @@ -15,915 +15,2766 @@ TEST_CASE("CellIntegrator", "[scheme]") { - SECTION("1D") + SECTION("scalar") { - using R1 = TinyVector<1>; - - const auto mesh = MeshDataBaseForTests::get().unordered1DMesh(); - auto f = [](const R1& x) -> double { return x[0] * x[0] + 1; }; - - Array<const double> int_f_per_cell = [=] { - Array<double> int_f(mesh->numberOfCells()); - auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { - auto cell_node_list = cell_to_node_matrix[cell_id]; - auto xr = mesh->xr(); - const double x_left = xr[cell_node_list[0]][0]; - const double x_right = xr[cell_node_list[1]][0]; - int_f[cell_id] = 1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left; - }); - - return int_f; - }(); - - SECTION("direct formula") + SECTION("1D") { - SECTION("all cells") + using R1 = TinyVector<1>; + + const auto mesh = MeshDataBaseForTests::get().unordered1DMesh(); + auto f = [](const R1& x) -> double { return x[0] * x[0] + 1; }; + + Array<const double> int_f_per_cell = [=] { + Array<double> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + const double x_left = xr[cell_node_list[0]][0]; + const double x_right = xr[cell_node_list[1]][0]; + int_f[cell_id] = 1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left; + }); + + return int_f; + }(); + + SECTION("direct formula") { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + + SECTION("tensorial formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + + SECTION("2D") + { + using R2 = TinyVector<2>; + + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + auto f = [](const R2& X) -> double { + const double x = X[0]; + const double y = X[1]; + return x * x + 2 * x * y + 3 * y * y + 2; + }; + + Array<const double> int_f_per_cell = [=] { + Array<double> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + double integral = 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]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + default: { + throw UnexpectedError("invalid cell type in 2d"); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("direct formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + + SECTION("tensorial formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> double { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + SECTION(mesh_name) + { + SECTION("direct formula") + { + Array<const double> int_f_per_cell = [=] { + Array<double> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + double integral = 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]]); + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + { // top tetrahedron + TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[3]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi)); + } + } + { // bottom tetrahedron + TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi)); + } + } + } else if (cell_node_list.size() == 6) { + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + { // top pyramid + PyramidTransformation T0(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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom pyramid + PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); + } + break; + } + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == 0); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = + CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == 0); + } + } + } + + SECTION("tensorial formula") + { + Array<const double> int_f_per_cell = [=] { + Array<double> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4)); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + double integral = 0; + + switch (cell_type[cell_id]) { + case CellType::Tetrahedron: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[3]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Pyramid: { + 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[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Prism: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + { // top tetrahedron + CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom tetrahedron + CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else if (cell_node_list.size() == 6) { + { // top pyramid + CubeTransformation T0(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[5]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom pyramid + CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); + } + break; + } + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<double> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + auto cell_type = mesh->connectivity().cellType(); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<double> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<double> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + } + } + } + + SECTION("vector") + { + using R2 = TinyVector<2>; + + SECTION("1D") + { + using R1 = TinyVector<1>; + + const auto mesh = MeshDataBaseForTests::get().unordered1DMesh(); + auto f = [](const R1& x) -> R2 { return R2{x[0] * x[0] + 1, 2 * x[0]}; }; + + Array<const R2> int_f_per_cell = [=] { + Array<R2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + const double x_left = xr[cell_node_list[0]][0]; + const double x_right = xr[cell_node_list[1]][0]; + int_f[cell_id] = R2{1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left, + x_right * x_right - x_left * x_left}; + }); + + return int_f; + }(); + + SECTION("direct formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + + SECTION("tensorial formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + + SECTION("2D") + { + using R2 = TinyVector<2>; + + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + auto f = [](const R2& X) -> R2 { + const double x = X[0]; + const double y = X[1]; + return R2{x * x + 2 * x * y + 3 * y * y + 2, 2 * x + 3 * y * y}; + }; + + Array<const R2> int_f_per_cell = [=] { + Array<R2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2 integral = zero; + + 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]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + default: { + throw UnexpectedError("invalid cell type in 2d"); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("direct formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + + SECTION("tensorial formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> R2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x * y - 3 * y * z + 2 * x}; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + SECTION(mesh_name) + { + SECTION("direct formula") + { + Array<const R2> int_f_per_cell = [=] { + Array<R2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2 integral = zero; + + 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]]); + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + { // top tetrahedron + TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[3]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi)); + } + } + { // bottom tetrahedron + TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi)); + } + } + } else if (cell_node_list.size() == 6) { + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + { // top pyramid + PyramidTransformation T0(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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom pyramid + PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); + } + break; + } + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == 0); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == 0); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == 0); + } + } + } + + SECTION("tensorial formula") + { + Array<const R2> int_f_per_cell = [=] { + Array<R2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4)); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2 integral = zero; + + switch (cell_type[cell_id]) { + case CellType::Tetrahedron: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[3]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Pyramid: { + 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[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Prism: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; + } + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + { // top tetrahedron + CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom tetrahedron + CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else if (cell_node_list.size() == 6) { + { // top pyramid + CubeTransformation T0(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[5]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom pyramid + CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); + } + break; + } + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); + + return int_f; + }(); + + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + auto cell_type = mesh->connectivity().cellType(); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfCells()); + + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("cell list") + { + SECTION("Array") + { + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + Array<R2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); + } + + SmallArray<R2> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } + } + } + } + } + + SECTION("matrix") + { + using R2x2 = TinyMatrix<2>; + auto l2Norm = [](const R2x2& A) -> double { + return std::sqrt(A(0, 0) * A(0, 0) + A(1, 0) * A(1, 0) + A(0, 1) * A(0, 1) + A(1, 1) * A(1, 1)); + }; + + SECTION("1D") + { + using R1 = TinyVector<1>; + + const auto mesh = MeshDataBaseForTests::get().unordered1DMesh(); + auto f = [](const R1& x) -> R2x2 { return R2x2{x[0] * x[0] + 1, 2 * x[0], 3 - x[0], 3 * x[0] - 1}; }; + + Array<const R2x2> int_f_per_cell = [=] { + Array<R2x2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + const double x_left = xr[cell_node_list[0]][0]; + const double x_right = xr[cell_node_list[1]][0]; + int_f[cell_id] = R2x2{1. / 3 * (x_right * x_right * x_right - x_left * x_left * x_left) + x_right - x_left, + x_right * x_right - x_left * x_left, + -0.5 * (x_right * x_right - x_left * x_left) + 3 * (x_right - x_left), + 3. / 2 * (x_right * x_right - x_left * x_left) - (x_right - x_left)}; + }); + + return int_f; + }(); + + SECTION("direct formula") + { + SECTION("all cells") + { + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); - Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); - SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } - } - SECTION("tensorial formula") - { - SECTION("all cells") + SECTION("tensorial formula") { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, - values); + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R1 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); - Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); - SmallArray<double> values = - CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } } - } - SECTION("2D") - { - using R2 = TinyVector<2>; - - const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); - - auto f = [](const R2& X) -> double { - const double x = X[0]; - const double y = X[1]; - return x * x + 2 * x * y + 3 * y * y + 2; - }; - - Array<const double> int_f_per_cell = [=] { - Array<double> int_f(mesh->numberOfCells()); - auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); - auto cell_type = mesh->connectivity().cellType(); - - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { - auto cell_node_list = cell_to_node_matrix[cell_id]; - auto xr = mesh->xr(); - double integral = 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]]); - auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4)); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + SECTION("2D") + { + using R2 = TinyVector<2>; + + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + auto f = [](const R2& X) -> R2x2 { + const double x = X[0]; + const double y = X[1]; + return R2x2{x * x + 2 * x * y + 3 * y * y + 2, 2 * x + 3 * y * y, 2 * x - 3 * y, 3 * x * x - 2 * y}; + }; + + Array<const R2x2> int_f_per_cell = [=] { + Array<R2x2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2x2 integral = zero; + + 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]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + break; } - 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]]); - auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(4)); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + 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]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; } - break; - } - default: { - throw UnexpectedError("invalid cell type in 2d"); - } - } - int_f[cell_id] = integral; - }); + default: { + throw UnexpectedError("invalid cell type in 2d"); + } + } + int_f[cell_id] = integral; + }); - return int_f; - }(); + return int_f; + }(); - SECTION("direct formula") - { - SECTION("all cells") + SECTION("direct formula") { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); - Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); - SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } - } - SECTION("tensorial formula") - { - SECTION("all cells") + SECTION("tensorial formula") { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, - values); + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); - Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); - SmallArray<double> values = - CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } } - } - SECTION("3D") - { - using R3 = TinyVector<3>; + SECTION("3D") + { + using R3 = TinyVector<3>; - auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); - auto f = [](const R3& X) -> double { - const double x = X[0]; - const double y = X[1]; - const double z = X[2]; - return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; - }; + auto f = [](const R3& X) -> R2x2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2x2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x * y - 3 * y * z + 2 * x, + x * x - 2 * x * z - 3 * y, 3 * z + x * y}; + }; - std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; - mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); - mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); - for (auto mesh_info : mesh_list) { - auto mesh_name = mesh_info.first; - auto mesh = mesh_info.second; + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; - SECTION(mesh_name) - { - SECTION("direct formula") + SECTION(mesh_name) { - Array<const double> int_f_per_cell = [=] { - Array<double> int_f(mesh->numberOfCells()); - auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); - auto cell_type = mesh->connectivity().cellType(); - - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { - auto cell_node_list = cell_to_node_matrix[cell_id]; - auto xr = mesh->xr(); - double integral = 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]]); - auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + SECTION("direct formula") + { + Array<const R2x2> int_f_per_cell = [=] { + Array<R2x2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2x2 integral = zero; + + 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]]); + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi)); + } + break; } - 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]]); - auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + 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]]); + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; } - 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]]); - auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4)); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + 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]]); + auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; } - 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]]); - auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4)); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + 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]]); + auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + break; } - break; - } - case CellType::Diamond: { - if (cell_node_list.size() == 5) { - auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); - { // top tetrahedron - TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]], - xr[cell_node_list[4]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi)); + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4)); + { // top tetrahedron + TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[3]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi)); + } } - } - { // bottom tetrahedron - TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]], - xr[cell_node_list[0]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi)); + { // bottom tetrahedron + TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi)); + } } - } - } else if (cell_node_list.size() == 6) { - auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); - { // top pyramid - PyramidTransformation T0(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]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } else if (cell_node_list.size() == 6) { + auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4)); + { // top pyramid + PyramidTransformation T0(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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } } - } - { // bottom pyramid - PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], - xr[cell_node_list[1]], xr[cell_node_list[0]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + { // bottom pyramid + PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); } - } else { - INFO("Diamond cells with more than 6 vertices are not tested"); - REQUIRE(false); + break; } - break; - } - default: { - INFO("Diamond cells not tested yet"); - REQUIRE(cell_type[cell_id] != CellType::Diamond); - } - } - int_f[cell_id] = integral; - }); + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); - return int_f; - }(); + return int_f; + }(); - SECTION("all cells") - { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, - values); + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == 0); } - REQUIRE(error == 0); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == 0); } - REQUIRE(error == 0); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == 0); } - - REQUIRE(error == 0); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); - Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == 0); } - REQUIRE(error == 0); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); - SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == 0); } - - REQUIRE(error == 0); } } - } - SECTION("tensorial formula") - { - Array<const double> int_f_per_cell = [=] { - Array<double> int_f(mesh->numberOfCells()); - auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); - auto cell_type = mesh->connectivity().cellType(); - - auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4)); - - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { - auto cell_node_list = cell_to_node_matrix[cell_id]; - auto xr = mesh->xr(); - double integral = 0; - - switch (cell_type[cell_id]) { - case CellType::Tetrahedron: { - CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], - xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]], - xr[cell_node_list[3]], xr[cell_node_list[3]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); - } - break; - } - case CellType::Pyramid: { - 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[4]], - xr[cell_node_list[4]], xr[cell_node_list[4]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); - } - break; - } - case CellType::Prism: { - CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], - xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]], - xr[cell_node_list[5]], xr[cell_node_list[5]]); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); - } - 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]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); - } - break; - } - case CellType::Diamond: { - if (cell_node_list.size() == 5) { - { // top tetrahedron - CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]], - xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]], - xr[cell_node_list[4]], xr[cell_node_list[4]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); - } + SECTION("tensorial formula") + { + Array<const R2x2> int_f_per_cell = [=] { + Array<R2x2> int_f(mesh->numberOfCells()); + auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + auto cell_type = mesh->connectivity().cellType(); + + auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4)); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto cell_node_list = cell_to_node_matrix[cell_id]; + auto xr = mesh->xr(); + R2x2 integral = zero; + + switch (cell_type[cell_id]) { + case CellType::Tetrahedron: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[3]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); } - { // bottom tetrahedron - CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]], - xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], - xr[cell_node_list[0]], xr[cell_node_list[0]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); - } + break; + } + case CellType::Pyramid: { + 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[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); } - } else if (cell_node_list.size() == 6) { - { // top pyramid - CubeTransformation T0(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[5]], - xr[cell_node_list[5]], xr[cell_node_list[5]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); - } + break; + } + case CellType::Prism: { + CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]], + xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); + } + 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]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi)); } - { // bottom pyramid - CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], - xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], - xr[cell_node_list[0]], xr[cell_node_list[0]]); - - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + break; + } + case CellType::Diamond: { + if (cell_node_list.size() == 5) { + { // top tetrahedron + CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]], + xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]], + xr[cell_node_list[4]], xr[cell_node_list[4]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } } + { // bottom tetrahedron + CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else if (cell_node_list.size() == 6) { + { // top pyramid + CubeTransformation T0(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[5]], + xr[cell_node_list[5]], xr[cell_node_list[5]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi)); + } + } + { // bottom pyramid + CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]], + xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]], + xr[cell_node_list[0]], xr[cell_node_list[0]]); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi)); + } + } + } else { + INFO("Diamond cells with more than 6 vertices are not tested"); + REQUIRE(false); } - } else { - INFO("Diamond cells with more than 6 vertices are not tested"); - REQUIRE(false); + break; } - break; - } - default: { - INFO("Diamond cells not tested yet"); - REQUIRE(cell_type[cell_id] != CellType::Diamond); - } - } - int_f[cell_id] = integral; - }); + default: { + INFO("Diamond cells not tested yet"); + REQUIRE(cell_type[cell_id] != CellType::Diamond); + } + } + int_f[cell_id] = integral; + }); - return int_f; - }(); + return int_f; + }(); - SECTION("all cells") - { - SECTION("CellValue") + SECTION("all cells") { - CellValue<double> values(mesh->connectivity()); - CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), - *mesh, values); + SECTION("CellValue") + { + CellValue<R2x2> values(mesh->connectivity()); + CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + auto cell_type = mesh->connectivity().cellType(); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - auto cell_type = mesh->connectivity().cellType(); - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfCells()); - SECTION("Array") - { - Array<double> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfCells()); + CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfCells()); - CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + error += l2Norm(int_f_per_cell[cell_id] - values[cell_id]); + } - double error = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - error += std::abs(int_f_per_cell[cell_id] - values[cell_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - } - SECTION("cell list") - { - SECTION("Array") + SECTION("cell list") { - Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - + SECTION("Array") { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } + + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + Array<R2x2> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); - Array<double> values = - CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("SmallArray") + { + SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; - SECTION("SmallArray") - { - SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2}; + { + size_t k = 0; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { + cell_list[k] = cell_id; + } - { - size_t k = 0; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) { - cell_list[k] = cell_id; + REQUIRE(k == cell_list.size()); } - REQUIRE(k == cell_list.size()); - } + SmallArray<R2x2> values = + CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); - SmallArray<double> values = - CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list); + double error = 0; + for (size_t i = 0; i < cell_list.size(); ++i) { + error += l2Norm(int_f_per_cell[cell_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < cell_list.size(); ++i) { - error += std::abs(int_f_per_cell[cell_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - - REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } } diff --git a/tests/test_CellType.cpp b/tests/test_CellType.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0e4a42331cf6606c9d8c0edd60f372676e6bb4c0 --- /dev/null +++ b/tests/test_CellType.cpp @@ -0,0 +1,25 @@ +#include <catch2/catch_test_macros.hpp> + +#include <mesh/CellType.hpp> + +#include <limits> +#include <type_traits> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("CellType", "[mesh]") +{ + REQUIRE(name(CellType::Line) == "line"); + + REQUIRE(name(CellType::Triangle) == "triangle"); + REQUIRE(name(CellType::Quadrangle) == "quadrangle"); + + REQUIRE(name(CellType::Diamond) == "diamond"); + REQUIRE(name(CellType::Hexahedron) == "hexahedron"); + REQUIRE(name(CellType::Prism) == "prism"); + REQUIRE(name(CellType::Pyramid) == "pyramid"); + REQUIRE(name(CellType::Tetrahedron) == "tetrahedron"); + + REQUIRE(name(static_cast<CellType>(std::numeric_limits<std::underlying_type_t<CellType>>::max())) == + "unknown cell type"); +} diff --git a/tests/test_CubeGaussQuadrature.cpp b/tests/test_CubeGaussQuadrature.cpp index 1db30eb178be24be223f63b5f6a68aa441261909..d4221fe322dc6dc4f898ac2563dc14d463b2e71f 100644 --- a/tests/test_CubeGaussQuadrature.cpp +++ b/tests/test_CubeGaussQuadrature.cpp @@ -488,6 +488,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]") SECTION("max implemented degree") { REQUIRE(QuadratureManager::instance().maxCubeDegree(QuadratureType::Gauss) == CubeGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula( + GaussQuadratureDescriptor(CubeGaussQuadrature ::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(CubeGaussQuadrature ::max_degree) + " on cubes"); } SECTION("Access functions") diff --git a/tests/test_EdgeIntegrator.cpp b/tests/test_EdgeIntegrator.cpp index f8ea09cee6bdeef7352281f207e2dfcc3b66e408..5dfccbc6c7427a2bd27634c49957209c61c6f3ce 100644 --- a/tests/test_EdgeIntegrator.cpp +++ b/tests/test_EdgeIntegrator.cpp @@ -1,6 +1,7 @@ #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> +#include <algebra/TinyMatrix.hpp> #include <analysis/GaussLegendreQuadratureDescriptor.hpp> #include <analysis/GaussLobattoQuadratureDescriptor.hpp> #include <analysis/GaussQuadratureDescriptor.hpp> @@ -15,242 +16,745 @@ TEST_CASE("EdgeIntegrator", "[scheme]") { - SECTION("3D") + SECTION("scalar") { - using R3 = TinyVector<3>; + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> double { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const double> int_f_per_edge = [=] { + Array<double> int_f(mesh->numberOfEdges()); + auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix(); + + parallel_for( + mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) { + auto edge_node_list = edge_to_node_matrix[edge_id]; + auto xr = mesh->xr(); + double integral = 0; + + LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); + } - auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + int_f[edge_id] = integral; + }); - auto f = [](const R3& X) -> double { - const double x = X[0]; - const double y = X[1]; - const double z = X[2]; - return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; - }; + return int_f; + }(); - std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; - mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); - mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); - - for (auto mesh_info : mesh_list) { - auto mesh_name = mesh_info.first; - auto mesh = mesh_info.second; - - Array<const double> int_f_per_edge = [=] { - Array<double> int_f(mesh->numberOfEdges()); - auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix(); - - parallel_for( - mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) { - auto edge_node_list = edge_to_node_matrix[edge_id]; - auto xr = mesh->xr(); - double integral = 0; - - LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]); - auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); - } + SECTION(mesh_name) + { + SECTION("direct formula") + { + SECTION("all edges") + { + SECTION("EdgeValue") + { + EdgeValue<double> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); - int_f[edge_id] = integral; - }); + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } - return int_f; - }(); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - SECTION(mesh_name) - { - SECTION("direct formula") - { - SECTION("all edges") - { - SECTION("EdgeValue") + SECTION("Array") + { + Array<double> values(mesh->numberOfEdges()); + + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("edge list") { - EdgeValue<double> values(mesh->connectivity()); - EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, - values); + SECTION("Array") + { + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + REQUIRE(k == edge_list.size()); + } + + Array<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<double> values = + EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } + } - SECTION("Array") + SECTION("tensorial formula") + { + SECTION("all edges") { - Array<double> values(mesh->numberOfEdges()); + SECTION("EdgeValue") + { + EdgeValue<double> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); - EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("Array") + { + Array<double> values(mesh->numberOfEdges()); + + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("edge list") { - SmallArray<double> values(mesh->numberOfEdges()); - EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + SECTION("Array") + { + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + REQUIRE(k == edge_list.size()); + } + + Array<double> values = + EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<double> values = + EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } + } + } + } + } + + SECTION("vector") + { + using R2 = TinyVector<2>; + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> R2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x + 3 * y - 1}; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const R2> int_f_per_edge = [=] { + Array<R2> int_f(mesh->numberOfEdges()); + auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix(); + + parallel_for( + mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) { + auto edge_node_list = edge_to_node_matrix[edge_id]; + auto xr = mesh->xr(); + R2 integral = zero; + + LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); + } + + int_f[edge_id] = integral; + }); - SECTION("edge list") + return int_f; + }(); + + SECTION(mesh_name) + { + SECTION("direct formula") { - SECTION("Array") + SECTION("all edges") { - Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + SECTION("EdgeValue") + { + EdgeValue<R2> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfEdges()); + + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + SECTION("SmallArray") { - size_t k = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { - edge_list[k] = edge_id; + SmallArray<R2> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); } - REQUIRE(k == edge_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } + } - Array<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + SECTION("edge list") + { + SECTION("Array") + { + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + Array<R2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < edge_list.size(); ++i) { - error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<R2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } + } - SECTION("SmallArray") + SECTION("tensorial formula") + { + SECTION("all edges") { - SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + SECTION("EdgeValue") + { + EdgeValue<R2> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") { - size_t k = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { - edge_list[k] = edge_id; + Array<R2> values(mesh->numberOfEdges()); + + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); } - REQUIRE(k == edge_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - SmallArray<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } - double error = 0; - for (size_t i = 0; i < edge_list.size(); ++i) { - error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } + } + + SECTION("edge list") + { + SECTION("Array") + { + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + REQUIRE(k == edge_list.size()); + } + + Array<R2> values = EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<R2> values = + EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } } + } + } + } - SECTION("tensorial formula") + SECTION("matrix") + { + using R2x2 = TinyMatrix<2>; + auto l2Norm = [](const R2x2& A) -> double { + return std::sqrt(A(0, 0) * A(0, 0) + A(1, 0) * A(1, 0) + A(0, 1) * A(0, 1) + A(1, 1) * A(1, 1)); + }; + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> R2x2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2x2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x + 3 * y - 1, 2 * x * x + 3 * x * y - z * z, + 3 - 2 * x + 3 * y}; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const R2x2> int_f_per_edge = [=] { + Array<R2x2> int_f(mesh->numberOfEdges()); + auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix(); + + parallel_for( + mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) { + auto edge_node_list = edge_to_node_matrix[edge_id]; + auto xr = mesh->xr(); + R2x2 integral = zero; + + LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); + } + + int_f[edge_id] = integral; + }); + + return int_f; + }(); + + SECTION(mesh_name) { - SECTION("all edges") + SECTION("direct formula") { - SECTION("EdgeValue") + SECTION("all edges") { - EdgeValue<double> values(mesh->connectivity()); - EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), - *mesh, values); + SECTION("EdgeValue") + { + EdgeValue<R2x2> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfEdges()); - SECTION("Array") - { - Array<double> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); - EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("edge list") { - SmallArray<double> values(mesh->numberOfEdges()); - EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + SECTION("Array") + { + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } - double error = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { - error += std::abs(int_f_per_edge[edge_id] - values[edge_id]); + REQUIRE(k == edge_list.size()); + } + + Array<R2x2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<R2x2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } - SECTION("edge list") + SECTION("tensorial formula") { - SECTION("Array") + SECTION("all edges") { - Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; - + SECTION("EdgeValue") { - size_t k = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { - edge_list[k] = edge_id; + EdgeValue<R2x2> values(mesh->connectivity()); + EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); } - REQUIRE(k == edge_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - Array<double> values = - EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfEdges()); + + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } - double error = 0; - for (size_t i = 0; i < edge_list.size(); ++i) { - error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfEdges()); + EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) { + error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("edge list") { - SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; - + SECTION("Array") { - size_t k = 0; - for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { - edge_list[k] = edge_id; + Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); } - REQUIRE(k == edge_list.size()); - } + Array<R2x2> values = + EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); - SmallArray<double> values = - EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < edge_list.size(); ++i) { - error += std::abs(int_f_per_edge[edge_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2}; + + { + size_t k = 0; + for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) { + edge_list[k] = edge_id; + } + + REQUIRE(k == edge_list.size()); + } + + SmallArray<R2x2> values = + EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list); + + double error = 0; + for (size_t i = 0; i < edge_list.size(); ++i) { + error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } } diff --git a/tests/test_FaceIntegrator.cpp b/tests/test_FaceIntegrator.cpp index 184f1a8f97db137f99f638f5feef048f4d78ec04..945d765a6d0a530c7c7babb5d8b5800d5af033d3 100644 --- a/tests/test_FaceIntegrator.cpp +++ b/tests/test_FaceIntegrator.cpp @@ -15,487 +15,1477 @@ TEST_CASE("FaceIntegrator", "[scheme]") { - SECTION("2D") + SECTION("scalar") { - using R2 = TinyVector<2>; + SECTION("2D") + { + using R2 = TinyVector<2>; - const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); - auto f = [](const R2& X) -> double { - const double x = X[0]; - const double y = X[1]; - return x * x + 2 * x * y + 3 * y * y + 2; - }; + auto f = [](const R2& X) -> double { + const double x = X[0]; + const double y = X[1]; + return x * x + 2 * x * y + 3 * y * y + 2; + }; - Array<const double> int_f_per_face = [=] { - Array<double> int_f(mesh->numberOfFaces()); - auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + Array<const double> int_f_per_face = [=] { + Array<double> int_f(mesh->numberOfFaces()); + auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); - parallel_for( - mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { - auto face_node_list = face_to_node_matrix[face_id]; - auto xr = mesh->xr(); - double integral = 0; - LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]); - auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { + auto face_node_list = face_to_node_matrix[face_id]; + auto xr = mesh->xr(); + double integral = 0; + LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); - } + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); + } - int_f[face_id] = integral; - }); + int_f[face_id] = integral; + }); - return int_f; - }(); + return int_f; + }(); - SECTION("direct formula") - { - SECTION("all faces") + SECTION("direct formula") { - SECTION("FaceValue") + SECTION("all faces") { - FaceValue<double> values(mesh->connectivity()); - FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + SECTION("FaceValue") + { + FaceValue<double> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<double> values(mesh->numberOfFaces()); - SECTION("Array") - { - Array<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<double> values(mesh->numberOfFaces()); - FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } - SECTION("face list") + SECTION("tensorial formula") { - SECTION("Array") + SECTION("all faces") { - Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - + SECTION("FaceValue") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + FaceValue<double> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); } - REQUIRE(k == face_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + SECTION("Array") + { + Array<double> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + REQUIRE(k == face_list.size()); + } + + Array<double> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<double> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); } - REQUIRE(k == face_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } + } + } + } + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> double { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const double> int_f_per_face = [=] { + Array<double> int_f(mesh->numberOfFaces()); + auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { + auto face_node_list = face_to_node_matrix[face_id]; + auto xr = mesh->xr(); + double integral = 0; + + switch (face_node_list.size()) { + case 3: { + TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi)); + } + break; + } + case 4: { + SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]], + xr[face_node_list[3]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi)); + } + break; + } + default: { + throw UnexpectedError("invalid face (node number must be 3 or 4)"); + } + } + int_f[face_id] = integral; + }); + + return int_f; + }(); + + SECTION(mesh_name) + { + SECTION("direct formula") + { + SECTION("all faces") + { + SECTION("FaceValue") + { + FaceValue<double> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); - SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<double> values = + FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("tensorial formula") + { + SECTION("all faces") + { + SECTION("FaceValue") + { + FaceValue<double> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<double> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += std::abs(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<double> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<double> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += std::abs(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } } } } + } - SECTION("tensorial formula") + SECTION("R^d") + { + SECTION("2D") { - SECTION("all faces") + using R2 = TinyVector<2>; + + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + auto f = [](const R2& X) -> R2 { + const double x = X[0]; + const double y = X[1]; + return R2{x * x + 2 * x * y + 3 * y * y + 2, 3 * x - 2 * y}; + }; + + Array<const R2> int_f_per_face = [=] { + Array<R2> int_f(mesh->numberOfFaces()); + auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { + auto face_node_list = face_to_node_matrix[face_id]; + auto xr = mesh->xr(); + R2 integral = zero; + LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); + } + + int_f[face_id] = integral; + }); + + return int_f; + }(); + + SECTION("direct formula") { - SECTION("FaceValue") + SECTION("all faces") { - FaceValue<double> values(mesh->connectivity()); - FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, - values); + SECTION("FaceValue") + { + FaceValue<R2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2& x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2> values(mesh->numberOfFaces()); - SECTION("Array") - { - Array<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); - FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<double> values(mesh->numberOfFaces()); - FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } - SECTION("face list") + SECTION("tensorial formula") { - SECTION("Array") + SECTION("all faces") { - Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - + SECTION("FaceValue") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + FaceValue<R2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); } - REQUIRE(k == face_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - Array<double> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + SECTION("Array") + { + Array<R2> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + Array<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); } - REQUIRE(k == face_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } + } + } + } + + SECTION("3D") + { + using R2 = TinyVector<2>; + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> R2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 3 * x - 2 * y + z}; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const R2> int_f_per_face = [=] { + Array<R2> int_f(mesh->numberOfFaces()); + auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { + auto face_node_list = face_to_node_matrix[face_id]; + auto xr = mesh->xr(); + R2 integral = zero; + + switch (face_node_list.size()) { + case 3: { + TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi)); + } + break; + } + case 4: { + SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]], + xr[face_node_list[3]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi)); + } + break; + } + default: { + throw UnexpectedError("invalid face (node number must be 3 or 4)"); + } + } + int_f[face_id] = integral; + }); - SmallArray<double> values = - FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + return int_f; + }(); + + SECTION(mesh_name) + { + SECTION("direct formula") + { + SECTION("all faces") + { + SECTION("FaceValue") + { + FaceValue<R2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("tensorial formula") + { + SECTION("all faces") + { + SECTION("FaceValue") + { + FaceValue<R2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("Array") + { + Array<R2> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + Array<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + } } } } } - SECTION("3D") + SECTION("R^dxd") { - using R3 = TinyVector<3>; - - auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); - - auto f = [](const R3& X) -> double { - const double x = X[0]; - const double y = X[1]; - const double z = X[2]; - return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1; + using R2x2 = TinyMatrix<2>; + auto l2Norm = [](const R2x2& A) -> double { + return std::sqrt(A(0, 0) * A(0, 0) + A(1, 0) * A(1, 0) + A(0, 1) * A(0, 1) + A(1, 1) * A(1, 1)); }; - std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; - mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); - mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + SECTION("2D") + { + using R2 = TinyVector<2>; + + const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh(); - for (auto mesh_info : mesh_list) { - auto mesh_name = mesh_info.first; - auto mesh = mesh_info.second; + auto f = [](const R2& X) -> R2x2 { + const double x = X[0]; + const double y = X[1]; + return R2x2{x * x + 2 * x * y + 3 * y * y + 2, 3 * x - 2 * y, 2 * x * x - 2 * y, 2 + x * y}; + }; - Array<const double> int_f_per_face = [=] { - Array<double> int_f(mesh->numberOfFaces()); + Array<const R2x2> int_f_per_face = [=] { + Array<R2x2> int_f(mesh->numberOfFaces()); auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); parallel_for( mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { auto face_node_list = face_to_node_matrix[face_id]; auto xr = mesh->xr(); - double integral = 0; + R2x2 integral = zero; + LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]); + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); - switch (face_node_list.size()) { - case 3: { - TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]); - auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi)); - } - break; - } - case 4: { - SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]], - xr[face_node_list[3]]); - auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2)); - for (size_t i = 0; i < qf.numberOfPoints(); ++i) { - const auto& xi = qf.point(i); - integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi)); - } - break; - } - default: { - throw UnexpectedError("invalid face (node number must be 3 or 4)"); - } + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.velocityNorm() * f(T(xi)); } + int_f[face_id] = integral; }); return int_f; }(); - SECTION(mesh_name) + SECTION("direct formula") { - SECTION("direct formula") + SECTION("all faces") { - SECTION("all faces") + SECTION("FaceValue") { - SECTION("FaceValue") - { - FaceValue<double> values(mesh->connectivity()); - FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, - values); + FaceValue<R2x2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2& x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values); - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); - } + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); } - SECTION("Array") - { - Array<double> values(mesh->numberOfFaces()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } - FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + REQUIRE(k == face_list.size()); } - SECTION("SmallArray") - { - SmallArray<double> values(mesh->numberOfFaces()); - FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + Array<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + REQUIRE(k == face_list.size()); + } + + SmallArray<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } + } + } - SECTION("face list") + SECTION("tensorial formula") + { + SECTION("all faces") + { + SECTION("FaceValue") { - SECTION("Array") - { - Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + FaceValue<R2x2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh, + values); - { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; - } + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - REQUIRE(k == face_list.size()); - } + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfFaces()); - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); - } + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); } - SECTION("SmallArray") - { - SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; - } + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } + } + + SECTION("face list") + { + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - REQUIRE(k == face_list.size()); + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; } - SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + REQUIRE(k == face_list.size()); + } + + Array<R2x2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + REQUIRE(k == face_list.size()); } + + SmallArray<R2x2> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } } + } + } + + SECTION("3D") + { + using R3 = TinyVector<3>; + + auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + auto f = [](const R3& X) -> R2x2 { + const double x = X[0]; + const double y = X[1]; + const double z = X[2]; + return R2x2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 3 * x - 2 * y + z, 2 * x - y * z, 3 * z - x * x}; + }; + + std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list; + mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh)); + mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh))); + + for (auto mesh_info : mesh_list) { + auto mesh_name = mesh_info.first; + auto mesh = mesh_info.second; + + Array<const R2x2> int_f_per_face = [=] { + Array<R2x2> int_f(mesh->numberOfFaces()); + auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) { + auto face_node_list = face_to_node_matrix[face_id]; + auto xr = mesh->xr(); + R2x2 integral = zero; + + switch (face_node_list.size()) { + case 3: { + TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]); + auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi)); + } + break; + } + case 4: { + SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]], + xr[face_node_list[3]]); + auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2)); + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& xi = qf.point(i); + integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi)); + } + break; + } + default: { + throw UnexpectedError("invalid face (node number must be 3 or 4)"); + } + } + int_f[face_id] = integral; + }); + + return int_f; + }(); - SECTION("tensorial formula") + SECTION(mesh_name) { - SECTION("all faces") + SECTION("direct formula") { - SECTION("FaceValue") + SECTION("all faces") { - FaceValue<double> values(mesh->connectivity()); - FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), - *mesh, values); + SECTION("FaceValue") + { + FaceValue<R2x2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, + values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); - } + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfFaces()); - SECTION("Array") - { - Array<double> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); - FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<double> values(mesh->numberOfFaces()); - FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + SECTION("Array") + { + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } - double error = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - error += std::abs(int_f_per_face[face_id] - values[face_id]); + REQUIRE(k == face_list.size()); + } + + Array<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } - SECTION("face list") + SECTION("tensorial formula") { - SECTION("Array") + SECTION("all faces") { - Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - + SECTION("FaceValue") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + FaceValue<R2x2> values(mesh->connectivity()); + FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10), + *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); } - REQUIRE(k == face_list.size()); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - Array<double> values = - FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + SECTION("Array") + { + Array<R2x2> values(mesh->numberOfFaces()); + + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<R2x2> values(mesh->numberOfFaces()); + FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values); + + double error = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + error += l2Norm(int_f_per_face[face_id] - values[face_id]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } - SECTION("SmallArray") + SECTION("face list") { - SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; - + SECTION("Array") { - size_t k = 0; - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { - face_list[k] = face_id; + Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); } - REQUIRE(k == face_list.size()); - } + Array<R2x2> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); - SmallArray<double> values = - FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } - double error = 0; - for (size_t i = 0; i < face_list.size(); ++i) { - error += std::abs(int_f_per_face[face_list[i]] - values[i]); + REQUIRE(error == Catch::Approx(0).margin(1E-10)); } - REQUIRE(error == Catch::Approx(0).margin(1E-10)); + SECTION("SmallArray") + { + SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2}; + + { + size_t k = 0; + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) { + face_list[k] = face_id; + } + + REQUIRE(k == face_list.size()); + } + + SmallArray<R2x2> values = + FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list); + + double error = 0; + for (size_t i = 0; i < face_list.size(); ++i) { + error += l2Norm(int_f_per_face[face_list[i]] - values[i]); + } + + REQUIRE(error == Catch::Approx(0).margin(1E-10)); + } } } } diff --git a/tests/test_GaussLegendreQuadratureDescriptor.cpp b/tests/test_GaussLegendreQuadratureDescriptor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8c28ee825b734b484a4f4ce7c1da425365342454 --- /dev/null +++ b/tests/test_GaussLegendreQuadratureDescriptor.cpp @@ -0,0 +1,16 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <analysis/GaussLegendreQuadratureDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("GaussLegendreQuadratureDescriptor", "[analysis]") +{ + GaussLegendreQuadratureDescriptor quadrature_descriptor(3); + + REQUIRE(quadrature_descriptor.isTensorial() == true); + REQUIRE(quadrature_descriptor.type() == QuadratureType::GaussLegendre); + REQUIRE(quadrature_descriptor.degree() == 3); + REQUIRE(quadrature_descriptor.name() == ::name(QuadratureType::GaussLegendre) + "(3)"); +} diff --git a/tests/test_GaussLobattoQuadratureDescriptor.cpp b/tests/test_GaussLobattoQuadratureDescriptor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fe5430735bd2cf2596ba63fd1857ec931f8a1bb3 --- /dev/null +++ b/tests/test_GaussLobattoQuadratureDescriptor.cpp @@ -0,0 +1,16 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <analysis/GaussLobattoQuadratureDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("GaussLobattoQuadratureDescriptor", "[analysis]") +{ + GaussLobattoQuadratureDescriptor quadrature_descriptor(3); + + REQUIRE(quadrature_descriptor.isTensorial() == true); + REQUIRE(quadrature_descriptor.type() == QuadratureType::GaussLobatto); + REQUIRE(quadrature_descriptor.degree() == 3); + REQUIRE(quadrature_descriptor.name() == ::name(QuadratureType::GaussLobatto) + "(3)"); +} diff --git a/tests/test_GaussQuadratureDescriptor.cpp b/tests/test_GaussQuadratureDescriptor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..033ac7b4a6ebc2d78f0777f254ddf2a2eff659bb --- /dev/null +++ b/tests/test_GaussQuadratureDescriptor.cpp @@ -0,0 +1,16 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <analysis/GaussQuadratureDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("GaussQuadratureDescriptor", "[analysis]") +{ + GaussQuadratureDescriptor quadrature_descriptor(3); + + REQUIRE(quadrature_descriptor.isTensorial() == false); + REQUIRE(quadrature_descriptor.type() == QuadratureType::Gauss); + REQUIRE(quadrature_descriptor.degree() == 3); + REQUIRE(quadrature_descriptor.name() == ::name(QuadratureType::Gauss) + "(3)"); +} diff --git a/tests/test_ItemId.cpp b/tests/test_ItemId.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cad3328b06c75d5e57a468b5d7f7efe8a9b36215 --- /dev/null +++ b/tests/test_ItemId.cpp @@ -0,0 +1,143 @@ +#include <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <mesh/ItemId.hpp> + +// Instantiate to ensure full coverage is performed +template class ItemIdT<ItemType::node>; +template class ItemIdT<ItemType::edge>; +template class ItemIdT<ItemType::face>; +template class ItemIdT<ItemType::cell>; + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("ItemId", "[mesh]") +{ + SECTION("NodeId") + { + NodeId node_id0(13); + REQUIRE(node_id0 == 13); + + NodeId node_id1 = 0; + REQUIRE(node_id1 == 0); + + REQUIRE(++node_id1 == 1); + + REQUIRE(node_id1++ == 1); + REQUIRE(node_id1 == 2); + + NodeId node_id2{node_id1}; + REQUIRE(node_id2 == 2); + + node_id2 = 17; + REQUIRE(node_id2 == 17); + + node_id2 = node_id1; + REQUIRE(node_id2 == 2); + + NodeId node_id3; + { + NodeId tmp_node_id{12}; + NodeId node_id4{std::move(tmp_node_id)}; + REQUIRE(node_id4 == 12); + node_id3 = std::move(node_id4); + } + REQUIRE(node_id3 == 12); + } + + SECTION("EdgeId") + { + EdgeId edge_id0(13); + REQUIRE(edge_id0 == 13); + + EdgeId edge_id1 = 0; + REQUIRE(edge_id1 == 0); + + REQUIRE(++edge_id1 == 1); + + REQUIRE(edge_id1++ == 1); + REQUIRE(edge_id1 == 2); + + EdgeId edge_id2{edge_id1}; + REQUIRE(edge_id2 == 2); + + edge_id2 = 17; + REQUIRE(edge_id2 == 17); + + edge_id2 = edge_id1; + REQUIRE(edge_id2 == 2); + + EdgeId edge_id3; + { + EdgeId tmp_edge_id{12}; + EdgeId edge_id4{std::move(tmp_edge_id)}; + REQUIRE(edge_id4 == 12); + edge_id3 = std::move(edge_id4); + } + REQUIRE(edge_id3 == 12); + } + + SECTION("FaceId") + { + FaceId face_id0(13); + REQUIRE(face_id0 == 13); + + FaceId face_id1 = 0; + REQUIRE(face_id1 == 0); + + REQUIRE(++face_id1 == 1); + + REQUIRE(face_id1++ == 1); + REQUIRE(face_id1 == 2); + + FaceId face_id2{face_id1}; + REQUIRE(face_id2 == 2); + + face_id2 = 17; + REQUIRE(face_id2 == 17); + + face_id2 = face_id1; + REQUIRE(face_id2 == 2); + + FaceId face_id3; + { + FaceId tmp_face_id{12}; + FaceId face_id4{std::move(tmp_face_id)}; + REQUIRE(face_id4 == 12); + face_id3 = std::move(face_id4); + } + REQUIRE(face_id3 == 12); + } + + SECTION("CellId") + { + CellId cell_id0(13); + REQUIRE(cell_id0 == 13); + + CellId cell_id1 = 0; + REQUIRE(cell_id1 == 0); + + REQUIRE(++cell_id1 == 1); + + REQUIRE(cell_id1++ == 1); + REQUIRE(cell_id1 == 2); + + CellId cell_id2{cell_id1}; + REQUIRE(cell_id2 == 2); + + cell_id2 = 17; + REQUIRE(cell_id2 == 17); + + cell_id2 = cell_id1; + REQUIRE(cell_id2 == 2); + + CellId cell_id3; + { + CellId tmp_cell_id{12}; + CellId cell_id4{std::move(tmp_cell_id)}; + REQUIRE(cell_id4 == 12); + cell_id3 = std::move(cell_id4); + } + REQUIRE(cell_id3 == 12); + } +} diff --git a/tests/test_PrismGaussQuadrature.cpp b/tests/test_PrismGaussQuadrature.cpp index ceff287e5ac5c823c9f427ba6801babfeaa47481..000edeed6a40556d702bd440a26aab90f37b7e5d 100644 --- a/tests/test_PrismGaussQuadrature.cpp +++ b/tests/test_PrismGaussQuadrature.cpp @@ -682,5 +682,9 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]") SECTION("max implemented degree") { REQUIRE(QuadratureManager::instance().maxPrismDegree(QuadratureType::Gauss) == PrismGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getPrismFormula( + GaussQuadratureDescriptor(PrismGaussQuadrature ::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(PrismGaussQuadrature ::max_degree) + " on prisms"); } } diff --git a/tests/test_PyramidGaussQuadrature.cpp b/tests/test_PyramidGaussQuadrature.cpp index e648c96fb13df144945e7c73cdf381545ff1169b..ea022b1360908ccd3f270d3afda53b710a486488 100644 --- a/tests/test_PyramidGaussQuadrature.cpp +++ b/tests/test_PyramidGaussQuadrature.cpp @@ -566,5 +566,9 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]") { REQUIRE(QuadratureManager::instance().maxPyramidDegree(QuadratureType::Gauss) == PyramidGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getPyramidFormula( + GaussQuadratureDescriptor(PyramidGaussQuadrature ::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(PyramidGaussQuadrature ::max_degree) + " on pyramids"); } } diff --git a/tests/test_QuadratureType.cpp b/tests/test_QuadratureType.cpp new file mode 100644 index 0000000000000000000000000000000000000000..89519aa97562e6b73647f826930ffc92add1ad76 --- /dev/null +++ b/tests/test_QuadratureType.cpp @@ -0,0 +1,16 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <analysis/QuadratureType.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("QuadratureType", "[analysis]") +{ + SECTION("name") + { + REQUIRE(name(QuadratureType::Gauss) == "Gauss"); + REQUIRE(name(QuadratureType::GaussLegendre) == "Gauss-Legendre"); + REQUIRE(name(QuadratureType::GaussLobatto) == "Gauss-Lobatto"); + } +} diff --git a/tests/test_RefId.cpp b/tests/test_RefId.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f14f42826ff1a4234136f72b2332b8734f272ba0 --- /dev/null +++ b/tests/test_RefId.cpp @@ -0,0 +1,69 @@ +#include <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <mesh/RefId.hpp> + +#include <sstream> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("RefId", "[mesh]") +{ + SECTION("tag with name") + { + RefId ref_id0(3, "tag_name"); + REQUIRE(ref_id0.tagNumber() == 3); + REQUIRE(ref_id0.tagName() == "tag_name"); + + RefId ref_id1(ref_id0); + REQUIRE(ref_id1.tagNumber() == 3); + REQUIRE(ref_id1.tagName() == "tag_name"); + REQUIRE(ref_id0 == ref_id1); + + RefId ref_id2(2, "another_tag_name"); + REQUIRE(ref_id2 < ref_id1); + + RefId ref_id3; + ref_id3 = ref_id0; + REQUIRE(ref_id3 == ref_id0); + + RefId ref_id4 = std::move(ref_id3); + REQUIRE(ref_id4 == ref_id0); + + RefId ref_id5(std::move(ref_id4)); + REQUIRE(ref_id5 == ref_id0); + + std::stringstream os; + os << ref_id0; + REQUIRE(os.str() == "tag_name(3)"); + } + + SECTION("tag with no name") + { + RefId ref_id0(3); + REQUIRE(ref_id0.tagNumber() == 3); + REQUIRE(ref_id0.tagName() == "3"); + + RefId ref_id1(ref_id0); + REQUIRE(ref_id1.tagNumber() == 3); + REQUIRE(ref_id1.tagName() == "3"); + REQUIRE(ref_id0 == ref_id1); + + RefId ref_id2(2); + REQUIRE(ref_id2 < ref_id1); + + RefId ref_id3; + ref_id3 = ref_id0; + REQUIRE(ref_id3 == ref_id0); + + RefId ref_id4 = std::move(ref_id3); + REQUIRE(ref_id4 == ref_id0); + + RefId ref_id5(std::move(ref_id4)); + REQUIRE(ref_id5 == ref_id0); + + std::stringstream os; + os << ref_id0; + REQUIRE(os.str() == "3"); + } +} diff --git a/tests/test_RefItemList.cpp b/tests/test_RefItemList.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4a9a7b7b2dc3a7aad46fd29703769d4550fc67dd --- /dev/null +++ b/tests/test_RefItemList.cpp @@ -0,0 +1,137 @@ +#include <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <mesh/RefItemList.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("RefItemList", "[mesh]") +{ + SECTION("nodes") + { + const Array<NodeId> node_id_array = convert_to_array(std::vector<NodeId>{1, 3, 7, 2, 4, 11}); + const RefId ref_id{3, "my_reference"}; + + RefItemList<ItemType::node> ref_node_list{ref_id, node_id_array}; + REQUIRE(ref_node_list.refId() == ref_id); + REQUIRE(ref_node_list.list().size() == node_id_array.size()); + REQUIRE(&(ref_node_list.list()[0]) == &(node_id_array[0])); + + { + RefItemList copy_ref_node_list{ref_node_list}; + REQUIRE(copy_ref_node_list.refId() == ref_id); + REQUIRE(copy_ref_node_list.list().size() == node_id_array.size()); + REQUIRE(&(copy_ref_node_list.list()[0]) == &(node_id_array[0])); + } + + { + RefItemList<ItemType::node> affect_ref_node_list; + affect_ref_node_list = ref_node_list; + REQUIRE(affect_ref_node_list.refId() == ref_id); + REQUIRE(affect_ref_node_list.list().size() == node_id_array.size()); + REQUIRE(&(affect_ref_node_list.list()[0]) == &(node_id_array[0])); + + RefItemList<ItemType::node> move_ref_node_list; + move_ref_node_list = std::move(affect_ref_node_list); + REQUIRE(move_ref_node_list.refId() == ref_id); + REQUIRE(move_ref_node_list.list().size() == node_id_array.size()); + REQUIRE(&(move_ref_node_list.list()[0]) == &(node_id_array[0])); + } + } + + SECTION("edges") + { + const Array<EdgeId> edge_id_array = convert_to_array(std::vector<EdgeId>{1, 3, 7, 2, 4, 11}); + const RefId ref_id{3, "my_reference"}; + + RefItemList<ItemType::edge> ref_edge_list{ref_id, edge_id_array}; + REQUIRE(ref_edge_list.refId() == ref_id); + REQUIRE(ref_edge_list.list().size() == edge_id_array.size()); + REQUIRE(&(ref_edge_list.list()[0]) == &(edge_id_array[0])); + + { + RefItemList copy_ref_edge_list{ref_edge_list}; + REQUIRE(copy_ref_edge_list.refId() == ref_id); + REQUIRE(copy_ref_edge_list.list().size() == edge_id_array.size()); + REQUIRE(&(copy_ref_edge_list.list()[0]) == &(edge_id_array[0])); + } + + { + RefItemList<ItemType::edge> affect_ref_edge_list; + affect_ref_edge_list = ref_edge_list; + REQUIRE(affect_ref_edge_list.refId() == ref_id); + REQUIRE(affect_ref_edge_list.list().size() == edge_id_array.size()); + REQUIRE(&(affect_ref_edge_list.list()[0]) == &(edge_id_array[0])); + + RefItemList<ItemType::edge> move_ref_edge_list; + move_ref_edge_list = std::move(affect_ref_edge_list); + REQUIRE(move_ref_edge_list.refId() == ref_id); + REQUIRE(move_ref_edge_list.list().size() == edge_id_array.size()); + REQUIRE(&(move_ref_edge_list.list()[0]) == &(edge_id_array[0])); + } + } + + SECTION("faces") + { + const Array<FaceId> face_id_array = convert_to_array(std::vector<FaceId>{1, 3, 7, 2, 4, 11}); + const RefId ref_id{3, "my_reference"}; + + RefItemList<ItemType::face> ref_face_list{ref_id, face_id_array}; + REQUIRE(ref_face_list.refId() == ref_id); + REQUIRE(ref_face_list.list().size() == face_id_array.size()); + REQUIRE(&(ref_face_list.list()[0]) == &(face_id_array[0])); + + { + RefItemList copy_ref_face_list{ref_face_list}; + REQUIRE(copy_ref_face_list.refId() == ref_id); + REQUIRE(copy_ref_face_list.list().size() == face_id_array.size()); + REQUIRE(&(copy_ref_face_list.list()[0]) == &(face_id_array[0])); + } + + { + RefItemList<ItemType::face> affect_ref_face_list; + affect_ref_face_list = ref_face_list; + REQUIRE(affect_ref_face_list.refId() == ref_id); + REQUIRE(affect_ref_face_list.list().size() == face_id_array.size()); + REQUIRE(&(affect_ref_face_list.list()[0]) == &(face_id_array[0])); + + RefItemList<ItemType::face> move_ref_face_list; + move_ref_face_list = std::move(affect_ref_face_list); + REQUIRE(move_ref_face_list.refId() == ref_id); + REQUIRE(move_ref_face_list.list().size() == face_id_array.size()); + REQUIRE(&(move_ref_face_list.list()[0]) == &(face_id_array[0])); + } + } + + SECTION("cells") + { + const Array<CellId> cell_id_array = convert_to_array(std::vector<CellId>{1, 3, 7, 2, 4, 11}); + const RefId ref_id{3, "my_reference"}; + + RefItemList<ItemType::cell> ref_cell_list{ref_id, cell_id_array}; + REQUIRE(ref_cell_list.refId() == ref_id); + REQUIRE(ref_cell_list.list().size() == cell_id_array.size()); + REQUIRE(&(ref_cell_list.list()[0]) == &(cell_id_array[0])); + + { + RefItemList copy_ref_cell_list{ref_cell_list}; + REQUIRE(copy_ref_cell_list.refId() == ref_id); + REQUIRE(copy_ref_cell_list.list().size() == cell_id_array.size()); + REQUIRE(&(copy_ref_cell_list.list()[0]) == &(cell_id_array[0])); + } + + { + RefItemList<ItemType::cell> affect_ref_cell_list; + affect_ref_cell_list = ref_cell_list; + REQUIRE(affect_ref_cell_list.refId() == ref_id); + REQUIRE(affect_ref_cell_list.list().size() == cell_id_array.size()); + REQUIRE(&(affect_ref_cell_list.list()[0]) == &(cell_id_array[0])); + + RefItemList<ItemType::cell> move_ref_cell_list; + move_ref_cell_list = std::move(affect_ref_cell_list); + REQUIRE(move_ref_cell_list.refId() == ref_id); + REQUIRE(move_ref_cell_list.list().size() == cell_id_array.size()); + REQUIRE(&(move_ref_cell_list.list()[0]) == &(cell_id_array[0])); + } + } +} diff --git a/tests/test_SocketModule.cpp b/tests/test_SocketModule.cpp index 879595965742b5aaa100784037506e58db22150f..8a3c3a3edecd3c3f28e2bb6506a8c05205a7dce0 100644 --- a/tests/test_SocketModule.cpp +++ b/tests/test_SocketModule.cpp @@ -48,7 +48,9 @@ TEST_CASE("SocketModule", "[language]") auto connect_socket_server = [&name_builtin_function](const std::string& hostname, int64_t port) { auto i_function = name_builtin_function.find("connectSocketServer:string*N"); - REQUIRE(i_function != name_builtin_function.end()); + if (i_function == name_builtin_function.end()) { + FAIL_CHECK("Cannot connect to server. This should NEVER happen"); + } DataVariant hostname_variant = hostname; DataVariant port_variant = port; @@ -56,8 +58,6 @@ TEST_CASE("SocketModule", "[language]") IBuiltinFunctionEmbedder& function_embedder = *i_function->second; DataVariant result_variant = function_embedder.apply({hostname_variant, port_variant}); - REQUIRE_NOTHROW(dynamic_cast<const DataHandler<const Socket>&>(std::get<EmbeddedData>(result_variant).get())); - return dynamic_cast<const DataHandler<const Socket>&>(std::get<EmbeddedData>(result_variant).get()).data_ptr(); }; diff --git a/tests/test_SquareGaussQuadrature.cpp b/tests/test_SquareGaussQuadrature.cpp index 5e6b9d6e5ee066cd8db48c6a2297bfc3d12cfcd8..74f442a3154e23e4436c9b170a9649e72a61d4ae 100644 --- a/tests/test_SquareGaussQuadrature.cpp +++ b/tests/test_SquareGaussQuadrature.cpp @@ -459,6 +459,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]") SECTION("max implemented degree") { REQUIRE(QuadratureManager::instance().maxSquareDegree(QuadratureType::Gauss) == SquareGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula( + GaussQuadratureDescriptor(SquareGaussQuadrature ::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(SquareGaussQuadrature ::max_degree) + " on squares"); } SECTION("Access functions") diff --git a/tests/test_TensorialGaussLegendreQuadrature.cpp b/tests/test_TensorialGaussLegendreQuadrature.cpp index aa8708ec610e0c43fa9ad2ded1f50c78e5461b91..342a4fbaa9ef8c54e274f7ad65800aee30986559 100644 --- a/tests/test_TensorialGaussLegendreQuadrature.cpp +++ b/tests/test_TensorialGaussLegendreQuadrature.cpp @@ -3,6 +3,7 @@ #include <catch2/matchers/catch_matchers_all.hpp> #include <analysis/GaussLegendreQuadratureDescriptor.hpp> +#include <analysis/GaussQuadratureDescriptor.hpp> #include <analysis/QuadratureManager.hpp> #include <analysis/TensorialGaussLegendreQuadrature.hpp> #include <utils/Exceptions.hpp> @@ -519,6 +520,15 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]") { REQUIRE(QuadratureManager::instance().maxLineDegree(QuadratureType::GaussLegendre) == TensorialGaussLegendreQuadrature<1>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula( + GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<1>::max_degree + 1)), + "error: Gauss-Legendre quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines"); + + REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula( + GaussQuadratureDescriptor(TensorialGaussLegendreQuadrature<1>::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines"); } } @@ -569,6 +579,16 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]") REQUIRE(integrate(px7y6, l7, 0, 1, 0.2, 0.8) == Catch::Approx(std::pow(0.8, 7) - std::pow(0.2, 7))); } + + SECTION("max implemented degree") + { + REQUIRE(QuadratureManager::instance().maxSquareDegree(QuadratureType::GaussLegendre) == + TensorialGaussLegendreQuadrature<2>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula( + GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<2>::max_degree + 1)), + "error: Gauss-Legendre quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLegendreQuadrature<2>::max_degree) + " on squares"); + } } SECTION("3D") @@ -624,5 +644,15 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]") Catch::Approx((std::pow(0.8, 7) - std::pow(0.2, 7)) * (0.7 - -0.1) + (0.8 - 0.2) * (std::pow(0.7, 8) - std::pow(-0.1, 8)))); } + + SECTION("max implemented degree") + { + REQUIRE(QuadratureManager::instance().maxCubeDegree(QuadratureType::GaussLegendre) == + TensorialGaussLegendreQuadrature<3>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula( + GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<3>::max_degree + 1)), + "error: Gauss-Legendre quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLegendreQuadrature<3>::max_degree) + " on cubes"); + } } } diff --git a/tests/test_TensorialGaussLobattoQuadrature.cpp b/tests/test_TensorialGaussLobattoQuadrature.cpp index e038b692ea6dce069e037640307e9d57864fffef..9ee4028cf0aa9b58d9f7478e91f2146f7f77a675 100644 --- a/tests/test_TensorialGaussLobattoQuadrature.cpp +++ b/tests/test_TensorialGaussLobattoQuadrature.cpp @@ -294,6 +294,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]") { REQUIRE(QuadratureManager::instance().maxLineDegree(QuadratureType::GaussLobatto) == TensorialGaussLobattoQuadrature<1>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula( + GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<1>::max_degree + 1)), + "error: Gauss-Lobatto quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree) + " on lines"); } } @@ -344,6 +348,16 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]") REQUIRE(integrate(px7y6, l7, 0, 1, 0.2, 0.8) == Catch::Approx(std::pow(0.8, 7) - std::pow(0.2, 7))); } + + SECTION("max implemented degree") + { + REQUIRE(QuadratureManager::instance().maxSquareDegree(QuadratureType::GaussLobatto) == + TensorialGaussLobattoQuadrature<2>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula( + GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<2>::max_degree + 1)), + "error: Gauss-Lobatto quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLobattoQuadrature<2>::max_degree) + " on squares"); + } } SECTION("3D") @@ -399,6 +413,16 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]") Catch::Approx((std::pow(0.8, 7) - std::pow(0.2, 7)) * (0.7 - -0.1) + (0.8 - 0.2) * (std::pow(0.7, 8) - std::pow(-0.1, 8)))); } + + SECTION("max implemented degree") + { + REQUIRE(QuadratureManager::instance().maxCubeDegree(QuadratureType::GaussLobatto) == + TensorialGaussLobattoQuadrature<3>::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula( + GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<3>::max_degree + 1)), + "error: Gauss-Lobatto quadrature formulae handle degrees up to " + + std::to_string(TensorialGaussLobattoQuadrature<3>::max_degree) + " on cubes"); + } } SECTION("Access functions") diff --git a/tests/test_TetrahedronGaussQuadrature.cpp b/tests/test_TetrahedronGaussQuadrature.cpp index 605d1dc1fa7476719fe16ac595ff55491f5bd8e8..a412c13688596e37916fd50a522b00a6aaae8485 100644 --- a/tests/test_TetrahedronGaussQuadrature.cpp +++ b/tests/test_TetrahedronGaussQuadrature.cpp @@ -679,5 +679,9 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]") { REQUIRE(QuadratureManager::instance().maxTetrahedronDegree(QuadratureType::Gauss) == TetrahedronGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getTetrahedronFormula( + GaussQuadratureDescriptor(TetrahedronGaussQuadrature ::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(TetrahedronGaussQuadrature ::max_degree) + " on tetrahedra"); } } diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp index 3d7b5dffaf2f0290d5b15dcb6796188512f62cf8..4ed18d60284861241861c45dd744a8e1b38c794b 100644 --- a/tests/test_TriangleGaussQuadrature.cpp +++ b/tests/test_TriangleGaussQuadrature.cpp @@ -624,5 +624,9 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]") { REQUIRE(QuadratureManager::instance().maxTriangleDegree(QuadratureType::Gauss) == TriangleGaussQuadrature::max_degree); + REQUIRE_THROWS_WITH(QuadratureManager::instance().getTriangleFormula( + GaussQuadratureDescriptor(TriangleGaussQuadrature::max_degree + 1)), + "error: Gauss quadrature formulae handle degrees up to " + + std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles"); } }