diff --git a/doc/userdoc.org b/doc/userdoc.org
index 6538215304ff3da89bd385e8b7918569490cc85e..9295c530bbf2773a6fcc2fe67e1ec9d07c52fd15 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -3406,6 +3406,19 @@ and the return value is an $\mathbb{R}^{3\times3}$ matrix.
   return value is $$\sum_{j\in\mathcal{J}} A_j.$$
 
 
+***** Utilities
+
+****** ~sum_to_Vh: Vh -> Vh~
+
+This function it computes the sum of all the components of
+$\vec{\mathbb{P}}_0(\mathbb{R})$ function and stores the result into a
+$\mathbb{P}_0(\mathbb{R})$.
+
+****** ~vectorize: (Vh) -> Vh~
+
+This function creates a $\vec{\mathbb{P}}_0(\mathbb{R})$ function using
+a tuple of $\mathbb{P}_0(\mathbb{R})$.
+
 ***** Interpolation and integration functions
 
 These functions are ways to define discrete functions from analytic
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index b15e9796f39f71c266741616db2afe1df2d5b0e1..fcf3b31a73da2c47bb59e27bc169195255effb30 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -323,6 +323,23 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
 
                                                      ));
 
+  scheme_module._addBuiltinFunction("sum_of_Vh",
+                                    std::function(
+
+                                      [](std::shared_ptr<const IDiscreteFunction> a)
+                                        -> std::shared_ptr<const IDiscreteFunction> { return sum_of_Vh_components(a); }
+
+                                      ));
+
+  scheme_module
+    ._addBuiltinFunction("vectorize",
+                         std::function(
+
+                           [](const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+                             -> std::shared_ptr<const IDiscreteFunction> { return vectorize(discrete_function_list); }
+
+                           ));
+
   scheme_module._addBuiltinFunction("integral_of_R", std::function(
 
                                                        [](std::shared_ptr<const IDiscreteFunction> a) -> double {
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index f49c4c033251bf46d7cd91e7c524af4d146db9d3..f83bd418592c29b5cded8b77bda6c717ff3af804 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -776,6 +776,125 @@ sum_of(const std::shared_ptr<const IDiscreteFunction>& f)
   }
 }
 
+std::shared_ptr<const IDiscreteFunction>
+sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+  if (f->descriptor().type() == DiscreteFunctionType::P0Vector) {
+    switch (f->mesh()->dimension()) {
+    case 1: {
+      constexpr size_t Dimension = 1;
+      using DiscreteFunctionType = DiscreteFunctionP0Vector<Dimension, double>;
+      return std::make_shared<const DiscreteFunctionP0<Dimension, double>>(
+        sumOfComponents(dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 2: {
+      constexpr size_t Dimension = 2;
+      using DiscreteFunctionType = DiscreteFunctionP0Vector<Dimension, double>;
+      return std::make_shared<const DiscreteFunctionP0<Dimension, double>>(
+        sumOfComponents(dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+    case 3: {
+      constexpr size_t Dimension = 3;
+      using DiscreteFunctionType = DiscreteFunctionP0Vector<Dimension, double>;
+      return std::make_shared<const DiscreteFunctionP0<Dimension, double>>(
+        sumOfComponents(dynamic_cast<const DiscreteFunctionType&>(*f)));
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("invalid mesh dimension");
+    }
+      // LCOV_EXCL_STOP
+    }
+  } else {
+    throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+vectorize(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+{
+  Assert(discrete_function_list.size() > 0);
+  std::shared_ptr p_i_mesh = [&] {
+    auto i = discrete_function_list.begin();
+
+    std::shared_ptr<const IMesh> i_mesh = (*i)->mesh();
+    ++i;
+    for (; i != discrete_function_list.end(); ++i) {
+      if ((*i)->mesh() != i_mesh) {
+        throw NormalError("discrete functions are not defined on the same mesh");
+      }
+    }
+
+    return i_mesh;
+  }();
+
+  for (auto&& i_discrete_function : discrete_function_list) {
+    if ((i_discrete_function->descriptor().type() != DiscreteFunctionType::P0) or
+        (i_discrete_function->dataType() != ASTNodeDataType::double_t)) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(i_discrete_function));
+    }
+  }
+
+  switch (p_i_mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension       = 1;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+  case 2: {
+    constexpr size_t Dimension       = 2;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+  case 3: {
+    constexpr size_t Dimension       = 3;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+    // LCOV_EXCL_STOP
+  }
+}
+
 template double sum_of<double>(const std::shared_ptr<const IDiscreteFunction>&);
 
 template TinyVector<1> sum_of<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 582e63d8551842e9f93e3962a08669ba34baf756..2369429a6267a8fdbc803d75b4213703dd970a48 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -65,6 +65,11 @@ template <size_t VectorDimension>
 std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
                                              const std::shared_ptr<const IDiscreteFunction>&);
 
+std::shared_ptr<const IDiscreteFunction> sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> vectorize(
+  const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list);
+
 double min(const std::shared_ptr<const IDiscreteFunction>&);
 
 std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&,
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 8eaa5d0a7580ddc65cb17e82d6f5551f86114e1c..97593e1b3c5a9a599fc0c17f7edca7a839f2b5ef 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -195,6 +195,26 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
     return product;
   }
 
+  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+  sumOfComponents(const DiscreteFunctionP0Vector& f)
+  {
+    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+
+    parallel_for(
+      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& f_cell_id = f[cell_id];
+
+        double sum = 0;
+        for (size_t i = 0; i < f.size(); ++i) {
+          sum += f_cell_id[i];
+        }
+
+        result[cell_id] = sum;
+      });
+
+    return result;
+  }
+
   PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
   dot(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
   {
diff --git a/tests/test_DiscreteFunctionP0Vector.cpp b/tests/test_DiscreteFunctionP0Vector.cpp
index 6e351c404a7b5f57d7b39e5a3f70d43cf62a5635..119a7a58a93444e2dba7d21a67d0897c9e47a271 100644
--- a/tests/test_DiscreteFunctionP0Vector.cpp
+++ b/tests/test_DiscreteFunctionP0Vector.cpp
@@ -535,6 +535,258 @@ TEST_CASE("DiscreteFunctionP0Vector", "[scheme]")
     }
   }
 
+  SECTION("functions")
+  {
+    SECTION("1D")
+    {
+      const size_t size = 3;
+
+      constexpr size_t Dimension = 1;
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh = named_mesh.mesh();
+
+          auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id][0]  = 2 * x + 1;
+              f[cell_id][1]  = x * x - 1;
+              f[cell_id][2]  = 2 + x;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              g[cell_id][0]  = (x + 1) * (x - 2) + 1;
+              g[cell_id][1]  = 3 * (x + 2) - 1;
+              g[cell_id][2]  = (x + 3) * 5;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          auto same_data = [](const auto& f, const auto& g) {
+            const size_t number_of_cells = g.size();
+            for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+              if (f[cell_id] != g[cell_id]) {
+                return false;
+              }
+            }
+            return true;
+          };
+
+          SECTION("dot")
+          {
+            Array<double> dot_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i] * g[cell_id][i];
+                }
+                dot_values[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(dot(f, g), dot_values));
+            REQUIRE(same_data(dot(const_f, g), dot_values));
+            REQUIRE(same_data(dot(f, const_g), dot_values));
+            REQUIRE(same_data(dot(const_f, const_g), dot_values));
+          }
+
+          SECTION("sumOfComponents")
+          {
+            Array<double> sum_of_components{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i];
+                }
+                sum_of_components[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(sumOfComponents(f), sum_of_components));
+            REQUIRE(same_data(sumOfComponents(const_f), sum_of_components));
+          }
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      const size_t size = 3;
+
+      constexpr size_t Dimension = 2;
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh = named_mesh.mesh();
+
+          auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id][0]  = 2 * x + 1;
+              f[cell_id][1]  = x * x - 1;
+              f[cell_id][2]  = 2 + x;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              g[cell_id][0]  = (x + 1) * (x - 2) + 1;
+              g[cell_id][1]  = 3 * (x + 2) - 1;
+              g[cell_id][2]  = (x + 3) * 5;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          auto same_data = [](const auto& f, const auto& g) {
+            const size_t number_of_cells = g.size();
+            for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+              if (f[cell_id] != g[cell_id]) {
+                return false;
+              }
+            }
+            return true;
+          };
+
+          SECTION("dot")
+          {
+            Array<double> dot_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i] * g[cell_id][i];
+                }
+                dot_values[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(dot(f, g), dot_values));
+            REQUIRE(same_data(dot(const_f, g), dot_values));
+            REQUIRE(same_data(dot(f, const_g), dot_values));
+            REQUIRE(same_data(dot(const_f, const_g), dot_values));
+          }
+
+          SECTION("sumOfComponents")
+          {
+            Array<double> sum_of_components{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i];
+                }
+                sum_of_components[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(sumOfComponents(f), sum_of_components));
+            REQUIRE(same_data(sumOfComponents(const_f), sum_of_components));
+          }
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      const size_t size = 3;
+
+      constexpr size_t Dimension = 3;
+
+      std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh = named_mesh.mesh();
+
+          auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();
+
+          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              f[cell_id][0]  = 2 * x + 1;
+              f[cell_id][1]  = x * x - 1;
+              f[cell_id][2]  = 2 + x;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
+          parallel_for(
+            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+              const double x = xj[cell_id][0];
+              g[cell_id][0]  = (x + 1) * (x - 2) + 1;
+              g[cell_id][1]  = 3 * (x + 2) - 1;
+              g[cell_id][2]  = (x + 3) * 5;
+            });
+
+          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
+          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};
+
+          auto same_data = [](const auto& f, const auto& g) {
+            const size_t number_of_cells = g.size();
+            for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
+              if (f[cell_id] != g[cell_id]) {
+                return false;
+              }
+            }
+            return true;
+          };
+
+          SECTION("dot")
+          {
+            Array<double> dot_values{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i] * g[cell_id][i];
+                }
+                dot_values[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(dot(f, g), dot_values));
+            REQUIRE(same_data(dot(const_f, g), dot_values));
+            REQUIRE(same_data(dot(f, const_g), dot_values));
+            REQUIRE(same_data(dot(const_f, const_g), dot_values));
+          }
+
+          SECTION("sumOfComponents")
+          {
+            Array<double> sum_of_components{mesh->numberOfCells()};
+            parallel_for(
+              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                double sum = 0;
+                for (size_t i = 0; i < size; ++i) {
+                  sum += f[cell_id][i];
+                }
+                sum_of_components[cell_id] = sum;
+              });
+
+            REQUIRE(same_data(sumOfComponents(f), sum_of_components));
+            REQUIRE(same_data(sumOfComponents(const_f), sum_of_components));
+          }
+        }
+      }
+    }
+  }
+
   SECTION("binary operators")
   {
     SECTION("1D")
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index a7f0cc0d74eb79f0897e2ea2a2a4b7212dbc306e..5386d9ae69f9b7b4fbd093d82e1daaef2cbc81c7 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -538,6 +538,68 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("sum_of_Vh Vh -> Vh")
+        {
+          {
+            auto p_sum_components = sum_of_Vh_components(p_Vector3_u);
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0<Dimension, double>& sum_components =
+              dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_sum_components);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              double sum = 0;
+              for (size_t i = 0; i < p_Vector3_u->size(); ++i) {
+                sum += (*p_Vector3_u)[cell_id][i];
+              }
+
+              direct_sum[cell_id] = sum;
+            }
+
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              if (sum_components[cell_id] != direct_sum[cell_id]) {
+                is_same = false;
+                break;
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1048,6 +1110,66 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("sum_of_Vh Vh -> Vh")
+        {
+          {
+            auto p_sum_components = sum_of_Vh_components(p_Vector3_u);
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0<Dimension, double>& sum_components =
+              dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_sum_components);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              double sum = 0;
+              for (size_t i = 0; i < p_Vector3_u->size(); ++i) {
+                sum += (*p_Vector3_u)[cell_id][i];
+              }
+
+              direct_sum[cell_id] = sum;
+            }
+
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              if (sum_components[cell_id] != direct_sum[cell_id]) {
+                is_same = false;
+                break;
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1558,6 +1680,68 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
         }
 
+        SECTION("sum_of_Vh Vh -> Vh")
+        {
+          {
+            auto p_sum_components = sum_of_Vh_components(p_Vector3_u);
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0<Dimension, double>& sum_components =
+              dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_sum_components);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              double sum = 0;
+              for (size_t i = 0; i < p_Vector3_u->size(); ++i) {
+                sum += (*p_Vector3_u)[cell_id][i];
+              }
+
+              direct_sum[cell_id] = sum;
+            }
+
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              if (sum_components[cell_id] != direct_sum[cell_id]) {
+                is_same = false;
+                break;
+              }
+            }
+
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+        }
+
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);