From 362d243fa57816ff2bb714bde85e6284b6eeee57 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 29 Mar 2005 23:32:01 +0000 Subject: [PATCH] Remove global operator* and replace it by a member operator that is at least 2 orders of magnitude faster. git-svn-id: https://svn.dealii.org/trunk@10311 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/symmetric_tensor.h | 410 +++++++++---------- 1 file changed, 193 insertions(+), 217 deletions(-) diff --git a/deal.II/base/include/base/symmetric_tensor.h b/deal.II/base/include/base/symmetric_tensor.h index 49e90a7320..b942334ebe 100644 --- a/deal.II/base/include/base/symmetric_tensor.h +++ b/deal.II/base/include/base/symmetric_tensor.h @@ -91,7 +91,50 @@ namespace internal } Assert (false, ExcInternalError()); return TableIndices<4>(); - } + } + + + /** + * Typedef template magic + * denoting the result of a + * double contraction between two + * tensors or ranks rank1 and + * rank2. In general, this is a + * tensor of rank + * rank1+rank2-4, but if + * this is zero it is a single + * scalar double. For this case, + * we have a specialization. + * + * @author Wolfgang Bangerth, 2005 + */ + template + struct double_contraction_result + { + typedef SymmetricTensor type; + }; + + + /** + * Typedef template magic + * denoting the result of a + * double contraction between two + * tensors or ranks rank1 and + * rank2. In general, this is a + * tensor of rank + * rank1+rank2-4, but if + * this is zero it is a single + * scalar double. For this case, + * we have a specialization. + * + * @author Wolfgang Bangerth, 2005 + */ + template + struct double_contraction_result<2,2,dim> + { + typedef double type; + }; + /** @@ -628,16 +671,37 @@ class SymmetricTensor SymmetricTensor operator - () const; /** - * Scalar product between two symmetric - * tensors. It is the contraction - * aijbij - * over all indices i,j. While - * it is possible to define other scalar - * products (and associated induced - * norms), this one seems to be the most - * appropriate one. + * Product between the present + * symmetric tensor and a tensor + * of rank 2. For example, if the + * present object is also a + * rank-2 tensor, then this is + * the scalar-product contraction + * aijbij + * over all indices + * i,j. In this case, + * the return value evaluates to + * a single scalar. While it is + * possible to define other + * scalar product (and associated + * induced norms), this one seems + * to be the most appropriate + * one. + * + * If the present object is a + * rank-4 tensor, the the result + * is a rank-2 tensor, the + * operation contracts over the + * last two indices of the + * present object and the indices + * of the argument, and the + * result is a tensor of rank 2. */ - double operator * (const SymmetricTensor &s) const; + typename internal::SymmetricTensorAccessors::double_contraction_result::type + operator * (const SymmetricTensor<2,dim> &s) const; + + typename internal::SymmetricTensorAccessors::double_contraction_result::type + operator * (const SymmetricTensor<4,dim> &s) const; /** * Return a read-write reference @@ -719,6 +783,11 @@ class SymmetricTensor * Data storage for a symmetric tensor. */ typename internal::SymmetricTensorAccessors::StorageType::base_tensor_type data; + + /** + * Make all other symmetric tensors friends. + */ + template friend class SymmetricTensor; }; @@ -944,7 +1013,7 @@ SymmetricTensor::memory_consumption () template <> -double +internal::SymmetricTensorAccessors::double_contraction_result<2,2,1>::type SymmetricTensor<2,1>::operator * (const SymmetricTensor<2,1> &s) const { return data[0] * s.data[0]; @@ -953,7 +1022,7 @@ SymmetricTensor<2,1>::operator * (const SymmetricTensor<2,1> &s) const template <> -double +internal::SymmetricTensorAccessors::double_contraction_result<2,2,2>::type SymmetricTensor<2,2>::operator * (const SymmetricTensor<2,2> &s) const { return (data[0] * s.data[0] + @@ -964,7 +1033,7 @@ SymmetricTensor<2,2>::operator * (const SymmetricTensor<2,2> &s) const template <> -double +internal::SymmetricTensorAccessors::double_contraction_result<2,2,3>::type SymmetricTensor<2,3>::operator * (const SymmetricTensor<2,3> &s) const { return (data[0] * s.data[0] + @@ -977,55 +1046,105 @@ SymmetricTensor<2,3>::operator * (const SymmetricTensor<2,3> &s) const +internal::SymmetricTensorAccessors::double_contraction_result<4,2,1>::type +SymmetricTensor<4,1>:: +operator * (const SymmetricTensor<2,1> &s) const +{ + const unsigned int dim = 1; + SymmetricTensor<2,dim> tmp; + tmp.data[0] = data[0][0] * s.data[0]; + return tmp; +} + + + template <> -double -SymmetricTensor<4,1>::operator * (const SymmetricTensor<4,1> &s) const +internal::SymmetricTensorAccessors::double_contraction_result<2,4,1>::type +SymmetricTensor<2,1>:: +operator * (const SymmetricTensor<4,1> &s) const { - return data[0][0] * s.data[0][0]; + const unsigned int dim = 1; + SymmetricTensor<2,dim> tmp; + tmp[0][0] = data[0] * s.data[0][0]; + return tmp; } template <> -double -SymmetricTensor<4,2>::operator * (const SymmetricTensor<4,2> &s) const +internal::SymmetricTensorAccessors::double_contraction_result<4,2,2>::type +SymmetricTensor<4,2>:: +operator * (const SymmetricTensor<2,2> &s) const { const unsigned int dim = 2; + SymmetricTensor<2,dim> tmp; - // this is not really efficient and - // could be improved by counting - // how often each tensor entry is - // accessed, but this isn't a - // really frequent operation anyway - double t = 0; - for (unsigned int i=0; i -double -SymmetricTensor<4,3>::operator * (const SymmetricTensor<4,3> &s) const +internal::SymmetricTensorAccessors::double_contraction_result<2,4,2>::type +SymmetricTensor<2,2>:: +operator * (const SymmetricTensor<4,2> &s) const +{ + const unsigned int dim = 2; + SymmetricTensor<2,dim> tmp; + + for (unsigned int i=0; i<3; ++i) + tmp.data[i] = data[0] * s.data[0][i] + + data[1] * s.data[1][i] + + 2 * data[2] * s.data[2][i]; + + return tmp; +} + + + +template <> +internal::SymmetricTensorAccessors::double_contraction_result<4,2,3>::type +SymmetricTensor<4,3>:: +operator * (const SymmetricTensor<2,3> &s) const { const unsigned int dim = 3; + SymmetricTensor<2,dim> tmp; - // this is not really efficient and - // could be improved by counting - // how often each tensor entry is - // accessed, but this isn't a - // really frequent operation anyway - double t = 0; - for (unsigned int i=0; i +internal::SymmetricTensorAccessors::double_contraction_result<2,4,3>::type +SymmetricTensor<2,3>:: +operator * (const SymmetricTensor<4,3> &s) const +{ + const unsigned int dim = 3; + SymmetricTensor<2,dim> tmp; + + for (unsigned int i=0; i<6; ++i) + tmp.data[i] = data[0] * s.data[0][i] + + data[1] * s.data[1][i] + + data[2] * s.data[2][i] + + 2 * data[3] * s.data[3][i] + + 2 * data[4] * s.data[4][i] + + 2 * data[5] * s.data[5][i]; + + return tmp; } @@ -1651,155 +1770,6 @@ operator / (const SymmetricTensor &t, } -/** - * Double contraction between a rank-4 and a rank-2 symmetric tensor, - * resulting in a symmetric tensor of rank 2. This operation is the - * symmetric tensor analogon of a matrix-vector multiplication. - * - * @related SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -SymmetricTensor<2,1> -operator * (const SymmetricTensor<4,1> &t, - const SymmetricTensor<2,1> &s) -{ - const unsigned int dim = 1; - SymmetricTensor<2,dim> tmp; - tmp[0][0] = t[0][0][0][0] * s[0][0]; - return tmp; -} - - - -/** - * Double contraction between a rank-4 and a rank-2 symmetric tensor, - * resulting in a symmetric tensor of rank 2. This operation is the - * symmetric tensor analogon of a matrix-vector multiplication. - * - * @related SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -SymmetricTensor<2,1> -operator * (const SymmetricTensor<2,1> &s, - const SymmetricTensor<4,1> &t) -{ - const unsigned int dim = 1; - SymmetricTensor<2,dim> tmp; - tmp[0][0] = t[0][0][0][0] * s[0][0]; - return tmp; -} - - - -/** - * Double contraction between a rank-4 and a rank-2 symmetric tensor, - * resulting in a symmetric tensor of rank 2. This operation is the - * symmetric tensor analogon of a matrix-vector multiplication. - * - * @related SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -SymmetricTensor<2,2> -operator * (const SymmetricTensor<4,2> &t, - const SymmetricTensor<2,2> &s) -{ - const unsigned int dim = 2; - SymmetricTensor<2,dim> tmp; - - for (unsigned int i=0; i -operator * (const SymmetricTensor<2,2> &s, - const SymmetricTensor<4,2> &t) -{ - const unsigned int dim = 2; - SymmetricTensor<2,dim> tmp; - - for (unsigned int i=0; i -operator * (const SymmetricTensor<4,3> &t, - const SymmetricTensor<2,3> &s) -{ - const unsigned int dim = 3; - SymmetricTensor<2,dim> tmp; - - for (unsigned int i=0; i -operator * (const SymmetricTensor<2,3> &s, - const SymmetricTensor<4,3> &t) -{ - const unsigned int dim = 3; - SymmetricTensor<2,dim> tmp; - - for (unsigned int i=0; i &s, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 @@ -1832,10 +1803,11 @@ double_contract (SymmetricTensor<2,1> &tmp, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 @@ -1856,10 +1828,11 @@ double_contract (SymmetricTensor<2,1> &tmp, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor @author Wolfgang Bangerth, 2005 */ @@ -1885,10 +1858,11 @@ double_contract (SymmetricTensor<2,2> &tmp, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 @@ -1915,10 +1889,11 @@ double_contract (SymmetricTensor<2,2> &tmp, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 @@ -1948,10 +1923,11 @@ double_contract (SymmetricTensor<2,3> &tmp, * argument to this function. This operation is the symmetric tensor * analogon of a matrix-vector multiplication. * - * This function does the same as a respective operator*, but it avoid - * a temporary object (that the compiler can optimize away in many - * cases, however). This function mostly exists for compatibility - * purposes with the general tensor class. + * This function does the same as the member operator* of the + * SymmetricTensor class. It should not be used, however, since the + * member operator has knowledge of the actual data storage format and + * is at least 2 orders of magnitude faster. This function mostly + * exists for compatibility purposes with the general tensor class. * * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 -- 2.39.5