diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 898f2a71c15b07e8bbc8946795f2450a9f315957..d6f9227f4b66e1f1abc442ca36f2cc21a8b9b04b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -10,6 +10,7 @@
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
 #include <scheme/DiscreteFunctionInterpoler.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
@@ -44,6 +45,15 @@ SchemeModule::SchemeModule()
 
                                     ));
 
+  this->_addBuiltinFunction("P0Vector",
+                            std::make_shared<
+                              BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunctionDescriptor>()>>(
+                              []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
+                                return std::make_shared<DiscreteFunctionDescriptorP0Vector>();
+                              }
+
+                              ));
+
   this->_addBuiltinFunction(
     "interpolate",
     std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
@@ -65,7 +75,17 @@ SchemeModule::SchemeModule()
       [](std::shared_ptr<const IMesh> mesh,
          std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
          const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
-        return DiscreteFunctionInterpoler{mesh, discrete_function_descriptor, function_id}.interpolate();
+        switch (discrete_function_descriptor->type()) {
+        case DiscreteFunctionType::P0: {
+          return DiscreteFunctionInterpoler{mesh, discrete_function_descriptor, function_id}.interpolate();
+        }
+        case DiscreteFunctionType::P0Vector: {
+          return DiscreteFunctionVectorInterpoler{mesh, discrete_function_descriptor, {function_id}}.interpolate();
+        }
+        default: {
+          throw NormalError("invalid function descriptor type");
+        }
+        }
       }
 
       ));
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index 8b2d50b12fe85092d837c699e62973bdec43da08..4827345a5864ccae5425f6463e00163d455a4aa8 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -49,6 +49,10 @@ DiscreteFunctionVectorInterpoler::_interpolate() const
 std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionVectorInterpoler::interpolate() const
 {
+  if (m_discrete_function_descriptor->type() != DiscreteFunctionType::P0Vector) {
+    throw NormalError("invalid discrete function type for vector interpolation");
+  }
+
   std::shared_ptr<IDiscreteFunction> discrete_function;
   switch (m_mesh->dimension()) {
   case 1: {
@@ -64,5 +68,4 @@ DiscreteFunctionVectorInterpoler::interpolate() const
     throw UnexpectedError("invalid dimension");
   }
   }
-  return nullptr;
 }