From 978b8dcb770c2845874335d6bf1611f692dff2b4 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 11 Sep 2015 21:36:59 -0500 Subject: [PATCH] Tensor - Implement operator* with TensorAccessors::contract --- include/deal.II/base/tensor.h | 473 +++++++++++----------------------- 1 file changed, 157 insertions(+), 316 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 7d231bc49d..b86ea17fb9 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -244,13 +244,13 @@ public: template void serialize(Archive &ar, const unsigned int version); -private: /** * Internal type declaration that is used to specialize the return type * of operator[]() for Tensor<1,dim,Number> */ typedef Number tensor_type; +private: /** * The value of this scalar object. */ @@ -557,13 +557,13 @@ public: << arg1 << ", but this is not possible for tensors of the current type."); -private: /** * Internal type declaration that is used to specialize the return type * of operator[]() for Tensor<1,dim,Number> */ typedef Tensor tensor_type; +private: /** * Array of tensors holding the subelements. */ @@ -1164,21 +1164,6 @@ Tensor::serialize(Archive &ar, const unsigned int) /* ----------------- Non-member functions operating on tensors. ------------ */ - -#ifndef DEAL_II_WITH_CXX11 -template -struct ProductType > -{ - typedef Tensor::type> type; -}; - -template -struct ProductType,U> -{ - typedef Tensor::type> type; -}; -#endif - /** * @name Output functions for Tensor objects */ @@ -1227,35 +1212,80 @@ std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p) */ //@{ + +#ifndef DEAL_II_WITH_CXX11 +template +struct ProductType > +{ + typedef Tensor::type> type; +}; + +template +struct ProductType,U> +{ + typedef Tensor::type> type; +}; +#endif + + + /** - * Scalar multiplication of a tensor of rank 0 with a scalar from the left. + * Scalar multiplication of a tensor of rank 0 with an object from the + * left. + * + * This function unwraps the underlying @p Number stored in the Tensor and + * multiplies @p object with it. * * @relates Tensor<0,dim,Number> * @relates EnableIfScalar */ -template +template inline -Tensor<0,dim,typename ProductType::type, Number>::type> -operator * (const OtherNumber factor, +typename ProductType::type +operator * (const Other object, const Tensor<0,dim,Number> &t) { - return factor * static_cast(t); + return object * static_cast(t); } /** - * Scalar multiplication of a tensor of rank 0 with a scalar from the right. + * Scalar multiplication of a tensor of rank 0 with an object from the + * right. + * + * This function unwraps the underlying @p Number stored in the Tensor and + * multiplies @p object with it. * * @relates Tensor<0,dim,Number> * @relates EnableIfScalar */ -template +template inline -Tensor<0,dim,typename ProductType::type>::type> +typename ProductType::type operator * (const Tensor<0,dim,Number> &t, - const OtherNumber factor) + const Other object) { - return static_cast(t) * factor; + return static_cast(t) * object; +} + + +/** + * Scalar multiplication of two tensors of rank 0. + * + * This function unwraps the underlying objects of type @p Number and @p + * OtherNumber that are stored within the Tensor and multiplies them. + * It returns an unwrapped number of product type. + * + * @relates Tensor<0,dim,Number> + */ +template +inline +typename ProductType::type // FIXME: TEST! +operator * (const Tensor<0, dim, Number> &src1, + const Tensor<0, dim, OtherNumber> &src2) +{ + return static_cast(src1) * + static_cast(src2); } @@ -1307,24 +1337,9 @@ operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) * Multiplication of a tensor of general rank 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 - * std::vector@ 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 - * Tensor@<1,dim,double@> by std::complex@, - * then the result will be a - * Tensor@<1,dim,std::complex@@>. 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. + * Only multiplication with a scalar number type (i.e., a floating point + * number, a complex floating point number, etc.), see the documentation of + * EnableIfScalar for details. * * @relates Tensor * @relates EnableIfScalar @@ -1450,20 +1465,54 @@ operator- (const Tensor &p, const Tensor */ //@{ + /** - * Returns the contraction of two Tensors of rank 0. + * The dot product (single contraction) for tensors: Return a tensor of + * rank $(\text{rank\_1} + \text{rank\_2} - 2)$ that is the contraction of + * the last index of a tensor @p src1 of rank @p rank_1 with the first + * index of a tensor @p src2 of rank @p rank_2: + * @f[ + * \text{result}_{i_1,..,i_{r1},j_1,..,j_{r2}} + * = \sum_{k} + * \text{left}_{i_1,..,i_{r1}, k} + * \text{right}_{j_1,..,j_{r2}, k} + * @f] + * + * @note For the Tensor class, the multiplication operator only + * performs a contraction over a single pair of indices. This is in + * contrast to the multiplication operator for symmetric tensors, which + * does the double contraction. + * + * @note In case the contraction yields tensor of rank 0 the scalar + * number is returned as an unwrapped number type * - * @relates Tensor<0,dim,Number> + * @relates Tensor */ -template +template ::type, + typename = typename std::enable_if::type> + inline -typename ProductType::type -operator* (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) +typename Tensor::type>::tensor_type +operator * (const Tensor &src1, + const Tensor &src2) { - return static_cast(p) * static_cast(q); + typename Tensor::type>::tensor_type result; + + TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor > + reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2); + TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered); + + return result; } + //@} +/** + * @name To be refactored + */ +//@{ /** @@ -1487,30 +1536,6 @@ contract (const Tensor<1,dim,Number> &src1, } -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - */ -template -inline -typename ProductType::type -operator * (const Tensor<1,dim,Number> &src1, - const Tensor<1,dim,OtherNumber> &src2) -{ - return contract(src1, src2); -} - - /** * Double contract two tensors of rank 2, thus computing the Frobenius inner * product sumi,j src1[i][j]*src2[i][j]. @@ -1553,37 +1578,6 @@ void contract (Tensor<1,dim,Number> &dest, } -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - * @author Wolfgang Bangerth, 2005 - */ -template -Tensor<1,dim,Number> -operator * (const Tensor<2,dim,Number> &src1, - const Tensor<1,dim,Number> &src2) -{ - Tensor<1,dim,Number> dest; - for (unsigned int i=0; idest[i] = sum_j src1[j] src2[j][i]. @@ -1606,38 +1600,6 @@ void contract (Tensor<1,dim,Number> &dest, } -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - * @author Wolfgang Bangerth, 2005 - */ -template -inline -Tensor<1,dim,Number> -operator * (const Tensor<1,dim,Number> &src1, - const Tensor<2,dim,Number> &src2) -{ - Tensor<1,dim,Number> dest; - for (unsigned int i=0; idest[i][k] = sum_j src1[i][j] src2[j][k]. @@ -1661,37 +1623,6 @@ void contract (Tensor<2,dim,Number> &dest, } - -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - * @author Wolfgang Bangerth, 2005 - */ -template -inline -Tensor<2,dim,Number> -operator * (const Tensor<2,dim,Number> &src1, - const Tensor<2,dim,Number> &src2) -{ - Tensor<2,dim,Number> dest; - for (unsigned int i=0; iindex1 of the first tensor, and @@ -1934,37 +1865,6 @@ void contract (Tensor<3,dim,Number> &dest, } -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - * @author Wolfgang Bangerth, 2005 - */ -template -inline -Tensor<3,dim,Number> -operator * (const Tensor<3,dim,Number> &src1, - const Tensor<2,dim,Number> &src2) -{ - Tensor<3,dim,Number> dest; - for (unsigned int i=0; idest[i][j][l] = sum_k src1[i][k] src2[k][j][l]. @@ -1987,61 +1887,6 @@ void contract (Tensor<3,dim,Number> &dest, } -/** - * Multiplication operator performing a contraction of the last index of the - * first argument and the first index of the second argument. This function - * therefore does the same as the corresponding contract function, - * but returns the result as a return value, rather than writing it into the - * reference given as the first argument to the contract function. - * - * Note that for the Tensor class, the multiplication operator only - * performs a contraction over a single pair of indices. This is in contrast - * to the multiplication operator for symmetric tensors, which does the double - * contraction. - * - * @relates Tensor - * @author Wolfgang Bangerth, 2005 - */ -template -inline -Tensor<3,dim,Number> -operator * (const Tensor<2,dim,Number> &src1, - const Tensor<3,dim,Number> &src2) -{ - Tensor<3,dim,Number> dest; - for (unsigned int i=0; idest[i][j][k][l] = sum_m src1[i][j][m] src2[m][k][l]. - * - * @relates Tensor - * @author Wolfgang Bangerth, 1998 - */ -template -inline -Tensor<4,dim,Number> -operator * (const Tensor<3,dim,Number> &src1, - const Tensor<3,dim,Number> &src2) -{ - Tensor<4,dim,Number> dest; - for (unsigned int i=0; isrc1 with the two indices * src2, creating a rank-2 tensor. This is the matrix-vector product @@ -2233,7 +2078,6 @@ void outer_product (Tensor<1,dim,Number> &dst, } - /** * Form the outer product of two tensors of rank 1 and 0, i.e. dst[i] = * src1[i] * src2. Of course, this is only a scaling of src1, @@ -2254,55 +2098,6 @@ void outer_product (Tensor<1,dim,Number> &dst, } -/** - * Cross-product in 2d. This is just a rotation by 90 degrees clockwise to - * compute the outer normal from a tangential vector. This function is defined - * for all space dimensions to allow for dimension independent programming - * (e.g. within switches over the space dimension), but may only be called if - * the actual dimension of the arguments is two (e.g. from the dim==2 - * case in the switch). - * - * @relates Tensor - * @author Guido Kanschat, 2001 - */ -template -inline -void -cross_product (Tensor<1,dim,Number> &dst, - const Tensor<1,dim,Number> &src) -{ - Assert (dim==2, ExcInternalError()); - - dst[0] = src[1]; - dst[1] = -src[0]; -} - - -/** - * Cross-product of 2 vectors in 3d. This function is defined for all space - * dimensions to allow for dimension independent programming (e.g. within - * switches over the space dimension), but may only be called if the actual - * dimension of the arguments is three (e.g. from the dim==3 case in - * the switch). - * - * @relates Tensor - * @author Guido Kanschat, 2001 - */ -template -inline -void -cross_product (Tensor<1,dim,Number> &dst, - const Tensor<1,dim,Number> &src1, - const Tensor<1,dim,Number> &src2) -{ - Assert (dim==3, ExcInternalError()); - - dst[0] = src1[1]*src2[2] - src1[2]*src2[1]; - dst[1] = src1[2]*src2[0] - src1[0]*src2[2]; - dst[2] = src1[0]*src2[1] - src1[1]*src2[0]; -} - - /** * 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 operator* for this @@ -2342,7 +2137,6 @@ Number determinant (const Tensor &t) } - /** * Compute the determinant of a tensor of rank one and dimension one. Since * this is a number, the return value is, of course, the number itself. @@ -2358,6 +2152,62 @@ Number determinant (const Tensor<1,1,Number> &t) } +//@} +/** + * @name Special operations on tensors of rank 1 + */ +//@{ + + +/** + * Cross-product in 2d. This is just a rotation by 90 degrees clockwise to + * compute the outer normal from a tangential vector. This function is defined + * for all space dimensions to allow for dimension independent programming + * (e.g. within switches over the space dimension), but may only be called if + * the actual dimension of the arguments is two (e.g. from the dim==2 + * case in the switch). + * + * @relates Tensor + * @author Guido Kanschat, 2001 + */ +template +inline +void +cross_product (Tensor<1,dim,Number> &dst, + const Tensor<1,dim,Number> &src) +{ + Assert (dim==2, ExcInternalError()); + + dst[0] = src[1]; + dst[1] = -src[0]; +} + + +/** + * Cross-product of 2 vectors in 3d. This function is defined for all space + * dimensions to allow for dimension independent programming (e.g. within + * switches over the space dimension), but may only be called if the actual + * dimension of the arguments is three (e.g. from the dim==3 case in + * the switch). + * + * @relates Tensor + * @author Guido Kanschat, 2001 + */ +template +inline +void +cross_product (Tensor<1,dim,Number> &dst, + const Tensor<1,dim,Number> &src1, + const Tensor<1,dim,Number> &src2) +{ + Assert (dim==3, ExcInternalError()); + + dst[0] = src1[1]*src2[2] - src1[2]*src2[1]; + dst[1] = src1[2]*src2[0] - src1[0]*src2[2]; + dst[2] = src1[0]*src2[1] - src1[1]*src2[0]; +} + + /** * Compute the determinant of a tensor of rank two and dimension one. Since * this is a number, the return value is, of course, the number itself. @@ -2373,7 +2223,6 @@ Number determinant (const Tensor<2,1,Number> &t) } - /** * Compute the determinant of a tensor or rank 2, here for dim==2. * @@ -2453,7 +2302,6 @@ Number determinant (const Tensor<2,dim,Number> &t) } - /** * Compute and return the trace of a tensor of rank 2, i.e. the sum of its * diagonal entries. @@ -2471,7 +2319,6 @@ Number trace (const Tensor<2,dim,Number> &d) } - /** * Compute and return the inverse of the given tensor. Since the compiler can * perform the return value optimization, and since the size of the return @@ -2540,7 +2387,6 @@ invert (const Tensor<2,dim,Number> &t) } - /** * Return the transpose of the given tensor. Since the compiler can perform * the return value optimization, and since the size of the return object is @@ -2587,8 +2433,6 @@ transpose (const Tensor<2,1,Number> &t) } - - /** * Return the transpose of the given tensor. This is the specialization of the * general template for dim==2. @@ -2606,8 +2450,6 @@ transpose (const Tensor<2,2,Number> &t) } - - /** * Return the transpose of the given tensor. This is the specialization of the * general template for dim==3. @@ -2657,7 +2499,6 @@ l1_norm (const Tensor<2,dim,Number> &t) } - /** * Return the $l_\infty$ norm of the given rank-2 tensor, where $||t||_\infty * = \max_i \sum_j |t_{ij}|$ (maximum of the sums over rows). @@ -2684,7 +2525,7 @@ linfty_norm (const Tensor<2,dim,Number> &t) return max; } - +//@} DEAL_II_NAMESPACE_CLOSE -- 2.39.5