#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)");
      }
    }
  }
}