cachingquadrature.hh

#ifndef DUNE_FEM_CACHINGQUADRATURE_HH
#define DUNE_FEM_CACHINGQUADRATURE_HH
 
#if __GNUC__ >= 13
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdangling-reference"
#endif
 
//- Dune includes
#include <dune/common/math.hh>
 
//- Local includes
#include <dune/fem/quadrature/elementquadrature.hh>
#include <dune/fem/quadrature/caching/twistutility.hh>
#include <dune/fem/quadrature/caching/pointmapper.hh>
#include <dune/fem/quadrature/caching/cacheprovider.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    class CachingInterface
    {
    protected:
      // do not create instances of this class
      CachingInterface ()
      {
      }
 
    public:
      static constexpr bool twisted () { return false; }
 
      inline int twistId () const { return 0; }
 
      inline size_t cachingPoint( const size_t quadraturePoint ) const
      {
        DUNE_THROW( NotImplemented,
                    "CachingInterface :: cachingPoint must be overloaded!" );
      }
 
      inline size_t interpolationPoint( const size_t quadraturePoint ) const
      {
        DUNE_THROW( NotImplemented,
                    "CachingInterface :: interpolationPoint must be overloaded!" );
      }
 
      inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
      {
        DUNE_THROW( NotImplemented,
                    "CachingInterface :: isInterpolationQuadrature must be overloaded!" );
      }
    };
 
 
 
    template< class GridPartImp, int codim, class IntegrationTraits, bool isQuadrature >
    class CachingQuadratureImpl;
 
 
    template< class GridPartImp, class IntegrationTraits, bool isQuadrature >
    class CachingQuadratureImpl< GridPartImp, 0, IntegrationTraits, isQuadrature >
    : public ElementPointListBase< GridPartImp, 0, IntegrationTraits >,
      public CachingInterface
    {
      typedef CachingQuadratureImpl< GridPartImp, 0, IntegrationTraits, isQuadrature > This;
      typedef ElementPointListBase< GridPartImp, 0, IntegrationTraits > Base;
 
    public:
      typedef typename Base::GridPartType GridPartType;
      typedef typename Base::EntityType   EntityType;
      static const int codimension = Base::codimension;
 
      typedef typename Base::CoordinateType CoordinateType;
 
      typedef QuadraturePointWrapper< This > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< This > IteratorType;
 
      static const int pointSetId = SelectQuadraturePointSetId<
          typename IntegrationTraits::IntegrationPointListType::Traits > :: value;
 
    protected:
      using Base::quadImp;
 
    public:
      using CachingInterface::twisted;
      using CachingInterface::twistId;
      using Base::localPoint;
      using Base::nop;
 
      template <class QuadratureKeyType>
      CachingQuadratureImpl( const EntityType& entity, const QuadratureKeyType& quadKey )
      : Base( entity.type(), quadKey )
      {
        CacheProvider< GridPartType, codimension >::registerQuadrature( quadImp() );
      }
 
      template <class QuadratureKeyType>
      CachingQuadratureImpl( const GeometryType &geometry, const QuadratureKeyType& quadKey )
      : Base( geometry, quadKey )
      {
        CacheProvider< GridPartType, codimension >::registerQuadrature( quadImp() );
      }
 
      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 ( std::size_t i ) const
      {
        return quadImp().weight( i );
      }
 
      inline size_t cachingPoint( const size_t quadraturePoint ) const
      {
        return quadraturePoint;
      }
 
      inline size_t interpolationPoint( const size_t quadraturePoint ) const
      {
        return quadraturePoint;
      }
 
      inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
      {
        // if pointSetId is not negative then we have an interpolation
        // quadrature if the number of point are equal to number of shape functions
        return (pointSetId >= 0) ? (nop() == numShapeFunctions) : false;
      }
    };
 
 
 
    template< typename GridPartImp, class IntegrationTraits, bool isQuadrature >
    class CachingQuadratureImpl< GridPartImp, 1, IntegrationTraits, isQuadrature >
    : public ElementPointListBase< GridPartImp, 1, IntegrationTraits >,
      public CachingInterface
    {
      typedef CachingQuadratureImpl< GridPartImp, 1, IntegrationTraits, isQuadrature > This;
      typedef ElementPointListBase< GridPartImp, 1, IntegrationTraits > Base;
 
    public:
      typedef GridPartImp GridPartType;
 
      typedef typename Base::RealType RealType;
      static const int dimension = Base::dimension;
      static const int codimension = Base::codimension;
 
      typedef typename Base::CoordinateType CoordinateType;
 
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
 
      typedef QuadraturePointWrapper< This > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< This > IteratorType;
 
      typedef ElementQuadratureImpl< GridPartType, codimension, IntegrationTraits, isQuadrature >
        NonConformingQuadratureType;
 
 
      // for compatibility
      typedef typename GridPartType::TwistUtilityType  TwistUtilityType;
      typedef IntersectionIteratorType IntersectionIterator;
 
    private:
      static const int quadPointSetId =
        SelectQuadraturePointSetId< typename IntegrationTraits::IntegrationPointListType::Traits > :: value;
 
    public:
      // Note: we also exclude GaussLegendre(0) here, because on faces it is not
      //       an interpolation rule
      static const int pointSetId = (quadPointSetId > 0) ? quadPointSetId :
                  SelectQuadraturePointSetId< void > :: value; // default value
 
    protected:
      typedef typename CachingTraits< RealType, dimension >::MapperPairType  MapperPairType;
      typedef typename CachingTraits< RealType, dimension >::PointVectorType PointVectorType;
 
      typedef CacheProvider< GridPartType, codimension >            CacheProviderType;
      typedef PointProvider< RealType, dimension, codimension>  PointProviderType;
 
      using Base::quadImp;
 
    public:
      using Base::localFaceIndex;
      using Base::elementGeometry;
      using Base::nop;
 
      template <class QuadratureKeyType>
      CachingQuadratureImpl ( const GridPartType &gridPart,
                              const IntersectionType &intersection,
                              const QuadratureKeyType& quadKey, const typename Base :: Side side )
        : Base( getPointList( intersection, quadKey, side ) ),
          side_(side),
          twist_( getTwist( gridPart, intersection, side ) ),
          mapper_( CacheProviderType::getMapper( quadImp(), elementGeometry(), localFaceIndex(), twist_) ),
          points_( PointProviderType::getPoints( quadImp().ipList().id(), elementGeometry() ) ),
          intersection_(intersection)
      {
      }
 
      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() ); }
 
      typename Base :: Side side() const { return side_; }
      bool isInside() const { return side_ == Base::INSIDE; }
 
      const CoordinateType &point ( const size_t i ) const
      {
        return points_[ cachingPoint( i ) ];
      }
 
      auto weight ( std::size_t i ) const
      {
        return quadImp().weight( i );
      }
 
      const IntersectionType &intersection() const
      {
        return intersection_;
      }
 
      static constexpr bool twisted() { return true; }
 
      inline int twistId () const { return twist_ + 4; }
 
      inline size_t cachingPoint( const size_t quadraturePoint ) const
      {
        assert( quadraturePoint < (size_t)nop() );
        return mapper_.first[ quadraturePoint ];
      }
 
      inline size_t interpolationPoint( const size_t quadraturePoint ) const
      {
        assert( quadraturePoint < mapper_.second.size() );
        return mapper_.second[ quadraturePoint ];
      }
 
      inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
      {
        // if pointSetId is not negative then we have an interpolation
        // quadrature if the number of point are equal to number of shape functions
        return (pointSetId < 0) ? false :
          quadImp().ipList().isFaceInterpolationQuadrature( numShapeFunctions );
      }
 
      // return local caching point
      // for debugging issues only
      size_t localCachingPoint ( const size_t i ) const
      {
        const auto& mapper = mapper_.first;
 
        assert( i < (size_t)nop() );
 
        assert( mapper[ i ] >= 0 );
        int faceIndex = localFaceIndex();
        unsigned int point = mapper[ i ] - faceIndex * mapper.size();
        assert( point < nop() );
 
        return point;
      }
 
    protected:
      template <class QuadratureKeyType>
      Base getPointList ( const IntersectionType &intersection,
                          const QuadratureKeyType& key,
                          const typename Base :: Side side )
      {
        switch( side )
        {
          case Base :: INSIDE:
            return Base( TwistUtilityType::elementGeometry( intersection, true ),
                         intersection.indexInInside(), key );
 
          case Base :: OUTSIDE:
            return Base( TwistUtilityType::elementGeometry( intersection, false ),
                         intersection.indexInOutside(), key );
 
          default:
            DUNE_THROW( InvalidStateException, "ElementIntegrationPointList: side must either be INSIDE or OUTSIDE." );
        }
      }
 
      int getTwist ( const GridPartType &gridPart,
                     const IntersectionType &intersection,
                     const typename Base :: Side side )
      {
        switch( side )
        {
          case Base :: INSIDE:
            return TwistUtilityType::twistInSelf( gridPart.grid(), intersection );
 
          case Base :: OUTSIDE:
            return TwistUtilityType::twistInNeighbor( gridPart.grid(), intersection );
 
          default:
            DUNE_THROW( InvalidStateException, "ElementIntegrationPointList: side must either be INSIDE or OUTSIDE." );
        }
      }
 
    private:
      const typename Base :: Side side_;
      const int twist_;
      const MapperPairType &mapper_;
      const PointVectorType &points_;
      const IntersectionType &intersection_;
    };
 
 
    template< typename GridPart, int codim, class IntegrationTraits >
    using CachingPointList = CachingQuadratureImpl< GridPart, codim, IntegrationTraits, false>;
 
    template< typename GridPart, int codim, template <class, int> class QuadratureTraits = DefaultQuadratureTraits >
    using CachingQuadrature = CachingQuadratureImpl< GridPart, codim, ElementQuadratureTraits< GridPart, codim, QuadratureTraits >, true>;
 
    template<class GridPart, class Entity>
    static inline auto cachingQuadrature(const GridPart& gridPart, const Entity& entity, unsigned quadOrder)
    {
      using Quadrature = Dune::Fem::CachingQuadrature<GridPart, 0>;
      return Quadrature(entity, quadOrder);
    }
 
  } //namespace Fem
 
} //namespace Dune
 
#if __GNUC__ >= 13
#pragma GCC diagnostic pop
#endif
 
#endif // #ifndef DUNE_FEM_CACHINGQUADRATURE_HH