operatorpartsoperations.hh

#ifndef  __DUNE_ACFEM_OPERATORPARTSOPERATIONS_HH__
#define __DUNE_ACFEM_OPERATORPARTSOPERATIONS_HH__
 
#include "../../expressions/expressionoperations.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    template<class Operation>
    struct OperatorPartsExpressionOperation
    {};
 
    template<class Operation>
    struct UnaryOperatorPartsExpressionOperation
      : public OperatorPartsExpressionOperation<Operation>
    {};
 
    template<class Operation>
    struct BinaryOperatorPartsExpressionOperation
      : public OperatorPartsExpressionOperation<Operation>
    {};
 
    //
    // Addition
 
    template<>
    struct BinaryOperatorPartsExpressionOperation<PlusOperation>
      : public OperatorPartsExpressionOperation<PlusOperation>
    {
      //flux
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void flux(const LeftOperand& f, const RightOperand& g,
                       const Entity &entity, const PointType& x,
                       const RangeType &value, const JacobianRangeType& jacobian,
                       JacobianRangeType& ret)
      {
        JacobianRangeType tmp;
 
        f.flux(entity, x, value, jacobian, ret);
        if (RightOperand::hasFlux) {
          g.flux(entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //linearizedFlux
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedFlux (const LeftOperand& f, const RightOperand& g,
                                  const RangeType& uBar, const JacobianRangeType& DuBar,
                                  const Entity &entity,const PointType& x,
                                  const RangeType  &value, const JacobianRangeType& jacobian,
                                  JacobianRangeType& ret)
      {
        JacobianRangeType tmp;
 
        f.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, ret);
        if (RightOperand::hasFlux) {
          g.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //source
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType, class Entity>
      static void source(const LeftOperand& f, const RightOperand& g,
                         const Entity &entity, const PointType& x,
                         const RangeType& value, const JacobianRangeType& jacobian,
                         RangeType& ret)
      {
        RangeType tmp;
 
        f.source(entity, x, value, jacobian, ret);
        if (RightOperand::hasSources) {
          g.source(entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //linearizedSource
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedSource(const LeftOperand& f, const RightOperand& g,
                                   const RangeType& uBar, const JacobianRangeType& DuBar,
                                   const Entity &entity,const PointType& x,
                                   const RangeType& value, const JacobianRangeType& jacobian,
                                   RangeType& ret)
      {
        RangeType tmp;
 
        f.linearizedSource(uBar, DuBar, entity, x, value, jacobian, ret);
        if (RightOperand::hasSources) {
          g.linearizedSource(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      // fluxDivergence
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType, class HessianRangeType,
               class Entity>
      static void fluxDivergence(const LeftOperand& f, const RightOperand& g,
                                 const Entity &entity, const PointType& x,
                                 const RangeType& value, const JacobianRangeType& jacobian,
                                 const HessianRangeType& hessian,
                                 RangeType& ret)
      {
        RangeType tmp;
 
        f.fluxDivergence(entity, x, value, jacobian, hessian, ret);
        if (RightOperand::hasFlux) {
          g.fluxDivergence(entity, x, value, jacobian, hessian, tmp);
          ret += tmp;
        }
      }
 
      // Robin b.c., non-linear
      template<class LeftOperand, class RightOperand,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void robinFlux(const LeftOperand& left, const RightOperand& right,
                            const Intersection& intersection,
                            const Point& x,
                            const DomainType& unitOuterNormal,
                            const RangeType& value,
                            RangeType& result)
      {
        RangeType tmp;
 
        left.robinFlux(intersection, x, unitOuterNormal, value, result);
        if (RightOperand::hasRobinFlux) {
          right.robinFlux(intersection, x, unitOuterNormal, value, tmp);
          result += tmp;
        }
      }
 
      // Robin b.c., linearization
      template<class LeftOperand, class RightOperand,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void linearizedRobinFlux(const LeftOperand& left,
                                      const RightOperand& right,
                                      const RangeType& uBar,
                                      const Intersection& intersection,
                                      const Point& x,
                                      const DomainType& unitOuterNormal,
                                      const RangeType& value,
                                      RangeType& result)
      {
        RangeType tmp;
 
        left.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        if (RightOperand::hasRobinFlux) {
          right.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, tmp);
          result += tmp;
        }
      }
    };
 
 
    //
    // Subtraction
 
    template<>
    struct BinaryOperatorPartsExpressionOperation<MinusOperation>
      : public OperatorPartsExpressionOperation<MinusOperation>
    {
      //flux
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void flux(const LeftOperand& f, const RightOperand& g,
                       const Entity &entity, const PointType& x,
                       const RangeType& value, const JacobianRangeType& jacobian,
                       JacobianRangeType& ret)
      {
        JacobianRangeType tmp;
 
        f.flux(entity, x, value, jacobian, ret);
        if (RightOperand::hasFlux) {
          g.flux(entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //linearizedFlux
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedFlux(const LeftOperand& f, const RightOperand& g,
                                 const RangeType& uBar, const JacobianRangeType& DuBar,
                                 const Entity &entity, const PointType& x,
                                 const RangeType  &value, const JacobianRangeType& jacobian,
                                 JacobianRangeType& ret)
      {
        JacobianRangeType tmp;
 
        f.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, ret);
        if (RightOperand::hasFlux) {
          g.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //source
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void source(const LeftOperand& f, const RightOperand& g,
                         const Entity &entity, const PointType& x,
                         const RangeType& value, const JacobianRangeType& jacobian,
                         RangeType& ret)
      {
        RangeType tmp;
 
        f.source(entity, x, value, jacobian, ret);
        if (RightOperand::hasSources) {
          g.source(entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //linearizedSource
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedSource(const LeftOperand& f, const RightOperand& g,
                                   const RangeType uBar, const JacobianRangeType DuBar,
                                   const Entity &entity, const PointType& x,
                                   const RangeType &value, const JacobianRangeType& jacobian,
                                   RangeType& ret)
      {
        RangeType tmp;
 
        f.linearizedSource(uBar, DuBar, entity, x, value, jacobian, ret);
        if (RightOperand::hasSources) {
          g.linearizedSource(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //fluxDivergence
      template<class LeftOperand, class RightOperand,
               class PointType, class RangeType, class JacobianRangeType, class HessianRangeType,
               class Entity>
      static void fluxDivergence(const LeftOperand& f, const RightOperand& g,
                                 const Entity &entity, const PointType& x,
                                 const RangeType &value, const JacobianRangeType& jacobian,
                                 const HessianRangeType& hessian,
                                 RangeType& ret)
      {
        RangeType tmp;
 
        f.fluxDivergence(entity, x, value, jacobian, hessian, ret);
        if (RightOperand::hasFlux) {
          g.fluxDivergence(entity, x, value, jacobian, hessian, tmp);
          ret -= tmp;
        }
      }
 
      // Robin b.c., non-linear
      template<class LeftOperand, class RightOperand,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void robinFlux(const LeftOperand& left, const RightOperand& right,
                            const Intersection& intersection,
                            const Point& x,
                            const DomainType& unitOuterNormal,
                            const RangeType& value,
                            RangeType& result)
      {
        RangeType tmp;
 
        left.robinFlux(intersection, x, unitOuterNormal, value, result);
        if (RightOperand::hasRobinFlux) {
          // If left is also Zero then the result is already correct.
          right.robinFlux(intersection, x, unitOuterNormal, value, tmp);
          result -= tmp;
        }
      }
 
      // Robin b.c., linearization
      template<class LeftOperand, class RightOperand,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void linearizedRobinFlux(const LeftOperand& left,
                                      const RightOperand& right,
                                      const RangeType& uBar,
                                      const Intersection& intersection,
                                      const Point& x,
                                      const DomainType& unitOuterNormal,
                                      const RangeType& value,
                                      RangeType& result)
      {
        RangeType tmp;
 
        left.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        if (RightOperand::hasRobinFlux) {
          right.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, tmp);
          result -= tmp;
        }
      }
    };
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_OPERATORPARTSOPERATIONS_HH__