Numerical quadrature#

Element integrals (Variational form and the discretised equations eq. (3)) are evaluated on the reference element \(\hat\Omega\) after the isoparametric mapping (Isoparametric mapping eq. (2)). Each kernel uses a Gauss-style quadrature rule chosen for the shape-function family of that element: enough points to integrate the bilinear form :math:`mathbf{B}^{!top}mathbf{C},mathbf{B}` exactly on the constant-Jacobian limit, and just enough to suppress numerical locking on a curved element. This chapter walks the rules in use across femorph-solver’s element kernels.

Notation#

A quadrature rule on \(\hat\Omega \subset \mathbb{R}^{d}\) is a finite set of pairs

\[\mathcal{Q} = \{ (\boldsymbol{\xi}_q,\, w_q) \}_{q = 1}^{n_q}, \qquad \int_{\hat\Omega} f(\boldsymbol{\xi}) \,\mathrm{d}\boldsymbol{\xi} \;\approx\; \sum_{q = 1}^{n_q} w_q\, f(\boldsymbol{\xi}_q),\]

with the property that the sum is exact for polynomials of degree \(\le p\), the rule’s degree of precision. The post-Jacobian element integral becomes

\[\int_{\Omega_e} f(\mathbf{x}) \,\mathrm{d}V \;=\; \int_{\hat\Omega} f(\mathbf{x}(\boldsymbol{\xi}))\, \det \mathbf{J}(\boldsymbol{\xi}) \,\mathrm{d}\boldsymbol{\xi} \;\approx\; \sum_{q = 1}^{n_q} w_q\, f(\mathbf{x}(\boldsymbol{\xi}_q))\, \det \mathbf{J}(\boldsymbol{\xi}_q).\]

Tensor-product Gauss-Legendre on the reference cube / square#

The 1D \(n\)-point Gauss-Legendre rule on \([-1, +1]\) is exact for polynomials of degree \(2n - 1\) — its abscissae are the roots of the Legendre polynomial \(P_n\) and the weights follow from the standard Lagrange-interpolant integral (Davis & Rabinowitz 1984 §2.7). Tensor-product extensions land directly on the reference hypercube:

Kernel

Rule

Degree of precision

HEX8

2 × 2 × 2 Gauss-Legendre (8 pts)

exact for trilinear \(\mathbf{B}^{\!\top} \mathbf{C}\, \mathbf{B}\) on a constant-Jacobian element (Zienkiewicz & Taylor Table 5.3).

HEX20

2 × 2 × 2 Gauss-Legendre (8 pts) for \(\mathbf{K}_e\); 14-point Irons rule for \(\mathbf{M}_e\)

reduced (Bedrosian 1992) on stiffness — same 8 points as HEX8. See Irons rule (HEX20 mass) below for the mass rule.

QUAD4_PLANE

2 × 2 Gauss-Legendre (4 pts)

exact bilinear on a regular quad (Hughes 2000 §3.6).

QUAD4_SHELL

2 × 2 Gauss + 1 × 1 selective-reduced shear

Malkus–Hughes 1978 cures shear-locking in thin-bending (CMAME 1978).

Reduced integration on HEX20 (8 points instead of the full 27-point \(3 \times 3 \times 3\) rule) is the standard practitioner choice for the serendipity hex (Bedrosian 1992) — the bending-mode response stays exact while spurious hourglass modes are absent because the 8 points still couple all 60 stiffness rows.

Irons rule (HEX20 mass)#

The 14-point Irons (1971) rule integrates a polynomial of degree \(\le 5\) exactly on the cube \([-1, +1]^{3}\) and is not a tensor product — Irons constructed it by hand as a degree-5 cubature with the minimum point count compatible with the cube’s symmetry group. It is preferred over the 27-point \(3 \times 3 \times 3\) rule for the consistent-mass integral \(\int \rho\, \mathbf{N}^{\!\top}\, \mathbf{N}\,\mathrm{d}V\) on HEX20 because it preserves the mass-matrix sparsity pattern without overshooting the polynomial degree of the integrand.

  • Reference: Irons, B. M. (1971) “Quadrature rules for brick based finite elements,” International Journal for Numerical Methods in Engineering 3 (2), 293–294.

Keast rule (TET10)#

The Keast (1986) family supplies symmetric cubatures on the unit tetrahedron with explicit rationality of the abscissae in volume coordinates \(L_0 + L_1 + L_2 + L_3 = 1\). femorph-solver’s TET10 kernel uses the 4-point Keast rule, exact for polynomials of degree \(\le 2\) in volume coordinates — sufficient for the quadratic-tet \(\mathbf{B}^{\!\top}\, \mathbf{C}\, \mathbf{B}\) integrand under a constant-Jacobian map.

The four points sit on the lines from each corner to the centroid at fractional coordinate \((\tfrac{5 - \sqrt{5}}{20},\, \tfrac{5 - \sqrt{5}}{20},\, \tfrac{5 - \sqrt{5}}{20},\, \tfrac{5 + 3\sqrt{5}}{20})\) (and symmetric permutations); each weight is \(1/24\).

  • Reference: Keast, P. (1986) “Moderate-degree tetrahedral quadrature formulas,” Computer Methods in Applied Mechanics and Engineering 55 (3), 339–348.

  • Reference: Stroud, A. H. (1971) Approximate Calculation of Multiple Integrals, Prentice-Hall — original derivations of low-order tet rules.

Triangle rules — Dunavant (WEDGE15)#

WEDGE15’s reference element factors as a triangle × \([-1, +1]\). The triangle integrals use the Dunavant (1985) family — symmetric cubatures with rational weights on the unit triangle \(L_0 + L_1 + L_2 = 1\). femorph-solver’s wedge stiffness uses a 3-point rule (degree-2 exact) tensored with a 3-point Gauss-Legendre rule (degree-5 exact in \(\zeta\)), giving a 9-point composite. The mass uses a 7-point triangle rule × 3-pt Legendre = 21 pts, matching the higher polynomial degree of \(\mathbf{N}^{\!\top}\, \mathbf{N}\).

  • Reference: Dunavant, D. A. (1985) “High degree efficient symmetrical Gaussian quadrature rules for the triangle,” International Journal for Numerical Methods in Engineering 21 (6), 1129–1148.

Pyramid rule — Duffy (PYR13)#

The pyramid is the most awkward of the standard linear-elastic reference shapes because the pyramidal limit of a hex is rational, not polynomial — Bedrosian’s 1992 pyramid shape functions explicitly carry a \((1 - \zeta)^{-1}\) factor. The Duffy (1982) transformation maps the pyramid to a collapsed cube (one vertex of the apex maps to a face of the cube) and the resulting integrand is finite at every Gauss point. femorph-solver’s PYR13 kernel uses a 2 × 2 × 2 Duffy collapsed-hex rule (8 points) on the transformed cube — degree \(\le 3\) exact in the rational basis.

  • Reference: Duffy, M. G. (1982) “Quadrature over a pyramid or cube of integrands with a singularity at a vertex,” SIAM Journal on Numerical Analysis 19 (6), 1260–1262.

  • Reference: Bedrosian, G. (1992) “Shape functions and integration formulas for three-dimensional finite element analysis,” International Journal for Numerical Methods in Engineering 35 (1), 95–108 — original PYR13 / WEDGE15 serendipity shape functions.

BEAM2 — analytic (no quadrature)#

BEAM2 (slender Euler-Bernoulli) integrates its 12 × 12 stiffness blocks in closed form (Cook Tables 2.5 and 16.3-1). The Hermite cubic basis lets the bilinear form \(\int E I\, \kappa^{\!\top}\, \kappa\, \mathrm{d}x\) collapse to the canonical \((EI / L^{3})\, [12, 6L, -12, 6L; \dots]\) block — no quadrature error to track, no integration points to choose.

How femorph-solver evaluates a quadrature rule#

Each element kernel’s ke / me / eel_batch routine:

  1. Iterates over the rule’s \(n_q\) points.

  2. At each point, computes the shape functions \(\mathbf{N}(\boldsymbol{\xi}_q)\) and their reference-frame derivatives \(\partial \mathbf{N}/\partial \boldsymbol{\xi}_q\).

  3. Forms \(\mathbf{J}(\boldsymbol{\xi}_q)\) from the nodal geometry, then \(\mathbf{B}(\boldsymbol{\xi}_q)\) via (3).

  4. Accumulates \(w_q\, \det \mathbf{J}(\boldsymbol{\xi}_q) \, \mathbf{B}^{\!\top}\, \mathbf{C}\, \mathbf{B}\) into the element stiffness (resp. \(w_q\, \det \mathbf{J}(\boldsymbol{\xi}_q)\, \rho\, \mathbf{N}^{\!\top}\, \mathbf{N}\) for mass).

The Python kernels in femorph_solver.elements each embed the rule’s abscissae / weights as module-level constants so the loop runs through nanobind into a C++ batch routine; the equations are still those above.

References#

  • Davis, P. J. and Rabinowitz, P. (1984) Methods of Numerical Integration, 2nd ed., Academic Press, §2.7 (Gauss-Legendre).

  • Stroud, A. H. (1971) Approximate Calculation of Multiple Integrals, Prentice-Hall (foundational reference for tetrahedron / pyramid cubature).

  • Zienkiewicz, O. C. and Taylor, R. L. (2013) The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, §5 (numerical integration in element formulation), Table 5.3 (Gauss-Legendre weights).

  • Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §6.7 (numerical integration).

  • Hughes, T. J. R. (2000) The Finite Element Method — Linear Static and Dynamic Finite Element Analysis, Dover, §3.6 (Gauss-Legendre on the bilinear quad).

  • Irons, B. M. (1971) “Quadrature rules for brick based finite elements,” IJNME 3 (2), 293–294.

  • Keast, P. (1986) “Moderate-degree tetrahedral quadrature formulas,” CMAME 55 (3), 339–348.

  • Dunavant, D. A. (1985) “High degree efficient symmetrical Gaussian quadrature rules for the triangle,” IJNME 21 (6), 1129–1148.

  • Duffy, M. G. (1982) “Quadrature over a pyramid or cube of integrands with a singularity at a vertex,” SIAM J. Numer. Anal. 19 (6), 1260–1262.

  • Bedrosian, G. (1992) “Shape functions and integration formulas for three-dimensional finite element analysis,” IJNME 35 (1), 95–108.

  • Malkus, D. S. and Hughes, T. J. R. (1978) “Mixed finite element methods — reduced and selective integration techniques: a unification of concepts,” CMAME 15 (1), 63–81 (selective-reduced shear for shells).