Numerical quadrature
====================

Element integrals (:doc:`variational_form` eq. :eq:`ke-me`) are
evaluated on the reference element :math:`\hat\Omega` after the
isoparametric mapping (:doc:`isoparametric` eq.
:eq:`jacobian`).  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 :math:`\hat\Omega \subset \mathbb{R}^{d}`
is a finite set of pairs

.. math::

    \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 :math:`\le p`, the rule's *degree of precision*.  The
post-Jacobian element integral becomes

.. math::

    \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 :math:`n`-point Gauss-Legendre rule on
:math:`[-1, +1]` is exact for polynomials of degree
:math:`2n - 1` — its abscissae are the roots of the Legendre
polynomial :math:`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:

.. list-table::
    :header-rows: 1
    :widths: 22 28 50

    * - Kernel
      - Rule
      - Degree of precision
    * - ``HEX8``
      - 2 × 2 × 2 Gauss-Legendre (8 pts)
      - exact for trilinear :math:`\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
        :math:`\mathbf{K}_e`; 14-point Irons rule for
        :math:`\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 :math:`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 :math:`\le 5` exactly on the cube
:math:`[-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 :math:`3 \times 3 \times 3`
rule for the consistent-mass integral
:math:`\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 :math:`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 :math:`\le 2` in volume
coordinates — sufficient for the quadratic-tet
:math:`\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
:math:`(\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 :math:`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 ×
:math:`[-1, +1]`.  The triangle integrals use the
Dunavant (1985) family — symmetric cubatures with rational
weights on the unit triangle :math:`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 :math:`\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
:math:`\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 :math:`(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
:math:`\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
:math:`\int E I\, \kappa^{\!\top}\, \kappa\, \mathrm{d}x`
collapse to the canonical
:math:`(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 :math:`n_q` points.
2. At each point, computes the shape functions
   :math:`\mathbf{N}(\boldsymbol{\xi}_q)` and their
   reference-frame derivatives
   :math:`\partial \mathbf{N}/\partial \boldsymbol{\xi}_q`.
3. Forms :math:`\mathbf{J}(\boldsymbol{\xi}_q)` from the nodal
   geometry, then :math:`\mathbf{B}(\boldsymbol{\xi}_q)` via
   :eq:`dn-dx`.
4. Accumulates :math:`w_q\, \det \mathbf{J}(\boldsymbol{\xi}_q)
   \, \mathbf{B}^{\!\top}\, \mathbf{C}\, \mathbf{B}` into the
   element stiffness (resp.
   :math:`w_q\, \det \mathbf{J}(\boldsymbol{\xi}_q)\, \rho\,
   \mathbf{N}^{\!\top}\, \mathbf{N}` for mass).

The Python kernels in :mod:`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).