expressions.hh

#ifndef __DUNE_ACFEM_TENSORS_EXPRESSIONS_HH__
#define __DUNE_ACFEM_TENSORS_EXPRESSIONS_HH__
 
#include "../common/literals.hh"
#include "../expressions/storage.hh"
#include "../expressions/expressionoperations.hh"
#include "../expressions/constantoperations.hh"
#include "../expressions/operandpromotion.hh"
 
#include "tensorbase.hh"
#include "expressiontraits.hh"
#include "bindings.hh"
//#include "optimization.hh"
#include "modules.hh"
#include "operations/einsum.hh"
#include "operations/reciprocal.hh"
#include "operationtraits.hh"
 
#include "binaryexpression.hh"
#include "unaryexpression.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace Tensor {
 
// Suck in using declarations. This is needed as we define operator
// functions of the same name in our own namespace.
#include "../expressions/usingstd.hh"
#include "../expressions/usingtyped.hh"
#include "../expressions/usingpromotion.hh"
#include "../mpl/usingcomparisons.hh"
 
      template<class F, class T, std::enable_if_t<IsTensor<T>::value, int> = 0>
      auto zero(F&&, T&&)
      {
        return zeros(TensorTraits<T>::signature(), Disclosure{});
      }
 
      template<class T0, class T1, std::enable_if_t<IsTensor<T0>::value && IsTensor<T1>::value, int> = 0>
      auto zero(MinusFunctor, T0&&, T1&&)
      {
        return zeros(TensorTraits<T0>::signature(), Disclosure{});
      }
 
      template<class Seq0, class Seq1, class Dims, class T0, class T1>
      auto zero(OperationTraits<EinsumOperation<Seq0, Seq1, Dims> >, T0&&, T1&&)
      {
        return zeros(EinsumSignature<T0, Seq0, T1, Seq1>{}, Disclosure{});
      }
 
      template<class Seq0, class Seq1, class Dims, class T0, class T1>
      auto zero(OperationTraits<TensorProductOperation<Seq0, Seq1, Dims> >, T0&&, T1&&)
      {
        return zeros(ProductSignature<T0, Seq0, T1, Seq1>{}, Disclosure{});
      }
 
      template<class T, std::enable_if_t<Tensor::IsTensor<T>::value, int> = 0>
      auto one(T&&)
      {
        return Tensor::ConstantTensor<IntFraction<1>, Seq<> >{};
      }
 
      template<class T0, class T1, std::enable_if_t<IsTensor<T0>::value && IsTensor<T1>::value, int> = 0>
      constexpr auto multiplyScalars(T0&& t0, T1&& t1)
      {
        return operate(ScalarEinsumFunctor{}, std::forward<T0>(t0), std::forward<T1>(t1));
      }
 
 
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr decltype(auto) operator+(T1&& t1, T2&& t2)
      {
        return finalize<PlusOperation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr decltype(auto) operator-(T1&& t1, T2&& t2)
      {
        return finalize<MinusOperation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      constexpr decltype(auto) operator-(T&& t)
      {
        return finalize<MinusOperation>(std::forward<T>(t));
      }
 
      template<std::size_t N, class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr auto contractInner(T1&& t1, T2&& t2, IndexConstant<N> = IndexConstant<N>{})
      {
        static_assert(TailPart<N, typename TensorTraits<T1>::Signature>{}
                      ==
                      HeadPart<N, typename TensorTraits<T2>::Signature>{},
                      "Tensor signatures have to match for inner product.");
        using LeftPos = MakeIndexSequence<N, TensorTraits<T1>::rank - N>;
        using RightPos = MakeIndexSequence<N>;
        return einsum<LeftPos, RightPos>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto contractInner(T1&& t1, T2&& t2)
      {
        return contractInner<1>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto outer(T1&& t1, T2&& t2)
      {
        return einsum<Seq<>, Seq<> >(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      constexpr auto trace(T&& t)
      {
        using IndexPositions = MakeIndexSequence<TensorTraits<T>::rank>;
 
        return einsum<IndexPositions, IndexPositions>(std::forward<T>(t), eye(TensorTraits<T>::signature()));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto inner(T1&& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::signature() == TensorTraits<T2>::signature(),
                      "Signature of both operands have to match for inner product.");
 
        using IndexPositions = MakeIndexSequence<TensorTraits<T1>::rank>;
 
        return einsum<IndexPositions, IndexPositions>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class Arg, std::enable_if_t<IsProperTensor<Arg>::value, int> = 0>
      constexpr auto frobenius2(Arg&& arg)
      {
        using Seq = MakeIndexSequence<TensorTraits<Arg>::rank>;
        return einsum<Seq, Seq>(std::forward<Arg>(arg), std::forward<Arg>(arg));
      }
 
      template<class T1, class T2,
               std::enable_if_t<(AreProperTensors<T1, T2>::value
                                 && TensorTraits<T1>::rank > 0
                                 && TensorTraits<T2>::rank > 0
                                 && TensorTraits<T1>::template dim<TensorTraits<T1>::rank - 1>() == TensorTraits<T2>::template dim<0>()
                 ), int> = 0>
      auto operator*(T1&& t1, T2&& t2)
      {
        return contractInner<1>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2,
               std::enable_if_t<(AreProperTensors<T1, T2>::value
                                 && TensorTraits<T1>::rank * TensorTraits<T2>::rank == 0
        ), int> = 0>
      auto operator*(T1&& t1, T2&& t2)
      {
        return einsum<Seq<>, Seq<> >(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto operator&(T1&& t1, T2&& t2)
      {
        return einsum(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2,
               std::enable_if_t<(AreProperTensors<T1, T2>::value
                                 && TensorTraits<T2>::rank == 0
                 ), int> = 0>
      auto operator/(T1&& t1, T2&& t2)
      {
        return std::forward<T1>(t1) * reciprocal(std::forward<T2>(t2));
      }
 
      template<class T1, class T2,
               std::enable_if_t<(IsTensorOperand<T1>::value
                                 && IsTensorOperand<T2>::value
                                 && (IsTensor<T1>::value || IsTensor<T2>::value)
        ), int> = 0>
      T1& operator+=(T1& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::hasMutableComponents,
                      "Destination of assignment must have assignable components.");
 
        decltype(auto) t1_ = tensor(std::forward<T1&>(t1));
        assign(t1_, tensor(std::forward<T2>(t2)),
               [](auto& dst, auto&& src) {
                 dst += src;
               });
        return t1;
      }
 
      template<class T1, class T2,
               std::enable_if_t<(IsTensorOperand<T1>::value
                                 && IsTensorOperand<T2>::value
                                 && (IsTensor<T1>::value || IsTensor<T2>::value)
        ), int> = 0>
      T1& operator-=(T1& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::hasMutableComponents,
                      "Destination of assignment must have assignable components.");
 
        decltype(auto) t1_ = tensor(std::forward<T1&>(t1));
        assign(t1_, tensor(std::forward<T2>(t2)),
               [](auto& dst, auto&& src) {
                 dst -= src;
               });
        return t1;
      }
 
      template<class T1, class T2,
               std::enable_if_t<(IsTensorOperand<T1>::value
                                 && IsTensorOperand<T2>::value
                                 && (IsTensor<T1>::value || IsTensor<T2>::value)
                                 && TensorTraits<T2>::rank == 0
                 ), int> = 0>
      T1& operator*=(T1& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::hasMutableComponents,
                      "Destination of assignment must have assignable components.");
 
        decltype(auto) t1_ = tensor(std::forward<T1&>(t1));
        assign(t1_, tensor(std::forward<T2>(t2)),
               [](auto& dst, auto&& src) {
                 dst *= src;
               });
        return t1;
      }
 
      template<class T1, class T2, class SFINAE = void>
      struct AreRuntimeComparable
        : FalseType
      {};
 
      template<class T1, class T2>
      struct AreRuntimeComparable<
        T1, T2, std::enable_if_t<(AreProperTensors<T1, T2>::value
                                  && !AreRuntimeEqualOperands<T1, T2>::value
                                  && std::is_same<typename std::decay_t<T1>::Signature, typename std::decay_t<T2>::Signature>::value
        )> >
        : TrueType
      {};
 
      template<class T1, class T2,
               std::enable_if_t<(AreProperTensors<T1, T2>::value
                                 && !std::is_same<typename std::decay_t<T1>::Signature,
                                                  typename std::decay_t<T2>::Signature>::value
                 ), int> = 0>
      constexpr bool operator==(const T1&, const T2&)
      {
        return false;
      }
 
      template<class T1, class T2,
               std::enable_if_t<(AreProperTensors<T1, T2>::value
                                 && AreRuntimeEqualOperands<T1, T2>::value
                 ), int> = 0>
      constexpr bool operator==(const T1&, const T2&)
      {
        return true;
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr bool operator!=(T1&& t1, T2&& t2)
      {
        return !(t1 == t2);
      }
 
      namespace {
 
        template<class T1, class T2, class F,
                 std::enable_if_t<(std::decay_t<T1>::dim() > Policy::TemplateForEachLimit::value), int> = 0>
        constexpr bool compareHelper(T1&& t1, T2&& t2, F&& f)
        {
          for (unsigned i = 0; i < TensorTraits<T1>::dim(); ++i) {
            if (!f(tensorValue(t1, multiIndex(i, TensorTraits<T1>::signature())), tensorValue(t2, multiIndex(i, TensorTraits<T2>::signature())))) {
              return false;
            }
          }
          return true;
        }
 
        template<class T1, class T2, class F,
                 std::enable_if_t<std::decay_t<T1>::dim() <= Policy::TemplateForEachLimit::value, int> = 0>
        constexpr bool compareHelper(T1&& t1, T2&& t2, F&& f)
        {
          return forEachWhile(MakeIndexSequence<TensorTraits<T1>::dim()>{}, [&](auto i) {
              using I = decltype(i);
              return f(tensorValue(t1, multiIndex<I::value>(TensorTraits<T1>::signature())),
                       tensorValue(t2, multiIndex<I::value>(TensorTraits<T2>::signature())));
            });
        }
 
      } // NS anonymous
 
      template<class T1, class T2, std::enable_if_t<AreRuntimeComparable<T1, T2>::value, int> = 0>
      constexpr bool operator==(T1&& t1, T2&& t2)
      {
        return compareHelper(std::forward<T1>(t1), std::forward<T2>(t2), [](auto&& a, auto&& b) {
            return a == b;
          });
      }
 
      template<class T1, class T2, std::enable_if_t<AreRuntimeComparable<T1, T2>::value, int> = 0>
      constexpr bool operator<(T1&& t1, T2&& t2)
      {
        return compareHelper(std::forward<T1>(t1), std::forward<T2>(t2), [](auto&& a, auto&& b) {
            return a < b;
          });
      }
 
      template<class T1, class T2, std::enable_if_t<AreRuntimeComparable<T1, T2>::value, int> = 0>
      constexpr bool operator<=(T1&& t1, T2&& t2)
      {
        return compareHelper(std::forward<T1>(t1), std::forward<T2>(t2), [](auto&& a, auto&& b) {
            return a <= b;
          });
      }
 
      template<class T1, class T2, std::enable_if_t<AreRuntimeComparable<T1, T2>::value, int> = 0>
      constexpr bool operator>(T1&& t1, T2&& t2)
      {
        return compareHelper(std::forward<T1>(t1), std::forward<T2>(t2), [](auto&& a, auto&& b) {
            return a > b;
          });
      }
 
      template<class T1, class T2, std::enable_if_t<AreRuntimeComparable<T1, T2>::value, int> = 0>
      constexpr bool operator>=(T1&& t1, T2&& t2)
      {
        return compareHelper(std::forward<T1>(t1), std::forward<T2>(t2), [](auto&& a, auto&& b) {
            return a >= b;
          });
      }
 
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto sqrt(T&& t)
      {
        return finalize<SqrtOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto sqr(T&& t)
      {
        return finalize<SquareOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto cos(T&& t)
      {
        return finalize<CosOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto sin(T&& t)
      {
        return finalize<SinOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto tan(T&& t)
      {
        return finalize<TanOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto acos(T&& t)
      {
        return finalize<AcosOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto asin(T&& t)
      {
        return finalize<AsinOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto atan(T&& t)
      {
        return finalize<AtanOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto erf(T&& t)
      {
        return finalize<ErfOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto lgamma(T&& t)
      {
        return finalize<LGammaOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto tgamma(T&& t)
      {
        return finalize<TGammaOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto exp(T&& t)
      {
        return finalize<ExpOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto log(T&& t)
      {
        return finalize<LogOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto cosh(T&& t)
      {
        return finalize<CoshOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto sinh(T&& t)
      {
        return finalize<SinhOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto tanh(T&& t)
      {
        return finalize<TanhOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto acosh(T&& t)
      {
        return finalize<AcoshOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto asinh(T&& t)
      {
        return finalize<AsinhOperation>(std::forward<T>(t));
      }
 
      template<class T, std::enable_if_t<IsProperTensor<T>::value, int> = 0>
      auto atanh(T&& t)
      {
        return finalize<AtanhOperation>(std::forward<T>(t));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto atan2(T1&& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::signature() == TensorTraits<T2>::signature(),
                      "Tensor signatures must coincide for component-wise operations.");
 
        return finalize<Atan2Operation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<
        class T1, class T2,
        std::enable_if_t<(AreProperTensors<T1, T2>::value
                          && TensorTraits<T1>::rank != 0
                          && TensorTraits<T2>::rank == 0
        ), int> = 0>
      auto pow(T1&& t1, T2&& t2)
      {
        return finalize<PowOperation>(
          std::forward<T1>(t1),
          operate(
            ScalarEinsumFunctor{},
            std::forward<T2>(t2),
            asExpression(ones<T1>()))
          );
      }
 
      template<
        class T1, class T2,
        std::enable_if_t<(AreProperTensors<T1, T2>::value
                          && TensorTraits<T1>::rank == TensorTraits<T2>::rank
        ), int> = 0>
      auto pow(T1&& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::signature() == TensorTraits<T2>::signature(),
                      "Tensor signatures must coincide for component-wise operations.");
        return finalize<PowOperation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto min(T1&& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::signature() == TensorTraits<T2>::signature(),
                      "Tensor signatures must coincide for component-wise operations.");
 
        return finalize<MinOperation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2, std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      auto max(T1&& t1, T2&& t2)
      {
        static_assert(TensorTraits<T1>::signature() == TensorTraits<T2>::signature(),
                      "Tensor signatures must coincide for component-wise operations.");
 
        return finalize<MaxOperation>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      //If, then else construct
      template<class T1, std::enable_if_t<IsProperTensor<T1>::value, int> = 0>
      auto ternary(T1&& t1)
      {
        return finalize<TernaryOperation>(std::forward<T1>(t1));
      }
 
 
    } // NS Tensor
 
#if 1
    using Tensor::operator-;
    using Tensor::operator+;
    using Tensor::operator*;
    using Tensor::operator/;
    using Tensor::operator&;
 
    using Tensor::operator==;
    using Tensor::operator!=;
    using Tensor::operator>=;
    using Tensor::operator<=;
    using Tensor::operator>;
    using Tensor::operator<;
#endif
 
 
 
  } // NS ACFem
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_EXPRESSIONS_HH__