diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
index 105dd0ff9ef05a9af8bf6389ac886f810421f6d1..cfa65b1787090a06c6de0ac6bf774b3771a4a1eb 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -109,7 +109,7 @@ DiscreteFunctionVectorIntegrator::_integrate() const
     default: {
       std::ostringstream os;
       os << "vector functions require scalar value type.\n"
-         << "Invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+         << "Invalid value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
       throw NormalError(os.str());
     }
     }
@@ -121,7 +121,7 @@ std::shared_ptr<IDiscreteFunction>
 DiscreteFunctionVectorIntegrator::integrate() const
 {
   if (m_discrete_function_descriptor->type() != DiscreteFunctionType::P0Vector) {
-    throw NormalError("invalid discrete function type for vector interpolation");
+    throw NormalError("invalid discrete function type for vector integration");
   }
 
   switch (m_mesh->dimension()) {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5035890f67986977626aaf13dbae968269aa66b6..4e39374a30ca179ae879bf522b604c5ddc96598b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -163,6 +163,7 @@ add_executable (mpi_unit_tests
   test_DiscreteFunctionInterpolerByZone.cpp
   test_DiscreteFunctionP0.cpp
   test_DiscreteFunctionP0Vector.cpp
+  test_DiscreteFunctionVectorIntegrator.cpp
   test_DiscreteFunctionVectorInterpoler.cpp
   test_DiscreteFunctionVectorInterpolerByZone.cpp
   test_EmbeddedIDiscreteFunctionMathFunctions.cpp
diff --git a/tests/test_DiscreteFunctionVectorIntegrator.cpp b/tests/test_DiscreteFunctionVectorIntegrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d0a9d6a78fd260c0b259aa630d9e79f85271ae89
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorIntegrator.cpp
@@ -0,0 +1,410 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/IntegrateCellValue.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorIntegrator", "[scheme]")
+{
+  auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+    for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+      if (fi[cell_id] != f[cell_id][i]) {
+        return false;
+      }
+    }
+
+    return true;
+  };
+
+  auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+                              const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+                              std::vector<FunctionSymbolId>& function_id_list) {
+    auto [i_symbol, found] = symbol_table->find(name, position);
+    REQUIRE(found);
+    REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+    FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+    function_id_list.push_back(function_symbol_id);
+  };
+
+  SECTION("1D")
+  {
+    constexpr size_t Dimension = 1;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_1d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_1d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    constexpr size_t Dimension = 2;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + x[1]) + 3;
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_2d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    constexpr size_t Dimension = 3;
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * x[2];
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        std::vector<FunctionSymbolId> function_id_list;
+        register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+        register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+        DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                    std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                    function_id_list);
+        std::shared_ptr discrete_function = integrator.integrate();
+
+        size_t i = 0;
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        {
+          CellValue<double> cell_value =
+            IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d);
+
+          REQUIRE(
+            same_cell_value(cell_value, i++,
+                            dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+        }
+
+        REQUIRE(i == function_id_list.size());
+      }
+    }
+  }
+
+  SECTION("errors")
+  {
+    std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+    std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+    for (auto named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_3d = named_mesh.mesh();
+
+        std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0] + x[1]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * x[2];
+let R2_scalar_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0] + x[1]) + 3 * x[2], x[0] - x[1]);
+)";
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        SECTION("invalid function type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+          register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                      std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: invalid function type
+note: expecting R^3 -> R
+note: provided function B_scalar_non_linear_2d: R^2 -> B)";
+
+          REQUIRE_THROWS_WITH(integrator.integrate(), error_msg);
+        }
+
+        SECTION("invalid value type")
+        {
+          std::vector<FunctionSymbolId> function_id_list;
+          register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+          register_function(position, symbol_table, "R2_scalar_non_linear_3d", function_id_list);
+
+          DiscreteFunctionVectorIntegrator integrator(mesh_3d, quadrature_descriptor,
+                                                      std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+                                                      function_id_list);
+
+          const std::string error_msg = R"(error: vector functions require scalar value type.
+Invalid value type: R^2)";
+
+          REQUIRE_THROWS_WITH(integrator.integrate(), error_msg);
+        }
+      }
+    }
+  }
+}