diff --git a/doc/userdoc.org b/doc/userdoc.org
index 6538215304ff3da89bd385e8b7918569490cc85e..da3ac7c0192a1c3fb214a7811a7525bc95e5c888 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 b15e9796f39f71c266741616db2afe1df2d5b0e1..99d9bfcec67af3819bf38236b9f0348149f80baa 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 f49c4c033251bf46d7cd91e7c524af4d146db9d3..25ec92c298c21052e61003d13fe1edaa13106ce1 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 582e63d8551842e9f93e3962a08669ba34baf756..360175e3ca416f6e0058b77076748823759bb346 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 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..5370595face3f88f7dfbffdfe186c53d3b5aa612 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);