<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
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,
--- /dev/null
+// ---------------------------------------------------------------------
+// $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;
+}
--- /dev/null
+
+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