General concepts » Quadrature formulas

Exactness ranges and default quadrature dispatch.

Contract

A quadrature formula on a reference polytope $ \widehat K $ has points $ \widehat x_q $ and weights $ w_q $ such that

\[ \int_{\widehat K} f(\widehat x)\,d\widehat x \simeq \sum_q w_q f(\widehat x_q). \]

In Rodin, the requested order $ p $ denotes the polynomial degree of exactness, not the number of quadrature points. The default dispatcher returns a rule exact on $ \mathbb P_p(\widehat K) $ with positive weights and points inside $ \widehat K $ .

This exactness choice is part of the assembled mathematical problem, not just a performance setting. See Numerical contracts for how quadrature interacts with interpolation, projection, manufactured solutions, and specialized integrators.

Default dispatch

Reference polytopeTabulated ruleFallback
PointTrivial ruleNone
SegmentGauss–LegendreSame rule
TriangleXiao–Gimbutas, degrees 1–50Positive conical product
QuadrilateralGauss–Legendre tensor productSame rule
TetrahedronXiao–Gimbutas, degrees 1–15Positive conical product
HexahedronGauss–Legendre tensor productSame rule
WedgeXG triangle times GL, degrees 1–50Conical triangle times GL
PyramidWitherden–Vincent, degrees 1–10Positive conical product

The complete directly constructible table ranges are:

FamilyTriangleQuadrilateralTetrahedronHexahedronWedgePyramid
Xiao–Gimbutas1–50–1–15–––
Witherden–Vincent1–201–201–101–101–101–10

The Xiao–Gimbutas [14] and Witherden–Vincent [13] rules are published positive-interior formulas. Rodin vendors their coefficients with the corresponding attribution licenses; no numerical search occurs at runtime.

Conical-product fallback

Let

\[ n_p=\left\lceil\frac{p+1}{2}\right\rceil. \]

Each one-dimensional Gaussian factor uses $ n_p $ points. The simplex rules follow the Duffy transformation [4] and the associated conical-product construction [11].

For the triangle,

\[ (r,s)=\bigl(u,(1-u)v\bigr), \qquad dr\,ds=(1-u)\,du\,dv, \]

so the rule is a Gauss–Jacobi rule with weight $ 1-u $ times a Gauss–Legendre rule. For the tetrahedron,

\[ (r,s,t)=\bigl(u,(1-u)v,(1-u)(1-v)w\bigr), \]

and the Jacobian $ (1-u)^2(1-v) $ is carried by two Gauss–Jacobi factors. For the pyramid,

\[ (x,y,z)=\bigl((1-z)u,(1-z)v,z\bigr), \]

and the Jacobian $ (1-z)^2 $ is carried by the vertical Gauss–Jacobi factor.

These fallbacks are deterministic, positive, interior, and available at arbitrary order. They are not fully symmetric and can use more points than optimized tabulated formulas.

API and caching

using namespace Rodin;

const auto& cached = QF::PolytopeQuadratureFormula::get(
  8, Geometry::Polytope::Type::Tetrahedron);

QF::PolytopeQuadratureFormula owned(
  8, Geometry::Polytope::Type::Tetrahedron);

PolytopeQuadratureFormula::get returns a process-wide canonical formula. Repeated per-thread requests are served by a small thread-local cache. PolytopeQuadratureFormula::build and the constructor return owned formulas. XiaoGimbutas, WitherdenVincent, GaussLegendre, GaussLobatto, and GrundmannMoller can also be selected explicitly. Grundmann–Moller is not a default because its weights are signed above degree one.

Curved elements

The rule order applies on the reference polytope. On a physical cell, the mapped integrand also contains the determinant of the geometric Jacobian. Non-affine transformations can therefore require an order above the degree suggested by the finite-element fields alone. A quadrature point must remain inside the reference polytope because a curved transformation is not required to extend beyond it.

See Also