Implement the operator* for mixed types also for SymmetricTensor.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 39e06367169af47bd0f56cb554cfcfc347182457..b449e56bd7e7343786a501ffc08e954bb346385c 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -355,7 +355,8 @@ inconvenience this causes.
   to store its elements as complex numbers, rather than as double
   variables. This is now changed: The result of the product of a Tensor
   and a scalar number is now a Tensor that stores its elements in a data
-  type appropriate for this product.
+  type appropriate for this product. The same approach is taken for the
+  SymmetricTensor class.
   <br>
   (Wolfgang Bangerth, 2015/02/11)
   </li>

diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h
index 63f4a025cca5cc489729ab4cd3b0e0583fced1ec..52b922de1351640053aa59cd51d674ccd72aeceb 100644 (file)
--- a/include/deal.II/base/symmetric_tensor.h
+++ b/include/deal.II/base/symmetric_tensor.h
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2005 - 2014 by the deal.II authors
+// Copyright (C) 2005 - 2015 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -518,7 +518,7 @@ public:
       n_independent_components;
 
   /**
-   * Default constructor. Creates a zero tensor.
+   * Default constructor. Creates a tensor with all entries equal to zero.
    */
   SymmetricTensor ();
 
@@ -551,6 +551,15 @@ public:
    */
   SymmetricTensor (const Number (&array) [n_independent_components]);
 
+  /**
+   * Copy constructor from tensors with different underlying scalar
+   * type. This obviously requires that the @p OtherNumber type is
+   * convertible to @p Number.
+   */
+  template <typename OtherNumber>
+  explicit
+  SymmetricTensor (const SymmetricTensor<rank,dim,OtherNumber> &initializer);
+
   /**
    * Assignment operator.
    */
@@ -955,6 +964,19 @@ SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Tensor<2,dim,Number> &t
 
 
 
+template <int rank, int dim, typename Number>
+template <typename OtherNumber>
+inline
+SymmetricTensor<rank,dim,Number>::
+SymmetricTensor (const SymmetricTensor<rank,dim,OtherNumber> &initializer)
+{
+  for (unsigned int i=0; i<n_independent_components; ++i)
+    data[i] = initializer.data[i];
+}
+
+
+
+
 template <int rank, int dim, typename Number>
 inline
 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
@@ -2644,8 +2666,10 @@ symmetrize (const Tensor<2,3,Number> &t)
 
 
 /**
- * Multiplication of a symmetric tensor of general rank with a scalar from the
- * right.
+ * Multiplication of a symmetric tensor of general rank with a scalar
+ * from the right. This version of the operator is used if the scalar
+ * has the same data type as is used to store the elements of the
+ * symmetric tensor.
  *
  * @relates SymmetricTensor
  */
@@ -2664,7 +2688,9 @@ operator * (const SymmetricTensor<rank,dim,Number> &t,
 
 /**
  * Multiplication of a symmetric tensor of general rank with a scalar from the
- * left.
+ * left. This version of the operator is used if the scalar
+ * has the same data type as is used to store the elements of the
+ * symmetric tensor.
  *
  * @relates SymmetricTensor
  */
@@ -2674,13 +2700,74 @@ SymmetricTensor<rank,dim,Number>
 operator * (const Number                            factor,
             const SymmetricTensor<rank,dim,Number> &t)
 {
-  SymmetricTensor<rank,dim,Number> tt = t;
-  tt *= factor;
+  // simply forward to the other operator
+  return t*factor;
+}
+
+
+/**
+ * Multiplication of a symmetric tensor with a scalar number from the right.
+ *
+ * The purpose of this operator is to enable only multiplication of a tensor
+ * by a scalar number (i.e., a floating point number, a complex floating point
+ * number, etc.). The function is written in a way that only allows the
+ * compiler to consider the function if the second argument is indeed a scalar
+ * number -- in other words, @p OtherNumber will not match, for example
+ * <code>std::vector@<double@></code> as the product of a tensor and a vector
+ * clearly would make no sense. The mechanism by which the compiler is
+ * prohibited of considering this operator for multiplication with non-scalar
+ * types are explained in the documentation of the EnableIfScalar class.
+ *
+ * The return type of the function is chosen so that it matches the types
+ * of both the tensor and the scalar argument. For example, if you multiply
+ * a <code>SymmetricTensor@<2,dim,double@></code> by <code>std::complex@<double@></code>,
+ * then the result will be a <code>SymmetricTensor@<2,dim,std::complex@<double@>@></code>.
+ * In other words, the type with which the returned tensor stores its
+ * components equals the type you would get if you multiplied an individual
+ * component of the input tensor by the scalar factor.
+ *
+ * @relates SymmetricTensor
+ * @relates EnableIfScalar
+ */
+template <int rank, int dim, typename Number, typename OtherNumber>
+inline
+SymmetricTensor<rank,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
+operator * (const SymmetricTensor<rank,dim,Number> &t,
+            const OtherNumber                    factor)
+{
+  // form the product. we have to convert the two factors into the final
+  // type via explicit casts because, for awkward reasons, the C++
+  // standard committee saw it fit to not define an
+  //   operator*(float,std::complex<double>)
+  // (as well as with switched arguments and double<->float).
+  typedef typename ProductType<Number,OtherNumber>::type product_type;
+  SymmetricTensor<rank,dim,product_type> tt(t);
+  tt *= product_type(factor);
   return tt;
 }
 
 
 
+/**
+ * Multiplication of a symmetric tensor with a scalar number from the left.
+ * See the discussion with the operator with switched arguments for more
+ * information about template arguments and the return type.
+ *
+ * @relates SymmetricTensor
+ * @relates EnableIfScalar
+ */
+template <int rank, int dim, typename Number, typename OtherNumber>
+inline
+SymmetricTensor<rank,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
+operator * (const Number                     factor,
+            const SymmetricTensor<rank,dim,OtherNumber> &t)
+{
+  // simply forward to the other operator with switched arguments
+  return (t*factor);
+}
+
+
+
 /**
  * Division of a symmetric tensor of general rank by a scalar.
  *

diff --git a/tests/base/product_type_01-symmetric.cc b/tests/base/product_type_01-symmetric.cc
new file mode 100644 (file)
index 0000000..c5ee98a
--- /dev/null
+++ b/tests/base/product_type_01-symmetric.cc
@@ -0,0 +1,62 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+// like _01 but for SymmetricTensor
+
+#include "../tests.h"
+#include <iomanip>
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+#include <typeinfo>
+#include <complex>
+
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/symmetric_tensor.h>
+
+
+template <typename T, typename U, typename CompareType>
+void check()
+{
+  Assert (typeid(typename ProductType<T,U>::type) == typeid(CompareType),
+         ExcInternalError());
+  Assert (typeid(typename ProductType<T,U>::type) == typeid(T() * U()),
+         ExcInternalError());
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  // check product with SymmetricTensor<2,dim>
+  check<SymmetricTensor<2,2,double>,double,SymmetricTensor<2,2,double> >();
+  check<SymmetricTensor<2,2,float>,double,SymmetricTensor<2,2,double> >();
+  check<double,SymmetricTensor<2,2,float>,SymmetricTensor<2,2,double> >();
+
+  // check product with std::complex. rather annoyingly, there is no
+  // product between std::complex<double> and float, or the other way
+  // around, so stay within the same type system
+  check<std::complex<double>,double,std::complex<double> >();
+  check<std::complex<float>,float,std::complex<float> >();
+  check<SymmetricTensor<2,2>,std::complex<double>,SymmetricTensor<2,2,std::complex<double> > >();
+  check<std::complex<double>,SymmetricTensor<2,2>,SymmetricTensor<2,2,std::complex<double> > >();
+
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/base/product_type_01-symmetric.output b/tests/base/product_type_01-symmetric.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
+++ b/tests/base/product_type_01-symmetric.output
@@ -0,0 +1,2 @@
+
+DEAL::OK

diff --git a/tests/base/product_type_03-symmetric.cc b/tests/base/product_type_03-symmetric.cc
new file mode 100644 (file)
index 0000000..a7a1dad
--- /dev/null
+++ b/tests/base/product_type_03-symmetric.cc
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+// like _03 but for SymmetricTensor
+
+#include "../tests.h"
+#include <iomanip>
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+#include <typeinfo>
+#include <complex>
+
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/symmetric_tensor.h>
+
+
+template <typename T, typename U, typename CompareType>
+void check()
+{
+  Assert (typeid(T() * U()) == typeid(CompareType),
+         ExcInternalError());
+  Assert (typeid(T() * U()) == typeid(CompareType),
+         ExcInternalError());
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  // check product of scalars
+  check<SymmetricTensor<2,1,double>,double,SymmetricTensor<2,1,double> >();
+  check<double,SymmetricTensor<2,1,double>,SymmetricTensor<2,1,double> >();
+
+  check<SymmetricTensor<2,1,double>,float,SymmetricTensor<2,1,double> >();
+  check<float,SymmetricTensor<2,1,double>,SymmetricTensor<2,1,double> >();
+
+  check<SymmetricTensor<2,1,double>,std::complex<double>,SymmetricTensor<2,1,std::complex<double> > >();
+  check<std::complex<double>,SymmetricTensor<2,1,double>,SymmetricTensor<2,1,std::complex<double> > >();
+
+  check<SymmetricTensor<2,1,double>,std::complex<float>,SymmetricTensor<2,1,std::complex<double> > >();
+  check<std::complex<float>,SymmetricTensor<2,1,double>,SymmetricTensor<2,1,std::complex<double> > >();
+
+  check<SymmetricTensor<2,1,float>,std::complex<double>,SymmetricTensor<2,1,std::complex<double> > >();
+  check<std::complex<double>,SymmetricTensor<2,1,float>,SymmetricTensor<2,1,std::complex<double> > >();
+
+  check<SymmetricTensor<2,1,float>,std::complex<float>,SymmetricTensor<2,1,std::complex<float> > >();
+  check<std::complex<float>,SymmetricTensor<2,1,float>,SymmetricTensor<2,1,std::complex<float> > >();
+  
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/base/product_type_03-symmetric.output b/tests/base/product_type_03-symmetric.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
+++ b/tests/base/product_type_03-symmetric.output
@@ -0,0 +1,2 @@
+
+DEAL::OK