diff --git a/tests/test_EmbeddedIDiscreteFunctionOperators.cpp b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
index dc413efce603bf8512d2b1fcaaee35a7a3170fa0..5232e8fe8d05aca29d5d56440071c0718a206627 100644
--- a/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
@@ -90,12 +90,15 @@
                                                                                            \
     const auto& fuv = dynamic_cast<const DiscreteFunctionType&>(*p_fuv);                   \
                                                                                            \
-    auto lhs_values = P_LHS->cellArrays();                                                 \
-    auto rhs_values = P_RHS->cellArrays();                                                 \
+    auto lhs_arrays = P_LHS->cellArrays();                                                 \
+    auto rhs_arrays = P_RHS->cellArrays();                                                 \
     bool is_same    = true;                                                                \
-    for (CellId cell_id = 0; cell_id < lhs_values.numberOfItems(); ++cell_id) {            \
+    REQUIRE(rhs_arrays.sizeOfArrays() > 0);                                                \
+    REQUIRE(lhs_arrays.sizeOfArrays() == rhs_arrays.sizeOfArrays());                       \
+    REQUIRE(lhs_arrays.sizeOfArrays() == fuv.size());                                      \
+    for (CellId cell_id = 0; cell_id < lhs_arrays.numberOfItems(); ++cell_id) {            \
       for (size_t i = 0; i < fuv.size(); ++i) {                                            \
-        if (fuv[cell_id][i] != (lhs_values[cell_id][i] OPERATOR rhs_values[cell_id][i])) { \
+        if (fuv[cell_id][i] != (lhs_arrays[cell_id][i] OPERATOR rhs_arrays[cell_id][i])) { \
           is_same = false;                                                                 \
           break;                                                                           \
         }                                                                                  \
@@ -116,11 +119,13 @@
                                                                                      \
     const auto& fuv = dynamic_cast<const DiscreteFunctionType&>(*p_fuv);             \
                                                                                      \
-    auto rhs_values = P_RHS->cellArrays();                                           \
+    auto rhs_arrays = P_RHS->cellArrays();                                           \
     bool is_same    = true;                                                          \
-    for (CellId cell_id = 0; cell_id < rhs_values.numberOfItems(); ++cell_id) {      \
+    REQUIRE(rhs_arrays.sizeOfArrays() > 0);                                          \
+    REQUIRE(fuv.size() == rhs_arrays.sizeOfArrays());                                \
+    for (CellId cell_id = 0; cell_id < rhs_arrays.numberOfItems(); ++cell_id) {      \
       for (size_t i = 0; i < fuv.size(); ++i) {                                      \
-        if (fuv[cell_id][i] != (LHS OPERATOR rhs_values[cell_id][i])) {              \
+        if (fuv[cell_id][i] != (LHS OPERATOR rhs_arrays[cell_id][i])) {              \
           is_same = false;                                                           \
           break;                                                                     \
         }                                                                            \
@@ -130,6 +135,56 @@
     REQUIRE(is_same);                                                                \
   }
 
+#define CHECK_SCALAR_VH_TO_VH(OPERATOR, P_RHS)                                   \
+  {                                                                              \
+    using DiscreteFunctionType = const std::decay_t<decltype(OPERATOR * P_RHS)>; \
+                                                                                 \
+    std::shared_ptr p_fu = OPERATOR P_RHS;                                       \
+                                                                                 \
+    REQUIRE(p_fu.use_count() > 0);                                               \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fu));           \
+                                                                                 \
+    const auto& fu = dynamic_cast<const DiscreteFunctionType&>(*p_fu);           \
+                                                                                 \
+    auto rhs_values = P_RHS->cellValues();                                       \
+    bool is_same    = true;                                                      \
+    for (CellId cell_id = 0; cell_id < rhs_values.numberOfItems(); ++cell_id) {  \
+      if (fu[cell_id] != (OPERATOR rhs_values[cell_id])) {                       \
+        is_same = false;                                                         \
+        break;                                                                   \
+      }                                                                          \
+    }                                                                            \
+                                                                                 \
+    REQUIRE(is_same);                                                            \
+  }
+
+#define CHECK_VECTOR_VH_TO_VH(OPERATOR, P_RHS)                                   \
+  {                                                                              \
+    using DiscreteFunctionType = const std::decay_t<decltype(OPERATOR * P_RHS)>; \
+                                                                                 \
+    std::shared_ptr p_fu = OPERATOR P_RHS;                                       \
+                                                                                 \
+    REQUIRE(p_fu.use_count() > 0);                                               \
+    REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionType&>(*p_fu));           \
+                                                                                 \
+    const auto& fu = dynamic_cast<const DiscreteFunctionType&>(*p_fu);           \
+                                                                                 \
+    auto rhs_arrays = P_RHS->cellArrays();                                       \
+    REQUIRE(rhs_arrays.sizeOfArrays() > 0);                                      \
+    REQUIRE(rhs_arrays.sizeOfArrays() == fu.size());                             \
+    bool is_same = true;                                                         \
+    for (CellId cell_id = 0; cell_id < rhs_arrays.numberOfItems(); ++cell_id) {  \
+      for (size_t i = 0; i < rhs_arrays.sizeOfArrays(); ++i) {                   \
+        if (fu[cell_id][i] != (OPERATOR rhs_arrays[cell_id][i])) {               \
+          is_same = false;                                                       \
+          break;                                                                 \
+        }                                                                        \
+      }                                                                          \
+    }                                                                            \
+                                                                                 \
+    REQUIRE(is_same);                                                            \
+  }
+
 TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 {
   SECTION("binary operators")
@@ -2150,4 +2205,418 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
       }
     }
   }
+
+  SECTION("unary operators")
+  {
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+
+      using Rd = TinyVector<Dimension>;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();
+
+      CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      CellValue<double> u_R_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;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_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 DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+      }();
+
+      constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+      }();
+
+      constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1], x[0] + x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+      }();
+
+      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] = {2 * xj[cell_id][0] + 1}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+      }();
+
+      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] = {2 * x[0] + 1, 1 - x[1],   //
+                           2 * x[1], -x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+      }();
+
+      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] = {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 DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+      }();
+
+      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 DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      SECTION("unary minus")
+      {
+        SECTION("- Vh -> Vh")
+        {
+          CHECK_SCALAR_VH_TO_VH(-, p_R_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1x1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2x2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3x3_u);
+
+          CHECK_VECTOR_VH_TO_VH(-, p_Vector3_u);
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+
+      using Rd = TinyVector<Dimension>;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();
+
+      CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      CellValue<double> u_R_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;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_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 DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+      }();
+
+      constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+      }();
+
+      constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1], x[0] + x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+      }();
+
+      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] = {2 * xj[cell_id][0] + 1}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+      }();
+
+      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] = {2 * x[0] + 1, 1 - x[1],   //
+                           2 * x[1], -x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+      }();
+
+      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] = {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 DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+      }();
+
+      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 DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      SECTION("unary minus")
+      {
+        SECTION("- Vh -> Vh")
+        {
+          CHECK_SCALAR_VH_TO_VH(-, p_R_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1x1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2x2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3x3_u);
+
+          CHECK_VECTOR_VH_TO_VH(-, p_Vector3_u);
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+
+      using Rd = TinyVector<Dimension>;
+
+      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();
+
+      CellValue<const Rd> xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+      CellValue<double> u_R_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;
+      }();
+
+      std::shared_ptr p_R_u = std::make_shared<const DiscreteFunctionP0<Dimension, double>>(mesh, u_R_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 DiscreteFunctionP0<Dimension, TinyVector<1>>>(mesh, uj);
+      }();
+
+      constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
+      }();
+
+      constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
+        if constexpr (Dimension == 1) {
+          return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+        } else if constexpr (Dimension == 2) {
+          return {x[0], x[1], x[0] + x[1]};
+        } else if constexpr (Dimension == 3) {
+          return {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]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
+      }();
+
+      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] = {2 * xj[cell_id][0] + 1}; });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
+      }();
+
+      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] = {2 * x[0] + 1, 1 - x[1],   //
+                           2 * x[1], -x[0]};
+          });
+
+        return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
+      }();
+
+      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] = {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 DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
+      }();
+
+      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 DiscreteFunctionP0Vector<Dimension, double>>(mesh, uj_vector);
+      }();
+
+      SECTION("unary minus")
+      {
+        SECTION("- Vh -> Vh")
+        {
+          CHECK_SCALAR_VH_TO_VH(-, p_R_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3_u);
+
+          CHECK_SCALAR_VH_TO_VH(-, p_R1x1_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R2x2_u);
+          CHECK_SCALAR_VH_TO_VH(-, p_R3x3_u);
+
+          CHECK_VECTOR_VH_TO_VH(-, p_Vector3_u);
+        }
+      }
+    }
+  }
 }