General concepts » Quadrature formulas

Contract

A quadrature formula on a reference polytope has points and weights such that

In Rodin, the requested order denotes the polynomial degree of exactness, not the number of quadrature points. The default dispatcher returns a rule exact on with positive weights and points inside .

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.

Formulas, mapped points, and basis tables

QuadratureFormulaBase provides read-only access to reference points and weights. Concrete families such as GaussLegendre can be constructed locally for immediate calculations. The interface carries no polymorphic-copy or cache-identity contract.

QuadratureFormula owns immutable reference samples retained until process termination. Its factory constructs a concrete family on demand:

const auto& formula = QF::QuadratureFormula::get<QF::GaussLegendre>(
  geometry, size_t{4});

The family and constructor argument value types select a shared pool. Integer and enum values, and references to published formulas, form its tuple key. The pool uses FlatMap with ordinary tuple comparison; formula references use total pointer ordering internally. A bounded thread-local lookup avoids repeated locking. Construction occurs outside the pool mutex; insertion publishes one immutable object even when several threads miss concurrently. Each pool entry constructs its concrete family and published interface in one allocation. The interface reads the retained rule directly; no second point or weight array is created. Optional read-only sample views provide direct point access and strided weight access, including uniform weights, without per-sample virtual dispatch. Pointwise-only families remain supported. Generators must own their samples or borrow only storage with process lifetime, must not retain caller-owned constructor arguments, and must be determined by their arguments and safe to construct concurrently. Counts and orders should be passed as size_t.

QuadratureFormula selects a family from geometry and requested exactness. Its get operation, and get, return published formulas through the same factory. Changing an integrator's order chooses another formula; it never modifies an existing formula.

PolytopeQuadrature owns points mapped to one physical polytope. It accepts only a published formula, so its borrowed source cannot expire. Mesh-owned quadratures remain valid until geometry flush; assembly rules instead own quadratures for their currently bound polytope.

ReferenceTabulation stores basis values and reference derivatives at those reference samples. Persistent table lookup accepts only a published formula. Tables belong to a bounded thread-local cache: formula permanence does not extend a borrowed table's lifetime across eviction. Immediate element.tabulate(table, localFormula) accepts a caller-owned QuadratureFormulaBase and leaves the result in caller-owned numerical storage.

QuadratureRule binds the integrand and physical polytope, selects the reference formula, and consumes mapped points and basis values during integration. Geometry-specific quantities remain separate from reusable reference-element values. A copied geometric point owns its reference coordinates and may retain a permanent formula pointer and index for lazy coordinate evaluation; if its basis table is absent, direct evaluation applies.

Default dispatch

The complete directly constructible table ranges are:

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

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

For the triangle,

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

and the Jacobian is carried by two Gauss–Jacobi factors. For the pyramid,

and the Jacobian 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::QuadratureFormula::get(
  8, Geometry::Polytope::Type::Tetrahedron);

const auto& explicitFamily = QF::QuadratureFormula::get<QF::XiaoGimbutas>(
  size_t{8}, Geometry::Polytope::Type::Tetrahedron);

QuadratureFormula::get returns a process-wide canonical formula. Repeated per-thread requests are served by a small thread-local cache. Published formulas are immutable and cannot be copied or moved. Concrete family constructors remain available for locally 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

• Numerical contracts
• Variational formulations