Add a test for TensorAccessors::contract

diff --git a/tests/base/tensor_accessors_03.cc b/tests/base/tensor_accessors_03.cc
new file mode 100644 (file)
index 0000000..d878bcf
--- /dev/null
+++ b/tests/base/tensor_accessors_03.cc
@@ -0,0 +1,115 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test TensorAccessors::contract
+
+#include "../tests.h"
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/tensor_accessors.h>
+#include <deal.II/base/table_indices.h>
+
+
+int main()
+{
+  initlog();
+
+  int c_left1 [5] = {0, 1, 2, 3, 4};
+  Tensor<1, 5, int> left1 (c_left1);
+  int c_right1 [5] = {0, 10, 20, 30, 40};
+  Tensor<1, 5, int> right1 (c_right1);
+
+  deallog << "left1: " << left1 << std::endl;
+  deallog << "right1: " << right1 << std::endl;
+
+  // Apply contract with no_contr == 0, this is a plain outer product:
+  {
+    Tensor<2, 5, int> result;
+    TensorAccessors::contract<0, 1, 1, 5>(result, left1, right1);
+
+    deallog << "Outer product of left1 and right1:" << std::endl;
+    deallog << result << std::endl;
+  }
+
+  // Apply contract with no_contr == 1, this is a scalar product:
+  {
+    int result = 0;
+    TensorAccessors::contract<1, 1, 1, 5>(result, left1, right1);
+
+    deallog << "Scalar product of left1 and right1:" << std::endl;
+    deallog << result << std::endl;
+  }
+
+  Tensor<2, 5, int> left2;
+  left2[0] = left1;
+  left2[1] = 2 * left1;
+  left2[2] = 4 * left1;
+  left2[3] = 8 * left1;
+  left2[4] = 16 * left1;
+
+  Tensor<2, 5, int> right2;
+  right2[0] = right1;
+  right2[1] = 2 * right1;
+  right2[2] = 4 * right1;
+  right2[3] = 8 * right1;
+  right2[4] = 16 * right1;
+
+  deallog << "left2: " << left2 << std::endl;
+  deallog << "right2: " << right2 << std::endl;
+
+  // Apply contract with no_contr == 0, this is a plain outer product:
+  {
+    Tensor<4, 5, int> result;
+    TensorAccessors::contract<0, 2, 2, 5>(result, left2, right2);
+
+    deallog << "Outer product of left2 and right2:" << std::endl;
+    deallog << result << std::endl;
+    // The result is verified by hand:
+    // for (unsigned int i = 0; i < 4; ++i)
+    //   for (unsigned int j = 0; j < 4; ++j)
+    //     for (unsigned int k = 0; k < 4; ++k)
+    //       for (unsigned int l = 0; l < 4; ++l)
+    //         {
+    //           deallog << "(" << i << "," << j << "," << k << "," << l << "):  ";
+    //           deallog << result[i][j][k][l] << " = " << left2[i][j] << " * " << right2[k][l] << std::endl;
+    //         }
+  }
+
+  // Apply contract with no_contr == 1, and switch indices of right2. This
+  // corresponds to a contraction of the last index of left2 with the first
+  // index of right2:
+  {
+    dealii::TensorAccessors::internal::ReorderedIndexView<0, 2, dealii::Tensor<2, 5, int> >
+    reordered = TensorAccessors::reordered_index_view<0, 2>(right2);
+
+    Tensor<2, 5, int> result;
+    TensorAccessors::contract<1, 2, 2, 5>(result, left2, reordered);
+
+    deallog << "Contract the last index of left2 with the first index of right2:" << std::endl;
+    deallog << result << std::endl;
+    // Verified to be the same as the old implementation
+    // deallog << contract(left2, right2) << std::endl;
+  }
+
+  // Apply contract with no_contr == 2, this is a double contraction.
+  {
+    int result;
+    TensorAccessors::contract<2, 2, 2, 5>(result, left2, right2);
+
+    deallog << "Double contraction:" << std::endl;
+    deallog << result << std::endl;
+    // Verified to be the same as the old implementation
+    // deallog << double_contract(left2, right2) << std::endl;
+  }
+}

diff --git a/tests/base/tensor_accessors_03.output b/tests/base/tensor_accessors_03.output
new file mode 100644 (file)
index 0000000..8a8ac55
--- /dev/null
+++ b/tests/base/tensor_accessors_03.output
@@ -0,0 +1,15 @@
+
+DEAL::left1: 0 1 2 3 4
+DEAL::right1: 0 10 20 30 40
+DEAL::Outer product of left1 and right1:
+DEAL::0 0 0 0 0 0 10 20 30 40 0 20 40 60 80 0 30 60 90 120 0 40 80 120 160
+DEAL::Scalar product of left1 and right1:
+DEAL::300
+DEAL::left2: 0 1 2 3 4 0 2 4 6 8 0 4 8 12 16 0 8 16 24 32 0 16 32 48 64
+DEAL::right2: 0 10 20 30 40 0 20 40 60 80 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640
+DEAL::Outer product of left2 and right2:
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 20 30 40 0 20 40 60 80 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 20 40 60 80 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 30 60 90 120 0 60 120 180 240 0 120 240 360 480 0 240 480 720 960 0 480 960 1440 1920 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 40 60 80 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 60 120 180 240 0 120 240 360 480 0 240 480 720 960 0 480 960 1440 1920 0 960 1920 2880 3840 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 40 80 120 160 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 120 240 360 480 0 240 480 720 960 0 480 960 1440 1920 0 960 1920 2880 3840 0 1920 3840 5760 7680 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 80 160 240 320 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 240 480 720 960 0 480 960 1440 1920 0 960 1920 2880 3840 0 1920 3840 5760 7680 0 3840 7680 11520 15360 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 5120 10240 15360 20480 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 160 320 480 640 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 320 640 960 1280 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 5120 10240 15360 20480 0 480 960 1440 1920 0 960 1920 2880 3840 0 1920 3840 5760 7680 0 3840 7680 11520 15360 0 7680 15360 23040 30720 0 640 1280 1920 2560 0 1280 2560 3840 5120 0 2560 5120 7680 10240 0 5120 10240 15360 20480 0 10240 20480 30720 40960
+DEAL::Contract the last index of left2 with the first index of right2:
+DEAL::0 980 1960 2940 3920 0 1960 3920 5880 7840 0 3920 7840 11760 15680 0 7840 15680 23520 31360 0 15680 31360 47040 62720
+DEAL::Double contraction:
+DEAL::102300