tensor.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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_tensor_h
#define dealii_tensor_h

#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.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_accessors.h>
#include <deal.II/base/utilities.h>

#ifdef DEAL_II_WITH_ADOLC
#  include <adolc/adouble.h> // Taped double
#endif

#include <cmath>
#include <ostream>
#include <vector>


DEAL_II_NAMESPACE_OPEN

// Forward declarations:

template <int dim, typename Number>
class Point;
template <int rank_, int dim, typename Number = double>
class Tensor;
template <typename Number>
class Vector;
template <typename Number>
class VectorizedArray;

#ifndef DOXYGEN
// Overload invalid tensor types of negative rank that come up during
// overload resolution of operator* and related contraction variants.
template <int dim, typename Number>
class Tensor<-2, dim, Number>
{};

template <int dim, typename Number>
class Tensor<-1, dim, Number>
{};
#endif /* DOXYGEN */


template <int dim, typename Number>
class Tensor<0, dim, Number>
{
public:
  static constexpr unsigned int dimension = dim;

  static constexpr unsigned int rank = 0;

  static constexpr unsigned int n_independent_components = 1;

  using real_type = typename numbers::NumberTraits<Number>::real_type;

  using value_type = Number;

  using array_type = Number;

  DEAL_II_CUDA_HOST_DEV
  Tensor();

  template <typename OtherNumber>
  Tensor(const Tensor<0, dim, OtherNumber> &initializer);

  template <typename OtherNumber>
  Tensor(const OtherNumber &initializer);

  Number *
  begin_raw();

  constexpr const Number *
  begin_raw() const;

  Number *
  end_raw();

  constexpr const Number *
  end_raw() const;

  DEAL_II_CUDA_HOST_DEV operator Number &();

  DEAL_II_CUDA_HOST_DEV operator const Number &() const;

  template <typename OtherNumber>
  Tensor &
  operator=(const Tensor<0, dim, OtherNumber> &rhs);

#ifdef __INTEL_COMPILER

  Tensor &
  operator=(const Tensor<0, dim, Number> &rhs);
#endif

  template <typename OtherNumber>
  Tensor &
  operator=(const OtherNumber &d);

  template <typename OtherNumber>
  bool
  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;

  template <typename OtherNumber>
  constexpr bool
  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;

  template <typename OtherNumber>
  Tensor &
  operator+=(const Tensor<0, dim, OtherNumber> &rhs);

  template <typename OtherNumber>
  Tensor &
  operator-=(const Tensor<0, dim, OtherNumber> &rhs);

  template <typename OtherNumber>
  DEAL_II_CUDA_HOST_DEV Tensor &
                        operator*=(const OtherNumber &factor);

  template <typename OtherNumber>
  Tensor &
  operator/=(const OtherNumber &factor);

  constexpr Tensor
  operator-() const;

  void
  clear();

  real_type
  norm() const;

  DEAL_II_CUDA_HOST_DEV real_type
                        norm_square() const;

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

  using tensor_type = Number;

private:
  Number value;

  template <typename OtherNumber>
  void
  unroll_recursion(Vector<OtherNumber> &result,
                   unsigned int &       start_index) const;

  template <int, int, typename>
  friend class Tensor;
};



template <int rank_, int dim, typename Number>
class Tensor
{
public:
  static constexpr unsigned int dimension = dim;

  static constexpr unsigned int rank = rank_;

  static constexpr unsigned int n_independent_components =
    Tensor<rank_ - 1, dim>::n_independent_components * dim;

  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;

  using array_type =
    typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];

  constexpr DEAL_II_CUDA_HOST_DEV
  Tensor() = default;

  explicit Tensor(const array_type &initializer);

  template <typename OtherNumber>
  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);

  template <typename OtherNumber>
  Tensor(
    const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);

  template <typename OtherNumber>
  constexpr
  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;

  DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);

  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
                                        operator[](const unsigned int i) const;

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

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

  Number *
  begin_raw();

  constexpr const Number *
  begin_raw() const;

  Number *
  end_raw();

  constexpr const Number *
  end_raw() const;

  template <typename OtherNumber>
  Tensor &
  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);

  Tensor &
  operator=(const Number &d);

  template <typename OtherNumber>
  bool
  operator==(const Tensor<rank_, dim, OtherNumber> &) const;

  template <typename OtherNumber>
  constexpr bool
  operator!=(const Tensor<rank_, dim, OtherNumber> &) const;

  template <typename OtherNumber>
  Tensor &
  operator+=(const Tensor<rank_, dim, OtherNumber> &);

  template <typename OtherNumber>
  Tensor &
  operator-=(const Tensor<rank_, dim, OtherNumber> &);

  template <typename OtherNumber>
  DEAL_II_CUDA_HOST_DEV Tensor &
                        operator*=(const OtherNumber &factor);

  template <typename OtherNumber>
  Tensor &
  operator/=(const OtherNumber &factor);

  Tensor
  operator-() const;

  void
  clear();

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

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

  template <typename OtherNumber>
  void
  unroll(Vector<OtherNumber> &result) const;

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

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

  static constexpr std::size_t
  memory_consumption();

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

  using tensor_type = Tensor<rank_, dim, Number>;

private:
  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
  // ... avoid a compiler warning in case of dim == 0 and ensure that the
  // array always has positive size.

  template <typename OtherNumber>
  void
  unroll_recursion(Vector<OtherNumber> &result,
                   unsigned int &       start_index) const;

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

  friend class Point<dim, Number>;
};


namespace internal
{
  template <int rank, int dim, typename T>
  struct NumberType<Tensor<rank, dim, T>>
  {
    static constexpr const Tensor<rank, dim, T> &
    value(const Tensor<rank, dim, T> &t)
    {
      return t;
    }

    static Tensor<rank, dim, T>
    value(const T &t)
    {
      Tensor<rank, dim, T> tmp;
      tmp = t;
      return tmp;
    }
  };

  template <int rank, int dim, typename T>
  struct NumberType<Tensor<rank, dim, VectorizedArray<T>>>
  {
    static constexpr const Tensor<rank, dim, VectorizedArray<T>> &
    value(const Tensor<rank, dim, VectorizedArray<T>> &t)
    {
      return t;
    }

    static Tensor<rank, dim, VectorizedArray<T>>
    value(const T &t)
    {
      Tensor<rank, dim, VectorizedArray<T>> tmp;
      tmp = internal::NumberType<VectorizedArray<T>>::value(t);
      return tmp;
    }

    static Tensor<rank, dim, VectorizedArray<T>>
    value(const VectorizedArray<T> &t)
    {
      Tensor<rank, dim, VectorizedArray<T>> tmp;
      tmp = t;
      return tmp;
    }
  };
} // namespace internal


/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/


template <int dim, typename Number>
DEAL_II_CUDA_HOST_DEV inline Tensor<0, dim, Number>::Tensor()
  // Some auto-differentiable numbers need explicit
  // zero initialization.
  : value(internal::NumberType<Number>::value(0.0))
{}



template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
{
  value = internal::NumberType<Number>::value(initializer);
}



template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p)
{
  value = p.value;
}



template <int dim, typename Number>
inline Number *
Tensor<0, dim, Number>::begin_raw()
{
  return std::addressof(value);
}



template <int dim, typename Number>
constexpr const Number *
Tensor<0, dim, Number>::begin_raw() const
{
  return std::addressof(value);
}



template <int dim, typename Number>
inline Number *
Tensor<0, dim, Number>::end_raw()
{
  return begin_raw() + n_independent_components;
}



template <int dim, typename Number>
constexpr const Number *
Tensor<0, dim, Number>::end_raw() const
{
  return begin_raw() + n_independent_components;
}



template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
{
  // We cannot use Assert inside a CUDA kernel
#ifndef __CUDA_ARCH__
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
#endif
  return value;
}


template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
                             operator const Number &() const
{
  // We cannot use Assert inside a CUDA kernel
#ifndef __CUDA_ARCH__
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
#endif
  return value;
}


template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
{
  value = internal::NumberType<Number>::value(p);
  return *this;
}


#ifdef __INTEL_COMPILER
template <int dim, typename Number>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p)
{
  value = p.value;
  return *this;
}
#endif


template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator=(const OtherNumber &d)
{
  value = internal::NumberType<Number>::value(d);
  return *this;
}


template <int dim, typename Number>
template <typename OtherNumber>
inline bool
Tensor<0, dim, Number>::operator==(const Tensor<0, dim, OtherNumber> &p) const
{
#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
  Assert(!(std::is_same<Number, adouble>::value ||
           std::is_same<OtherNumber, adouble>::value),
         ExcMessage(
           "The Tensor equality operator for ADOL-C taped numbers has not yet "
           "been extended to support advanced branching."));
#endif

  return numbers::values_are_equal(value, p.value);
}


template <int dim, typename Number>
template <typename OtherNumber>
constexpr bool
Tensor<0, dim, Number>::operator!=(const Tensor<0, dim, OtherNumber> &p) const
{
  return !((*this) == p);
}


template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p)
{
  value += p.value;
  return *this;
}


template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator-=(const Tensor<0, dim, OtherNumber> &p)
{
  value -= p.value;
  return *this;
}



namespace internal
{
  namespace ComplexWorkaround
  {
    template <typename Number, typename OtherNumber>
    inline DEAL_II_CUDA_HOST_DEV void
    multiply_assign_scalar(Number &val, const OtherNumber &s)
    {
      val *= s;
    }

#ifdef __CUDA_ARCH__
    template <typename Number, typename OtherNumber>
    inline DEAL_II_CUDA_HOST_DEV void
    multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
    {
      printf("This function is not implemented for std::complex<Number>!\n");
      assert(false);
    }
#endif
  } // namespace ComplexWorkaround
} // namespace internal


template <int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
{
  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
  return *this;
}



template <int dim, typename Number>
template <typename OtherNumber>
inline Tensor<0, dim, Number> &
Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
{
  value /= s;
  return *this;
}


template <int dim, typename Number>
constexpr Tensor<0, dim, Number>
Tensor<0, dim, Number>::operator-() const
{
  return -value;
}


template <int dim, typename Number>
inline typename Tensor<0, dim, Number>::real_type
Tensor<0, dim, Number>::norm() const
{
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
  return numbers::NumberTraits<Number>::abs(value);
}


template <int dim, typename Number>
inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV
                                                  Tensor<0, dim, Number>::norm_square() const
{
  // We cannot use Assert inside a CUDA kernel
#ifndef __CUDA_ARCH__
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
#endif
  return numbers::NumberTraits<Number>::abs_square(value);
}


template <int dim, typename Number>
template <typename OtherNumber>
inline void
Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
                                         unsigned int &       index) const
{
  Assert(dim != 0,
         ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
  result[index] = value;
  ++index;
}


template <int dim, typename Number>
inline void
Tensor<0, dim, Number>::clear()
{
  // Some auto-differentiable numbers need explicit
  // zero initialization.
  value = internal::NumberType<Number>::value(0.0);
}


template <int dim, typename Number>
template <class Archive>
inline void
Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
{
  ar &value;
}


/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/


template <int rank_, int dim, typename Number>
inline DEAL_II_ALWAYS_INLINE
Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
Tensor<rank_, dim, Number>::Tensor(
  const Tensor<rank_, dim, OtherNumber> &initializer)
{
  for (unsigned int i = 0; i != dim; ++i)
    values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
Tensor<rank_, dim, Number>::Tensor(
  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] = initializer[i];
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
                                operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
{
  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
}



namespace internal
{
  namespace TensorSubscriptor
  {
    template <typename ArrayElementType, int dim>
    inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType &
                                                       subscript(ArrayElementType * values,
                                                                 const unsigned int i,
                                                                 std::integral_constant<int, dim>)
    {
      // We cannot use Assert in a CUDA kernel
#ifndef __CUDA_ARCH__
      Assert(i < dim, ExcIndexRange(i, 0, dim));
#endif
      return values[i];
    }


    template <typename ArrayElementType>
    ArrayElementType &
    subscript(ArrayElementType *,
              const unsigned int,
              std::integral_constant<int, 0>)
    {
      Assert(
        false,
        ExcMessage(
          "Cannot access elements of an object of type Tensor<rank,0,Number>."));
      static ArrayElementType t;
      return t;
    }
  } // namespace TensorSubscriptor
} // namespace internal


template <int rank_, int dim, typename Number>
inline DEAL_II_ALWAYS_INLINE                       DEAL_II_CUDA_HOST_DEV
  typename Tensor<rank_, dim, Number>::value_type &Tensor<rank_, dim, Number>::
                                                   operator[](const unsigned int i)
{
  return ::internal::TensorSubscriptor::subscript(
    values, i, std::integral_constant<int, dim>());
}


template <int rank_, int dim, typename Number>
constexpr DEAL_II_ALWAYS_INLINE
    DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
    Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
{
  return ::internal::TensorSubscriptor::subscript(
    values, i, std::integral_constant<int, dim>());
}


template <int rank_, int dim, typename Number>
inline const Number &Tensor<rank_, dim, Number>::
                     operator[](const TableIndices<rank_> &indices) const
{
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));

  return TensorAccessors::extract<rank_>(*this, indices);
}



template <int rank_, int dim, typename Number>
inline Number &Tensor<rank_, dim, Number>::
               operator[](const TableIndices<rank_> &indices)
{
  Assert(dim != 0,
         ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));

  return TensorAccessors::extract<rank_>(*this, indices);
}



template <int rank_, int dim, typename Number>
inline Number *
Tensor<rank_, dim, Number>::begin_raw()
{
  return std::addressof(
    this->operator[](this->unrolled_to_component_indices(0)));
}



template <int rank_, int dim, typename Number>
constexpr const Number *
Tensor<rank_, dim, Number>::begin_raw() const
{
  return std::addressof(
    this->operator[](this->unrolled_to_component_indices(0)));
}



template <int rank_, int dim, typename Number>
inline Number *
Tensor<rank_, dim, Number>::end_raw()
{
  return begin_raw() + n_independent_components;
}



template <int rank_, int dim, typename Number>
constexpr const Number *
Tensor<rank_, dim, Number>::end_raw() const
{
  return begin_raw() + n_independent_components;
}



template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator=(const Tensor<rank_, dim, OtherNumber> &t)
{
  if (dim > 0)
    std::copy(&t.values[0], &t.values[0] + dim, &values[0]);
  return *this;
}


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

  for (unsigned int i = 0; i < dim; ++i)
    values[i] = internal::NumberType<Number>::value(0.0);
  return *this;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline bool
Tensor<rank_, dim, Number>::
operator==(const Tensor<rank_, dim, OtherNumber> &p) const
{
  for (unsigned int i = 0; i < dim; ++i)
    if (values[i] != p.values[i])
      return false;
  return true;
}


// At some places in the library, we have Point<0> for formal reasons
// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
// in the above function that the loop end check always fails, we
// implement this function here
template <>
template <>
inline bool
Tensor<1, 0, double>::operator==(const Tensor<1, 0, double> &) const
{
  return true;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
constexpr bool
Tensor<rank_, dim, Number>::
operator!=(const Tensor<rank_, dim, OtherNumber> &p) const
{
  return !((*this) == p);
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator+=(const Tensor<rank_, dim, OtherNumber> &p)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] += p.values[i];
  return *this;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator-=(const Tensor<rank_, dim, OtherNumber> &p)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] -= p.values[i];
  return *this;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] *= s;
  return *this;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline Tensor<rank_, dim, Number> &
Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] /= s;
  return *this;
}


template <int rank_, int dim, typename Number>
inline Tensor<rank_, dim, Number>
Tensor<rank_, dim, Number>::operator-() const
{
  Tensor<rank_, dim, Number> tmp;

  for (unsigned int i = 0; i < dim; ++i)
    tmp.values[i] = -values[i];

  return tmp;
}


template <int rank_, int dim, typename Number>
inline typename numbers::NumberTraits<Number>::real_type
Tensor<rank_, dim, Number>::norm() const
{
  return std::sqrt(norm_square());
}


template <int rank_, int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
Tensor<rank_, dim, Number>::norm_square() const
{
  typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
    typename numbers::NumberTraits<Number>::real_type>::value(0.0);
  for (unsigned int i = 0; i < dim; ++i)
    s += values[i].norm_square();

  return s;
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline void
Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
{
  AssertDimension(result.size(),
                  (Utilities::fixed_power<rank_, unsigned int>(dim)));

  unsigned int index = 0;
  unroll_recursion(result, index);
}


template <int rank_, int dim, typename Number>
template <typename OtherNumber>
inline void
Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
                                             unsigned int &       index) const
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i].unroll_recursion(result, index);
}


template <int rank_, int dim, typename Number>
inline unsigned int
Tensor<rank_, dim, Number>::component_to_unrolled_index(
  const TableIndices<rank_> &indices)
{
  unsigned int index = 0;
  for (int r = 0; r < rank_; ++r)
    index = index * dim + indices[r];

  return index;
}



namespace internal
{
  // unrolled_to_component_indices is instantiated from DataOut for dim==0
  // and rank=2. Make sure we don't have compiler warnings.

  template <int dim>
  inline unsigned int
  mod(const unsigned int x)
  {
    return x % dim;
  }

  template <>
  inline unsigned int
  mod<0>(const unsigned int x)
  {
    Assert(false, ExcInternalError());
    return x;
  }

  template <int dim>
  inline unsigned int
  div(const unsigned int x)
  {
    return x / dim;
  }

  template <>
  inline unsigned int
  div<0>(const unsigned int x)
  {
    Assert(false, ExcInternalError());
    return x;
  }

} // namespace internal



template <int rank_, int dim, typename Number>
inline TableIndices<rank_>
Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
{
  Assert(i < n_independent_components,
         ExcIndexRange(i, 0, n_independent_components));

  TableIndices<rank_> indices;

  unsigned int remainder = i;
  for (int r = rank_ - 1; r >= 0; --r)
    {
      indices[r] = internal::mod<dim>(remainder);
      remainder  = internal::div<dim>(remainder);
    }
  Assert(remainder == 0, ExcInternalError());

  return indices;
}


template <int rank_, int dim, typename Number>
inline void
Tensor<rank_, dim, Number>::clear()
{
  for (unsigned int i = 0; i < dim; ++i)
    values[i] = internal::NumberType<Number>::value(0.0);
}


template <int rank_, int dim, typename Number>
constexpr std::size_t
Tensor<rank_, dim, Number>::memory_consumption()
{
  return sizeof(Tensor<rank_, dim, Number>);
}


template <int rank_, int dim, typename Number>
template <class Archive>
inline void
Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
{
  ar &values;
}


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


template <int rank_, int dim, typename Number>
inline std::ostream &
operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
{
  for (unsigned int i = 0; i < dim; ++i)
    {
      out << p[i];
      if (i != dim - 1)
        out << ' ';
    }

  return out;
}


template <int dim, typename Number>
inline std::ostream &
operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
{
  out << static_cast<const Number &>(p);
  return out;
}





template <int dim, typename Number, typename Other>
constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Other, Number>::type
operator*(const Other &object, const Tensor<0, dim, Number> &t)
{
  return object * static_cast<const Number &>(t);
}



template <int dim, typename Number, typename Other>
constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, Other>::type
operator*(const Tensor<0, dim, Number> &t, const Other &object)
{
  return static_cast<const Number &>(t) * object;
}


template <int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
operator*(const Tensor<0, dim, Number> &     src1,
          const Tensor<0, dim, OtherNumber> &src2)
{
  return static_cast<const Number &>(src1) *
         static_cast<const OtherNumber &>(src2);
}


template <int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  Tensor<0,
         dim,
         typename ProductType<Number,
                              typename EnableIfScalar<OtherNumber>::type>::type>
  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
{
  return static_cast<const Number &>(t) / factor;
}


template <int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
  operator+(const Tensor<0, dim, Number> &     p,
            const Tensor<0, dim, OtherNumber> &q)
{
  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
}


template <int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
  operator-(const Tensor<0, dim, Number> &     p,
            const Tensor<0, dim, OtherNumber> &q)
{
  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
  Tensor<rank,
         dim,
         typename ProductType<Number,
                              typename EnableIfScalar<OtherNumber>::type>::type>
  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
{
  // recurse over the base objects
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tt;
  for (unsigned int d = 0; d < dim; ++d)
    tt[d] = t[d] * factor;
  return tt;
}


template <int rank, int dim, typename Number, typename OtherNumber>
constexpr DEAL_II_ALWAYS_INLINE
  Tensor<rank,
         dim,
         typename ProductType<typename EnableIfScalar<Number>::type,
                              OtherNumber>::type>
  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
{
  // simply forward to the operator above
  return t * factor;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline Tensor<
  rank,
  dim,
  typename ProductType<Number,
                       typename EnableIfScalar<OtherNumber>::type>::type>
operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
{
  // recurse over the base objects
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tt;
  for (unsigned int d = 0; d < dim; ++d)
    tt[d] = t[d] / factor;
  return tt;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
  operator+(const Tensor<rank, dim, Number> &     p,
            const Tensor<rank, dim, OtherNumber> &q)
{
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tmp(p);

  for (unsigned int i = 0; i < dim; ++i)
    tmp[i] += q[i];

  return tmp;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
  operator-(const Tensor<rank, dim, Number> &     p,
            const Tensor<rank, dim, OtherNumber> &q)
{
  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tmp(p);

  for (unsigned int i = 0; i < dim; ++i)
    tmp[i] -= q[i];

  return tmp;
}





template <int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
inline 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 Tensor<rank_2, dim, OtherNumber> &src2)
{
  typename Tensor<rank_1 + rank_2 - 2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
    result;

  TensorAccessors::internal::
    ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
      reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);

  return result;
}


template <int index_1,
          int index_2,
          int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
  typename Tensor<rank_1 + rank_2 - 2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
  contract(const Tensor<rank_1, dim, Number> &     src1,
           const Tensor<rank_2, dim, OtherNumber> &src2)
{
  Assert(0 <= index_1 && index_1 < rank_1,
         ExcMessage(
           "The specified index_1 must lie within the range [0,rank_1)"));
  Assert(0 <= index_2 && index_2 < rank_2,
         ExcMessage(
           "The specified index_2 must lie within the range [0,rank_2)"));

  using namespace TensorAccessors;
  using namespace TensorAccessors::internal;

  // Reorder index_1 to the end of src1:
  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
    reord_01 = reordered_index_view<index_1, rank_1>(src1);

  // Reorder index_2 to the end of src2:
  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
    reord_02 = reordered_index_view<index_2, rank_2>(src2);

  typename Tensor<rank_1 + rank_2 - 2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
    result;
  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
  return result;
}


template <int index_1,
          int index_2,
          int index_3,
          int index_4,
          int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
inline
  typename Tensor<rank_1 + rank_2 - 4,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
  double_contract(const Tensor<rank_1, dim, Number> &     src1,
                  const Tensor<rank_2, dim, OtherNumber> &src2)
{
  Assert(0 <= index_1 && index_1 < rank_1,
         ExcMessage(
           "The specified index_1 must lie within the range [0,rank_1)"));
  Assert(0 <= index_3 && index_3 < rank_1,
         ExcMessage(
           "The specified index_3 must lie within the range [0,rank_1)"));
  Assert(index_1 != index_3,
         ExcMessage("index_1 and index_3 must not be the same"));
  Assert(0 <= index_2 && index_2 < rank_2,
         ExcMessage(
           "The specified index_2 must lie within the range [0,rank_2)"));
  Assert(0 <= index_4 && index_4 < rank_2,
         ExcMessage(
           "The specified index_4 must lie within the range [0,rank_2)"));
  Assert(index_2 != index_4,
         ExcMessage("index_2 and index_4 must not be the same"));

  using namespace TensorAccessors;
  using namespace TensorAccessors::internal;

  // Reorder index_1 to the end of src1:
  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
    reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);

  // Reorder index_2 to the end of src2:
  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
    reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);

  // Now, reorder index_3 to the end of src1. We have to make sure to
  // preserve the original ordering: index_1 has been removed. If
  // index_3 > index_1, we have to use (index_3 - 1) instead:
  ReorderedIndexView<
    (index_3 < index_1 ? index_3 : index_3 - 1),
    rank_1,
    ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
    reord_3 =
      TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
                                                                  index_3 - 1,
    rank_1 > (reord_1);

  // Now, reorder index_4 to the end of src2. We have to make sure to
  // preserve the original ordering: index_2 has been removed. If
  // index_4 > index_2, we have to use (index_4 - 1) instead:
  ReorderedIndexView<
    (index_4 < index_2 ? index_4 : index_4 - 1),
    rank_2,
    ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
    reord_4 =
      TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
                                                                  index_4 - 1,
    rank_2 > (reord_2);

  typename Tensor<rank_1 + rank_2 - 4,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
    result;
  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
  return result;
}


template <int rank, int dim, typename Number, typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
scalar_product(const Tensor<rank, dim, Number> &     left,
               const Tensor<rank, dim, OtherNumber> &right)
{
  typename ProductType<Number, OtherNumber>::type result;
  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
  return result;
}


template <template <int, int, typename> class TensorT1,
          template <int, int, typename> class TensorT2,
          template <int, int, typename> class TensorT3,
          int rank_1,
          int rank_2,
          int dim,
          typename T1,
          typename T2,
          typename T3>
typename ProductType<T1, typename ProductType<T2, T3>::type>::type
contract3(const TensorT1<rank_1, dim, T1> &         left,
          const TensorT2<rank_1 + rank_2, dim, T2> &middle,
          const TensorT3<rank_2, dim, T3> &         right)
{
  using return_type =
    typename ProductType<T1, typename ProductType<T2, T3>::type>::type;
  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
                                                                      middle,
                                                                      right);
}


template <int rank_1,
          int rank_2,
          int dim,
          typename Number,
          typename OtherNumber>
inline DEAL_II_ALWAYS_INLINE
  Tensor<rank_1 + rank_2, dim, typename ProductType<Number, OtherNumber>::type>
  outer_product(const Tensor<rank_1, dim, Number> &     src1,
                const Tensor<rank_2, dim, OtherNumber> &src2)
{
  typename Tensor<rank_1 + rank_2,
                  dim,
                  typename ProductType<Number, OtherNumber>::type>::tensor_type
    result;
  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
  return result;
}





template <int dim, typename Number>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
                             cross_product_2d(const Tensor<1, dim, Number> &src)
{
  Assert(dim == 2, ExcInternalError());

  Tensor<1, dim, Number> result;

  result[0] = src[1];
  result[1] = -src[0];

  return result;
}


template <int dim, typename Number>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
                             cross_product_3d(const Tensor<1, dim, Number> &src1,
                                              const Tensor<1, dim, Number> &src2)
{
  Assert(dim == 3, ExcInternalError());

  Tensor<1, dim, Number> result;

  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
  result[2] = src1[0] * src2[1] - src1[1] * src2[0];

  return result;
}





template <int dim, typename Number>
inline Number
determinant(const Tensor<2, dim, Number> &t)
{
  // Compute the determinant using the Laplace expansion of the
  // determinant. We expand along the last row.
  Number det = internal::NumberType<Number>::value(0.0);

  for (unsigned int k = 0; k < dim; ++k)
    {
      Tensor<2, dim - 1, Number> minor;
      for (unsigned int i = 0; i < dim - 1; ++i)
        for (unsigned int j = 0; j < dim - 1; ++j)
          minor[i][j] = t[i][j < k ? j : j + 1];

      const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);

      det += t[dim - 1][k] * cofactor;
    }

  return ((dim % 2 == 0) ? 1. : -1.) * det;
}

template <typename Number>
constexpr Number
determinant(const Tensor<2, 1, Number> &t)
{
  return t[0][0];
}


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


template <int dim, typename Number>
inline Tensor<2, dim, Number>
invert(const Tensor<2, dim, Number> &)
{
  Number return_tensor[dim][dim];

  // if desired, take over the
  // inversion of a 4x4 tensor
  // from the FullMatrix
  AssertThrow(false, ExcNotImplemented());

  return Tensor<2, dim, Number>(return_tensor);
}


#ifndef DOXYGEN

template <typename Number>
inline Tensor<2, 1, Number>
invert(const Tensor<2, 1, Number> &t)
{
  Number return_tensor[1][1];

  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);

  return Tensor<2, 1, Number>(return_tensor);
}


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

  // this is Maple output,
  // thus a bit unstructured
  const Number inv_det_t = internal::NumberType<Number>::value(
    1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
  return_tensor[0][0] = t[1][1];
  return_tensor[0][1] = -t[0][1];
  return_tensor[1][0] = -t[1][0];
  return_tensor[1][1] = t[0][0];
  return_tensor *= inv_det_t;

  return return_tensor;
}


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

  const Number t4  = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
               t6  = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
               t8  = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
               t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
               t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
               t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
               inv_det_t = internal::NumberType<Number>::value(
                 1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
                        t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
                        internal::NumberType<Number>::value(t[1][2] * t[2][1]);
  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
                        internal::NumberType<Number>::value(t[0][1] * t[2][2]);
  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
                        internal::NumberType<Number>::value(t[0][2] * t[1][1]);
  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
                        internal::NumberType<Number>::value(t[1][0] * t[2][2]);
  return_tensor[1][1] =
    internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
  return_tensor[1][2] = t00 - t6;
  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
                        internal::NumberType<Number>::value(t[1][1] * t[2][0]);
  return_tensor[2][1] =
    t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
  return_tensor *= inv_det_t;

  return return_tensor;
}

#endif /* DOXYGEN */


template <int dim, typename Number>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
                             transpose(const Tensor<2, dim, Number> &t)
{
  Tensor<2, dim, Number> tt;
  for (unsigned int i = 0; i < dim; ++i)
    {
      tt[i][i] = t[i][i];
      for (unsigned int j = i + 1; j < dim; ++j)
        {
          tt[i][j] = t[j][i];
          tt[j][i] = t[i][j];
        };
    }
  return tt;
}


template <int dim, typename Number>
constexpr Tensor<2, dim, Number>
adjugate(const Tensor<2, dim, Number> &t)
{
  return determinant(t) * invert(t);
}


template <int dim, typename Number>
constexpr Tensor<2, dim, Number>
cofactor(const Tensor<2, dim, Number> &t)
{
  return transpose(adjugate(t));
}


template <int dim, typename Number>
inline Number
l1_norm(const Tensor<2, dim, Number> &t)
{
  Number max = internal::NumberType<Number>::value(0.0);
  for (unsigned int j = 0; j < dim; ++j)
    {
      Number sum = internal::NumberType<Number>::value(0.0);
      for (unsigned int i = 0; i < dim; ++i)
        sum += std::fabs(t[i][j]);

      if (sum > max)
        max = sum;
    }

  return max;
}


template <int dim, typename Number>
inline Number
linfty_norm(const Tensor<2, dim, Number> &t)
{
  Number max = internal::NumberType<Number>::value(0.0);
  for (unsigned int i = 0; i < dim; ++i)
    {
      Number sum = internal::NumberType<Number>::value(0.0);
      for (unsigned int j = 0; j < dim; ++j)
        sum += std::fabs(t[i][j]);

      if (sum > max)
        max = sum;
    }

  return max;
}



#ifndef DOXYGEN


#  ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING

// Specialization of functions for ADOL-C number types when
// the advanced branching feature is used
template <int dim>
inline adouble
l1_norm(const Tensor<2, dim, adouble> &t)
{
  adouble max = internal::NumberType<adouble>::value(0.0);
  for (unsigned int j = 0; j < dim; ++j)
    {
      adouble sum = internal::NumberType<adouble>::value(0.0);
      for (unsigned int i = 0; i < dim; ++i)
        sum += std::fabs(t[i][j]);

      condassign(max, (sum > max), sum, max);
    }

  return max;
}


template <int dim>
inline adouble
linfty_norm(const Tensor<2, dim, adouble> &t)
{
  adouble max = internal::NumberType<adouble>::value(0.0);
  for (unsigned int i = 0; i < dim; ++i)
    {
      adouble sum = internal::NumberType<adouble>::value(0.0);
      for (unsigned int j = 0; j < dim; ++j)
        sum += std::fabs(t[i][j]);

      condassign(max, (sum > max), sum, max);
    }

  return max;
}

#  endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

// include deprecated non-member functions operating on Tensor
#include <deal.II/base/tensor_deprecated.h>

#endif