restriction.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONS_RESTRICTION_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONS_RESTRICTION_HH__
 
#include "../../expressions/storage.hh"
#include "../../expressions/operationtraits.hh"
#include "../../mpl/insertat.hh"
#include "../../mpl/toarray.hh"
#include "../tensorbase.hh"
#include "../expressiontraits.hh"
#include "restrictiondetail.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // Forward.
    template<class Dims, class Pivot>
    struct RestrictionOperation;
 
    namespace Tensor {
 
      template<class Tensor, class Dims, class Indices = Seq<> >
      class Restriction;
 
      template<class Tensor, std::size_t... IndexPositions, std::size_t... PivotIndices>
      class Restriction<Tensor, Seq<IndexPositions...>, Seq<PivotIndices...> >
        : public TensorBase<typename TensorTraits<Tensor>::FieldType,
                            ACFem::SubSequenceComplement<typename TensorTraits<Tensor>::Signature, IndexPositions...>,
                            Restriction<Tensor, Seq<IndexPositions...>, Seq<PivotIndices...> > >
        , public Expressions::Storage<OperationTraits<RestrictionOperation<Seq<IndexPositions...>, Seq<PivotIndices...> > >, Tensor>
      {
        static_assert(isSorted(Seq<IndexPositions...>{}),
                      "Index positions must be sorted.");
 
       public:
        using DefectPositions = Seq<IndexPositions...>;
        using PivotSequence = Seq<PivotIndices...>;
        using FunctorType = OperationTraits<RestrictionOperation<DefectPositions, PivotSequence> >;
       private:
        using ArgType = Tensor;
        using HostType = std::decay_t<Tensor>;
        using BaseType = TensorBase<typename HostType::FieldType,
                                    ACFem::SubSequenceComplement<typename HostType::Signature, IndexPositions...>,
                                    Restriction<Tensor, DefectPositions, PivotSequence> >;
        using StorageType = Expressions::Storage<FunctorType, Tensor>;
       public:
        using DefectSignature = SequenceSlice<typename HostType::Signature, DefectPositions>;
        static constexpr std::size_t defectRank_ = DefectSignature::size();
       public:
        using BaseType::rank;
        using typename BaseType::Signature;
        using typename BaseType::FieldType;
        using StorageType::operand;
        using StorageType::operation;
        using BaseType::operator[];
 
        template<class Arg,
                 std::enable_if_t<std::is_constructible<ArgType, Arg>::value, int> = 0>
        Restriction(Arg&& host, PivotSequence = PivotSequence{})
          : StorageType(std::forward<Arg>(host))
        {}
 
        template<
          class... Dummy,
          std::enable_if_t<(sizeof...(Dummy) == 0
                            && IsTypedValue<ArgType>::value), int> = 0>
        Restriction(Dummy&&...)
          : StorageType(ArgType{})
        {}
 
        auto operator[](std::size_t i)
        {
          using Positions = JoinedDefects<DefectPositions, Seq<0> >;
          using Injections = JoinedInjections<DefectPositions, Seq<0> >;
          using ResultType = Restriction<Tensor, Positions, Seq<> >;
          return ResultType(operand(0_c), insertAt(std::make_tuple(PivotIndices...), i, Head<Injections>::value));
        }
 
        auto operator[](std::size_t i) const
        {
          using Positions = JoinedDefects<DefectPositions, Seq<0> >;
          using Injections = JoinedInjections<DefectPositions, Seq<0> >;
          using ResultType = Restriction<Tensor, Positions, Seq<> >;
          return ResultType(operand(0_c), insertAt(std::make_tuple(PivotIndices...), i, Head<Injections>::value));
        }
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices)
        {
          return tensorValue(operand(0_c), insertAt<IndexPositions...>(std::forward_as_tuple(indices...), toArray(PivotSequence{})));
        }
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices) const
        {
          return tensorValue(operand(0_c), insertAt<IndexPositions...>(std::forward_as_tuple(indices...), toArray(PivotSequence{})));
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
          ), int> = 0>
        decltype(auto) constexpr operator()(Seq<Indices...>)
        {
          return operand(0_c)(InsertAt<Seq<Indices...>, PivotSequence, DefectPositions>{});
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
          ), int> = 0>
        decltype(auto) constexpr operator()(Seq<Indices...>) const
        {
          return operand(0_c)(InsertAt<Seq<Indices...>, PivotSequence, DefectPositions>{});
        }
 
        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 and sort it, otherwise
          // the JoinedDefects stuff will not work properly
          using TruncatedPos = HeadPart<sizeof...(Indices), Pos>;
          using RealPos = typename SortSequence<TruncatedPos>::Result;
          using Permutation = typename SortSequence<TruncatedPos>::Permutation;
          using RealIndices = PermuteSequence<Seq<Indices...>, Permutation>;
 
          using Positions = JoinedDefects<DefectPositions, RealPos>;
          using Injections = JoinedInjections<DefectPositions, RealPos>;
 
          return HostType::isZero(InsertAt<PivotSequence, RealIndices, Injections>{}, Positions{});
        }
 
        static constexpr auto lookAt()
        {
          return PivotSequence{};
        }
 
        std::string name() const
        {
          using T = Tensor;
          std::string pfx = std::is_reference<T>::value ? (RefersConst<T>::value ? "cref" : "ref") : "";
          return operationName(operation(), pfx+operand(0_c).name());
        }
 
      };
 
      template<class Tensor, std::size_t... IndexPositions>
      class Restriction<Tensor, Seq<IndexPositions...>, Seq<> >
        : public TensorBase<typename TensorTraits<Tensor>::FieldType,
                            ACFem::SubSequenceComplement<typename TensorTraits<Tensor>::Signature, IndexPositions...>,
                            Restriction<Tensor, Seq<IndexPositions...>, Seq<> > >
        , public Expressions::Storage<OperationTraits<RestrictionOperation<Seq<IndexPositions...>, Seq<> > >, Tensor>
      {
       public:
        using DefectPositions = Seq<IndexPositions...>;
        using PivotSequence = Seq<>;
        using FunctorType = OperationTraits<RestrictionOperation<Seq<IndexPositions...>, Seq<> > >;
       private:
        using ArgType = Tensor;
        using HostType = std::decay_t<Tensor>;
        using BaseType = TensorBase<typename HostType::FieldType,
                                    ACFem::SubSequenceComplement<typename HostType::Signature, IndexPositions...>,
                                    Restriction<Tensor, Seq<IndexPositions...>, Seq<> > >;
        using StorageType = Expressions::Storage<FunctorType, Tensor>;
        using StorageType::functor_;
       public:
        using StorageType::operation;
        using StorageType::operand;
        using DefectSignature = SequenceSlice<typename HostType::Signature, DefectPositions>;
        static constexpr std::size_t defectRank_ = DefectSignature::size();
       public:
        using BaseType::rank;
        using typename BaseType::Signature;
        using typename BaseType::FieldType;
        using PivotIndexType = typename FunctorType::IndexType;
 
       private:
        template<class Arg, class T, std::size_t... Idx,
                 std::enable_if_t<std::is_constructible<ArgType, Arg>::value, int> = 0>
        Restriction(Arg&& host, T&& t, IndexSequence<Idx...>)
          : StorageType(std::forward<Arg>(host), FunctorType({{ get<Idx>(std::forward<T>(t))... }}))
        {}
 
       public:
        template<class Arg, class T,
                 std::enable_if_t<(std::is_constructible<ArgType, Arg>::value
                                   && IsTupleLike<T>::value
                                   && size<T>() == sizeof...(IndexPositions)
                                   && IsIntegralTuple<T>::value)
                                , int> = 0>
        Restriction(Arg&& host, const T& tuple)
          : StorageType(FunctorType(tuple), std::forward<Arg>(host))
        {}
 
        template<class Arg, class T, std::size_t N,
                 std::enable_if_t<(std::is_constructible<ArgType, Arg>::value
                                   && IsTupleLike<T>::value
                                   && (N == sizeof...(IndexPositions) || N == 0)
                                   && std::is_integral<T>::value)
                                , int> = 0>
        Restriction(Arg&& host, const T (&t)[N])
          : StorageType(std::forward<Arg>(host), FunctorType(t))
        {}
 
        template<class Arg, class T,
                 std::enable_if_t<(std::is_constructible<ArgType, Arg>::value
                                   && IsTupleLike<T>::value
                                   && size<T>() == 0
                                   && IsIntegralTuple<T>::value)
                                , int> = 0>
        Restriction(Arg&& host, T&& tuple)
          : StorageType(std::forward<Arg>(host))
        {}
 
        template<class Arg, class... Indices,
                 std::enable_if_t<(std::is_constructible<ArgType, Arg>::value
                                   && sizeof...(Indices) != 0
                                   && sizeof...(Indices) == sizeof...(IndexPositions)
                                   && IsIntegralPack<Indices...>::value
          ), int> = 0>
        Restriction(Arg&& host, Indices... indices)
          : StorageType(std::forward<Arg>(host), FunctorType({{{ (PivotIndexType)indices... }}}))
        {}
 
        template<class Arg, std::size_t... Pivots,
                 std::enable_if_t<(std::is_constructible<ArgType, Arg>::value
                                   && sizeof...(Pivots) != 0
          ), int> = 0>
        Restriction(Arg&& host, IndexSequence<Pivots...>)
          : Restriction(std::forward<Arg>(host), Pivots...)
        {}
 
        template<class Arg, std::enable_if_t<std::is_constructible<ArgType, Arg>::value, int> = 0>
        Restriction(ArgType&& host, Seq<> = Seq<>{})
          : StorageType(std::forward<Arg>(host))
        {}
 
        auto operator[](std::size_t i)
        {
          using Positions = JoinedDefects<DefectPositions, Seq<0> >;
          using Injections = JoinedInjections<DefectPositions, Seq<0> >;
 
          return Restriction<Tensor, Positions, Seq<> >(operand(0_c), insertAt(lookAt(), i, Head<Injections>{}));
        }
 
        auto operator[](std::size_t i) const
        {
          using Positions = JoinedDefects<DefectPositions, Seq<0> >;
          using Injections = JoinedInjections<DefectPositions, Seq<0> >;
 
          return Restriction<Tensor, Positions, Seq<> >(operand(0_c), insertAt(lookAt(), i, Head<Injections>{}));
        }
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices)
        {
          return tensorValue(operand(0_c), insertAt<IndexPositions...>(std::forward_as_tuple(indices...), lookAt()));
        }
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices) const
        {
          return tensorValue(operand(0_c), insertAt<IndexPositions...>(std::forward_as_tuple(indices...), lookAt()));
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
          ), int> = 0>
        decltype(auto) constexpr operator()(Seq<Indices...>)
        {
          return (*this)(Indices...);
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
          ), int> = 0>
        decltype(auto) constexpr operator()(Seq<Indices...>) const
        {
          return (*this)(Indices...);
        }
 
       private:
        template<class Indices, class Positions, class Injections, std::size_t... N>
        static bool constexpr isZeroExpander(Indices, Positions, Injections, Seq<N...>)
        {
          return (... && HostType::isZero(InsertAt<MultiIndex<N, DefectSignature>, Indices, Injections>{}, Positions{}));
        }
 
       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 and sort it, otherwise
          // the JoinedDefects stuff will not work properly
          using TruncatedPos = HeadPart<sizeof...(Indices), Pos>;
          using RealPos = typename SortSequence<TruncatedPos>::Result;
          using Permutation = typename SortSequence<TruncatedPos>::Permutation;
          using RealIndices = PermuteSequence<Seq<Indices...>, Permutation>;
 
          using Positions = JoinedDefects<DefectPositions, RealPos>;
          using Injections = JoinedInjections<DefectPositions, RealPos>;
 
          //std::clog << Seq<1,2,3>{} << std::endl;
          //std::clog << Injections{} << std::endl;
 
          return isZeroExpander(RealIndices{}, Positions{}, Injections{},
                                MakeIndexSequence<multiDim(DefectSignature{})>{});
        }
 
        template<class... Indices,
                 std::enable_if_t<((sizeof...(Indices) > 0)
                                   && IsIntegralPack<Indices...>::value)
                 , int> = 0>
        void lookAt(Indices... indices)
        {
          static_assert(sizeof...(Indices) <= defectRank_,
                        "The number of pivot indices must match the co-dimension.");
          auto indexList = { (std::size_t)indices... };
          std::copy(indexList.begin(), indexList.end(), functor_.pivotIndices_.end() - indexList.size());
        }
 
        template<class T,
                 std::enable_if_t<(IsTupleLike<T>::value)
                                  , int> = 0>
        void lookAt(T&& tuple)
        {
          static_assert(size<T>() <= defectRank_,
                        "The number of pivot indices must match the co-dimension.");
          std::copy(tuple.begin(), tuple.end(), functor_.pivotIndices_.end() - tuple.size());
        }
 
        const auto& lookAt() const
        {
          return functor_.pivotIndices_;
        }
 
        auto& lookAt()
        {
          return functor_.pivotIndices_;
        }
 
        std::string name() const
        {
          using T = Tensor;
          std::string pfx = std::is_reference<T>::value ? (RefersConst<T>::value ? "cref" : "ref") : "";
          return operationName(operation(), pfx+operand(0_c).name());
        }
 
      };
 
      template<class T>
      struct IsRestriction
        : FalseType
      {};
 
      template<class T>
      struct IsRestriction<T&>
        : IsRestriction<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsRestriction<T&&>
        : IsRestriction<std::decay_t<T> >
      {};
 
      template<class T, class DefectPositions, class PivotIndices>
      struct IsRestriction<Restriction<T, DefectPositions, PivotIndices> >
        : TrueType
      {};
 
      template<class T, class SFINAE = void>
      struct IsConstRestriction
        : FalseType
      {};
 
      template<class T>
      struct IsConstRestriction<T, std::enable_if_t<!IsDecay<T>::value> >
        : IsConstRestriction<std::decay_t<T> >
      {};
 
      template<class T, class Pos, class Pivots>
      struct IsConstRestriction<Restriction<T, Pos, Pivots>,
                                std::enable_if_t<Pos::size() == Pivots::size()> >
        : TrueType
      {};
 
      template<class T>
      struct IsDynamicRestriction
        : FalseType
      {};
 
      template<class T>
      struct IsDynamicRestriction<T&>
        : IsDynamicRestriction<std::decay_t<T> >
      {};
 
      template<class T>
      struct IsDynamicRestriction<T&&>
        : IsDynamicRestriction<std::decay_t<T> >
      {};
 
      template<class T, std::size_t P0, std::size_t... PRest>
      struct IsDynamicRestriction<Restriction<T, Seq<P0, PRest...>, Seq<> > >
        : TrueType
      {};
 
      template<class Pos, class Ind, class T>
      constexpr auto operate(Expressions::DontOptimize, const OperationTraits<RestrictionOperation<Pos, Ind> >& f, T&& t)
      {
        DUNE_ACFEM_RECORD_OPTIMIZATION;
 
        return Restriction<T, Pos, Ind>(std::forward<T>(t), f.pivotIndices_);
      }
 
      using Expressions::finalize;
 
      template<std::size_t... IndexPositions, class T, std::size_t... Indices,
               std::enable_if_t<(IsProperTensor<T>::value
                                 && sizeof...(IndexPositions) == sizeof...(Indices)
        ), int> = 0>
      auto restriction(T&& t, Seq<Indices...>, Seq<IndexPositions...> = Seq<IndexPositions...>{})
      {
        using Operation = RestrictionOperation<Seq<IndexPositions...>, Seq<Indices...> >;
        return finalize<Operation>(std::forward<T>(t));
      }
 
      template<std::size_t... IndexPositions, class T, class Indices,
               std::enable_if_t<(IsTupleLike<Indices>::value
                                 && IsProperTensor<T>::value
        ), int> = 0>
      auto restriction(T&& t, Indices&& indices, Seq<IndexPositions...> = Seq<IndexPositions...>{})
      {
        static_assert(sizeof...(IndexPositions) == size<Indices>() || size<Indices>() == 0,
                      "The number of pivot indices must match the number of defect positions.");
 
        using Functor = OperationTraits<RestrictionOperation<Seq<IndexPositions...>, Seq<> > >;
 
        return finalize(Functor(indices), std::forward<T>(t));
      }
 
      template<std::size_t... IndexPositions, class T, class I, std::size_t N, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto restriction(T&& t, const I (&l)[N], Seq<IndexPositions...> = Seq<IndexPositions...>{})
      {
        static_assert(sizeof...(IndexPositions) == N || N == 0,
                      "Either all or no pivot indices must be given.");
 
        using Functor = OperationTraits<RestrictionOperation<Seq<IndexPositions...>, Seq<> > >;
 
        return finalize(Functor(l), std::forward<T>(t));
      }
 
      template<std::size_t... IndexPositions, class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto restriction(T&& t)
      {
        using Operation = RestrictionOperation<Seq<IndexPositions...>, Seq<> >;
        return finalize<Operation>(std::forward<T>(t));
      }
 
      template<class T, class DefectPositions, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto restriction(T&& t, DefectPositions)
      {
        using Operation = RestrictionOperation<DefectPositions, Seq<> >;
        return finalize<Operation>(std::forward<T>(t));
      }
 
    } // NS Tensor
 
    using Tensor::restriction;
 
    namespace Expressions {
 
      template<class T>
      constexpr inline std::size_t WeightV<T, std::enable_if_t<Tensor::IsConstRestriction<T>::value> > =
        ExpressionWeightV<Operand<0, T> > + prod(typename std::decay_t<T>::DefectPositions{}) + prod(std::decay_t<T>::lookAt());
 
    }
 
  } // NS ACFem
 
  template<class Tensor, class Pos, class Indices>
  struct FieldTraits<ACFem::Tensor::Restriction<Tensor, Pos, Indices> >
    : FieldTraits<Tensor>
  {};
 
 
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONS_RESTRICTION_HH__