l2boundaryfunctional.hh

#ifndef __DUNE_ACFEM_L2_SCALAR_PRODUCT_FUNCTIONAL_HH__
#define __DUNE_ACFEM_L2_SCALAR_PRODUCT_FUNCTIONAL_HH__
 
#include "linearfunctional.hh"
#include "../../models/boundaryindicator.hh"
#include "../../common/quadrature.hh"
#include "../../functions/gridfunctionexpression.hh"
#include "../../functions/boundaryfunctionexpression.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             class BndryIndicator = EntireBoundaryIndicatorType,
             template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class L2BoundaryFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             class BndryIndicator,
             template<class> class QuadratureTraits>
    class LocalL2BoundaryFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             class BndryIndicator,
             template<class> class QuadratureTraits>
    struct L2BoundaryFunctionalTraits
      : public
    LinearFunctionalTraitsDefault<
      L2BoundaryFunctional<DiscreteFunctionSpace,
                           GridFunction,
                           BndryIndicator,
                           QuadratureTraits>,
      LocalL2BoundaryFunctional<DiscreteFunctionSpace,
                                GridFunction,
                                BndryIndicator,
                                QuadratureTraits> >
    {
     private:
      typedef L2BoundaryFunctionalTraits ThisType;
     public:
      /*There are two possibilities for zero-expressions:
       *
       * - GridFunction is a zero-expression (either zero-function of zero-supported)
       * - BndryIndicator is a zero-expression
       */
      typedef typename
      std::conditional<ExpressionTraits<GridFunction>::isZero || ExpressionTraits<BndryIndicator>::isZero,
                       ZeroFunctionalExpression<DiscreteFunctionSpace, ThisType>,
                       DiscreteLinearFunctionalExpression<DiscreteFunctionSpace, ThisType>
                       >::type
        ExpressionBaseType;
    };
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             class BndryIndicator,
             template<class> class QuadratureTraits>
    class L2BoundaryFunctional
      : public L2BoundaryFunctionalTraits<DiscreteFunctionSpace,
                                          GridFunction,
                                          BndryIndicator,
                                          QuadratureTraits>::ExpressionBaseType
    {
      typedef L2BoundaryFunctional ThisType;
     public:
      typedef
      L2BoundaryFunctionalTraits<DiscreteFunctionSpace, GridFunction, BndryIndicator, QuadratureTraits>
      TraitsType;
      typedef GridFunction L2FunctionType;
      typedef BndryIndicator OuterIndicatorType;
     protected:
      typedef typename TraitsType::ExpressionBaseType BaseType;
     public:
      typedef
      decltype(std::declval<OuterIndicatorType>() * std::declval<L2FunctionType>())
        BndryFunctionType;
      typedef typename BndryFunctionType::IndicatorType IndicatorType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
      typedef
      Fem::Function<typename L2FunctionType::FunctionSpaceType, L2FunctionType>
      L2FunctionInterfaceType;
      typedef
      Fem::DiscreteFunctionSpaceInterface<typename DiscreteFunctionSpaceType::Traits>
      DiscreteFunctionSpaceInterfaceType;
      typedef
      BoundaryIndicatorInterface<OuterIndicatorType>
      IndicatorInterfaceType;
 
      L2BoundaryFunctional(const DiscreteFunctionSpaceInterfaceType& space,
                           const L2FunctionInterfaceType& L2Fct,
                           const IndicatorInterfaceType& indicator = OuterIndicatorType(),
                           unsigned degree = 2*DiscreteFunctionSpaceType::polynomialOrder)
        : BaseType(static_cast<const DiscreteFunctionSpaceType&>(space)),
          bndryFct_(asImp(indicator) * asImp(L2Fct)),
          degree_(degree)
      {}
 
      using BaseType::operator();
      using BaseType::coefficients;
      using BaseType::space;
      using BaseType::localFunctional;
 
      std::string name() const
      {
        return "[" + bndryFct_().name() + ", . ]";
      }
 
      const BndryFunctionType& function() const
      {
        return bndryFct_();
      }
 
      unsigned quadratureDegree() const
      {
        return degree_;
      }
 
      const IndicatorType& indicator() const
      {
        return function().indicator();
      }
 
      IndicatorType& indicator()
      {
        return function().indicator();
      }
 
     protected:
      ExpressionStorage<BndryFunctionType> bndryFct_;
      const unsigned degree_;
    };
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             class BndryIndicator,
             template<class> class QuadratureTraits>
    class LocalL2BoundaryFunctional
      : public LocalLinearFunctionalDefault<DiscreteFunctionSpace,
                                            L2BoundaryFunctionalTraits<DiscreteFunctionSpace,
                                                                       GridFunction,
                                                                       BndryIndicator,
                                                                       QuadratureTraits> >
    {
     public:
      typedef
      L2BoundaryFunctionalTraits<DiscreteFunctionSpace,
                                 GridFunction,
                                 BndryIndicator,
                                 QuadratureTraits>
      TraitsType;
     private:
      typedef LocalL2BoundaryFunctional ThisType;
      typedef LocalLinearFunctionalDefault<DiscreteFunctionSpace, TraitsType> BaseType;
     public:
      typedef
      L2BoundaryFunctional<DiscreteFunctionSpace, GridFunction, BndryIndicator, QuadratureTraits>
      FunctionalType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType FieldType;
      typedef FieldType RangeType;
     protected:
      typedef typename FunctionalType::IndicatorType IndicatorType;
      typedef typename FunctionalType::BndryFunctionType BndryFunctionType;
      typedef typename BndryFunctionType::LocalFunctionType LocalBndryFunctionType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef QuadratureTraits<GridPartType> QuadratureTraitsType;
      typedef typename QuadratureTraitsType::BulkQuadratureType QuadratureType;
      typedef typename QuadratureTraitsType::BulkMassQuadratureType MassQuadratureType;
      typedef typename QuadratureTraitsType::FaceQuadratureType FaceQuadratureType;
      typedef typename QuadratureTraitsType::FaceMassQuadratureType FaceMassQuadratureType;
 
      typedef Dune::Fem::IntersectionQuadrature<FaceMassQuadratureType, true> IntersectionMassQuadratureType;
      typedef typename IntersectionMassQuadratureType::FaceQuadratureType     BndryMassQuadratureType;
 
     public:
      LocalL2BoundaryFunctional(const FunctionalType& phi)
        : BaseType(phi),
          bndryLocal_(phi.function()),
          indicator_(phi.function().indicator()),
          gridPart_(phi.space().gridPart()),
          degree_(phi.quadratureDegree())
      {}
 
      LocalL2BoundaryFunctional(const EntityType& entity, const FunctionalType& phi)
        : BaseType(entity, phi),
          bndryLocal_(phi.function()),
          indicator_(phi.function().indicator()),
          gridPart_(phi.space().gridPart()),
          degree_(phi.quadratureDegree())
      {}
 
     public:
 
      template<class LocalFunction>
      RangeType operator()(const LocalFunction& wLocal) const
      {
        if (ExpressionTraits<IndicatorType>::isZero) {
          // nothing to do
          return 0;
        }
 
        const auto& entity = wLocal.entity();
 
        if (!entity.hasBoundaryIntersections()) {
          return 0;
        }
 
        return BaseType::operator()(wLocal);
      }
 
      template<class LocalFunction>
      void coefficients(const RangeType& c, LocalFunction& coeffs) const
      {
        if (ExpressionTraits<IndicatorType>::isZero) {
          return;
        }
 
        const auto& entity = coeffs.entity();
 
        if (!entity.hasBoundaryIntersections()) {
          return;
        }
 
        const auto iend = gridPart_.iend(entity);
        for (auto it = gridPart_.ibegin(entity); it != iend; ++it) {
 
          const auto &intersection = *it;
 
          if (intersection.neighbor() || !intersection.boundary()) {
            continue;
          }
 
          if (!indicator_.applies(intersection)) {
            continue;
          }
 
          bndryLocal_.init(entity, intersection); // assume this is cheap ... in this context
 
          const int quadOrder = degree_;
 
          IntersectionMassQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
          const BndryMassQuadratureType &bndryQuad = intersectionQuad.inside();
          const int numQuadraturePoints = bndryQuad.nop();
 
          for (int pt = 0; pt < numQuadraturePoints; ++pt) {
 
            // quadrature weight
            const double weight = bndryQuad.weight(pt);
 
            // integration element
            const double integrationElement = intersection.geometry().integrationElement(bndryQuad.localPoint(pt));
 
            typename LocalBndryFunctionType::RangeType f;
            bndryLocal_.evaluate(bndryQuad[pt], f);
            f *= weight * integrationElement * c;
            coeffs.axpy(bndryQuad[pt], f);
 
          } // end boundary quadrature loop
 
        } // end Intersection loop
 
      }
 
      template<class LocalFunction>
      void coefficients(LocalFunction& coeffs) const
      {
        BaseType::coefficients(coeffs);
      }
 
      using BaseType::operator();
      using BaseType::functional;
      using BaseType::entity;
 
     protected:
      using BaseType::entity_;
      mutable LocalBndryFunctionType bndryLocal_;
      const IndicatorType& indicator_;
      const GridPartType& gridPart_;
      const unsigned degree_;
    };
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_L2_SCALAR_PRODUCT_FUNCTIONAL_HH__