diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c9577ef8cc653035ecb36e6a278290a99c3fb549..913789852ec8b7f12b02a1dc0e1f138abe7dfc6e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -305,10 +305,10 @@ SchemeModule::SchemeModule()
 
   this->_addBuiltinFunction("glace_fluxes", std::function(
 
-                                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                                                 const std::shared_ptr<const IDiscreteFunction>& u,
-                                                 const std::shared_ptr<const IDiscreteFunction>& c,
-                                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                                    bc_descriptor_list)
                                                 -> std::tuple<std::shared_ptr<const ItemValueVariant>,
@@ -324,17 +324,17 @@ SchemeModule::SchemeModule()
   this->_addBuiltinFunction("glace_solver",
                             std::function(
 
-                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                                 const std::shared_ptr<const IDiscreteFunction>& u,
-                                 const std::shared_ptr<const IDiscreteFunction>& E,
-                                 const std::shared_ptr<const IDiscreteFunction>& c,
-                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
-                                 const double& dt)
-                                -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>> {
+                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                   .solver()
                                   .apply(AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
@@ -346,10 +346,10 @@ SchemeModule::SchemeModule()
   this->_addBuiltinFunction("eucclhyd_fluxes",
                             std::function(
 
-                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                                 const std::shared_ptr<const IDiscreteFunction>& u,
-                                 const std::shared_ptr<const IDiscreteFunction>& c,
-                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list)
                                 -> std::tuple<std::shared_ptr<const ItemValueVariant>,
@@ -365,17 +365,17 @@ SchemeModule::SchemeModule()
   this->_addBuiltinFunction("eucclhyd_solver",
                             std::function(
 
-                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
-                                 const std::shared_ptr<const IDiscreteFunction>& u,
-                                 const std::shared_ptr<const IDiscreteFunction>& E,
-                                 const std::shared_ptr<const IDiscreteFunction>& c,
-                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
-                                 const double& dt)
-                                -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>> {
+                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                   .solver()
                                   .apply(AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
@@ -387,15 +387,15 @@ SchemeModule::SchemeModule()
   this->_addBuiltinFunction("apply_acoustic_fluxes",
                             std::function(
 
-                              [](const std::shared_ptr<const IDiscreteFunction>& rho,            //
-                                 const std::shared_ptr<const IDiscreteFunction>& u,              //
-                                 const std::shared_ptr<const IDiscreteFunction>& E,              //
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,      //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,        //
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,        //
                                  const std::shared_ptr<const ItemValueVariant>& ur,              //
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
-                                 const double& dt)
-                                -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>,
-                                              std::shared_ptr<const IDiscreteFunction>> {
+                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return AcousticSolverHandler{getCommonMesh({rho, u, E})}   //
                                   .solver()
                                   .apply_fluxes(dt, rho, u, E, ur, Fjr);
@@ -412,12 +412,13 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("acoustic_dt",
-                            std::function(
+  this->_addBuiltinFunction("acoustic_dt", std::function(
 
-                              [](const std::shared_ptr<const IDiscreteFunction>& c) -> double { return acoustic_dt(c); }
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+                                               return acoustic_dt(c);
+                                             }
 
-                              ));
+                                             ));
 
   this->_addBuiltinFunction("cell_volume",
                             std::function(
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 589b1d380d6ba1fbac10a61290552fe85ab1d193..98e26a68959865636e73356fa442e79357dc8251 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -13,7 +13,6 @@
 #include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <utils/Socket.hpp>
@@ -23,7 +22,7 @@
 
 template <size_t Dimension>
 double
-acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
+acoustic_dt(const DiscreteFunctionP0<Dimension, const double>& c)
 {
   const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh());
 
@@ -38,23 +37,19 @@ acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
 }
 
 double
-acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c)
+acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c)
 {
-  if ((c->descriptor().type() != DiscreteFunctionType::P0) or (c->dataType() != ASTNodeDataType::double_t)) {
-    throw NormalError("invalid discrete function type");
-  }
-
-  std::shared_ptr mesh = c->mesh();
+  std::shared_ptr mesh = getCommonMesh({c});
 
   switch (mesh->dimension()) {
   case 1: {
-    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*c));
+    return acoustic_dt(c->get<DiscreteFunctionP0<1, const double>>());
   }
   case 2: {
-    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*c));
+    return acoustic_dt(c->get<DiscreteFunctionP0<2, const double>>());
   }
   case 3: {
-    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*c));
+    return acoustic_dt(c->get<DiscreteFunctionP0<2, const double>>());
   }
   default: {
     throw UnexpectedError("invalid mesh dimension");
@@ -72,8 +67,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   using MeshType     = Mesh<Connectivity<Dimension>>;
   using MeshDataType = MeshData<Dimension>;
 
-  using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
-  using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
+  using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
+  using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
 
   class FixedBoundaryCondition;
   class PressureBoundaryCondition;
@@ -381,26 +376,26 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
-                 const std::shared_ptr<const IDiscreteFunction>& i_rho,
-                 const std::shared_ptr<const IDiscreteFunction>& i_c,
-                 const std::shared_ptr<const IDiscreteFunction>& i_u,
-                 const std::shared_ptr<const IDiscreteFunction>& i_p,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({i_rho, i_c, i_u, i_p});
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
     if (not i_mesh) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
-    if (not checkDiscretizationType({i_rho, i_c, i_u, i_p}, DiscreteFunctionType::P0)) {
+    if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
       throw NormalError("acoustic solver expects P0 functions");
     }
 
     const MeshType& mesh              = dynamic_cast<const MeshType&>(*i_mesh);
-    const DiscreteScalarFunction& rho = dynamic_cast<const DiscreteScalarFunction&>(*i_rho);
-    const DiscreteScalarFunction& c   = dynamic_cast<const DiscreteScalarFunction&>(*i_c);
-    const DiscreteVectorFunction& u   = dynamic_cast<const DiscreteVectorFunction&>(*i_u);
-    const DiscreteScalarFunction& p   = dynamic_cast<const DiscreteScalarFunction&>(*i_p);
+    const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c   = c_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u   = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p   = p_v->get<DiscreteScalarFunction>();
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
 
@@ -422,14 +417,14 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   }
 
   std::tuple<std::shared_ptr<const IMesh>,
-             std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
-             std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
-             std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
   apply_fluxes(const double& dt,
                const MeshType& mesh,
-               const DiscreteFunctionP0<Dimension, double>& rho,
-               const DiscreteFunctionP0<Dimension, Rd>& u,
-               const DiscreteFunctionP0<Dimension, double>& E,
+               const DiscreteScalarFunction& rho,
+               const DiscreteVectorFunction& u,
+               const DiscreteScalarFunction& E,
                const NodeValue<const Rd>& ur,
                const NodeValuePerCell<const Rd>& Fjr) const
   {
@@ -473,51 +468,51 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     parallel_for(
       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
 
-    return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
-            std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
-            std::make_shared<DiscreteScalarFunction>(new_mesh, new_E)};
+    return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
   }
 
   std::tuple<std::shared_ptr<const IMesh>,
-             std::shared_ptr<const IDiscreteFunction>,
-             std::shared_ptr<const IDiscreteFunction>,
-             std::shared_ptr<const IDiscreteFunction>>
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
   apply_fluxes(const double& dt,
-               const std::shared_ptr<const IDiscreteFunction>& rho,
-               const std::shared_ptr<const IDiscreteFunction>& u,
-               const std::shared_ptr<const IDiscreteFunction>& E,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
   {
-    std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
     if (not i_mesh) {
       throw NormalError("discrete functions are not defined on the same mesh");
     }
 
-    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+    if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
       throw NormalError("acoustic solver expects P0 functions");
     }
 
-    return this->apply_fluxes(dt,                                                  //
-                              dynamic_cast<const MeshType&>(*i_mesh),              //
-                              dynamic_cast<const DiscreteScalarFunction&>(*rho),   //
-                              dynamic_cast<const DiscreteVectorFunction&>(*u),     //
-                              dynamic_cast<const DiscreteScalarFunction&>(*E),     //
-                              ur->get<NodeValue<const Rd>>(),                      //
+    return this->apply_fluxes(dt,                                       //
+                              dynamic_cast<const MeshType&>(*i_mesh),   //
+                              rho_v->get<DiscreteScalarFunction>(),     //
+                              u_v->get<DiscreteVectorFunction>(),       //
+                              E_v->get<DiscreteScalarFunction>(),       //
+                              ur->get<NodeValue<const Rd>>(),           //
                               Fjr->get<NodeValuePerCell<const Rd>>());
   }
 
   std::tuple<std::shared_ptr<const IMesh>,
-             std::shared_ptr<const IDiscreteFunction>,
-             std::shared_ptr<const IDiscreteFunction>,
-             std::shared_ptr<const IDiscreteFunction>>
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
   apply(const SolverType& solver_type,
         const double& dt,
-        const std::shared_ptr<const IDiscreteFunction>& rho,
-        const std::shared_ptr<const IDiscreteFunction>& u,
-        const std::shared_ptr<const IDiscreteFunction>& E,
-        const std::shared_ptr<const IDiscreteFunction>& c,
-        const std::shared_ptr<const IDiscreteFunction>& p,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E,
+        const std::shared_ptr<const DiscreteFunctionVariant>& c,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 489b8e795d843f9677daf562214bd0266f803baa..36b62a32a1442b18cbe014a9e90ee28ffad08c8a 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -5,13 +5,13 @@
 #include <tuple>
 #include <vector>
 
-class IDiscreteFunction;
+class DiscreteFunctionVariant;
 class IBoundaryConditionDescriptor;
 class IMesh;
 class ItemValueVariant;
 class SubItemValuePerItemVariant;
 
-double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
+double acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c);
 
 class AcousticSolverHandler
 {
@@ -29,34 +29,34 @@ class AcousticSolverHandler
                        const std::shared_ptr<const SubItemValuePerItemVariant>>
     compute_fluxes(
       const SolverType& solver_type,
-      const std::shared_ptr<const IDiscreteFunction>& rho,
-      const std::shared_ptr<const IDiscreteFunction>& c,
-      const std::shared_ptr<const IDiscreteFunction>& u,
-      const std::shared_ptr<const IDiscreteFunction>& p,
+      const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+      const std::shared_ptr<const DiscreteFunctionVariant>& c,
+      const std::shared_ptr<const DiscreteFunctionVariant>& u,
+      const std::shared_ptr<const DiscreteFunctionVariant>& p,
       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
     virtual std::tuple<std::shared_ptr<const IMesh>,
-                       std::shared_ptr<const IDiscreteFunction>,
-                       std::shared_ptr<const IDiscreteFunction>,
-                       std::shared_ptr<const IDiscreteFunction>>
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
     apply_fluxes(const double& dt,
-                 const std::shared_ptr<const IDiscreteFunction>& rho,
-                 const std::shared_ptr<const IDiscreteFunction>& u,
-                 const std::shared_ptr<const IDiscreteFunction>& E,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
                  const std::shared_ptr<const ItemValueVariant>& ur,
                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
 
     virtual std::tuple<std::shared_ptr<const IMesh>,
-                       std::shared_ptr<const IDiscreteFunction>,
-                       std::shared_ptr<const IDiscreteFunction>,
-                       std::shared_ptr<const IDiscreteFunction>>
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
     apply(const SolverType& solver_type,
           const double& dt,
-          const std::shared_ptr<const IDiscreteFunction>& rho,
-          const std::shared_ptr<const IDiscreteFunction>& u,
-          const std::shared_ptr<const IDiscreteFunction>& E,
-          const std::shared_ptr<const IDiscreteFunction>& c,
-          const std::shared_ptr<const IDiscreteFunction>& p,
+          const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+          const std::shared_ptr<const DiscreteFunctionVariant>& u,
+          const std::shared_ptr<const DiscreteFunctionVariant>& E,
+          const std::shared_ptr<const DiscreteFunctionVariant>& c,
+          const std::shared_ptr<const DiscreteFunctionVariant>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
     IAcousticSolver()                  = default;
diff --git a/src/scheme/DiscreteFunctionUtils.hpp b/src/scheme/DiscreteFunctionUtils.hpp
index ed75fe94c3ceb9cfe4b292c33d1ec432cf744ede..a75394cef80e862a134a68c38c45be96859369ce 100644
--- a/src/scheme/DiscreteFunctionUtils.hpp
+++ b/src/scheme/DiscreteFunctionUtils.hpp
@@ -10,11 +10,12 @@
 
 PUGS_INLINE
 bool
-checkDiscretizationType(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list,
+checkDiscretizationType(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_list,
                         const DiscreteFunctionType& discrete_function_type)
 {
-  for (const auto& discrete_function : discrete_function_list) {
-    if (discrete_function->descriptor().type() != discrete_function_type) {
+  for (const auto& discrete_function_variant : discrete_function_list) {
+    if (not std::visit([&](auto&& f) { return f.descriptor().type() == discrete_function_type; },
+                       discrete_function_variant->discreteFunction())) {
       return false;
     }
   }
diff --git a/tests/test_DiscreteFunctionUtils.cpp b/tests/test_DiscreteFunctionUtils.cpp
index f6762403b2637b764ee9f565a16eefe5aacba008..b206c64062be67e450bfd4df1a06ea779acb9bc1 100644
--- a/tests/test_DiscreteFunctionUtils.cpp
+++ b/tests/test_DiscreteFunctionUtils.cpp
@@ -40,22 +40,27 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("check discretization type")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-          std::shared_ptr vh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-
-          std::shared_ptr qh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh_copy);
-
-          std::shared_ptr Uh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-          std::shared_ptr Vh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-
-          REQUIRE(checkDiscretizationType({uh}, DiscreteFunctionType::P0));
-          REQUIRE(checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({Uh}, DiscreteFunctionType::P0));
+          DiscreteFunctionP0<Dimension, double> uh(mesh);
+          DiscreteFunctionP0<Dimension, double> vh(mesh);
+          DiscreteFunctionP0<Dimension, double> qh(mesh_copy);
+
+          DiscreteFunctionP0Vector<Dimension, double> Uh(mesh, 3);
+          DiscreteFunctionP0Vector<Dimension, double> Vh(mesh, 3);
+
+          auto uh_v = std::make_shared<DiscreteFunctionVariant>(uh);
+          auto vh_v = std::make_shared<DiscreteFunctionVariant>(vh);
+          auto qh_v = std::make_shared<DiscreteFunctionVariant>(qh);
+          auto Uh_v = std::make_shared<DiscreteFunctionVariant>(Uh);
+          auto Vh_v = std::make_shared<DiscreteFunctionVariant>(Vh);
+
+          REQUIRE(checkDiscretizationType({uh_v}, DiscreteFunctionType::P0));
+          REQUIRE(checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0));
         }
 
         SECTION("scalar function shallow copy")
@@ -216,22 +221,27 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("check discretization type")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-          std::shared_ptr vh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-
-          std::shared_ptr qh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh_copy);
-
-          std::shared_ptr Uh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-          std::shared_ptr Vh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-
-          REQUIRE(checkDiscretizationType({uh}, DiscreteFunctionType::P0));
-          REQUIRE(checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({Uh}, DiscreteFunctionType::P0));
+          DiscreteFunctionP0<Dimension, double> uh(mesh);
+          DiscreteFunctionP0<Dimension, double> vh(mesh);
+          DiscreteFunctionP0<Dimension, double> qh(mesh_copy);
+
+          DiscreteFunctionP0Vector<Dimension, double> Uh(mesh, 3);
+          DiscreteFunctionP0Vector<Dimension, double> Vh(mesh, 3);
+
+          auto uh_v = std::make_shared<DiscreteFunctionVariant>(uh);
+          auto vh_v = std::make_shared<DiscreteFunctionVariant>(vh);
+          auto qh_v = std::make_shared<DiscreteFunctionVariant>(qh);
+          auto Uh_v = std::make_shared<DiscreteFunctionVariant>(Uh);
+          auto Vh_v = std::make_shared<DiscreteFunctionVariant>(Vh);
+
+          REQUIRE(checkDiscretizationType({uh_v}, DiscreteFunctionType::P0));
+          REQUIRE(checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0));
         }
 
         SECTION("scalar function shallow copy")
@@ -392,22 +402,27 @@ TEST_CASE("DiscreteFunctionUtils", "[scheme]")
 
         SECTION("check discretization type")
         {
-          std::shared_ptr uh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-          std::shared_ptr vh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-
-          std::shared_ptr qh = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh_copy);
-
-          std::shared_ptr Uh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-          std::shared_ptr Vh = std::make_shared<DiscreteFunctionP0Vector<Dimension, double>>(mesh, 3);
-
-          REQUIRE(checkDiscretizationType({uh}, DiscreteFunctionType::P0));
-          REQUIRE(checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({uh, vh, qh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0Vector));
-          REQUIRE(not checkDiscretizationType({Uh, Vh}, DiscreteFunctionType::P0));
-          REQUIRE(not checkDiscretizationType({Uh}, DiscreteFunctionType::P0));
+          DiscreteFunctionP0<Dimension, double> uh(mesh);
+          DiscreteFunctionP0<Dimension, double> vh(mesh);
+          DiscreteFunctionP0<Dimension, double> qh(mesh_copy);
+
+          DiscreteFunctionP0Vector<Dimension, double> Uh(mesh, 3);
+          DiscreteFunctionP0Vector<Dimension, double> Vh(mesh, 3);
+
+          auto uh_v = std::make_shared<DiscreteFunctionVariant>(uh);
+          auto vh_v = std::make_shared<DiscreteFunctionVariant>(vh);
+          auto qh_v = std::make_shared<DiscreteFunctionVariant>(qh);
+          auto Uh_v = std::make_shared<DiscreteFunctionVariant>(Uh);
+          auto Vh_v = std::make_shared<DiscreteFunctionVariant>(Vh);
+
+          REQUIRE(checkDiscretizationType({uh_v}, DiscreteFunctionType::P0));
+          REQUIRE(checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({uh_v, vh_v, qh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0Vector));
+          REQUIRE(not checkDiscretizationType({Uh_v, Vh_v}, DiscreteFunctionType::P0));
+          REQUIRE(not checkDiscretizationType({Uh_v}, DiscreteFunctionType::P0));
         }
 
         SECTION("scalar function shallow copy")