elementquadrature.hh

#ifndef DUNE_FEM_ELEMENTQUADRATURE_HH
#define DUNE_FEM_ELEMENTQUADRATURE_HH
 
#include <dune/fem/quadrature/elementpointlistbase.hh>
#include <dune/fem/quadrature/caching/twistutility.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    template< typename GridPartImp, class IntegrationPointList >
    class Agglomeration;
 
    template< class GridPartImp, int codim, template< class, int > class QuadratureTraits >
    struct ElementQuadratureTraits
    {
      // type of single coordinate
      typedef typename GridPartImp :: ctype ctype;
 
      // dimension of quadrature
      enum { dimension = GridPartImp ::  dimension };
 
      // codimension of quadrature
      enum { codimension = codim };
 
      // type of used integration point list
      typedef Quadrature< ctype, dimension-codim, QuadratureTraits > IntegrationPointListType;
 
      // type of local coordinate (with respect to the codim-0 entity)
      typedef typename Quadrature< ctype, dimension, QuadratureTraits > :: CoordinateType
        CoordinateType;
    };
 
    template< class GridPartImp, int codim, class IntegrationTraits, bool isQuadrature >
    class ElementQuadratureImpl;
 
 
    template< class GridPartImp, class IntegrationTraits, bool isQuadrature >
    class ElementQuadratureImpl< GridPartImp, 0, IntegrationTraits, isQuadrature >
    : public ElementPointListBase< GridPartImp, 0, IntegrationTraits >
    {
      typedef ElementQuadratureImpl< GridPartImp, 0, IntegrationTraits, isQuadrature > This;
      typedef ElementPointListBase< GridPartImp, 0, IntegrationTraits >  Base;
 
    public:
      typedef typename Base :: GridPartType  GridPartType;
 
      typedef typename Base::CoordinateType  CoordinateType;
 
      typedef QuadraturePointWrapper< This > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< This > IteratorType;
 
      // for compatibility
      typedef typename Base :: EntityType     EntityType;
 
      typedef typename Base :: IntegrationPointListType  IntegrationPointListType;
 
    protected:
      template <class QuadratureKeyType>
      IntegrationPointListType createQuadrature(  const EntityType &entity, const QuadratureKeyType& quadKey, const bool checkGeomType )
      {
        const GeometryType geomType = entity.type();
        if( checkGeomType && ! geomType.isNone() )
        {
          // return default element quadratures for given geometry type
          return IntegrationPointListType( geomType, quadKey );
        }
        else // compute weights and points based on sub-triangulation
        {
          typedef Agglomeration< GridPartType, IntegrationPointListType > AgglomerationType;
          return AgglomerationType::computeQuadrature( entity, quadKey );
        }
      }
 
      using Base::quadImp;
 
    public:
      using Base::localPoint;
      using Base::nop;
 
      template <class QuadratureKeyType>
      ElementQuadratureImpl( const EntityType &entity, const QuadratureKeyType& quadKey, const bool checkGeomType = isQuadrature )
      : Base( createQuadrature( entity, quadKey, checkGeomType ) )
      {}
 
      template <class QuadratureKeyType>
      ElementQuadratureImpl( const GeometryType &geometry, const QuadratureKeyType& quadKey )
      : Base( geometry, quadKey )
      {}
 
      const QuadraturePointWrapperType operator[] ( const size_t i ) const
      {
        return QuadraturePointWrapperType( *this, i );
      }
 
      IteratorType begin () const noexcept { return IteratorType( *this, 0 ); }
      IteratorType end () const noexcept { return IteratorType( *this, nop() ); }
 
      const CoordinateType &point ( const size_t i ) const
      {
        return localPoint( i );
      }
 
      auto weight( size_t i ) const
      {
        return quadImp().weight( i );
      }
    };
 
 
 
    template< class GridPartImp, class IntegrationTraits, bool isQuadrature >
    class ElementQuadratureImpl< GridPartImp, 1, IntegrationTraits, isQuadrature >
    : public ElementPointListBase< GridPartImp, 1, IntegrationTraits >
    {
      typedef ElementQuadratureImpl< GridPartImp, 1, IntegrationTraits, isQuadrature > This;
      typedef ElementPointListBase< GridPartImp, 1, IntegrationTraits > Base;
 
    public:
      typedef GridPartImp GridPartType;
 
      static const int dimension = Base::dimension;
 
      typedef typename Base::CoordinateType CoordinateType;
 
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
 
      typedef QuadraturePointWrapper< This > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< This > IteratorType;
 
      typedef This NonConformingQuadratureType;
 
 
      // for compatibility
      typedef typename GridPartType::TwistUtilityType  TwistUtilityType;
      typedef IntersectionIteratorType IntersectionIterator;
 
 
      using Base::localPoint;
      using Base::elementGeometry;
      using Base::nop;
 
      template <class QuadratureKeyType>
      ElementQuadratureImpl( const GridPartType &gridPart,
                                       const IntersectionType &intersection,
                                       const QuadratureKeyType& quadKey,
                                       const typename Base :: Side side )
      : Base( getPointList( intersection, quadKey, side ) ),
        side_(side),
        intersection_(intersection),
        referenceGeometry_( side == Base::INSIDE ?  intersection.geometryInInside() : intersection.geometryInOutside())
      {}
 
      const QuadraturePointWrapperType operator[] ( size_t i ) const
      {
        return QuadraturePointWrapperType( *this, i );
      }
 
      IteratorType begin () const noexcept { return IteratorType( *this, 0 ); }
      IteratorType end () const noexcept { return IteratorType( *this, nop() ); }
      typename Base :: Side side() const { return side_; }
      bool isInside() const { return side_ == Base::INSIDE; }
 
      const CoordinateType &point ( size_t i ) const
      {
        dummy_ = referenceGeometry_.global( localPoint( i ) );
        return dummy_;
      }
 
      auto weight( size_t i ) const
      {
        return quadImp().weight( i );
      }
 
      using Base::localFaceIndex;
 
      const IntersectionType &intersection() const
      {
        return intersection_;
      }
 
    protected:
      using Base::quadImp;
 
      Base getPointList ( const IntersectionType &intersection, const int order,
                          const typename Base :: Side side )
      {
        switch( side )
        {
          case Base :: INSIDE:
            return Base( TwistUtilityType::elementGeometry( intersection, true ),
                         intersection.type(), intersection.indexInInside(), order );
 
          case Base ::OUTSIDE:
            return Base( TwistUtilityType::elementGeometry( intersection, false ),
                         intersection.type(), intersection.indexInOutside(), order );
 
          default:
            DUNE_THROW( InvalidStateException, "ElementIntegrationPointList: side must either be INSIDE or OUTSIDE." );
        }
      }
 
    private:
      typedef typename IntersectionIteratorType::Intersection::LocalGeometry ReferenceGeometry;
 
      const typename Base :: Side side_;
      const IntersectionType &intersection_;
      ReferenceGeometry referenceGeometry_;
      mutable CoordinateType dummy_;
    };
 
    template< class GridPart, int codim, template <class, int> class IntegrationTraits = DefaultQuadratureTraits >
    using ElementIntegrationPointList = ElementQuadratureImpl< GridPart, codim, ElementQuadratureTraits< GridPart, codim, IntegrationTraits>, false >;
 
    template< class GridPart, int codim, template <class, int> class IntegrationTraits = DefaultQuadratureTraits >
    using ElementQuadrature = ElementQuadratureImpl< GridPart, codim, ElementQuadratureTraits< GridPart, codim, IntegrationTraits>, true >;
 
    template<class GridPart, class Entity>
    static inline auto elementQuadrature(const GridPart& gridPart, const Entity& entity, unsigned quadOrder)
    {
      using Quadrature = Dune::Fem::ElementQuadrature<GridPart, 0>;
      return Quadrature(entity, quadOrder);
    }
 
  } // namespace Fem
 
} // namespace Dune
 
#include "agglomerationquadrature.hh"
 
#endif // #ifndef DUNE_FEM_ELEMENTQUADRATURE_HH