blockeye.hh

#ifndef __DUNE_ACFEM_TENSORS_MODULES_BLOCKEYE_HH__
#define __DUNE_ACFEM_TENSORS_MODULES_BLOCKEYE_HH__
 
#include "../tensorbase.hh"
#include "constant.hh"
#include "eye.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace Tensor {
 
      template<std::size_t BlockRank, class BlockSignature, class Field = IntFraction<1> >
      class BlockEye;
 
      template<std::size_t BlockRank, class BlockSignature, class Field>
      struct HasSpecializedExpressionTraits<BlockEye<BlockRank, BlockSignature, Field> >
        : TrueType
      {};
 
      template<class T>
      struct IsBlockEye
        : FalseType
      {};
 
      template<class T>
      struct IsBlockEye<T&>
        : IsBlockEye<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsBlockEye<T&&>
        : IsBlockEye<std::decay_t<T> >
      {};
 
      template<std::size_t BlockRank, class BlockSignature, class Field>
      struct IsBlockEye<BlockEye<BlockRank, BlockSignature, Field> >
        : TrueType
      {};
 
      template<class T, class SFINAE = void>
      struct BlockEyeTraits
      {
        using ValueType = void;
        static constexpr std::size_t blockRank_ = 0;
        using BlockSignature = Seq<>;
        static constexpr std::size_t blockSize_ = 0;
      };
 
      template<class T>
      struct BlockEyeTraits<T, std::enable_if_t<!IsDecay<T>::value> >
        : BlockEyeTraits<std::decay_t<T> >
      {};
 
      template<std::size_t Rank, class Signature, class Field>
      struct BlockEyeTraits<BlockEye<Rank, Signature, Field> >
      {
        using ValueType = Field;
        static constexpr std::size_t blockRank_ = Rank;
        using BlockSignature = Signature;
        static constexpr std::size_t blockSize_ = BlockSignature::size();
      };
 
      template<class Pos, std::size_t BlockRank, class Dims, class SFINAE = void>
      constexpr inline bool IsBlockEyePositionsV = false;
 
      template<std::size_t... P, std::size_t BlockRank, class Dims>
      constexpr inline bool IsBlockEyePositionsV<Seq<P...>, BlockRank, Dims, std::enable_if_t<(sizeof...(P) % Dims::size() == 0)> > = std::is_same<
        typename SortSequence<Seq<(P % Dims::size())...> >::Result,
        MakeIndexSequence<Dims::size(), 0, 1, sizeof...(P) / Dims::size()> >::value;
 
      template<std::size_t BlockRank, std::size_t... Dimensions, class Field>
      class BlockEye<BlockRank, Seq<Dimensions...>, Field>
        : public TensorBase<FloatingPointClosure<Field>, SequenceProd<BlockRank, Seq<Dimensions...> >,
                            BlockEye<BlockRank, Seq<Dimensions...>, Field> >
        , public Expressions::SelfExpression<BlockEye<BlockRank, Seq<Dimensions...>, Field > >
        , public MPL::UniqueTags<ConstantExpression,
                                 ConditionalType<IsTypedValue<Field>::value, TypedValueExpression, void>
                                 >
      {
        using BaseType = TensorBase<FloatingPointClosure<Field>, SequenceProd<BlockRank, Seq<Dimensions...> >,
                                    BlockEye<BlockRank, Seq<Dimensions...>, Field> >;
        static constexpr std::size_t blockSize_ = sizeof...(Dimensions);
       public:
        using BaseType::rank;
        using typename BaseType::FieldType;
        using typename BaseType::Signature;
        using ValueType = Field;
        using BlockSignature = Seq<Dimensions...>;
        static constexpr std::size_t blockRank_ = BlockRank;
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        FieldType operator()(Dims... indices) const
        {
          bool isConstant = true;
#if __GNUC__ <= 7 && !defined(__clang__)
          const auto indexArray = std::array<std::size_t, sizeof...(Dims)>({{ (std::size_t)indices... }});
          for (unsigned i = 0; i < blockSize_; ++i) {
            const auto first = indexArray[i];
            for (unsigned j = 1; j < BlockRank; ++j) {
              const auto next = indexArray[i+j*blockSize_];
              if (next != first) {
                return FieldType(0);
              }
            }
          }
#else
          // GCC does not like the following and catches an internal compiler error:
          forLoop<blockSize_>([&](auto i) {
              using I = decltype(i);
              std::size_t first = get<I::value>(std::forward_as_tuple(indices...));
              forEach(MakeIndexSequence<BlockRank-1, 1>{}, [&](auto j) {
                  using J = decltype(j);
                  std::size_t next = std::get<I::value+J::value*blockSize_>(std::forward_as_tuple(indices...));
                  if (next != first) {
                    isConstant = false;
                  }
                });
            });
#endif
          if (!isConstant) {
            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<class Indices, class Pos>
        struct IsZeroFunctor
        {
          // Get sequence of positions equal to index nr. N.*/
          template<std::size_t N>
          using PositionsOf = TransformSequence<
              Pos, Seq<>, IndexFunctor, AcceptModuloFunctor<std::size_t, blockSize_, N> >;
 
          // The following must yield a constant sequence.
          template<std::size_t N>
          using IndexBlock = SequenceSlice<Indices, PositionsOf<N> >;
 
          template<std::size_t N>
          using ZeroMatch = BoolConstant<!isConstant(IndexBlock<N>{})>;
        };
 
        template<class Indices, class Pos, std::size_t... N>
        static bool constexpr isZeroExpander(Indices, Pos, Seq<N...>)
        {
          return (... || (IsZeroFunctor<Indices, Pos>::template ZeroMatch<N>::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{}, MakeSequenceFor<RealPos>{});
        }
 
        std::string name() const
        {
          return "blockEye<"+std::to_string(blockRank_)+"|"+toString(BlockSignature{})+">";
        }
      };
 
      // blockEye<2>(TensorTraits<T>::signature())
 
      template<std::size_t Rank, std::size_t... Dimensions, class F = Expressions::Closure>
      auto blockEye(Seq<Dimensions...> = Seq<Dimensions...>{}, F closure = F{})
      {
        if constexpr (Rank > 0) {
          if constexpr (sizeof...(Dimensions) == 1) {
            return eye(MakeConstantSequence<Rank, Head<Seq<Dimensions...> >::value>{}, F{});
          } else {
            return closure(BlockEye<Rank, Seq<Dimensions...>, IntFraction<1> >{});
          }
        } else {
          return constantTensor(1_f, Seq<>{}, F{});
        }
      }
 
      template<std::size_t Rank, class T, class F = Expressions::Closure,
               std::enable_if_t<!IsSequence<T>::value, int> = 0>
      auto blockEye(T&&, F = F{})
      {
        return blockEye<Rank>(typename TensorTraits<T>::Signature{}, F{});
      }
 
 
 
    } // NS Tensor
 
    template<class Field, std::size_t BlockRank, class BlockSignature>
    struct ExpressionTraits<Tensor::BlockEye<BlockRank, BlockSignature, Field> >
      : ExpressionTraits<Field>
    {
     private:
      using BaseType = ExpressionTraits<Field>;
     public:
      using ExpressionType = Tensor::BlockEye<BlockRank, BlockSignature, 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;
      using Sign = ExpressionSign<isNonZero, BaseType::isSemiPositive, BaseType::isSemiNegative>;
    };
 
  } // NS ACFem
 
  template<class Field, std::size_t N, class Signature>
  struct FieldTraits<ACFem::Tensor::BlockEye<N, Signature, Field> >
    : FieldTraits<std::decay_t<Field> >
  {};
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_MODULES_BLOCKEYE_HH__