diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 2e1efdbe2a9736804f2c7c89409f0ed72113af73..c9577ef8cc653035ecb36e6a278290a99c3fb549 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -407,8 +407,8 @@ SchemeModule::SchemeModule()
                             std::function(
 
                               [](const std::shared_ptr<const IMesh>& mesh,
-                                 const std::shared_ptr<const IDiscreteFunction>& v)
-                                -> std::shared_ptr<const IDiscreteFunction> { return shallowCopy(mesh, v); }
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& v)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> { return shallowCopy(mesh, v); }
 
                               ));
 
diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp
index f9b6873406adc409a438ba2164b6ac637e1fa2d4..225cdaf15e2ca1d87f78093b8707a26ca564126f 100644
--- a/src/scheme/DiscreteFunctionUtils.cpp
+++ b/src/scheme/DiscreteFunctionUtils.cpp
@@ -65,114 +65,57 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
   return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh, discrete_function->cellValues());
 }
 
-template <size_t Dimension>
-std::shared_ptr<const IDiscreteFunction>
-shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
-            const std::shared_ptr<const IDiscreteFunction>& discrete_function)
+template <typename MeshType, typename DiscreteFunctionT>
+std::shared_ptr<const DiscreteFunctionVariant>
+shallowCopy(const std::shared_ptr<const MeshType>& mesh, const DiscreteFunctionT& f)
 {
-  const std::shared_ptr function_mesh =
-    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(discrete_function->mesh());
+  const std::shared_ptr function_mesh = std::dynamic_pointer_cast<const MeshType>(f.mesh());
 
   if (mesh->shared_connectivity() != function_mesh->shared_connectivity()) {
     throw NormalError("cannot shallow copy when connectivity changes");
   }
 
-  switch (discrete_function->descriptor().type()) {
-  case DiscreteFunctionType::P0: {
-    switch (discrete_function->dataType()) {
-    case ASTNodeDataType::double_t: {
-      return shallowCopy(mesh,
-                         std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, double>>(discrete_function));
+  if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
+    if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+      return std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionT(mesh, f.cellValues()));
+    } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+      return std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionT(mesh, f.cellArrays()));
+    } else {
+      throw UnexpectedError("invalid discrete function type");
     }
-    case ASTNodeDataType::vector_t: {
-      switch (discrete_function->dataType().dimension()) {
-      case 1: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(
-                                   discrete_function));
-      }
-      case 2: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(
-                                   discrete_function));
-      }
-      case 3: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(
-                                   discrete_function));
-      }
-        // LCOV_EXCL_START
-      default: {
-        throw UnexpectedError("invalid data vector dimension: " + stringify(discrete_function->dataType().dimension()));
-      }
-        // LCOV_EXCL_STOP
-      }
-    }
-    case ASTNodeDataType::matrix_t: {
-      if (discrete_function->dataType().numberOfRows() != discrete_function->dataType().numberOfColumns()) {
-        // LCOV_EXCL_START
-        throw UnexpectedError(
-          "invalid data matrix dimensions: " + stringify(discrete_function->dataType().numberOfRows()) + "x" +
-          stringify(discrete_function->dataType().numberOfColumns()));
-        // LCOV_EXCL_STOP
+  } else {
+    throw UnexpectedError("invalid mesh types");
+  }
+}
+
+std::shared_ptr<const DiscreteFunctionVariant>
+shallowCopy(const std::shared_ptr<const IMesh>& mesh,
+            const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_variant)
+{
+  return std::visit(
+    [&](auto&& f) {
+      if (mesh == f.mesh()) {
+        return discrete_function_variant;
+      } else if (mesh->dimension() != f.mesh()->dimension()) {
+        throw NormalError("incompatible mesh dimensions");
       }
-      switch (discrete_function->dataType().numberOfRows()) {
+
+      switch (mesh->dimension()) {
       case 1: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(
-                                   discrete_function));
+        return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), f);
       }
       case 2: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(
-                                   discrete_function));
+        return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), f);
       }
       case 3: {
-        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(
-                                   discrete_function));
+        return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), f);
       }
         // LCOV_EXCL_START
       default: {
-        throw UnexpectedError(
-          "invalid data matrix dimensions: " + stringify(discrete_function->dataType().numberOfRows()) + "x" +
-          stringify(discrete_function->dataType().numberOfColumns()));
+        throw UnexpectedError("invalid mesh dimension");
       }
         // LCOV_EXCL_STOP
       }
-    }
-      // LCOV_EXCL_START
-    default: {
-      throw UnexpectedError("invalid kind of P0 function: invalid data type");
-    }
-      // LCOV_EXCL_STOP
-    }
-  }
-    // LCOV_EXCL_START
-  default: {
-    throw UnexpectedError("invalid discretization type");
-  }
-    // LCOV_EXCL_STOP
-  }
-}
-
-std::shared_ptr<const IDiscreteFunction>
-shallowCopy(const std::shared_ptr<const IMesh>& mesh, const std::shared_ptr<const IDiscreteFunction>& discrete_function)
-{
-  if (mesh == discrete_function->mesh()) {
-    return discrete_function;
-  } else if (mesh->dimension() != discrete_function->mesh()->dimension()) {
-    throw NormalError("incompatible mesh dimensions");
-  }
-
-  switch (mesh->dimension()) {
-  case 1: {
-    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), discrete_function);
-  }
-  case 2: {
-    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), discrete_function);
-  }
-  case 3: {
-    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), discrete_function);
-  }
-    // LCOV_EXCL_START
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-    // LCOV_EXCL_STOP
-  }
+    },
+    discrete_function_variant->discreteFunction());
 }
diff --git a/src/scheme/DiscreteFunctionUtils.hpp b/src/scheme/DiscreteFunctionUtils.hpp
index c788e3768fc11266601e9e90614893a0ac7b7388..ed75fe94c3ceb9cfe4b292c33d1ec432cf744ede 100644
--- a/src/scheme/DiscreteFunctionUtils.hpp
+++ b/src/scheme/DiscreteFunctionUtils.hpp
@@ -43,7 +43,8 @@ std::shared_ptr<const IMesh> getCommonMesh(
 
 bool hasSameMesh(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list);
 
-std::shared_ptr<const IDiscreteFunction> shallowCopy(const std::shared_ptr<const IMesh>& mesh,
-                                                     const std::shared_ptr<const IDiscreteFunction>& discrete_function);
+std::shared_ptr<const DiscreteFunctionVariant> shallowCopy(
+  const std::shared_ptr<const IMesh>& mesh,
+  const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function);
 
 #endif   // DISCRETE_FUNCTION_UTILS_HPP
diff --git a/tests/test_DiscreteFunctionUtils.cpp b/tests/test_DiscreteFunctionUtils.cpp
index c1cb07afea956f0ae9f70f6497d08bfda06edc71..f6762403b2637b764ee9f565a16eefe5aacba008 100644
--- a/tests/test_DiscreteFunctionUtils.cpp
+++ b/tests/test_DiscreteFunctionUtils.cpp
@@ -60,7 +60,8 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("scalar function shallow copy")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const double>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, double>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -68,14 +69,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1 function shallow copy")
         {
-          using DataType     = TinyVector<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -83,14 +85,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2 function shallow copy")
         {
-          using DataType     = TinyVector<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -98,14 +101,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3 function shallow copy")
         {
-          using DataType     = TinyVector<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -113,14 +117,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1x1 function shallow copy")
         {
-          using DataType     = TinyMatrix<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -128,14 +133,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2x2 function shallow copy")
         {
-          using DataType     = TinyMatrix<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -143,14 +149,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3x3 function shallow copy")
         {
-          using DataType     = TinyMatrix<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -158,8 +165,24 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
+        }
+
+        SECTION("P0Vector function shallow copy")
+        {
+          using DiscreteFunctionT = DiscreteFunctionP0Vector<Dimension, const double>;
+          std::shared_ptr uh =
+            std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0Vector<Dimension, double>(mesh, 2));
+          std::shared_ptr vh = shallowCopy(mesh, uh);
+
+          REQUIRE(uh == vh);
+
+          std::shared_ptr wh = shallowCopy(mesh_copy, uh);
+
+          REQUIRE(uh != wh);
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]) ==
+                  &(wh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]));
         }
       }
     }
@@ -213,7 +236,8 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("scalar function shallow copy")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const double>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, double>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -221,14 +245,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1 function shallow copy")
         {
-          using DataType     = TinyVector<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -236,14 +261,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2 function shallow copy")
         {
-          using DataType     = TinyVector<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -251,14 +277,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3 function shallow copy")
         {
-          using DataType     = TinyVector<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -266,14 +293,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1x1 function shallow copy")
         {
-          using DataType     = TinyMatrix<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -281,14 +309,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2x2 function shallow copy")
         {
-          using DataType     = TinyMatrix<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -296,14 +325,31 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3x3 function shallow copy")
         {
-          using DataType     = TinyMatrix<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
+          std::shared_ptr vh = shallowCopy(mesh, uh);
+
+          REQUIRE(uh == vh);
+
+          std::shared_ptr wh = shallowCopy(mesh_copy, uh);
+
+          REQUIRE(uh != wh);
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
+        }
+
+        SECTION("P0Vector function shallow copy")
+        {
+          using DiscreteFunctionT = DiscreteFunctionP0Vector<Dimension, const double>;
+          std::shared_ptr uh =
+            std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0Vector<Dimension, double>(mesh, 2));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -311,8 +357,8 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]) ==
+                  &(wh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]));
         }
       }
     }
@@ -366,7 +412,8 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("scalar function shallow copy")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const double>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, double>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -374,14 +421,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1 function shallow copy")
         {
-          using DataType     = TinyVector<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -389,14 +437,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2 function shallow copy")
         {
-          using DataType     = TinyVector<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -404,14 +453,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3 function shallow copy")
         {
-          using DataType     = TinyVector<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyVector<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -419,14 +469,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^1x1 function shallow copy")
         {
-          using DataType     = TinyMatrix<1>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<1>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -434,14 +485,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^2x2 function shallow copy")
         {
-          using DataType     = TinyMatrix<2>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<2>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -449,14 +501,15 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
         }
 
         SECTION("R^3x3 function shallow copy")
         {
-          using DataType     = TinyMatrix<3>;
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh);
+          using DataType          = TinyMatrix<3>;
+          using DiscreteFunctionT = DiscreteFunctionP0<Dimension, const DataType>;
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, DataType>(mesh));
           std::shared_ptr vh = shallowCopy(mesh, uh);
 
           REQUIRE(uh == vh);
@@ -464,8 +517,24 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr wh = shallowCopy(mesh_copy, uh);
 
           REQUIRE(uh != wh);
-          REQUIRE(&(uh->cellValues()[CellId{0}]) ==
-                  &(dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*wh).cellValues()[CellId{0}]));
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellValues()[CellId{0}]) ==
+                  &(wh->get<DiscreteFunctionT>().cellValues()[CellId{0}]));
+        }
+
+        SECTION("P0Vector function shallow copy")
+        {
+          using DiscreteFunctionT = DiscreteFunctionP0Vector<Dimension, const double>;
+          std::shared_ptr uh =
+            std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0Vector<Dimension, double>(mesh, 2));
+          std::shared_ptr vh = shallowCopy(mesh, uh);
+
+          REQUIRE(uh == vh);
+
+          std::shared_ptr wh = shallowCopy(mesh_copy, uh);
+
+          REQUIRE(uh != wh);
+          REQUIRE(&(uh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]) ==
+                  &(wh->get<DiscreteFunctionT>().cellArrays()[CellId{0}][0]));
         }
       }
     }
@@ -487,7 +556,7 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
           std::shared_ptr other_mesh =
             CartesianMeshBuilder{TinyVector<1>{-1}, TinyVector<1>{3}, TinyVector<1, size_t>{19}}.mesh();
 
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+          std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, double>(mesh));
 
           REQUIRE_THROWS_WITH(shallowCopy(other_mesh, uh), "error: cannot shallow copy when connectivity changes");
         }
@@ -501,7 +570,7 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
       std::shared_ptr mesh_1d = MeshDataBaseForTests::get().cartesian1DMesh();
       std::shared_ptr mesh_2d = MeshDataBaseForTests::get().cartesian2DMesh();
 
-      std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh_1d);
+      std::shared_ptr uh = std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, double>(mesh_1d));
 
       REQUIRE_THROWS_WITH(shallowCopy(mesh_2d, uh), "error: incompatible mesh dimensions");
     }