divergencefunctional.hh

#ifndef __DUNE_ACFEM_DIVERGENCE_FUNCTIONAL_HH__
#define __DUNE_ACFEM_DIVERGENCE_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 DivergenceFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             template<class> class QuadratureTraits>
    class LocalDivergenceFunctional;
 
    template<class DiscreteFunctionSpace,
             class GridFunction,
             template<class> class QuadratureTraits>
    struct DivergenceFunctionalTraits
      : public
    LinearFunctionalTraitsDefault<
      DivergenceFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits>,
      LocalDivergenceFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits> >
    {
     private:
      typedef DivergenceFunctionalTraits 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 DivergenceFunctional
      : public DivergenceFunctionalTraits<DiscreteFunctionSpace, GridFunction, QuadratureTraits>::ExpressionBaseType
    {
      typedef DivergenceFunctional ThisType;
     public:
      typedef
      DivergenceFunctionalTraits<DiscreteFunctionSpace, GridFunction, QuadratureTraits>
      TraitsType;
     private:
      typedef typename TraitsType::ExpressionBaseType BaseType;
      typedef ExpressionStorage<GridFunction> DivergenceFunctionStorage;
     public:
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef GridFunction DivergenceFunctionType;
      typedef
      Fem::Function<typename DivergenceFunctionType::FunctionSpaceType, DivergenceFunctionType>
      DivergenceFunctionInterfaceType;
      typedef
      Fem::DiscreteFunctionSpaceInterface<typename DiscreteFunctionSpaceType::Traits>
      DiscreteFunctionSpaceInterfaceType;
 
      static_assert(DivergenceFunctionType::DiscreteFunctionSpaceType::FunctionSpaceType::dimRange
                    ==
                    DivergenceFunctionType::DiscreteFunctionSpaceType::FunctionSpaceType::dimDomain,
                    "must have dimRange == dimDomain for divergence");
 
      static_assert(DiscreteFunctionSpaceType::FunctionSpaceType::dimRange == 1,
                    "must have dimRange == 1 for test space, div is a scalar quantity");
 
      DivergenceFunctional(const DiscreteFunctionSpaceInterfaceType& space,
                           const DivergenceFunctionInterfaceType& divFct,
                           unsigned degree = (DiscreteFunctionSpaceType::polynomialOrder
                                              +
                                              DivergenceFunctionType::DiscreteFunctionSpaceType::polynomialOrder - 1))
        : BaseType(static_cast<const DiscreteFunctionSpaceType&>(space)),
          divFct_(static_cast<const DivergenceFunctionType&>(divFct)),
          degree_(degree > 111 ? (2*space.order() - 1) : degree)
      {}
 
      using BaseType::operator();
      using BaseType::coefficients;
      using BaseType::space;
      using BaseType::localFunctional;
 
      const DivergenceFunctionType& function() const
      {
        return divFct_();
      }
 
      unsigned quadratureDegree() const
      {
        return degree_;
      }
 
      std::string name() const
      {
        return "(div " + divFct_().name() + ", . )";
      }
 
     protected:
      DivergenceFunctionStorage divFct_;
      const unsigned degree_;
    };
 
    template<class DiscreteFunctionSpace, class GridFunction, template<class> class QuadratureTraits>
    class LocalDivergenceFunctional
      : public
    LocalLinearFunctionalDefault<DiscreteFunctionSpace,
                                 DivergenceFunctionalTraits<DiscreteFunctionSpace,
                                                            GridFunction,
                                                            QuadratureTraits> >
    {
     public:
      typedef
      DivergenceFunctionalTraits<DiscreteFunctionSpace, GridFunction, QuadratureTraits>
      TraitsType;
     private:
      typedef LocalDivergenceFunctional ThisType;
      typedef LocalLinearFunctionalDefault<DiscreteFunctionSpace, TraitsType> BaseType;
     public:
      typedef
      DivergenceFunctional<DiscreteFunctionSpace, GridFunction, QuadratureTraits>
      FunctionalType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType FieldType;
      typedef FieldType RangeType;
     protected:
      typedef typename FunctionalType::DivergenceFunctionType DivergenceFunctionType;
      typedef typename DivergenceFunctionType::LocalFunctionType LocalDivergenceFunctionType;
      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:
      LocalDivergenceFunctional(const FunctionalType& phi)
        : BaseType(phi),
          divLocal_(phi.function()),
          degree_(phi.quadratureDegree())
      {}
 
      LocalDivergenceFunctional(const EntityType& entity, const FunctionalType& phi)
        : BaseType(entity, phi),
          divLocal_(phi.function()),
          degree_(phi.quadratureDegree())
      {}
 
     public:
 
      void init(const EntityType& entity)
      {
        divLocal_.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 = divLocal.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 LocalDivergenceFunctionType::JacobianRangeType Df;
          divLocal_.jacobian(quadrature[pt], Df);
 
          typename DiscreteFunctionSpaceType::RangeType divf(0);
          for (int i = 0; i < LocalDivergenceFunctionType::dimDomain; ++i) {
            divf[0] += Df[i][i];
          }
 
          // add to local function
          divf *= weight * c;
          coeffs.axpy(quadrature[pt], divf);
        }  // 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_;
      LocalDivergenceFunctionType divLocal_;
      const unsigned degree_;
    };
 
    template<class DiscreteFunction, class GridFunction>
    static inline
    double
    divergenceFunctional(const
                         Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& u,
                         const
                         Fem::DiscreteFunctionInterface<DiscreteFunction>& p)
    {
      DivergenceFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType, GridFunction>
        functional(p.space(), u);
      return functional(p);
    }
 
    template<class DiscreteSpaceTraits, class GridFunction>
    static inline
    DivergenceFunctional<typename DiscreteSpaceTraits::DiscreteFunctionSpaceType, GridFunction>
    divergenceFunctional(const
                         Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& u,
                         const
                         Fem::DiscreteFunctionSpaceInterface<DiscreteSpaceTraits>& space)
    {
      typedef
        DivergenceFunctional<typename DiscreteSpaceTraits::DiscreteFunctionSpaceType, GridFunction>
        FunctionalType;
 
      return FunctionalType(space, u);
    }
 
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_DIVERGENCE_FUNCTIONAL_HH__