tensor_accessors.h

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

#ifndef dealii_tensor_accessors_h
#define dealii_tensor_accessors_h

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

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


DEAL_II_NAMESPACE_OPEN

namespace TensorAccessors
{
  // forward declarations
  namespace internal
  {
    template <int index, int rank, typename T>
    class ReorderedIndexView;
    template <int position, int rank>
    struct ExtractHelper;
    template <int no_contr, int rank_1, int rank_2, int dim>
    class Contract;
    template <int rank_1, int rank_2, int dim>
    class Contract3;
  } // namespace internal


  template <typename T>
  struct ValueType
  {
    using value_type = typename T::value_type;
  };

  template <typename T>
  struct ValueType<const T>
  {
    using value_type = const typename T::value_type;
  };

  template <typename T, std::size_t N>
  struct ValueType<T[N]>
  {
    using value_type = T;
  };

  template <typename T, std::size_t N>
  struct ValueType<const T[N]>
  {
    using value_type = const T;
  };


  template <int deref_steps, typename T>
  struct ReturnType
  {
    using value_type =
      typename ReturnType<deref_steps - 1,
                          typename ValueType<T>::value_type>::value_type;
  };

  template <typename T>
  struct ReturnType<0, T>
  {
    using value_type = T;
  };


  template <int index, int rank, typename T>
  inline DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView<index, rank, T>
                               reordered_index_view(T &t)
  {
    static_assert(0 <= index && index < rank,
                  "The specified index must lie within the range [0,rank)");

    return internal::ReorderedIndexView<index, rank, T>(t);
  }


  template <int rank, typename T, typename ArrayType>
  typename ReturnType<rank, T>::value_type &
  extract(T &t, const ArrayType &indices)
  {
    return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(
      t, indices);
  }


  template <int no_contr,
            int rank_1,
            int rank_2,
            int dim,
            typename T1,
            typename T2,
            typename T3>
  inline DEAL_II_ALWAYS_INLINE void
  contract(T1 &result, const T2 &left, const T3 &right)
  {
    static_assert(rank_1 >= no_contr,
                  "The rank of the left tensor must be "
                  "equal or greater than the number of "
                  "contractions");
    static_assert(rank_2 >= no_contr,
                  "The rank of the right tensor must be "
                  "equal or greater than the number of "
                  "contractions");

    internal::Contract<no_contr, rank_1, rank_2, dim>::
      template contract<T1, T2, T3>(result, left, right);
  }


  template <int rank_1,
            int rank_2,
            int dim,
            typename T1,
            typename T2,
            typename T3,
            typename T4>
  T1
  contract3(const T2 &left, const T3 &middle, const T4 &right)
  {
    return internal::Contract3<rank_1, rank_2, dim>::
      template contract3<T1, T2, T3, T4>(left, middle, right);
  }


  namespace internal
  {
    // -------------------------------------------------------------------------
    // Forward declarations and type traits
    // -------------------------------------------------------------------------

    template <int rank, typename S>
    class StoreIndex;
    template <typename T>
    class Identity;
    template <int no_contr, int dim>
    class Contract2;

    template <typename T>
    struct ReferenceType
    {
      using type = T &;
    };

    template <int rank, typename S>
    struct ReferenceType<StoreIndex<rank, S>>
    {
      using type = StoreIndex<rank, S>;
    };

    template <int index, int rank, typename T>
    struct ReferenceType<ReorderedIndexView<index, rank, T>>
    {
      using type = ReorderedIndexView<index, rank, T>;
    };


    // TODO: Is there a possibility to just have the following block of
    // explanation on an internal page in doxygen? If, yes. Doxygen
    // wizards, your call!

    // -------------------------------------------------------------------------
    // Implementation of helper classes for reordered_index_view
    // -------------------------------------------------------------------------

    // OK. This is utterly brutal template magic. Therefore, we will not
    // comment on the individual internal helper classes, because this is
    // of not much value, but explain the general recursion procedure.
    //
    // (In order of appearance)
    //
    // Our task is to reorder access to a tensor object where a specified
    // index is moved to the end. Thus we want to construct an object
    // <code>reordered</code> out of a <code>tensor</code> where the
    // following access patterns are equivalent:
    // @code
    //   tensor    [i_0]...[i_index-1][i_index][i_index+1]...[i_n]
    //   reordered [i_0]...[i_index_1][i_index+1]...[i_n][i_index]
    // @endcode
    //
    // The first task is to get rid of the application of
    // [i_0]...[i_index-1]. This is a classical recursion pattern - relay
    // the task from <index, rank> to <index-1, rank-1> by accessing the
    // subtensor object:

    template <int index, int rank, typename T>
    class ReorderedIndexView
    {
    public:
      ReorderedIndexView(typename ReferenceType<T>::type t)
        : t_(t)
      {}

      using value_type = ReorderedIndexView<index - 1,
                                            rank - 1,
                                            typename ValueType<T>::value_type>;

      // Recurse by applying index j directly:
      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
      {
        return value_type(t_[j]);
      }

    private:
      typename ReferenceType<T>::type t_;
    };

    // At some point we hit the condition index == 0 and rank > 1, i.e.,
    // the first index should be reordered to the end.
    //
    // At this point we cannot be lazy any more and have to start storing
    // indices because we get them in the wrong order. The user supplies
    //   [i_0][i_1]...[i_{rank - 1}]
    // but we have to call the subtensor object with
    //   [i_{rank - 1}[i_0][i_1]...[i_{rank-2}]
    //
    // So give up and relay the task to the StoreIndex class:

    template <int rank, typename T>
    class ReorderedIndexView<0, rank, T>
    {
    public:
      ReorderedIndexView(typename ReferenceType<T>::type t)
        : t_(t)
      {}

      using value_type = StoreIndex<rank - 1, internal::Identity<T>>;

      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
      {
        return value_type(Identity<T>(t_), j);
      }

    private:
      typename ReferenceType<T>::type t_;
    };

    // Sometimes, we're lucky and don't have to do anything. In this case
    // just return the original tensor.

    template <typename T>
    class ReorderedIndexView<0, 1, T>
    {
    public:
      ReorderedIndexView(typename ReferenceType<T>::type t)
        : t_(t)
      {}

      using value_type =
        typename ReferenceType<typename ValueType<T>::value_type>::type;

      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
      {
        return t_[j];
      }

    private:
      typename ReferenceType<T>::type t_;
    };

    // Here, Identity is a helper class to ground the recursion in
    // StoreIndex. Its implementation is easy - we haven't stored any
    // indices yet. So, we just provide a function apply that returns the
    // application of an index j to the stored tensor t_:

    template <typename T>
    class Identity
    {
    public:
      Identity(typename ReferenceType<T>::type t)
        : t_(t)
      {}

      using return_type = typename ValueType<T>::value_type;

      inline DEAL_II_ALWAYS_INLINE typename ReferenceType<return_type>::type
      apply(unsigned int j) const
      {
        return t_[j];
      }

    private:
      typename ReferenceType<T>::type t_;
    };

    // StoreIndex is a class that stores an index recursively with every
    // invocation of operator[](unsigned int j): We do this by recursively
    // creating a new StoreIndex class of lower rank that stores the
    // supplied index j and holds a copy of the current class (with all
    // other stored indices). Again, we provide an apply member function
    // that knows how to apply an index on the highest rank and all
    // subsequently stored indices:

    template <int rank, typename S>
    class StoreIndex
    {
    public:
      StoreIndex(S s, int i)
        : s_(s)
        , i_(i)
      {}

      using value_type = StoreIndex<rank - 1, StoreIndex<rank, S>>;

      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
      {
        return value_type(*this, j);
      }

      using return_type =
        typename ValueType<typename S::return_type>::value_type;

      inline typename ReferenceType<return_type>::type
      apply(unsigned int j) const
      {
        return s_.apply(j)[i_];
      }

    private:
      const S   s_;
      const int i_;
    };

    // We have to store indices until we hit rank == 1. Then, upon the next
    // invocation of operator[](unsigned int j) we have all necessary
    // information available to return the actual object.

    template <typename S>
    class StoreIndex<1, S>
    {
    public:
      StoreIndex(S s, int i)
        : s_(s)
        , i_(i)
      {}

      using return_type =
        typename ValueType<typename S::return_type>::value_type;
      using value_type = return_type;

      inline DEAL_II_ALWAYS_INLINE return_type &operator[](unsigned int j) const
      {
        return s_.apply(j)[i_];
      }

    private:
      const S   s_;
      const int i_;
    };


    // -------------------------------------------------------------------------
    // Implementation of helper classes for extract
    // -------------------------------------------------------------------------

    // Straightforward recursion implemented by specializing ExtractHelper
    // for position == rank. We use the type trait ReturnType<rank, T> to
    // have an idea what the final type will be.
    template <int position, int rank>
    struct ExtractHelper
    {
      template <typename T, typename ArrayType>
      inline static typename ReturnType<rank - position, T>::value_type &
      extract(T &t, const ArrayType &indices)
      {
        return ExtractHelper<position + 1, rank>::template extract<
          typename ValueType<T>::value_type,
          ArrayType>(t[indices[position]], indices);
      }
    };

    // For position == rank there is nothing to extract, just return the
    // object.
    template <int rank>
    struct ExtractHelper<rank, rank>
    {
      template <typename T, typename ArrayType>
      inline static T &
      extract(T &t, const ArrayType &)
      {
        return t;
      }
    };


    // -------------------------------------------------------------------------
    // Implementation of helper classes for contract
    // -------------------------------------------------------------------------

    // Straightforward recursive pattern:
    //
    // As long as rank_1 > no_contr, assign indices from the left tensor to
    // result. This builds up the first part of the nested outer loops:
    //
    // for(unsigned int i_0; i_0 < dim; ++i_0)
    //   ...
    //     for(i_; i_ < dim; ++i_)
    //       [...]
    //         result[i_0]..[i_] ... left[i_0]..[i_] ...

    template <int no_contr, int rank_1, int rank_2, int dim>
    class Contract
    {
    public:
      template <typename T1, typename T2, typename T3>
      inline DEAL_II_ALWAYS_INLINE static void
      contract(T1 &result, const T2 &left, const T3 &right)
      {
        for (unsigned int i = 0; i < dim; ++i)
          Contract<no_contr, rank_1 - 1, rank_2, dim>::contract(result[i],
                                                                left[i],
                                                                right);
      }
    };

    // If rank_1 == no_contr leave out the remaining no_contr indices for
    // the contraction and assign indices from the right tensor to the
    // result. This builds up the second part of the nested loops:
    //
    //  for(unsigned int i_0 = 0; i_0 < dim; ++i_0)
    //    ...
    //      for(unsigned int i_ = 0; i_ < dim; ++i_)
    //        for(unsigned int j_0 = 0; j_0 < dim; ++j_0)
    //          ...
    //            for(unsigned int j_ = 0; j_ < dim; ++j_)
    //             [...]
    //               result[i_0]..[i_][j_0]..[j_] ... left[i_0]..[i_] ...
    //               right[j_0]..[j_]
    //

    template <int no_contr, int rank_2, int dim>
    class Contract<no_contr, no_contr, rank_2, dim>
    {
    public:
      template <typename T1, typename T2, typename T3>
      inline DEAL_II_ALWAYS_INLINE static void
      contract(T1 &result, const T2 &left, const T3 &right)
      {
        for (unsigned int i = 0; i < dim; ++i)
          Contract<no_contr, no_contr, rank_2 - 1, dim>::contract(result[i],
                                                                  left,
                                                                  right[i]);
      }
    };

    // If rank_1 == rank_2 == no_contr we have built up all of the outer
    // loop. Now, it is time to do the actual contraction:
    //
    // [...]
    //   {
    //     result[i_0]..[i_][j_0]..[j_] = 0.;
    //     for(unsigned int k_0 = 0; k_0 < dim; ++k_0)
    //       ...
    //         for(unsigned int k_ = 0; k_ < dim; ++k_)
    //           result[i_0]..[i_][j_0]..[j_] += left[i_0]..[i_][k_0]..[k_] *
    //           right[j_0]..[j_][k_0]..[k_];
    //   }
    //
    //  Relay this summation to another helper class.

    template <int no_contr, int dim>
    class Contract<no_contr, no_contr, no_contr, dim>
    {
    public:
      template <typename T1, typename T2, typename T3>
      inline DEAL_II_ALWAYS_INLINE static void
      contract(T1 &result, const T2 &left, const T3 &right)
      {
        result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
      }
    };

    // Straightforward recursion:
    //
    // Contract leftmost index and recurse one down.

    template <int no_contr, int dim>
    class Contract2
    {
    public:
      template <typename T1, typename T2, typename T3>
      inline DEAL_II_ALWAYS_INLINE static T1
      contract2(const T2 &left, const T3 &right)
      {
        // Some auto-differentiable numbers need explicit
        // zero initialization.
        T1 result = ::internal::NumberType<T1>::value(0.0);
        for (unsigned int i = 0; i < dim; ++i)
          result +=
            Contract2<no_contr - 1, dim>::template contract2<T1>(left[i],
                                                                 right[i]);
        return result;
      }
    };

    // A contraction of two objects of order 0 is just a scalar
    // multiplication:

    template <int dim>
    class Contract2<0, dim>
    {
    public:
      template <typename T1, typename T2, typename T3>
      inline DEAL_II_ALWAYS_INLINE static T1
      contract2(const T2 &left, const T3 &right)
      {
        return left * right;
      }
    };


    // -------------------------------------------------------------------------
    // Implementation of helper classes for contract3
    // -------------------------------------------------------------------------

    // Fully contract three tensorial objects
    //
    // As long as rank_1 > 0, recurse over left and middle:
    //
    // for(unsigned int i_0; i_0 < dim; ++i_0)
    //   ...
    //     for(i_; i_ < dim; ++i_)
    //       [...]
    //         left[i_0]..[i_] ... middle[i_0]..[i_] ... right

    template <int rank_1, int rank_2, int dim>
    class Contract3
    {
    public:
      template <typename T1, typename T2, typename T3, typename T4>
      static inline T1
      contract3(const T2 &left, const T3 &middle, const T4 &right)
      {
        // Some auto-differentiable numbers need explicit
        // zero initialization.
        T1 result = ::internal::NumberType<T1>::value(0.0);
        for (unsigned int i = 0; i < dim; ++i)
          result += Contract3<rank_1 - 1, rank_2, dim>::template contract3<T1>(
            left[i], middle[i], right);
        return result;
      }
    };

    // If rank_1 ==0, continue to recurse over middle and right:
    //
    // for(unsigned int i_0; i_0 < dim; ++i_0)
    //   ...
    //     for(i_; i_ < dim; ++i_)
    //       for(unsigned int j_0; j_0 < dim; ++j_0)
    //         ...
    //           for(j_; j_ < dim; ++j_)
    //             [...]
    //               left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ...
    //               right[j_0]..[j_]

    template <int rank_2, int dim>
    class Contract3<0, rank_2, dim>
    {
    public:
      template <typename T1, typename T2, typename T3, typename T4>
      static inline T1
      contract3(const T2 &left, const T3 &middle, const T4 &right)
      {
        // Some auto-differentiable numbers need explicit
        // zero initialization.
        T1 result = ::internal::NumberType<T1>::value(0.0);
        for (unsigned int i = 0; i < dim; ++i)
          result +=
            Contract3<0, rank_2 - 1, dim>::template contract3<T1>(left,
                                                                  middle[i],
                                                                  right[i]);
        return result;
      }
    };

    // Contraction of three tensorial objects of rank 0 is just a scalar
    // multiplication.

    template <int dim>
    class Contract3<0, 0, dim>
    {
    public:
      template <typename T1, typename T2, typename T3, typename T4>
      static inline T1
      contract3(const T2 &left, const T3 &middle, const T4 &right)
      {
        return left * middle * right;
      }
    };

    // -------------------------------------------------------------------------

  } /* namespace internal */
} /* namespace TensorAccessors */

DEAL_II_NAMESPACE_CLOSE

#endif /* dealii_tensor_accessors_h */