robinmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_ROBINMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_ROBINMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace, class Indicator = EntireBoundaryIndicatorType>
    class RobinBoundaryOperatorParts
      : public OperatorPartsExpression<RobinBoundaryOperatorParts<FunctionSpace, Indicator> >
    {
      typedef RobinBoundaryOperatorParts ThisType;
      typedef OperatorPartsTraits<ThisType> TraitsType;
     public:
      typedef Indicator IndicatorType;
      typedef FunctionSpace FunctionSpaceType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
 
      RobinBoundaryOperatorParts(const BoundaryIndicatorInterface<IndicatorType>& indicator = IndicatorType(),
                                 const std::string& name = "")
        : indicator_(indicator),
          name_(name == "" ? "Robin(1)" : name),
          supported_(false)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        supported_ = isDomainBoundary(intersection) && indicator_().applies(intersection);
        return supported_;
      }
 
      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 = supported_ ? value : RangeType(0.); // it's that simple in this case ...
      }
 
     protected:
      ExpressionStorage<IndicatorType> indicator_;
      const std::string name_;
      mutable bool supported_;
    };
 
    template<class FunctionSpace, class Indicator>
    struct OperatorPartsTraits<RobinBoundaryOperatorParts<FunctionSpace, Indicator> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = true,
        isSymmetric = true,
        isSemiDefinite = true
      };
 
      enum ConstituentFlags {
        hasFlux = false,
        hasSources = false,
        hasRobinFlux = true,
      };
    };
 
    template<class FunctionSpace>
    class RobinBoundaryOperatorParts<FunctionSpace, EmptyBoundaryIndicatorType>
      : public ZeroOperatorParts<FunctionSpace>,
        public ExpressionTemplate<RobinBoundaryOperatorParts<FunctionSpace, EmptyBoundaryIndicatorType> >
    {
      typedef ZeroOperatorParts<FunctionSpace> BaseType;
     public:
      typedef FunctionSpace FunctionSpaceType;
      typedef EmptyBoundaryIndicatorType IndicatorType;
 
      RobinBoundaryOperatorParts(const BoundaryIndicatorInterface<IndicatorType>& indicator = IndicatorType(),
                                 const std::string& name = "")
        : BaseType(name == "" ? ("Robin(1)") : name)
      {}
    };
 
 
    template<class Object, class Indicator = EntireBoundaryIndicatorType>
    static inline
    RobinBoundaryOperatorParts<typename Object::FunctionSpaceType, Indicator>
    robinOperatorParts(const Object& object,
                       const BoundaryIndicatorInterface<Indicator>& where = Indicator(),
                       const std::string& name = "")
    {
      typedef RobinBoundaryOperatorParts<typename Object::FunctionSpaceType, Indicator> OperatorPartsType;
      return OperatorPartsType(where, name);
    }
 
    template<class Object, class Indicator = EntireBoundaryIndicatorType>
    OperatorPartsAdapterModel<RobinBoundaryOperatorParts<typename Object::FunctionSpaceType, Indicator>,
                              typename Object::GridPartType>
    robinZeroModel(const Object& object,
                   const BoundaryIndicatorInterface<Indicator>& where = Indicator())
    {
      typedef RobinBoundaryOperatorParts<typename Object::FunctionSpaceType, Indicator> OperatorPartsType;
      typedef typename Object::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(where));
    }
 
    template<class GridFunction, class Indicator = EntireBoundaryIndicatorType>
    auto
    robinBoundaryModel(const Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& values,
                       const BoundaryIndicatorInterface<Indicator>& where = Indicator())
      -> decltype(robinZeroModel(asImp(values),
                                 where * neumannBoundaryModel(asImp(values), where).neumannIndicator())
                  +
                  neumannBoundaryModel(asImp(values), where))
    {
      auto neumannModel = neumannBoundaryModel(values, where);
      auto robinIndicator = where * neumannModel.neumannIndicator();
 
      return robinZeroModel(asImp(values), robinIndicator) + neumannModel;
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_ROBINMODEL_HH__