6#ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
7#define DUNE_GEOMETRY_QUADRATURERULES_HH
16#include <dune/common/fvector.hh>
17#include <dune/common/exceptions.hh>
18#include <dune/common/stdstreams.hh>
19#include <dune/common/stdthread.hh>
20#include <dune/common/visibility.hh>
33 template<
typename ct,
int dim>
34 class QuadraturePoint;
39 template<
typename ct,
int dim>
40 struct tuple_size<
Dune::QuadraturePoint<ct,dim>> :
public std::integral_constant<std::size_t,2> {};
42 template<
typename ct,
int dim>
43 struct tuple_element<0,
Dune::QuadraturePoint<ct,dim>> {
using type = Dune::FieldVector<ct, dim>; };
45 template<
typename ct,
int dim>
46 struct tuple_element<1,
Dune::QuadraturePoint<ct,dim>> {
using type = ct; };
65 template<
typename ct,
int dim>
75 typedef Dune::FieldVector<ct,dim>
Vector;
111 template<std::
size_t index, std::enable_if_t<(index<=1),
int> = 0>
112 std::tuple_element_t<index, QuadraturePo
int<ct, dim>> get() const
114 if constexpr (index == 0) {
123 FieldVector<ct, dim> local;
130 namespace QuadratureType {
212 template<
typename ct,
int dim>
231 constexpr static int d = dim;
237 virtual int order ()
const {
return delivered_order; }
245 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
259 template<
typename ctype,
int dim>
266 using QuadratureOrderVector = std::vector<std::pair<std::once_flag, QuadratureRule> >;
269 using GeometryTypeVector = std::vector<std::pair<std::once_flag, QuadratureOrderVector> >;
272 using QuadratureCacheVector = std::vector<std::pair<std::once_flag, GeometryTypeVector> >;
277 assert(t.
dim()==dim);
279 DUNE_ASSERT_CALL_ONCE();
281 static QuadratureCacheVector quadratureCache(QuadratureType::size);
283 auto& [ onceFlagQuadratureType, geometryTypes ] = quadratureCache[qt];
285 std::call_once(onceFlagQuadratureType, [&types = geometryTypes]{
286 types = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
289 auto& [ onceFlagGeometryType, quadratureOrders ] = geometryTypes[LocalGeometryTypeIndex::index(t)];
291 std::call_once(onceFlagGeometryType, [&, &orders = quadratureOrders]{
294 orders = QuadratureOrderVector(numRules);
298 auto& [ onceFlagQuadratureOrder, quadratureRule ] = quadratureOrders[dim == 0 ? 0 : p];
300 std::call_once(onceFlagQuadratureOrder, [&, &rule = quadratureRule]{
304 return quadratureRule;
328 return instance()._rule(t,p,qt);
335 return instance()._rule(gt,p,qt);
341#define DUNE_INCLUDING_IMPLEMENTATION
344#include "quadraturerules/pointquadrature.hh"
346#include "quadraturerules/gausslobattoquadrature.hh"
347#include "quadraturerules/gaussquadrature.hh"
348#include "quadraturerules/gaussradauleftquadrature.hh"
349#include "quadraturerules/gaussradaurightquadrature.hh"
350#include "quadraturerules/jacobi1quadrature.hh"
351#include "quadraturerules/jacobi2quadrature.hh"
352#include "quadraturerules/jacobiNquadrature.hh"
354#include "quadraturerules/prismquadrature.hh"
356#include "quadraturerules/simplexquadrature.hh"
357#include "quadraturerules/tensorproductquadrature.hh"
359#undef DUNE_INCLUDING_IMPLEMENTATION
369 template<
typename ctype,
int dim>
375 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
379 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
383 template<
typename ctype>
386 constexpr static int dim = 0;
392 return std::numeric_limits<int>::max();
394 DUNE_THROW(Exception,
"Unknown GeometryType");
400 return PointQuadratureRule<ctype>();
402 DUNE_THROW(Exception,
"Unknown GeometryType");
406 template<
typename ctype>
409 constexpr static int dim = 1;
417 return GaussQuadratureRule1D<ctype>::highest_order;
419 return Jacobi1QuadratureRule1D<ctype>::highest_order;
421 return Jacobi2QuadratureRule1D<ctype>::highest_order;
423 return GaussLobattoQuadratureRule1D<ctype>::highest_order;
425 return JacobiNQuadratureRule1D<ctype>::maxOrder();
427 return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
429 return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
431 DUNE_THROW(Exception,
"Unknown QuadratureType");
434 DUNE_THROW(Exception,
"Unknown GeometryType");
442 return GaussQuadratureRule1D<ctype>(p);
444 return Jacobi1QuadratureRule1D<ctype>(p);
446 return Jacobi2QuadratureRule1D<ctype>(p);
448 return GaussLobattoQuadratureRule1D<ctype>(p);
450 return JacobiNQuadratureRule1D<ctype>(p);
452 return GaussRadauLeftQuadratureRule1D<ctype>(p);
454 return GaussRadauRightQuadratureRule1D<ctype>(p);
456 DUNE_THROW(Exception,
"Unknown QuadratureType");
459 DUNE_THROW(Exception,
"Unknown GeometryType");
463 template<
typename ctype>
466 constexpr static int dim = 2;
471 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
474 (order,
static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
481 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
483 return SimplexQuadratureRule<ctype,dim>(p);
485 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
489 template<
typename ctype>
492 constexpr static int dim = 3;
497 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
500 (order,
static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
503 (order,
static_cast<unsigned>(PrismQuadratureRule<ctype,dim>::highest_order));
511 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
513 return SimplexQuadratureRule<ctype,dim>(p);
517 && p <= PrismQuadratureRule<ctype,dim>::highest_order)
519 return PrismQuadratureRule<ctype,dim>(p);
521 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
525#ifndef DUNE_NO_EXTERN_QUADRATURERULES
526 extern template class GaussLobattoQuadratureRule<double, 1>;
527 extern template class GaussQuadratureRule<double, 1>;
528 extern template class GaussRadauLeftQuadratureRule<double, 1>;
529 extern template class GaussRadauRightQuadratureRule<double, 1>;
530 extern template class Jacobi1QuadratureRule<double, 1>;
531 extern template class Jacobi2QuadratureRule<double, 1>;
532 extern template class JacobiNQuadratureRule<double, 1>;
533 extern template class PrismQuadratureRule<double, 3>;
534 extern template class SimplexQuadratureRule<double, 2>;
535 extern template class SimplexQuadratureRule<double, 3>;
Helper classes to provide indices for geometrytypes for use in a vector.
A unique label for each type of element that can occur in a grid.
Definition affinegeometry.hh:21
Enum
Definition quadraturerules.hh:131
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:168
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:155
@ GaussRadauRight
Gauss-Radau rules including the right endpoint.
Definition quadraturerules.hh:193
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition quadraturerules.hh:148
@ GaussLobatto
Gauss-Lobatto rules.
Definition quadraturerules.hh:176
@ GaussRadauLeft
Gauss-Radau rules including the left endpoint.
Definition quadraturerules.hh:184
@ GaussLegendre
Gauss-Legendre rules (default)
Definition quadraturerules.hh:141
Single evaluation point in a quadrature rule.
Definition quadraturerules.hh:66
const Vector & position() const
return local coordinates of integration point i
Definition quadraturerules.hh:82
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition quadraturerules.hh:75
ct Field
Number type used for coordinates and quadrature weights.
Definition quadraturerules.hh:72
const ct & weight() const
return weight associated with integration point i
Definition quadraturerules.hh:88
ct weight_
Definition quadraturerules.hh:124
static constexpr int dimension
Dimension of the integration domain.
Definition quadraturerules.hh:69
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition quadraturerules.hh:78
FieldVector< ct, dim > local
Definition quadraturerules.hh:123
Dune::FieldVector< ct, dim > type
Definition quadraturerules.hh:43
ct type
Definition quadraturerules.hh:46
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition quadraturerules.hh:55
Abstract base class for quadrature rules.
Definition quadraturerules.hh:214
virtual ~QuadratureRule()
Definition quadraturerules.hh:241
virtual GeometryType type() const
return type of element
Definition quadraturerules.hh:240
int delivered_order
Definition quadraturerules.hh:249
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition quadraturerules.hh:228
GeometryType geometry_type
Definition quadraturerules.hh:248
ct CoordType
The type used for coordinates.
Definition quadraturerules.hh:234
QuadratureRule()
Default constructor.
Definition quadraturerules.hh:221
virtual int order() const
return order
Definition quadraturerules.hh:237
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid
Definition quadraturerules.hh:225
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition quadraturerules.hh:245
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition quadraturerules.hh:370
A container for all quadrature rules of dimension dim
Definition quadraturerules.hh:260
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition quadraturerules.hh:319
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:332
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:326
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114
constexpr bool isPrism() const
Return true if entity is a prism.
Definition type.hh:309
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition type.hh:279
constexpr unsigned int dim() const
Return dimension of the type.
Definition type.hh:360
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition type.hh:120
constexpr bool isLine() const
Return true if entity is a line segment.
Definition type.hh:284
constexpr unsigned int id() const
Return the topology id of the type.
Definition type.hh:365
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition type.hh:319