From 21dc89dfbeb1b8944a72701f47548837cc823d8a Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Sun, 18 Aug 2019 21:01:51 -0600
Subject: [PATCH] Add a test.

---
 tests/tensors/tensor_30.cc     | 77 ++++++++++++++++++++++++++++++++++
 tests/tensors/tensor_30.output |  5 +++
 2 files changed, 82 insertions(+)
 create mode 100644 tests/tensors/tensor_30.cc
 create mode 100644 tests/tensors/tensor_30.output

diff --git a/tests/tensors/tensor_30.cc b/tests/tensors/tensor_30.cc
new file mode 100644
index 0000000000..72bfb85199
--- /dev/null
+++ b/tests/tensors/tensor_30.cc
@@ -0,0 +1,77 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Additional tests for tensor contraction
+//
+// A mysteriously failing pull request suggested that the contraction
+// of a rank-3 tensor and a rank-1 tensor might be using the wrong
+// indices. Check this.
+
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+  Tensor<3, dim> a;
+  Tensor<1, dim> b;
+
+  for (unsigned int i = 0; i < dim; ++i)
+    for (unsigned int j = 0; j < dim; ++j)
+      for (unsigned int k = 0; k < dim; ++k)
+        a[i][j][k] = i + 2 * j + 3 * k;
+
+  for (unsigned int k = 0; k < dim; ++k)
+    b[k] = k;
+
+  // Compute
+  //   c(i,j) = sum_k a(i,j,k)*b(k)
+  //          = sum_k (i + 2j + 3k)k
+  //          = (i + 2j) sum k  + 3 sum k^2
+  // which in 2d equates to
+  //          = (i + 2j) + 3
+  // and in 3d equates to
+  //          = 3(i + 2j) + 3*5
+  //
+  // In other words, in 2d this is the 2x2 matrix
+  //    [[3,4],[4,6]]
+  // and in 3d is the 3x3 matrix
+  //    [[15,21,27],[18,24,30],[21,27,33]]
+  Tensor<2, dim> c = a * b;
+
+  deallog << c << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/tensors/tensor_30.output b/tests/tensors/tensor_30.output
new file mode 100644
index 0000000000..ad198cbaf6
--- /dev/null
+++ b/tests/tensors/tensor_30.output
@@ -0,0 +1,5 @@
+
+DEAL::0.00000
+DEAL::3.00000 5.00000 4.00000 6.00000
+DEAL::15.0000 21.0000 27.0000 18.0000 24.0000 30.0000 21.0000 27.0000 33.0000
+DEAL::OK
-- 
2.39.5