From e3aaea5900a7084db2bd36f09ffaffc708a95d64 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 28 Mar 2005 17:31:03 +0000 Subject: [PATCH] Add a class implementing symmetric tensors of rank 2. git-svn-id: https://svn.dealii.org/trunk@10249 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/symmetric_tensor.h | 625 +++++++++++++++++++ 1 file changed, 625 insertions(+) create mode 100644 deal.II/base/include/base/symmetric_tensor.h diff --git a/deal.II/base/include/base/symmetric_tensor.h b/deal.II/base/include/base/symmetric_tensor.h new file mode 100644 index 0000000000..3af4dcd575 --- /dev/null +++ b/deal.II/base/include/base/symmetric_tensor.h @@ -0,0 +1,625 @@ +//---------------------------- symmetric_tensor.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 by the deal authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- symmetric_tensor.h --------------------------- +#ifndef __deal2__symmetric_tensor_h +#define __deal2__symmetric_tensor_h + + +template class SymmetricTensor; + + + +/** + * 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. + * + * 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. + * + * @author Wolfgang Bangerth, 2005 + */ +template +class SymmetricTensor<2,dim> +{ + public: + /** + * Provide a way to get the + * dimension of an object without + * explicit knowledge of it's + * data type. Implementation is + * this way instead of providing + * a function dimension() + * because now it is possible to + * get the dimension at compile + * time without the expansion and + * preevaluation of an inlined + * function; the compiler may + * therefore produce more + * efficient code and you may use + * this value to declare other + * 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 + * tensor. + */ + SymmetricTensor (); + + /** + * Constructor. Generate a symmetric + * tensor from a general one. Assumes + * that @p t is already symmetric, but + * this is not checked: we simply copy + * only a subset of elements. + */ + SymmetricTensor (const Tensor<2,dim> &t); + + /** + * Assignment operator. + */ + SymmetricTensor & operator = (const SymmetricTensor &); + + /** + * Test for equality of two tensors. + */ + bool operator == (const SymmetricTensor &) const; + + /** + * Test for inequality of two tensors. + */ + bool operator != (const SymmetricTensor &) const; + + /** + * Add another tensor. + */ + SymmetricTensor & operator += (const SymmetricTensor &); + + /** + * Subtract another tensor. + */ + SymmetricTensor & operator -= (const SymmetricTensor &); + + /** + * Scale the tensor by factor, + * i.e. multiply all components by + * factor. + */ + SymmetricTensor & operator *= (const double factor); + + /** + * Scale the vector by + * 1/factor. + */ + SymmetricTensor & operator /= (const double factor); + + /** + * Add two tensors. If possible, you + * should use operator += + * instead since this does not need the + * creation of a temporary. + */ + SymmetricTensor operator + (const SymmetricTensor &) const; + + /** + * Subtract two tensors. If possible, + * you should use operator -= + * instead since this does not need the + * creation of a temporary. + */ + SymmetricTensor operator - (const SymmetricTensor &) const; + + /** + * Unary minus operator. Negate all + * entries of a tensor. + */ + SymmetricTensor operator - () const; + + /** + * Return the Frobenius-norm of a tensor, + * i.e. the square root of the sum of + * squares of all entries. + */ + double norm () const; + + /** + * Reset all values to zero. + * + * Note that this is partly inconsistent + * with the semantics of the @p clear() + * member functions of the STL and of + * several other classes within deal.II + * which not only reset the values of + * stored elements to zero, but release + * all memory and return the object into + * a virginial state. However, since the + * size of objects of the present type is + * determined by its template parameters, + * resizing is not an option, and indeed + * the state where all elements have a + * zero value is the state right after + * construction of such an object. + */ + void clear (); + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this + * object. + */ + static unsigned int memory_consumption (); + + + private: + /** + * Number of independent components of a + * symmetric tensor of rank 2. We store + * only the upper right half of it. + */ + static const unsigned int + n_tensor_components = (dim*dim + dim)/2; + + /** + * Declare the type in which we actually + * store the data. + */ + typedef Tensor<1,n_tensor_components> StorageType; + + /** + * Data storage for a symmetric tensor. + */ + StorageType data; +}; + + + +// ------------------------- inline functions ------------------------ + +template +inline +SymmetricTensor<2,dim>::SymmetricTensor () +{} + + + +template <> +inline +SymmetricTensor<2,2>::SymmetricTensor (const Tensor<2,2> &t) +{ + Assert (t[0][1] == t[1][0], ExcInternalError()); + + data[0] = t[0][0]; + data[1] = t[1][1]; + data[2] = t[0][1]; +} + + + +template <> +inline +SymmetricTensor<2,3>::SymmetricTensor (const Tensor<2,3> &t) +{ + Assert (t[0][1] == t[1][0], ExcInternalError()); + Assert (t[0][2] == t[2][0], ExcInternalError()); + Assert (t[1][2] == t[2][1], ExcInternalError()); + + data[0] = t[0][0]; + data[1] = t[1][1]; + data[2] = t[2][2]; + data[3] = t[0][1]; + data[4] = t[0][2]; + data[5] = t[1][2]; +} + + +template +inline +SymmetricTensor<2,dim> & +SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t) +{ + data = t.data; + return *this; +} + + + +template +inline +bool +SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const +{ + return data == t.data; +} + + + +template +inline +bool +SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const +{ + return data != t.data; +} + + + +template +inline +SymmetricTensor<2,dim> & +SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t) +{ + data += t.data; + return *this; +} + + + +template +inline +SymmetricTensor<2,dim> & +SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t) +{ + data -= t.data; + return *this; +} + + + +template +inline +SymmetricTensor<2,dim> & +SymmetricTensor<2,dim>::operator *= (const double d) +{ + data *= d; + return *this; +} + + + +template +inline +SymmetricTensor<2,dim> & +SymmetricTensor<2,dim>::operator /= (const double d) +{ + data /= d; + return *this; +} + + + +template +inline +SymmetricTensor<2,dim> +SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const +{ + SymmetricTensor tmp = *this; + tmp.data += t.data; + return tmp; +} + + + +template +inline +SymmetricTensor<2,dim> +SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const +{ + SymmetricTensor tmp = *this; + tmp.data -= t.data; + return tmp; +} + + + +template +inline +SymmetricTensor<2,dim> +SymmetricTensor<2,dim>::operator - () const +{ + SymmetricTensor tmp = *this; + tmp.data = -tmp.data; + return tmp; +} + + + +template +inline +void +SymmetricTensor<2,dim>::clear () +{ + data.clear (); +} + + + +template +inline +unsigned int +SymmetricTensor<2,dim>::memory_consumption () +{ + return StorageType::memory_consumption (); +} + + + +template <> +double +SymmetricTensor<2,1>::norm () const +{ + return std::sqrt(data[0]*data[0]); +} + + + +template <> +double +SymmetricTensor<2,2>::norm () const +{ + return std::sqrt(data[0]*data[0] + data[1]*data[1] + 2*data[2]*data[2]); +} + + + +template <> +double +SymmetricTensor<2,3>::norm () const +{ + return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2] + + 2*data[3]*data[3] + 2*data[4]*data[4] + 2*data[5]*data[5]); +} + + +/* ----------------- Non-member functions operating on tensors. ------------ */ + +/** + * 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. + * + * @relates Tensor + * @author Wolfgang Bangerth, 2005 + */ +inline +double determinant (const SymmetricTensor<1,1> &t) +{ + Assert (false, ExcNotImplemented()); + return 0; + +// return t[0]; +} + + + +/** + * Compute the determinant of a tensor or rank 2, here for dim==2. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +inline +double determinant (const SymmetricTensor<2,2> &t) +{ + Assert (false, ExcNotImplemented()); + return 0; + +// return ((t[0][0] * t[1][1]) - +// (t[1][0] * t[0][1])); +} + + + + +/** + * Compute the determinant of a tensor or rank 2, here for dim==3. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +inline +double determinant (const SymmetricTensor<2,3> &t) +{ + Assert (false, ExcNotImplemented()); + return 0; + +// // get this using Maple: +// // with(linalg); +// // a := matrix(3,3); +// // x := det(a); +// // readlib(C); +// // C(x, optimized); +// return ( t[0][0]*t[1][1]*t[2][2] +// -t[0][0]*t[1][2]*t[2][1] +// -t[1][0]*t[0][1]*t[2][2] +// +t[1][0]*t[0][2]*t[2][1] +// +t[2][0]*t[0][1]*t[1][2] +// -t[2][0]*t[0][2]*t[1][1] ); +} + + + +/** + * Compute and return the trace of a tensor of rank 2, i.e. the sum of + * its diagonal entries. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +template +double trace (const SymmetricTensor<2,dim> &d) +{ + Assert (false, ExcNotImplemented()); + return 0; + +// double t=0; +// for (unsigned int i=0; i +inline +SymmetricTensor<2,dim> +invert (const SymmetricTensor<2,dim> &t) +{ + Assert (false, ExcNotImplemented()); + return SymmetricTensor<2,dim>(); + +// SymmetricTensor<2,dim> return_tensor; +// switch (dim) +// { +// case 1: +// return_tensor[0][0] = 1.0/t[0][0]; +// return return_tensor; +// case 2: +// // this is Maple output, +// // thus a bit unstructured +// { +// const double t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]); +// return_tensor[0][0] = t[1][1]*t4; +// return_tensor[0][1] = -t[0][1]*t4; +// return_tensor[1][0] = -t[1][0]*t4; +// return_tensor[1][1] = t[0][0]*t4; +// return return_tensor; +// }; + +// case 3: +// { +// const double t4 = t[0][0]*t[1][1], +// t6 = t[0][0]*t[1][2], +// t8 = t[0][1]*t[1][0], +// t00 = t[0][2]*t[1][0], +// t01 = t[0][1]*t[2][0], +// t04 = t[0][2]*t[2][0], +// t07 = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+ +// t00*t[2][1]+t01*t[1][2]-t04*t[1][1]); +// return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07; +// return_tensor[0][1] = -(t[0][1]*t[2][2]-t[0][2]*t[2][1])*t07; +// return_tensor[0][2] = -(-t[0][1]*t[1][2]+t[0][2]*t[1][1])*t07; +// return_tensor[1][0] = -(t[1][0]*t[2][2]-t[1][2]*t[2][0])*t07; +// return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07; +// return_tensor[1][2] = -(t6-t00)*t07; +// return_tensor[2][0] = -(-t[1][0]*t[2][1]+t[1][1]*t[2][0])*t07; +// return_tensor[2][1] = -(t[0][0]*t[2][1]-t01)*t07; +// return_tensor[2][2] = (t4-t8)*t07; +// return return_tensor; +// }; + +// // if desired, take over the +// // inversion of a 4x4 tensor +// // from the FullMatrix + +// default: +// AssertThrow (false, ExcNotImplemented()); +// }; +// return return_tensor; +} + + + +/** + * Return the transpose of the given symmetric tensor. Since we are working + * with symmetric objects, the transpose is of course the same as the original + * tensor. This function mainly exists for compatibility with the Tensor + * class. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +template +inline +SymmetricTensor<2,dim> +transpose (const SymmetricTensor<2,dim> &t) +{ + return t; +} + + +/** + * Multiplication of a symmetric tensor of general rank with a scalar double + * from the right. + * + * @relates SymmetricTensor + */ +template +inline +SymmetricTensor +operator * (const SymmetricTensor &t, + const double factor) +{ + SymmetricTensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Multiplication of a symmetric tensor of general rank with a scalar double + * from the left. + * + * @relates SymmetricTensor + */ +template +inline +SymmetricTensor +operator * (const double factor, + const SymmetricTensor &t) +{ + SymmetricTensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Division of a symmetric tensor of general rank by a scalar double. + * + * @relates SymmetricTensor + */ +template +inline +SymmetricTensor +operator / (const SymmetricTensor &t, + const double factor) +{ + SymmetricTensor tt = t; + tt /= factor; + return tt; +} + + + +#endif -- 2.39.5