Quadrature formulas
Exactness ranges and default quadrature dispatch.
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 Numericalcontracts" for how quadrature interacts with interpolation, projection,
manufactured solutions, and specialized integrators.
@section guides-quadrature-dispatch Default dispatch
<table>
<tr><th>Reference polytope</th><th>Tabulated rule</th><th>Fallback</th></tr>
<tr><td>Point</td><td>Trivial rule</td><td>None</td></tr>
<tr><td>Segment</td><td>Gauss–Legendre</td><td>Same rule</td></tr>
<tr><td>Triangle</td><td>Xiao–Gimbutas, degrees 1–50</td><td>Positive conical product</td></tr>
<tr><td>Quadrilateral</td><td>Gauss–Legendre tensor product</td><td>Same rule</td></tr>
<tr><td>Tetrahedron</td><td>Xiao–Gimbutas, degrees 1–15</td><td>Positive conical product</td></tr>
<tr><td>Hexahedron</td><td>Gauss–Legendre tensor product</td><td>Same rule</td></tr>
<tr><td>Wedge</td><td>XG triangle times GL, degrees 1–50</td><td>Conical triangle times GL</td></tr>
<tr><td>Pyramid</td><td>Witherden–Vincent, degrees 1–10</td><td>Positive conical product</td></tr>
</table>
The complete directly constructible table ranges are:
<table>
<tr><th>Family</th><th>Triangle</th><th>Quadrilateral</th><th>Tetrahedron</th><th>Hexahedron</th><th>Wedge</th><th>Pyramid</th></tr>
<tr><td>Xiao–Gimbutas</td><td>1–50</td><td>–</td><td>1–15</td><td>–</td><td>–</td><td>–</td></tr>
<tr><td>Witherden–Vincent</td><td>1–20</td><td>1–20</td><td>1–10</td><td>1–10</td><td>1–10</td><td>1–10</td></tr>
</table>
The Xiao–Gimbutas @cite xiao2010numerical and Witherden–Vincent
@cite witherden2015identification rules are published positive-interior
formulas.
Rodin vendors their coefficients with the corresponding attribution licenses;
no numerical search occurs at runtime.
@section guides-quadrature-fallback Conical-product fallback
Let
\iline 59 \iline 60 \_form#2198@_fakenl@_fakenl
Each one-dimensional Gaussian factor uses \_form#2199 points. The simplex
rules follow the Duffy transformation @cite duffy1982quadrature and the
associated conical-product construction @cite stroud1971approximate.
For the triangle,
\iline 67 \iline 68 \iline 69 \_form#2200@_fakenl@_fakenl@_fakenl
so the rule is a Gauss–Jacobi rule with weight \_form#2201 times a
Gauss–Legendre rule. For the tetrahedron,
\iline 73 \iline 74 \_form#2202@_fakenl@_fakenl
and the Jacobian \_form#592 is carried by two Gauss–Jacobi
factors. For the pyramid,
\iline 78 \iline 79 \_form#2203@_fakenl@_fakenl
and the Jacobian \_form#595 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.
@section guides-quadrature-api API and caching
@code{.cpp}
using namespace Rodin;
const auto& cached = QF::PolytopeQuadratureFormula::get(
8, Geometry::Polytope::Type::Tetrahedron);
QF::PolytopeQuadratureFormula owned(
8, Geometry::Polytope::Type::Tetrahedron);
@endcode
<tt>PolytopeQuadratureFormula::get</tt> returns a process-wide canonical formula.
Repeated per-thread requests are served by a small thread-local cache.
<tt>PolytopeQuadratureFormula::build</tt> and the constructor return owned
formulas. <tt>XiaoGimbutas</tt>, <tt>WitherdenVincent</tt>, <tt>GaussLegendre</tt>,
<tt>GaussLobatto</tt>, and <tt>GrundmannMoller</tt> can also be selected explicitly.
Grundmann–Moller is not a default because its weights are signed above
degree one.
@section guides-quadrature-curved 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.
@section guides-quadrature-see-also See Also
- @ref guides-numerical-contracts "Numerical contracts"
- @ref guides-variational "Variational formulations"
@xmlonly <mcss:footernavigation xmlns:mcss="http:/