Plug doubleDot for discrete P0 functions of Rdxd to the language

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)

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>&);
 

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,   //