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>
20 #include <dune/geometry/quadraturerules/nocopyvector.hh>
42 template<
typename ct,
int dim>
52 typedef Dune::FieldVector<ct,dim>
Vector;
80 namespace QuadratureType {
95 template<
typename ct,
int dim>
128 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
142 template<
typename ctype,
int dim>
155 typedef NoCopyVector<std::pair<std::once_flag, QuadratureRule> >
156 QuadratureOrderVector;
159 static void initQuadratureOrderVector(QuadratureOrderVector *qov,
170 typedef NoCopyVector<std::pair<std::once_flag, QuadratureOrderVector> >
174 static void initGeometryTypeVector(GeometryTypeVector *gtv)
182 assert(t.
dim()==dim);
184 DUNE_ASSERT_CALL_ONCE();
186 static NoCopyVector<std::pair<
191 auto & quadratureTypeLevel = quadratureCache[qt];
192 std::call_once(quadratureTypeLevel.first, initGeometryTypeVector,
193 &quadratureTypeLevel.second);
195 auto & geometryTypeLevel =
197 std::call_once(geometryTypeLevel.first, initQuadratureOrderVector,
198 &geometryTypeLevel.second, qt, t);
201 auto & quadratureOrderLevel = geometryTypeLevel.second[dim == 0 ? 0 : p];
202 std::call_once(quadratureOrderLevel.first, initQuadratureRule,
203 &quadratureOrderLevel.second, qt, t, p);
205 return quadratureOrderLevel.second;
227 return instance()._rule(t,p,qt);
233 return instance()._rule(gt,p,qt);
239 #include "quadraturerules/pointquadrature.hh"
244 template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
246 template<
typename ct>
248 static void init(
int p,
249 std::vector< FieldVector<ct, 1> > & _points,
250 std::vector< ct > & _weight,
251 int & delivered_order);
253 template<
typename ct>
255 static void init(
int p,
256 std::vector< FieldVector<ct, 1> > & _points,
257 std::vector< ct > & _weight,
258 int & delivered_order);
262 template<
typename ct>
278 std::vector< FieldVector<ct, dim> > _points;
279 std::vector< ct > _weight;
284 assert(_points.size() == _weight.size());
285 for (
size_t i = 0; i < _points.size(); i++)
290 extern template GaussQuadratureRule1D<float>::GaussQuadratureRule1D(
int);
291 extern template GaussQuadratureRule1D<double>::GaussQuadratureRule1D(
int);
295 #define DUNE_INCLUDING_IMPLEMENTATION
296 #include "quadraturerules/gauss_imp.hh"
301 template<
typename ct,
302 bool fundamental = std::numeric_limits<ct>::is_specialized>
304 template<
typename ct>
306 static void init(
int p,
307 std::vector< FieldVector<ct, 1> > & _points,
308 std::vector< ct > & _weight,
309 int & delivered_order);
311 template<
typename ct>
313 static void init(
int p,
314 std::vector< FieldVector<ct, 1> > & _points,
315 std::vector< ct > & _weight,
316 int & delivered_order);
322 template<
typename ct>
331 enum { highest_order=61 };
340 std::vector< FieldVector<ct, dim> > _points;
341 std::vector< ct > _weight;
346 (p, _points, _weight, delivered_order);
347 this->delivered_order = delivered_order;
348 assert(_points.size() == _weight.size());
349 for (
size_t i = 0; i < _points.size(); i++)
355 extern template Jacobi1QuadratureRule1D<float>::Jacobi1QuadratureRule1D(
int);
356 extern template Jacobi1QuadratureRule1D<double>::Jacobi1QuadratureRule1D(
int);
361 #define DUNE_INCLUDING_IMPLEMENTATION
362 #include "quadraturerules/jacobi_1_0_imp.hh"
367 template<
typename ct,
368 bool fundamental = std::numeric_limits<ct>::is_specialized>
370 template<
typename ct>
372 static void init(
int p,
373 std::vector< FieldVector<ct, 1> > & _points,
374 std::vector< ct > & _weight,
375 int & delivered_order);
377 template<
typename ct>
379 static void init(
int p,
380 std::vector< FieldVector<ct, 1> > & _points,
381 std::vector< ct > & _weight,
382 int & delivered_order);
388 template<
typename ct>
397 enum { highest_order=61 };
406 std::vector< FieldVector<ct, dim> > _points;
407 std::vector< ct > _weight;
412 (p, _points, _weight, delivered_order);
414 this->delivered_order = delivered_order;
415 assert(_points.size() == _weight.size());
416 for (
size_t i = 0; i < _points.size(); i++)
422 extern template Jacobi2QuadratureRule1D<float>::Jacobi2QuadratureRule1D(
int);
423 extern template Jacobi2QuadratureRule1D<double>::Jacobi2QuadratureRule1D(
int);
428 #define DUNE_INCLUDING_IMPLEMENTATION
429 #include "quadraturerules/jacobi_2_0_imp.hh"
434 template<
typename ct,
435 bool fundamental = std::numeric_limits<ct>::is_specialized>
437 template<
typename ct>
439 static void init(
int p,
440 std::vector< FieldVector<ct, 1> > & _points,
441 std::vector< ct > & _weight,
442 int & delivered_order);
444 template<
typename ct>
446 static void init(
int p,
447 std::vector< FieldVector<ct, 1> > & _points,
448 std::vector< ct > & _weight,
449 int & delivered_order);
455 template<
typename ct>
464 enum { highest_order=31 };
473 std::vector< FieldVector<ct, dim> > _points;
474 std::vector< ct > _weight;
479 (p, _points, _weight, delivered_order);
481 this->delivered_order = delivered_order;
482 assert(_points.size() == _weight.size());
483 for (
size_t i = 0; i < _points.size(); i++)
489 extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(
int);
490 extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(
int);
495 #define DUNE_INCLUDING_IMPLEMENTATION
496 #include "quadraturerules/gausslobatto_imp.hh"
498 #include "quadraturerules/tensorproductquadrature.hh"
500 #include "quadraturerules/simplexquadrature.hh"
518 enum { highest_order=2 };
552 W[m][0] = 0.16666666666666666 / 2.0;
553 W[m][1] = 0.16666666666666666 / 2.0;
554 W[m][2] = 0.16666666666666666 / 2.0;
555 W[m][3] = 0.16666666666666666 / 2.0;
556 W[m][4] = 0.16666666666666666 / 2.0;
557 W[m][5] = 0.16666666666666666 / 2.0;
564 G[m][0][0] =0.66666666666666666 ;
565 G[m][0][1] =0.16666666666666666 ;
566 G[m][0][2] =0.211324865405187 ;
568 G[m][1][0] = 0.16666666666666666;
569 G[m][1][1] =0.66666666666666666 ;
570 G[m][1][2] = 0.211324865405187;
572 G[m][2][0] = 0.16666666666666666;
573 G[m][2][1] = 0.16666666666666666;
574 G[m][2][2] = 0.211324865405187;
576 G[m][3][0] = 0.66666666666666666;
577 G[m][3][1] = 0.16666666666666666;
578 G[m][3][2] = 0.788675134594813;
580 G[m][4][0] = 0.16666666666666666;
581 G[m][4][1] = 0.66666666666666666;
582 G[m][4][2] = 0.788675134594813;
584 G[m][5][0] = 0.16666666666666666;
585 G[m][5][1] = 0.16666666666666666;
586 G[m][5][2] = 0.788675134594813;
588 W[m][0] = 0.16666666666666666 / 2.0;
589 W[m][1] = 0.16666666666666666 / 2.0;
590 W[m][2] = 0.16666666666666666 / 2.0;
591 W[m][3] = 0.16666666666666666 / 2.0;
592 W[m][4] = 0.16666666666666666 / 2.0;
593 W[m][5] = 0.16666666666666666 / 2.0;
600 FieldVector<double, 3>
point(
int m,
int i)
618 FieldVector<double, 3> G[MAXP+1][MAXP];
620 double W[MAXP+1][MAXP];
644 template<
typename ct,
int dim>
650 template<
typename ct>
659 enum { highest_order = 2 };
668 for(
int i=0; i<m; ++i)
670 FieldVector<ct,3> local;
671 for (
int k=0; k<d; k++)
687 template<
typename ctype,
int dim>
688 class QuadratureRuleFactory {
693 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
697 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
701 template<
typename ctype>
710 return std::numeric_limits<int>::max();
712 DUNE_THROW(Exception,
"Unknown GeometryType");
718 return PointQuadratureRule<ctype>();
720 DUNE_THROW(Exception,
"Unknown GeometryType");
724 template<
typename ctype>
743 DUNE_THROW(Exception,
"Unknown QuadratureType");
746 DUNE_THROW(Exception,
"Unknown GeometryType");
762 DUNE_THROW(Exception,
"Unknown QuadratureType");
765 DUNE_THROW(Exception,
"Unknown GeometryType");
769 template<
typename ctype>
777 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
780 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
787 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
789 return SimplexQuadratureRule<ctype,dim>(p);
791 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
795 template<
typename ctype>
803 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
806 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
816 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
818 return SimplexQuadratureRule<ctype,dim>(p);
822 && p <= PrismQuadratureRule<ctype,dim>::highest_order)
824 return PrismQuadratureRule<ctype,dim>(p);
826 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
832 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH
Definition: quadraturerules.hh:88
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:225
int delivered_order
Definition: quadraturerules.hh:132
unsigned int id() const
Return the topology id the type.
Definition: type.hh:326
Definition: quadraturerules.hh:514
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:29
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:61
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:128
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:67
Definition: affinegeometry.hh:18
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:96
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:35
Definition: quadraturerules.hh:303
FieldVector< ct, dim > local
Definition: quadraturerules.hh:73
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition: quadraturerules.hh:218
bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:271
~GaussLobattoQuadratureRule1D()
Definition: quadraturerules.hh:466
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:55
Enum
Definition: quadraturerules.hh:81
Definition: quadraturerules.hh:87
Jacobi-Gauss quadrature for alpha=1, beta=0.
Definition: quadraturerules.hh:323
Quadrature rules for prisms.
Definition: quadraturerules.hh:645
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:630
Gauss quadrature rule in 1D.
Definition: quadraturerules.hh:263
bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:296
Definition: quadraturerules.hh:84
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid.
Definition: quadraturerules.hh:108
Definition: quadraturerules.hh:725
Definition: quadraturerules.hh:85
ct weight_
Definition: quadraturerules.hh:74
~Jacobi1QuadratureRule1D()
Definition: quadraturerules.hh:333
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:143
~GaussQuadratureRule1D()
Definition: quadraturerules.hh:271
virtual ~QuadratureRule()
Definition: quadraturerules.hh:124
Definition: quadraturerules.hh:114
Definition: quadraturerules.hh:46
Helper classes to provide indices for geometrytypes for use in a vector.
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:637
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:104
bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:266
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:389
PrismQuadraturePoints()
initialize quadrature points on the interval for all orders
Definition: quadraturerules.hh:521
int order(int m)
Definition: quadraturerules.hh:612
Definition: quadraturerules.hh:702
~Jacobi2QuadratureRule1D()
Definition: quadraturerules.hh:399
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:123
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition: quadraturerules.hh:111
ct CoordType
The type used for coordinates.
Definition: quadraturerules.hh:117
Definition: quadraturerules.hh:245
Definition: quadraturerules.hh:268
A unique label for each type of element that can occur in a grid.
static std::size_t index(const GeometryType >)
Compute the index for the given geometry type within its dimension.
Definition: typeindex.hh:70
Definition: quadraturerules.hh:796
Definition: quadraturerules.hh:82
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:638
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:230
GeometryType geometry_type
Definition: quadraturerules.hh:131
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:629
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
Definition: quadraturerules.hh:269
unsigned int dim() const
Return dimension of the type.
Definition: type.hh:321
Definition: quadraturerules.hh:510
Definition: quadraturerules.hh:369
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:456
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition: quadraturerules.hh:52
FieldVector< double, 3 > point(int m, int i)
Definition: quadraturerules.hh:600
double weight(int m, int i)
Definition: quadraturerules.hh:606
Quadrature rules for prisms.
Definition: quadraturerules.hh:651
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType...
Definition: quadraturerules.hh:137
virtual int order() const
return order
Definition: quadraturerules.hh:120
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:43
bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:306
~PrismQuadratureRule()
Definition: quadraturerules.hh:661
Definition: quadraturerules.hh:436
static DUNE_CONSTEXPR std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition: typeindex.hh:58
ct Field
Number type used for coordinates and quadrature weights.
Definition: quadraturerules.hh:49
Definition: quadraturerules.hh:770