modeloperations.hh

#ifndef  __DUNE_ACFEM_MODELOPERATIONS_HH__
#define __DUNE_ACFEM_MODELOPERATIONS_HH__
 
#include "../expressions/expressionoperations.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    template<class Operation>
    struct ModelExpressionOperation
    {};
 
    template<class Operation>
    struct UnaryModelExpressionOperation
      : public ModelExpressionOperation<Operation>
    {};
 
    template<class Operation>
    struct BinaryModelExpressionOperation
      : public ModelExpressionOperation<Operation>
    {};
 
    //
    // Addition of Models
 
    template<>
    struct BinaryModelExpressionOperation<PlusOperation>
      : public ModelExpressionOperation<PlusOperation>
    {
      //flux
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void flux(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.flux(entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //linearizedFlux
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedFlux (const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //source
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType, class Entity>
      static void source(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasSources) {
          g.source(entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      //linearizedSource
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedSource(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasSources) {
          g.linearizedSource(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret += tmp;
        }
      }
 
      // fluxDivergence
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType, class HessianRangeType,
               class Entity>
      static void fluxDivergence(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.fluxDivergence(entity, x, value, jacobian, hessian, tmp);
          ret += tmp;
        }
      }
 
      // Robin b.c., non-linear
      template<class LeftModel, class RightModel,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void robinFlux(const LeftModel& lm, const RightModel& rm,
                            const Intersection& intersection,
                            const Point& x,
                            const DomainType& unitOuterNormal,
                            const RangeType& value,
                            RangeType& result)
      {
        RangeType tmp;
 
        lm.robinFlux(intersection, x, unitOuterNormal, value, result);
        if (ModelConstituents<RightModel>::hasRobinBoundary) {
          rm.robinFlux(intersection, x, unitOuterNormal, value, tmp);
          result += tmp;
        }
      }
 
      // Robin b.c., linearization
      template<class LeftModel, class RightModel,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void linearizedRobinFlux(const LeftModel& lm,
                                      const RightModel& rm,
                                      const RangeType& uBar,
                                      const Intersection& intersection,
                                      const Point& x,
                                      const DomainType& unitOuterNormal,
                                      const RangeType& value,
                                      RangeType& result)
      {
        RangeType tmp;
 
        lm.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        if (ModelConstituents<RightModel>::hasRobinBoundary) {
          rm.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, tmp);
          result += tmp;
        }
      }
    };
 
 
    //
    // Subtraction of Models
 
    template<>
    struct BinaryModelExpressionOperation<MinusOperation>
      : public ModelExpressionOperation<MinusOperation>
    {
      //flux
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void flux(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.flux(entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //linearizedFlux
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedFlux(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.linearizedFlux(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //source
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void source(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasSources) {
          g.source(entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //linearizedSource
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType,
               class Entity>
      static void linearizedSource(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasSources) {
          g.linearizedSource(uBar, DuBar, entity, x, value, jacobian, tmp);
          ret -= tmp;
        }
      }
 
      //fluxDivergence
      template<class LeftModel, class RightModel,
               class PointType, class RangeType, class JacobianRangeType, class HessianRangeType,
               class Entity>
      static void fluxDivergence(const LeftModel& f, const RightModel& 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 (ModelConstituents<RightModel>::hasFlux) {
          g.fluxDivergence(entity, x, value, jacobian, hessian, tmp);
          ret -= tmp;
        }
      }
 
      // Robin b.c., non-linear
      template<class LeftModel, class RightModel,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void robinFlux(const LeftModel& lm, const RightModel& rm,
                            const Intersection& intersection,
                            const Point& x,
                            const DomainType& unitOuterNormal,
                            const RangeType& value,
                            RangeType& result)
      {
        RangeType tmp;
 
        lm.robinFlux(intersection, x, unitOuterNormal, value, result);
        if (ModelConstituents<RightModel>::hasRobinBoundary) {
          // If left is also Zero then the result is already correct.
          rm.robinFlux(intersection, x, unitOuterNormal, value, tmp);
          result -= tmp;
        }
      }
 
      // Robin b.c., linearization
      template<class LeftModel, class RightModel,
               class Intersection, class Point,
               class DomainType, class RangeType>
      static void linearizedRobinFlux(const LeftModel& lm,
                                      const RightModel& rm,
                                      const RangeType& uBar,
                                      const Intersection& intersection,
                                      const Point& x,
                                      const DomainType& unitOuterNormal,
                                      const RangeType& value,
                                      RangeType& result)
      {
        RangeType tmp;
 
        lm.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, result);
        if (ModelConstituents<RightModel>::hasRobinBoundary) {
          rm.linearizedRobinFlux(uBar, intersection, x, unitOuterNormal, value, tmp);
          result -= tmp;
        }
      }
    };
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_MODELOPERATIONS_HH__