lumpingquadrature.hh

#ifndef DUNE_FEM_LUMPING_QUADRATURE_HH
#define DUNE_FEM_LUMPING_QUADRATURE_HH
 
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/space/lagrange/lagrangepoints.hh>
 
namespace Dune {
 
namespace Fem {
 
template<class FieldImp, Dune::GeometryType::Id geometryId>
class LumpingQuadrature
  : public QuadratureImp<FieldImp, Dune::GeometryType(geometryId).dim()>
{
 public:
  typedef FieldImp FieldType;
  static constexpr auto dimension = Dune::GeometryType(geometryId).dim();
 
 private:
  typedef LumpingQuadrature<FieldType, geometryId> ThisType;
  typedef QuadratureImp<FieldType, dimension> BaseType;
 
  static const int maxOrder_ = 6;
  int order_ ;
 
 public:
  typedef typename BaseType::CoordinateType CoordinateType;
 
  LumpingQuadrature(const GeometryType& gt, int order, int id)
    : BaseType(id),
      order_( order )
  {
    assert( order_ <= maxOrder_ && order_ >= 1 );
 
    // create Lagrange points to define quadrature
    typedef LagrangePointListImplementation< FieldType, Dune::GeometryType(geometryId).id(), dimension, maxOrder_ > LP;
    LP pts( gt, order_, id );
 
    const auto &refElement = Dune::ReferenceElements< FieldType, dimension >::general( gt );
    //const std::size_t numCorners = refElement.size( dimension );
    const std::size_t numCorners = pts.nop();
    const FieldType weight = refElement.volume() / FieldType(numCorners);
    for( std::size_t i = 0 ; i < numCorners; ++i )
    {
      //auto pt = refElement.position( i, dimension );
      auto pt =  pts.point( i );
      this->addQuadraturePoint( pt, weight );
      //std::cout << "Point " << i << ": (" << pt << ") " << weight << std::endl;
    }
  }
 
  virtual GeometryType geometryType() const { return Dune::GeometryType(geometryId); }
  virtual int order () const { return order_; }
 
  static std::size_t maxOrder () { return maxOrder_; }
};
 
template<class FieldType, int dimension>
struct DefaultLumpingQuadratureTraits
{
  typedef QuadratureImp<FieldType, dimension> IntegrationPointListType;
 
  static constexpr Dune::GeometryType::Id simplexId = Dune::GeometryTypes::simplex(dimension);
  static constexpr Dune::GeometryType::Id cubeId    = Dune::GeometryTypes::cube(dimension);
  static constexpr Dune::GeometryType::Id prismId   = Dune::GeometryTypes::prism ;
  static constexpr Dune::GeometryType::Id pyramidId = Dune::GeometryTypes::pyramid;
 
  typedef LumpingQuadrature<FieldType, simplexId>  SimplexQuadratureType;
  typedef LumpingQuadrature<FieldType, cubeId   >  CubeQuadratureType;
  typedef LumpingQuadrature<FieldType, prismId  >  PrismQuadratureType;
  typedef LumpingQuadrature<FieldType, pyramidId>  PyramidQuadratureType;
 
  typedef SimplexQuadratureType LineQuadratureType;
  typedef SimplexQuadratureType PointQuadratureType;
 
  typedef int QuadratureKeyType;
};
 
// LumpingQuadrature uses CachingQuadrature with a different traits class for creating the quadratures.
template<class GridPartImp, int codim>
using CachingLumpingQuadrature = CachingQuadrature< GridPartImp, codim, DefaultLumpingQuadratureTraits >;
 
} // Fem
 
} // Dune
 
 
#endif // DUNE_FEM_LUMPING_QUADRATURE_HH