operatorfunctional.hh

#ifndef __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_OPERATORFUNCTIONAL_HH__
#define __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_OPERATORFUNCTIONAL_HH__
 
#include <dune/fem/operator/common/operator.hh>
 
#include "../../../common/literals.hh"
#include "../../../expressions/terminal.hh"
 
#include "../linearfunctional.hh"
#include "../functionaltraits.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace LinearFunctional {
 
      using namespace Literals;
 
      template<class DomainFunction, class Operator>
      class OperatorFunctional
        : public LinearFunctionalBase<typename std::decay_t<Operator>::RangeFunctionType::DiscreteFunctionSpaceType>
        , public Expressions::SelfExpression<OperatorFunctional<DomainFunction, Operator> >
      {
       public:
        using OperatorType = std::decay_t<Operator>;
        using DiscreteFunctionSpaceType = typename OperatorType::RangeFunctionType::DiscreteFunctionSpaceType;
       private:
        using BaseType = LinearFunctionalBase<DiscreteFunctionSpaceType>;
       public:
        using typename BaseType::RangeFieldType;
        using DomainFunctionType = typename OperatorType::DomainFunctionType;
 
        static_assert(std::is_constructible<DomainFunctionType, DomainFunction>::value,
                      "Incompatible DomainFunction");
 
        template<class DomFct, class Op,
                 std::enable_if_t<(std::is_constructible<Operator, Op>::value
                                   && std::is_constructible<DomainFunction, DomFct>::value
          ), int> = 0>
        OperatorFunctional(Op&& op, DomFct&& domFct, const DiscreteFunctionSpaceType& space)
          : BaseType(space)
          , domainFunction_(std::forward<DomFct>(domFct)), operator_(std::forward<Op>(op))
        {}
 
        template<class FieldArg, class DiscreteFunction,
                 std::enable_if_t<std::is_constructible<RangeFieldType, FieldArg>::value, int> = 0>
        void coefficients(FieldArg&& s, DiscreteFunction& values) const
        {
          if constexpr (ExpressionTraits<FieldArg>::isZero) {
            values.clear();
          } else {
            operator_(domainFunction_, values);
            if constexpr (!ExpressionTraits<FieldArg>::isOne) {
              values *= s;
            }
          }
        }
 
        template<class DiscreteFunction>
        void coefficients(DiscreteFunction& df) const
        {
          coefficients(1_f, df);
        }
 
        template<class DiscreteFunction>
        auto operator()(const DiscreteFunction& df) const
        {
          return inner(*this, df);
        }
 
        std::string name() const
        {
          return "Op(" + domainFunction_.name() + ",.)";
        }
 
       protected:
        DomainFunction domainFunction_;
        Operator operator_;
      };
 
      template<class PDEOperator, class DomainFunction>
      auto operatorFunctional(PDEOperator&& op,
                              DomainFunction&& domFct,
                              const typename std::decay_t<PDEOperator>::RangeFunctionType::DiscreteFunctionSpaceType& space)
      {
        static_assert(IsDiscreteFunction<typename std::decay_t<PDEOperator>::RangeFunctionType>::value
                      && std::is_constructible<typename std::decay_t<PDEOperator>::DomainFunctionType, DomainFunction>::value,
          "DomainFunction type is incompatibel with given PDEOperator type.");
 
        return OperatorFunctional<DomainFunction, PDEOperator>(
          std::forward<PDEOperator>(op),
          std::forward<DomainFunction>(domFct),
          space
          );
      }
 
    } // LinearFunctional::
 
    using LinearFunctional::operatorFunctional;
 
  } // ACFem::
 
} // Dune::
 
#endif // __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_OPERATORFUNCTIONAL_HH__