From 94d6b7d4a1c75566874e258171f6d6a1c09f1172 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Sun, 4 Dec 2022 15:57:33 +0100
Subject: [PATCH] Add sum_of_Vh function

This function computes the sum of the components of a
DiscreteFunctionP0Vector into a DiscreteFunctionP0
---
 doc/userdoc.org                               |   5 +
 .../modules/MathFunctionRegisterForVh.cpp     |   8 +
 ...EmbeddedIDiscreteFunctionMathFunctions.cpp |  34 +++
 ...EmbeddedIDiscreteFunctionMathFunctions.hpp |   2 +
 src/scheme/DiscreteFunctionP0Vector.hpp       |  20 ++
 tests/test_DiscreteFunctionP0Vector.cpp       | 252 ++++++++++++++++++
 ...EmbeddedIDiscreteFunctionMathFunctions.cpp |  99 +++++++
 7 files changed, 420 insertions(+)

diff --git a/doc/userdoc.org b/doc/userdoc.org
index 653821530..da3ac7c01 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -3218,6 +3218,11 @@ Here is the list of the functions
 - ~tan: Vh -> Vh~
 - ~tanh: Vh -> Vh~
 
+Additionally the following function is defined from
+$\vec{\mathbb{P}}_0(\mathbb{R})$ to $\mathbb{P}_0(\mathbb{R})$
+- ~sum_to_Vh: Vh -> Vh~
+it computes the sum of all the components of the vector function.
+
 ******  ~Vh*Vh -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index b15e9796f..99d9bfcec 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -323,6 +323,14 @@ 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("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 f49c4c033..25ec92c29 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -776,6 +776,40 @@ 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));
+  }
+}
+
 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 582e63d85..360175e3c 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -65,6 +65,8 @@ 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>&);
+
 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 8eaa5d0a7..97593e1b3 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 6e351c404..119a7a58a 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 a7f0cc0d7..5370595fa 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -538,6 +538,39 @@ 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("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1048,6 +1081,39 @@ 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("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1558,6 +1624,39 @@ 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("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
-- 
GitLab