kronecker.hh

#ifndef __DUNE_ACFEM_TENSORS_MODULES_KRONECKER_HH__
#define __DUNE_ACFEM_TENSORS_MODULES_KRONECKER_HH__
 
#include "../../mpl/assign.hh"
#include "../tensorbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace Tensor {
 
      template<class Signature, class PivotIndices, class Field>
      class KroneckerDelta;
 
      template<class Signature, class PivotIndices, class Field>
      struct HasSpecializedExpressionTraits<KroneckerDelta<Signature, PivotIndices, Field> >
        : TrueType
      {};
 
      template<class Field, std::size_t... Dimensions, std::size_t... PivotIndices>
      class KroneckerDelta<Seq<Dimensions...>, Seq<PivotIndices...>, Field>
        : public TensorBase<FloatingPointClosure<Field>, Seq<Dimensions...>, KroneckerDelta<Seq<Dimensions...>, Seq<PivotIndices...>, Field> >
        , public Expressions::SelfExpression<KroneckerDelta<Seq<Dimensions...>, Seq<PivotIndices...>, Field> >
        , public MPL::UniqueTags<ConstantExpression,
                                 ConditionalType<IsTypedValue<Field>::value, TypedValueExpression, void> >
      {
        using BaseType = TensorBase<FloatingPointClosure<Field>, Seq<Dimensions...>, KroneckerDelta<Seq<Dimensions...>, Seq<PivotIndices...>, Field> >;
       public:
        using BaseType::rank;
        using typename BaseType::Signature;
        using typename BaseType::FieldType;
        using PivotSequence = Seq<PivotIndices...>;
        using ValueType = Field;
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        FieldType operator()(Dims... indices) const
        {
          if (std::make_tuple((std::size_t)indices...) != std::make_tuple(PivotIndices...)) {
            return FieldType(0);
          } else if (IsTypedValue<ValueType>::value) {
            return FieldType(ValueType{});
          } else {
            return FieldType(1);
          }
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<sizeof...(Indices) == rank, int> = 0>
        decltype(auto) operator()(Seq<Indices...>) const
        {
          using ReturnType = ConditionalType<isZero<Indices...>(),
                                             IntFraction<0>,
                                             ConditionalType<IsTypedValue<ValueType>::value,
                                                             ValueType,
                                                             IntFraction<1> > >;
          return ReturnType{};
        }
 
       private:
        template<std::size_t... Indices, std::size_t... Pos>
        static bool constexpr isZeroExpander(Seq<Indices...>, Seq<Pos...>)
        {
          return !(... && (Indices == Get<Pos, Seq<PivotIndices...> >::value));
        }
 
       public:
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          // Truncate Pos to size of Indices...
          using RealPos = HeadPart<sizeof...(Indices), Pos>;
 
          return isZeroExpander(Seq<Indices...>{}, RealPos{});
        }
 
        std::string name() const
        {
          return "delta<"+toString(Signature{})+"@["+toString(PivotSequence{})+"]>";
        }
 
        static constexpr auto lookAt()
        {
          return PivotSequence{};
        }
      };
 
      template<class Field, std::size_t D0, std::size_t... RestDims>
      class KroneckerDelta<Seq<D0, RestDims...>, Seq<>, Field>
        : public TensorBase<FloatingPointClosure<Field>, Seq<D0, RestDims...>, KroneckerDelta<Seq<D0, RestDims...>, Seq<>, Field> >
        , public TerminalExpression
        , public ConstantExpression
      {
        using BaseType = TensorBase<FloatingPointClosure<Field>, Seq<D0, RestDims...>, KroneckerDelta<Seq<D0, RestDims...>, Seq<>, Field> >;
       public:
        using BaseType::rank;
        using typename BaseType::Signature;
        using typename BaseType::FieldType;
        using ValueType = Field;
 
       private:
        template<std::size_t... I, class T>
        KroneckerDelta(T&& tuple, Seq<I...>&&)
          : pivotIndices_({{ (std::size_t)std::get<I>(tuple)... }})
        {}
 
       public:
        template<class T,
                 std::enable_if_t<(IsTupleLike<T>::value
                                   && size<T>() == rank
                                   && IsIntegralTuple<T>::value
                                  ), int> = 0>
        KroneckerDelta(T&& tuple)
          : KroneckerDelta(std::forward<T>(tuple), MakeSequenceFor<T>{})
        {}
 
        template<class... Indices,
                 std::enable_if_t<(sizeof...(Indices) != 0
                                   &&
                                   sizeof...(Indices) == rank
                                   &&
                                   IsIntegralPack<Indices...>::value
                                  ), int> = 0>
        KroneckerDelta(Indices... indices)
          : pivotIndices_({{ (std::size_t)indices... }})
        {}
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        FieldType operator()(Dims... indices) const
        {
          if (std::array<std::size_t, rank>({{ (std::size_t)indices... }}) != pivotIndices_) {
            return FieldType(0);
          } else if (IsTypedValue<ValueType>::value) {
            return FieldType(ValueType{});
          } else {
            return FieldType(1);
          }
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<sizeof...(Indices) == rank, int> = 0>
        decltype(auto) operator()(Seq<Indices...>) const
        {
          return (*this)(Indices...);
        }
 
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          return false; // we just do not know ...
        }
 
        template<class... Indices,
                 std::enable_if_t<((sizeof...(Indices) > 0)
                                   && IsIntegralPack<Indices...>::value)
                 , int> = 0>
        void lookAt(Indices... indices)
        {
          static_assert(sizeof...(Indices) == rank,
                        "The number of pivot indices must match the co-dimension.");
          pivotIndices_ = {{ (std::size_t)indices... }};
        }
 
        template<class T,
                 std::enable_if_t<(IsTupleLike<T>::value)
                                  , int> = 0>
        void lookAt(T&& tuple)
        {
          static_assert(size<T>() == rank,
                        "The number of pivot indices must match the co-dimension.");
          ACFem::assign(pivotIndices_, tuple);
        }
 
        const std::array<std::size_t, rank>& lookAt() const&
        {
          return pivotIndices_;
        }
 
        std::string name() const
        {
          return "delta<"+toString(Signature{})+"@["+toString(pivotIndices_)+"]>";
        }
 
       private:
        std::array<std::size_t, rank> pivotIndices_;
      };
 
      template<class T, class SFINAE = void>
      struct IsKroneckerDelta
        : FalseType
      {};
 
      template<class T>
      struct IsKroneckerDelta<T, std::enable_if_t<!IsDecay<T>::value> >
        : IsKroneckerDelta<std::decay_t<T> >
      {};
 
      template<std::size_t... Pos, class PivotIndices, class Field>
      struct IsKroneckerDelta<KroneckerDelta<Seq<Pos...>, PivotIndices, Field> >
        : TrueType
      {};
 
      template<class T>
      struct IsConstKroneckerDelta
        : FalseType
      {};
 
      template<class T>
      struct IsConstKroneckerDelta<T&>
        : IsConstKroneckerDelta<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsConstKroneckerDelta<T&&>
        : IsConstKroneckerDelta<std::decay_t<T> >
      {};
 
      template<std::size_t... Pos, std::size_t... Pivot, class Field>
      struct IsConstKroneckerDelta<KroneckerDelta<Seq<Pos...>, Seq<Pivot...>, Field> >
        : BoolConstant<(sizeof...(Pos) == sizeof...(Pivot))>
      {};
 
      template<class T>
      struct IsDiagonalKroneckerDelta
        : FalseType
      {};
 
      template<class T>
      struct IsDiagonalKroneckerDelta<T&>
        : IsDiagonalKroneckerDelta<T>
      {};
 
      template<class T>
      struct IsDiagonalKroneckerDelta<T&&>
        : IsDiagonalKroneckerDelta<T>
      {};
 
      template<std::size_t... Pos, std::size_t Pivot0, std::size_t... Pivot, class Field>
      struct IsDiagonalKroneckerDelta<KroneckerDelta<Seq<Pos...>, Seq<Pivot0, Pivot...>, Field> >
        : BoolConstant<(... && (Pivot0 == Pivot))>
      {};
 
      template<class T>
      using PivotSequence = typename std::decay_t<T>::PivotSequence;
 
      template<std::size_t... Index, class Dims, class F = Expressions::Closure,
               std::enable_if_t<IsSequence<Dims>::value, int> = 0>
      auto kroneckerDelta(Dims = Dims{}, Seq<Index...> = Seq<Index...>{}, F closure = F{})
      {
        static_assert(Dims::size() == sizeof...(Index),
                      "Kronecker-Delta needs an index tuple of the size of the rank of the tensor.");
        static_assert(Seq<Index...>{} < Dims{},
                      "Index tuple not in range of tensor dimensions");
 
        return closure(KroneckerDelta<Dims, Seq<Index...>, IntFraction<1> >{});
      }
 
      template<std::size_t... Dimensions, std::size_t... Index, class F = Expressions::Closure>
      auto kroneckerDelta(Seq<Dimensions...>, Seq<Index...>, F closure = F{})
      {
        static_assert(sizeof...(Dimensions) == sizeof...(Index),
                      "Kronecker-Delta needs an index tuple of the size of the rank of the tensor.");
        static_assert(Seq<Index...>{} < Seq<Dimensions...>{},
                      "Index tuple not in range of tensor dimensions");
 
        return closure(KroneckerDelta<Seq<Dimensions...>, Seq<Index...>, IntFraction<1> >{});
      }
 
      template<class T, std::size_t... Index, class F = Expressions::Closure,
               std::enable_if_t<IsTensorOperand<T>::value, int> = 0>
      auto kroneckerDelta(T&&, Seq<Index...>, F = F{})
      {
        return kroneckerDelta(typename TensorTraits<T>::Signature{}, Seq<Index...>{}, F{});
      }
 
      template<std::size_t... Dimensions, class Tuple, class F = Expressions::Closure,
               std::enable_if_t<IsTupleLike<Tuple>::value, int> = 0>
      auto kroneckerDelta(Seq<Dimensions...>, const Tuple& tuple, F closure = F{})
      {
        static_assert(sizeof...(Dimensions) == size<Tuple>(),
                      "Kronecker-Delta needs an index tuple of the size of the rank of the tensor.");
        return closure(KroneckerDelta<Seq<Dimensions...>, Seq<>, IntFraction<1> >(tuple));
      }
 
      template<std::size_t... Dimensions, class... Indices,
               std::enable_if_t<IsIntegralPack<Indices...>::value, int> = 0>
      auto kroneckerDelta(Seq<Dimensions...>, Indices... indices)
      {
        static_assert(sizeof...(Dimensions) == sizeof...(Indices),
                      "Kronecker-Delta needs an index tuple of the size of the rank of the tensor.");
        return expressionClosure(KroneckerDelta<Seq<Dimensions...>, Seq<>, IntFraction<1> >(indices...));
      }
 
 
 
    } // NS Tensor
 
    using Tensor::kroneckerDelta;
 
    template<class Field, class Signature, class Indices>
    struct ExpressionTraits<Tensor::KroneckerDelta<Signature, Indices, Field> >
      : ExpressionTraits<Field>
    {
     private:
      using BaseType = ExpressionTraits<Field>;
     public:
      using ExpressionType = Tensor::KroneckerDelta<Signature, Indices, Field>;
      static constexpr bool isOne = BaseType::isOne && ExpressionType::rank == 0;
      static constexpr bool isMinusOne = BaseType::isMinusOne && ExpressionType::rank == 0;
      static constexpr bool isNonZero = BaseType::isNonZero && ExpressionType::rank == 0;
      static constexpr bool isPositive = BaseType::isPositive && ExpressionType::rank == 0;
      static constexpr bool isNegative = BaseType::isNegative && ExpressionType::rank == 0;
      static constexpr bool isNonSingular = ExpressionType::rank == 0;
      using Sign = ExpressionSign<isNonZero, BaseType::isSemiPositive, BaseType::isSemiNegative>;
 
      static constexpr bool isIndependent = BaseType::isIndependent;
      static constexpr bool isConstant = BaseType::isConstant && !std::is_same<Indices, Tensor::Seq<> >::value;
      static constexpr bool isTypedValue = BaseType::isTypedValue && !std::is_same<Indices, Tensor::Seq<> >::value;
    };
 
    namespace Expressions {
 
      template<class Field, class Signature, std::size_t Pivot0, std::size_t... Pivot>
      constexpr inline std::size_t WeightV<Tensor::KroneckerDelta<Signature, Tensor::Seq<Pivot0, Pivot...>, Field> > = (Pivot0 + ... + Pivot) + 1000*(... && (Pivot0 != Pivot));
 
    }
 
  } // NS ACFem
 
  template<class Field, class Signature, class Indices>
  struct FieldTraits<ACFem::Tensor::KroneckerDelta<Signature, Indices, Field> >
    : FieldTraits<std::decay_t<Field> >
  {};
 
} // NS Dune
 
#endif // __DUNE_VECTOR_EXPRESSIONS_KRONECKER_HH__