From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Fri, 13 Feb 2015 16:07:13 +0000 (-0600) Subject: Implement the operator* for mixed types also for SymmetricTensor. X-Git-Tag: v8.3.0-rc1~465^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d6008a8252fd558cbfe548816faad1b34210f025;p=dealii.git Implement the operator* for mixed types also for SymmetricTensor. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 39e0636716..b449e56bd7 100644 --- 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 63f4a025cc..52b922de13 100644 --- 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 index 0000000000..c5ee98a6a7 --- /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 index 0000000000..0fd8fc12f0 --- /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 index 0000000000..a7a1dada2c --- /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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/product_type_03-symmetric.output @@ -0,0 +1,2 @@ + +DEAL::OK