symmetric_tensor.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__symmetric_tensor_h
#define dealii__symmetric_tensor_h


#include <deal.II/base/tensor.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/table_indices.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/vectorization.h>

DEAL_II_NAMESPACE_OPEN

template <int rank, int dim, typename Number=double> class SymmetricTensor;

template <int dim, typename Number> SymmetricTensor<2,dim,Number>
unit_symmetric_tensor ();
template <int dim, typename Number> SymmetricTensor<4,dim,Number>
deviator_tensor ();
template <int dim, typename Number> SymmetricTensor<4,dim,Number>
identity_tensor ();
template <int dim, typename Number> SymmetricTensor<2,dim,Number>
invert (const SymmetricTensor<2,dim,Number> &);
template <int dim, typename Number> SymmetricTensor<4,dim,Number>
invert (const SymmetricTensor<4,dim,Number> &);
template <int dim2, typename Number> Number
trace (const SymmetricTensor<2,dim2,Number> &);

template <int dim, typename Number> SymmetricTensor<2,dim,Number>
deviator (const SymmetricTensor<2,dim,Number> &);
template <int dim, typename Number> Number
determinant (const SymmetricTensor<2,dim,Number> &);



namespace internal
{
  namespace SymmetricTensorAccessors
  {
    inline
    TableIndices<2> merge (const TableIndices<2> &previous_indices,
                           const unsigned int     new_index,
                           const unsigned int     position)
    {
      Assert (position < 2, ExcIndexRange (position, 0, 2));

      if (position == 0)
        return TableIndices<2>(new_index);
      else
        return TableIndices<2>(previous_indices[0], new_index);
    }



    inline
    TableIndices<4> merge (const TableIndices<4> &previous_indices,
                           const unsigned int     new_index,
                           const unsigned int     position)
    {
      Assert (position < 4, ExcIndexRange (position, 0, 4));

      switch (position)
        {
        case 0:
          return TableIndices<4>(new_index);
        case 1:
          return TableIndices<4>(previous_indices[0],
                                 new_index);
        case 2:
          return TableIndices<4>(previous_indices[0],
                                 previous_indices[1],
                                 new_index);
        case 3:
          return TableIndices<4>(previous_indices[0],
                                 previous_indices[1],
                                 previous_indices[2],
                                 new_index);
        }
      Assert (false, ExcInternalError());
      return TableIndices<4>();
    }


    template <int rank1, int rank2, int dim, typename Number>
    struct double_contraction_result
    {
      typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
    };


    template <int dim, typename Number>
    struct double_contraction_result<2,2,dim,Number>
    {
      typedef Number type;
    };



    template <int rank, int dim, typename Number>
    struct StorageType;

    template <int dim, typename Number>
    struct StorageType<2,dim,Number>
    {
      static const unsigned int
      n_independent_components = (dim*dim + dim)/2;

      typedef Tensor<1,n_independent_components,Number> base_tensor_type;
    };



    template <int dim, typename Number>
    struct StorageType<4,dim,Number>
    {
      static const unsigned int
      n_rank2_components = (dim*dim + dim)/2;

      static const unsigned int
      n_independent_components = (n_rank2_components *
                                  StorageType<2,dim,Number>::n_independent_components);

      typedef Tensor<2,n_rank2_components,Number> base_tensor_type;
    };



    template <int rank, int dim, bool constness, typename Number>
    struct AccessorTypes;

    template <int rank, int dim, typename Number>
    struct AccessorTypes<rank,dim,true,Number>
    {
      typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;

      typedef Number reference;
    };

    template <int rank, int dim, typename Number>
    struct AccessorTypes<rank,dim,false,Number>
    {
      typedef ::SymmetricTensor<rank,dim,Number> tensor_type;

      typedef Number &reference;
    };


    template <int rank, int dim, bool constness, int P, typename Number>
    class Accessor
    {
    public:
      typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
      typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;

    private:
      Accessor (tensor_type              &tensor,
                const TableIndices<rank> &previous_indices);

      Accessor (const Accessor &a);

    public:

      Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);

      Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;

    private:
      tensor_type             &tensor;
      const TableIndices<rank> previous_indices;

      // declare some other classes
      // as friends. make sure to
      // work around bugs in some
      // compilers
      template <int,int,typename> friend class ::SymmetricTensor;
      template <int,int,bool,int,typename>
      friend class Accessor;
#  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
      friend class ::SymmetricTensor<rank,dim,Number>;
      friend class Accessor<rank,dim,constness,P+1,Number>;
#  endif
    };



    template <int rank, int dim, bool constness, typename Number>
    class Accessor<rank,dim,constness,1,Number>
    {
    public:
      typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
      typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;

    private:
      Accessor (tensor_type              &tensor,
                const TableIndices<rank> &previous_indices);

      Accessor ();

      Accessor (const Accessor &a);

    public:

      reference operator [] (const unsigned int);

      reference operator [] (const unsigned int) const;

    private:
      tensor_type             &tensor;
      const TableIndices<rank> previous_indices;

      // declare some other classes
      // as friends. make sure to
      // work around bugs in some
      // compilers
      template <int,int,typename> friend class ::SymmetricTensor;
      template <int,int,bool,int,typename>
      friend class SymmetricTensorAccessors::Accessor;
#  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
      friend class ::SymmetricTensor<rank,dim,Number>;
      friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
#  endif
    };
  }
}



template <int rank, int dim, typename Number>
class SymmetricTensor
{
public:
  static const unsigned int dimension = dim;

  static const unsigned int n_independent_components
    = internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::
      n_independent_components;

  SymmetricTensor ();

  SymmetricTensor (const Tensor<2,dim,Number> &t);

  SymmetricTensor (const Number (&array) [n_independent_components]);

  template <typename OtherNumber>
  explicit
  SymmetricTensor (const SymmetricTensor<rank,dim,OtherNumber> &initializer);

  SymmetricTensor &operator = (const Number d);

  operator Tensor<rank,dim,Number> () const;

  bool operator == (const SymmetricTensor &) const;

  bool operator != (const SymmetricTensor &) const;

  SymmetricTensor &operator += (const SymmetricTensor &);

  SymmetricTensor &operator -= (const SymmetricTensor &);

  SymmetricTensor &operator *= (const Number factor);

  SymmetricTensor &operator /= (const Number factor);

  SymmetricTensor   operator + (const SymmetricTensor &s) const;

  SymmetricTensor   operator - (const SymmetricTensor &s) const;

  SymmetricTensor   operator - () const;

  typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
  operator * (const SymmetricTensor<2,dim,Number> &s) const;

  typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
  operator * (const SymmetricTensor<4,dim,Number> &s) const;

  Number &operator() (const TableIndices<rank> &indices);

  Number operator() (const TableIndices<rank> &indices) const;

  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
  operator [] (const unsigned int row) const;

  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
  operator [] (const unsigned int row);

  Number
  operator [] (const TableIndices<rank> &indices) const;

  Number &
  operator [] (const TableIndices<rank> &indices);

  Number
  access_raw_entry (const unsigned int unrolled_index) const;

  Number &
  access_raw_entry (const unsigned int unrolled_index);

  typename numbers::NumberTraits<Number>::real_type
  norm () const;

  static
  unsigned int
  component_to_unrolled_index (const TableIndices<rank> &indices);

  static
  TableIndices<rank>
  unrolled_to_component_indices (const unsigned int i);

  void clear ();

  static std::size_t memory_consumption ();

  template <class Archive>
  void serialize(Archive &ar, const unsigned int version);

private:
  typedef
  internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>
  base_tensor_descriptor;

  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;

  base_tensor_type data;

  template <int, int, typename> friend class SymmetricTensor;

  template <int dim2, typename Number2>
  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);

  template <int dim2, typename Number2>
  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);

  template <int dim2, typename Number2>
  friend SymmetricTensor<2,dim2,Number2>
  deviator (const SymmetricTensor<2,dim2,Number2> &t);

  template <int dim2, typename Number2>
  friend SymmetricTensor<2,dim2,Number2> unit_symmetric_tensor ();

  template <int dim2, typename Number2>
  friend SymmetricTensor<4,dim2,Number2> deviator_tensor ();

  template <int dim2, typename Number2>
  friend SymmetricTensor<4,dim2,Number2> identity_tensor ();

  template <int dim2, typename Number2>
  friend SymmetricTensor<2,dim2,Number2> invert (const SymmetricTensor<2,dim2,Number2> &);

  template <int dim2, typename Number2>
  friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &);
};



// ------------------------- inline functions ------------------------

#ifndef DOXYGEN

namespace internal
{
  namespace SymmetricTensorAccessors
  {
    template <int rank, int dim, bool constness, int P, typename Number>
    Accessor<rank,dim,constness,P,Number>::
    Accessor (tensor_type              &tensor,
              const TableIndices<rank> &previous_indices)
      :
      tensor (tensor),
      previous_indices (previous_indices)
    {}


    template <int rank, int dim, bool constness, int P, typename Number>
    Accessor<rank,dim,constness,P,Number>::
    Accessor (const Accessor &a)
      :
      tensor (a.tensor),
      previous_indices (a.previous_indices)
    {}



    template <int rank, int dim, bool constness, int P, typename Number>
    Accessor<rank,dim,constness,P-1,Number>
    Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
    {
      return Accessor<rank,dim,constness,P-1,Number> (tensor,
                                                      merge (previous_indices, i, rank-P));
    }



    template <int rank, int dim, bool constness, int P, typename Number>
    Accessor<rank,dim,constness,P-1,Number>
    Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
    {
      return Accessor<rank,dim,constness,P-1,Number> (tensor,
                                                      merge (previous_indices, i, rank-P));
    }



    template <int rank, int dim, bool constness, typename Number>
    Accessor<rank,dim,constness,1,Number>::
    Accessor (tensor_type              &tensor,
              const TableIndices<rank> &previous_indices)
      :
      tensor (tensor),
      previous_indices (previous_indices)
    {}



    template <int rank, int dim, bool constness, typename Number>
    Accessor<rank,dim,constness,1,Number>::
    Accessor (const Accessor &a)
      :
      tensor (a.tensor),
      previous_indices (a.previous_indices)
    {}



    template <int rank, int dim, bool constness, typename Number>
    typename Accessor<rank,dim,constness,1,Number>::reference
    Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
    {
      return tensor(merge (previous_indices, i, rank-1));
    }


    template <int rank, int dim, bool constness, typename Number>
    typename Accessor<rank,dim,constness,1,Number>::reference
    Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
    {
      return tensor(merge (previous_indices, i, rank-1));
    }
  }
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>::SymmetricTensor ()
{}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Tensor<2,dim,Number> &t)
{
  Assert (rank == 2, ExcNotImplemented());
  switch (dim)
    {
    case 2:
      Assert (t[0][1] == t[1][0], ExcInternalError());

      data[0] = t[0][0];
      data[1] = t[1][1];
      data[2] = t[0][1];

      break;
    case 3:
      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];

      break;
    default:
      for (unsigned int d=0; d<dim; ++d)
        for (unsigned int e=0; e<d; ++e)
          Assert(t[d][e] == t[e][d], ExcInternalError());

      for (unsigned int d=0; d<dim; ++d)
        data[d] = t[d][d];

      for (unsigned int d=0, c=0; d<dim; ++d)
        for (unsigned int e=d+1; e<dim; ++e, ++c)
          data[dim+c] = t[d][e];
    }
}



template <int rank, int dim, typename Number>
template <typename OtherNumber>
inline
SymmetricTensor<rank,dim,Number>::
SymmetricTensor (const SymmetricTensor<rank,dim,OtherNumber> &initializer)
{
  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
    data[i] = initializer.data[i];
}




template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
  :
  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
{
  // ensure that the reinterpret_cast above actually works
  Assert (sizeof(typename base_tensor_type::array_type)
          == sizeof(array),
          ExcInternalError());
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number> &
SymmetricTensor<rank,dim,Number>::operator = (const Number d)
{
  Assert (d==Number(), ExcMessage ("Only assignment with zero is allowed"));
  (void) d;

  data = 0;

  return *this;
}


namespace internal
{
  namespace SymmetricTensor
  {
    template <int dim, typename Number>
    ::Tensor<2,dim,Number>
    convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
    {
      Number t[dim][dim];

      // diagonal entries are stored first
      for (unsigned int d=0; d<dim; ++d)
        t[d][d] = s.access_raw_entry(d);

      // off-diagonal entries come next, row by row
      for (unsigned int d=0, c=0; d<dim; ++d)
        for (unsigned int e=d+1; e<dim; ++e, ++c)
          {
            t[d][e] = s.access_raw_entry(dim+c);
            t[e][d] = s.access_raw_entry(dim+c);
          }
      return ::Tensor<2,dim,Number>(t);
    }


    template <int dim, typename Number>
    ::Tensor<4,dim,Number>
    convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
    {
      // utilize the symmetry properties of SymmetricTensor<4,dim>
      // discussed in the class documentation to avoid accessing all
      // independent elements of the input tensor more than once
      ::Tensor<4,dim,Number> t;

      for (unsigned int i=0; i<dim; ++i)
        for (unsigned int j=i; j<dim; ++j)
          for (unsigned int k=0; k<dim; ++k)
            for (unsigned int l=k; l<dim; ++l)
              t[TableIndices<4>(i,j,k,l)]
                = t[TableIndices<4>(i,j,l,k)]
                  = t[TableIndices<4>(j,i,k,l)]
                    = t[TableIndices<4>(j,i,l,k)]
                      = st[TableIndices<4>(i,j,k,l)];

      return t;
    }
  }
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>::
operator Tensor<rank,dim,Number> () const
{
  return internal::SymmetricTensor::convert_to_tensor (*this);
}



template <int rank, int dim, typename Number>
inline
bool
SymmetricTensor<rank,dim,Number>::operator ==
(const SymmetricTensor<rank,dim,Number> &t) const
{
  return data == t.data;
}



template <int rank, int dim, typename Number>
inline
bool
SymmetricTensor<rank,dim,Number>::operator !=
(const SymmetricTensor<rank,dim,Number> &t) const
{
  return data != t.data;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number> &
SymmetricTensor<rank,dim,Number>::operator +=
(const SymmetricTensor<rank,dim,Number> &t)
{
  data += t.data;
  return *this;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number> &
SymmetricTensor<rank,dim,Number>::operator -=
(const SymmetricTensor<rank,dim,Number> &t)
{
  data -= t.data;
  return *this;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number> &
SymmetricTensor<rank,dim,Number>::operator *= (const Number d)
{
  data *= d;
  return *this;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number> &
SymmetricTensor<rank,dim,Number>::operator /= (const Number d)
{
  data /= d;
  return *this;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
SymmetricTensor<rank,dim,Number>::operator + (const SymmetricTensor &t) const
{
  SymmetricTensor tmp = *this;
  tmp.data += t.data;
  return tmp;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
SymmetricTensor<rank,dim,Number>::operator - (const SymmetricTensor &t) const
{
  SymmetricTensor tmp = *this;
  tmp.data -= t.data;
  return tmp;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
SymmetricTensor<rank,dim,Number>::operator - () const
{
  SymmetricTensor tmp = *this;
  tmp.data = -tmp.data;
  return tmp;
}



template <int rank, int dim, typename Number>
inline
void
SymmetricTensor<rank,dim,Number>::clear ()
{
  data.clear ();
}



template <int rank, int dim, typename Number>
inline
std::size_t
SymmetricTensor<rank,dim,Number>::memory_consumption ()
{
  // all memory consists of statically allocated memory of the current
  // object, no pointers
  return sizeof(SymmetricTensor<rank,dim,Number>);
}



namespace internal
{

  template <int dim, typename Number>
  inline
  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
                              const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
  {
    switch (dim)
      {
      case 1:
        return data[0] * sdata[0];
      default:
        // Start with the non-diagonal part to avoid some multiplications by
        // 2.
        Number sum = data[dim] * sdata[dim];
        for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
          sum += data[d] * sdata[d];
        sum += sum; // sum = sum * 2.;

        // Now add the contributions from the diagonal
        for (unsigned int d=0; d<dim; ++d)
          sum += data[d] * sdata[d];
        return sum;
      }
  }



  template <int dim, typename Number>
  inline
  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
                              const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
  {
    const unsigned int data_dim =
      SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components;
    Number tmp [data_dim];
    for (unsigned int i=0; i<data_dim; ++i)
      tmp[i] = perform_double_contraction<dim,Number>(data[i], sdata);
    return ::SymmetricTensor<2,dim,Number>(tmp);
  }



  template <int dim, typename Number>
  inline
  typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type
  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
                              const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
  {
    typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type tmp;
    for (unsigned int i=0; i<tmp.dimension; ++i)
      {
        // Start with the non-diagonal part
        Number sum = data[dim] * sdata[dim][i];
        for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
          sum += data[d] * sdata[d][i];
        sum += sum; // sum = sum * 2.;

        // Now add the contributions from the diagonal
        for (unsigned int d=0; d<dim; ++d)
          sum += data[d] * sdata[d][i];
        tmp[i] = sum;
      }
    return tmp;
  }



  template <int dim, typename Number>
  inline
  typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type
  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
                              const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
  {
    const unsigned int data_dim =
      SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components;
    typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type tmp;
    for (unsigned int i=0; i<data_dim; ++i)
      for (unsigned int j=0; j<data_dim; ++j)
        {
          // Start with the non-diagonal part
          for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
            tmp[i][j] += data[i][d] * sdata[d][j];
          tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;

          // Now add the contributions from the diagonal
          for (unsigned int d=0; d<dim; ++d)
            tmp[i][j] += data[i][d] * sdata[d][j];
        }
    return tmp;
  }

} // end of namespace internal



template <int rank, int dim, typename Number>
inline
typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<2,dim,Number> &s) const
{
  // need to have two different function calls
  // because a scalar and rank-2 tensor are not
  // the same data type (see internal function
  // above)
  return internal::perform_double_contraction<dim,Number> (data, s.data);
}



template <int rank, int dim, typename Number>
inline
typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<4,dim,Number> &s) const
{
  typename internal::SymmetricTensorAccessors::
  double_contraction_result<rank,4,dim,Number>::type tmp;
  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
  return tmp;
}



// internal namespace to switch between the
// access of different tensors. There used to
// be explicit instantiations before for
// different ranks and dimensions, but since
// we now allow for templates on the data
// type, and since we cannot partially
// specialize the implementation, this got
// into a separate namespace
namespace internal
{
  template <int dim, typename Number>
  inline
  Number &
  symmetric_tensor_access (const TableIndices<2> &indices,
                           typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
  {
    // 1d is very simple and done first
    if (dim == 1)
      return data[0];

    // 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.
    switch (dim)
      {
      case 2:
        // 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];

      default:
        // to do the rest, sort our indices before comparing
      {
        TableIndices<2> sorted_indices (indices);
        sorted_indices.sort ();

        for (unsigned int d=0, c=0; d<dim; ++d)
          for (unsigned int e=d+1; e<dim; ++e, ++c)
            if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
              return data[dim+c];
        Assert (false, ExcInternalError());
      }
      }

    static Number dummy_but_referenceable = Number();
    return dummy_but_referenceable;
  }



  template <int dim, typename Number>
  inline
  Number
  symmetric_tensor_access (const TableIndices<2> &indices,
                           const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
  {
    // 1d is very simple and done first
    if (dim == 1)
      return data[0];

    // 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.
    switch (dim)
      {
      case 2:
        // 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];

      default:
        // to do the rest, sort our indices before comparing
      {
        TableIndices<2> sorted_indices (indices);
        sorted_indices.sort ();

        for (unsigned int d=0, c=0; d<dim; ++d)
          for (unsigned int e=d+1; e<dim; ++e, ++c)
            if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
              return data[dim+c];
        Assert (false, ExcInternalError());
      }
      }

    static Number dummy_but_referenceable = Number();
    return dummy_but_referenceable;
  }



  template <int dim, typename Number>
  inline
  Number &
  symmetric_tensor_access (const TableIndices<4> &indices,
                           typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
  {
    switch (dim)
      {
      case 1:
        return data[0][0];

      case 2:
        // each entry of the tensor can be
        // thought of as an entry in a
        // matrix that maps the rolled-out
        // rank-2 tensors into rolled-out
        // rank-2 tensors. this is the
        // format in which we store rank-4
        // tensors. determine which
        // position the present entry is
        // stored in
      {
        unsigned int base_index[2] ;
        if ((indices[0] == 0) && (indices[1] == 0))
          base_index[0] = 0;
        else if ((indices[0] == 1) && (indices[1] == 1))
          base_index[0] = 1;
        else
          base_index[0] = 2;

        if ((indices[2] == 0) && (indices[3] == 0))
          base_index[1] = 0;
        else if ((indices[2] == 1) && (indices[3] == 1))
          base_index[1] = 1;
        else
          base_index[1] = 2;

        return data[base_index[0]][base_index[1]];
      }

      case 3:
        // each entry of the tensor can be
        // thought of as an entry in a
        // matrix that maps the rolled-out
        // rank-2 tensors into rolled-out
        // rank-2 tensors. this is the
        // format in which we store rank-4
        // tensors. determine which
        // position the present entry is
        // stored in
      {
        unsigned int base_index[2] ;
        if ((indices[0] == 0) && (indices[1] == 0))
          base_index[0] = 0;
        else if ((indices[0] == 1) && (indices[1] == 1))
          base_index[0] = 1;
        else if ((indices[0] == 2) && (indices[1] == 2))
          base_index[0] = 2;
        else if (((indices[0] == 0) && (indices[1] == 1)) ||
                 ((indices[0] == 1) && (indices[1] == 0)))
          base_index[0] = 3;
        else if (((indices[0] == 0) && (indices[1] == 2)) ||
                 ((indices[0] == 2) && (indices[1] == 0)))
          base_index[0] = 4;
        else
          {
            Assert (((indices[0] == 1) && (indices[1] == 2)) ||
                    ((indices[0] == 2) && (indices[1] == 1)),
                    ExcInternalError());
            base_index[0] = 5;
          }

        if ((indices[2] == 0) && (indices[3] == 0))
          base_index[1] = 0;
        else if ((indices[2] == 1) && (indices[3] == 1))
          base_index[1] = 1;
        else if ((indices[2] == 2) && (indices[3] == 2))
          base_index[1] = 2;
        else if (((indices[2] == 0) && (indices[3] == 1)) ||
                 ((indices[2] == 1) && (indices[3] == 0)))
          base_index[1] = 3;
        else if (((indices[2] == 0) && (indices[3] == 2)) ||
                 ((indices[2] == 2) && (indices[3] == 0)))
          base_index[1] = 4;
        else
          {
            Assert (((indices[2] == 1) && (indices[3] == 2)) ||
                    ((indices[2] == 2) && (indices[3] == 1)),
                    ExcInternalError());
            base_index[1] = 5;
          }

        return data[base_index[0]][base_index[1]];
      }

      default:
        Assert (false, ExcNotImplemented());
      }

    static Number dummy;
    return dummy;
  }


  template <int dim, typename Number>
  inline
  Number
  symmetric_tensor_access (const TableIndices<4> &indices,
                           const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
  {
    switch (dim)
      {
      case 1:
        return data[0][0];

      case 2:
        // each entry of the tensor can be
        // thought of as an entry in a
        // matrix that maps the rolled-out
        // rank-2 tensors into rolled-out
        // rank-2 tensors. this is the
        // format in which we store rank-4
        // tensors. determine which
        // position the present entry is
        // stored in
      {
        unsigned int base_index[2] ;
        if ((indices[0] == 0) && (indices[1] == 0))
          base_index[0] = 0;
        else if ((indices[0] == 1) && (indices[1] == 1))
          base_index[0] = 1;
        else
          base_index[0] = 2;

        if ((indices[2] == 0) && (indices[3] == 0))
          base_index[1] = 0;
        else if ((indices[2] == 1) && (indices[3] == 1))
          base_index[1] = 1;
        else
          base_index[1] = 2;

        return data[base_index[0]][base_index[1]];
      }

      case 3:
        // each entry of the tensor can be
        // thought of as an entry in a
        // matrix that maps the rolled-out
        // rank-2 tensors into rolled-out
        // rank-2 tensors. this is the
        // format in which we store rank-4
        // tensors. determine which
        // position the present entry is
        // stored in
      {
        unsigned int base_index[2] ;
        if ((indices[0] == 0) && (indices[1] == 0))
          base_index[0] = 0;
        else if ((indices[0] == 1) && (indices[1] == 1))
          base_index[0] = 1;
        else if ((indices[0] == 2) && (indices[1] == 2))
          base_index[0] = 2;
        else if (((indices[0] == 0) && (indices[1] == 1)) ||
                 ((indices[0] == 1) && (indices[1] == 0)))
          base_index[0] = 3;
        else if (((indices[0] == 0) && (indices[1] == 2)) ||
                 ((indices[0] == 2) && (indices[1] == 0)))
          base_index[0] = 4;
        else
          {
            Assert (((indices[0] == 1) && (indices[1] == 2)) ||
                    ((indices[0] == 2) && (indices[1] == 1)),
                    ExcInternalError());
            base_index[0] = 5;
          }

        if ((indices[2] == 0) && (indices[3] == 0))
          base_index[1] = 0;
        else if ((indices[2] == 1) && (indices[3] == 1))
          base_index[1] = 1;
        else if ((indices[2] == 2) && (indices[3] == 2))
          base_index[1] = 2;
        else if (((indices[2] == 0) && (indices[3] == 1)) ||
                 ((indices[2] == 1) && (indices[3] == 0)))
          base_index[1] = 3;
        else if (((indices[2] == 0) && (indices[3] == 2)) ||
                 ((indices[2] == 2) && (indices[3] == 0)))
          base_index[1] = 4;
        else
          {
            Assert (((indices[2] == 1) && (indices[3] == 2)) ||
                    ((indices[2] == 2) && (indices[3] == 1)),
                    ExcInternalError());
            base_index[1] = 5;
          }

        return data[base_index[0]][base_index[1]];
      }

      default:
        Assert (false, ExcNotImplemented());
      }

    static Number dummy;
    return dummy;
  }

} // end of namespace internal



template <int rank, int dim, typename Number>
inline
Number &
SymmetricTensor<rank,dim,Number>::operator () (const TableIndices<rank> &indices)
{
  for (unsigned int r=0; r<rank; ++r)
    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
  return internal::symmetric_tensor_access<dim,Number> (indices, data);
}



template <int rank, int dim, typename Number>
inline
Number
SymmetricTensor<rank,dim,Number>::operator ()
(const TableIndices<rank> &indices) const
{
  for (unsigned int r=0; r<rank; ++r)
    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
  return internal::symmetric_tensor_access<dim,Number> (indices, data);
}



template <int rank, int dim, typename Number>
internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
{
  return
    internal::SymmetricTensorAccessors::
    Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
}



template <int rank, int dim, typename Number>
internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
{
  return
    internal::SymmetricTensorAccessors::
    Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
}



template <int rank, int dim, typename Number>
inline
Number
SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices) const
{
  return operator()(indices);
}



template <int rank, int dim, typename Number>
inline
Number &
SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices)
{
  return operator()(indices);
}




namespace internal
{
  namespace SymmetricTensor
  {
    template <int dim, typename Number>
    unsigned int
    entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
                      const unsigned int index)
    {
      return index;
    }


    template <int dim, typename Number>
    ::TableIndices<2>
    entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
                      const unsigned int index)
    {
      return
        internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type::
        unrolled_to_component_indices(index);
    }

  }
}



template <int rank, int dim, typename Number>
inline
Number
SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
{
  AssertIndexRange (index, n_independent_components);
  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
}



template <int rank, int dim, typename Number>
inline
Number &
SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
{
  AssertIndexRange (index, n_independent_components);
  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
}



namespace internal
{
  template <int dim, typename Number>
  inline
  typename numbers::NumberTraits<Number>::real_type
  compute_norm (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
  {
    switch (dim)
      {
      case 1:
        return numbers::NumberTraits<Number>::abs(data[0]);

      case 2:
        return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
                         numbers::NumberTraits<Number>::abs_square(data[1]) +
                         2. * numbers::NumberTraits<Number>::abs_square(data[2]));

      case 3:
        return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
                         numbers::NumberTraits<Number>::abs_square(data[1]) +
                         numbers::NumberTraits<Number>::abs_square(data[2]) +
                         2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
                         2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
                         2. * numbers::NumberTraits<Number>::abs_square(data[5]));

      default:
      {
        typename numbers::NumberTraits<Number>::real_type return_value
          = typename numbers::NumberTraits<Number>::real_type();

        for (unsigned int d=0; d<dim; ++d)
          return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
        for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
          return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);

        return std::sqrt(return_value);
      }
      }
  }



  template <int dim, typename Number>
  inline
  typename numbers::NumberTraits<Number>::real_type
  compute_norm (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
  {
    switch (dim)
      {
      case 1:
        return numbers::NumberTraits<Number>::abs (data[0][0]);

      default:
      {
        typename numbers::NumberTraits<Number>::real_type return_value
          = typename numbers::NumberTraits<Number>::real_type();

        const unsigned int n_independent_components = data.dimension;

        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=0; j<dim; ++j)
            return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=dim; j<n_independent_components; ++j)
            return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
        for (unsigned int i=dim; i<n_independent_components; ++i)
          for (unsigned int j=0; j<dim; ++j)
            return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
        for (unsigned int i=dim; i<n_independent_components; ++i)
          for (unsigned int j=dim; j<n_independent_components; ++j)
            return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);

        return std::sqrt(return_value);
      }
      }
  }

} // end of namespace internal



template <int rank, int dim, typename Number>
inline
typename numbers::NumberTraits<Number>::real_type
SymmetricTensor<rank,dim,Number>::norm () const
{
  return internal::compute_norm<dim,Number> (data);
}



namespace internal
{
  namespace SymmetricTensor
  {
    namespace
    {
      // a function to do the unrolling from a set of indices to a
      // scalar index into the array in which we store the elements of
      // a symmetric tensor
      //
      // this function is for rank-2 tensors
      template <int dim>
      inline
      unsigned int
      component_to_unrolled_index
      (const TableIndices<2> &indices)
      {
        Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
        Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));

        switch (dim)
          {
          case 1:
          {
            return 0;
          }

          case 2:
          {
            static const unsigned int table[2][2] = {{0, 2},
              {2, 1}
            };
            return table[indices[0]][indices[1]];
          }

          case 3:
          {
            static const unsigned int table[3][3] = {{0, 3, 4},
              {3, 1, 5},
              {4, 5, 2}
            };
            return table[indices[0]][indices[1]];
          }

          case 4:
          {
            static const unsigned int table[4][4] = {{0, 4, 5, 6},
              {4, 1, 7, 8},
              {5, 7, 2, 9},
              {6, 8, 9, 3}
            };
            return table[indices[0]][indices[1]];
          }

          default:
            // for the remainder, manually figure out the numbering
          {
            if (indices[0] == indices[1])
              return indices[0];

            TableIndices<2> sorted_indices (indices);
            sorted_indices.sort ();

            for (unsigned int d=0, c=0; d<dim; ++d)
              for (unsigned int e=d+1; e<dim; ++e, ++c)
                if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
                  return dim+c;

            // should never get here:
            Assert(false, ExcInternalError());
            return 0;
          }
          }
      }

      // a function to do the unrolling from a set of indices to a
      // scalar index into the array in which we store the elements of
      // a symmetric tensor
      //
      // this function is for tensors of ranks not already handled
      // above
      template <int dim, int rank>
      inline
      unsigned int
      component_to_unrolled_index
      (const TableIndices<rank> &indices)
      {
        (void)indices;
        Assert (false, ExcNotImplemented());
        return numbers::invalid_unsigned_int;
      }
    }
  }
}


template <int rank, int dim, typename Number>
inline
unsigned int
SymmetricTensor<rank,dim,Number>::component_to_unrolled_index
(const TableIndices<rank> &indices)
{
  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
}



namespace internal
{
  namespace SymmetricTensor
  {
    namespace
    {
      // a function to do the inverse of the unrolling from a set of
      // indices to a scalar index into the array in which we store
      // the elements of a symmetric tensor. in other words, it goes
      // from the scalar index into the array to a set of indices of
      // the tensor
      //
      // this function is for rank-2 tensors
      template <int dim>
      inline
      TableIndices<2>
      unrolled_to_component_indices
      (const unsigned int i,
       const int2type<2> &)
      {
        Assert ((i < ::SymmetricTensor<2,dim,double>::n_independent_components),
                ExcIndexRange(i, 0, ::SymmetricTensor<2,dim,double>::n_independent_components));
        switch (dim)
          {
          case 1:
          {
            return TableIndices<2>(0,0);
          }

          case 2:
          {
            const TableIndices<2> table[3] =
            {
              TableIndices<2> (0,0),
              TableIndices<2> (1,1),
              TableIndices<2> (0,1)
            };
            return table[i];
          }

          case 3:
          {
            const TableIndices<2> table[6] =
            {
              TableIndices<2> (0,0),
              TableIndices<2> (1,1),
              TableIndices<2> (2,2),
              TableIndices<2> (0,1),
              TableIndices<2> (0,2),
              TableIndices<2> (1,2)
            };
            return table[i];
          }

          default:
            if (i<dim)
              return TableIndices<2> (i,i);

            for (unsigned int d=0, c=0; d<dim; ++d)
              for (unsigned int e=d+1; e<dim; ++e, ++c)
                if (c==i)
                  return TableIndices<2>(d,e);

            // should never get here:
            Assert(false, ExcInternalError());
            return TableIndices<2>(0, 0);
          }
      }

      // a function to do the inverse of the unrolling from a set of
      // indices to a scalar index into the array in which we store
      // the elements of a symmetric tensor. in other words, it goes
      // from the scalar index into the array to a set of indices of
      // the tensor
      //
      // this function is for tensors of a rank not already handled
      // above
      template <int dim, int rank>
      inline
      TableIndices<rank>
      unrolled_to_component_indices
      (const unsigned int i,
       const int2type<rank> &)
      {
        (void)i;
        Assert ((i < ::SymmetricTensor<rank,dim,double>::n_independent_components),
                ExcIndexRange(i, 0, ::SymmetricTensor<rank,dim,double>::n_independent_components));
        Assert (false, ExcNotImplemented());
        return TableIndices<rank>();
      }

    }
  }
}

template <int rank, int dim, typename Number>
inline
TableIndices<rank>
SymmetricTensor<rank,dim,Number>::unrolled_to_component_indices
(const unsigned int i)
{
  return
    internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
        internal::int2type<rank>());
}



template <int rank, int dim, typename Number>
template <class Archive>
inline
void
SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
{
  ar &data;
}


#endif // DOXYGEN

/* ----------------- Non-member functions operating on tensors. ------------ */


template <int rank, int dim, typename Number, typename OtherNumber>
inline
Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
operator+(const SymmetricTensor<rank, dim, Number> &left,
          const Tensor<rank, dim, OtherNumber> &right)
{
  return Tensor<rank, dim, Number>(left) + right;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline
Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
operator+(const Tensor<rank, dim, Number> &left,
          const SymmetricTensor<rank, dim, OtherNumber> &right)
{
  return left + Tensor<rank, dim, OtherNumber>(right);
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline
Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
operator-(const SymmetricTensor<rank, dim, Number> &left,
          const Tensor<rank, dim, OtherNumber> &right)
{
  return Tensor<rank, dim, Number>(left) - right;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline
Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
operator-(const Tensor<rank, dim, Number> &left,
          const SymmetricTensor<rank, dim, OtherNumber> &right)
{
  return left - Tensor<rank, dim, OtherNumber>(right);
}



template <int dim, typename Number>
inline
Number determinant (const SymmetricTensor<2,dim,Number> &t)
{
  switch (dim)
    {
    case 1:
      return t.data[0];
    case 2:
      return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
    case 3:
    {
      // in analogy to general tensors, but
      // there's something to be simplified for
      // the present case
      const Number tmp = t.data[3]*t.data[4]*t.data[5];
      return ( tmp + tmp
               +t.data[0]*t.data[1]*t.data[2]
               -t.data[0]*t.data[5]*t.data[5]
               -t.data[1]*t.data[4]*t.data[4]
               -t.data[2]*t.data[3]*t.data[3]);
    }
    default:
      Assert (false, ExcNotImplemented());
      return 0;
    }
}



template <int dim, typename Number>
inline
Number third_invariant (const SymmetricTensor<2,dim,Number> &t)
{
  return determinant (t);
}



template <int dim, typename Number>
Number trace (const SymmetricTensor<2,dim,Number> &d)
{
  Number t = d.data[0];
  for (unsigned int i=1; i<dim; ++i)
    t += d.data[i];
  return t;
}


template <int dim, typename Number>
inline
Number first_invariant (const SymmetricTensor<2,dim,Number> &t)
{
  return trace (t);
}


template <typename Number>
inline
Number second_invariant (const SymmetricTensor<2,1,Number> &)
{
  return 0;
}



template <typename Number>
inline
Number second_invariant (const SymmetricTensor<2,2,Number> &t)
{
  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
}



template <typename Number>
inline
Number second_invariant (const SymmetricTensor<2,3,Number> &t)
{
  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
          - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
}




template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
transpose (const SymmetricTensor<rank,dim,Number> &t)
{
  return t;
}



template <int dim, typename Number>
inline
SymmetricTensor<2,dim,Number>
deviator (const SymmetricTensor<2,dim,Number> &t)
{
  SymmetricTensor<2,dim,Number> tmp = t;

  // subtract scaled trace from the diagonal
  const Number tr = trace(t) / dim;
  for (unsigned int i=0; i<dim; ++i)
    tmp.data[i] -= tr;

  return tmp;
}



template <int dim, typename Number>
inline
SymmetricTensor<2,dim,Number>
unit_symmetric_tensor ()
{
  // create a default constructed matrix filled with
  // zeros, then set the diagonal elements to one
  SymmetricTensor<2,dim,Number> tmp;
  switch (dim)
    {
    case 1:
      tmp.data[0] = 1;
      break;
    case 2:
      tmp.data[0] = tmp.data[1] = 1;
      break;
    case 3:
      tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
      break;
    default:
      for (unsigned int d=0; d<dim; ++d)
        tmp.data[d] = 1;
    }
  return tmp;
}



template <int dim>
inline
SymmetricTensor<2,dim>
unit_symmetric_tensor ()
{
  return unit_symmetric_tensor<dim,double>();
}



template <int dim, typename Number>
inline
SymmetricTensor<4,dim,Number>
deviator_tensor ()
{
  SymmetricTensor<4,dim,Number> tmp;

  // fill the elements treating the diagonal
  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=0; j<dim; ++j)
      tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;

  // then fill the ones that copy over the
  // non-diagonal elements. note that during
  // the double-contraction, we handle the
  // off-diagonal elements twice, so simply
  // copying requires a weight of 1/2
  for (unsigned int i=dim;
       i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
       ++i)
    tmp.data[i][i] = 0.5;

  return tmp;
}



template <int dim>
inline
SymmetricTensor<4,dim>
deviator_tensor ()
{
  return deviator_tensor<dim,double>();
}



template <int dim, typename Number>
inline
SymmetricTensor<4,dim,Number>
identity_tensor ()
{
  SymmetricTensor<4,dim,Number> tmp;

  // fill the elements treating the diagonal
  for (unsigned int i=0; i<dim; ++i)
    tmp.data[i][i] = 1;

  // then fill the ones that copy over the
  // non-diagonal elements. note that during
  // the double-contraction, we handle the
  // off-diagonal elements twice, so simply
  // copying requires a weight of 1/2
  for (unsigned int i=dim;
       i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
       ++i)
    tmp.data[i][i] = 0.5;

  return tmp;
}



template <int dim>
inline
SymmetricTensor<4,dim>
identity_tensor ()
{
  return identity_tensor<dim,double>();
}



template <int dim, typename Number>
inline
SymmetricTensor<2,dim,Number>
invert (const SymmetricTensor<2,dim,Number> &)
{
  // if desired, take over the
  // inversion of a 4x4 tensor
  // from the FullMatrix
  AssertThrow (false, ExcNotImplemented());

  return SymmetricTensor<2,dim,Number>();
}



#ifndef DOXYGEN

template <typename Number>
inline
SymmetricTensor<2,1,Number>
invert (const SymmetricTensor<2,1,Number> &t)
{
  SymmetricTensor<2,1,Number> tmp;

  tmp[0][0] = 1.0/t[0][0];

  return tmp;
}



template <typename Number>
inline
SymmetricTensor<2,2,Number>
invert (const SymmetricTensor<2,2,Number> &t)
{
  SymmetricTensor<2,2,Number> tmp;

  // Sympy result: ([
  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
  // [-t01/(t00*t11 - t01**2),  t00/(t00*t11 - t01**2)]  ])
  const TableIndices<2> idx_00 (0,0);
  const TableIndices<2> idx_01 (0,1);
  const TableIndices<2> idx_11 (1,1);
  const Number inv_det_t
    = 1.0/(t[idx_00]*t[idx_11]
           - t[idx_01]*t[idx_01]);
  tmp[idx_00] = t[idx_11];
  tmp[idx_01] = -t[idx_01];
  tmp[idx_11] = t[idx_00];
  tmp *= inv_det_t;

  return tmp;
}



template <typename Number>
inline
SymmetricTensor<2,3,Number>
invert (const SymmetricTensor<2,3,Number> &t)
{
  SymmetricTensor<2,3,Number> tmp;

  // Sympy result: ([
  // [  (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
  //    (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
  //    (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
  // [  (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
  //    (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
  //    (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
  // [  (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
  //    (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
  //    (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)]  ])
  const TableIndices<2> idx_00 (0,0);
  const TableIndices<2> idx_01 (0,1);
  const TableIndices<2> idx_02 (0,2);
  const TableIndices<2> idx_11 (1,1);
  const TableIndices<2> idx_12 (1,2);
  const TableIndices<2> idx_22 (2,2);
  const Number inv_det_t
    = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
           - t[idx_00]*t[idx_12]*t[idx_12]
           - t[idx_01]*t[idx_01]*t[idx_22]
           + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
           - t[idx_02]*t[idx_02]*t[idx_11]);
  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
  tmp *= inv_det_t;

  return tmp;
}

#endif /* DOXYGEN */



template <int dim, typename Number>
inline
SymmetricTensor<4,dim,Number>
invert (const SymmetricTensor<4,dim,Number> &t)
{
  SymmetricTensor<4,dim,Number> tmp;
  switch (dim)
    {
    case 1:
      tmp.data[0][0] = 1./t.data[0][0];
      break;
    case 2:

      // inverting this tensor is a little more
      // complicated than necessary, since we
      // store the data of 't' as a 3x3 matrix
      // t.data, but the product between a rank-4
      // and a rank-2 tensor is really not the
      // product between this matrix and the
      // 3-vector of a rhs, but rather
      //
      // B.vec = t.data * mult * A.vec
      //
      // where mult is a 3x3 matrix with
      // entries [[1,0,0],[0,1,0],[0,0,2]] to
      // capture the fact that we need to add up
      // both the c_ij12*a_12 and the c_ij21*a_21
      // terms
      //
      // in addition, in this scheme, the
      // identity tensor has the matrix
      // representation mult^-1.
      //
      // the inverse of 't' therefore has the
      // matrix representation
      //
      // inv.data = mult^-1 * t.data^-1 * mult^-1
      //
      // in order to compute it, let's first
      // compute the inverse of t.data and put it
      // into tmp.data; at the end of the
      // function we then scale the last row and
      // column of the inverse by 1/2,
      // corresponding to the left and right
      // multiplication with mult^-1
    {
      const Number t4 = t.data[0][0]*t.data[1][1],
                   t6 = t.data[0][0]*t.data[1][2],
                   t8 = t.data[0][1]*t.data[1][0],
                   t00 = t.data[0][2]*t.data[1][0],
                   t01 = t.data[0][1]*t.data[2][0],
                   t04 = t.data[0][2]*t.data[2][0],
                   t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
                              t8*t.data[2][2]+t00*t.data[2][1]+
                              t01*t.data[1][2]-t04*t.data[1][1]);
      tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
      tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
      tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
      tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
      tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
      tmp.data[1][2] = -(t6-t00)*t07;
      tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
      tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
      tmp.data[2][2] = (t4-t8)*t07;

      // scale last row and column as mentioned
      // above
      tmp.data[2][0] /= 2;
      tmp.data[2][1] /= 2;
      tmp.data[0][2] /= 2;
      tmp.data[1][2] /= 2;
      tmp.data[2][2] /= 4;
    }
    break;
    default:
      Assert (false, ExcNotImplemented());
    }
  return tmp;
}



template <>
SymmetricTensor<4,3,double>
invert (const SymmetricTensor<4,3,double> &t);
// this function is implemented in the .cc file for double data types



template <int dim, typename Number>
inline
SymmetricTensor<4,dim,Number>
outer_product (const SymmetricTensor<2,dim,Number> &t1,
               const SymmetricTensor<2,dim,Number> &t2)
{
  SymmetricTensor<4,dim,Number> tmp;

  // fill only the elements really needed
  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=i; j<dim; ++j)
      for (unsigned int k=0; k<dim; ++k)
        for (unsigned int l=k; l<dim; ++l)
          tmp[i][j][k][l] = t1[i][j] * t2[k][l];

  return tmp;
}



template <int dim,typename Number>
inline
SymmetricTensor<2,dim,Number>
symmetrize (const Tensor<2,dim,Number> &t)
{
  Number array[(dim*dim+dim)/2];
  for (unsigned int d=0; d<dim; ++d)
    array[d] = t[d][d];
  for (unsigned int d=0, c=0; d<dim; ++d)
    for (unsigned int e=d+1; e<dim; ++e, ++c)
      array[dim+c] = (t[d][e]+t[e][d])*0.5;
  return SymmetricTensor<2,dim,Number>(array);
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
operator * (const SymmetricTensor<rank,dim,Number> &t,
            const Number                            factor)
{
  SymmetricTensor<rank,dim,Number> tt = t;
  tt *= factor;
  return tt;
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
operator * (const Number                            factor,
            const SymmetricTensor<rank,dim,Number> &t)
{
  // simply forward to the other operator
  return t*factor;
}


#ifndef DEAL_II_WITH_CXX11

template <typename T, typename U, int rank, int dim>
struct ProductType<T,SymmetricTensor<rank,dim,U> >
{
  typedef SymmetricTensor<rank,dim,typename ProductType<T,U>::type> type;
};

template <typename T, typename U, int rank, int dim>
struct ProductType<SymmetricTensor<rank,dim,T>,U>
{
  typedef SymmetricTensor<rank,dim,typename ProductType<T,U>::type> type;
};

#endif



template <int rank, int dim, typename Number, typename OtherNumber>
inline
SymmetricTensor<rank,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
operator * (const SymmetricTensor<rank,dim,Number> &t,
            const OtherNumber                    factor)
{
  // form the product. we have to convert the two factors into the final
  // type via explicit casts because, for awkward reasons, the C++
  // standard committee saw it fit to not define an
  //   operator*(float,std::complex<double>)
  // (as well as with switched arguments and double<->float).
  typedef typename ProductType<Number,OtherNumber>::type product_type;
  SymmetricTensor<rank,dim,product_type> tt(t);
  // we used to shorten the following by 'tt *= product_type(factor);'
  // which requires that a converting constructor
  // 'product_type::product_type(const OtherNumber) is defined.
  // however, a user-defined constructor is not allowed for aggregates,
  // e.g. VectorizedArray. therefore, we work around this issue using a
  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
  // which we assume to be defined.
  product_type new_factor;
  new_factor = factor;
  tt *= new_factor;
  return tt;
}



template <int rank, int dim, typename Number, typename OtherNumber>
inline
SymmetricTensor<rank,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
operator * (const Number                     factor,
            const SymmetricTensor<rank,dim,OtherNumber> &t)
{
  // simply forward to the other operator with switched arguments
  return (t*factor);
}



template <int rank, int dim, typename Number>
inline
SymmetricTensor<rank,dim,Number>
operator / (const SymmetricTensor<rank,dim,Number> &t,
            const Number                            factor)
{
  SymmetricTensor<rank,dim,Number> tt = t;
  tt /= factor;
  return tt;
}



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;
}



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;
}



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;
}

template <int dim, typename Number>
inline
Number
scalar_product (const SymmetricTensor<2,dim,Number> &t1,
                const SymmetricTensor<2,dim,Number> &t2)
{
  return (t1*t2);
}


template <int dim, typename Number>
inline
Number
scalar_product (const SymmetricTensor<2,dim,Number> &t1,
                const Tensor<2,dim,Number> &t2)
{
  Number s = 0;
  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=0; j<dim; ++j)
      s += t1[i][j] * t2[i][j];
  return s;
}


template <int dim, typename Number>
inline
Number
scalar_product (const Tensor<2,dim,Number> &t1,
                const SymmetricTensor<2,dim,Number> &t2)
{
  return scalar_product(t2, t1);
}


template <typename Number>
inline
void
double_contract (SymmetricTensor<2,1,Number> &tmp,
                 const SymmetricTensor<4,1,Number> &t,
                 const SymmetricTensor<2,1,Number> &s)
{
  tmp[0][0] = t[0][0][0][0] * s[0][0];
}



template <typename Number>
inline
void
double_contract (SymmetricTensor<2,1,Number> &tmp,
                 const SymmetricTensor<2,1,Number> &s,
                 const SymmetricTensor<4,1,Number> &t)
{
  tmp[0][0] = t[0][0][0][0] * s[0][0];
}



template <typename Number>
inline
void
double_contract (SymmetricTensor<2,2,Number> &tmp,
                 const SymmetricTensor<4,2,Number> &t,
                 const SymmetricTensor<2,2,Number> &s)
{
  const unsigned int dim = 2;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=i; j<dim; ++j)
      tmp[i][j] = t[i][j][0][0] * s[0][0] +
                  t[i][j][1][1] * s[1][1] +
                  2 * t[i][j][0][1] * s[0][1];
}



template <typename Number>
inline
void
double_contract (SymmetricTensor<2,2,Number> &tmp,
                 const SymmetricTensor<2,2,Number> &s,
                 const SymmetricTensor<4,2,Number> &t)
{
  const unsigned int dim = 2;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=i; j<dim; ++j)
      tmp[i][j] = s[0][0] * t[0][0][i][j] * +
                  s[1][1] * t[1][1][i][j] +
                  2 * s[0][1] * t[0][1][i][j];
}



template <typename Number>
inline
void
double_contract (SymmetricTensor<2,3,Number> &tmp,
                 const SymmetricTensor<4,3,Number> &t,
                 const SymmetricTensor<2,3,Number> &s)
{
  const unsigned int dim = 3;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=i; j<dim; ++j)
      tmp[i][j] = t[i][j][0][0] * s[0][0] +
                  t[i][j][1][1] * s[1][1] +
                  t[i][j][2][2] * s[2][2] +
                  2 * t[i][j][0][1] * s[0][1] +
                  2 * t[i][j][0][2] * s[0][2] +
                  2 * t[i][j][1][2] * s[1][2];
}



template <typename Number>
inline
void
double_contract (SymmetricTensor<2,3,Number> &tmp,
                 const SymmetricTensor<2,3,Number> &s,
                 const SymmetricTensor<4,3,Number> &t)
{
  const unsigned int dim = 3;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=i; j<dim; ++j)
      tmp[i][j] = s[0][0] * t[0][0][i][j] +
                  s[1][1] * t[1][1][i][j] +
                  s[2][2] * t[2][2][i][j] +
                  2 * s[0][1] * t[0][1][i][j] +
                  2 * s[0][2] * t[0][2][i][j] +
                  2 * s[1][2] * t[1][2][i][j];
}



template <int dim, typename Number>
Tensor<1,dim,Number>
operator * (const SymmetricTensor<2,dim,Number> &src1,
            const Tensor<1,dim,Number> &src2)
{
  Tensor<1,dim,Number> dest;
  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=0; j<dim; ++j)
      dest[i] += src1[i][j] * src2[j];
  return dest;
}


template <int dim, typename Number>
Tensor<1,dim,Number>
operator * (const Tensor<1,dim,Number> &src1,
            const SymmetricTensor<2,dim,Number> &src2)
{
  // this is easy for symmetric tensors:
  return src2 * src1;
}


template <int dim, typename Number>
inline
std::ostream &operator << (std::ostream &out,
                           const SymmetricTensor<2,dim,Number> &t)
{
  //make out lives a bit simpler by outputing
  //the tensor through the operator for the
  //general Tensor class
  Tensor<2,dim,Number> tt;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=0; j<dim; ++j)
      tt[i][j] = t[i][j];

  return out << tt;
}



template <int dim, typename Number>
inline
std::ostream &operator << (std::ostream &out,
                           const SymmetricTensor<4,dim,Number> &t)
{
  //make out lives a bit simpler by outputing
  //the tensor through the operator for the
  //general Tensor class
  Tensor<4,dim,Number> tt;

  for (unsigned int i=0; i<dim; ++i)
    for (unsigned int j=0; j<dim; ++j)
      for (unsigned int k=0; k<dim; ++k)
        for (unsigned int l=0; l<dim; ++l)
          tt[i][j][k][l] = t[i][j][k][l];

  return out << tt;
}


DEAL_II_NAMESPACE_CLOSE

#endif