]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Scott Miller (and testcases by me): Add scalar_product functions between...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Aug 2013 02:52:42 +0000 (02:52 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Aug 2013 02:52:42 +0000 (02:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@30315 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/base/symmetric_tensor.h
tests/base/symmetric_tensor_27.cc [new file with mode: 0644]
tests/base/symmetric_tensor_27/cmp/generic [new file with mode: 0644]

index 89bdcf89845ffd226d82f70daa79dc2eef431e5a..0b615ad489ff8e4ec7c4aba05dac8f6588074ca4 100644 (file)
@@ -56,6 +56,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  New: There are now global functions <code>scalar_product</code>
+  that compute the scalar product (double contraction) between
+  tensors of rank 2.
+  <br>
+  (Scott Miller, 2013/08/14)
+  </li>
+
   <li>
   Fixed: Creating objects of type MappingQ was previously only possible
   for low order polynomials. For orders higher than around 6, one ran
index 473792c92b63927e619ce4f29d9be59587822fec..93ffc6296a2c3e9b1f639e424da5ced3952fc312 100644 (file)
@@ -2842,8 +2842,68 @@ operator / (const SymmetricTensor<rank,dim> &t,
   return tt;
 }
 
+/**
+ * Compute the scalar product $a:b=\sum_{i,j} a_{ij}b_{ij}$ between two
+ * tensors $a,b$ of rank 2. In the current case where both arguments are
+ * symmetric tensors, this is equivalent to calling the expression
+ * <code>t1*t2</code> which uses the overloaded <code>operator*</code>
+ * between two symmetric tensors of rank 2.
+ *
+ * @relates SymmetricTensor
+ */
+template <int dim, typename Number>
+inline
+Number
+scalar_product (const SymmetricTensor<2,dim,Number> &t1,
+                const SymmetricTensor<2,dim,Number> &t2)
+{
+  return (t1*t2);
+}
+
+
+/**
+ * Compute the scalar product $a:b=\sum_{i,j} a_{ij}b_{ij}$ between two
+ * tensors $a,b$ of rank 2. We don't use <code>operator*</code> for this
+ * operation since the product between two tensors is usually assumed to be
+ * the contraction over the last index of the first tensor and the first index
+ * of the second tensor, for example $(a\cdot b)_{ij}=\sum_k a_{ik}b_{kj}$.
+ *
+ * @relates Tensor
+ * @relates SymmetricTensor
+ */
+template <int dim, typename Number>
+inline
+Number
+scalar_product (const SymmetricTensor<2,dim,Number> &t1,
+                const Tensor<2,dim,Number> &t2)
+{
+  Number s = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      s += t1[i][j] * t2[i][j];
+  return s;
+}
 
 
+/**
+ * Compute the scalar product $a:b=\sum_{i,j} a_{ij}b_{ij}$ between two
+ * tensors $a,b$ of rank 2. We don't use <code>operator*</code> for this
+ * operation since the product between two tensors is usually assumed to be
+ * the contraction over the last index of the first tensor and the first index
+ * of the second tensor, for example $(a\cdot b)_{ij}=\sum_k a_{ik}b_{kj}$.
+ *
+ * @relates Tensor
+ * @relates SymmetricTensor
+ */
+template <int dim, typename Number>
+inline
+Number
+scalar_product (const Tensor<2,dim,Number> &t1,
+                const SymmetricTensor<2,dim,Number> &t2)
+{
+  return scalar_product(t2, t1);
+}
+
 
 /**
  * Double contraction between a rank-4 and a rank-2 symmetric tensor,
diff --git a/tests/base/symmetric_tensor_27.cc b/tests/base/symmetric_tensor_27.cc
new file mode 100644 (file)
index 0000000..4d03de1
--- /dev/null
@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2005 - 2013 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 scalar_product between tensors and symmetric tensors
+
+#include "../tests.h"
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/logstream.h>
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+void test ()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  // initialize symmetric and non-symmetric tensors. in the former case, we
+  // only need to initialize one half
+  SymmetricTensor<2,dim> s;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=i; j<dim; ++j)
+      s[i][j] = (1.+(i+1)*(j*2));
+
+  Tensor<2,dim> t;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      t[i][j] = (1.+(i+1)*(j*2));
+
+  deallog << scalar_product(s,s) << std::endl;
+  Assert (scalar_product(s,s) == s*s, ExcInternalError());
+  
+  deallog << scalar_product(s,t) << std::endl;
+  Assert (scalar_product(s,t) == s*symmetrize(t), ExcInternalError());
+  Assert (scalar_product(s,t) == symmetrize(t)*s, ExcInternalError());
+  
+  deallog << scalar_product(t,s) << std::endl;
+  Assert (scalar_product(t,s) == s*symmetrize(t), ExcInternalError());
+  Assert (scalar_product(t,s) == symmetrize(t)*s, ExcInternalError());
+}
+
+
+
+
+int main ()
+{
+  std::ofstream logfile("symmetric_tensor_27/output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/base/symmetric_tensor_27/cmp/generic b/tests/base/symmetric_tensor_27/cmp/generic
new file mode 100644 (file)
index 0000000..48633ff
--- /dev/null
@@ -0,0 +1,14 @@
+
+DEAL::dim=1
+DEAL::1.00
+DEAL::1.00
+DEAL::1.00
+DEAL::dim=2
+DEAL::44.0
+DEAL::38.0
+DEAL::38.0
+DEAL::dim=3
+DEAL::425.
+DEAL::381.
+DEAL::381.
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.