4 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
5 #define DUNE_GEOMETRY_QUADRATURERULES_HH
14 #include <dune/common/fvector.hh>
15 #include <dune/common/exceptions.hh>
16 #include <dune/common/stdstreams.hh>
17 #include <dune/common/stdthread.hh>
18 #include <dune/common/visibility.hh>
41 template<
typename ct,
int dim>
51 typedef Dune::FieldVector<ct,dim>
Vector;
79 namespace QuadratureType {
124 template<
typename ct,
int dim>
157 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
171 template<
typename ctype,
int dim>
184 typedef std::vector<std::pair<std::once_flag, QuadratureRule> >
185 QuadratureOrderVector;
188 static void initQuadratureOrderVector(QuadratureOrderVector *qov,
194 *qov = QuadratureOrderVector(1);
199 typedef std::vector<std::pair<std::once_flag, QuadratureOrderVector> >
203 static void initGeometryTypeVector(GeometryTypeVector *gtv)
211 assert(t.
dim()==dim);
213 DUNE_ASSERT_CALL_ONCE();
215 static std::vector<std::pair<
220 auto & quadratureTypeLevel = quadratureCache[qt];
221 std::call_once(quadratureTypeLevel.first, initGeometryTypeVector,
222 &quadratureTypeLevel.second);
224 auto & geometryTypeLevel =
226 std::call_once(geometryTypeLevel.first, initQuadratureOrderVector,
227 &geometryTypeLevel.second, qt, t);
230 auto & quadratureOrderLevel = geometryTypeLevel.second[dim == 0 ? 0 : p];
231 std::call_once(quadratureOrderLevel.first, initQuadratureRule,
232 &quadratureOrderLevel.second, qt, t, p);
234 return quadratureOrderLevel.second;
256 return instance()._rule(t,p,qt);
259 DUNE_NO_DEPRECATED_BEGIN
264 return instance()._rule(gt,p,qt);
266 DUNE_NO_DEPRECATED_END
271 #include "quadraturerules/pointquadrature.hh"
276 template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
278 template<
typename ct>
281 std::vector< FieldVector<ct, 1> > & _points,
282 std::vector< ct > & _weight,
283 int & delivered_order);
285 template<
typename ct>
288 std::vector< FieldVector<ct, 1> > & _points,
289 std::vector< ct > & _weight,
290 int & delivered_order);
294 template<
typename ct>
310 std::vector< FieldVector<ct, dim> > _points;
311 std::vector< ct > _weight;
316 assert(_points.size() == _weight.size());
317 for (
size_t i = 0; i < _points.size(); i++)
322 extern template GaussQuadratureRule1D<float>::GaussQuadratureRule1D(
int);
323 extern template GaussQuadratureRule1D<double>::GaussQuadratureRule1D(
int);
327 #define DUNE_INCLUDING_IMPLEMENTATION
328 #include "quadraturerules/gauss_imp.hh"
333 template<
typename ct,
334 bool fundamental = std::numeric_limits<ct>::is_specialized>
336 template<
typename ct>
339 std::vector< FieldVector<ct, 1> > & _points,
340 std::vector< ct > & _weight,
341 int & delivered_order);
343 template<
typename ct>
346 std::vector< FieldVector<ct, 1> > & _points,
347 std::vector< ct > & _weight,
348 int & delivered_order);
354 template<
typename ct>
363 enum { highest_order=61 };
372 std::vector< FieldVector<ct, dim> > _points;
373 std::vector< ct > _weight;
378 (p, _points, _weight, deliveredOrder_);
379 this->delivered_order = deliveredOrder_;
380 assert(_points.size() == _weight.size());
381 for (
size_t i = 0; i < _points.size(); i++)
387 extern template Jacobi1QuadratureRule1D<float>::Jacobi1QuadratureRule1D(
int);
388 extern template Jacobi1QuadratureRule1D<double>::Jacobi1QuadratureRule1D(
int);
393 #define DUNE_INCLUDING_IMPLEMENTATION
394 #include "quadraturerules/jacobi_1_0_imp.hh"
399 template<
typename ct,
400 bool fundamental = std::numeric_limits<ct>::is_specialized>
402 template<
typename ct>
405 std::vector< FieldVector<ct, 1> > & _points,
406 std::vector< ct > & _weight,
407 int & delivered_order);
409 template<
typename ct>
412 std::vector< FieldVector<ct, 1> > & _points,
413 std::vector< ct > & _weight,
414 int & delivered_order);
420 template<
typename ct>
429 enum { highest_order=61 };
438 std::vector< FieldVector<ct, dim> > _points;
439 std::vector< ct > _weight;
444 (p, _points, _weight, deliveredOrder_);
446 this->delivered_order = deliveredOrder_;
447 assert(_points.size() == _weight.size());
448 for (
size_t i = 0; i < _points.size(); i++)
454 extern template Jacobi2QuadratureRule1D<float>::Jacobi2QuadratureRule1D(
int);
455 extern template Jacobi2QuadratureRule1D<double>::Jacobi2QuadratureRule1D(
int);
460 #define DUNE_INCLUDING_IMPLEMENTATION
461 #include "quadraturerules/jacobi_2_0_imp.hh"
466 template<
typename ct,
467 bool fundamental = std::numeric_limits<ct>::is_specialized>
469 template<
typename ct>
472 std::vector< FieldVector<ct, 1> > & _points,
473 std::vector< ct > & _weight,
474 int & delivered_order);
476 template<
typename ct>
479 std::vector< FieldVector<ct, 1> > & _points,
480 std::vector< ct > & _weight,
481 int & delivered_order);
487 template<
typename ct>
496 enum { highest_order=31 };
505 std::vector< FieldVector<ct, dim> > _points;
506 std::vector< ct > _weight;
511 (p, _points, _weight, deliveredOrder_);
513 this->delivered_order = deliveredOrder_;
514 assert(_points.size() == _weight.size());
515 for (
size_t i = 0; i < _points.size(); i++)
521 extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(
int);
522 extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(
int);
527 #define DUNE_INCLUDING_IMPLEMENTATION
528 #include "quadraturerules/gausslobatto_imp.hh"
530 #include "quadraturerules/tensorproductquadrature.hh"
532 #include "quadraturerules/simplexquadrature.hh"
550 enum { highest_order=2 };
584 W[m][0] = 0.16666666666666666 / 2.0;
585 W[m][1] = 0.16666666666666666 / 2.0;
586 W[m][2] = 0.16666666666666666 / 2.0;
587 W[m][3] = 0.16666666666666666 / 2.0;
588 W[m][4] = 0.16666666666666666 / 2.0;
589 W[m][5] = 0.16666666666666666 / 2.0;
596 G[m][0][0] =0.66666666666666666 ;
597 G[m][0][1] =0.16666666666666666 ;
598 G[m][0][2] =0.211324865405187 ;
600 G[m][1][0] = 0.16666666666666666;
601 G[m][1][1] =0.66666666666666666 ;
602 G[m][1][2] = 0.211324865405187;
604 G[m][2][0] = 0.16666666666666666;
605 G[m][2][1] = 0.16666666666666666;
606 G[m][2][2] = 0.211324865405187;
608 G[m][3][0] = 0.66666666666666666;
609 G[m][3][1] = 0.16666666666666666;
610 G[m][3][2] = 0.788675134594813;
612 G[m][4][0] = 0.16666666666666666;
613 G[m][4][1] = 0.66666666666666666;
614 G[m][4][2] = 0.788675134594813;
616 G[m][5][0] = 0.16666666666666666;
617 G[m][5][1] = 0.16666666666666666;
618 G[m][5][2] = 0.788675134594813;
620 W[m][0] = 0.16666666666666666 / 2.0;
621 W[m][1] = 0.16666666666666666 / 2.0;
622 W[m][2] = 0.16666666666666666 / 2.0;
623 W[m][3] = 0.16666666666666666 / 2.0;
624 W[m][4] = 0.16666666666666666 / 2.0;
625 W[m][5] = 0.16666666666666666 / 2.0;
632 FieldVector<double, 3>
point(
int m,
int i)
650 FieldVector<double, 3> G[MAXP+1][MAXP];
652 double W[MAXP+1][MAXP];
676 template<
typename ct,
int dim>
682 template<
typename ct>
691 enum { highest_order = 2 };
700 for(
int i=0; i<m; ++i)
702 FieldVector<ct,3> local;
703 for (
int k=0; k<d; k++)
719 template<
typename ctype,
int dim>
725 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
729 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
733 template<
typename ctype>
742 return std::numeric_limits<int>::max();
744 DUNE_THROW(Exception,
"Unknown GeometryType");
750 return PointQuadratureRule<ctype>();
752 DUNE_THROW(Exception,
"Unknown GeometryType");
756 template<
typename ctype>
775 return JacobiNQuadratureRule1D<ctype>::maxOrder();
777 DUNE_THROW(Exception,
"Unknown QuadratureType");
780 DUNE_THROW(Exception,
"Unknown GeometryType");
796 return JacobiNQuadratureRule1D<ctype>(p);
798 DUNE_THROW(Exception,
"Unknown QuadratureType");
801 DUNE_THROW(Exception,
"Unknown GeometryType");
805 template<
typename ctype>
813 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
816 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
823 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
825 return SimplexQuadratureRule<ctype,dim>(p);
827 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
831 template<
typename ctype>
839 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
842 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
853 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
855 return SimplexQuadratureRule<ctype,dim>(p);
863 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
A unique label for each type of element that can occur in a grid.
Helper classes to provide indices for geometrytypes for use in a vector.
constexpr GeometryType line
GeometryType representing a line.
Definition: type.hh:803
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:833
Definition: affinegeometry.hh:19
Enum
Definition: quadraturerules.hh:80
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:115
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition: quadraturerules.hh:102
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition: quadraturerules.hh:95
@ size
Definition: quadraturerules.hh:117
@ GaussLobatto
Definition: quadraturerules.hh:116
@ GaussLegendre
Gauss-Legendre rules (default)
Definition: quadraturerules.hh:88
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:34
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:42
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition: quadraturerules.hh:51
ct Field
Number type used for coordinates and quadrature weights.
Definition: quadraturerules.hh:48
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:60
@ dimension
Definition: quadraturerules.hh:45
ct weight_
Definition: quadraturerules.hh:73
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:66
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:54
FieldVector< ct, dim > local
Definition: quadraturerules.hh:72
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:126
virtual ~QuadratureRule()
Definition: quadraturerules.hh:153
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:152
int delivered_order
Definition: quadraturerules.hh:161
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:140
GeometryType geometry_type
Definition: quadraturerules.hh:160
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:146
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:133
virtual int order() const
return order
Definition: quadraturerules.hh:149
@ d
Definition: quadraturerules.hh:143
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid
Definition: quadraturerules.hh:137
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:157
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition: quadraturerules.hh:720
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:172
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition: quadraturerules.hh:247
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:254
static DUNE_NO_DEPRECATED_BEGIN const QuadratureRule & rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:261
Definition: quadraturerules.hh:277
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
Gauss quadrature rule in 1D.
Definition: quadraturerules.hh:297
@ dim
Definition: quadraturerules.hh:300
@ highest_order
Definition: quadraturerules.hh:301
~GaussQuadratureRule1D()
Definition: quadraturerules.hh:303
Definition: quadraturerules.hh:335
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
Jacobi-Gauss quadrature for alpha=1, beta=0.
Definition: quadraturerules.hh:357
Definition: quadraturerules.hh:401
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:423
Definition: quadraturerules.hh:468
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
static void init(int p, std::vector< FieldVector< ct, 1 > > &_points, std::vector< ct > &_weight, int &delivered_order)
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:490
Definition: quadraturerules.hh:542
Definition: quadraturerules.hh:547
double weight(int m, int i)
Definition: quadraturerules.hh:638
int order(int m)
Definition: quadraturerules.hh:644
PrismQuadraturePoints()
initialize quadrature points on the interval for all orders
Definition: quadraturerules.hh:553
FieldVector< double, 3 > point(int m, int i)
Definition: quadraturerules.hh:632
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:661
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:662
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:669
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:670
Quadrature rules for prisms.
Definition: quadraturerules.hh:677
Quadrature rules for prisms.
Definition: quadraturerules.hh:684
Definition: quadraturerules.hh:734
Definition: quadraturerules.hh:757
Definition: quadraturerules.hh:806
Definition: quadraturerules.hh:832
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
constexpr bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:619
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:589
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:644
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:285
constexpr bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:594
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:649
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:629
static constexpr std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition: typeindex.hh:56
static constexpr std::size_t index(const GeometryType >)
Compute the index for the given geometry type within its dimension.
Definition: typeindex.hh:68