operatorpartsexpression.hh

#ifndef __OPERATORPARTSEXPRESSION_HH__
#define __OPERATORPARTSEXPRESSION_HH__
 
#include "operatorparts.hh"
#include "operatorpartsexpressionbase.hh"
#include "operatorpartsoperations.hh"
#include "zerooperatorparts.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    template<class BinOp, class LeftOperand, class RightOperand>
    class BinaryOperatorPartsExpression;
 
    template<class UnOp, class Operand>
    class UnaryOperatorPartsExpression;
 
    template<class Operand>
    class UnaryOperatorPartsExpression<IdentityOperation, Operand>
      : public OperatorPartsExpression<UnaryOperatorPartsExpression<IdentityOperation, Operand> >
    {
      typedef Operand                                  OperandType;
      typedef IdentityOperation                        OpType;
      typedef UnaryOperatorPartsExpression<OpType, OperandType> ThisType;
      typedef OperatorPartsExpression<ThisType>         BaseType;
      typedef OperatorPartsInterface<ThisType>                   InterfaceType;
      typedef OperatorPartsTraits<ThisType>             TraitsType;
    public:
      typedef OperandType ContainedExpressionType;
 
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename InterfaceType::DomainFieldType DomainFieldType ;
      typedef typename InterfaceType::RangeFieldType RangeFieldType ;
      typedef typename InterfaceType::DomainType DomainType ;
      typedef typename InterfaceType::RangeType RangeType ;
      typedef typename InterfaceType::JacobianRangeType JacobianRangeType;
      typedef typename InterfaceType::HessianRangeType HessianRangeType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain, dimRange = FunctionSpaceType::dimRange };
 
    public:
 
      // constructor
      UnaryOperatorPartsExpression(const OperandType& op)
        : operand_(op)
      {}
 
      // copy constructor
      UnaryOperatorPartsExpression(const ThisType &other)
        : operand_( other.operand_ )
      {}
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        operand_().setEntity(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        return operand_().setIntersection(intersection); // recurse ...
      }
 
      std::string name() const
      {
        return OpType::name() + "(" + operand_().name() + ")";
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &flux ) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().flux(entity, x, value, jacobian, flux);
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          JacobianRangeType &flux) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().linearizedFlux(uBar, DuBar, entity, x, value, jacobian, flux);
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType &value,
                  const JacobianRangeType &jacobian,
                  RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().source(entity, x, value, jacobian, result);
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType &uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType &value,
                            const JacobianRangeType &jacobian,
                            RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().linearizedSource(uBar, DuBar, entity, x, value, jacobian, result);
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().robinFlux(intersection, x, unitOuterNormal, value, result);
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          const HessianRangeType &hessian,
                          RangeType &result) const
      {
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
        operand_().fluxDivergence(entity, x, value, jacobian, hessian, result);
      }
 
      const ContainedExpressionType& containedExpression() const
      {
        return operand_();
      }
 
    protected:
      ExpressionStorage<OperandType> operand_;
    };
 
    //traits
    template<class Operand>
    struct OperatorPartsTraits<UnaryOperatorPartsExpression<IdentityOperation, Operand> >
      : public DefaultOperatorPartsTraits<typename Operand::FunctionSpaceType>
    {
    protected:
      typedef Operand OperandType;
    public:
      enum StructureFlags {
        isLinear = OperandType::isLinear,
        isSymmetric = OperandType::isSymmetric,
        isSemiDefinite = OperandType::isSemiDefinite
      };
 
      enum ConstituentFlags {
        hasFlux = OperandType::hasFlux,
        hasSources = OperandType::hasSources,
        hasRobinFlux = OperandType::hasRobinFlux
      };
    };
 
    template<class Operand>
    class UnaryOperatorPartsExpression<MinusOperation, Operand>
      : public OperatorPartsExpression<UnaryOperatorPartsExpression<MinusOperation, Operand> >
    {
      typedef Operand                           OperandType;
      typedef MinusOperation                    OpType;
      typedef UnaryOperatorPartsExpression<OpType, OperandType> ThisType;
      typedef OperatorPartsExpression<ThisType> BaseType;
      typedef OperatorPartsInterface<ThisType>           InterfaceType;
      typedef OperatorPartsTraits<ThisType>     TraitsType;
    public:
      typedef OperandType ContainedExpressionType;
 
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename InterfaceType::DomainFieldType DomainFieldType ;
      typedef typename InterfaceType::RangeFieldType RangeFieldType ;
      typedef typename InterfaceType::DomainType DomainType ;
      typedef typename InterfaceType::RangeType RangeType ;
      typedef typename InterfaceType::JacobianRangeType JacobianRangeType;
      typedef typename InterfaceType::HessianRangeType HessianRangeType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain, dimRange = FunctionSpaceType::dimRange };
 
    public:
 
      // constructor
      UnaryOperatorPartsExpression(const OperandType& f)
        : operand_(f)
      {}
 
      // copy constructor
      UnaryOperatorPartsExpression(const ThisType &other)
        : operand_( other.operand_ )
      {}
 
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        operand_().setEntity(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        return operand_().setIntersection(intersection); // recurse ...
      }
 
      std::string name() const
      {
        return OpType::name() + "(" + operand_().name() + ")";
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &flux) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().flux(entity, x, value, jacobian, flux);
        flux *= -1.;
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          JacobianRangeType &flux ) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().linearizedFlux(uBar, DuBar, entity, x, value, jacobian, flux);
        flux *= -1.;
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType &value,
                  const JacobianRangeType &jacobian,
                  RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().source(entity, x, value, jacobian, result);
        result *= -1.;
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType &uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType &value,
                            const JacobianRangeType &jacobian,
                            RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().linearizedSource(uBar, DuBar, entity, x, value, jacobian, result);
        result *= -1.;
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().robinFlux(intersection, x, unitOuterNormal, value, result);
        result *= -1.;
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        result *= -1.0;
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          const HessianRangeType &hessian,
                          RangeType &result) const
      {
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
        operand_().fluxDivergence(entity, x, value, jacobian, hessian, result);
        result *= -1.;
      }
 
      const ContainedExpressionType& containedExpression() const
      {
        return operand_();
      }
 
    protected:
      ExpressionStorage<OperandType> operand_;
    };
 
    //traits
    template<class Operand>
    struct OperatorPartsTraits<UnaryOperatorPartsExpression<MinusOperation, Operand> >
      : public DefaultOperatorPartsTraits<typename Operand::FunctionSpaceType>
    {
    protected:
      typedef Operand OperandType;
    public:
      enum StructureFlags {
        isLinear = OperandType::isLinear,
        isSymmetric = OperandType::isSymmetric,
        isSemiDefinite = OperandType::isSemiDefinite
      };
 
      enum ConstituentFlags {
        hasFlux = OperandType::hasFlux,
        hasSources = OperandType::hasSources,
        hasRobinFlux = OperandType::hasRobinFlux
      };
 
      typedef typename OperandType::FunctionSpaceType FunctionSpaceType;
    };
 
    // BinaryOperatorPartsExpression
    // -------------------
 
    template<class BinOp, class LeftOperand, class RightOperand>
    class BinaryOperatorPartsExpression
      : public OperatorPartsExpression<BinaryOperatorPartsExpression<BinOp, LeftOperand, RightOperand> >
    {
      typedef BinOp                               BinOpType;
      typedef BinaryOperatorPartsExpression<BinOpType, LeftOperand, RightOperand> ThisType;
      typedef OperatorPartsExpression<ThisType>   BaseType;
      typedef OperatorPartsInterface<ThisType>             InterfaceType;
      typedef BinaryOperatorPartsExpressionOperation<BinOpType> OperationType;
      typedef OperatorPartsTraits<ThisType>       TraitsType;
    public:
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename InterfaceType::DomainFieldType DomainFieldType ;
      typedef typename InterfaceType::RangeFieldType RangeFieldType ;
      typedef typename InterfaceType::DomainType DomainType ;
      typedef typename InterfaceType::RangeType RangeType ;
      typedef typename InterfaceType::JacobianRangeType JacobianRangeType;
      typedef typename InterfaceType::HessianRangeType HessianRangeType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain, dimRange = FunctionSpaceType::dimRange };
 
    public:
 
      // constructor
      BinaryOperatorPartsExpression(const LeftOperand &f,
                                    const RightOperand &g)
        : leftOperand_( f ), rightOperand_( g )
      {}
 
      // copy constructor
      BinaryOperatorPartsExpression(const ThisType &other )
        : leftOperand_( other.leftOperand_ ),
          rightOperand_( other.rightOperand_ )
      {}
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        leftOperand_().setEntity(entity);
        rightOperand_().setEntity(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        // This corresponds to the combined Robin-indicator ...  Do
        // not wrap the calls to setIntersection() into the return
        // value, as then the second one is not executed if the first
        // one yields true.
        bool leftRobin  = leftOperand_().setIntersection(intersection);
        bool rightRobin = rightOperand_().setIntersection(intersection);
 
        return leftRobin || rightRobin;
      }
 
      std::string name() const
      {
        return "(" + leftOperand_().name() + " " + BinOpType::name() + " " + rightOperand_().name() + ")";
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &flux ) const
      {
        OperationType::flux(leftOperand_(), rightOperand_(), entity, x, value, jacobian, flux);
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          JacobianRangeType &flux ) const
      {
        OperationType::linearizedFlux(leftOperand_(), rightOperand_(), uBar, DuBar, entity, x, value, jacobian, flux);
      }
 
      template<class Entity, class Point>
      void source( const Entity& entity,
                   const Point &x,
                   const RangeType &value,
                   const JacobianRangeType &jacobian,
                   RangeType &result) const
      {
        OperationType::source(leftOperand_(), rightOperand_(), entity, x, value, jacobian, result);
      }
 
      template<class Entity, class Point>
      void linearizedSource( const RangeType &uBar,
                             const JacobianRangeType& DuBar,
                             const Entity& entity,
                             const Point &x,
                             const RangeType &value,
                             const JacobianRangeType &jacobian,
                             RangeType &result) const
      {
        OperationType::linearizedSource(leftOperand_(), rightOperand_(), uBar, DuBar, entity, x, value, jacobian, result);
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        OperationType::robinFlux(leftOperand_(), rightOperand_(), intersection, x, unitOuterNormal, value, result);
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        OperationType::linearizedRobinFlux(leftOperand_(), rightOperand_(), uBar, intersection, x, unitOuterNormal, value, result);
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          const HessianRangeType &hessian,
                          RangeType &result) const
      {
        OperationType::fluxDivergence(leftOperand_(), rightOperand_(), entity, x, value, jacobian, hessian, result);
      }
 
    protected:
      ExpressionStorage<LeftOperand>  leftOperand_;
      ExpressionStorage<RightOperand> rightOperand_;
    };
 
    template<class BinOp, class LeftOperand, class RightOperand>
    struct OperatorPartsTraits<BinaryOperatorPartsExpression<BinOp, LeftOperand, RightOperand> >
      : public DefaultOperatorPartsTraits<typename LeftOperand::FunctionSpaceType>
    {
    private:
      typedef LeftOperand LeftOperandType;
      typedef RightOperand RightOperandType;
      typedef typename LeftOperand::FunctionSpaceType LeftFunctionSpaceType;
      typedef typename RightOperandType::FunctionSpaceType RightFunctionSpaceType;
 
      typedef BinaryOperatorPartsExpression<BinOp, LeftOperandType, RightOperandType> OperandType;
 
      static_assert((std::is_same<LeftFunctionSpaceType, RightFunctionSpaceType>::value),
                    "Function spaces must agree, cannot combine them");
    public:
      enum StructureFlags {
        isLinear = LeftOperandType::isLinear && RightOperandType::isLinear,
        isSymmetric = LeftOperandType::isSymmetric && RightOperandType::isSymmetric,
        isSemiDefinite = LeftOperandType::isSemiDefinite && RightOperandType::isSemiDefinite // this is not true
      };
 
      enum ConstituentFlags {
        hasFlux = LeftOperandType::hasFlux || RightOperandType::hasFlux,
        hasSources = LeftOperandType::hasSources || RightOperandType::hasSources,
        hasRobinFlux = LeftOperandType::hasRobinFlux || RightOperandType::hasRobinFlux,
      };
 
      typedef LeftFunctionSpaceType FunctionSpaceType;
    };
 
    //
    // Addition
 
    template<class LeftOperand, class RightOperand>
    BinaryOperatorPartsExpression<PlusOperation, LeftOperand, RightOperand>
    operator+(const OperatorPartsInterface<LeftOperand>& f_,
              const OperatorPartsInterface<RightOperand>& g_)
    {
      const LeftOperand& f(static_cast<const LeftOperand&>(f_));
      const RightOperand& g(static_cast<const RightOperand&>(g_));
 
      typedef BinaryOperatorPartsExpression<PlusOperation, LeftOperand, RightOperand> ExpressionType;
 
      return ExpressionType(f, g);
    }
 
 
    //
    // Subtraction
 
    template<class LeftOperand, class RightOperand>
    BinaryOperatorPartsExpression<MinusOperation, LeftOperand, RightOperand>
    operator-(const OperatorPartsInterface<LeftOperand>& f_,
              const OperatorPartsInterface<RightOperand>& g_)
    {
      const LeftOperand& f(static_cast<const LeftOperand&>(f_));
      const RightOperand& g(static_cast<const RightOperand&>(g_));
 
      typedef BinaryOperatorPartsExpression<MinusOperation, LeftOperand, RightOperand> ExpressionType;
 
      return ExpressionType(f, g);
    }
 
    template<class Factor, class OperandType>
    class BinaryOperatorPartsExpression<SMultiplyOperation, Factor, OperandType>
      : public OperatorPartsExpression<BinaryOperatorPartsExpression<SMultiplyOperation, Factor, OperandType> >
    {
      typedef SMultiplyOperation                BinOpType;
      typedef BinaryOperatorPartsExpression<BinOpType, Factor, OperandType> ThisType;
      typedef OperatorPartsExpression<ThisType> BaseType;
      typedef OperatorPartsInterface<ThisType>           InterfaceType;
      typedef OperatorPartsTraits<ThisType>     TraitsType;
    public:
      typedef Factor FactorType;
 
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename InterfaceType::DomainFieldType DomainFieldType ;
      typedef typename InterfaceType::RangeFieldType RangeFieldType ;
      typedef typename InterfaceType::DomainType DomainType ;
      typedef typename InterfaceType::RangeType RangeType ;
      typedef typename InterfaceType::JacobianRangeType JacobianRangeType;
      typedef typename InterfaceType::HessianRangeType HessianRangeType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain, dimRange = FunctionSpaceType::dimRange };
 
    public:
 
      // constructor
      BinaryOperatorPartsExpression(const FactorType& s, const OperandType& f)
        : scalar_(s), operand_(f)
      {}
 
      // copy constructor
      BinaryOperatorPartsExpression(const ThisType &other)
        : scalar_( other.scalar_ ), operand_( other.operand_ )
      {}
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        operand_().setEntity(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        return operand_().setIntersection(intersection); // recurse ...
      }
 
      std::string name() const
      {
        return "(" + std::to_string(parameterValue(scalar_())) + " " + BinOpType::name() + " " + operand_().name() + ")";
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &flux) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().flux(entity, x, value, jacobian, flux);
        flux *= parameterValue(scalar_());
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          JacobianRangeType &flux) const
      {
        if (!OperandType::hasFlux) {
          flux = 0.;
          return;
        }
        operand_().linearizedFlux(uBar, DuBar, entity, x, value, jacobian, flux);
        flux *= parameterValue(scalar_());
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType &value,
                  const JacobianRangeType &jacobian,
                  RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().source(entity, x, value, jacobian, result);
        result *= parameterValue(scalar_());
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType &uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType &value,
                            const JacobianRangeType &jacobian,
                            RangeType &result) const
      {
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
        operand_().linearizedSource(uBar, DuBar, entity, x, value, jacobian, result);
        result *= parameterValue(scalar_());
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().robinFlux(intersection, x, unitOuterNormal, value, result);
        result *= parameterValue(scalar_());
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
        operand_().linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        result *= parameterValue(scalar_());
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          const HessianRangeType &hessian,
                          RangeType &result) const
      {
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
        operand_().fluxDivergence(entity, x, value, jacobian, hessian, result);
        result *= parameterValue(scalar_());
      }
 
    public:
      const FactorType& scalar() const { return scalar_(); }
      const OperandType& operand() const { return operand_(); }
    protected:
      ExpressionStorage<FactorType> scalar_;
      ExpressionStorage<OperandType> operand_;
    };
 
    //traits
    template<class Factor, class Operand>
    struct OperatorPartsTraits<BinaryOperatorPartsExpression<SMultiplyOperation, Factor, Operand> >
      : public DefaultOperatorPartsTraits<typename Operand::FunctionSpaceType>
    {
    protected:
      typedef Operand OperandType;
      typedef Factor FactorType;
      typedef typename OperandType::FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
    public:
      enum StructureFlags {
        isLinear = OperandType::isLinear,
        isSymmetric = OperandType::isSymmetric,
        isSemiDefinite = OperandType::isSemiDefinite // this is not true when multiplying by negative factors ...
      };
 
      enum ConstituentFlags {
        hasFlux = OperandType::hasFlux,
        hasSources = OperandType::hasSources
      };
 
      typedef SMultiplyOperation BinOpType;
      typedef typename OperandType::FunctionSpaceType FunctionSpaceType;
    };
 
 
    template<class OperandType>
    BinaryOperatorPartsExpression<SMultiplyOperation, typename OperandType::RangeFieldType, OperandType>
    operator*(const typename OperandType::RangeFieldType& s_,
              const OperatorPartsInterface<OperandType>& f_)
    {
      typedef typename OperandType::RangeFieldType RangeFieldType;
      const OperandType& f(static_cast<const OperandType&>(f_));
 
      typedef BinaryOperatorPartsExpression<SMultiplyOperation, RangeFieldType, OperandType> ExpressionType;
 
      return ExpressionType(s_, f);
    }
 
    template<class OperandType>
    auto operator*(const OperatorPartsInterface<OperandType>& op_,
                   const typename OperandType::RangeFieldType& s_)
      -> decltype(s_ * asImp(op_))
    {
      return s_ * asImp(op_);
    }
 
    template<class Parameter, class OperandType>
    BinaryOperatorPartsExpression<SMultiplyOperation, Parameter, OperandType>
    operator*(const ParameterInterface<Parameter>& s_,
              const OperatorPartsInterface<OperandType>& f_)
    {
      typedef Parameter ParameterType;
 
      const ParameterType& s(static_cast<const ParameterType&>(s_));
      const OperandType& f(static_cast<const OperandType&>(f_));
 
      typedef BinaryOperatorPartsExpression<SMultiplyOperation, ParameterType, OperandType> ExpressionType;
 
      return ExpressionType(s, f);
    }
 
    template<class Parameter, class OperandType>
    auto operator*(const OperatorPartsInterface<OperandType>& op_,
                   const ParameterInterface<Parameter>& s_)
      -> decltype(asImp(s_) * asImp(op_))
    {
      return asImp(s_) * asImp(op_);
    }
 
#if 0
    // We have real unary minus by now.
    // unary minus operator
    template<class OperandType>
    BinaryOperatorPartsExpression<SMultiplyOperation, typename OperandType::RangeFieldType, OperandType>
    operator*(const OperatorPartsInterface<OperandType>& f_)
    {
      typedef typename OperandType::RangeFieldType RangeFieldType;
      const OperandType& f(static_cast<const OperandType&>(f_));
 
      typedef BinaryOperatorPartsExpression<SMultiplyOperation, RangeFieldType, OperandType> ExpressionType;
 
      return ExpressionType(-1.0, f);
    }
#endif
 
    template<class FactorFunction, class OperandType>
    class BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, OperandType>
      : public OperatorPartsExpression<BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, OperandType> >
    {
      typedef SMultiplyOperation                     BinOpType;
      typedef BinaryOperatorPartsExpression          ThisType;
      typedef OperatorPartsExpression<ThisType>      BaseType;
      typedef OperatorPartsInterface<ThisType>                InterfaceType;
      typedef OperatorPartsTraits<ThisType>          TraitsType;
      typedef FactorFunction                         FactorFunctionType;
      typedef typename FactorFunctionType::LocalFunctionType FactorLocalFunctionType;
      typedef typename FactorFunctionType::RangeType         FactorRangeType;
      typedef typename FactorFunctionType::JacobianRangeType FactorJacobianRangeType;
 
    public:
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename InterfaceType::DomainFieldType DomainFieldType ;
      typedef typename InterfaceType::RangeFieldType RangeFieldType ;
      typedef typename InterfaceType::DomainType DomainType ;
      typedef typename InterfaceType::RangeType RangeType ;
      typedef typename InterfaceType::JacobianRangeType JacobianRangeType;
      typedef typename InterfaceType::HessianRangeType HessianRangeType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain, dimRange = FunctionSpaceType::dimRange };
      enum { isLinear = TraitsType::isLinear,
             isSymmetric = TraitsType::isSymmetric,
             isSemiDefinite = TraitsType::isSemiDefinite
      };
 
    public:
 
      // constructor
      BinaryOperatorPartsExpression(const FactorFunctionType& f, const OperandType& op)
        : factor_(f), operand_(op), localFactor_(factor_())
      {}
 
      // copy constructor
      BinaryOperatorPartsExpression(const ThisType &other)
        : factor_(other.factor_), operand_(other.operand_), localFactor_(factor_())
      {}
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        operand_().setEntity(entity);
        localFactor_.init(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        // Should we also initialize localFactor_ here? For the moment
        // we assume that any high/level code first computes the
        // bulk-contributions, so localFactor_ should alreay be
        // initialized here
        //
        // Would the following work? Probably not
        // assert(&localFactor_.entity() == &*intersection.inside());
        return operand_().setIntersection(intersection); // recurse ...
      }
 
      std::string name() const
      {
        return "(" + factor_().name() + " " + BinOpType::name() + " " + operand_().name() + ")";
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &result) const
      {
        assert(&entity == &localFactor_.entity());
 
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
        operand_().flux(entity, x, value, jacobian, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          JacobianRangeType &result) const
      {
        assert(&entity == &localFactor_.entity());
 
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
 
        operand_().linearizedFlux(uBar, DuBar, entity, x, value, jacobian, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType &value,
                  const JacobianRangeType &jacobian,
                  RangeType &result) const
      {
        assert(&entity == &localFactor_.entity());
 
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
 
        operand_().source(entity, x, value, jacobian, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType &uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType &value,
                            const JacobianRangeType &jacobian,
                            RangeType &result) const
      {
        assert(&entity == &localFactor_.entity());
 
        if (!OperandType::hasSources) {
          result = 0.;
          return;
        }
 
        operand_().linearizedSource(uBar, DuBar, entity, x, value, jacobian, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
 
        operand_().robinFlux(intersection, x, unitOuterNormal, value, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        if (!OperandType::hasRobinFlux) {
          result = 0.;
          return;
        }
 
        operand_().linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        result *= factorValue;
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType &value,
                          const JacobianRangeType &jacobian,
                          const HessianRangeType &hessian,
                          RangeType &result) const
      {
        if (!OperandType::hasFlux) {
          result = 0.;
          return;
        }
 
        // Product rule
 
        operand_().fluxDivergence(entity, x, value, jacobian, hessian, result);
 
        FactorRangeType factorValue;
        localFactor_.evaluate(x, factorValue);
 
        JacobianRangeType fluxValue;
        FactorJacobianRangeType factorJacobian;
 
        operand_().flux(entity, x, value, jacobian, fluxValue);
        localFactor_.jacobian(x, factorJacobian);
 
        for (int i = 0; i < dimRange; ++i) {
          result += factorJacobian[i] * fluxValue[i];
        }
      }
 
    protected:
      ExpressionStorage<FactorFunctionType> factor_;
      ExpressionStorage<OperandType> operand_;
      mutable FactorLocalFunctionType localFactor_;
    };
 
    //traits
    template<class FactorFunction, class Operand>
    struct OperatorPartsTraits<BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, Operand> >
      : public DefaultOperatorPartsTraits<typename Operand::FunctionSpaceType>
    {
    protected:
      typedef MultiplyOperation BinOpType;
      typedef FactorFunction FactorFunctionType;
      typedef Operand FactorOperandType;
 
    public:
 
      enum StructureFlags
        {
          isLinear = FactorOperandType::isLinear,
          isSymmetric = FactorOperandType::isSymmetric,
          isSemiDefinite = FactorOperandType::isSemiDefinite // this is not necessarily true
        };
 
      enum ConstituentFlags {
        hasFlux = FactorOperandType::hasFlux,
        hasSources = FactorOperandType::hasSources,
        hasRobinFlux = FactorOperandType::hasRobinFlux
      };
 
      typedef typename FactorOperandType::FunctionSpaceType FunctionSpaceType;
 
    };
 
    template<class FactorFunction, class OperandType>
    BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, OperandType>
    operator*(const Fem::Function<typename FactorFunction::FunctionSpaceType::ScalarFunctionSpaceType, FactorFunction>& f_,
              const OperatorPartsInterface<OperandType>& op_)
    {
      const FactorFunction& f(static_cast<const FactorFunction&>(f_));
      const OperandType& op(static_cast<const OperandType&>(op_));
 
      typedef BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, OperandType> ExpressionType;
 
      return ExpressionType(f, op);
    }
 
    template<class FactorFunction, class OperandType>
    BinaryOperatorPartsExpression<MultiplyOperation, FactorFunction, OperandType>
    operator*(const OperatorPartsInterface<OperandType>& op_,
              const Fem::Function<typename FactorFunction::FunctionSpaceType::ScalarFunctionSpaceType, FactorFunction>& f_)
    {
      const FactorFunction& f(static_cast<const FactorFunction&>(f_));
      const OperandType& op(static_cast<const OperandType&>(op_));
 
      return f * op;
    }
 
    template<class Operand>
    class UnaryOperatorPartsExpression<MinusOperation, Operand>
    operator-(const OperatorPartsInterface<Operand>& op_)
    {
      const Operand& op(static_cast<const Operand&>(op_));
 
      typedef UnaryOperatorPartsExpression<MinusOperation, Operand> ExpressionType;
 
      return ExpressionType(op);
    }
 
    template<class Operand>
    class UnaryOperatorPartsExpression<IdentityOperation, Operand>
    operator*(const OperatorPartsInterface<Operand>& op_)
    {
      const Operand& op(static_cast<const Operand&>(op_));
 
      typedef UnaryOperatorPartsExpression<IdentityOperation, Operand> ExpressionType;
 
      return ExpressionType(op);
    }
 
    template<class Operand>
    auto
    operator+(const OperatorPartsInterface<Operand>& op_)
      -> decltype(*op_)
    {
      return *op_;
    }
 
    template<class Expression>
    Expression
    operator*(const OperatorPartsExpression<Expression>& op_)
    {
      return op_.expression();
    }
 
    template<class Operand>
    static inline
    auto operator*(const typename Operand::RangeFieldType& s_,
                   const BinaryOperatorPartsExpression<SMultiplyOperation, typename Operand::RangeFieldType, Operand>& op_)
      -> decltype((s_ * op_.scalar()) * op_.operand())
    {
      return (s_ * op_.scalar()) * op_.operand();
    }
 
#if 0
    template<class Parameter, class Field, class Operand>
    static inline
    auto operator*(const ParameterInterface<Parameter>& p_,
                   const BinaryOperatorPartsExpression<SMultiplyOperation, Field, Operand>& op_)
      -> decltype(op_.scalar() * (asImp(p_) * op_.operand()))
    {
      return op_.scalar() * (asImp(p_) * op_.operand());
    }
#else
    template<class Parameter, class Operand>
    static inline
    auto operator*(const ParameterInterface<Parameter>& p_,
                   const OperatorPartsInterface<BinaryOperatorPartsExpression<SMultiplyOperation, typename Operand::RangeFieldType, Operand> >& op_)
      -> decltype(asImp(op_).scalar() * (asImp(p_) * asImp(op_).operand()))
    {
      return asImp(op_).scalar() * (asImp(p_) * asImp(op_).operand());
    }
#endif
 
    template<class Field, class Parameter, class Operand>
    static inline
    auto operator*(const BinaryParameterExpression<SMultiplyOperation, Field, Parameter>& p_,
                   const OperatorPartsInterface<Operand>& op_)
      -> decltype(p_.left() * (p_.right() * asImp(op_)))
    {
      return p_.left() * (p_.right() * asImp(op_));
    }
 
    template<class Function, class Field, class Operand>
    static inline
    auto operator*(const Fem::Function<typename Function::FunctionSpaceType, Function>& f_,
                   const BinaryOperatorPartsExpression<SMultiplyOperation, Field, Operand>& op_)
      -> decltype(op_.scalar() * (asImp(f_) * op_.operand()))
    {
      return op_.scalar() * (asImp(f_) * op_.operand());
    }
 
    template<class Field, class Function, class Operand>
    static inline
    auto operator*(const BinaryGridFunctionExpression<SMultiplyOperation, Field, Function>& f_,
                   const OperatorPartsInterface<Operand>& op_)
      -> decltype(f_.scalar() * (f_.function() * asImp(op_)))
    {
      return f_.scalar() * (f_.function() * asImp(op_));
    }
 
    template<class Operand>
    auto
    operator-(const UnaryOperatorPartsExpression<MinusOperation, Operand>& op_)
      -> decltype(*op_.containedExpression())
    {
      return *op_.containedExpression();
    }
 
    template<class Operand>
    auto
    operator+(const OperatorPartsInterface<Operand>& op,
              const ZeroOperatorParts<typename Operand::FunctionSpaceType>& z)
      -> decltype(*asImp(op))
    {
      return *asImp(op);
    }
 
    template<class Operand>
    auto
    operator+(const ZeroOperatorParts<typename Operand::FunctionSpaceType>& z,
              const OperatorPartsInterface<Operand>& op)
      -> decltype(*asImp(op))
    {
      return *asImp(op);
    }
 
    template<class FunctionSpace>
    auto
    operator+(const ZeroOperatorParts<FunctionSpace>& z1,
              const ZeroOperatorParts<FunctionSpace>& z2)
      -> decltype(*z1)
    {
      return *z1;
    }
 
    template<class Operand>
    auto
    operator-(const OperatorPartsInterface<Operand>& op,
              const ZeroOperatorParts<typename Operand::FunctionSpaceType>& z)
      -> decltype(*asImp(op))
    {
      return *asImp(op);
    }
 
    template<class Operand>
    auto
    operator-(const ZeroOperatorParts<typename Operand::FunctionSpaceType>& z,
              const OperatorPartsInterface<Operand>& op)
      -> decltype(-op)
    {
      return -op;
    }
 
    template<class FunctionSpace>
    auto
    operator-(const ZeroOperatorParts<FunctionSpace>& z1,
              const ZeroOperatorParts<FunctionSpace>& z2)
      -> decltype(*z1)
    {
      return *z1;
    }
 
    template<class FunctionSpace>
    auto
    operator*(const typename FunctionSpace::RangeFieldType& s_,
              const ZeroOperatorParts<FunctionSpace>& z)
      ->decltype(*z)
    {
      return *z;
    }
 
    template<class FunctionSpace>
    auto
    operator*(const ZeroOperatorParts<FunctionSpace>& z,
              const typename FunctionSpace::RangeFieldType& s)
      ->decltype(*z)
    {
      return *z;
    }
 
    template<class Parameter, class FunctionSpace>
    auto
    operator*(const ParameterInterface<Parameter>& s,
              const ZeroOperatorParts<FunctionSpace>& z)
      ->decltype(*z)
    {
      return *z;
    }
 
    template<class Parameter, class FunctionSpace>
    auto
    operator*(const ZeroOperatorParts<FunctionSpace>& z,
              const ParameterInterface<Parameter>& s)
      ->decltype(*z)
    {
      return *z;
    }
 
    template<class Function, class ScalarSpace, class FunctionSpace>
    auto
    operator*(const Fem::Function<ScalarSpace, Function>& f,
              const ZeroOperatorParts<FunctionSpace>& z)
      ->decltype(*z)
    {
      static_assert(std::is_same<ScalarSpace, typename FunctionSpace::ScalarFunctionSpaceType>::value,
                    "Incompatible function spaces.");
      return *z;
    }
 
    template<class Function, class ScalarSpace, class FunctionSpace>
    auto
    operator*(const ZeroOperatorParts<FunctionSpace>& z,
              const Fem::Function<ScalarSpace, Function>& f)
      ->decltype(*z)
    {
      static_assert(std::is_same<ScalarSpace, typename FunctionSpace::ScalarFunctionSpaceType>::value,
                    "Incompatible function spaces.");
      return *z;
    }
 
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __OPERATORPARTSEXPRESSION_HH__