boundarysupportedfunction.hh

#ifndef __DUNE_ACFEM_BOUNDARY_SUPPORTED_FUNCTION_HH__
#define __DUNE_ACFEM_BOUNDARY_SUPPORTED_FUNCTION_HH__
 
#include <dune/fem/function/common/function.hh>
 
#include "../models/boundaryindicator.hh"
#include "../common/stringhelper.hh"
#include "../functions/gridfunctionexpressionbase.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    struct HasBoundarySupport
    {};
 
    struct HasNoBoundarySupport
    {};
 
    // forward
    template<class GridFunction, class Indicator = EntireBoundaryIndicatorType>
    class BoundarySupportedFunction;
 
    template<class F, bool bndrySupport = std::is_base_of<HasBoundarySupport, F>::value>
    struct BoundaryFunctionTraits;
 
    template<class F>
    struct BoundaryFunctionTraits<F, true>
    {
      enum {
        hasBoundarySupport = true 
      };
      typedef F FunctionType;
      typedef typename F::FunctionSpaceType FunctionSpaceType;
      typedef typename F::IndicatorType IndicatorType;
      typedef typename F::LocalFunctionType LocalFunctionType;
      typedef typename F::GridPartType GridPartType;
      typedef typename F::EntityType EntityType;
      typedef typename GridPartType::IntersectionType IntersectionType;
 
      enum {
        emptySupport = ExpressionTraits<IndicatorType>::isZero,
        globalSupport = std::is_same<IndicatorType, EntireBoundaryIndicatorType>::value
      };
 
      template<class OtherIndicator>
      struct ProductIndicator
      {
        typedef
        typename std::decay<decltype(*std::declval<IndicatorType>()
                                     &&
                                     *std::declval<OtherIndicator>())>::type
        Type;
 
        // form zero expressions, when appropriate
        enum {
          emptySupport = ExpressionTraits<Type>::isZero,
          globalSupport = std::is_same<Type, EntireBoundaryIndicatorType>::value
        };
 
        typedef BoundarySupportedFunction<FunctionType, OtherIndicator> BoundaryFunctionType;
        typedef typename
        std::conditional<emptySupport,
                         GridFunctionExpression<FunctionSpaceType, BoundaryFunctionType, ZeroExpression>,
                         GridFunctionExpression<FunctionSpaceType, BoundaryFunctionType, ExpressionTemplate>
                         >::type
        ExpressionBaseType;
      };
 
      static void init(LocalFunctionType& lf,
                       const EntityType& entity,
                       const IntersectionType& intersection)
      {
        lf.init(entity, intersection);
      }
 
      static bool supported(const LocalFunctionType& lf)
      {
        return globalSupport || (!emptySupport && lf.supported());
      }
 
      static typename F::IndicatorType indicator(const F& f)
      {
          return f.indicator();
      }
    };
 
    template<class F>
    struct BoundaryFunctionTraits<F, false>
    {
      enum {
        hasBoundarySupport = false 
      };
      typedef F FunctionType;
      typedef typename F::FunctionSpaceType FunctionSpaceType;
      typedef EntireBoundaryIndicatorType IndicatorType;
      typedef typename F::LocalFunctionType LocalFunctionType;
      typedef typename F::GridPartType GridPartType;
 
      typedef typename F::EntityType EntityType;
 
      typedef typename GridPartType::IntersectionType IntersectionType;
 
      enum {
        emptySupport = false, 
        globalSupport = true  
      };
 
      template<class OtherIndicator>
      struct ProductIndicator
      {
        typedef OtherIndicator Type;
 
        // Convenience
        enum {
          emptySupport = std::is_same<Type, EmptyBoundaryIndicatorType>::value,
          globalSupport = std::is_same<Type, EntireBoundaryIndicatorType>::value
        };
        typedef BoundarySupportedFunction<FunctionType, OtherIndicator> BoundaryFunctionType;
        typedef typename
        std::conditional<emptySupport,
                         GridFunctionExpression<FunctionSpaceType, BoundaryFunctionType, ZeroExpression>,
                         GridFunctionExpression<FunctionSpaceType, BoundaryFunctionType, ExpressionTemplate>
                         >::type
        ExpressionBaseType;
      };
 
      static void init(LocalFunctionType& lf,
                       const EntityType& entity,
                       const IntersectionType& intersection)
      {
        lf.init(entity);
      }
 
      static bool supported(const LocalFunctionType& lf)
      {
        return true;
      }
 
      static IndicatorType indicator(const F& ignored)
      {
        return EntireBoundaryIndicatorType();
      }
    };
 
 
    template<class GridFunction, class Indicator>
    class BoundarySupportedFunction
      : public BoundaryFunctionTraits<GridFunction>::template ProductIndicator<Indicator>::ExpressionBaseType,
        public HasBoundarySupport
    {
      typedef BoundaryFunctionTraits<GridFunction> BoundaryTraitsType;
      typedef BoundarySupportedFunction ThisType;
      typedef typename
      BoundaryTraitsType::template ProductIndicator<Indicator>::ExpressionBaseType
      ExpressionBaseType;
      class LocalFunction;
      using ExpressionBaseType::expressionName_;
 
      typedef typename BoundaryTraitsType::IndicatorType InnerIndicatorType;
     public:
      typedef LocalFunction LocalFunctionType;
 
      typedef
      typename std::decay<decltype(*std::declval<Indicator>()
                                   &&
                                   *std::declval<InnerIndicatorType>())>::type
      IndicatorType;
      typedef Indicator OuterIndicatorType;
 
      enum {
        emptySupport = std::is_same<IndicatorType, EmptyBoundaryIndicatorType>::value,
        globalSupport = std::is_same<IndicatorType, EntireBoundaryIndicatorType>::value
      };
 
      typedef GridFunction GridFunctionType;
 
      typedef ThisType DiscreteFunctionType;
      typedef typename GridFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename GridFunctionType::FunctionSpaceType FunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename GridFunctionType::EntityType EntityType;
      typedef typename GridPartType::IntersectionType IntersectionType;
 
      // Do not define a default argument for indicator, otherwise
      // this could result in automatic type-casts from any function into this type.
      BoundarySupportedFunction(const Fem::Function<FunctionSpaceType, GridFunctionType>& function,
                                const BoundaryIndicatorInterface<Indicator>& indicator)
        : gridFunction_(function),
          outerIndicator_(indicator),
          indicator_(*asImp(outerIndicator()) && *BoundaryTraitsType::indicator(gridFunction_()))
      {
        expressionName_ = makeName(gridFunction_().name());
      }
 
      void evaluate(const DomainType &global, RangeType &result) const
      {
        if (emptySupport) {
          result = 0.;
          return;
        }
        function().evaluate(global, result);
      }
 
      void jacobian(const DomainType &global, JacobianRangeType &result) const
      {
        if (emptySupport) {
          result = 0.;
          return;
        }
        function().jacobian(global, result);
      }
 
      void hessian(const DomainType &global, HessianRangeType &result) const
      {
        if (emptySupport) {
          result = 0.;
          return;
        }
        function().hessian(global, result);
      }
 
      const LocalFunctionType localFunction(const EntityType &entity) const
      {
        return LocalFunctionType(entity, *this);
      }
 
      LocalFunctionType localFunction(const EntityType &entity)
      {
        return LocalFunctionType(entity, *this);
      }
 
      const LocalFunctionType localFunction(const EntityType& entity,
                                            const IntersectionType &intersection) const
      {
        return LocalFunctionType(entity, intersection, *this);
      }
 
      LocalFunctionType localFunction(const EntityType& entity,
                                      const IntersectionType &intersection)
      {
        return LocalFunctionType(entity, intersection, *this);
      }
 
      const DiscreteFunctionSpaceType& space() const
      {
        return function().space();
      }
 
      const GridPartType& gridPart() const
      {
        return function().gridPart();
      }
 
      const GridFunctionType& function() const
      {
        return gridFunction_();
      }
 
      const IndicatorType& indicator() const
      {
        return indicator_();
      }
 
      const OuterIndicatorType& outerIndicator() const
      {
        return outerIndicator_();
      }
 
     protected:
      // Remove any outer redundant parenthesis and wrap into "B(.)"
      static std::string makeName(const std::string& name)
      {
        std::string result(name);
        trimParenthesis(result);
        result = "B(" + result + ")";
        return result;
      }
 
      ExpressionStorage<GridFunctionType> gridFunction_;
      ExpressionStorage<OuterIndicatorType> outerIndicator_;
      ExpressionStorage<IndicatorType> indicator_;
    };
 
    template<class GridFunction, class Indicator>
    class BoundarySupportedFunction<GridFunction, Indicator>::LocalFunction
    {
     public:
      typedef GridFunction GridFunctionType;
      typedef typename GridFunctionType::FunctionSpaceType FunctionSpaceType;
      typedef typename GridFunctionType::LocalFunctionType BulkLocalFunctionType;
      typedef Indicator IndicatorType;
      typedef BoundarySupportedFunction<GridFunction, Indicator> DiscreteFunctionType;
 
      enum {
        emptySupport = DiscreteFunctionType::emptySupport,
        globalSupport = DiscreteFunctionType::globalSupport
      };
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionType::GridPartType GridPartType;
      typedef typename DiscreteFunctionType::EntityType EntityType;
      typedef typename DiscreteFunctionType::IntersectionType IntersectionType;
 
      enum {
        dimRange = DiscreteFunctionSpaceType::dimRange,
        dimDomain = GridPartType::dimensionworld,
        dimension = GridPartType::dimension
      };
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
 
      LocalFunction(const EntityType& entity, const DiscreteFunctionType& df)
        : intersection_(0),
          df_(df),
          localFunction_(df_.function().localFunction(entity)),
          supported_(true)
      {}
 
      LocalFunction(const DiscreteFunctionType& df)
        : intersection_(0),
          df_(df),
          localFunction_(df_.function()),
          supported_(true)
      {}
 
      LocalFunction(const EntityType& entity,
                    const IntersectionType& intersection,
                    const DiscreteFunctionType& df)
        : df_(df),
          localFunction_(BulkLocalFunctionType(df_.function()))
      {
        init(entity, intersection);
      }
 
      template<class PointType>
      void evaluate(const PointType &x, RangeType &ret) const
      {
        if (!supported()) {
          ret = 0;
        } else {
          localFunction_.evaluate(x, ret);
        }
      }
 
      template<class PointType>
      void jacobian(const PointType &x, JacobianRangeType &ret) const
      {
        if (!supported()) {
          ret = 0;
        } else {
          localFunction_.jacobian(x, ret);
        }
      }
 
      template<class PointType>
      void hessian(const PointType &x, HessianRangeType &ret) const
      {
        if (!supported()) {
          ret = 0;
        } else {
          localFunction_.hessian(x, ret);
        }
      }
 
      template<class QuadratureType, class VectorType>
      void evaluateQuadrature(const QuadratureType& quadrature, VectorType& values) const
      {
        assert(values.size() == quadrature.nop());
        DiscreteFunctionType::evaluateQuadratureImp(*this, quadrature, values, values[0]);
      }
 
      int order() const
      {
        return localFunction_.order();
      }
 
      void init(const EntityType &entity)
      {
        localFunction_.init(entity);
        supported_ = true;
      }
 
      void init(const EntityType& entity, const IntersectionType& isect)
      {
        intersection_ = &isect;
        BoundaryTraitsType::init(localFunction_, entity, intersection());
        supported_ = globalSupport || (!emptySupport && df_.indicator().applies(intersection()));
      }
 
      const EntityType& entity() const
      {
        localFunction_.entity();
      }
 
      const IntersectionType& intersection() const
      {
        assert(intersection_ != 0);
        return *intersection_;
      }
 
      bool supported() const
      {
        return globalSupport || (!emptySupport && supported_);
      }
 
     protected:
      const IntersectionType* intersection_;
      const DiscreteFunctionType& df_;
      BulkLocalFunctionType localFunction_;
      bool supported_; 
    };
 
    template<class GridFunction, class Indicator = EntireBoundaryIndicatorType>
    BoundarySupportedFunction<GridFunction, Indicator>
    boundarySupportedFunction(const Fem::Function<typename GridFunction::FunctionSpaceType,
                                                  GridFunction>& f,
                              const BoundaryIndicatorInterface<Indicator>&
                              c = EntireBoundaryIndicatorType())
    {
      typedef BoundarySupportedFunction<GridFunction, Indicator> FunctionType;
 
      return FunctionType(f, c);
    }
 
    template<class GridFunction, class Indicator>
    BoundarySupportedFunction<GridFunction, Indicator>
    boundarySupportedFunction(const BoundarySupportedFunction<GridFunction, Indicator>& f,
                              const EntireBoundaryIndicatorType& = EntireBoundaryIndicatorType())
    {
      return f;
    }
 
  } // ACFem::
 
} // Dune::
 
#endif // __DUNE_ACFEM_BOUNDARY_SUPPORTED_FUNCTION_HH__