Elementwise multiplication for Tensor<rank, dim>

diff --git a/doc/news/changes/minor/20190610RolandRichter b/doc/news/changes/minor/20190610RolandRichter
new file mode 100644 (file)
index 0000000..47913c8
--- /dev/null
+++ b/doc/news/changes/minor/20190610RolandRichter
@@ -0,0 +1,4 @@
+Added: The function "Schur-product" was added for the Tensor-class, 
+allowing the entrywise multiplication of tensor objects of general rank
+<br>
+(Roland Richter, 2019/06/10)

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 140b2e1df47d8f435e95dccfb11a4842cbb4fa13..2531af9fad817ff406bffab6f90084e29289beb2 100644 (file)
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -1817,6 +1817,73 @@ DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
   return tmp;
 }
 
+/**
+ * Entrywise multiplication of two tensor objects of rank 0 (i.e. the
+ * multiplication of two scalar values).
+ *
+ * @relatesalso Tensor
+ */
+template <int dim, typename Number, typename OtherNumber>
+inline DEAL_II_ALWAYS_INLINE
+  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
+  schur_product(const Tensor<0, dim, Number> &     src1,
+                const Tensor<0, dim, OtherNumber> &src2)
+{
+  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type> tmp(src1);
+
+  tmp *= src2;
+
+  return tmp;
+}
+
+/**
+ * Entrywise multiplication of two tensor objects of rank 1.
+ *
+ * @relatesalso Tensor
+ */
+template <int dim, typename Number, typename OtherNumber>
+inline DEAL_II_ALWAYS_INLINE
+  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
+  schur_product(const Tensor<1, dim, Number> &     src1,
+                const Tensor<1, dim, OtherNumber> &src2)
+{
+  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type> tmp(src1);
+
+  for (unsigned int i = 0; i < dim; ++i)
+    tmp[i] *= src2[i];
+
+  return tmp;
+}
+
+/**
+ * Entrywise multiplication of two tensor objects of general rank.
+ *
+ * This multiplication is also called "Hadamard-product" (c.f.
+ * https://en.wikipedia.org/wiki/Hadamard_product_(matrices)), and generates a
+ * new tensor of size <rank, dim>:
+ * @f[
+ *   \text{result}_{i, j}
+ *   = \text{left}_{i, j}\cdot
+ *     \text{right}_{i, j}
+ * @f]
+ *
+ * @tparam rank The rank of both tensors.
+ *
+ * @relatesalso Tensor
+ */
+template <int rank, int dim, typename Number, typename OtherNumber>
+inline DEAL_II_ALWAYS_INLINE
+  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
+  schur_product(const Tensor<rank, dim, Number> &     src1,
+                const Tensor<rank, dim, OtherNumber> &src2)
+{
+  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tmp(src1);
+
+  for (unsigned int i = 0; i < dim; ++i)
+    tmp[i] = schur_product(src1[i], src2[i]);
+
+  return tmp;
+}
 
 //@}
 /**

diff --git a/tests/tensors/tensor_10.cc b/tests/tensors/tensor_10.cc
new file mode 100644 (file)
index 0000000..93b5085
--- /dev/null
+++ b/tests/tensors/tensor_10.cc
@@ -0,0 +1,98 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// check schur_product for Tensor<rank, dim> with rank = 0, 1, 2, 3 and dim = 2
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+template <int rank, int dim>
+void
+test_same(const Tensor<rank, dim> &t, const Tensor<rank, dim> &compare)
+{
+  deallog << "Constant rank " << rank << " and dim " << dim << '\t'
+          << schur_product(t, t) << " compare " << compare << std::endl;
+}
+
+template <int rank, int dim>
+void
+test_different(const Tensor<rank, dim> &t,
+               const Tensor<rank, dim> &t2,
+               const Tensor<rank, dim> &compare)
+{
+  deallog << "Constant rank " << rank << " and dim " << dim << '\t'
+          << schur_product(t, t2) << " compare " << compare << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    Tensor<0, 2> a, b, compare;
+    a       = 2;
+    b       = 2;
+    compare = 4;
+    test_same(a, compare);
+    test_different(a, b, compare);
+  }
+
+  {
+    Tensor<1, 2> a, b, compare;
+    for (size_t d = 0; d < 2; ++d)
+      {
+        a[d]       = 2;
+        b[d]       = 2;
+        compare[d] = 4;
+      }
+    test_same(a, compare);
+    test_different(a, b, compare);
+  }
+
+  {
+    Tensor<2, 2> a, b, compare;
+    for (size_t c = 0; c < 2; ++c)
+      for (size_t d = 0; d < 2; ++d)
+        {
+          a[c][d]       = 2;
+          b[c][d]       = 2;
+          compare[c][d] = 4;
+        }
+    test_same(a, compare);
+    test_different(a, b, compare);
+  }
+
+  {
+    Tensor<3, 2> a, b, compare;
+
+    for (size_t c = 0; c < 2; ++c)
+      for (size_t d = 0; d < 2; ++d)
+        for (size_t e = 0; e < 2; ++e)
+          {
+            a[c][d][e]       = 2;
+            b[c][d][e]       = 2;
+            compare[c][d][e] = 4;
+          }
+    test_same(a, compare);
+    test_different(a, b, compare);
+  }
+
+  return 0;
+}

diff --git a/tests/tensors/tensor_10.output b/tests/tensors/tensor_10.output
new file mode 100644 (file)
index 0000000..baa36bc
--- /dev/null
+++ b/tests/tensors/tensor_10.output
@@ -0,0 +1,9 @@
+
+DEAL::Constant rank 0 and dim 2        4.00000 compare 4.00000
+DEAL::Constant rank 0 and dim 2        4.00000 compare 4.00000
+DEAL::Constant rank 1 and dim 2        4.00000 4.00000 compare 4.00000 4.00000
+DEAL::Constant rank 1 and dim 2        4.00000 4.00000 compare 4.00000 4.00000
+DEAL::Constant rank 2 and dim 2        4.00000 4.00000 4.00000 4.00000 compare 4.00000 4.00000 4.00000 4.00000
+DEAL::Constant rank 2 and dim 2        4.00000 4.00000 4.00000 4.00000 compare 4.00000 4.00000 4.00000 4.00000
+DEAL::Constant rank 3 and dim 2        4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 compare 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000
+DEAL::Constant rank 3 and dim 2        4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 compare 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000 4.00000