--- /dev/null
+//---------------------------- 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 <int rank, int dim> 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 <int dim>
+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 <tt>dimension()</tt>
+ * 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 <tt>factor</tt>,
+ * i.e. multiply all components by
+ * <tt>factor</tt>.
+ */
+ SymmetricTensor & operator *= (const double factor);
+
+ /**
+ * Scale the vector by
+ * <tt>1/factor</tt>.
+ */
+ SymmetricTensor & operator /= (const double factor);
+
+ /**
+ * Add two tensors. If possible, you
+ * should use <tt>operator +=</tt>
+ * instead since this does not need the
+ * creation of a temporary.
+ */
+ SymmetricTensor operator + (const SymmetricTensor &) const;
+
+ /**
+ * Subtract two tensors. If possible,
+ * you should use <tt>operator -=</tt>
+ * 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 <int dim>
+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 <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t)
+{
+ data = t.data;
+ return *this;
+}
+
+
+
+template <int dim>
+inline
+bool
+SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const
+{
+ return data == t.data;
+}
+
+
+
+template <int dim>
+inline
+bool
+SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const
+{
+ return data != t.data;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t)
+{
+ data += t.data;
+ return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t)
+{
+ data -= t.data;
+ return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator *= (const double d)
+{
+ data *= d;
+ return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator /= (const double d)
+{
+ data /= d;
+ return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const
+{
+ SymmetricTensor tmp = *this;
+ tmp.data += t.data;
+ return tmp;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const
+{
+ SymmetricTensor tmp = *this;
+ tmp.data -= t.data;
+ return tmp;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator - () const
+{
+ SymmetricTensor tmp = *this;
+ tmp.data = -tmp.data;
+ return tmp;
+}
+
+
+
+template <int dim>
+inline
+void
+SymmetricTensor<2,dim>::clear ()
+{
+ data.clear ();
+}
+
+
+
+template <int dim>
+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 <tt>dim==2</tt>.
+ *
+ * @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 <tt>dim==3</tt>.
+ *
+ * @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 <int dim>
+double trace (const SymmetricTensor<2,dim> &d)
+{
+ Assert (false, ExcNotImplemented());
+ return 0;
+
+// double t=0;
+// for (unsigned int i=0; i<dim; ++i)
+// t += d[i][i];
+// return t;
+}
+
+
+
+/**
+ * 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 object is known, it is acceptable to return the
+ * result by value, rather than by reference as a parameter.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+template <int dim>
+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 <int dim>
+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 <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator * (const SymmetricTensor<rank,dim> &t,
+ const double factor)
+{
+ SymmetricTensor<rank,dim> tt = t;
+ tt *= factor;
+ return tt;
+}
+
+
+
+/**
+ * Multiplication of a symmetric tensor of general rank with a scalar double
+ * from the left.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator * (const double factor,
+ const SymmetricTensor<rank,dim> &t)
+{
+ SymmetricTensor<rank,dim> tt = t;
+ tt *= factor;
+ return tt;
+}
+
+
+
+/**
+ * Division of a symmetric tensor of general rank by a scalar double.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator / (const SymmetricTensor<rank,dim> &t,
+ const double factor)
+{
+ SymmetricTensor<rank,dim> tt = t;
+ tt /= factor;
+ return tt;
+}
+
+
+
+#endif