operatorparts.hh

#ifndef __DUNE_ACFEM_OPERATORPARTS_HH__
#define __DUNE_ACFEM_OPERATORPARTS_HH__
 
#include <dune/fem/misc/bartonnackmaninterface.hh>
 
namespace Dune {
 
  namespace ACFem {
 
    template<class PartsImpl>
    struct OperatorPartsTraits;
 
    template<class FunctionSpace>
    struct DefaultOperatorPartsTraits
    {
      typedef FunctionSpace FunctionSpaceType;
 
      typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
 
      enum StructureFlags {
        isLinear = false, 
        isSymmetric = false, 
        isSemiDefinite = false 
      };
 
      enum ConstituentFlags {
        hasFlux = true,      
        hasSources = true,   
        hasRobinFlux = true, 
      };
 
    };
 
    template<class PartsImpl>
    class OperatorPartsInterface
      : public Fem::BartonNackmanInterface<OperatorPartsInterface<PartsImpl>, PartsImpl>
    {
      typedef OperatorPartsInterface ThisType;
      typedef PartsImpl ImplementationType;
      typedef OperatorPartsTraits<ImplementationType> TraitsType;
      typedef Fem::BartonNackmanInterface<ThisType, ImplementationType> BaseType;
      using BaseType::asImp;
    protected:
      OperatorPartsInterface() {}
      OperatorPartsInterface(const OperatorPartsInterface& other) {}
      OperatorPartsInterface(OperatorPartsInterface&& other) {}
      OperatorPartsInterface& operator=(const OperatorPartsInterface& other) {}
    public:
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain,
             dimRange = FunctionSpaceType::dimRange };
 
      enum StructureFlags {
        isLinear = TraitsType::isLinear, 
        isSymmetric = TraitsType::isSymmetric, 
        isSemiDefinite = TraitsType::isSemiDefinite 
      };
 
      enum ConstituentFlags {
        hasFlux = TraitsType::hasFlux,
        hasSources = TraitsType::hasSources,
        hasRobinFlux = TraitsType::hasRobinFlux,
      };
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().setEntity(entity));
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        bool result;
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(result = asImp().setIntersection(intersection));
        return result;
      }
 
      std::string name() const
      {
        return asImp().name();
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType &value,
                const JacobianRangeType &jacobian,
                JacobianRangeType &flux) const
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().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
      {
        CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(asImp().fluxDivergence(entity, x, value, jacobian, hessian, result));
      }
    };
 
    template<class PartsImpl>
    class DefaultOperatorParts
      : public OperatorPartsInterface<PartsImpl>
    {
      typedef PartsImpl ImplementationType;
      typedef OperatorPartsInterface<ImplementationType> BaseType;
      typedef DefaultOperatorParts<ImplementationType> ThisType;
      typedef OperatorPartsTraits<ImplementationType> TraitsType;
      typedef DefaultOperatorPartsTraits<typename TraitsType::FunctionSpaceType>
      DefaultTraitsType;
    protected:
      DefaultOperatorParts() {}
    public:
      typedef typename DefaultTraitsType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      enum {
        dimDomain = FunctionSpaceType::dimDomain,
        dimRange = FunctionSpaceType::dimRange
      };
 
      enum StructureFlags {
        isLinear = TraitsType::isLinear, 
        isSymmetric = TraitsType::isSymmetric, 
        isSemiDefinite = TraitsType::isSemiDefinite 
      };
 
      enum ConstituentFlags {
        hasFlux = TraitsType::hasFlux,
        hasSources = TraitsType::hasSources,
        hasRobinFlux = TraitsType::hasRobinFlux,
      };
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        // nothing
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        // this is false as the default indicator classifies the
        // empty set.
        return hasRobinFlux;
      }
 
      std::string name() const
      {
        return "Model";
      }
 
      template<class Entity, class Point>
      void flux (const Entity& entity,
                 const Point &x,
                 const RangeType &value,
                 const JacobianRangeType &jacobian,
                 JacobianRangeType &flux) const
      {
        BaseType::linearizedFlux(value, jacobian, 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
      {
        flux = JacobianRangeType(0.);
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType &value,
                  const JacobianRangeType &jacobian,
                  RangeType &result) const
      {
        BaseType::linearizedSource(value, jacobian, 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
      {
        result = RangeType(0.);
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        BaseType::linearizedRobinFlux(value, 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
      {
        result = RangeType(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
      {
        result = RangeType(0.);
      }
 
    };
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
#endif  // __DUNE_ACFEM_OPERATORPARTS_HH__