From d22cc36758f8c6026ecbe60e31bfd7ff2a71a671 Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Sun, 4 Dec 2022 20:24:30 +0100 Subject: [PATCH] Add the vectorize function This function builds a discrete function vector from a tuple of scalar discrete functions --- doc/userdoc.org | 18 ++-- .../modules/MathFunctionRegisterForVh.cpp | 9 ++ ...EmbeddedIDiscreteFunctionMathFunctions.cpp | 85 +++++++++++++++++++ ...EmbeddedIDiscreteFunctionMathFunctions.hpp | 3 + ...EmbeddedIDiscreteFunctionMathFunctions.cpp | 85 +++++++++++++++++++ 5 files changed, 195 insertions(+), 5 deletions(-) diff --git a/doc/userdoc.org b/doc/userdoc.org index da3ac7c01..9295c530b 100644 --- a/doc/userdoc.org +++ b/doc/userdoc.org @@ -3218,11 +3218,6 @@ 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 @@ -3411,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 99d9bfcec..fcf3b31a7 100644 --- a/src/language/modules/MathFunctionRegisterForVh.cpp +++ b/src/language/modules/MathFunctionRegisterForVh.cpp @@ -331,6 +331,15 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module )); + 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 25ec92c29..f83bd4185 100644 --- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp +++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp @@ -810,6 +810,91 @@ sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>& 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 360175e3c..2369429a6 100644 --- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp +++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp @@ -67,6 +67,9 @@ std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&, 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/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp index 5370595fa..5386d9ae6 100644 --- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp +++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp @@ -571,6 +571,35 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]") 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); @@ -1114,6 +1143,33 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]") 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); @@ -1657,6 +1713,35 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]") 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); -- GitLab