symmetric_tensor.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2021 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_symmetric_tensor_h
#define dealii_symmetric_tensor_h
 
 
#include <deal.II/base/config.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/tensor.h>
 
#include <algorithm>
#include <array>
#include <functional>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <int rank, int dim, typename Number = double>
class SymmetricTensor;
#endif
 
template <int dim, typename Number>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
                                               unit_symmetric_tensor();
 
template <int dim, typename Number>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
                                               deviator_tensor();
 
template <int dim, typename Number>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
                                               identity_tensor();
 
template <int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
                                invert(const SymmetricTensor<2, dim, Number> &);
 
template <int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
                                invert(const SymmetricTensor<4, dim, Number> &);
 
template <int dim2, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE Number
                                       trace(const SymmetricTensor<2, dim2, Number> &);
 
template <int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
                                       deviator(const SymmetricTensor<2, dim, Number> &);
 
template <int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE Number
                                       determinant(const SymmetricTensor<2, dim, Number> &);
 
 
 
namespace internal
{
  // Workaround: The following 4 overloads are necessary to be able to
  // compile the library with Apple Clang 8 and older. We should remove
  // these overloads again when we bump the minimal required version to
  // something later than clang-3.6 / Apple Clang 6.3.
  template <int rank, int dim, typename T, typename U>
  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
  {
    using type =
      SymmetricTensor<rank,
                      dim,
                      std::complex<typename ProductType<T, U>::type>>;
  };
 
  template <int rank, int dim, typename T, typename U>
  struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
                         std::complex<U>>
  {
    using type =
      SymmetricTensor<rank,
                      dim,
                      std::complex<typename ProductType<T, U>::type>>;
  };
 
  template <typename T, int rank, int dim, typename U>
  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
  {
    using type =
      SymmetricTensor<rank,
                      dim,
                      std::complex<typename ProductType<T, U>::type>>;
  };
 
  template <int rank, int dim, typename T, typename U>
  struct ProductTypeImpl<std::complex<T>,
                         SymmetricTensor<rank, dim, std::complex<U>>>
  {
    using type =
      SymmetricTensor<rank,
                      dim,
                      std::complex<typename ProductType<T, U>::type>>;
  };
  // end workaround
 
  namespace SymmetricTensorImplementation
  {
    template <int rank, int dim, typename Number>
    struct Inverse;
  } // namespace SymmetricTensorImplementation
 
  namespace SymmetricTensorAccessors
  {
    DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
                                                   merge(const TableIndices<2> &previous_indices,
                                                         const unsigned int     new_index,
                                                         const unsigned int     position)
    {
      AssertIndexRange(position, 2);
 
      if (position == 0)
        return {new_index, numbers::invalid_unsigned_int};
      else
        return {previous_indices[0], new_index};
    }
 
 
 
    DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<4>
                                                   merge(const TableIndices<4> &previous_indices,
                                                         const unsigned int     new_index,
                                                         const unsigned int     position)
    {
      AssertIndexRange(position, 4);
 
      switch (position)
        {
          case 0:
            return {new_index,
                    numbers::invalid_unsigned_int,
                    numbers::invalid_unsigned_int,
                    numbers::invalid_unsigned_int};
          case 1:
            return {previous_indices[0],
                    new_index,
                    numbers::invalid_unsigned_int,
                    numbers::invalid_unsigned_int};
          case 2:
            return {previous_indices[0],
                    previous_indices[1],
                    new_index,
                    numbers::invalid_unsigned_int};
          case 3:
            return {previous_indices[0],
                    previous_indices[1],
                    previous_indices[2],
                    new_index};
          default:
            Assert(false, ExcInternalError());
            return {};
        }
    }
 
 
    template <int rank1,
              int rank2,
              int dim,
              typename Number,
              typename OtherNumber = Number>
    struct double_contraction_result
    {
      using value_type = typename ProductType<Number, OtherNumber>::type;
      using type =
        ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
    };
 
 
    template <int dim, typename Number, typename OtherNumber>
    struct double_contraction_result<2, 2, dim, Number, OtherNumber>
    {
      using type = typename ProductType<Number, OtherNumber>::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;
 
      using base_tensor_type = Tensor<1, n_independent_components, Number>;
    };
 
 
 
    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);
 
      using base_tensor_type = Tensor<2, n_rank2_components, Number>;
    };
 
 
 
    template <int rank, int dim, bool constness, typename Number>
    struct AccessorTypes;
 
    template <int rank, int dim, typename Number>
    struct AccessorTypes<rank, dim, true, Number>
    {
      using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
 
      using reference = Number;
    };
 
    template <int rank, int dim, typename Number>
    struct AccessorTypes<rank, dim, false, Number>
    {
      using tensor_type = ::SymmetricTensor<rank, dim, Number>;
 
      using reference = Number &;
    };
 
 
    template <int rank, int dim, bool constness, int P, typename Number>
    class Accessor
    {
    public:
      using reference =
        typename AccessorTypes<rank, dim, constness, Number>::reference;
      using tensor_type =
        typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
 
    private:
      constexpr Accessor(tensor_type &             tensor,
                         const TableIndices<rank> &previous_indices);
 
      constexpr DEAL_II_ALWAYS_INLINE
      Accessor(const Accessor &) = default;
 
    public:
      constexpr Accessor<rank, dim, constness, P - 1, Number>
      operator[](const unsigned int i);
 
      constexpr 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;
      friend class ::SymmetricTensor<rank, dim, Number>;
      friend class Accessor<rank, dim, constness, P + 1, Number>;
    };
 
 
 
    template <int rank, int dim, bool constness, typename Number>
    class Accessor<rank, dim, constness, 1, Number>
    {
    public:
      using reference =
        typename AccessorTypes<rank, dim, constness, Number>::reference;
      using tensor_type =
        typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
 
    private:
      constexpr Accessor(tensor_type &             tensor,
                         const TableIndices<rank> &previous_indices);
 
      constexpr DEAL_II_ALWAYS_INLINE
      Accessor(const Accessor &) = default;
 
    public:
      constexpr reference operator[](const unsigned int);
 
      constexpr 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;
      friend class ::SymmetricTensor<rank, dim, Number>;
      friend class SymmetricTensorAccessors::
        Accessor<rank, dim, constness, 2, Number>;
    };
  } // namespace SymmetricTensorAccessors
} // namespace internal
 
 
 
template <int rank_, int dim, typename Number>
class SymmetricTensor
{
public:
  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
 
  static const unsigned int dimension = dim;
 
  static const unsigned int rank = rank_;
 
  static constexpr unsigned int n_independent_components =
    internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>::
      n_independent_components;
 
  constexpr DEAL_II_ALWAYS_INLINE
  SymmetricTensor() = default;
 
  template <typename OtherNumber>
  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
 
  constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
 
  template <typename OtherNumber>
  constexpr explicit SymmetricTensor(
    const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
 
  Number *
  begin_raw();
 
  const Number *
  begin_raw() const;
 
  Number *
  end_raw();
 
  const Number *
  end_raw() const;
 
  template <typename OtherNumber>
  constexpr SymmetricTensor &
  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
 
  constexpr SymmetricTensor &
  operator=(const Number &d);
 
  constexpr operator Tensor<rank_, dim, Number>() const;
 
  constexpr bool
  operator==(const SymmetricTensor &) const;
 
  constexpr bool
  operator!=(const SymmetricTensor &) const;
 
  template <typename OtherNumber>
  constexpr SymmetricTensor &
  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
 
  template <typename OtherNumber>
  constexpr SymmetricTensor &
  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
 
  template <typename OtherNumber>
  constexpr SymmetricTensor &
  operator*=(const OtherNumber &factor);
 
  template <typename OtherNumber>
  constexpr SymmetricTensor &
  operator/=(const OtherNumber &factor);
 
  constexpr SymmetricTensor
  operator-() const;
 
  template <typename OtherNumber>
  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
    operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const;
 
  template <typename OtherNumber>
  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
    operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const;
 
  constexpr Number &
  operator()(const TableIndices<rank_> &indices);
 
  constexpr const Number &
  operator()(const TableIndices<rank_> &indices) const;
 
  constexpr internal::SymmetricTensorAccessors::
    Accessor<rank_, dim, true, rank_ - 1, Number>
    operator[](const unsigned int row) const;
 
  constexpr internal::SymmetricTensorAccessors::
    Accessor<rank_, dim, false, rank_ - 1, Number>
    operator[](const unsigned int row);
 
  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
 
  constexpr Number &operator[](const TableIndices<rank_> &indices);
 
  constexpr const Number &
  access_raw_entry(const unsigned int unrolled_index) const;
 
  constexpr Number &
  access_raw_entry(const unsigned int unrolled_index);
 
  constexpr typename numbers::NumberTraits<Number>::real_type
  norm() const;
 
  static constexpr unsigned int
  component_to_unrolled_index(const TableIndices<rank_> &indices);
 
  static constexpr TableIndices<rank_>
  unrolled_to_component_indices(const unsigned int i);
 
  constexpr void
  clear();
 
  static constexpr std::size_t
  memory_consumption();
 
  template <class Archive>
  void
  serialize(Archive &ar, const unsigned int version);
 
private:
  using base_tensor_descriptor =
    internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>;
 
  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
 
  base_tensor_type data;
 
  // Make all other symmetric tensors friends.
  template <int, int, typename>
  friend class SymmetricTensor;
 
  // Make a few more functions friends.
  template <int dim2, typename Number2>
  friend constexpr Number2
  trace(const SymmetricTensor<2, dim2, Number2> &d);
 
  template <int dim2, typename Number2>
  friend constexpr Number2
  determinant(const SymmetricTensor<2, dim2, Number2> &t);
 
  template <int dim2, typename Number2>
  friend constexpr SymmetricTensor<2, dim2, Number2>
  deviator(const SymmetricTensor<2, dim2, Number2> &t);
 
  template <int dim2, typename Number2>
  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
                           unit_symmetric_tensor();
 
  template <int dim2, typename Number2>
  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
                           deviator_tensor();
 
  template <int dim2, typename Number2>
  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
                           identity_tensor();
 
 
  // Make a few helper classes friends as well.
  friend struct internal::SymmetricTensorImplementation::
    Inverse<2, dim, Number>;
 
  friend struct internal::SymmetricTensorImplementation::
    Inverse<4, dim, Number>;
};
 
 
 
// ------------------------- inline functions ------------------------
 
#ifndef DOXYGEN
 
// provide declarations for static members
template <int rank, int dim, typename Number>
const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
 
template <int rank_, int dim, typename Number>
constexpr unsigned int
  SymmetricTensor<rank_, dim, Number>::n_independent_components;
 
namespace internal
{
  namespace SymmetricTensorAccessors
  {
    template <int rank_, int dim, bool constness, int P, typename Number>
    constexpr DEAL_II_ALWAYS_INLINE
    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>
    constexpr inline DEAL_II_ALWAYS_INLINE
        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>
    constexpr DEAL_II_ALWAYS_INLINE
        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>
    constexpr DEAL_II_ALWAYS_INLINE
    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>
    constexpr inline DEAL_II_ALWAYS_INLINE
      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>
    constexpr DEAL_II_ALWAYS_INLINE
      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));
    }
  } // namespace SymmetricTensorAccessors
} // namespace internal
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
  const Tensor<2, dim, OtherNumber> &t)
{
  static_assert(rank == 2, "This function is only implemented for rank==2");
  for (unsigned int d = 0; d < dim; ++d)
    for (unsigned int e = 0; e < d; ++e)
      Assert(t[d][e] == t[e][d],
             ExcMessage("The incoming Tensor must be exactly symmetric."));
 
  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>
constexpr DEAL_II_ALWAYS_INLINE
SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
  : data(initializer.data)
{}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_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>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
                                       SymmetricTensor<rank_, dim, Number>::
                                       operator=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
{
  data = t.data;
  return *this;
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
SymmetricTensor<rank_, dim, Number>::operator=(const Number &d)
{
  Assert(numbers::value_is_zero(d),
         ExcMessage("Only assignment with zero is allowed"));
  (void)d;
 
  data = internal::NumberType<Number>::value(0.0);
 
  return *this;
}
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    template <int dim, typename Number>
    constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
                                           convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
    {
      ::Tensor<2, dim, Number> t;
 
      // 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 t;
    }
 
 
    template <int dim, typename Number>
    constexpr ::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 <typename Number>
    struct Inverse<2, 1, Number>
    {
      constexpr static inline DEAL_II_ALWAYS_INLINE
        ::SymmetricTensor<2, 1, Number>
        value(const ::SymmetricTensor<2, 1, Number> &t)
      {
        ::SymmetricTensor<2, 1, Number> tmp;
 
        tmp[0][0] = 1.0 / t[0][0];
 
        return tmp;
      }
    };
 
 
    template <typename Number>
    struct Inverse<2, 2, Number>
    {
      constexpr static inline DEAL_II_ALWAYS_INLINE
        ::SymmetricTensor<2, 2, Number>
        value(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>
    struct Inverse<2, 3, Number>
    {
      constexpr static ::SymmetricTensor<2, 3, Number>
      value(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)]  ])
        //
        // =
        //
        // [  (t11*t22 - t12**2)/det_t,
        //    (-t01*t22 + t02*t12)/det_t,
        //    (t01*t12 - t02*t11)/det_t],
        // [  (-t01*t22 + t02*t12)/det_t,
        //    (t00*t22 - t02**2)/det_t,
        //    (-t00*t12 + t01*t02)/det_t],
        // [  (t01*t12 - t02*t11)/det_t,
        //    (-t00*t12 + t01*t02)/det_t,
        //    (t00*t11 - t01**2)/det_t]   ])
        //
        // with det_t = (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;
      }
    };
 
 
    template <typename Number>
    struct Inverse<4, 1, Number>
    {
      constexpr static inline ::SymmetricTensor<4, 1, Number>
      value(const ::SymmetricTensor<4, 1, Number> &t)
      {
        ::SymmetricTensor<4, 1, Number> tmp;
        tmp.data[0][0] = 1.0 / t.data[0][0];
        return tmp;
      }
    };
 
 
    template <typename Number>
    struct Inverse<4, 2, Number>
    {
      constexpr static inline ::SymmetricTensor<4, 2, Number>
      value(const ::SymmetricTensor<4, 2, Number> &t)
      {
        ::SymmetricTensor<4, 2, Number> tmp;
 
        // 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;
 
        return tmp;
      }
    };
 
 
    template <typename Number>
    struct Inverse<4, 3, Number>
    {
      static ::SymmetricTensor<4, 3, Number>
      value(const ::SymmetricTensor<4, 3, Number> &t)
      {
        ::SymmetricTensor<4, 3, Number> tmp = t;
 
        // This function follows the exact same scheme as the 2d case, except
        // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
        // Instead, we use the Gauss-Jordan algorithm implemented for
        // FullMatrix. For historical reasons the following code is copied from
        // there, with the tangential benefit that we do not need to copy the
        // tensor entries to and from the FullMatrix.
        const unsigned int N = 6;
 
        // First get an estimate of the size of the elements of this matrix,
        // for later checks whether the pivot element is large enough, or
        // whether we have to fear that the matrix is not regular.
        Number diagonal_sum = internal::NumberType<Number>::value(0.0);
        for (unsigned int i = 0; i < N; ++i)
          diagonal_sum += std::fabs(tmp.data[i][i]);
        const Number typical_diagonal_element =
          diagonal_sum / static_cast<double>(N);
        (void)typical_diagonal_element;
 
        unsigned int p[N];
        for (unsigned int i = 0; i < N; ++i)
          p[i] = i;
 
        for (unsigned int j = 0; j < N; ++j)
          {
            // Pivot search: search that part of the line on and right of the
            // diagonal for the largest element.
            Number       max = std::fabs(tmp.data[j][j]);
            unsigned int r   = j;
            for (unsigned int i = j + 1; i < N; ++i)
              if (std::fabs(tmp.data[i][j]) > max)
                {
                  max = std::fabs(tmp.data[i][j]);
                  r   = i;
                }
 
            // Check whether the pivot is too small
            Assert(max > 1.e-16 * typical_diagonal_element,
                   ExcMessage("This tensor seems to be noninvertible"));
 
            // Row interchange
            if (r > j)
              {
                for (unsigned int k = 0; k < N; ++k)
                  std::swap(tmp.data[j][k], tmp.data[r][k]);
 
                std::swap(p[j], p[r]);
              }
 
            // Transformation
            const Number hr = 1. / tmp.data[j][j];
            tmp.data[j][j]  = hr;
            for (unsigned int k = 0; k < N; ++k)
              {
                if (k == j)
                  continue;
                for (unsigned int i = 0; i < N; ++i)
                  {
                    if (i == j)
                      continue;
                    tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
                  }
              }
            for (unsigned int i = 0; i < N; ++i)
              {
                tmp.data[i][j] *= hr;
                tmp.data[j][i] *= -hr;
              }
            tmp.data[j][j] = hr;
          }
 
        // Column interchange
        Number hv[N];
        for (unsigned int i = 0; i < N; ++i)
          {
            for (unsigned int k = 0; k < N; ++k)
              hv[p[k]] = tmp.data[i][k];
            for (unsigned int k = 0; k < N; ++k)
              tmp.data[i][k] = hv[k];
          }
 
        // Scale rows and columns. The mult matrix
        // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
        for (unsigned int i = 3; i < 6; ++i)
          for (unsigned int j = 0; j < 3; ++j)
            tmp.data[i][j] /= 2;
 
        for (unsigned int i = 0; i < 3; ++i)
          for (unsigned int j = 3; j < 6; ++j)
            tmp.data[i][j] /= 2;
 
        for (unsigned int i = 3; i < 6; ++i)
          for (unsigned int j = 3; j < 6; ++j)
            tmp.data[i][j] /= 4;
 
        return tmp;
      }
    };
 
  } // namespace SymmetricTensorImplementation
} // namespace internal
 
 
 
template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
                                operator Tensor<rank_, dim, Number>() const
{
  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr bool
SymmetricTensor<rank_, dim, Number>::
operator==(const SymmetricTensor<rank_, dim, Number> &t) const
{
  return data == t.data;
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr bool
SymmetricTensor<rank_, dim, Number>::
operator!=(const SymmetricTensor<rank_, dim, Number> &t) const
{
  return data != t.data;
}
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
                                       SymmetricTensor<rank_, dim, Number>::
                                       operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
{
  data += t.data;
  return *this;
}
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
                                       SymmetricTensor<rank_, dim, Number>::
                                       operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
{
  data -= t.data;
  return *this;
}
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
SymmetricTensor<rank_, dim, Number>::operator*=(const OtherNumber &d)
{
  data *= d;
  return *this;
}
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
SymmetricTensor<rank_, dim, Number>::operator/=(const OtherNumber &d)
{
  data /= d;
  return *this;
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_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>
constexpr inline DEAL_II_ALWAYS_INLINE void
SymmetricTensor<rank_, dim, Number>::clear()
{
  data.clear();
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr 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, typename OtherNumber = Number>
  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
    typename SymmetricTensorAccessors::
      double_contraction_result<2, 2, dim, Number, OtherNumber>::type
      perform_double_contraction(
        const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
          base_tensor_type &data,
        const typename SymmetricTensorAccessors::
          StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
  {
    using result_type = typename SymmetricTensorAccessors::
      double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
 
    switch (dim)
      {
        case 1:
          return data[0] * sdata[0];
        default:
          // Start with the non-diagonal part to avoid some multiplications by
          // 2.
 
          result_type 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, typename OtherNumber = Number>
  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
    typename SymmetricTensorAccessors::
      double_contraction_result<4, 2, dim, Number, OtherNumber>::type
      perform_double_contraction(
        const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
          base_tensor_type &data,
        const typename SymmetricTensorAccessors::
          StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
  {
    using result_type = typename SymmetricTensorAccessors::
      double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
    using value_type = typename SymmetricTensorAccessors::
      double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
 
    const unsigned int data_dim = SymmetricTensorAccessors::
      StorageType<2, dim, value_type>::n_independent_components;
    value_type tmp[data_dim]{};
    for (unsigned int i = 0; i < data_dim; ++i)
      tmp[i] =
        perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
    return result_type(tmp);
  }
 
 
 
  template <int dim, typename Number, typename OtherNumber = Number>
  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
    typename SymmetricTensorAccessors::StorageType<
      2,
      dim,
      typename SymmetricTensorAccessors::
        double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
      base_tensor_type
      perform_double_contraction(
        const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
          base_tensor_type &data,
        const typename SymmetricTensorAccessors::
          StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
  {
    using value_type = typename SymmetricTensorAccessors::
      double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
    using base_tensor_type = typename SymmetricTensorAccessors::
      StorageType<2, dim, value_type>::base_tensor_type;
 
    base_tensor_type tmp;
    for (unsigned int i = 0; i < tmp.dimension; ++i)
      {
        // Start with the non-diagonal part
        value_type 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, typename OtherNumber = Number>
  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
    typename SymmetricTensorAccessors::StorageType<
      4,
      dim,
      typename SymmetricTensorAccessors::
        double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
      base_tensor_type
      perform_double_contraction(
        const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
          base_tensor_type &data,
        const typename SymmetricTensorAccessors::
          StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
  {
    using value_type = typename SymmetricTensorAccessors::
      double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
    using base_tensor_type = typename SymmetricTensorAccessors::
      StorageType<4, dim, value_type>::base_tensor_type;
 
    const unsigned int data_dim = SymmetricTensorAccessors::
      StorageType<2, dim, value_type>::n_independent_components;
    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>
template <typename OtherNumber>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
  typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
      SymmetricTensor<rank_, dim, Number>::
      operator*(const SymmetricTensor<2, dim, OtherNumber> &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, OtherNumber>(data,
                                                                        s.data);
}
 
 
 
template <int rank_, int dim, typename Number>
template <typename OtherNumber>
DEAL_II_CONSTEXPR inline typename internal::SymmetricTensorAccessors::
  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
    SymmetricTensor<rank_, dim, Number>::
    operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const
{
  typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
  tmp.data =
    internal::perform_double_contraction<dim, Number, OtherNumber>(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
{
  // The variables within this struct will be referenced in the next functions.
  // It is a workaround that allows returning a reference to a static variable
  // while allowing constexpr evaluation of the function.
  // It has to be defined outside the function because constexpr functions
  // cannot define static variables.
  // A similar struct has also been defined in tensor.h
  template <typename Type>
  struct Uninitialized
  {
    static Type value;
  };
 
  template <typename Type>
  Type Uninitialized<Type>::value;
 
  template <int dim, typename Number>
  constexpr inline DEAL_II_ALWAYS_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(std::min(indices[0], indices[1]),
                                           std::max(indices[0], indices[1]));
            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());
          }
      }
 
    // The code should never reach there.
    // Returns a dummy reference to a dummy variable just to make the
    // compiler happy.
    return Uninitialized<Number>::value;
  }
 
 
 
  template <int dim, typename Number>
  constexpr inline DEAL_II_ALWAYS_INLINE const 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(std::min(indices[0], indices[1]),
                                           std::max(indices[0], indices[1]));
            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());
          }
      }
 
    // The code should never reach there.
    // Returns a dummy reference to a dummy variable just to make the
    // compiler happy.
    return Uninitialized<Number>::value;
  }
 
 
 
  template <int dim, typename Number>
  constexpr 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
          {
            constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
            return data[base_index[indices[0]][indices[1]]]
                       [base_index[indices[2]][indices[3]]];
          }
        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
          {
            constexpr std::size_t base_index[3][3] = {{0, 3, 4},
                                                      {3, 1, 5},
                                                      {4, 5, 2}};
            return data[base_index[indices[0]][indices[1]]]
                       [base_index[indices[2]][indices[3]]];
          }
 
        default:
          Assert(false, ExcNotImplemented());
      }
 
    // The code should never reach there.
    // Returns a dummy reference to a dummy variable just to make the
    // compiler happy.
    return Uninitialized<Number>::value;
  }
 
 
  template <int dim, typename Number>
  constexpr inline DEAL_II_ALWAYS_INLINE const 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
          {
            constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
            return data[base_index[indices[0]][indices[1]]]
                       [base_index[indices[2]][indices[3]]];
          }
        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
          {
            constexpr std::size_t base_index[3][3] = {{0, 3, 4},
                                                      {3, 1, 5},
                                                      {4, 5, 2}};
            return data[base_index[indices[0]][indices[1]]]
                       [base_index[indices[2]][indices[3]]];
          }
 
        default:
          Assert(false, ExcNotImplemented());
      }
 
    // The code should never reach there.
    // Returns a dummy reference to a dummy variable just to make the
    // compiler happy.
    return Uninitialized<Number>::value;
  }
 
} // end of namespace internal
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE Number &
                                       SymmetricTensor<rank_, dim, Number>::
                                       operator()(const TableIndices<rank_> &indices)
{
  for (unsigned int r = 0; r < rank; ++r)
    AssertIndexRange(indices[r], dimension);
  return internal::symmetric_tensor_access<dim, Number>(indices, data);
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE const Number &
                                             SymmetricTensor<rank_, dim, Number>::
                                             operator()(const TableIndices<rank_> &indices) const
{
  for (unsigned int r = 0; r < rank; ++r)
    AssertIndexRange(indices[r], dimension);
  return internal::symmetric_tensor_access<dim, Number>(indices, data);
}
 
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    template <int rank_>
    constexpr TableIndices<rank_>
    get_partially_filled_indices(const unsigned int row,
                                 const std::integral_constant<int, 2> &)
    {
      return TableIndices<rank_>(row, numbers::invalid_unsigned_int);
    }
 
 
    template <int rank_>
    constexpr TableIndices<rank_>
    get_partially_filled_indices(const unsigned int row,
                                 const std::integral_constant<int, 4> &)
    {
      return TableIndices<rank_>(row,
                                 numbers::invalid_unsigned_int,
                                 numbers::invalid_unsigned_int,
                                 numbers::invalid_unsigned_int);
    }
  } // namespace SymmetricTensorImplementation
} // namespace internal
 
 
template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE 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,
      internal::SymmetricTensorImplementation::get_partially_filled_indices<
        rank_>(row, std::integral_constant<int, rank_>()));
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE 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,
      internal::SymmetricTensorImplementation::get_partially_filled_indices<
        rank_>(row, std::integral_constant<int, rank_>()));
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE const Number &
                                      SymmetricTensor<rank_, dim, Number>::
                                      operator[](const TableIndices<rank_> &indices) const
{
  return operator()(indices);
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE Number &
                                       SymmetricTensor<rank_, dim, Number>::
                                       operator[](const TableIndices<rank_> &indices)
{
  return operator()(indices);
}
 
 
 
template <int rank_, int dim, typename Number>
inline Number *
SymmetricTensor<rank_, dim, Number>::begin_raw()
{
  return std::addressof(this->access_raw_entry(0));
}
 
 
 
template <int rank_, int dim, typename Number>
inline const Number *
SymmetricTensor<rank_, dim, Number>::begin_raw() const
{
  return std::addressof(this->access_raw_entry(0));
}
 
 
 
template <int rank_, int dim, typename Number>
inline Number *
SymmetricTensor<rank_, dim, Number>::end_raw()
{
  return begin_raw() + n_independent_components;
}
 
 
 
template <int rank_, int dim, typename Number>
inline const Number *
SymmetricTensor<rank_, dim, Number>::end_raw() const
{
  return begin_raw() + n_independent_components;
}
 
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    template <int dim, typename Number>
    constexpr unsigned int
    entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
                     const unsigned int index)
    {
      return index;
    }
 
 
    template <int dim, typename Number>
    constexpr ::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);
    }
 
  } // namespace SymmetricTensorImplementation
} // namespace internal
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline const Number &
SymmetricTensor<rank_, dim, Number>::access_raw_entry(
  const unsigned int index) const
{
  AssertIndexRange(index, n_independent_components);
  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
                                                                        index)];
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline Number &
SymmetricTensor<rank_, dim, Number>::access_raw_entry(const unsigned int index)
{
  AssertIndexRange(index, n_independent_components);
  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
                                                                        index)];
}
 
 
 
namespace internal
{
  template <int dim, typename Number>
  constexpr 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>
  constexpr 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>
constexpr typename numbers::NumberTraits<Number>::real_type
SymmetricTensor<rank_, dim, Number>::norm() const
{
  return internal::compute_norm<dim, Number>(data);
}
 
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    // 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>
    constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
    component_to_unrolled_index(const TableIndices<2> &indices)
    {
      AssertIndexRange(indices[0], dim);
      AssertIndexRange(indices[1], dim);
 
      switch (dim)
        {
          case 1:
            {
              return 0;
            }
 
          case 2:
            {
              constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
              return table[indices[0]][indices[1]];
            }
 
          case 3:
            {
              constexpr unsigned int table[3][3] = {{0, 3, 4},
                                                    {3, 1, 5},
                                                    {4, 5, 2}};
              return table[indices[0]][indices[1]];
            }
 
          case 4:
            {
              constexpr 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_>
    constexpr inline unsigned int
    component_to_unrolled_index(const TableIndices<rank_> &indices)
    {
      (void)indices;
      Assert(false, ExcNotImplemented());
      return numbers::invalid_unsigned_int;
    }
  } // namespace SymmetricTensorImplementation
} // namespace internal
 
 
template <int rank_, int dim, typename Number>
constexpr unsigned int
SymmetricTensor<rank_, dim, Number>::component_to_unrolled_index(
  const TableIndices<rank_> &indices)
{
  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
    dim>(indices);
}
 
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    // 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>
    constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
                                           unrolled_to_component_indices(const unsigned int i,
                                                                         const std::integral_constant<int, 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 {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 {i, i};
 
            for (unsigned int d = 0, c = dim; d < dim; ++d)
              for (unsigned int e = d + 1; e < dim; ++e, ++c)
                if (c == i)
                  return {d, e};
 
            // should never get here:
            Assert(false, ExcInternalError());
            return {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_>
    constexpr inline
      typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
      unrolled_to_component_indices(const unsigned int i,
                                    const std::integral_constant<int, 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_>();
    }
 
  } // namespace SymmetricTensorImplementation
} // namespace internal
 
template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE TableIndices<rank_>
                                SymmetricTensor<rank_, dim, Number>::unrolled_to_component_indices(
  const unsigned int i)
{
  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
    dim>(i, std::integral_constant<int, 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>
constexpr inline DEAL_II_ALWAYS_INLINE
  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
  operator+(const SymmetricTensor<rank_, dim, Number> &     left,
            const SymmetricTensor<rank_, dim, OtherNumber> &right)
{
  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
    tmp = left;
  tmp += right;
  return tmp;
}
 
 
template <int rank_, int dim, typename Number, typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE
  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
  operator-(const SymmetricTensor<rank_, dim, Number> &     left,
            const SymmetricTensor<rank_, dim, OtherNumber> &right)
{
  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
    tmp = left;
  tmp -= right;
  return tmp;
}
 
 
template <int rank_, int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_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>
constexpr DEAL_II_ALWAYS_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>
constexpr DEAL_II_ALWAYS_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>
constexpr DEAL_II_ALWAYS_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>
constexpr inline DEAL_II_ALWAYS_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 internal::NumberType<Number>::value(0.0);
    }
}
 
 
 
template <int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE Number
                                third_invariant(const SymmetricTensor<2, dim, Number> &t)
{
  return determinant(t);
}
 
 
 
template <int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_INLINE 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>
constexpr Number
first_invariant(const SymmetricTensor<2, dim, Number> &t)
{
  return trace(t);
}
 
 
template <typename Number>
constexpr DEAL_II_ALWAYS_INLINE Number
                                second_invariant(const SymmetricTensor<2, 1, Number> &)
{
  return internal::NumberType<Number>::value(0.0);
}
 
 
 
template <typename Number>
constexpr DEAL_II_ALWAYS_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>
constexpr DEAL_II_ALWAYS_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 <typename Number>
std::array<Number, 1>
eigenvalues(const SymmetricTensor<2, 1, Number> &T);
 
 
 
template <typename Number>
std::array<Number, 2>
eigenvalues(const SymmetricTensor<2, 2, Number> &T);
 
 
 
template <typename Number>
std::array<Number, 3>
eigenvalues(const SymmetricTensor<2, 3, Number> &T);
 
 
 
namespace internal
{
  namespace SymmetricTensorImplementation
  {
    template <int dim, typename Number>
    void
    tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
                   ::Tensor<2, dim, Number> &               Q,
                   std::array<Number, dim> &                      d,
                   std::array<Number, dim - 1> &                  e);
 
 
 
    template <int dim, typename Number>
    std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
    ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
 
 
 
    template <int dim, typename Number>
    std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
      jacobi(::SymmetricTensor<2, dim, Number> A);
 
 
 
    template <typename Number>
    std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
    hybrid(const ::SymmetricTensor<2, 2, Number> &A);
 
 
 
    template <typename Number>
    std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
    hybrid(const ::SymmetricTensor<2, 3, Number> &A);
 
    template <int dim, typename Number>
    struct SortEigenValuesVectors
    {
      using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
      bool
      operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
      {
        return lhs.first > rhs.first;
      }
    };
 
  } // namespace SymmetricTensorImplementation
 
} // namespace internal
 
 
 
// The line below is to ensure that doxygen puts the full description
// of this global enumeration into the documentation
// See https://stackoverflow.com/a/1717984
enum struct SymmetricTensorEigenvectorMethod
{
  hybrid,
  ql_implicit_shifts,
  jacobi
};
 
 
 
template <int dim, typename Number>
std::array<std::pair<Number, Tensor<1, dim, Number>>,
           std::integral_constant<int, dim>::value>
eigenvectors(const SymmetricTensor<2, dim, Number> &T,
             const SymmetricTensorEigenvectorMethod method =
               SymmetricTensorEigenvectorMethod::ql_implicit_shifts);
 
 
 
template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
                                transpose(const SymmetricTensor<rank_, dim, Number> &t)
{
  return t;
}
 
 
 
template <int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_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>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_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] = internal::NumberType<Number>::value(1.);
        break;
      case 2:
        tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
        break;
      case 3:
        tmp.data[0] = tmp.data[1] = tmp.data[2] =
          internal::NumberType<Number>::value(1.);
        break;
      default:
        for (unsigned int d = 0; d < dim; ++d)
          tmp.data[d] = internal::NumberType<Number>::value(1.);
    }
  return tmp;
}
 
 
 
template <int dim>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim>
                                               unit_symmetric_tensor()
{
  return unit_symmetric_tensor<dim, double>();
}
 
 
 
template <int dim, typename Number>
DEAL_II_CONSTEXPR 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] =
        internal::NumberType<Number>::value((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] = internal::NumberType<Number>::value(0.5);
 
  return tmp;
}
 
 
 
template <int dim>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
                                               deviator_tensor()
{
  return deviator_tensor<dim, double>();
}
 
 
 
template <int dim, typename Number>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_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] = internal::NumberType<Number>::value(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] = internal::NumberType<Number>::value(0.5);
 
  return tmp;
}
 
 
 
template <int dim>
DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
                                               identity_tensor()
{
  return identity_tensor<dim, double>();
}
 
 
 
template <int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
                                invert(const SymmetricTensor<2, dim, Number> &t)
{
  return internal::SymmetricTensorImplementation::Inverse<2, dim, Number>::
    value(t);
}
 
 
 
template <int dim, typename Number>
constexpr SymmetricTensor<4, dim, Number>
invert(const SymmetricTensor<4, dim, Number> &t)
{
  return internal::SymmetricTensorImplementation::Inverse<4, dim, Number>::
    value(t);
}
 
 
 
template <int dim, typename Number>
constexpr 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>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
                                       symmetrize(const Tensor<2, dim, Number> &t)
{
  SymmetricTensor<2, dim, Number> result;
  for (unsigned int d = 0; d < dim; ++d)
    result[d][d] = t[d][d];
 
  const Number half = internal::NumberType<Number>::value(0.5);
  for (unsigned int d = 0; d < dim; ++d)
    for (unsigned int e = d + 1; e < dim; ++e)
      result[d][e] = (t[d][e] + t[e][d]) * half;
  return result;
}
 
 
 
template <int rank_, int dim, typename Number>
constexpr inline DEAL_II_ALWAYS_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>
constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
                                operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
{
  // simply forward to the other operator
  return t * factor;
}
 
 
 
template <int rank_, int dim, typename Number, typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_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).
  using product_type = typename ProductType<Number, OtherNumber>::type;
  SymmetricTensor<rank_, dim, product_type> tt(t);
  tt *= internal::NumberType<product_type>::value(factor);
  return tt;
}
 
 
 
template <int rank_, int dim, typename Number, typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
  rank_,
  dim,
  typename ProductType<OtherNumber,
                       typename EnableIfScalar<Number>::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, typename OtherNumber>
constexpr inline SymmetricTensor<
  rank_,
  dim,
  typename ProductType<Number,
                       typename EnableIfScalar<OtherNumber>::type>::type>
operator/(const SymmetricTensor<rank_, dim, Number> &t,
          const OtherNumber &                        factor)
{
  using product_type = typename ProductType<Number, OtherNumber>::type;
  SymmetricTensor<rank_, dim, product_type> tt(t);
  tt /= internal::NumberType<product_type>::value(factor);
  return tt;
}
 
 
 
template <int rank_, int dim>
constexpr inline DEAL_II_ALWAYS_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>
constexpr inline DEAL_II_ALWAYS_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>
constexpr 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, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
scalar_product(const SymmetricTensor<2, dim, Number> &     t1,
               const SymmetricTensor<2, dim, OtherNumber> &t2)
{
  return (t1 * t2);
}
 
 
template <int dim, typename Number, typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE
  typename ProductType<Number, OtherNumber>::type
  scalar_product(const SymmetricTensor<2, dim, Number> &t1,
                 const Tensor<2, dim, OtherNumber> &    t2)
{
  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
    typename ProductType<Number, OtherNumber>::type>::value(0.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, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
scalar_product(const Tensor<2, dim, Number> &              t1,
               const SymmetricTensor<2, dim, OtherNumber> &t2)
{
  return scalar_product(t2, t1);
}
 
 
template <typename Number, typename OtherNumber>
constexpr inline DEAL_II_ALWAYS_INLINE void double_contract(
  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<4, 1, Number> &                                   t,
  const SymmetricTensor<2, 1, OtherNumber> &                              s)
{
  tmp[0][0] = t[0][0][0][0] * s[0][0];
}
 
 
 
template <typename Number, typename OtherNumber>
constexpr inline void double_contract(
  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<2, 1, Number> &                                   s,
  const SymmetricTensor<4, 1, OtherNumber> &                              t)
{
  tmp[0][0] = t[0][0][0][0] * s[0][0];
}
 
 
 
template <typename Number, typename OtherNumber>
constexpr inline void double_contract(
  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<4, 2, Number> &                                   t,
  const SymmetricTensor<2, 2, OtherNumber> &                              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, typename OtherNumber>
constexpr inline void double_contract(
  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<2, 2, Number> &                                   s,
  const SymmetricTensor<4, 2, OtherNumber> &                              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, typename OtherNumber>
constexpr inline void double_contract(
  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<4, 3, Number> &                                   t,
  const SymmetricTensor<2, 3, OtherNumber> &                              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, typename OtherNumber>
constexpr inline void double_contract(
  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
  const SymmetricTensor<2, 3, Number> &                                   s,
  const SymmetricTensor<4, 3, OtherNumber> &                              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, typename OtherNumber>
constexpr Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
operator*(const SymmetricTensor<2, dim, Number> &src1,
          const Tensor<1, dim, OtherNumber> &    src2)
{
  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type> 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, typename OtherNumber>
constexpr Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
operator*(const Tensor<1, dim, Number> &              src1,
          const SymmetricTensor<2, dim, OtherNumber> &src2)
{
  // this is easy for symmetric tensors:
  return src2 * src1;
}
 
 
 
template <int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  typename Tensor<rank_1 + rank_2 - 2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
  operator*(const Tensor<rank_1, dim, Number> &              src1,
            const SymmetricTensor<rank_2, dim, OtherNumber> &src2)
{
  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
}
 
 
 
template <int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  typename Tensor<rank_1 + rank_2 - 2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
  operator*(const SymmetricTensor<rank_1, dim, Number> &src1,
            const Tensor<rank_2, dim, OtherNumber> &    src2)
{
  return Tensor<rank_1, dim, Number>(src1) * src2;
}
 
 
 
template <int dim, typename Number>
inline std::ostream &
operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
{
  // make our lives a bit simpler by outputting
  // 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 our lives a bit simpler by outputting
  // 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