From edbcf3d5c251cf522972e4d81e31ac85b4697a52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com> Date: Tue, 25 May 2021 22:01:29 +0200 Subject: [PATCH] Add min/max functions: Vh*Vh->Vh, R*Vh->Vh, Vh*R->Vh, Vh->R It is required that Vh has scalar values. --- .../modules/MathFunctionRegisterForVh.cpp | 54 ++++ ...EmbeddedIDiscreteFunctionMathFunctions.cpp | 246 +++++++++++++++++- ...EmbeddedIDiscreteFunctionMathFunctions.hpp | 18 ++ src/scheme/DiscreteFunctionP0.hpp | 93 +++++++ 4 files changed, 408 insertions(+), 3 deletions(-) diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp index da49a0383..ca47bb9ac 100644 --- a/src/language/modules/MathFunctionRegisterForVh.cpp +++ b/src/language/modules/MathFunctionRegisterForVh.cpp @@ -197,4 +197,58 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module const IDiscreteFunction>(const TinyVector<3>&, std::shared_ptr<const IDiscreteFunction>)>>( [](const TinyVector<3>& a, std::shared_ptr<const IDiscreteFunction> b) -> std::shared_ptr<const IDiscreteFunction> { return dot(a, b); })); + + scheme_module + ._addBuiltinFunction("min", + std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>( + [](std::shared_ptr<const IDiscreteFunction> a) -> double { return min(a); })); + + scheme_module + ._addBuiltinFunction("min", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>)>>( + [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b) + -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); })); + + scheme_module + ._addBuiltinFunction("min", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>( + [](double a, std::shared_ptr<const IDiscreteFunction> b) + -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); })); + + scheme_module + ._addBuiltinFunction("min", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>( + [](std::shared_ptr<const IDiscreteFunction> a, + double b) -> std::shared_ptr<const IDiscreteFunction> { return min(a, b); })); + + scheme_module + ._addBuiltinFunction("max", + std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>( + [](std::shared_ptr<const IDiscreteFunction> a) -> double { return max(a); })); + + scheme_module + ._addBuiltinFunction("max", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>)>>( + [](std::shared_ptr<const IDiscreteFunction> a, std::shared_ptr<const IDiscreteFunction> b) + -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); })); + + scheme_module + ._addBuiltinFunction("max", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(double, std::shared_ptr<const IDiscreteFunction>)>>( + [](double a, std::shared_ptr<const IDiscreteFunction> b) + -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); })); + + scheme_module + ._addBuiltinFunction("max", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>, double)>>( + [](std::shared_ptr<const IDiscreteFunction> a, + double b) -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); })); } diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp index e785033cb..e980c91fa 100644 --- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp +++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp @@ -111,7 +111,7 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<c } } else { std::stringstream os; - os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(g); + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g); throw NormalError(os.str()); } } @@ -139,7 +139,7 @@ atan2(const double a, const std::shared_ptr<const IDiscreteFunction>& f) } } else { std::stringstream os; - os << "incompatible operands type " << operand_type_name(a) << " and " << operand_type_name(f); + os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f); throw NormalError(os.str()); } } @@ -167,7 +167,7 @@ atan2(const std::shared_ptr<const IDiscreteFunction>& f, const double a) } } else { std::stringstream os; - os << "incompatible operands type " << operand_type_name(f) << " and " << operand_type_name(a); + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a); throw NormalError(os.str()); } } @@ -479,3 +479,243 @@ template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<2>&, template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<3>&, const std::shared_ptr<const IDiscreteFunction>&); + +double +min(const std::shared_ptr<const IDiscreteFunction>& f) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return min(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return min(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return min(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + throw NormalError("invalid operand type " + operand_type_name(f)); + } +} + +std::shared_ptr<const IDiscreteFunction> +min(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g) +{ + if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and + (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) { + std::shared_ptr mesh = getCommonMesh({f, g}); + + if (mesh.use_count() == 0) { + throw NormalError("operands are defined on different meshes"); + } + + switch (mesh->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>( + min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>( + min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>( + min(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g); + throw NormalError(os.str()); + } +} + +std::shared_ptr<const IDiscreteFunction> +min(const double a, const std::shared_ptr<const IDiscreteFunction>& f) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>(min(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f); + throw NormalError(os.str()); + } +} + +std::shared_ptr<const IDiscreteFunction> +min(const std::shared_ptr<const IDiscreteFunction>& f, const double a) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>(min(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a); + throw NormalError(os.str()); + } +} + +double +max(const std::shared_ptr<const IDiscreteFunction>& f) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return max(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return max(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return max(dynamic_cast<const DiscreteFunctionType&>(*f)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + throw NormalError("invalid operand type " + operand_type_name(f)); + } +} + +std::shared_ptr<const IDiscreteFunction> +max(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g) +{ + if ((f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) and + (g->dataType() == ASTNodeDataType::double_t and g->descriptor().type() == DiscreteFunctionType::P0)) { + std::shared_ptr mesh = getCommonMesh({f, g}); + + if (mesh.use_count() == 0) { + throw NormalError("operands are defined on different meshes"); + } + + switch (mesh->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>( + max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>( + max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>( + max(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g))); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(g); + throw NormalError(os.str()); + } +} + +std::shared_ptr<const IDiscreteFunction> +max(const double a, const std::shared_ptr<const IDiscreteFunction>& f) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>(max(a, dynamic_cast<const DiscreteFunctionType&>(*f))); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(a) << " and " << operand_type_name(f); + throw NormalError(os.str()); + } +} + +std::shared_ptr<const IDiscreteFunction> +max(const std::shared_ptr<const IDiscreteFunction>& f, const double a) +{ + if (f->dataType() == ASTNodeDataType::double_t and f->descriptor().type() == DiscreteFunctionType::P0) { + switch (f->mesh()->dimension()) { + case 1: { + using DiscreteFunctionType = DiscreteFunctionP0<1, double>; + return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + case 2: { + using DiscreteFunctionType = DiscreteFunctionP0<2, double>; + return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + case 3: { + using DiscreteFunctionType = DiscreteFunctionP0<3, double>; + return std::make_shared<const DiscreteFunctionType>(max(dynamic_cast<const DiscreteFunctionType&>(*f), a)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + } else { + std::stringstream os; + os << "incompatible operand types " << operand_type_name(f) << " and " << operand_type_name(a); + throw NormalError(os.str()); + } +} diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp index 51ece8b82..abcfc50b4 100644 --- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp +++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp @@ -65,4 +65,22 @@ template <size_t VectorDimension> std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&, 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>&, + const std::shared_ptr<const IDiscreteFunction>&); + +std::shared_ptr<const IDiscreteFunction> min(const double, const std::shared_ptr<const IDiscreteFunction>&); + +std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&, const double); + +double max(const std::shared_ptr<const IDiscreteFunction>&); + +std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&); + +std::shared_ptr<const IDiscreteFunction> max(const double, const std::shared_ptr<const IDiscreteFunction>&); + +std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&, const double); + #endif // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp index 1d9493a8b..943b3e6be 100644 --- a/src/scheme/DiscreteFunctionP0.hpp +++ b/src/scheme/DiscreteFunctionP0.hpp @@ -4,6 +4,7 @@ #include <scheme/IDiscreteFunction.hpp> #include <mesh/Connectivity.hpp> +#include <mesh/ItemValueUtils.hpp> #include <mesh/Mesh.hpp> #include <mesh/MeshData.hpp> #include <mesh/MeshDataManager.hpp> @@ -559,6 +560,98 @@ class DiscreteFunctionP0 : public IDiscreteFunction return result; } + PUGS_INLINE friend double + min(const DiscreteFunctionP0& f) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + + return min(f.m_cell_values); + } + + PUGS_INLINE friend DiscreteFunctionP0 + min(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt()); + Assert(f.m_mesh == g.m_mesh); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], g[cell_id]); }); + + return result; + } + + PUGS_INLINE friend DiscreteFunctionP0 + min(const double a, const DiscreteFunctionP0& f) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(a, f[cell_id]); }); + + return result; + } + + PUGS_INLINE friend DiscreteFunctionP0 + min(const DiscreteFunctionP0& f, const double a) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::min(f[cell_id], a); }); + + return result; + } + + PUGS_INLINE friend double + max(const DiscreteFunctionP0& f) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + + return max(f.m_cell_values); + } + + PUGS_INLINE friend DiscreteFunctionP0 + max(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt()); + Assert(f.m_mesh == g.m_mesh); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], g[cell_id]); }); + + return result; + } + + PUGS_INLINE friend DiscreteFunctionP0 + max(const double a, const DiscreteFunctionP0& f) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(a, f[cell_id]); }); + + return result; + } + + PUGS_INLINE friend DiscreteFunctionP0 + max(const DiscreteFunctionP0& f, const double a) + { + static_assert(std::is_arithmetic_v<DataType>); + Assert(f.m_cell_values.isBuilt()); + DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh}; + parallel_for( + f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::max(f[cell_id], a); }); + + return result; + } + DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} {} DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value) -- GitLab