Commit 0b7983ce authored by Simon Praetorius's avatar Simon Praetorius
Browse files

Use a container with flexible key type as cache for quadrature rules

parent 5521af30
Pipeline #37484 passed with stage
in 8 minutes and 16 seconds
......@@ -7,7 +7,7 @@
#include <algorithm>
#include <iostream>
#include <limits>
#include <mutex>
#include <map>
#include <utility>
#include <vector>
......@@ -195,98 +195,44 @@ namespace Dune {
\ingroup Quadrature
template<typename ctype, int dim>
class QuadratureRules {
/** \brief Internal short-hand notation for the type of quadrature rules this container contains */
typedef Dune::QuadratureRule<ctype, dim> QuadratureRule;
//! \brief a quadrature rule (for each quadrature order, geometry type,
//! and quadrature type)
static void initQuadratureRule(QuadratureRule *qr, QuadratureType::Enum qt,
const GeometryType &t, int p)
*qr = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
typedef std::vector<std::pair<std::once_flag, QuadratureRule> >
QuadratureOrderVector; // indexed by quadrature order
//! \brief initialize the vector indexed by the quadrature order (for each
//! geometry type and quadrature type)
static void initQuadratureOrderVector(QuadratureOrderVector *qov,
QuadratureType::Enum qt,
const GeometryType &t)
if(dim == 0)
// we only need one quadrature rule for points, not maxint
*qov = QuadratureOrderVector(1);
*qov = QuadratureOrderVector(QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt)+1);
typedef std::vector<std::pair<std::once_flag, QuadratureOrderVector> >
GeometryTypeVector; // indexed by geometry type
//! \brief initialize the vector indexed by the geometry type (for each
//! quadrature type)
static void initGeometryTypeVector(GeometryTypeVector *gtv)
*gtv = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
class QuadratureRules
//! type of quadrature rule provided by this cache
using QuadratureRule = Dune::QuadratureRule<ctype, dim>;
//! real rule creator
DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
//! identifier for quadrature rules, used in the cache container
struct QuadratureKey
unsigned int id; // topologyId
int p; // order
QuadratureType::Enum qt; // quadrature type
static std::vector<std::pair< // indexed by quadrature type
> > quadratureCache(QuadratureType::size);
auto & quadratureTypeLevel = quadratureCache[qt];
std::call_once(quadratureTypeLevel.first, initGeometryTypeVector,
auto & geometryTypeLevel =
std::call_once(geometryTypeLevel.first, initQuadratureOrderVector,
&geometryTypeLevel.second, qt, t);
// we only have one quadrature rule for points
auto & quadratureOrderLevel = geometryTypeLevel.second[dim == 0 ? 0 : p];
std::call_once(quadratureOrderLevel.first, initQuadratureRule,
&quadratureOrderLevel.second, qt, t, p);
friend bool operator< (const QuadratureKey& lhs, const QuadratureKey& rhs)
return std::tie(,lhs.p,lhs.qt) < std::tie(,rhs.p,rhs.qt);
return quadratureOrderLevel.second;
//! singleton provider
DUNE_EXPORT static QuadratureRules& instance()
static QuadratureRules instance;
return instance;
//! private constructor
QuadratureRules () {}
//! maximum quadrature order for given geometry type and quadrature type
static unsigned
maxOrder(const GeometryType& t,
QuadratureType::Enum qt=QuadratureType::GaussLegendre)
static unsigned maxOrder (const GeometryType& t,
QuadratureType::Enum qt = QuadratureType::GaussLegendre)
return QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt);
//! select the appropriate QuadratureRule for GeometryType t and order p
static const QuadratureRule& rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
static const QuadratureRule& rule (const GeometryType& t, int p,
QuadratureType::Enum qt = QuadratureType::GaussLegendre)
return instance()._rule(t,p,qt);
// Container to store the cached quadrature rules
thread_local std::map<QuadratureKey, QuadratureRule> quadCache;
//! @copydoc rule
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);
auto [it,inserted] = quadCache.try_emplace(QuadratureKey{,p,qt});
if (inserted)
it->second = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
return it->second;
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment