#include <catch2/catch_test_macros.hpp> #include <catch2/matchers/catch_matchers_all.hpp> #include <MeshDataBaseForTests.hpp> #include <scheme/DiscreteFunctionP0.hpp> #include <language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp> #include <scheme/DiscreteFunctionP0Vector.hpp> #include <scheme/DiscreteFunctionVariant.hpp> #include <test_EmbeddedDiscreteFunctionMathFunctions.hpp> // clazy:excludeall=non-pod-global-static TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]") { constexpr size_t Dimension = 1; using Rd = TinyVector<Dimension>; std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes(); using DiscreteFunctionR = DiscreteFunctionP0<const double>; using DiscreteFunctionR1 = DiscreteFunctionP0<const TinyVector<1>>; using DiscreteFunctionR2 = DiscreteFunctionP0<const TinyVector<2>>; using DiscreteFunctionR3 = DiscreteFunctionP0<const TinyVector<3>>; using DiscreteFunctionR1x1 = DiscreteFunctionP0<const TinyMatrix<1>>; using DiscreteFunctionR2x2 = DiscreteFunctionP0<const TinyMatrix<2>>; using DiscreteFunctionR3x3 = DiscreteFunctionP0<const TinyMatrix<3>>; using DiscreteFunctionVector = DiscreteFunctionP0Vector<const double>; for (const auto& named_mesh : mesh_list) { SECTION(named_mesh.name()) { auto mesh = named_mesh.mesh()->get<Mesh<Dimension>>(); std::shared_ptr other_mesh = std::make_shared<const Mesh<Dimension>>(mesh->shared_connectivity(), mesh->xr()); CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj(); CellValue<double> values = [=] { CellValue<double> build_values{mesh->connectivity()}; parallel_for( build_values.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.2 + std::cos(l2Norm(xj[cell_id])); }); return build_values; }(); CellValue<double> positive_values = [=] { CellValue<double> build_values{mesh->connectivity()}; parallel_for( build_values.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 2 + std::sin(l2Norm(xj[cell_id])); }); return build_values; }(); CellValue<double> bounded_values = [=] { CellValue<double> build_values{mesh->connectivity()}; parallel_for( build_values.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { build_values[cell_id] = 0.9 * std::sin(l2Norm(xj[cell_id])); }); return build_values; }(); std::shared_ptr p_u = std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR(mesh, values)); std::shared_ptr p_other_mesh_u = std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR(other_mesh, values)); std::shared_ptr p_positive_u = std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR(mesh, positive_values)); std::shared_ptr p_bounded_u = std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR(mesh, bounded_values)); std::shared_ptr p_R1_u = [=] { CellValue<TinyVector<1>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id][0] = 2 * xj[cell_id][0] + 1; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1(mesh, uj)); }(); std::shared_ptr p_R1_v = [=] { CellValue<TinyVector<1>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id][0] = xj[cell_id][0] * xj[cell_id][0] + 1; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1(mesh, vj)); }(); std::shared_ptr p_other_mesh_R1_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR1(other_mesh, p_R1_u->get<DiscreteFunctionR1>().cellValues())); constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> { if constexpr (Dimension == 1) { return TinyVector<2>{x[0], 1 + x[0] * x[0]}; } else if constexpr (Dimension == 2) { return TinyVector<2>{x[0], x[1]}; } else if constexpr (Dimension == 3) { return TinyVector<2>{x[0], x[1] + x[2]}; } }; std::shared_ptr p_R2_u = [=] { CellValue<TinyVector<2>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<2> x = to_2d(xj[cell_id]); uj[cell_id] = TinyVector<2>{2 * x[0] + 1, 1 - x[1]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2(mesh, uj)); }(); std::shared_ptr p_R2_v = [=] { CellValue<TinyVector<2>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<2> x = to_2d(xj[cell_id]); vj[cell_id] = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2(mesh, vj)); }(); std::shared_ptr p_other_mesh_R2_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR2(other_mesh, p_R2_u->get<DiscreteFunctionR2>().cellValues())); constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> { if constexpr (Dimension == 1) { return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]}; } else if constexpr (Dimension == 2) { return TinyVector<3>{x[0], x[1], x[0] + x[1]}; } else if constexpr (Dimension == 3) { return TinyVector<3>{x[0], x[1], x[2]}; } }; std::shared_ptr p_R3_u = [=] { CellValue<TinyVector<3>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); uj[cell_id] = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3(mesh, uj)); }(); std::shared_ptr p_R3_v = [=] { CellValue<TinyVector<3>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); vj[cell_id] = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3(mesh, vj)); }(); std::shared_ptr p_other_mesh_R3_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR3(other_mesh, p_R3_u->get<DiscreteFunctionR3>().cellValues())); std::shared_ptr p_R1x1_u = [=] { CellValue<TinyMatrix<1>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1x1(mesh, uj)); }(); std::shared_ptr p_other_mesh_R1x1_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR1x1(other_mesh, p_R1x1_u->get<DiscreteFunctionR1x1>().cellValues())); std::shared_ptr p_R1x1_v = [=] { CellValue<TinyMatrix<1>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = TinyMatrix<1>{3 * xj[cell_id][0] - 2}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1x1(mesh, vj)); }(); std::shared_ptr p_R2x2_u = [=] { CellValue<TinyMatrix<2>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<2> x = to_2d(xj[cell_id]); uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1], // 2 * x[1], -x[0]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2x2(mesh, uj)); }(); std::shared_ptr p_other_mesh_R2x2_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR2x2(other_mesh, p_R2x2_u->get<DiscreteFunctionR2x2>().cellValues())); std::shared_ptr p_R2x2_v = [=] { CellValue<TinyMatrix<2>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<2> x = to_2d(xj[cell_id]); vj[cell_id] = TinyMatrix<2>{3 * x[0] + 2, 2 - x[0], // 5 * x[1], x[1]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2x2(mesh, vj)); }(); std::shared_ptr p_R3x3_u = [=] { CellValue<TinyMatrix<3>> uj{mesh->connectivity()}; parallel_for( uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1, 1 - x[1], 3, // 2 * x[1], -x[0], x[0] - x[1], // 3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3x3(mesh, uj)); }(); std::shared_ptr p_other_mesh_R3x3_u = std::make_shared<const DiscreteFunctionVariant>( DiscreteFunctionR3x3(other_mesh, p_R3x3_u->get<DiscreteFunctionR3x3>().cellValues())); std::shared_ptr p_R3x3_v = [=] { CellValue<TinyMatrix<3>> vj{mesh->connectivity()}; parallel_for( vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); vj[cell_id] = TinyMatrix<3>{3 * x[0] + 1, 2 - x[1], 2, // 2 * x[0], x[2], x[1] - x[2], // 2 * x[1] - x[2], x[1] - 2 * x[2], x[0] - x[2]}; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3x3(mesh, vj)); }(); std::shared_ptr p_Vector3_u = [=] { CellArray<double> uj_vector{mesh->connectivity(), 3}; parallel_for( uj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); uj_vector[cell_id][0] = 2 * x[0] + 1; uj_vector[cell_id][1] = 1 - x[1] * x[2]; uj_vector[cell_id][2] = x[0] + x[2]; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionVector(mesh, uj_vector)); }(); std::shared_ptr p_Vector3_v = [=] { CellArray<double> vj_vector{mesh->connectivity(), 3}; parallel_for( vj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); vj_vector[cell_id][0] = x[0] * x[1] + 1; vj_vector[cell_id][1] = 2 * x[1]; vj_vector[cell_id][2] = x[2] * x[0]; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionVector(mesh, vj_vector)); }(); std::shared_ptr p_Vector2_w = [=] { CellArray<double> wj_vector{mesh->connectivity(), 2}; parallel_for( wj_vector.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { const TinyVector<3> x = to_3d(xj[cell_id]); wj_vector[cell_id][0] = x[0] + x[1] * 2; wj_vector[cell_id][1] = x[0] * x[1]; }); return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionVector(mesh, wj_vector)); }(); SECTION("sqrt Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_positive_u, sqrt, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(sqrt(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("abs Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, abs, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(abs(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("sin Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, sin, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(sin(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("cos Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, cos, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(cos(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("tan Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, tan, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(tan(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("asin Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_bounded_u, asin, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(asin(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("acos Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_bounded_u, acos, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(acos(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("atan Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_bounded_u, atan, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(atan(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("sinh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, sinh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(sinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("cosh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, cosh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(cosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("tanh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, tanh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(tanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("asinh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_positive_u, asinh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(asinh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("acosh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_positive_u, acosh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(acosh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("atanh Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_bounded_u, atanh, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(atanh(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("exp Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_u, exp, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(exp(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("log Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_positive_u, log, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(log(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("atan2 Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, atan2, // DiscreteFunctionR, DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(atan2(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(atan2(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)"); REQUIRE_THROWS_WITH(atan2(p_u, p_other_mesh_u), "error: operands are defined on different meshes"); } SECTION("atan2 Vh*R -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 3.6, atan2, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(atan2(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R"); } SECTION("atan2 R*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.4, p_u, atan2, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(atan2(2.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)"); } SECTION("min Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, min, // DiscreteFunctionR, DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(::min(p_u, p_other_mesh_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(::min(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(::min(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)"); } SECTION("min Vh*R -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, min, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(min(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R"); } SECTION("min R*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, min, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(min(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)"); } SECTION("min Vh -> R") { REQUIRE(min(p_u) == min(p_u->get<DiscreteFunctionR>().cellValues())); REQUIRE_THROWS_WITH(min(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("max Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_u, p_bounded_u, max, // DiscreteFunctionR, DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(::max(p_u, p_other_mesh_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(::max(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(::max(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)"); } SECTION("max Vh*R -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_u, 1.2, max, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(max(p_R1_u, 2.1), "error: incompatible operand types Vh(P0:R^1) and R"); } SECTION("max Vh -> R") { REQUIRE(max(p_u) == max(p_u->get<DiscreteFunctionR>().cellValues())); REQUIRE_THROWS_WITH(max(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); } SECTION("max R*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(0.4, p_u, max, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(max(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)"); } SECTION("pow Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_positive_u, p_bounded_u, pow, // DiscreteFunctionR, DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(pow(p_u, p_other_mesh_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(pow(p_u, p_R1_u), "error: incompatible operand types Vh(P0:R) and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(pow(p_R1_u, p_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R)"); } SECTION("pow Vh*R -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_positive_u, 3.3, pow, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(pow(p_R1_u, 3.1), "error: incompatible operand types Vh(P0:R^1) and R"); } SECTION("pow R*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION(2.1, p_u, pow, // DiscreteFunctionR, DiscreteFunctionR); REQUIRE_THROWS_WITH(pow(3.1, p_R1_u), "error: incompatible operand types R and Vh(P0:R^1)"); } SECTION("dot Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, dot, // DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, dot, // DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, dot, // DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR); { std::shared_ptr p_fuv = ::dot(p_Vector3_u, p_Vector3_v); REQUIRE(p_fuv.use_count() > 0); const DiscreteFunctionVector& u = p_Vector3_u->get<DiscreteFunctionVector>(); const DiscreteFunctionVector& v = p_Vector3_v->get<DiscreteFunctionVector>(); const DiscreteFunctionR& fuv = p_fuv->get<DiscreteFunctionR>(); bool is_same = true; auto u_arrays = u.cellArrays(); auto v_arrays = v.cellArrays(); for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) { using namespace std; double dot_u_v = [&](auto&& a, auto&& b) { double sum = 0; for (size_t i = 0; i < a.size(); ++i) { sum += a[i] * b[i]; } return sum; }(u_arrays[cell_id], v_arrays[cell_id]); if (fuv[cell_id] != dot_u_v) { is_same = false; break; } } REQUIRE(is_same); } REQUIRE_THROWS_WITH(dot(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(dot(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(dot(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(dot(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)"); REQUIRE_THROWS_WITH(dot(p_u, p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(dot(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)"); REQUIRE_THROWS_WITH(dot(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)"); REQUIRE_THROWS_WITH(dot(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)"); REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension"); } SECTION("doubleDot Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, p_R1x1_v, doubleDot, // DiscreteFunctionR1x1, DiscreteFunctionR1x1, DiscreteFunctionR); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, p_R2x2_v, doubleDot, // DiscreteFunctionR2x2, DiscreteFunctionR2x2, DiscreteFunctionR); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, p_R3x3_v, doubleDot, // DiscreteFunctionR3x3, DiscreteFunctionR3x3, DiscreteFunctionR); REQUIRE_THROWS_WITH(doubleDot(p_R1x1_u, p_other_mesh_R1x1_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(doubleDot(p_R2x2_u, p_other_mesh_R2x2_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(doubleDot(p_R3x3_u, p_other_mesh_R3x3_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(doubleDot(p_R1x1_u, p_R3x3_u), "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^3x3)"); REQUIRE_THROWS_WITH(doubleDot(p_u, p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(doubleDot(p_R1_u, p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(doubleDot(p_R2_u, p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(doubleDot(p_R3_u, p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(doubleDot(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("tensorProduct Vh*Vh -> Vh") { CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, // DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, // DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2); CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, // DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3); REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u), "error: operands are defined on different meshes"); REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u), "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)"); REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)"); REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)"); REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)"); } SECTION("det Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, // DiscreteFunctionR1x1, DiscreteFunctionR); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, det, // DiscreteFunctionR2x2, DiscreteFunctionR); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, det, // DiscreteFunctionR3x3, DiscreteFunctionR); REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(det(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("trace Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, trace, // DiscreteFunctionR1x1, DiscreteFunctionR); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, trace, // DiscreteFunctionR2x2, DiscreteFunctionR); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, trace, // DiscreteFunctionR3x3, DiscreteFunctionR); REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(trace(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("inverse Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, inverse, // DiscreteFunctionR1x1, DiscreteFunctionR1x1); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, inverse, // DiscreteFunctionR2x2, DiscreteFunctionR2x2); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, inverse, // DiscreteFunctionR3x3, DiscreteFunctionR3x3); REQUIRE_THROWS_WITH(inverse(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(inverse(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("transpose Vh -> Vh") { CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, transpose, // DiscreteFunctionR1x1, DiscreteFunctionR1x1); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, transpose, // DiscreteFunctionR2x2, DiscreteFunctionR2x2); CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, transpose, // DiscreteFunctionR3x3, DiscreteFunctionR3x3); REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(transpose(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("sum_of_Vh Vh -> Vh") { { auto p_sum_components = sum_of_Vh_components(p_Vector3_u); REQUIRE(p_sum_components.use_count() == 1); const DiscreteFunctionR& sum_components = p_sum_components->get<DiscreteFunctionR>(); const DiscreteFunctionVector& vector3_u = p_Vector3_u->get<DiscreteFunctionVector>(); DiscreteFunctionP0<double> direct_sum(mesh); for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { double sum = 0; for (size_t i = 0; i < vector3_u.size(); ++i) { sum += vector3_u[cell_id][i]; } direct_sum[cell_id] = sum; } bool is_same = true; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { if (sum_components[cell_id] != direct_sum[cell_id]) { is_same = false; break; } } REQUIRE(is_same); } REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); } SECTION("vectorize (Vh) -> Vh") { { std::shared_ptr p_vector3 = vectorize(std::vector{p_u, p_positive_u, p_bounded_u}); REQUIRE(p_vector3.use_count() == 1); const DiscreteFunctionVector vector3 = p_vector3->get<DiscreteFunctionVector>(); const DiscreteFunctionR& u = p_u->get<DiscreteFunctionR>(); const DiscreteFunctionR& positive_u = p_positive_u->get<DiscreteFunctionR>(); const DiscreteFunctionR& bounded_u = p_bounded_u->get<DiscreteFunctionR>(); REQUIRE(vector3.size() == 3); bool is_same = true; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { is_same &= (u[cell_id] == vector3[cell_id][0]); is_same &= (positive_u[cell_id] == vector3[cell_id][1]); is_same &= (bounded_u[cell_id] == vector3[cell_id][2]); } REQUIRE(is_same); } REQUIRE_THROWS_WITH(vectorize(std::vector{p_u, p_other_mesh_u}), "error: discrete functions are not defined on the same mesh"); REQUIRE_THROWS_WITH(vectorize(std::vector{p_R1_u}), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(vectorize(std::vector{p_Vector3_u}), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("dot Vh*Rd -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot, // DiscreteFunctionR1, DiscreteFunctionR); CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), dot, // DiscreteFunctionR2, DiscreteFunctionR); CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), dot, // DiscreteFunctionR3, DiscreteFunctionR); REQUIRE_THROWS_WITH(dot(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(dot(p_Vector3_u, TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0Vector:R) and R^2"); REQUIRE_THROWS_WITH(dot(p_R1_u, (TinyVector<2>{-6, 2})), "error: incompatible operand types Vh(P0:R^1) and R^2"); REQUIRE_THROWS_WITH(dot(p_R2_u, (TinyVector<3>{-1, 5, 2})), "error: incompatible operand types Vh(P0:R^2) and R^3"); REQUIRE_THROWS_WITH(dot(p_R3_u, (TinyVector<1>{-1})), "error: incompatible operand types Vh(P0:R^3) and R^1"); } SECTION("dot Rd*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, dot, // DiscreteFunctionR1, DiscreteFunctionR); CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, dot, // DiscreteFunctionR2, DiscreteFunctionR); CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, dot, // DiscreteFunctionR3, DiscreteFunctionR); REQUIRE_THROWS_WITH(dot(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(dot(TinyVector<2>{1, 2}, p_Vector3_u), "error: incompatible operand types R^2 and Vh(P0Vector:R)"); REQUIRE_THROWS_WITH(dot((TinyVector<2>{-6, 2}), p_R1_u), "error: incompatible operand types R^2 and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(dot((TinyVector<3>{-1, 5, 2}), p_R2_u), "error: incompatible operand types R^3 and Vh(P0:R^2)"); REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)"); } SECTION("tensorProduct Vh*Rd -> Vh") { CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, // DiscreteFunctionR1, DiscreteFunctionR1x1); CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, // DiscreteFunctionR2, DiscreteFunctionR2x2); CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, // DiscreteFunctionR3, DiscreteFunctionR3x3); REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}), "error: incompatible operand types Vh(P0Vector:R) and R^2"); REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})), "error: incompatible operand types Vh(P0:R^1) and R^2"); REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})), "error: incompatible operand types Vh(P0:R^2) and R^3"); REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})), "error: incompatible operand types Vh(P0:R^3) and R^1"); } SECTION("tensorProduct Rd*Vh -> Vh") { CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, // DiscreteFunctionR1, DiscreteFunctionR1x1); CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, // DiscreteFunctionR2, DiscreteFunctionR2x2); CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, // DiscreteFunctionR3, DiscreteFunctionR3x3); REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u), "error: incompatible operand types R^2 and Vh(P0Vector:R)"); REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u), "error: incompatible operand types R^2 and Vh(P0:R^1)"); REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u), "error: incompatible operand types R^3 and Vh(P0:R^2)"); REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)"); } SECTION("sum_of_R* Vh -> R*") { REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues())); REQUIRE(sum_of<TinyVector<1>>(p_R1_u) == sum(p_R1_u->get<DiscreteFunctionR1>().cellValues())); REQUIRE(sum_of<TinyVector<2>>(p_R2_u) == sum(p_R2_u->get<DiscreteFunctionR2>().cellValues())); REQUIRE(sum_of<TinyVector<3>>(p_R3_u) == sum(p_R3_u->get<DiscreteFunctionR3>().cellValues())); REQUIRE(sum_of<TinyMatrix<1>>(p_R1x1_u) == sum(p_R1x1_u->get<DiscreteFunctionR1x1>().cellValues())); REQUIRE(sum_of<TinyMatrix<2>>(p_R2x2_u) == sum(p_R2x2_u->get<DiscreteFunctionR2x2>().cellValues())); REQUIRE(sum_of<TinyMatrix<3>>(p_R3x3_u) == sum(p_R3x3_u->get<DiscreteFunctionR3x3>().cellValues())); REQUIRE_THROWS_WITH(sum_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)"); REQUIRE_THROWS_WITH(sum_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)"); REQUIRE_THROWS_WITH(sum_of<double>(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } SECTION("integral_of_R* Vh -> R*") { auto integrate_locally = [&](const auto& cell_values) { const auto& Vj = MeshDataManager::instance().getMeshData(*mesh).Vj(); using DataType = decltype(double{} * cell_values[CellId{0}]); CellValue<DataType> local_integral{mesh->connectivity()}; parallel_for( local_integral.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { local_integral[cell_id] = Vj[cell_id] * cell_values[cell_id]; }); return local_integral; }; REQUIRE(integral_of<double>(p_u) == sum(integrate_locally(p_u->get<DiscreteFunctionR>().cellValues()))); REQUIRE(integral_of<TinyVector<1>>(p_R1_u) == sum(integrate_locally(p_R1_u->get<DiscreteFunctionR1>().cellValues()))); REQUIRE(integral_of<TinyVector<2>>(p_R2_u) == sum(integrate_locally(p_R2_u->get<DiscreteFunctionR2>().cellValues()))); REQUIRE(integral_of<TinyVector<3>>(p_R3_u) == sum(integrate_locally(p_R3_u->get<DiscreteFunctionR3>().cellValues()))); REQUIRE(integral_of<TinyMatrix<1>>(p_R1x1_u) == sum(integrate_locally(p_R1x1_u->get<DiscreteFunctionR1x1>().cellValues()))); REQUIRE(integral_of<TinyMatrix<2>>(p_R2x2_u) == sum(integrate_locally(p_R2x2_u->get<DiscreteFunctionR2x2>().cellValues()))); REQUIRE(integral_of<TinyMatrix<3>>(p_R3x3_u) == sum(integrate_locally(p_R3x3_u->get<DiscreteFunctionR3x3>().cellValues()))); REQUIRE_THROWS_WITH(integral_of<TinyVector<1>>(p_u), "error: invalid operand type Vh(P0:R)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R1_u), "error: invalid operand type Vh(P0:R^1)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R2_u), "error: invalid operand type Vh(P0:R^2)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R3_u), "error: invalid operand type Vh(P0:R^3)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)"); REQUIRE_THROWS_WITH(integral_of<double>(p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)"); REQUIRE_THROWS_WITH(integral_of<double>(p_Vector3_u), "error: invalid operand type Vh(P0Vector:R)"); } } } }