4 #ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
5 #define DUNE_GEOMETRY_QUADRATURERULES_HH
12 #include <dune/common/fvector.hh>
13 #include <dune/common/exceptions.hh>
14 #include <dune/common/stdstreams.hh>
15 #include <dune/common/visibility.hh>
16 #include <dune/common/deprecated.hh>
38 template<
typename ct,
int dim>
42 enum {
d=dim } DUNE_DEPRECATED_MSG(
"Use 'dimension' instead");
43 typedef ct
CoordType DUNE_DEPRECATED_MSG(
"Use type 'Field' instead");
47 typedef Dune::FieldVector<ct,dim>
Vector;
75 namespace QuadratureType {
92 template<
typename ct,
int dim>
121 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
135 template<
typename ctype,
int dim>
141 typedef std::pair<GeometryType,int> QuadratureRuleKey;
149 static std::map<QuadratureRuleKey, QuadratureRule> _quadratureMap;
150 QuadratureRuleKey key(t,p);
151 if (_quadratureMap.find(key) == _quadratureMap.end()) {
159 _quadratureMap.insert(std::make_pair(key, rule));
161 return _quadratureMap.find(key)->second;
175 return instance()._rule(t,p,qt);
181 return instance()._rule(gt,p,qt);
187 #include "quadraturerules/pointquadrature.hh"
192 template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
194 template<
typename ct>
196 static void init(
int p,
197 std::vector< FieldVector<ct, 1> > & _points,
198 std::vector< ct > & _weight,
199 int & delivered_order);
201 template<
typename ct>
203 static void init(
int p,
204 std::vector< FieldVector<ct, 1> > & _points,
205 std::vector< ct > & _weight,
206 int & delivered_order);
210 template<
typename ct>
226 std::vector< FieldVector<ct, dim> > _points;
227 std::vector< ct > _weight;
232 assert(_points.size() == _weight.size());
233 for (
size_t i = 0; i < _points.size(); i++)
238 extern template GaussQuadratureRule1D<float>::GaussQuadratureRule1D(
int);
239 extern template GaussQuadratureRule1D<double>::GaussQuadratureRule1D(
int);
243 #define DUNE_INCLUDING_IMPLEMENTATION
244 #include "quadraturerules/gauss_imp.hh"
249 template<
typename ct,
250 bool fundamental = std::numeric_limits<ct>::is_specialized>
252 template<
typename ct>
254 static void init(
int p,
255 std::vector< FieldVector<ct, 1> > & _points,
256 std::vector< ct > & _weight,
257 int & delivered_order);
259 template<
typename ct>
261 static void init(
int p,
262 std::vector< FieldVector<ct, 1> > & _points,
263 std::vector< ct > & _weight,
264 int & delivered_order);
270 template<
typename ct>
279 enum { highest_order=61 };
288 std::vector< FieldVector<ct, dim> > _points;
289 std::vector< ct > _weight;
294 (p, _points, _weight, delivered_order);
295 this->delivered_order = delivered_order;
296 assert(_points.size() == _weight.size());
297 for (
size_t i = 0; i < _points.size(); i++)
303 extern template Jacobi1QuadratureRule1D<float>::Jacobi1QuadratureRule1D(
int);
304 extern template Jacobi1QuadratureRule1D<double>::Jacobi1QuadratureRule1D(
int);
309 #define DUNE_INCLUDING_IMPLEMENTATION
310 #include "quadraturerules/jacobi_1_0_imp.hh"
315 template<
typename ct,
316 bool fundamental = std::numeric_limits<ct>::is_specialized>
318 template<
typename ct>
320 static void init(
int p,
321 std::vector< FieldVector<ct, 1> > & _points,
322 std::vector< ct > & _weight,
323 int & delivered_order);
325 template<
typename ct>
327 static void init(
int p,
328 std::vector< FieldVector<ct, 1> > & _points,
329 std::vector< ct > & _weight,
330 int & delivered_order);
336 template<
typename ct>
345 enum { highest_order=61 };
354 std::vector< FieldVector<ct, dim> > _points;
355 std::vector< ct > _weight;
360 (p, _points, _weight, delivered_order);
362 this->delivered_order = delivered_order;
363 assert(_points.size() == _weight.size());
364 for (
size_t i = 0; i < _points.size(); i++)
370 extern template Jacobi2QuadratureRule1D<float>::Jacobi2QuadratureRule1D(
int);
371 extern template Jacobi2QuadratureRule1D<double>::Jacobi2QuadratureRule1D(
int);
376 #define DUNE_INCLUDING_IMPLEMENTATION
377 #include "quadraturerules/jacobi_2_0_imp.hh"
382 template<
typename ct,
383 bool fundamental = std::numeric_limits<ct>::is_specialized>
385 template<
typename ct>
387 static void init(
int p,
388 std::vector< FieldVector<ct, 1> > & _points,
389 std::vector< ct > & _weight,
390 int & delivered_order);
392 template<
typename ct>
394 static void init(
int p,
395 std::vector< FieldVector<ct, 1> > & _points,
396 std::vector< ct > & _weight,
397 int & delivered_order);
403 template<
typename ct>
412 enum { highest_order=61 };
421 std::vector< FieldVector<ct, dim> > _points;
422 std::vector< ct > _weight;
427 (p, _points, _weight, delivered_order);
429 this->delivered_order = delivered_order;
430 assert(_points.size() == _weight.size());
431 for (
size_t i = 0; i < _points.size(); i++)
437 extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(
int);
438 extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(
int);
443 #define DUNE_INCLUDING_IMPLEMENTATION
444 #include "quadraturerules/gausslobatto_imp.hh"
446 #include "quadraturerules/tensorproductquadrature.hh"
448 #include "quadraturerules/simplexquadrature.hh"
466 enum { highest_order=2 };
500 W[m][0] = 0.16666666666666666 / 2.0;
501 W[m][1] = 0.16666666666666666 / 2.0;
502 W[m][2] = 0.16666666666666666 / 2.0;
503 W[m][3] = 0.16666666666666666 / 2.0;
504 W[m][4] = 0.16666666666666666 / 2.0;
505 W[m][5] = 0.16666666666666666 / 2.0;
512 G[m][0][0] =0.66666666666666666 ;
513 G[m][0][1] =0.16666666666666666 ;
514 G[m][0][2] =0.211324865405187 ;
516 G[m][1][0] = 0.16666666666666666;
517 G[m][1][1] =0.66666666666666666 ;
518 G[m][1][2] = 0.211324865405187;
520 G[m][2][0] = 0.16666666666666666;
521 G[m][2][1] = 0.16666666666666666;
522 G[m][2][2] = 0.211324865405187;
524 G[m][3][0] = 0.66666666666666666;
525 G[m][3][1] = 0.16666666666666666;
526 G[m][3][2] = 0.788675134594813;
528 G[m][4][0] = 0.16666666666666666;
529 G[m][4][1] = 0.66666666666666666;
530 G[m][4][2] = 0.788675134594813;
532 G[m][5][0] = 0.16666666666666666;
533 G[m][5][1] = 0.16666666666666666;
534 G[m][5][2] = 0.788675134594813;
536 W[m][0] = 0.16666666666666666 / 2.0;
537 W[m][1] = 0.16666666666666666 / 2.0;
538 W[m][2] = 0.16666666666666666 / 2.0;
539 W[m][3] = 0.16666666666666666 / 2.0;
540 W[m][4] = 0.16666666666666666 / 2.0;
541 W[m][5] = 0.16666666666666666 / 2.0;
548 FieldVector<double, 3>
point(
int m,
int i)
566 FieldVector<double, 3> G[MAXP+1][MAXP];
568 double W[MAXP+1][MAXP];
592 template<
typename ct,
int dim>
598 template<
typename ct>
607 enum { highest_order = 2 };
616 for(
int i=0; i<m; ++i)
618 FieldVector<ct,3> local;
619 for (
int k=0; k<d; k++)
635 template<
typename ctype,
int dim>
636 class QuadratureRuleFactory {
641 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
645 template<
typename ctype>
654 return PointQuadratureRule<ctype>();
656 DUNE_THROW(Exception,
"Unknown GeometryType");
660 template<
typename ctype>
679 DUNE_THROW(Exception,
"Unknown QuadratureType");
682 DUNE_THROW(Exception,
"Unknown GeometryType");
686 template<
typename ctype>
693 if (t.
isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
695 return SimplexQuadratureRule<ctype,dim>(p);
697 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
701 template<
typename ctype>
708 if (t.
isSimplex() && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
710 return SimplexQuadratureRule<ctype,dim>(p);
712 if (t.
isPrism() && p <= PrismQuadratureRule<ctype,dim>::highest_order)
716 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
722 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH
Definition: quadraturerules.hh:78
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:116
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:110
Definition: quadraturerules.hh:193
Quadrature rules for prisms.
Definition: quadraturerules.hh:593
double weight(int m, int i)
Definition: quadraturerules.hh:554
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:31
Definition: quadraturerules.hh:45
Definition: quadraturerules.hh:77
A unique label for each type of element that can occur in a grid.
Definition: quadraturerules.hh:83
Dune::FieldVector< ct, dim > Vector
Definition: quadraturerules.hh:47
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType...
Definition: quadraturerules.hh:130
~Jacobi2QuadratureRule1D()
Definition: quadraturerules.hh:347
Definition: quadraturerules.hh:384
bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:307
Definition: quadraturerules.hh:458
bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:267
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:104
Definition: quadraturerules.hh:687
int delivered_order
Definition: quadraturerules.hh:125
~PrismQuadratureRule()
Definition: quadraturerules.hh:609
Definition: quadraturerules.hh:217
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:404
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:98
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:173
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:121
static 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:178
Definition: quadraturerules.hh:107
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:56
Definition: quadraturerules.hh:42
Enum
Definition: quadraturerules.hh:76
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:62
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:586
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:585
Definition: quadraturerules.hh:81
Definition: quadraturerules.hh:317
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:578
bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:272
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:337
unsigned int id() const
Return the topology id the type.
Definition: type.hh:327
Definition: quadraturerules.hh:216
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:93
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid.
Definition: quadraturerules.hh:101
Definition: quadraturerules.hh:251
bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:297
~GaussLobattoQuadratureRule1D()
Definition: quadraturerules.hh:414
Definition: quadraturerules.hh:661
GeometryType geometry_type
Definition: quadraturerules.hh:124
Definition: quadraturerules.hh:85
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:577
virtual ~QuadratureRule()
Definition: quadraturerules.hh:117
ct CoordType
Definition: quadraturerules.hh:43
ct Field
Definition: quadraturerules.hh:46
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:136
Quadrature rules for prisms.
Definition: quadraturerules.hh:599
FieldVector< double, 3 > point(int m, int i)
Definition: quadraturerules.hh:548
virtual int order() const
return order
Definition: quadraturerules.hh:113
FieldVector< ct, dim > local
Definition: quadraturerules.hh:68
PrismQuadraturePoints()
initialize quadrature points on the interval for all orders
Definition: quadraturerules.hh:469
~Jacobi1QuadratureRule1D()
Definition: quadraturerules.hh:281
int order(int m)
Definition: quadraturerules.hh:560
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:39
Jacobi-Gauss quadrature for alpha=1, beta=0.
Definition: quadraturerules.hh:271
~GaussQuadratureRule1D()
Definition: quadraturerules.hh:219
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
Definition: quadraturerules.hh:80
ct weight_
Definition: quadraturerules.hh:69
Definition: quadraturerules.hh:702
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:50
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:29
Definition: quadraturerules.hh:646
Gauss quadrature rule in 1D.
Definition: quadraturerules.hh:211
Definition: quadraturerules.hh:82
Definition: quadraturerules.hh:462