tensorbase.hh

#ifndef __DUNE_ACFEM_TENSORS_TENSORBASE_HH__
#define __DUNE_ACFEM_TENSORS_TENSORBASE_HH__
 
#include <dune/common/fvector.hh>
#include <dune/common/tupleutility.hh>
 
#include "../common/types.hh"
#include "../mpl/generators.hh"
#include "../mpl/subtuple.hh"
#include "../mpl/filter.hh"
#include "../mpl/compare.hh"
#include "../mpl/multiindex.hh"
#include "../expressions/constantoperations.hh"
#include "../expressions/interface.hh"
 
#include "operations/restrictionoperators.hh"
 
#define DUNE_ACFEM_TENSOR_WORKAROUND_GCC(MAJOR) DUNE_ACFEM_IS_GCC(MAJOR, MAJOR)
 
namespace Dune
{
  namespace ACFem
  {
    namespace Tensor
    {
      using namespace Literals;
      using namespace Expressions;
 
      template<std::size_t... Ints>
      using Seq = IndexSequence<Ints...>;
 
      template<std::size_t... Dimensions>
      using TensorSignature = Seq<Dimensions...>;
 
      template<class T, class SFINAE = void>
      struct IsTensor
        : FalseType
      {};
 
      template<class T>
      struct IsTensor<T, std::enable_if_t<!IsDecay<T>::value> >
        : IsTensor<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsTensor<T, std::enable_if_t<(IsDecay<T>::value
                                           && T::Signature::size() == T::rank
                                          )> >
        : TrueType
      {};
 
      template<class T,
               std::enable_if_t<(IsTensor<T>::value
                                 && !IsTypedValue<T>::value
                                 && !Expressions::IsClosure<std::decay_t<T> >::value
                 ), int> = 0>
      decltype(auto) tensor(T&& t)
      {
        return std::forward<T>(t);
      }
 
      template<class T,
               std::enable_if_t<(IsTensor<T>::value
                                 && IsTypedValue<T>::value
                                 && !Expressions::IsClosure<std::decay_t<T> >::value
                 ), int> = 0>
      std::decay_t<T> tensor(T&& t)
      {
        return std::forward<T>(t);
      }
 
      template<
        class T,
        std::enable_if_t<(IsTensor<T>::value
                          && Expressions::IsClosure<std::decay_t<T> >::value
                          && !IsTypedValue<Operand<0, T> >::value
          ), int> = 0>
      constexpr decltype(auto) tensor(T&& t)
      {
        Expressions::danglingReferenceCheck<0>(std::forward<T>(t));
        return operand<0>(std::forward<T>(t));
      }
 
      template<
        class T,
        std::enable_if_t<(IsTensor<T>::value
                          && Expressions::IsClosure<std::decay_t<T> >::value
                          && IsTypedValue<Operand<0, T> >::value
          ), int> = 0>
      std::decay_t<Operand<0, T> > tensor(T&& t)
      {
        return std::forward<T>(t).template operand<0>();
      }
 
      template<class T, class SFINAE = void>
      struct IsTensorOperand
        : FalseType
      {};
 
      template<class T>
      struct IsTensorOperand<T, std::enable_if_t<!IsDecay<T>::value> >
        : IsTensorOperand<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsTensorOperand<T, std::enable_if_t<(IsTensor<T>::value && IsDecay<T>::value)> >
        : TrueType
      {};
 
      // forward
      template<class T, class SFINAE = void>
      struct TensorTraits;
 
      template<class T, std::enable_if_t<(IsTensorOperand<T>::value && !IsTensor<T>::value), int> = 0>
      auto tensor(T&& t)
      {
        return TensorTraits<T>::toTensor(std::forward<T>(t));
      }
 
      template<class Field, class Seq, class Impl>
      struct TensorBase;
 
      // Forward
      template<class T>
      class RestrictionOperators;
 
      template<class Field, std::size_t... Dimensions, class Impl>
      struct TensorBase<Field, TensorSignature<Dimensions...>, Impl>
        : public RestrictionOperators<Impl>
      {
       public:
        using TensorType = Impl;
        using FieldType = Field;
        using Signature = TensorSignature<Dimensions...>;
 
        static constexpr std::size_t rank = sizeof...(Dimensions);
 
        template<std::size_t... Ind, std::enable_if_t<(sizeof...(Ind) <= rank), int> = 0>
        static std::size_t constexpr dim(Seq<Ind...> = Seq<Ind...>{})
        {
          using SubSignature = ConditionalType<sizeof...(Ind) == 0, Signature, SubSequence<Signature, Ind...> >;
          return multiDim(SubSignature{});
        }
 
        template<std::size_t... Ind, std::enable_if_t<(sizeof...(Ind) > rank), int> = 0>
        static std::size_t constexpr dim(Seq<Ind...> = Seq<Ind...>{})
        {
          return std::numeric_limits<std::size_t>::max();
        }
 
        static auto constexpr signature()
        {
          return Signature{};
        }
 
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          return false;
        }
 
        std::string name() const
        {
          return name(Signature{});
        }
 
        template<class T, std::enable_if_t<(rank == 0 && SameDecay<T, FieldType>::value), int> = 0>
        operator T () const
        {
          return static_cast<const Impl&>(*this)();
        }
 
       protected:
        std::string name(Seq<>) const
        {
          return std::to_string(static_cast<const Impl&>(*this)());
        }
 
        template<std::size_t D0, std::size_t... Rest>
        std::string name(Seq<D0, Rest...> sig) const
        {
          return "tensor<"+toString(sig)+">";
        }
 
      };
 
      template<class T, class SFINAE>
      struct TensorTraits
      {
        using Signature = Seq<>;
        static constexpr std::size_t rank = std::numeric_limits<std::size_t>::max();
      };
 
      template<class T>
      constexpr inline bool HasTensorTraitsV =
        TensorTraits<T>::Signature::size() == TensorTraits<T>::rank;
 
      template<class T>
      struct TensorTraitsBase
      {
        using TensorType = std::decay_t<T>;
        using FieldType = typename TensorType::FieldType;
        using Signature = typename TensorType::Signature;
        static constexpr std::size_t rank = TensorType::rank;
 
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          return TensorType::isZero(Seq<Indices...>{}, Pos{});
        }
 
        template<std::size_t... Ind, std::enable_if_t<(sizeof...(Ind) <= rank), int> = 0>
        static std::size_t constexpr dim(Seq<Ind...> = Seq<Ind...>{})
        {
          using SubSignature = ConditionalType<sizeof...(Ind) == 0, Signature, SubSequence<Signature, Ind...> >;
          return multiDim(SubSignature{});
        }
 
        template<std::size_t... Ind, std::enable_if_t<(sizeof...(Ind) > rank), int> = 0>
        static std::size_t constexpr dim(Seq<Ind...> = Seq<Ind...>{})
        {
          return std::numeric_limits<std::size_t>::max();
        }
 
        static Signature constexpr signature()
        {
          return Signature{};
        }
 
        using ZeroIndex = MakeIndexSequence<rank, 0, 0>;
       private:
        static TrueType  mutableArg(FieldType&);
        static FalseType mutableArg(const FieldType&);
        static FalseType mutableArg(FieldType&&);
#ifndef NDEBUG
        static_assert(ZeroIndex::size() == rank, "bug");
#endif
        template<class U, class SFINE = void>
        struct ElementZeroHelper
        {
          using Type = decltype(std::declval<const U&>()(ZeroIndex{}));
          static constexpr bool hasMutableCallOperator = false;
        };
 
        template<class U>
        struct ElementZeroHelper<U, VoidType<decltype(std::declval<U&>()(ZeroIndex{}))> >
        {
          using Type = decltype(std::declval<U&>()(ZeroIndex{}));
          static constexpr bool hasMutableCallOperator = true;
        };
       public:
        using ElementZero = typename ElementZeroHelper<T>::Type;
        static constexpr bool hasMutableComponents = decltype(mutableArg(std::declval<ElementZero>()))::value;
        static constexpr bool hasMutableCallOperator = ElementZeroHelper<T>::hasMutableCallOperator;
      };
 
      template<class T>
      struct TensorTraits<T, std::enable_if_t<!IsDecay<T>::value> >
        : TensorTraits<std::decay_t<T> >
      {};
 
      template<class T>
      struct TensorTraits<T, std::enable_if_t<IsTensor<T>::value && IsDecay<T>::value> >
        : TensorTraitsBase<T>
      {};
 
      template<class T, std::enable_if_t<IsTensorOperand<T>::value, int> = 0>
      constexpr auto rank(T&& t)
      {
        return TensorTraits<T>::rank;
      }
 
      template<class T, std::enable_if_t<IsTensorOperand<T>::value, int> = 0>
      constexpr auto signature(T&& t)
      {
        return typename TensorTraits<T>::Signature{};
      }
 
      template<class T, class SFINAE = void>
      struct IsTensorZero
        : FalseType
      {};
 
      template<class T>
      struct IsTensorZero<
        T,
        std::enable_if_t<(IsTensor<T>::value
                          && IsDecay<T>::value
                          && T::isZero()
        )> >
        : TrueType
      {};
 
      template<class T, class SFINAE = void>
      struct IsTensorNotZero
        : FalseType
      {};
 
      template<class T>
      struct IsTensorNotZero<
        T,
        std::enable_if_t<(IsTensor<T>::value
                          && IsDecay<T>::value
                          && !T::isZero()
        )> >
        : TrueType
      {};
 
      namespace {
 
        template<class Tensor, class Tuple, std::size_t... Unpack>
        decltype(auto) tensorValueHelper(Tensor&& t,
                                         Tuple&& indices,
                                         IndexSequence<Unpack...>&&)
        {
          return std::forward<Tensor>(t)(std::get<Unpack>(indices)...);
        }
      }
 
      template<
        class Tensor,
        class Tuple,
        std::enable_if_t<(IsTensor<Tensor>::value
                          && IsTupleLike<Tuple>::value
        ), int> = 0>
      decltype(auto) tensorValue(Tensor&& tensor, Tuple&& indices)
      {
        return tensorValueHelper(std::forward<Tensor>(tensor), std::forward<Tuple>(indices),
                                 MakeIndexSequence<TensorTraits<Tensor>::rank>{});
      }
 
      template<
        class Tensor,
        std::size_t... Indices,
        std::enable_if_t<(IsTensor<Tensor>::value
                          && sizeof...(Indices) == TensorTraits<Tensor>::rank
        ), int> = 0>
      decltype(auto) tensorValue(Tensor&& tensor, IndexSequence<Indices...>)
      {
        // could check whether Indices... is in range.
        return std::forward<Tensor>(tensor)(IndexSequence<Indices...>{});
      }
 
      template<
        class T, class Indices,
        std::enable_if_t<(IsTensorOperand<T>::value
                          && !IsTensor<T>::value
        ), int> = 0>
      decltype(auto) tensorValue(T&& t, Indices&& indices)
      {
        return tensorValue(tensor(std::forward<T>(t)), std::forward<Indices>(indices));
      }
 
      template<class T, std::enable_if_t<!Expressions::IsClosure<T>::value && IsTensorOperand<T>::value, int> = 0>
      constexpr auto name(T&& t)
      {
        return tensor(std::forward<T>(t)).name();
      }
 
      template<class T, std::enable_if_t<Expressions::IsClosure<T>::value && IsTensor<T>::value, int> = 0>
      constexpr auto name(T&& t)
      {
        return std::forward<T>(t).name();
      }
 
    } // NS Tensor
 
    using Tensor::HasTensorTraitsV;
    using Tensor::IsTensor;
    using Tensor::IsTensorOperand;
    using Tensor::TensorTraits;
    using Tensor::tensor;
 
  } // NS ACFem
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_TENSORBASE_HH__