bindabletensorfunction.hh

#ifndef __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
#define __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
 
#include <dune/fem/function/localfunction/bindable.hh>
 
#include "../common/ostream.hh"
#include "../common/quadraturepoint.hh"
#include "../expressions/examine.hh"
#include "../expressions/polynomialtraits.hh"
#include "../tensors/optimization/subexpression.hh"
#include "../tensors/bindings/dune/fieldvectortraits.hh"
#include "../tensors/ostream.hh"
 
#include "policy.hh"
#include "placeholders/placeholdertraits.hh"
#include "placeholders/localfunctionderivatives.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace GridFunction {
 
      using namespace Literals;
 
      using Tensor::indeterminate;
      using Tensor::TensorTraits;
      using Tensor::derivative;
      using Tensor::Seq;
      using Expressions::replaceIndeterminate;
      using Expressions::IsClosure;
      using Expressions::extract;
      using Expressions::extractPlaceholders;
      using Expressions::extractIndeterminates;
      using Expressions::PlaceholderTuple;
      using Expressions::ExplicitlyDependsOn;
      using Expressions::asExpression;
      using Expressions::finalize;
      using Expressions::EnclosedType;
      using Expressions::AreRuntimeEqual;
      using Expressions::ExamineOr;
 
      template<class Expr, class GridPart, std::size_t AutoDiffOrder, std::size_t IndeterminateId>
      class BindableTensorFunction;
 
      template<class Expr, class GridPart, std::size_t AutoDiffOrder = Policy::autoDiffLimit(), std::size_t IndeterminateId = Policy::indeterminateId()>
      class BindableTensorFunction
        : public Fem::BindableGridFunction<GridPart, Tensor::FieldVectorVectorType<FloatingPointClosure<typename TensorTraits<Expr>::FieldType>, typename TensorTraits<Expr>::Signature> >
        , public Expressions::SelfExpression<BindableTensorFunction<Expr, GridPart, AutoDiffOrder, IndeterminateId> >
      {
       public:
        using FieldType = FloatingPointClosure<typename TensorTraits<Expr>::FieldType>;
        using TensorSignature = typename TensorTraits<Expr>::Signature;
        using RangeType = typename Tensor::FieldVectorVectorType<FieldType, TensorSignature>;
       private:
        using BaseType = Fem::BindableGridFunction<GridPart, RangeType>;
       public:
        using DomainType = typename BaseType::DomainType;
        using JacobianRangeType = typename BaseType::JacobianRangeType;
        using HessianRangeType = typename BaseType::HessianRangeType;
        using EntityType = typename BaseType::EntityType;
 
        using BaseType::global;
        using BaseType::entity;
        using BaseType::gridPart;
 
        static constexpr std::size_t indeterminateId_ = IndeterminateId; 
 
        static_assert(!ExamineOr<Expr, Tensor::IsAutoDiffPlaceholder>::value,
                      "Expr should not contain AD placeholders.");
 
        static constexpr std::size_t requestedADOrder_ = AutoDiffOrder;
        static constexpr std::size_t autoDiffOrder_ =
          (hasConstantPolynomialDegreeAtMost<0, Expr, IndeterminateId>()
           ? 0
           : (hasConstantPolynomialDegreeAtMost<1, Expr, IndeterminateId>()
              ? std::min(AutoDiffOrder, 1UL)
              : AutoDiffOrder));
 
        using IndeterminateType =
          std::decay_t<decltype(asExpression(Tensor::indeterminate<IndeterminateId>(std::declval<DomainType>())))>;
 
        using ValuesType = decltype(
          replaceIndeterminate<IndeterminateId>(std::declval<Expr>(), std::declval<IndeterminateType&>())
          );
 
        static_assert(!std::is_rvalue_reference<ValuesType>::value,
                      "Expr must not be an rvalue reference.");
 
        using ValuesCSECascade = std::decay_t<decltype(
          Tensor::autoDiffCSECascade<autoDiffOrder_>(std::declval<ValuesType&&>(), std::declval<IndeterminateType>())
          )>;
 
        using ValuesCSEType = std::decay_t<decltype(
          injectSubExpressions(std::declval<ValuesType&&>(), std::declval<ValuesCSECascade>())
          )>;
 
        BindableTensorFunction(Expr&& t, const GridPart& gridPart)
          : BaseType(gridPart)
          , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType(0))))
          , expression_(replaceIndeterminate<indeterminateId_>(std::forward<Expr>(t), indeterminate_))
          , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
          , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
          , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
          , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
          , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
          , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
          , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
          , name_(expression_.name())
        {}
 
        BindableTensorFunction(const BindableTensorFunction& other)
          : BaseType(other)
          , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType())))
          , expression_(replaceIndeterminate<indeterminateId_>(const_cast<ValuesType&&>(other.expression_), indeterminate_))
          , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
          , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
          , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
          , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
          , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
          , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
          , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
          , name_(other.name_)
        {}
 
        BindableTensorFunction(BindableTensorFunction&& other)
          : BaseType(other)
          , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType())))
          , expression_(replaceIndeterminate<indeterminateId_>(const_cast<ValuesType&&>(other.expression_), indeterminate_))
          , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
          , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
          , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
          , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
          , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
          , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
          , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
          , name_(std::move(other.name_))
        {}
 
        void bind(const EntityType& entity)
        {
          BaseType::bind(entity);
          forEach(bindablePlaceholders_, [&](auto&& p) {
              p.bind(entity);
            });
        }
 
        void unbind()
        {
          forEach(unbindablePlaceholders_, [&](auto&& p) {
              p.unbind();
            });
          BaseType::unbind();
#ifndef NDEBUG
          // in order to catch bogus bind cascades
          BaseType::bind(EntityType());
#endif
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        const auto& operator()(const Point& point) const
        {
          initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
          initPlaceholders(point, valueSetValuePlaceholders_);
          initPlaceholders(valueIndeterminateSetValuePlaceholders_);
 
          valuesCSECascade_.evaluate(valuesCSECascade_.sExpTraits().upTo(0_c));
 
          return values_;
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void evaluate(const Point& point, RangeType& result) const
        {
          result = (*this)(point);
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void jacobian(const Point& point, JacobianRangeType& result) const
        {
          if constexpr (autoDiffOrder_ < 1) {
            result = 0;
          } else {
            initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
            initPlaceholders(point, valueSetValuePlaceholders_);
            initPlaceholders(valueIndeterminateSetValuePlaceholders_);
 
            const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(1_c);
            valuesCSECascade_.evaluate(evalTraits);
 
            // The call to derivative() will just peek out the pre-evaluated derivative
            assign(result, derivative(evalTraits(values_), indeterminate_));
          }
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void hessian(const Point& point, HessianRangeType& result) const
        {
          if constexpr (autoDiffOrder_ < 2) {
            result = 0;
          } else {
            static_assert(autoDiffOrder_ >= 2, "We were constructed without 2nd derivatives ...");
 
            initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
            initPlaceholders(point, valueSetValuePlaceholders_);
            initPlaceholders(valueIndeterminateSetValuePlaceholders_);
 
            const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(2_c);
            valuesCSECascade_.evaluate(evalTraits);
 
            // The call to derivative() will just peek out the pre-evaluated derivative
            assign(result, derivative<2>(evalTraits(values_), indeterminate_));
          }
        }
 
        ValuesType expression()&&
        {
          return expression_;
        }
 
        decltype(auto) expression()&
        {
          return expression_;
        }
 
        decltype(auto) expression() const&
        {
          return expression_;
        }
 
        template<std::size_t I>
        auto operator[](const IndexConstant<I>)&
        {
          return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
        }
 
        template<std::size_t I>
        auto operator[](const IndexConstant<I>)const&
        {
          return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
        }
 
        template<std::size_t I>
        auto operator[](const IndexConstant<I>)const&&
        {
          return gridFunction(
            gridPart(),
            restriction<0>(std::move(values_), Seq<I>{}),
            requestedADOrder(),
            indeterminateId()
            );
        }
 
        constexpr unsigned int order() const
        {
          return Expressions::polynomialDegree<IndeterminateId>(values_);
        }
 
        std::string expressionName() const
        {
          return expression_.name();
        }
 
        std::string name() const
        {
          return name_;
        }
 
        void setName(const std::string & name)
        {
          name_ = name;
        }
 
        static constexpr auto indeterminateId()
        {
          return IndexConstant<indeterminateId_>{};
        }
 
        static constexpr auto autoDiffOrder()
        {
          return IndexConstant<(hasConstantPolynomialDegreeOf<0, Expr, IndeterminateId>()
                                ? Policy::autoDiffLimit()
                                : (hasConstantPolynomialDegreeOf<1, Expr, IndeterminateId>() && autoDiffOrder_ >= 1
                                   ? Policy::autoDiffLimit()
                                   : autoDiffOrder_))>{};
        }
 
        static constexpr auto requestedADOrder()
        {
          return IndexConstant<requestedADOrder_>{};
        }
 
        decltype(auto) indeterminate() const& {
          return indeterminate_;
        }
 
        decltype(auto) indeterminate()& {
          return indeterminate_;
        }
 
        IndeterminateType indeterminate()&& {
          return indeterminate_;
        }
 
        static constexpr auto signature()
        {
          return typename TensorTraits<ValuesType>::Signature{};
        }
 
        void printInfo(std::ostream& out = std::clog) const
        {
          //std::clog << "Memory Footprint: " << sizeof(*this) << std::endl;
          std::clog << "IndeterminateType: " << typeString<IndeterminateType>() << std::endl;
          forLoop<3>([](auto i) {
              using I = decltype(i);
              std::clog << "Constant Degree of " << I::value << ": " << hasConstantPolynomialDegreeOf<I::value, ValuesType, IndeterminateId>() << std::endl;
                     });
          std::clog << "Estimated order: " << order() << std::endl;
          std::clog << "Value Placeholders: "
                    << " #" << size(valuePlaceholders_) << ": "
                    << typeString<decltype(valuePlaceholders_)>() << std::endl;
          std::clog << "  bindable:   "
                    << " #" << size(bindablePlaceholders_) << ": "
                    << typeString<decltype(bindablePlaceholders_)>() << std::endl;
          std::clog << "  unbindable: "
                    << " #" << size(unbindablePlaceholders_) << ": "
                    << typeString<decltype(unbindablePlaceholders_)>() << std::endl;
          std::clog << "  set value:  "
                    << " #" << size(valueSetValuePlaceholders_) << ": "
                    << typeString<decltype(valueSetValuePlaceholders_)>() << std::endl;
          std::clog << "  set indet:  "
                    << " #" << size(valueIndeterminateSetValuePlaceholders_) << ": "
                    << typeString<decltype(valueIndeterminateSetValuePlaceholders_)>() << std::endl;
          std::clog << "ValuesType:" << std::endl
                    << "  " << values_.name() << std::endl
                    << "  " << typeString<decltype(values_)>() << std::endl;
          std::clog << "  depends explicitly on " << indeterminate_.name() << ": " << ExplicitlyDependsOn<indeterminateId_, ValuesType>::value << std::endl;
          std::clog << "  is zero: " << ExpressionTraits<ValuesType>::isZero << std::endl;
          std::clog << "  tensor signature: " << TensorTraits<ValuesType>::signature() << std::endl;
          std::clog << "  fem signature:    " << TensorTraits<RangeType>::signature() << std::endl;
          std::clog << "  CSE expressions: " << std::endl;
          std::size_t count = 0;
          forLoop<ValuesCSECascade::depth()>([this,&count](auto i) {
              forEach(valuesCSECascade_.template get<i.value>(),[&count](const auto& t) {
                  ++count;
                  std::clog << t.name() << std::endl;
                });
            });
          std::clog << "  #CSE expressions: " << count << std::endl;
        }
 
       private:
        template<class FemResult, class TensorResult>
        struct NeedsScalarDisambiguation
          : BoolConstant<((TensorTraits<FemResult>::rank - TensorTraits<TensorResult>::rank) == 1
                          && TensorTraits<FemResult>::template dim<0>() == 1)>
        {};
 
        template<class FemResult, class TensorResult,
                 std::enable_if_t<(NeedsScalarDisambiguation<FemResult, TensorResult>::value), int> = 0>
        static void assign(FemResult& to, TensorResult&& from)
        {
          tensor(to[0]) = from;
        }
 
        template<class FemResult, class TensorResult,
                 std::enable_if_t<(!NeedsScalarDisambiguation<FemResult, TensorResult>::value), int> = 0>
        static void assign(FemResult& to, TensorResult&& from)
        {
          tensor(to) = from;
        }
 
        template<
          class Point, class IsDependent,
          std::enable_if_t<((IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value)
                            && IsDependent::value), int> = 0>
        void initIndeterminate(const Point& point, IsDependent) const
        {
          //assert(entity() != EntityType());
          indeterminate_.operand(0_c) = global(point);
        }
 
        template<
          class Point, class IsDependent,
          std::enable_if_t<((IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value)
                            && !IsDependent::value), int> = 0>
        void initIndeterminate(const Point& point, IsDependent) const
        {}
 
        template<class T>
        struct IsRelevantPlaceholder
          : BoolConstant<(IsPlaceholderExpression<T>::value
                          && !Tensor::IsAutoDiffPlaceholder<T>::value
                          && !RefersConst<T>::value
          )>
        {};
 
        template<class T, class Cascade>
        static auto extractPlaceholders(T&& expr, Cascade&& cascade)
        {
          return
            extract<IsRelevantPlaceholder>(
              std::forward<T>(expr),
              extract<IsRelevantPlaceholder>(
                std::forward<Cascade>(cascade),
                std::tuple<>{},
                FalseType{}, // recursive
                PredicateWrapper<AreRuntimeEqual>{}
                ),
              FalseType{}, // recursive
              PredicateWrapper<AreRuntimeEqual>{}
              );
        }
 
        template<class T, class Cascade>
        using PlaceholderTuple = std::decay_t<decltype(extractPlaceholders(std::declval<T>(), std::declval<Cascade>()))>;
 
        template<class Point, class Placeholders,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        static void initPlaceholders(const Point& point, const Placeholders& placeholders)
        {
          forLoop<size<Placeholders>()>([&](auto i) {
              get<i.value>(placeholders).setValue(point);
            });
        }
 
        template<class Placeholders>
        void initPlaceholders(const Placeholders& placeholders) const
        {
          forLoop<size<Placeholders>()>([&](auto i) {
              get<i.value>(placeholders).setValue(indeterminate_);
            });
        }
 
        mutable IndeterminateType indeterminate_;
 
        ValuesType expression_;
 
        mutable ValuesCSECascade valuesCSECascade_;
        ValuesCSEType values_;
 
        PlaceholderTuple<ValuesCSEType&, ValuesCSECascade&> valuePlaceholders_;
        FilteredTuple<HasBind, decltype(valuePlaceholders_)> bindablePlaceholders_;
        FilteredTuple<HasUnbind, decltype(valuePlaceholders_)> unbindablePlaceholders_;
 
        FilteredTuple<HasSetValue, decltype(valuePlaceholders_)> valueSetValuePlaceholders_;
 
        FilteredTuple<HasIndeterminateSetValue, decltype(valuePlaceholders_)> valueIndeterminateSetValuePlaceholders_;
        std::string name_;
      };
 
 
    } // NS GridFunction
 
  } // NS ACFem
 
} // NS Dune
 
#endif // __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__