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>
46 enum { dimension = dim };
52 typedef Dune::FieldVector<ct,dim>
Vector;
80 namespace QuadratureType {
95 template<
typename ct,
int dim>
120 virtual int order ()
const {
return delivered_order; }
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);
230 #pragma GCC diagnostic push 231 #pragma GCC diagnostic ignored "-Wdeprecated-declarations" 236 return instance()._rule(gt,p,qt);
239 #pragma GCC diagnostic pop 243 #include "quadraturerules/pointquadrature.hh" 248 template<typename ct, bool fundamental = std::numeric_limits<ct>::is_specialized>
250 template<
typename ct>
252 static void init(
int p,
253 std::vector< FieldVector<ct, 1> > & _points,
254 std::vector< ct > & _weight,
255 int & delivered_order);
257 template<
typename ct>
259 static void init(
int p,
260 std::vector< FieldVector<ct, 1> > & _points,
261 std::vector< ct > & _weight,
262 int & delivered_order);
266 template<
typename ct>
273 enum { highest_order=61 };
282 std::vector< FieldVector<ct, dim> > _points;
283 std::vector< ct > _weight;
286 (p, _points, _weight, this->delivered_order);
288 assert(_points.size() == _weight.size());
289 for (
size_t i = 0; i < _points.size(); i++)
299 #define DUNE_INCLUDING_IMPLEMENTATION 300 #include "quadraturerules/gauss_imp.hh" 305 template<
typename ct,
306 bool fundamental = std::numeric_limits<ct>::is_specialized>
308 template<
typename ct>
310 static void init(
int p,
311 std::vector< FieldVector<ct, 1> > & _points,
312 std::vector< ct > & _weight,
313 int & delivered_order);
315 template<
typename ct>
317 static void init(
int p,
318 std::vector< FieldVector<ct, 1> > & _points,
319 std::vector< ct > & _weight,
320 int & delivered_order);
326 template<
typename ct>
335 enum { highest_order=61 };
344 std::vector< FieldVector<ct, dim> > _points;
345 std::vector< ct > _weight;
350 (p, _points, _weight, deliveredOrder_);
351 this->delivered_order = deliveredOrder_;
352 assert(_points.size() == _weight.size());
353 for (
size_t i = 0; i < _points.size(); i++)
365 #define DUNE_INCLUDING_IMPLEMENTATION 366 #include "quadraturerules/jacobi_1_0_imp.hh" 371 template<
typename ct,
372 bool fundamental = std::numeric_limits<ct>::is_specialized>
374 template<
typename ct>
376 static void init(
int p,
377 std::vector< FieldVector<ct, 1> > & _points,
378 std::vector< ct > & _weight,
379 int & delivered_order);
381 template<
typename ct>
383 static void init(
int p,
384 std::vector< FieldVector<ct, 1> > & _points,
385 std::vector< ct > & _weight,
386 int & delivered_order);
392 template<
typename ct>
401 enum { highest_order=61 };
410 std::vector< FieldVector<ct, dim> > _points;
411 std::vector< ct > _weight;
416 (p, _points, _weight, deliveredOrder_);
418 this->delivered_order = deliveredOrder_;
419 assert(_points.size() == _weight.size());
420 for (
size_t i = 0; i < _points.size(); i++)
432 #define DUNE_INCLUDING_IMPLEMENTATION 433 #include "quadraturerules/jacobi_2_0_imp.hh" 438 template<
typename ct,
439 bool fundamental = std::numeric_limits<ct>::is_specialized>
441 template<
typename ct>
443 static void init(
int p,
444 std::vector< FieldVector<ct, 1> > & _points,
445 std::vector< ct > & _weight,
446 int & delivered_order);
448 template<
typename ct>
450 static void init(
int p,
451 std::vector< FieldVector<ct, 1> > & _points,
452 std::vector< ct > & _weight,
453 int & delivered_order);
459 template<
typename ct>
468 enum { highest_order=31 };
477 std::vector< FieldVector<ct, dim> > _points;
478 std::vector< ct > _weight;
483 (p, _points, _weight, deliveredOrder_);
485 this->delivered_order = deliveredOrder_;
486 assert(_points.size() == _weight.size());
487 for (
size_t i = 0; i < _points.size(); i++)
499 #define DUNE_INCLUDING_IMPLEMENTATION 500 #include "quadraturerules/gausslobatto_imp.hh" 502 #include "quadraturerules/tensorproductquadrature.hh" 504 #include "quadraturerules/simplexquadrature.hh" 522 enum { highest_order=2 };
556 W[m][0] = 0.16666666666666666 / 2.0;
557 W[m][1] = 0.16666666666666666 / 2.0;
558 W[m][2] = 0.16666666666666666 / 2.0;
559 W[m][3] = 0.16666666666666666 / 2.0;
560 W[m][4] = 0.16666666666666666 / 2.0;
561 W[m][5] = 0.16666666666666666 / 2.0;
568 G[m][0][0] =0.66666666666666666 ;
569 G[m][0][1] =0.16666666666666666 ;
570 G[m][0][2] =0.211324865405187 ;
572 G[m][1][0] = 0.16666666666666666;
573 G[m][1][1] =0.66666666666666666 ;
574 G[m][1][2] = 0.211324865405187;
576 G[m][2][0] = 0.16666666666666666;
577 G[m][2][1] = 0.16666666666666666;
578 G[m][2][2] = 0.211324865405187;
580 G[m][3][0] = 0.66666666666666666;
581 G[m][3][1] = 0.16666666666666666;
582 G[m][3][2] = 0.788675134594813;
584 G[m][4][0] = 0.16666666666666666;
585 G[m][4][1] = 0.66666666666666666;
586 G[m][4][2] = 0.788675134594813;
588 G[m][5][0] = 0.16666666666666666;
589 G[m][5][1] = 0.16666666666666666;
590 G[m][5][2] = 0.788675134594813;
592 W[m][0] = 0.16666666666666666 / 2.0;
593 W[m][1] = 0.16666666666666666 / 2.0;
594 W[m][2] = 0.16666666666666666 / 2.0;
595 W[m][3] = 0.16666666666666666 / 2.0;
596 W[m][4] = 0.16666666666666666 / 2.0;
597 W[m][5] = 0.16666666666666666 / 2.0;
604 FieldVector<double, 3>
point(
int m,
int i)
622 FieldVector<double, 3> G[MAXP+1][MAXP];
624 double W[MAXP+1][MAXP];
648 template<
typename ct,
int dim>
654 template<
typename ct>
663 enum { highest_order = 2 };
672 for(
int i=0; i<m; ++i)
674 FieldVector<ct,3> local;
675 for (
int k=0; k<d; k++)
691 template<
typename ctype,
int dim>
697 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
701 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
705 template<
typename ctype>
714 return std::numeric_limits<int>::max();
716 DUNE_THROW(Exception,
"Unknown GeometryType");
722 return PointQuadratureRule<ctype>();
724 DUNE_THROW(Exception,
"Unknown GeometryType");
728 template<
typename ctype>
747 DUNE_THROW(Exception,
"Unknown QuadratureType");
750 DUNE_THROW(Exception,
"Unknown GeometryType");
766 DUNE_THROW(Exception,
"Unknown QuadratureType");
769 DUNE_THROW(Exception,
"Unknown GeometryType");
773 template<
typename ctype>
781 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
784 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
791 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
793 return SimplexQuadratureRule<ctype,dim>(p);
795 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
799 template<
typename ctype>
807 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
810 (order,
unsigned(SimplexQuadratureRule<ctype,dim>::highest_order));
820 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
822 return SimplexQuadratureRule<ctype,dim>(p);
826 && p <= PrismQuadratureRule<ctype,dim>::highest_order)
830 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
836 #endif // DUNE_GEOMETRY_QUADRATURERULES_HH int delivered_order
Definition: quadraturerules.hh:132
~GaussLobattoQuadratureRule1D()
Definition: quadraturerules.hh:470
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:577
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:633
constexpr bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:522
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid.
Definition: quadraturerules.hh:108
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:572
Definition: quadraturerules.hh:85
~PrismQuadratureRule()
Definition: quadraturerules.hh:665
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition: quadraturerules.hh:35
Quadrature rules for prisms.
Definition: quadraturerules.hh:649
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition: quadraturerules.hh:52
const ct & weight() const
return weight associated with integration point i
Definition: quadraturerules.hh:67
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
Definition: quadraturerules.hh:87
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:763
Definition: quadraturerules.hh:249
Definition: quadraturerules.hh:84
Definition: quadraturerules.hh:800
ct Field
Number type used for coordinates and quadrature weights.
Definition: quadraturerules.hh:49
Enum
Definition: quadraturerules.hh:81
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
constexpr bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:547
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition: quadraturerules.hh:128
const Vector & position() const
return local coordinates of integration point i
Definition: quadraturerules.hh:61
virtual int order() const
return order
Definition: quadraturerules.hh:120
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition: type.hh:282
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:557
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:96
PrismQuadraturePoints()
initialize quadrature points on the interval for all orders
Definition: quadraturerules.hh:525
GeometryType geometry_type
Definition: quadraturerules.hh:131
Definition: quadraturerules.hh:88
Definition: quadraturerules.hh:706
~GaussQuadratureRule1D()
Definition: quadraturerules.hh:275
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:276
virtual ~QuadratureRule()
Definition: quadraturerules.hh:124
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:642
QuadratureRule()
Default constructor.
Definition: quadraturerules.hh:104
virtual GeometryType type() const
return type of element
Definition: quadraturerules.hh:123
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:460
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:143
static constexpr std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition: typeindex.hh:56
Definition: quadraturerules.hh:774
Single evaluation point in a quadrature rule.
Definition: quadraturerules.hh:43
Quadrature rules for prisms.
Definition: quadraturerules.hh:655
A unique label for each type of element that can occur in a grid.
Jacobi-Gauss quadrature for alpha=1, beta=0.
Definition: quadraturerules.hh:327
Definition: quadraturerules.hh:307
~Jacobi1QuadratureRule1D()
Definition: quadraturerules.hh:337
Helper classes to provide indices for geometrytypes for use in a vector.
static PrismQuadraturePoints< 3 > prqp
Definition: quadraturerules.hh:634
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType...
Definition: quadraturerules.hh:137
Definition: quadraturerules.hh:729
Definition: quadraturerules.hh:440
FieldVector< ct, dim > local
Definition: quadraturerules.hh:73
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:517
static constexpr std::size_t index(const GeometryType >)
Compute the index for the given geometry type within its dimension.
Definition: typeindex.hh:68
Gauss quadrature rule in 1D.
Definition: quadraturerules.hh:267
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
Definition: quadraturerules.hh:518
Definition: affinegeometry.hh:18
Jacobi-Gauss quadrature for alpha=2, beta=0.
Definition: quadraturerules.hh:393
ct weight_
Definition: quadraturerules.hh:74
Definition: quadraturerules.hh:373
FieldVector< double, 3 > point(int m, int i)
Definition: quadraturerules.hh:604
double weight(int m, int i)
Definition: quadraturerules.hh:610
constexpr GeometryType line
GeometryType representing a line.
Definition: type.hh:733
Singleton holding the Prism Quadrature points.
Definition: quadraturerules.hh:641
int order(int m)
Definition: quadraturerules.hh:616
Definition: quadraturerules.hh:82
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition: quadraturerules.hh:55
Definition: quadraturerules.hh:514
~Jacobi2QuadratureRule1D()
Definition: quadraturerules.hh:403