quadrature.hh

#ifndef DUNE_FEM_QUADRATURE_HH
#define DUNE_FEM_QUADRATURE_HH
 
#include <cassert>
#include <memory>
 
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
 
#include <dune/geometry/quadraturerules.hh>
 
#include <dune/grid/common/grid.hh>
 
#include <dune/fem/common/coordinate.hh>
#include <dune/fem/common/utility.hh>
#include <dune/fem/storage/envelope.hh>
#include <dune/fem/version.hh>
#include <dune/fem/quadrature/quadprovider.hh>
#include <dune/fem/quadrature/defaultquadratures.hh>
 
namespace Dune
{
  namespace Fem
  {
  template< class Quadrature >
  class QuadraturePointWrapper
  {
  public:
    typedef Quadrature QuadratureType;
 
    typedef typename QuadratureType::CoordinateType CoordinateType;
    typedef typename QuadratureType::RealType RealType;
    typedef typename QuadratureType::LocalCoordinateType LocalCoordinateType;
 
  private:
    typedef QuadraturePointWrapper< QuadratureType > ThisType;
 
  protected:
    const QuadratureType &quadrature_;
    const unsigned int index_;
 
  public:
    QuadraturePointWrapper ( const QuadratureType &quadrature, unsigned int index )
      : quadrature_( quadrature ), index_( index )
    {}
 
    ThisType &operator= ( const ThisType & ) = delete;
 
    inline const QuadratureType &quadrature () const
    {
      return quadrature_;
    }
 
    unsigned int index () const { return index_; }
    const CoordinateType &position () const { return quadrature().point( index() ); }
    auto weight () const { return quadrature().weight( index() ); }
    const LocalCoordinateType &localPosition () const { return quadrature().localPoint( index() ); }
  };
 
 
 
  template< class Quadrature >
  static inline typename Quadrature::CoordinateType
  coordinate ( const QuadraturePointWrapper< Quadrature > &x )
  {
    return x.position();
  }
 
  template< class Quadrature >
  class QuadraturePointIterator
  {
    typedef QuadraturePointIterator< Quadrature > ThisType;
 
  public:
    typedef std::forward_iterator_tag iterator_category;
    typedef QuadraturePointWrapper< Quadrature > value_type;
    typedef std::ptrdiff_t difference_type;
    typedef Envelope< QuadraturePointWrapper< Quadrature > > pointer;
    typedef QuadraturePointWrapper< Quadrature > reference;
 
    QuadraturePointIterator () noexcept = default;
 
    QuadraturePointIterator ( const Quadrature &quadrature, std::size_t point ) noexcept
      : quadrature_( &quadrature ), point_( point  )
    {}
 
    reference operator* () const noexcept { return value_type( quadrature(), point_ ); }
    pointer operator-> () const noexcept { return pointer( value_type( quadrature(), point_ ) ); }
 
    bool operator== ( const ThisType &other ) const noexcept { return (point_ == other.point_); }
    bool operator!= ( const ThisType &other ) const noexcept { return (point_ != other.point_); }
 
    ThisType &operator++ () noexcept { ++point_; return *this; }
    ThisType operator++ ( int ) noexcept { ThisType copy( *this ); ++(*this); return copy; }
 
    const Quadrature &quadrature () const noexcept { assert( quadrature_ ); return *quadrature_; }
 
  protected:
    const Quadrature *quadrature_ = nullptr;
    std::size_t point_ = 0;
  };
 
 
 
  template< typename FieldImp, int dim,
            template< class, int > class IntegrationTraits,
            bool isQuadrature = false >
  class IntegrationPointList
  {
  public:
    typedef FieldImp FieldType;
 
    enum { dimension = dim };
 
    typedef IntegrationTraits< FieldType, dimension > Traits;
 
  private:
    typedef IntegrationPointList< FieldType, dimension, IntegrationTraits > ThisType;
 
    typedef QuadratureProvider< FieldType, dimension, IntegrationTraits >
      QuadratureProviderType;
 
    typedef typename std::conditional< isQuadrature,
                             const FieldType&,
                             FieldType > :: type WeightReturnType;
  public:
    typedef typename Traits :: IntegrationPointListType IntegrationPointListType;
 
    typedef typename IntegrationPointListType :: CoordinateType CoordinateType;
 
    typedef QuadraturePointWrapper< ThisType > QuadraturePointWrapperType;
 
    typedef std::shared_ptr< const IntegrationPointListType >  IntegrationPointListStorageType;
 
    enum { codimension = 0 };
 
  protected:
    // when the integration point list is obtained from the
    // QuadratureProvider then it should not be deleted
    struct NoDelete
    {
      void operator()( const IntegrationPointListType* ) {}
    };
 
    IntegrationPointListStorageType ipListPtr_;
 
  public:
    template <class FactoryTraits>
    inline IntegrationPointList ( const FactoryTraits traits,
                                  const GeometryType &geometryType,
                                  const typename FactoryTraits::QuadratureKeyType& quadKey )
    : ipListPtr_( &QuadratureProviderType :: getQuadrature( geometryType, quadKey ), NoDelete() )
    //: ipListPtr_( &QuadratureProviderType :: getQuadrature( traits, geometryType, quadKey ), NoDelete() )
    {
    }
 
    template <class QuadratureKeyType>
    inline IntegrationPointList ( const GeometryType &geometryType,
                                  const QuadratureKeyType& quadKey )
    : IntegrationPointList( DefaultQuadratureTraits< FieldType, dimension > (),
                            geometryType, quadKey )
    {
    }
 
    template <class FactoryTraits>
    inline IntegrationPointList ( const FactoryTraits traits,
                                  const GeometryType &geometryType,
                                  const GeometryType &elementGeometry,
                                  const typename FactoryTraits::QuadratureKeyType& quadKey )
    //: ipListPtr_( &QuadratureProviderType :: getQuadrature( traits, geometryType, elementGeometry, quadKey ), NoDelete() )
    : ipListPtr_( &QuadratureProviderType :: getQuadrature( geometryType, elementGeometry, quadKey ), NoDelete() )
    {
    }
 
    template <class QuadratureKeyType>
    inline IntegrationPointList ( const GeometryType &geometryType,
                                  const GeometryType &elementGeometry,
                                  const QuadratureKeyType& quadKey )
    : IntegrationPointList( DefaultQuadratureTraits< FieldType, dimension > (),
                            geometryType, elementGeometry, quadKey )
    {
    }
 
    inline IntegrationPointList ( const IntegrationPointListType& ipList )
    : ipListPtr_( &ipList, NoDelete() )
    {
    }
 
    inline IntegrationPointList ( const IntegrationPointListStorageType& ipListPtr )
    : ipListPtr_( ipListPtr )
    {
    }
 
    inline IntegrationPointList ( const IntegrationPointList &org ) = default;
 
    const QuadraturePointWrapperType operator[] ( unsigned int i ) const
    {
      return QuadraturePointWrapperType( *this, i );
    }
 
    const IntegrationPointListType &ipList () const
    {
      assert( ipListPtr_ );
      return *ipListPtr_;
    }
 
    int nop () const
    {
      return ipList().nop();
    }
 
    const CoordinateType &point ( size_t i ) const
    {
      return ipList().point( i );
    }
 
    size_t id () const
    {
      return ipList().id();
    }
 
    int order () const
    {
      return ipList().order();
    }
 
    GeometryType geometryType () const
    {
      return ipList().geometryType();
    }
 
    auto interpolationPoints( const int reqDim ) const
    {
      return ipList().interpolationPoints( reqDim );
    }
 
    bool isFaceInterpolationQuadrature( const size_t numShapeFunctions ) const
    {
      return ipList().isFaceInterpolationQuadrature( numShapeFunctions );
    }
 
    WeightReturnType weight( size_t i ) const
    {
      if constexpr ( isQuadrature )
      {
        return ipList().weight( i );
      }
      else
        return 1.0;
    }
  };
 
 
  template< class FieldImp, int dim,
            template< class, int > class QuadratureTraits = DefaultQuadratureTraits >
  using Quadrature = IntegrationPointList< FieldImp, dim, QuadratureTraits, true >;
 
 
 
  template <class Quadrature>
  using SelectQuadraturePointSetId = detail::SelectPointSetId< Quadrature, -Dune::QuadratureType::size >;
 
  } //end namespace Fem
 
} //end namespace Dune
 
#endif