Loading...
Searching...
No Matches
quadraturerules.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
#define DUNE_GEOMETRY_QUADRATURERULES_HH
#include <algorithm>
#include <iostream>
#include <limits>
#include <mutex>
#include <utility>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/stdthread.hh>
#include <dune/common/visibility.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/typeindex.hh>
namespace Dune {
// forward declaration
template<typename ct, int dim>
class QuadraturePoint;
}
// class specialization of standard classes that allow to use structured bindings on QuadraturePoint
namespace std {
template<typename ct, int dim>
struct tuple_size<Dune::QuadraturePoint<ct,dim>> : public std::integral_constant<std::size_t,2> {};
template<typename ct, int dim>
struct tuple_element<0, Dune::QuadraturePoint<ct,dim>> { using type = Dune::FieldVector<ct, dim>; };
template<typename ct, int dim>
struct tuple_element<1, Dune::QuadraturePoint<ct,dim>> { using type = ct; };
}
namespace Dune {
class QuadratureOrderOutOfRange : public NotImplemented {};
template<typename ct, int dim>
class QuadraturePoint {
public:
constexpr static int dimension = dim;
typedef ct Field;
typedef Dune::FieldVector<ct,dim> Vector;
QuadraturePoint (const Vector& x, ct w) : local(x), weight_(w)
{}
const Vector& position () const
{
return local;
}
const ct &weight () const
{
return weight_;
}
template<std::size_t index, std::enable_if_t<(index<=1), int> = 0>
std::tuple_element_t<index, QuadraturePoint<ct, dim>> get() const
{
if constexpr (index == 0) {
return local;
}
else {
return weight_;
}
}
protected:
FieldVector<ct, dim> local;
ct weight_;
};
namespace QuadratureType {
enum Enum {
GaussLegendre = 0,
GaussJacobi_1_0 = 1,
GaussJacobi_2_0 = 2,
GaussJacobi_n_0 = 3,
GaussLobatto = 4,
GaussRadauLeft = 5,
GaussRadauRight = 6,
size
};
}
template<typename ct, int dim>
class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
{
public:
QuadratureRule() : delivered_order(-1) {}
protected:
QuadratureRule(GeometryType t) : geometry_type(t), delivered_order(-1) {}
QuadratureRule(GeometryType t, int order) : geometry_type(t), delivered_order(order) {}
public:
constexpr static int d = dim;
typedef ct CoordType;
virtual int order () const { return delivered_order; }
virtual GeometryType type () const { return geometry_type; }
virtual ~QuadratureRule(){}
typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;
protected:
GeometryType geometry_type;
int delivered_order;
};
// Forward declaration of the factory class,
// needed internally by the QuadratureRules container class.
template<typename ctype, int dim> class QuadratureRuleFactory;
template<typename ctype, int dim>
class QuadratureRules {
using QuadratureRule = Dune::QuadratureRule<ctype, dim>;
// indexed by quadrature order
using QuadratureOrderVector = std::vector<std::pair<std::once_flag, QuadratureRule> >;
// indexed by geometry type
using GeometryTypeVector = std::vector<std::pair<std::once_flag, QuadratureOrderVector> >;
// indexed by quadrature type enum
using QuadratureCacheVector = std::vector<std::pair<std::once_flag, GeometryTypeVector> >;
DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre) const
{
assert(t.dim()==dim);
DUNE_ASSERT_CALL_ONCE();
static QuadratureCacheVector quadratureCache(QuadratureType::size);
auto& [ onceFlagQuadratureType, geometryTypes ] = quadratureCache[qt];
// initialize geometry types for this quadrature type once
std::call_once(onceFlagQuadratureType, [&types = geometryTypes]{
types = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
});
auto& [ onceFlagGeometryType, quadratureOrders ] = geometryTypes[LocalGeometryTypeIndex::index(t)];
// initialize quadrature orders for this geometry type and quadrature type once
std::call_once(onceFlagGeometryType, [&, &orders = quadratureOrders]{
// we only need one quadrature rule for points, not maxint
const auto numRules = dim == 0 ? 1 : QuadratureRuleFactory<ctype,dim>::maxOrder(t, qt)+1;
orders = QuadratureOrderVector(numRules);
});
// we only have one quadrature rule for points
auto& [ onceFlagQuadratureOrder, quadratureRule ] = quadratureOrders[dim == 0 ? 0 : p];
// initialize quadrature rule once
std::call_once(onceFlagQuadratureOrder, [&, &rule = quadratureRule]{
rule = QuadratureRuleFactory<ctype,dim>::rule(t, p, qt);
});
return quadratureRule;
}
DUNE_EXPORT static QuadratureRules& instance()
{
static QuadratureRules instance;
return instance;
}
QuadratureRules () = default;
public:
static unsigned
maxOrder(const GeometryType& t,
QuadratureType::Enum qt=QuadratureType::GaussLegendre)
{
return QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt);
}
static const QuadratureRule& rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
{
return instance()._rule(t,p,qt);
}
static const QuadratureRule& rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
{
GeometryType gt(t,dim);
return instance()._rule(gt,p,qt);
}
};
} // end namespace Dune
#define DUNE_INCLUDING_IMPLEMENTATION
// 0d rules
#include "quadraturerules/pointquadrature.hh"
// 1d rules
#include "quadraturerules/gausslobattoquadrature.hh"
#include "quadraturerules/gaussquadrature.hh"
#include "quadraturerules/gaussradauleftquadrature.hh"
#include "quadraturerules/gaussradaurightquadrature.hh"
#include "quadraturerules/jacobi1quadrature.hh"
#include "quadraturerules/jacobi2quadrature.hh"
#include "quadraturerules/jacobiNquadrature.hh"
// 3d rules
#include "quadraturerules/prismquadrature.hh"
// general rules
#include "quadraturerules/simplexquadrature.hh"
#include "quadraturerules/tensorproductquadrature.hh"
#undef DUNE_INCLUDING_IMPLEMENTATION
namespace Dune {
template<typename ctype, int dim>
class QuadratureRuleFactory {
private:
friend class QuadratureRules<ctype, dim>;
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
{
return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
}
static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
{
return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
}
};
template<typename ctype>
class QuadratureRuleFactory<ctype, 0> {
private:
constexpr static int dim = 0;
friend class QuadratureRules<ctype, dim>;
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum)
{
if (t.isVertex())
{
return std::numeric_limits<int>::max();
}
DUNE_THROW(Exception, "Unknown GeometryType");
}
static QuadratureRule<ctype, dim> rule(const GeometryType& t, int , QuadratureType::Enum)
{
if (t.isVertex())
{
return PointQuadratureRule<ctype>();
}
DUNE_THROW(Exception, "Unknown GeometryType");
}
};
template<typename ctype>
class QuadratureRuleFactory<ctype, 1> {
private:
constexpr static int dim = 1;
friend class QuadratureRules<ctype, dim>;
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
{
if (t.isLine())
{
switch (qt) {
case QuadratureType::GaussLegendre :
return GaussQuadratureRule1D<ctype>::highest_order;
case QuadratureType::GaussJacobi_1_0 :
return Jacobi1QuadratureRule1D<ctype>::highest_order;
case QuadratureType::GaussJacobi_2_0 :
return Jacobi2QuadratureRule1D<ctype>::highest_order;
case QuadratureType::GaussLobatto :
return GaussLobattoQuadratureRule1D<ctype>::highest_order;
case QuadratureType::GaussJacobi_n_0 :
return JacobiNQuadratureRule1D<ctype>::maxOrder();
case QuadratureType::GaussRadauLeft :
return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
case QuadratureType::GaussRadauRight :
return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
default :
DUNE_THROW(Exception, "Unknown QuadratureType");
}
}
DUNE_THROW(Exception, "Unknown GeometryType");
}
static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
{
if (t.isLine())
{
switch (qt) {
case QuadratureType::GaussLegendre :
return GaussQuadratureRule1D<ctype>(p);
case QuadratureType::GaussJacobi_1_0 :
return Jacobi1QuadratureRule1D<ctype>(p);
case QuadratureType::GaussJacobi_2_0 :
return Jacobi2QuadratureRule1D<ctype>(p);
case QuadratureType::GaussLobatto :
return GaussLobattoQuadratureRule1D<ctype>(p);
case QuadratureType::GaussJacobi_n_0 :
return JacobiNQuadratureRule1D<ctype>(p);
case QuadratureType::GaussRadauLeft :
return GaussRadauLeftQuadratureRule1D<ctype>(p);
case QuadratureType::GaussRadauRight :
return GaussRadauRightQuadratureRule1D<ctype>(p);
default :
DUNE_THROW(Exception, "Unknown QuadratureType");
}
}
DUNE_THROW(Exception, "Unknown GeometryType");
}
};
template<typename ctype>
class QuadratureRuleFactory<ctype, 2> {
private:
constexpr static int dim = 2;
friend class QuadratureRules<ctype, dim>;
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
{
unsigned order =
TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
if (t.isSimplex())
order = std::max
(order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
return order;
}
static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
{
if (t.isSimplex()
&& ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
&& p <= SimplexQuadratureRule<ctype,dim>::highest_order)
{
return SimplexQuadratureRule<ctype,dim>(p);
}
return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
}
};
template<typename ctype>
class QuadratureRuleFactory<ctype, 3> {
private:
constexpr static int dim = 3;
friend class QuadratureRules<ctype, dim>;
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
{
unsigned order =
TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
if (t.isSimplex())
order = std::max
(order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
if (t.isPrism())
order = std::max
(order, static_cast<unsigned>(PrismQuadratureRule<ctype,dim>::highest_order));
return order;
}
static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
{
if (t.isSimplex()
&& ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
&& p <= SimplexQuadratureRule<ctype,dim>::highest_order)
{
return SimplexQuadratureRule<ctype,dim>(p);
}
if (t.isPrism()
&& qt == QuadratureType::GaussLegendre
&& p <= PrismQuadratureRule<ctype,dim>::highest_order)
{
return PrismQuadratureRule<ctype,dim>(p);
}
return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
}
};
#ifndef DUNE_NO_EXTERN_QUADRATURERULES
extern template class GaussLobattoQuadratureRule<double, 1>;
extern template class GaussQuadratureRule<double, 1>;
extern template class GaussRadauLeftQuadratureRule<double, 1>;
extern template class GaussRadauRightQuadratureRule<double, 1>;
extern template class Jacobi1QuadratureRule<double, 1>;
extern template class Jacobi2QuadratureRule<double, 1>;
extern template class JacobiNQuadratureRule<double, 1>;
extern template class PrismQuadratureRule<double, 3>;
extern template class SimplexQuadratureRule<double, 2>;
extern template class SimplexQuadratureRule<double, 3>;
#endif // !DUNE_NO_EXTERN_QUADRATURERULES
} // end namespace
#endif // DUNE_GEOMETRY_QUADRATURERULES_HH
typeindex.hh
Helper classes to provide indices for geometrytypes for use in a vector.
type.hh
A unique label for each type of element that can occur in a grid.
std
STL namespace.
Dune
Definition affinegeometry.hh:22
Dune::QuadratureType::Enum
Enum
Definition quadraturerules.hh:131
Dune::QuadratureType::GaussJacobi_n_0
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:168
Dune::QuadratureType::GaussJacobi_2_0
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:155
Dune::QuadratureType::GaussRadauRight
@ GaussRadauRight
Gauss-Radau rules including the right endpoint.
Definition quadraturerules.hh:193
Dune::QuadratureType::GaussJacobi_1_0
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition quadraturerules.hh:148
Dune::QuadratureType::size
@ size
Definition quadraturerules.hh:194
Dune::QuadratureType::GaussLobatto
@ GaussLobatto
Gauss-Lobatto rules.
Definition quadraturerules.hh:176
Dune::QuadratureType::GaussRadauLeft
@ GaussRadauLeft
Gauss-Radau rules including the left endpoint.
Definition quadraturerules.hh:184
Dune::QuadratureType::GaussLegendre
@ GaussLegendre
Gauss-Legendre rules (default).
Definition quadraturerules.hh:141
Dune::QuadraturePoint
Single evaluation point in a quadrature rule.
Definition quadraturerules.hh:66
Dune::QuadraturePoint::position
const Vector & position() const
return local coordinates of integration point i
Definition quadraturerules.hh:82
Dune::QuadraturePoint::Vector
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition quadraturerules.hh:75
Dune::QuadraturePoint::Field
ct Field
Number type used for coordinates and quadrature weights.
Definition quadraturerules.hh:72
Dune::QuadraturePoint::weight
const ct & weight() const
return weight associated with integration point i
Definition quadraturerules.hh:88
Dune::QuadraturePoint::weight_
ct weight_
Definition quadraturerules.hh:124
Dune::QuadraturePoint::dimension
static constexpr int dimension
Dimension of the integration domain.
Definition quadraturerules.hh:69
Dune::QuadraturePoint::QuadraturePoint
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition quadraturerules.hh:78
Dune::QuadraturePoint::local
FieldVector< ct, dim > local
Definition quadraturerules.hh:123
std::tuple_element< 0, Dune::QuadraturePoint< ct, dim > >::type
Dune::FieldVector< ct, dim > type
Definition quadraturerules.hh:43
std::tuple_element< 1, Dune::QuadraturePoint< ct, dim > >::type
ct type
Definition quadraturerules.hh:46
Dune::QuadratureOrderOutOfRange
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition quadraturerules.hh:55
Dune::QuadratureRule
Abstract base class for quadrature rules.
Definition quadraturerules.hh:214
Dune::QuadratureRule::~QuadratureRule
virtual ~QuadratureRule()
Definition quadraturerules.hh:241
Dune::QuadratureRule< ctype, dim >::d
static constexpr int d
Definition quadraturerules.hh:231
Dune::QuadratureRule::type
virtual GeometryType type() const
return type of element
Definition quadraturerules.hh:240
Dune::QuadratureRule< ctype, dim >::delivered_order
int delivered_order
Definition quadraturerules.hh:249
Dune::QuadratureRule::QuadratureRule
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition quadraturerules.hh:228
Dune::QuadratureRule< ctype, dim >::geometry_type
GeometryType geometry_type
Definition quadraturerules.hh:248
Dune::QuadratureRule::CoordType
ct CoordType
The type used for coordinates.
Definition quadraturerules.hh:234
Dune::QuadratureRule::QuadratureRule
QuadratureRule()
Default constructor.
Definition quadraturerules.hh:221
Dune::QuadratureRule< ctype, dim >::order
virtual int order() const
Definition quadraturerules.hh:237
Dune::QuadratureRule::QuadratureRule
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid.
Definition quadraturerules.hh:225
Dune::QuadratureRule::iterator
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition quadraturerules.hh:245
Dune::QuadratureRuleFactory
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition quadraturerules.hh:370
Dune::QuadratureRules
A container for all quadrature rules of dimension dim.
Definition quadraturerules.hh:260
Dune::QuadratureRules::maxOrder
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
Dune::QuadratureRules::rule
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
Dune::QuadratureRules::rule
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
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114
Dune::GeometryType::isPrism
constexpr bool isPrism() const
Return true if entity is a prism.
Definition type.hh:309
Dune::GeometryType::isVertex
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition type.hh:279
Dune::GeometryType::dim
constexpr unsigned int dim() const
Return dimension of the type.
Definition type.hh:360
Dune::GeometryType::BasicType
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition type.hh:120
Dune::GeometryType::isLine
constexpr bool isLine() const
Return true if entity is a line segment.
Definition type.hh:284
Dune::GeometryType::id
constexpr unsigned int id() const
Return the topology id of the type.
Definition type.hh:365
Dune::GeometryType::isSimplex
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition type.hh:319
Dune::LocalGeometryTypeIndex::size
static constexpr std::size_t size(std::size_t dim)
Compute total number of geometry types for the given dimension.
Definition typeindex.hh:61
Dune::LocalGeometryTypeIndex::index
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type within its dimension.
Definition typeindex.hh:73