elementquadrature.hh

#ifndef DUNE_FEM_ELEMENTQUADRATURE_HH
#define DUNE_FEM_ELEMENTQUADRATURE_HH
 
#include <dune/geometry/quadraturerules.hh>
 
#include <dune/fem/common/utility.hh>
 
#include "quadrature.hh"
#include "elementpointlist.hh"
 
namespace Dune
{
 
  namespace Fem
  {
 
    template< typename GridPartImp, class IntegrationPointList >
    class Agglomeration;
 
 
    template< typename GridPartImp, int codim, template <class, int> class QuadratureTraits = DefaultQuadratureTraits >
    class ElementQuadrature;
 
 
 
    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< typename GridPartImp, template< class, int > class QuadratureTraits >
    class ElementQuadrature< GridPartImp, 0, QuadratureTraits >
    : public ElementIntegrationPointList< GridPartImp, 0, ElementQuadratureTraits< GridPartImp, 0, QuadratureTraits > >
    {
      typedef ElementQuadrature< GridPartImp, 0, QuadratureTraits > ThisType;
 
    public:
      typedef ElementQuadratureTraits< GridPartImp, 0, QuadratureTraits > IntegrationTraits;
      typedef ElementIntegrationPointList< GridPartImp, 0, IntegrationTraits > BaseType;
 
      typedef GridPartImp GridPartType;
 
      enum { codimension = 0 };
 
      enum { dimension = GridPartType :: dimension };
 
      typedef typename GridPartType :: ctype RealType;
 
      typedef typename IntegrationTraits :: CoordinateType CoordinateType;
 
      typedef typename BaseType :: QuadratureKeyType QuadratureKeyType;
 
      typedef QuadraturePointWrapper< ThisType > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< ThisType > IteratorType;
 
      // for compatibility
      typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
 
      typedef typename BaseType :: IntegrationPointListType  IntegrationPointListType;
 
    protected:
      using BaseType::quadImp;
 
      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 );
        }
      }
 
    public:
      using BaseType::nop;
 
      ElementQuadrature( const EntityType &entity, const QuadratureKeyType& quadKey, const bool checkGeomType = true )
      : BaseType( createQuadrature( entity, quadKey, checkGeomType ) )
      {}
 
      ElementQuadrature( const GeometryType &type, const QuadratureKeyType& quadKey )
      : BaseType( type, quadKey )
      {
        // when type is none then entity has to be passed in order to create
        // the sub triangulation etc.
        // assert( ! type.isNone() );
      }
 
      ElementQuadrature( const ThisType &org )
      : BaseType( org )
      {}
 
      QuadraturePointWrapperType operator[] ( std::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 RealType &weight( size_t i ) const
      {
        return quadImp().weight( i );
      }
    };
 
 
 
    template< class GridPartImp, template< class, int > class QuadratureTraits >
    class ElementQuadrature< GridPartImp, 1, QuadratureTraits >
    : public ElementIntegrationPointList
      < GridPartImp, 1, ElementQuadratureTraits< GridPartImp, 1, QuadratureTraits > >
    {
    public:
      typedef GridPartImp GridPartType;
 
      enum { codimension = 1 };
 
    private:
      typedef ElementQuadratureTraits< GridPartType, codimension, QuadratureTraits > IntegrationTraits;
 
      typedef ElementQuadrature< GridPartType, codimension, QuadratureTraits > ThisType;
      typedef ElementIntegrationPointList< GridPartType, 1, IntegrationTraits >
        BaseType;
 
    protected:
      using BaseType :: quadImp;
 
    public:
      enum { dimension = GridPartType :: dimension };
 
      typedef typename GridPartType :: ctype RealType;
 
      typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType :: Intersection IntersectionType;
 
      typedef typename IntegrationTraits :: CoordinateType CoordinateType;
 
      typedef typename BaseType :: QuadratureKeyType QuadratureKeyType;
 
      typedef QuadraturePointWrapper< ThisType > QuadraturePointWrapperType;
      typedef QuadraturePointIterator< ThisType > IteratorType;
 
      typedef typename IntegrationTraits :: IntegrationPointListType :: CoordinateType
        LocalCoordinateType;
 
      typedef ThisType NonConformingQuadratureType;
 
    public:
      using BaseType::nop;
 
      ElementQuadrature ( const GridPartType &gridPart,
                          const IntersectionType &intersection,
                          const QuadratureKeyType& quadKey,
                          typename BaseType :: Side side )
      : BaseType( gridPart, intersection, quadKey, side )
      {}
 
      ElementQuadrature( const ElementQuadrature &org )
      : BaseType( org )
      {
      }
 
      QuadraturePointWrapperType operator[] ( std::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 RealType &weight( size_t i ) const
      {
        return quadImp().weight( i );
      }
    };
 
    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