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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_tensor_accessors_h
#define dealii_tensor_accessors_h

#include <deal.II/base/config.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/table_indices.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;
  }


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

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

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

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


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

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


  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
    {
      typedef T &type;
    };

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

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


    // 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) {}

      typedef ReorderedIndexView<index - 1, rank - 1, typename ValueType<T>::value_type>
      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) {}

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

      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) {}

      typedef typename ReferenceType<typename ValueType<T>::value_type>::type value_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) {}

      typedef typename ValueType<T>::value_type return_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) {}

      typedef StoreIndex<rank - 1, StoreIndex<rank, S> > value_type;

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

      typedef typename ValueType<typename S::return_type>::value_type return_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) {}

      typedef typename ValueType<typename S::return_type>::value_type return_type;
      typedef return_type value_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 */