diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
index 961cc0e4092ba57c001866a52d7aed113e6b0f5d..5d000813aeff1e670a2f547a69827a679cbfa699 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
@@ -315,6 +315,107 @@ 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>
+doubleDot(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_matrix_v<DataType>) {
+ return std::make_shared<DiscreteFunctionVariant>(doubleDot(f, g));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ }
+ },
+ f_v->discreteFunction(), g_v->discreteFunction());
+}
+
+template <size_t MatrixDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyMatrix<MatrixDimension>& 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_matrix_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyMatrix<MatrixDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(doubleDot(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 MatrixDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+doubleDot(const TinyMatrix<MatrixDimension>& 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_matrix_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyMatrix<MatrixDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(doubleDot(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> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyMatrix<1>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyMatrix<2>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyMatrix<3>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
+ const TinyMatrix<1>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
+ const TinyMatrix<2>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
+ const TinyMatrix<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)
diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
index bcda26113a4fac6d72e7b54edfe0f0e8d8d74b2d..35ad4f38a03b1f0dfc9962877fe169c8c3265f71 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> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template <size_t SquareMatrixNbRow>
+std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyMatrix<SquareMatrixNbRow>&);
+
+template <size_t SquareMatrixNbRow>
+std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const TinyMatrix<SquareMatrixNbRow>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>&,
const std::shared_ptr<const DiscreteFunctionVariant>&);
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
index 41f3c0d2233b1884051f0eafd5df637ce37ed7e0..8b0b3fbb7b2475ba885434331607977badc2ce1c 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
@@ -171,6 +171,18 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1x1(mesh, uj));
}();
+ std::shared_ptr p_other_mesh_R1x1_u = std::make_shared<const DiscreteFunctionVariant>(
+ DiscreteFunctionR1x1(other_mesh, p_R1x1_u->get<DiscreteFunctionR1x1>().cellValues()));
+
+ std::shared_ptr p_R1x1_v = [=] {
+ CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
+ parallel_for(
+ vj.numberOfItems(),
+ PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = TinyMatrix<1>{3 * xj[cell_id][0] - 2}; });
+
+ return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR1x1(mesh, vj));
+ }();
+
std::shared_ptr p_R2x2_u = [=] {
CellValue<TinyMatrix<2>> uj{mesh->connectivity()};
parallel_for(
@@ -184,6 +196,22 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2x2(mesh, uj));
}();
+ std::shared_ptr p_other_mesh_R2x2_u = std::make_shared<const DiscreteFunctionVariant>(
+ DiscreteFunctionR2x2(other_mesh, p_R2x2_u->get<DiscreteFunctionR2x2>().cellValues()));
+
+ std::shared_ptr p_R2x2_v = [=] {
+ CellValue<TinyMatrix<2>> vj{mesh->connectivity()};
+ parallel_for(
+ vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<2> x = to_2d(xj[cell_id]);
+
+ vj[cell_id] = TinyMatrix<2>{3 * x[0] + 2, 2 - x[0], //
+ 5 * x[1], x[1]};
+ });
+
+ return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR2x2(mesh, vj));
+ }();
+
std::shared_ptr p_R3x3_u = [=] {
CellValue<TinyMatrix<3>> uj{mesh->connectivity()};
parallel_for(
@@ -198,6 +226,25 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3x3(mesh, uj));
}();
+ std::shared_ptr p_other_mesh_R3x3_u = std::make_shared<const DiscreteFunctionVariant>(
+ DiscreteFunctionR3x3(other_mesh, p_R3x3_u->get<DiscreteFunctionR3x3>().cellValues()));
+
+ std::shared_ptr p_R3x3_v = [=] {
+ CellValue<TinyMatrix<3>> vj{mesh->connectivity()};
+ parallel_for(
+ vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<3> x = to_3d(xj[cell_id]);
+
+ vj[cell_id] = TinyMatrix<3>{3 * x[0] + 1, 2 - x[1], 2,
+ //
+ 2 * x[0], x[2], x[1] - x[2],
+ //
+ 2 * x[1] - x[2], x[1] - 2 * x[2], x[0] - x[2]};
+ });
+
+ return std::make_shared<const DiscreteFunctionVariant>(DiscreteFunctionR3x3(mesh, vj));
+ }();
+
std::shared_ptr p_Vector3_u = [=] {
CellArray<double> uj_vector{mesh->connectivity(), 3};
parallel_for(
@@ -502,6 +549,30 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("doubleDot Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, p_R1x1_v, doubleDot, //
+ DiscreteFunctionR1x1, DiscreteFunctionR1x1, DiscreteFunctionR);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2x2_u, p_R2x2_v, doubleDot, //
+ DiscreteFunctionR2x2, DiscreteFunctionR2x2, DiscreteFunctionR);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3x3_u, p_R3x3_v, doubleDot, //
+ DiscreteFunctionR3x3, DiscreteFunctionR3x3, DiscreteFunctionR);
+
+ REQUIRE_THROWS_WITH(doubleDot(p_R1x1_u, p_other_mesh_R1x1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(doubleDot(p_R2x2_u, p_other_mesh_R2x2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(doubleDot(p_R3x3_u, p_other_mesh_R3x3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(doubleDot(p_R1x1_u, p_R3x3_u),
+ "error: incompatible operand types Vh(P0:R^1x1) and Vh(P0:R^3x3)");
+ REQUIRE_THROWS_WITH(doubleDot(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(doubleDot(p_R1_u, p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(doubleDot(p_R2_u, p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(doubleDot(p_R3_u, p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(doubleDot(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)");
+ }
+
SECTION("tensorProduct Vh*Vh -> Vh")
{
CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //