divergence.hh

#ifndef __DUNE_ACFEM_TENSOR_OPERATIONS_DIVERGENCE_HH__
#define __DUNE_ACFEM_TENSOR_OPERATIONS_DIVERGENCE_HH__
 
#include "derivative.hh"
 
namespace Dune::ACFem {
 
  namespace Tensor {
 
    template<std::size_t IndeterminateId, std::size_t... Dims, class T,
             std::enable_if_t<IsTensorOperand<T>::value, int> = 0>
    constexpr decltype(auto) divergence(T&& t, Seq<Dims...> = Seq<Dims...>{}, IndexConstant<IndeterminateId> = IndexConstant<IndeterminateId>{})
    {
      using IndeterminateSignature = Seq<Dims...>;
      constexpr std::size_t indeterminateRank = sizeof...(Dims);
      using Signature = typename TensorTraits<T>::Signature;
      using DivTailSignature = TailPart<indeterminateRank, Signature>;
 
      static_assert(std::is_same<DivTailSignature, IndeterminateSignature>::value,
                    "Tail of signature does not match signature of indeterminate.");
 
      constexpr std::size_t offset = TensorTraits<T>::rank - indeterminateRank;
      using PosLeft = MakeIndexSequence<2*indeterminateRank, offset>;
      using PosRight = MakeIndexSequence<2*indeterminateRank>;
 
      return einsum<PosLeft, PosRight>(derivative<IndeterminateId, Dims...>(std::forward<T>(t)), blockEye<2>(IndeterminateSignature{}, Expressions::Disclosure{}));
    }
 
    template<
      class T, class X,
      std::enable_if_t<(IsTensorOperand<T>::value
                        && IsTensor<X>::value
                        && IsIndeterminateExpression<EnclosedType<X> >::value
      ), int> = 0>
    constexpr decltype(auto) divergence(T&& t, X&&)\
    {
      using XType = EnclosedType<X>;
      constexpr std::size_t id = IndeterminateTraits<XType>::id_;
      using Signature = typename TensorTraits<XType>::Signature;
 
      return divergence<id>(std::forward<T>(t), Signature{});
    }
 
  } // Tensor::
 
} // Dune::ACFEM::
 
#endif // __DUNE_ACFEM_TENSOR_OPERATIONS_DIVERGENCE_HH__