#include <base/tensor.h>
+#include <base/table_indices.h>
template <int rank, int dim> class SymmetricTensor;
-template <int dim> class SymmetricTensor<2,dim>;
namespace internal
template <int rank, int dim>
struct AccessorTypes<rank, dim,true>
{
- typedef
- const typename StorageType<rank,dim>::base_tensor_type
- base_tensor_type;
+ typedef const ::SymmetricTensor<rank,dim> tensor_type;
typedef double reference;
};
template <int rank, int dim>
struct AccessorTypes<rank,dim,false>
{
- typedef
- typename StorageType<rank,dim>::base_tensor_type
- base_tensor_type;
+ typedef ::SymmetricTensor<rank,dim> tensor_type;
typedef double &reference;
};
+
+ template <int rank, int dim, bool constness>
+ 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 <tt>operator[]</tt>;
+ * the accessor thereby points to a row
+ * of the tensor. Calling
+ * <tt>operator[]</tt> 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 <int dim, bool constness>
+ 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 <int dim, bool constness>
- 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 <int,int> 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 <int,int> 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.
+ * 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 <tt>Tensor</tt> 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
+ * <tt>SymmetricTensor<2,dim></tt>, <tt>SymmetricTensor<4,dim></tt>, etc, but
+ * <tt>SymmetricTensor<1,dim></tt> and <tt>SymmetricTensor<3,dim></tt> do not
+ * exist and their use will most likely lead to compiler errors.
*
* @author Wolfgang Bangerth, 2005
*/
-template <int dim>
-class SymmetricTensor<2,dim>
+template <int rank, int dim>
+class SymmetricTensor
{
public:
/**
* 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
*/
double operator * (const SymmetricTensor &s) const;
+ /**
+ * Return a read-write reference
+ * to the indicated element.
+ */
+ double & operator() (const TableIndices<rank> &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<rank> &indices) const;
+
/**
* Access the elements of a row of this
* symmetric tensor. This function is
* called for constant tensors.
*/
- internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true>
+ internal::SymmetricTensor::Accessor<rank,dim,true>
operator [] (const unsigned int row) const;
/**
* symmetric tensor. This function is
* called for non-constant tensors.
*/
- internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false>
+ internal::SymmetricTensor::Accessor<rank,dim,false>
operator [] (const unsigned int row);
/**
{
namespace SymmetricTensor
{
- namespace Rank2Accessors
+ template <int dim, bool constness>
+ Accessor<2,dim,constness>::
+ Accessor (tensor_type &tensor,
+ const unsigned int row)
+ :
+ tensor (tensor),
+ row (row)
+ {}
+
+
+
+ template <int dim, bool constness>
+ typename Accessor<2,dim,constness>::reference
+ Accessor<2,dim,constness>::
+ operator[] (const unsigned int column)
{
- template <int dim, bool constness>
- Accessor<dim,constness>::
- Accessor (base_tensor_type &base_tensor,
- const unsigned int row)
- :
- base_tensor (base_tensor),
- row (row)
- {
- Assert (row < dim, ExcIndexRange (row, 0, dim));
- }
-
-
-
- template <int dim, bool constness>
- typename Accessor<dim,constness>::reference
- Accessor<dim,constness>::
- 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 <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim>::SymmetricTensor ()
+SymmetricTensor<rank,dim>::SymmetricTensor ()
{}
}
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator = (const SymmetricTensor<rank,dim> &t)
{
data = t.data;
return *this;
-template <int dim>
+template <int rank, int dim>
inline
bool
-SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const
+SymmetricTensor<rank,dim>::operator == (const SymmetricTensor<rank,dim> &t) const
{
return data == t.data;
}
-template <int dim>
+template <int rank, int dim>
inline
bool
-SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const
+SymmetricTensor<rank,dim>::operator != (const SymmetricTensor<rank,dim> &t) const
{
return data != t.data;
}
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator += (const SymmetricTensor<rank,dim> &t)
{
data += t.data;
return *this;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator -= (const SymmetricTensor<rank,dim> &t)
{
data -= t.data;
return *this;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator *= (const double d)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator *= (const double d)
{
data *= d;
return *this;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator /= (const double d)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator /= (const double d)
{
data /= d;
return *this;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator + (const SymmetricTensor &t) const
{
SymmetricTensor tmp = *this;
tmp.data += t.data;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator - (const SymmetricTensor &t) const
{
SymmetricTensor tmp = *this;
tmp.data -= t.data;
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator - () const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator - () const
{
SymmetricTensor tmp = *this;
tmp.data = -tmp.data;
-template <int dim>
+template <int rank, int dim>
inline
void
-SymmetricTensor<2,dim>::clear ()
+SymmetricTensor<rank,dim>::clear ()
{
data.clear ();
}
-template <int dim>
+template <int rank, int dim>
inline
unsigned int
-SymmetricTensor<2,dim>::memory_consumption ()
+SymmetricTensor<rank,dim>::memory_consumption ()
{
- return internal::SymmetricTensor::StorageType<2,dim>::memory_consumption ();
+ return
+ internal::SymmetricTensor::StorageType<rank,dim>::memory_consumption ();
}
-template <int dim>
+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<dim; ++i)
- t += data[i] * s.data[i];
+ return data[0] * s.data[0];
+}
- for (; i<internal::SymmetricTensor::StorageType<2,dim>::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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ return data[0];
+}
+
+
+
+template <>
+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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ return data[0];
+}
+
+
+
+template <>
+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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ // first treat the main diagonal
+ // elements, which are stored
+ // consecutively at the beginning
+ if (indices[0] == indices[1])
+ return data[indices[0]];
+
+ // the rest is messier and requires a few
+ // switches. at least for the 2x2 case it
+ // is reasonably simple
+ Assert (((indices[0]==1) && (indices[1]==0)) ||
+ ((indices[0]==0) && (indices[1]==1)),
+ ExcInternalError());
+ return data[2];
+}
+
+
+
+template <>
+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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ // first treat the main diagonal
+ // elements, which are stored
+ // consecutively at the beginning
+ if (indices[0] == indices[1])
+ return data[indices[0]];
+
+ // the rest is messier and requires a few
+ // switches. at least for the 2x2 case it
+ // is reasonably simple
+ Assert (((indices[0]==1) && (indices[1]==0)) ||
+ ((indices[0]==0) && (indices[1]==1)),
+ ExcInternalError());
+ return data[2];
}
+
+
+
+template <>
+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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ // first treat the main diagonal
+ // elements, which are stored
+ // consecutively at the beginning
+ if (indices[0] == indices[1])
+ return data[indices[0]];
+
+ // the rest is messier and requires a few
+ // switches, but simpler if we just sort
+ // our indices
+ TableIndices<2> 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 <int dim>
-internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true>
-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<rank; ++r)
+ Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+ // first treat the main diagonal
+ // elements, which are stored
+ // consecutively at the beginning
+ if (indices[0] == indices[1])
+ return data[indices[0]];
+
+ // the rest is messier and requires a few
+ // switches, but simpler if we just sort
+ // our indices
+ TableIndices<2> 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 <int rank, int dim>
+internal::SymmetricTensor::Accessor<rank,dim,true>
+SymmetricTensor<rank,dim>::operator [] (const unsigned int row) const
{
return
- internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true> (data, row);
+ internal::SymmetricTensor::Accessor<rank,dim,true> (*this, row);
}
-template <int dim>
-internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false>
-SymmetricTensor<2,dim>::operator [] (const unsigned int row)
+template <int rank, int dim>
+internal::SymmetricTensor::Accessor<rank,dim,false>
+SymmetricTensor<rank,dim>::operator [] (const unsigned int row)
{
return
- internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false> (data, row);
+ internal::SymmetricTensor::Accessor<rank,dim,false> (*this, row);
}
* @relates SymmetricTensor
* @author Wolfgang Bangerth, 2005
*/
-template <int dim>
-double trace (const SymmetricTensor<2,dim> &d)
+template <int rank, int dim>
+double trace (const SymmetricTensor<rank,dim> &d)
{
double t=0;
for (unsigned int i=0; i<dim; ++i)
* @relates SymmetricTensor
* @author Wolfgang Bangerth, 2005
*/
-template <int dim>
+template <int rank, int dim>
inline
-SymmetricTensor<2,dim>
-transpose (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim>
+transpose (const SymmetricTensor<rank,dim> &t)
{
return t;
}