functionexpression.hh

#ifndef __DUNE_ACFEM_FUNCTIONEXPRESSION_HH__
#define __DUNE_ACFEM_FUNCTIONEXPRESSION_HH__
 
#include <dune/fem/version.hh>
#include <dune/fem/function/common/function.hh>
 
#include "../functions/constantfunction.hh"
#include "../functions/parameterfunction.hh"
#include "../functions/functionoperations.hh"
#include "../expressions/parameterinterface.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    template<class FunctionSpaceImp, class ExpressionImp>
    class FunctionExpression
      : public Fem::Function<FunctionSpaceImp, ExpressionImp>,
        public ExpressionTemplate<ExpressionImp>
    {};
 
    template<class UnOp, class FunctionType>
    class UnaryFunctionExpression;
 
    template<class BinOp, class LeftFunctionType, class RightFunctionType>
    class BinaryFunctionExpression;
 
 
    /******************************************************************************
     *
     * Traits and helpers
     *
     */
 
    /* Determine the resulting function space */
    template<class LeftSpace, class RightSpace>
    struct FunctionMultiplicationResultTraits
    {
      typedef LeftSpace LeftFunctionSpaceType;
      typedef RightSpace RightFunctionSpaceType;
 
      typedef typename LeftFunctionSpaceType::ScalarFunctionSpaceType LeftScalarSpaceType;
      typedef typename RightFunctionSpaceType::ScalarFunctionSpaceType RightScalarSpaceType;
 
      static_assert(std::is_same<LeftScalarSpaceType, RightScalarSpaceType>::value,
                    "Scalar functions spaces do not coincide");
 
      enum {
        sameSpace = std::is_same<LeftFunctionSpaceType, RightFunctionSpaceType>::value,
        leftScalar = std::is_same<LeftFunctionSpaceType, LeftScalarSpaceType>::value,
        rightScalar = std::is_same<RightFunctionSpaceType, RightScalarSpaceType>::value
      };
 
      static_assert(sameSpace || leftScalar || rightScalar,
                    "Only S-multiplication and scalar products are supported");
 
      // Scalar products result in scalars, s-mult results in vectors/
      typedef typename
      std::conditional<sameSpace,
                       LeftScalarSpaceType,
                       typename std::conditional<leftScalar,
                                                 RightFunctionSpaceType,
                                                 LeftFunctionSpaceType>::type>::type
      FunctionSpaceType;
    };
 
    template<class BinOp, class LeftFunctionType, class RightFunctionType>
    struct BinaryFunctionExpressionTraits;
 
    template<class LeftFunction, class RightFunction>
    struct BinaryFunctionExpressionTraits<MultiplyOperation, LeftFunction, RightFunction>
    {
      typedef LeftFunction LeftFunctionType;
      typedef RightFunction RightFunctionType;
      typedef typename LeftFunctionType::FunctionSpaceType LeftFunctionSpaceType;
      typedef typename RightFunctionType::FunctionSpaceType RightFunctionSpaceType;
 
      typedef typename
      FunctionMultiplicationResultTraits<LeftFunctionSpaceType, RightFunctionSpaceType>::FunctionSpaceType
      FunctionSpaceType;
    };
 
    template<class Factor, class Function>
    struct BinaryFunctionExpressionTraits<SMultiplyOperation, Factor, Function>
    {
      typedef Function FunctionType;
      typedef Factor FactorType;
      typedef typename FunctionType::FunctionSpaceType FunctionSpaceType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
    };
 
    template<class LeftFunction, class RightFunction>
    struct BinaryFunctionExpressionTraits<PlusOperation, LeftFunction, RightFunction>
    {
      typedef LeftFunction LeftFunctionType;
      typedef RightFunction RightFunctionType;
      typedef typename LeftFunctionType::FunctionSpaceType LeftFunctionSpaceType;
      typedef typename RightFunctionType::FunctionSpaceType RightFunctionSpaceType;
 
      static_assert(std::is_same<LeftFunctionSpaceType, RightFunctionSpaceType>::value,
                    "Function spaces must agree, cannot combine them");
 
      typedef LeftFunctionSpaceType FunctionSpaceType;
    };
 
    template<class LeftFunction, class RightFunction>
    struct BinaryFunctionExpressionTraits<MinusOperation, LeftFunction, RightFunction>
      : public BinaryFunctionExpressionTraits<PlusOperation, LeftFunction, RightFunction>
    {};
 
    // BinaryFunctionExpression
    // -------------------
 
    template<class BinOp, class LeftFunction, class RightFunction>
    class BinaryFunctionExpression
      : public FunctionExpression<typename
                                  BinaryFunctionExpressionTraits<BinOp, LeftFunction, RightFunction>::FunctionSpaceType,
                                  BinaryFunctionExpression<BinOp, LeftFunction, RightFunction> >
    {
      typedef LeftFunction LeftFunctionType;
      typedef RightFunction RightFunctionType;
     public:
      typedef BinaryFunctionExpressionTraits<BinOp, LeftFunctionType, RightFunctionType> TraitsType;
     private:
      typedef BinaryFunctionExpression<BinOp, LeftFunctionType, RightFunctionType> ThisType;
      typedef Fem::Function<typename TraitsType::FunctionSpaceType, ThisType> BaseType;
      typedef BinOp OperationTagType;
      typedef BinaryFunctionExpressionOperation<BinOp> OperationType;
     public:
      // type of discrete function space
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
     public:
      // reference to function this local belongs to
      BinaryFunctionExpression(const LeftFunctionType &f,
                               const RightFunctionType &g)
        : leftFunction_(f), rightFunction_(g)
      {}
 
      // We DO allow copying
      BinaryFunctionExpression(const ThisType &other)
        : leftFunction_(other.leftFunction_),
          rightFunction_(other.rightFunction_)
      {}
 
      void evaluate(const DomainType &global, RangeType &result) const
      {
        OperationType::evaluate(leftFunction_(), rightFunction_(), global, result);
      }
 
      void jacobian(const DomainType &global, JacobianRangeType &result) const
      {
        OperationType::jacobian(leftFunction_(), rightFunction_(), global, result);
      }
 
      void hessian(const DomainType &global, HessianRangeType &result) const
      {
        OperationType::jacobian(leftFunction_(), rightFunction_(), global, result);
      }
 
     public:
      const LeftFunctionType& leftFunction() const { return leftFunction_(); }
      const RightFunctionType& rightFunction() const { return rightFunction_(); }
 
     private:
      ExpressionStorage<LeftFunctionType> leftFunction_;
      ExpressionStorage<RightFunctionType> rightFunction_;
    };
 
    // BinaryFunctionExpression
    // -------------------
 
    template<class FactorType, class FunctionType>
    class BinaryFunctionExpression<SMultiplyOperation, FactorType, FunctionType>
      : public FunctionExpression<typename FunctionType::FunctionSpaceType,
                                  BinaryFunctionExpression<SMultiplyOperation, FactorType, FunctionType> >
    {
      typedef BinaryFunctionExpression<SMultiplyOperation, FactorType, FunctionType> ThisType;
      typedef Fem::Function<typename FunctionType::FunctionSpaceType, ThisType> BaseType;
      typedef SMultiplyOperation OperationTagType;
      typedef BinaryFunctionExpressionOperation<SMultiplyOperation> OperationType;
     public:
      typedef BinaryFunctionExpressionTraits<SMultiplyOperation, FactorType, FunctionType> TraitsType;
 
      // type of discrete function space
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
     public:
      // reference to function this local belongs to
      BinaryFunctionExpression(const FactorType &s, const FunctionType &f)
        : scalar_(s),
          function_(f)
      {}
 
      // We DO allow copying of expressions
      BinaryFunctionExpression(const ThisType &other)
        : scalar_(other.scalar_),
          function_(other.function_)
      {}
 
      void evaluate(const DomainType &global, RangeType &result) const
      {
        function_().evaluate(global, result);
        result *= parameterValue(scalar_());
      }
 
      void jacobian(const DomainType &global, JacobianRangeType &result) const
      {
        function_().jacobian(global, result);
        result *= parameterValue(scalar_());
      }
 
      void hessian(const DomainType &global, HessianRangeType &result) const
      {
        function_().hessian(global, result);
        result *= parameterValue(scalar_());
      }
 
     public:
      const FactorType& scalar() const { return scalar_(); }
      const FunctionType& function() const { return function_(); }
 
     protected:
      ExpressionStorage<FactorType>   scalar_;
      ExpressionStorage<FunctionType> function_;
    };
 
    //
    // Unary expressions
 
    // The result space may depend on the operation. The default is to
    // have the result space equal to the function space (inversion,
    // negation etc.)
    template<class UnOp, class FunctionType>
    struct UnaryFunctionOperationTraits
    {
      typedef typename FunctionType::FunctionSpaceType FunctionSpaceType;
    };
 
    template<class FunctionType>
    struct UnaryFunctionOperationTraits<SquareOperation, FunctionType>
    {
      typedef typename FunctionType::FunctionSpaceType::ScalarFunctionSpaceType FunctionSpaceType;
    };
 
    template <class UnOp, class FunctionType>
    struct UnaryFunctionExpressionTraits
    {
      typedef UnaryFunctionOperationTraits<UnOp, FunctionType> OperationTraitsType;
      typedef typename OperationTraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
    };
 
    // UnaryFunctionExpression
    // -------------------
 
    template<class UnOp, class FunctionType>
    class UnaryFunctionExpression
      : public FunctionExpression<typename UnaryFunctionExpressionTraits<UnOp, FunctionType>::FunctionSpaceType,
                                  UnaryFunctionExpression<UnOp, FunctionType> >
    {
     public:
      typedef UnaryFunctionExpressionTraits<UnOp, FunctionType> TraitsType;
     private:
      typedef UnaryFunctionExpression<UnOp, FunctionType> ThisType;
      typedef Fem::Function<typename TraitsType::FunctionSpaceType, ThisType> BaseType;
      typedef UnOp OperationTagType;
      typedef UnaryFunctionExpressionOperation<UnOp> OperationType;
     public:
      // type of discrete function space
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      // reference to function this local belongs to
      UnaryFunctionExpression(const FunctionType &f)
        : function_(f)
      {}
 
      // reference to function this local belongs to
      UnaryFunctionExpression(const ThisType &other)
        : function_(other.function_)
      {}
 
      void evaluate(const DomainType &global, RangeType &result) const
      {
        OperationType::evaluate(function_(), global, result);
      }
 
      void jacobian(const DomainType &global, JacobianRangeType &result) const
      {
        OperationType::jacobian(function_(), global, result);
      }
 
      void hessian(const DomainType &global, HessianRangeType &result) const
      {
        OperationType::hessian(function_(), global, result);
      }
 
     public:
      const FunctionType& function() const { return function_(); }
 
     private:
      ExpressionStorage<FunctionType> function_;
    };
 
    // UnaryFunctionExpression<IdentityOperation>
    // -------------------
 
    template<class FunctionType>
    class UnaryFunctionExpression<IdentityOperation, FunctionType>
      : public FunctionExpression<typename FunctionType::FunctionSpaceType,
                                  UnaryFunctionExpression<IdentityOperation, FunctionType> >
    {
      typedef UnaryFunctionExpression<IdentityOperation, FunctionType> ThisType;
      typedef Fem::Function<typename FunctionType::FunctionSpaceType, ThisType> BaseType;
      typedef IdentityOperation OperationTagType;
      typedef UnaryFunctionExpressionOperation<OperationTagType> OperationType;
     public:
      typedef UnaryFunctionExpressionTraits<OperationTagType, FunctionType> TraitsType;
 
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      // reference to function this local belongs to
      UnaryFunctionExpression(const FunctionType &f)
        : function_(f)
      {}
 
      // reference to function this local belongs to
      UnaryFunctionExpression(const ThisType &other)
        : function_(other.function_)
      {}
 
      void evaluate(const DomainType &global, RangeType &result) const
      {
        function_().evaluate(global, result);
      }
 
      void jacobian(const DomainType &global, JacobianRangeType &result) const
      {
        function_().jacobian(global, result);
      }
 
      void hessian(const DomainType &global, HessianRangeType &result) const
      {
        function_().hessian(global, result);
      }
     private:
      ExpressionStorage<FunctionType> function_;
    };
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
 
#endif // #ifndef __DUNE_ACFEM_FUNCTIONEXPRESSION_HH__