From 73100210d54e4c34fbd0ca96d3018ea43219a8d5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 29 Mar 2005 15:46:41 +0000 Subject: [PATCH] Move closer to general symmetric tensors. git-svn-id: https://svn.dealii.org/trunk@10288 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/symmetric_tensor.h | 587 ++++++++++++------- 1 file changed, 380 insertions(+), 207 deletions(-) diff --git a/deal.II/base/include/base/symmetric_tensor.h b/deal.II/base/include/base/symmetric_tensor.h index c52b1bd999..d8a3b40889 100644 --- a/deal.II/base/include/base/symmetric_tensor.h +++ b/deal.II/base/include/base/symmetric_tensor.h @@ -15,10 +15,10 @@ #include +#include template class SymmetricTensor; -template class SymmetricTensor<2,dim>; namespace internal @@ -122,9 +122,7 @@ namespace internal template struct AccessorTypes { - typedef - const typename StorageType::base_tensor_type - base_tensor_type; + typedef const ::SymmetricTensor tensor_type; typedef double reference; }; @@ -141,86 +139,95 @@ namespace internal template struct AccessorTypes { - typedef - typename StorageType::base_tensor_type - base_tensor_type; + typedef ::SymmetricTensor tensor_type; typedef double &reference; }; + + template + class Accessor; - namespace Rank2Accessors + /** + * Accessor class to access the elements + * of individual rows in a symmetric + * tensor of rank 2. Since the elements + * of symmetric tensors are not stored as + * in a table, the accessors are a little + * more involved. However, for tensors of + * rank 2 they are still relatively + * simple in that an accessor is created + * by the SymmetricTensor class with the + * first access to operator[]; + * the accessor thereby points to a row + * of the tensor. Calling + * operator[] on the accessor + * then selects an entry of this + * row. Note that if this entry is not + * actually stored, then the transpose + * entry is chosen as that is guaranteed + * to be stored. + * + * @author Wolfgang Bangerth, 2005 + */ + template + class Accessor<2,dim,constness> { + public: + /** + * Import which tensor we work on. + */ + typedef + typename AccessorTypes<2,dim,constness>::tensor_type + tensor_type; + + /** + * The type of a reference to an + * individual element of the + * symmetric tensor. If the tensor + * is constant, we can only return + * a value instead of a reference. + */ + typedef typename AccessorTypes<2,dim,constness>::reference reference; + + /** + * Constructor. Take the tensor to + * access as well as the row we + * point to as arguments. + */ + Accessor (tensor_type &tensor, + const unsigned int row); - /** - * Accessor class to access the - * elements of individual rows in a - * symmetric tensor. Since the elements - * of symmetric tensors are not stored - * as in a table, the accessors are a - * little more involved. - * - * @author Wolfgang Bangerth, 2005 - */ - template - class Accessor - { - public: - /** - * Import which tensor we work on. - */ - typedef - typename AccessorTypes<2,dim,constness>::base_tensor_type - base_tensor_type; - - /** - * The type of a reference to an - * individual element of the - * symmetric tensor. If the tensor - * is constant, we can only return - * a value instead of a reference. - */ - typedef typename AccessorTypes<2,dim,constness>::reference reference; - - /** - * Constructor. Take the tensor to - * access as well as the row we - * point to as arguments. - */ - Accessor (base_tensor_type &tensor, - const unsigned int row); - - /** - * Return a reference to an element - * of this row (if we point to a - * non-const tensor), or the value - * of the element (in case this is - * a constant tensor). - */ - reference operator[] (const unsigned int column); + /** + * Return a reference to an element + * of this row (if we point to a + * non-const tensor), or the value + * of the element (in case this is + * a constant tensor). + */ + reference operator[] (const unsigned int column); - private: - /** - * Reference to the tensor we - * access. - */ - base_tensor_type &base_tensor; - - /** - * Index of the row we access. - */ - const unsigned int row; - - /** - * Make the symmetric tensor - * classes a friend, since they are - * the only ones who can create - * objects like this. - */ - template class ::SymmetricTensor; - }; + private: + /** + * Reference to the tensor we + * access. + */ + tensor_type &tensor; + + /** + * Index of the row we access. + */ + const unsigned int row; + + /** + * Make the symmetric tensor + * classes a friend, since they are + * the only ones who can create + * objects like this. + */ + template class ::SymmetricTensor; + }; - } } } @@ -228,19 +235,42 @@ namespace internal /** - * Provide a class that stores symmetric tensors of rank 2 efficiently, - * i.e. only store half of the off-diagonal elements of the full tensor. + * Provide a class that stores symmetric tensors of rank 2,4,... efficiently, + * i.e. only store those off-diagonal elements of the full tensor that are not + * redundant. For example, for symmetric 2x2 tensors, this would be the + * elements 11, 22, and 12, while the element 21 is equal to the 12 element. + * + * Using this class for symmetric tensors of rank 2 has advantages over + * matrices in many cases since the dimension is known to the compiler as well + * as the location of the data. It is therefore possible to produce far more + * efficient code than for matrices with runtime-dependent dimension. It is + * also more efficient than using the more general Tensor class, + * since less elements are stored, and the class automatically makes sure that + * the tensor represents a symmetric object. + * + * For tensors of higher rank, the savings in storage are even higher. For + * example for the 3x3x3x3 tensors of rank 4, only 36 instead of the full 81 + * entries have to be stored. + * + * Tensors of rank 4 are considered symmetric if they are operators mapping + * symmetric rank-2 tensors onto symmetric rank-2 tensors. This entails + * certain symmetry properties on the elements in their 4-dimensional index + * space. + * + * Symmetric tensors are most often used in structural and fluid mechanics, + * where strains and stresses are usually symmetric tensors, and the + * stress-strain relationship is given by a symmetric rank-4 tensor. * - * Using this tensor class for objects of rank 2 has advantages over - * matrices in many cases since the dimension is known to the compiler - * as well as the location of the data. It is therefore possible to - * produce far more efficient code than for matrices with - * runtime-dependent dimension. + * Note that symmetric tensors only exist with even numbers of indices. In + * other words, the only objects that you can use are + * SymmetricTensor<2,dim>, SymmetricTensor<4,dim>, etc, but + * SymmetricTensor<1,dim> and SymmetricTensor<3,dim> do not + * exist and their use will most likely lead to compiler errors. * * @author Wolfgang Bangerth, 2005 */ -template -class SymmetricTensor<2,dim> +template +class SymmetricTensor { public: /** @@ -261,12 +291,6 @@ class SymmetricTensor<2,dim> * data types. */ static const unsigned int dimension = dim; - - /** - * Publish the rank of this tensor to - * the outside world. - */ - static const unsigned int rank = 2; /** * Default constructor. Creates a zero @@ -355,12 +379,33 @@ class SymmetricTensor<2,dim> */ double operator * (const SymmetricTensor &s) const; + /** + * Return a read-write reference + * to the indicated element. + */ + double & operator() (const TableIndices &indices); + + /** + * Return the value of the + * indicated element as a + * read-only reference. + * + * We return the requested value + * as a constant reference rather + * than by value since this + * object may hold data types + * that may be large, and we + * don't know here whether + * copying is expensive or not. + */ + double operator() (const TableIndices &indices) const; + /** * Access the elements of a row of this * symmetric tensor. This function is * called for constant tensors. */ - internal::SymmetricTensor::Rank2Accessors::Accessor + internal::SymmetricTensor::Accessor operator [] (const unsigned int row) const; /** @@ -368,7 +413,7 @@ class SymmetricTensor<2,dim> * symmetric tensor. This function is * called for non-constant tensors. */ - internal::SymmetricTensor::Rank2Accessors::Accessor + internal::SymmetricTensor::Accessor operator [] (const unsigned int row); /** @@ -424,74 +469,32 @@ namespace internal { namespace SymmetricTensor { - namespace Rank2Accessors + template + Accessor<2,dim,constness>:: + Accessor (tensor_type &tensor, + const unsigned int row) + : + tensor (tensor), + row (row) + {} + + + + template + typename Accessor<2,dim,constness>::reference + Accessor<2,dim,constness>:: + operator[] (const unsigned int column) { - template - Accessor:: - Accessor (base_tensor_type &base_tensor, - const unsigned int row) - : - base_tensor (base_tensor), - row (row) - { - Assert (row < dim, ExcIndexRange (row, 0, dim)); - } - - - - template - typename Accessor::reference - Accessor:: - operator[] (const unsigned int column) - { - Assert (column < dim, ExcIndexRange (column, 0, dim)); - - // first treat the main diagonal - // elements, which are stored - // consecutively at the beginning - if (row == column) - return base_tensor[row]; - - // the rest is messier and requires a - // few switches. if someone has a - // better idea, help is welcome - switch (dim) - { - case 2: - Assert (((row==1) && (column==0)) || ((row==0) && (column==1)), - ExcInternalError()); - return base_tensor[2]; - - case 3: - if (((row==0) && (column==1)) || - ((row==1) && (column==0))) - return base_tensor[3]; - else if (((row==0) && (column==2)) || - ((row==2) && (column==0))) - return base_tensor[4]; - else if (((row==1) && (column==2)) || - ((row==2) && (column==1))) - return base_tensor[5]; - else - Assert (false, ExcInternalError()); - - default: - Assert (false, ExcNotImplemented()); - } - - Assert (false, ExcInternalError()); - static double dummy_but_referenceable = 0; - return dummy_but_referenceable; - } + return tensor(TableIndices<2> (row, column)); } } } -template +template inline -SymmetricTensor<2,dim>::SymmetricTensor () +SymmetricTensor::SymmetricTensor () {} @@ -526,10 +529,10 @@ SymmetricTensor<2,3>::SymmetricTensor (const Tensor<2,3> &t) } -template +template inline -SymmetricTensor<2,dim> & -SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t) +SymmetricTensor & +SymmetricTensor::operator = (const SymmetricTensor &t) { data = t.data; return *this; @@ -537,30 +540,30 @@ SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t) -template +template inline bool -SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const +SymmetricTensor::operator == (const SymmetricTensor &t) const { return data == t.data; } -template +template inline bool -SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const +SymmetricTensor::operator != (const SymmetricTensor &t) const { return data != t.data; } -template +template inline -SymmetricTensor<2,dim> & -SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t) +SymmetricTensor & +SymmetricTensor::operator += (const SymmetricTensor &t) { data += t.data; return *this; @@ -568,10 +571,10 @@ SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t) -template +template inline -SymmetricTensor<2,dim> & -SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t) +SymmetricTensor & +SymmetricTensor::operator -= (const SymmetricTensor &t) { data -= t.data; return *this; @@ -579,10 +582,10 @@ SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t) -template +template inline -SymmetricTensor<2,dim> & -SymmetricTensor<2,dim>::operator *= (const double d) +SymmetricTensor & +SymmetricTensor::operator *= (const double d) { data *= d; return *this; @@ -590,10 +593,10 @@ SymmetricTensor<2,dim>::operator *= (const double d) -template +template inline -SymmetricTensor<2,dim> & -SymmetricTensor<2,dim>::operator /= (const double d) +SymmetricTensor & +SymmetricTensor::operator /= (const double d) { data /= d; return *this; @@ -601,10 +604,10 @@ SymmetricTensor<2,dim>::operator /= (const double d) -template +template inline -SymmetricTensor<2,dim> -SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const +SymmetricTensor +SymmetricTensor::operator + (const SymmetricTensor &t) const { SymmetricTensor tmp = *this; tmp.data += t.data; @@ -613,10 +616,10 @@ SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const -template +template inline -SymmetricTensor<2,dim> -SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const +SymmetricTensor +SymmetricTensor::operator - (const SymmetricTensor &t) const { SymmetricTensor tmp = *this; tmp.data -= t.data; @@ -625,10 +628,10 @@ SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const -template +template inline -SymmetricTensor<2,dim> -SymmetricTensor<2,dim>::operator - () const +SymmetricTensor +SymmetricTensor::operator - () const { SymmetricTensor tmp = *this; tmp.data = -tmp.data; @@ -637,59 +640,229 @@ SymmetricTensor<2,dim>::operator - () const -template +template inline void -SymmetricTensor<2,dim>::clear () +SymmetricTensor::clear () { data.clear (); } -template +template inline unsigned int -SymmetricTensor<2,dim>::memory_consumption () +SymmetricTensor::memory_consumption () { - return internal::SymmetricTensor::StorageType<2,dim>::memory_consumption (); + return + internal::SymmetricTensor::StorageType::memory_consumption (); } -template +template <> double -SymmetricTensor<2,dim>::operator * (const SymmetricTensor &s) const +SymmetricTensor<2,1>::operator * (const SymmetricTensor<2,1> &s) const { - double t = 0; - unsigned int i=0; - for (; i::n_independent_components; ++i) - t += 2 * data[i] * s.data[i]; - return t; + +template <> +double +SymmetricTensor<2,2>::operator * (const SymmetricTensor<2,2> &s) const +{ + return (data[0] * s.data[0] + + data[1] * s.data[1] + + 2*data[2] * s.data[2]); +} + + + +template <> +double +SymmetricTensor<2,3>::operator * (const SymmetricTensor<2,3> &s) const +{ + return (data[0] * s.data[0] + + data[1] * s.data[1] + + data[2] * s.data[2] + + 2*data[3] * s.data[3] + + 2*data[4] * s.data[4] + + 2*data[5] * s.data[5]); +} + + + +template <> +double & +SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices) +{ + const unsigned int rank = 2; + const unsigned int dim = 1; + for (unsigned int r=0; r +double +SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices) const +{ + const unsigned int rank = 2; + const unsigned int dim = 1; + for (unsigned int r=0; r +double & +SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices) +{ + const unsigned int rank = 2; + const unsigned int dim = 2; + for (unsigned int r=0; r +double +SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices) const +{ + const unsigned int rank = 2; + const unsigned int dim = 2; + for (unsigned int r=0; r +double & +SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices) +{ + const unsigned int rank = 2; + const unsigned int dim = 3; + for (unsigned int r=0; r sorted_indices (indices); + sorted_indices.sort (); + if ((sorted_indices[0]==0) && (sorted_indices[1]==1)) + return data[3]; + else if ((sorted_indices[0]==0) && (sorted_indices[1]==2)) + return data[4]; + else if ((sorted_indices[0]==1) && (sorted_indices[1]==2)) + return data[5]; + else + Assert (false, ExcInternalError()); + + static double dummy_but_referenceable = 0; + return dummy_but_referenceable; +} + -template -internal::SymmetricTensor::Rank2Accessors::Accessor -SymmetricTensor<2,dim>::operator [] (const unsigned int row) const +template <> +double +SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices) const +{ + const unsigned int rank = 2; + const unsigned int dim = 3; + for (unsigned int r=0; r sorted_indices (indices); + sorted_indices.sort (); + + if ((sorted_indices[0]==0) && (sorted_indices[1]==1)) + return data[3]; + else if ((sorted_indices[0]==0) && (sorted_indices[1]==2)) + return data[4]; + else if ((sorted_indices[0]==1) && (sorted_indices[1]==2)) + return data[5]; + else + Assert (false, ExcInternalError()); + + static double dummy_but_referenceable = 0; + return dummy_but_referenceable; +} + + + +template +internal::SymmetricTensor::Accessor +SymmetricTensor::operator [] (const unsigned int row) const { return - internal::SymmetricTensor::Rank2Accessors::Accessor (data, row); + internal::SymmetricTensor::Accessor (*this, row); } -template -internal::SymmetricTensor::Rank2Accessors::Accessor -SymmetricTensor<2,dim>::operator [] (const unsigned int row) +template +internal::SymmetricTensor::Accessor +SymmetricTensor::operator [] (const unsigned int row) { return - internal::SymmetricTensor::Rank2Accessors::Accessor (data, row); + internal::SymmetricTensor::Accessor (*this, row); } @@ -766,8 +939,8 @@ double determinant (const SymmetricTensor<2,3> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template -double trace (const SymmetricTensor<2,dim> &d) +template +double trace (const SymmetricTensor &d) { double t=0; for (unsigned int i=0; i &d) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template inline -SymmetricTensor<2,dim> -transpose (const SymmetricTensor<2,dim> &t) +SymmetricTensor +transpose (const SymmetricTensor &t) { return t; } -- 2.39.5