From d8a810211857be87cdce2bb9a460227d39358c66 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 7 Aug 2024 01:02:38 +0200
Subject: [PATCH] Add tensorProduct for DiscreteFunctionP0
The tensor product is performed cell by cell
Also give access to the tensor product in the language. In this
special case vectors have to be of the same dimension since non-square
matrices are not allowed there
Also add related tests
---
.../modules/MathFunctionRegisterForVh.cpp | 70 ++++++++++++
src/language/modules/MathModule.cpp | 15 +++
.../EmbeddedDiscreteFunctionMathFunctions.cpp | 104 ++++++++++++++++++
.../EmbeddedDiscreteFunctionMathFunctions.hpp | 11 ++
src/scheme/DiscreteFunctionP0.hpp | 45 ++++++++
tests/test_BuiltinFunctionProcessor.cpp | 34 ++++++
tests/test_DiscreteFunctionP0.cpp | 47 ++++++++
...mbeddedDiscreteFunctionMathFunctions1D.cpp | 65 +++++++++++
...mbeddedDiscreteFunctionMathFunctions2D.cpp | 66 +++++++++++
...mbeddedDiscreteFunctionMathFunctions3D.cpp | 66 +++++++++++
tests/test_MathModule.cpp | 59 +++++++++-
11 files changed, 581 insertions(+), 1 deletion(-)
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 38dc08b6e..b31fe47d0 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -212,6 +212,76 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
));
+ scheme_module._addBuiltinFunction("tensorProduct", std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<1> b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<2> b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<3>& b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<1> a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<2> a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<3>& a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
scheme_module._addBuiltinFunction("det", std::function(
[](std::shared_ptr<const DiscreteFunctionVariant> A)
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index 4a1e7b35b..af1dfeafb 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -71,6 +71,21 @@ MathModule::MathModule()
return dot(x, y);
}));
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<1> x, const TinyVector<1> y) -> TinyMatrix<1> {
+ return tensorProduct(x, y);
+ }));
+
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<2> x, const TinyVector<2> y) -> TinyMatrix<2> {
+ return tensorProduct(x, y);
+ }));
+
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<3>& x, const TinyVector<3>& y) -> TinyMatrix<3> {
+ return tensorProduct(x, y);
+ }));
+
this->_addBuiltinFunction("det", std::function([](const TinyMatrix<1> A) -> double { return det(A); }));
this->_addBuiltinFunction("det", std::function([](const TinyMatrix<2>& A) -> double { return det(A); }));
diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
index 715de1be2..961cc0e40 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
@@ -315,6 +315,110 @@ template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<2>&
template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<3>&,
const std::shared_ptr<const DiscreteFunctionVariant>&);
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
+{
+ if (not hasSameMesh({f_v, g_v})) {
+ throw NormalError("operands are defined on different meshes");
+ }
+
+ return std::visit(
+ [&](auto&& f, auto&& g) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using TypeOfF = std::decay_t<decltype(f)>;
+ using TypeOfG = std::decay_t<decltype(g)>;
+ if constexpr (not std::is_same_v<TypeOfF, TypeOfG>) {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
+ } else {
+ using DataType = std::decay_t<typename TypeOfF::data_type>;
+ if constexpr (is_discrete_function_P0_v<TypeOfF>) {
+ if constexpr (is_tiny_vector_v<DataType>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct(f, g));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ }
+ },
+ f_v->discreteFunction(), g_v->discreteFunction());
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyVector<VectorDimension>& a)
+{
+ return std::visit(
+ [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
+ using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
+ if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(discrete_function0, a));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
+ }
+ },
+ f->discreteFunction());
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const TinyVector<VectorDimension>& a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
+{
+ return std::visit(
+ [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
+ using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
+ if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(a, discrete_function0));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+ }
+ },
+ f->discreteFunction());
+}
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<1>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<2>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<3>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<1>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<2>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<3>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
std::shared_ptr<const DiscreteFunctionVariant>
det(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
index 66b784921..bcda26113 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
@@ -69,6 +69,17 @@ template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<VectorDimension>&,
const std::shared_ptr<const DiscreteFunctionVariant>&);
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<VectorDimension>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const TinyVector<VectorDimension>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
std::shared_ptr<const DiscreteFunctionVariant> det(const std::shared_ptr<const DiscreteFunctionVariant>&);
std::shared_ptr<const DiscreteFunctionVariant> trace(const std::shared_ptr<const DiscreteFunctionVariant>&);
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index df43f5290..ea0ec56ae 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -606,6 +606,51 @@ class DiscreteFunctionP0
return result;
}
+ template <typename LHSDataType>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(LHSDataType{}, DataType{}))>
+ tensorProduct(const LHSDataType& u, const DiscreteFunctionP0& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<LHSDataType>>);
+ Assert(v.m_cell_values.isBuilt());
+ DiscreteFunctionP0<decltype(tensorProduct(LHSDataType{}, DataType{}))> result{v.m_mesh_v};
+ parallel_for(
+ v.m_mesh_v->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u, v[cell_id]); });
+
+ return result;
+ }
+
+ template <typename RHSDataType>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(DataType{}, RHSDataType{}))>
+ tensorProduct(const DiscreteFunctionP0& u, const RHSDataType& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<RHSDataType>>);
+ Assert(u.m_cell_values.isBuilt());
+ DiscreteFunctionP0<decltype(tensorProduct(DataType{}, RHSDataType{}))> result{u.m_mesh_v};
+ parallel_for(
+ u.m_mesh_v->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u[cell_id], v); });
+
+ return result;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(DataType{}, DataType2{}))>
+ tensorProduct(const DiscreteFunctionP0& u, const DiscreteFunctionP0<DataType2>& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<DataType2>>);
+ Assert(u.m_cell_values.isBuilt());
+ Assert(v.m_cell_values.isBuilt());
+ Assert(u.m_mesh_v->id() == v.m_mesh_v->id());
+ DiscreteFunctionP0<decltype(tensorProduct(DataType{}, DataType2{}))> result{u.m_mesh_v};
+ parallel_for(
+ u.m_mesh_v->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u[cell_id], v[cell_id]); });
+
+ return result;
+ }
+
PUGS_INLINE friend double
min(const DiscreteFunctionP0& f)
{
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index b678612fe..f600bfa44 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -362,6 +362,40 @@ let s:R, s = dot(x,y);
CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
}
+ { // dot
+ tested_function_set.insert("tensorProduct:R^1*R^1");
+ std::string_view data = R"(
+import math;
+let x:R^1, x = -2;
+let y:R^1, y = 4;
+let s:R^1x1, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", TinyMatrix<1>{-2 * 4});
+ }
+
+ { // dot
+ tested_function_set.insert("tensorProduct:R^2*R^2");
+ std::string_view data = R"(
+import math;
+let x:R^2, x = [-2, 3];
+let y:R^2, y = [4, 3];
+let s:R^2x2, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", tensorProduct(TinyVector<2>{-2, 3}, TinyVector<2>{4, 3}));
+ }
+
+ { // dot
+ tested_function_set.insert("tensorProduct:R^3*R^3");
+ std::string_view data = R"(
+import math;
+let x:R^3, x = [-2, 3, 4];
+let y:R^3, y = [4, 3, 5];
+let s:R^3x3, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s",
+ tensorProduct(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
+ }
+
{ // det
tested_function_set.insert("det:R^1x1");
std::string_view data = R"(
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 5b1133b4e..14c96b0db 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -2729,6 +2729,53 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
}
+ SECTION("tensorProduct(uh,hv)")
+ {
+ DiscreteFunctionP0<TinyVector<2>> uh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ uh[cell_id] = TinyVector<2>{x + 1, 2 * x - 3};
+ });
+
+ DiscreteFunctionP0<TinyVector<2>> vh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ vh[cell_id] = TinyVector<2>{2.3 * x, 1 - x};
+ });
+
+ CHECK_STD_BINARY_MATH_FUNCTION(uh, vh, tensorProduct);
+ }
+
+ SECTION("tensorProduct(uh,v)")
+ {
+ DiscreteFunctionP0<TinyVector<2>> uh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ uh[cell_id] = TinyVector<2>{x + 1, 2 * x - 3};
+ });
+
+ const TinyVector<2> v{1, 2};
+
+ CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(uh, v, tensorProduct);
+ }
+
+ SECTION("tensorProduct(u,hv)")
+ {
+ const TinyVector<2> u{3, -2};
+
+ DiscreteFunctionP0<TinyVector<2>> vh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ vh[cell_id] = TinyVector<2>{2.3 * x, 1 - x};
+ });
+
+ CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, tensorProduct);
+ }
+
SECTION("scalar sum")
{
const CellValue<const double> cell_value = positive_function.cellValues();
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
index 5843a8ae7..41f3c0d22 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
@@ -502,6 +502,29 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -675,6 +698,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
index 11582b982..919d4f828 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
@@ -504,6 +504,30 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions2D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -677,6 +701,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions2D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
index edca0a4eb..99994e9d0 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
@@ -502,6 +502,30 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions3D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -675,6 +699,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions3D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index f299d1e36..89a480e1a 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -13,7 +13,7 @@ TEST_CASE("MathModule", "[language]")
MathModule math_module;
const auto& name_builtin_function = math_module.getNameBuiltinFunctionMap();
- REQUIRE(name_builtin_function.size() == 42);
+ REQUIRE(name_builtin_function.size() == 45);
SECTION("Z -> N")
{
@@ -458,6 +458,63 @@ TEST_CASE("MathModule", "[language]")
}
}
+ SECTION("(R^d, R^d) -> R^dxd")
+ {
+ SECTION("tensorProduct:R^1*R^1")
+ {
+ TinyVector<1> arg0{3};
+ TinyVector<1> arg1{2};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^1*R^1");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<1> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+ }
+
+ SECTION("tensorProduct:R^2*R^2")
+ {
+ TinyVector<2> arg0{3, 2};
+ TinyVector<2> arg1{-2, 5};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^2*R^2");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<2> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+ }
+
+ SECTION("tensorProduct:R^3*R^3")
+ {
+ TinyVector<3> arg0{3, 2, 4};
+ TinyVector<3> arg1{-2, 5, 2};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^3*R^3");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<3> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+ }
+ }
+
SECTION("R^dxd -> double")
{
SECTION("det:R^1x1")
--
GitLab