l2innerproductfunctional.hh

#ifndef __DUNE_ACFEM_L2_INNER_PRODUCT_FUNCTIONAL_HH__
#define __DUNE_ACFEM_L2_INNER_PRODUCT_FUNCTIONAL_HH__
 
#include "functionalexpressionbase.hh"
#include "../../common/quadrature.hh"
#include "../../functions/gridfunctionexpression.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class L2InnerProductFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             template<class> class QuadratureTraits>
    class LocalL2InnerProductFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             template<class> class QuadratureTraits>
    struct L2InnerProductFunctionalTraits
      : public
    LinearFunctionalTraitsDefault<
      L2InnerProductFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits>,
      LocalL2InnerProductFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits> >
    {
     private:
      typedef L2InnerProductFunctionalTraits ThisType;
     public:
      // if GridFunction inherits from ZeroExpression, then we are zero, too
      typedef typename
      std::conditional<ExpressionTraits<GridFunction>::isZero,
                       ZeroFunctionalExpression<DiscreteFunctionSpace, ThisType>,
                       DiscreteLinearFunctionalExpression<DiscreteFunctionSpace, ThisType>
                       >::type
      ExpressionBaseType;
    };
 
    template<class DiscreteFunctionSpace, class GridFunction, template<class> class QuadratureTraits>
    class L2InnerProductFunctional
      : public L2InnerProductFunctionalTraits<DiscreteFunctionSpace, GridFunction, QuadratureTraits>::ExpressionBaseType
    {
      typedef L2InnerProductFunctional ThisType;
     public:
      typedef
      L2InnerProductFunctionalTraits<DiscreteFunctionSpace, GridFunction, QuadratureTraits>
      TraitsType;
     private:
      typedef typename TraitsType::ExpressionBaseType BaseType;
      typedef ExpressionStorage<GridFunction> L2FunctionStorage;
     public:
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef GridFunction L2FunctionType;
      typedef
      Fem::Function<typename L2FunctionType::FunctionSpaceType, L2FunctionType>
      L2FunctionInterfaceType;
 
      L2InnerProductFunctional(const DiscreteFunctionSpaceType& space,
                               const L2FunctionInterfaceType& L2Fct,
                               unsigned degree = 2*DiscreteFunctionSpaceType::polynomialOrder)
        : BaseType(space),
          L2Fct_(static_cast<const L2FunctionType&>(L2Fct)),
          degree_(degree)
      {}
 
      using BaseType::operator();
      using BaseType::coefficients;
      using BaseType::space;
      using BaseType::localFunctional;
 
      const L2FunctionType& function() const
      {
        return L2Fct_();
      }
 
      unsigned quadratureDegree() const
      {
        return degree_;
      }
 
      std::string name() const
      {
        return "(" + L2Fct_().name() + ", . )";
      }
 
     protected:
      L2FunctionStorage L2Fct_;
      const unsigned degree_;
    };
 
    template<class DiscreteFunctionSpace, class GridFunction, template<class> class QuadratureTraits>
    class LocalL2InnerProductFunctional
      : public LocalLinearFunctionalDefault<DiscreteFunctionSpace,
                                            L2InnerProductFunctionalTraits<DiscreteFunctionSpace,
                                                                           GridFunction,
                                                                           QuadratureTraits> >
    {
     public:
      typedef
      L2InnerProductFunctionalTraits<DiscreteFunctionSpace,
                                     GridFunction,
                                     QuadratureTraits>
      TraitsType;
     private:
      typedef LocalL2InnerProductFunctional ThisType;
      typedef LocalLinearFunctionalDefault<DiscreteFunctionSpace, TraitsType> BaseType;
     public:
      typedef
      L2InnerProductFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits>
      FunctionalType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType FieldType;
      typedef FieldType RangeType;
     protected:
      typedef typename FunctionalType::L2FunctionType L2FunctionType;
      typedef typename L2FunctionType::LocalFunctionType LocalL2FunctionType;
      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;
 
     public:
      LocalL2InnerProductFunctional(const FunctionalType& phi)
        : BaseType(phi),
          l2Local_(phi.function()),
          degree_(phi.quadratureDegree())
      {}
 
      LocalL2InnerProductFunctional(const EntityType& entity, const FunctionalType& phi)
        : BaseType(entity, phi),
          l2Local_(phi.function()),
          degree_(phi.quadratureDegree())
      {}
 
     public:
 
      void init(const EntityType& entity)
      {
        l2Local_.init(entity);
        BaseType::init(entity);
      }
 
      template<class LocalFunction>
      void coefficients(const RangeType& c, LocalFunction& coeffs) const
      {
        const auto& entity = coeffs.entity();
        const auto& geometry = entity.geometry();
 
        // obtain quadrature order
        // const int quadOrder = l2Local.order() + coeffs.order();
        const int quadOrder = degree_;
 
        MassQuadratureType quadrature(entity, quadOrder);
        const size_t numQuadraturePoints = quadrature.nop();
        for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
 
          const typename MassQuadratureType::CoordinateType &x = quadrature.point(pt);
          const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 
          typename LocalL2FunctionType::RangeType f;
          l2Local_.evaluate(quadrature[pt], f);
 
          // add to local function
          f *= weight * c;
          coeffs.axpy(quadrature[pt], f);
        }  // end quad 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_;
      LocalL2FunctionType l2Local_;
      const unsigned degree_;
    };
 
    template<class DiscreteFunction, class GridFunction>
    static inline
    double
    l2InnerProductFunctional(const
                             Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& u,
                             const
                             Fem::DiscreteFunctionInterface<DiscreteFunction> &p)
    {
      L2InnerProductFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType, GridFunction>
        functional(p.space(), u);
      return functional(p);
    }
 
    template<class DiscreteSpaceTraits, class GridFunction>
    static inline
    L2InnerProductFunctional<typename DiscreteSpaceTraits::DiscreteFunctionSpaceType, GridFunction>
    l2InnerProductFunctional(const
                             Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& u,
                             const
                             Fem::DiscreteFunctionSpaceInterface<DiscreteSpaceTraits>& space)
    {
      typedef
        L2InnerProductFunctional<typename DiscreteSpaceTraits::DiscreteFunctionSpaceType, GridFunction>
        FunctionalType;
 
      return FunctionalType(space, u);
    }
 
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_L2_INNER_PRODUCT_FUNCTIONAL_HH__