Theory
======

FEM foundations used throughout femorph-solver.  Every formula
below is drawn from the open literature — Bathe, Zienkiewicz–Taylor,
Cook, Hughes, Belytschko, Timoshenko, Keast, Irons, Simo–Rifai — and
the per-kernel docstrings cite the specific section / equation the
code implements.

The chapters below expand the high-level summary that follows
into standalone derivations.  New readers should walk
:doc:`variational_form` and :doc:`isoparametric` first, then
visit the per-topic chapters from there.

.. toctree::
    :maxdepth: 1
    :caption: Theory chapters

    variational_form
    isoparametric
    quadrature
    assembly
    bc_elimination
    eigenproblem

Variational form and the discretised equations
----------------------------------------------

For a linear-elastic body :math:`\Omega` with displacement field
:math:`\mathbf{u}`, the principle of virtual work reduces to

.. math::

    \int_{\Omega} \boldsymbol{\varepsilon}(\delta\mathbf{u})^{\!\top}
    \, \mathbf{C}\, \boldsymbol{\varepsilon}(\mathbf{u})\,\mathrm{d}V
    \;=\;
    \int_{\Omega} \delta\mathbf{u}^{\!\top} \mathbf{b}\,\mathrm{d}V
    +
    \int_{\partial\Omega_t} \delta\mathbf{u}^{\!\top}\mathbf{t}\,\mathrm{d}A,

with :math:`\boldsymbol{\varepsilon}=\tfrac{1}{2}(\nabla\mathbf{u}+\nabla\mathbf{u}^{\!\top})`
the symmetric gradient and :math:`\mathbf{C}` the fourth-order elasticity
tensor.  After discretisation :math:`\mathbf{u}\approx\mathbf{N u}_e`
and :math:`\boldsymbol{\varepsilon}=\mathbf{B u}_e` the element
stiffness and mass take the standard form

.. math::

    \mathbf{K}_e = \int_{\Omega_e} \mathbf{B}^{\!\top}\mathbf{C}\mathbf{B}\,\mathrm{d}V,
    \qquad
    \mathbf{M}_e = \int_{\Omega_e} \rho\,\mathbf{N}^{\!\top}\mathbf{N}\,\mathrm{d}V.

**Reference:** Zienkiewicz & Taylor, *The Finite Element Method:
Its Basis and Fundamentals*, 7th ed., 2013, §2 + §8; Bathe, *Finite
Element Procedures*, 2nd ed., 2014, §4.  See
:doc:`variational_form` for the full derivation.

Isoparametric mapping
---------------------

Every solid element in femorph-solver is isoparametric: the same
shape functions that interpolate displacement interpolate the
geometry, so ``J(ξ) = ∂x/∂ξ = ∑ᵢ (∂Nᵢ/∂ξ) xᵢ`` and derivatives
of ``N`` with respect to physical coordinates use ``dN/dx = J⁻ᵀ
dN/dξ``.  The integrals above transform to the reference element
with ``dV = det(J) dξ`` and are evaluated by Gauss quadrature.

**Reference:** Zienkiewicz & Taylor §4 + §5.  Cook, Malkus, Plesha,
Witt, *Concepts and Applications of FEA*, 4th ed., 2002, §6.  See
:doc:`isoparametric` for the Jacobian transform, distortion
measures, and per-kernel reference-element table.

Quadrature rules in use
-----------------------

Each element's Gauss rule is fixed by its kernel module:

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

    * - Element kernel
      - Rule
      - Source
    * - ``HEX8``
      - 2×2×2 Gauss-Legendre
      - Zienkiewicz & Taylor Table 5.3
    * - ``HEX20``
      - 2×2×2 reduced (K) + 14-pt Irons (M)
      - Bedrosian 1992 (K); Irons 1971 (M)
    * - ``TET10``
      - 4-point Keast rule
      - Keast, CMAME 1986
    * - ``WEDGE15``
      - 9-pt (K) + 21-pt (M) triangle × line
      - Dunavant 1985 triangle × Gauss-Legendre
    * - ``PYR13``
      - 2×2×2 Duffy collapsed-hex
      - Duffy 1982; Bedrosian 1992 shapes
    * - ``QUAD4_PLANE`` / ``QUAD4_SHELL``
      - 2×2 Gauss + 1×1 selective-reduced shear (shell)
      - Malkus-Hughes, CMAME 1978
    * - ``BEAM2``
      - Analytic (Hermite cubic closed-form)
      - Cook Tables 2.5 / 16.3-1

Assembly and boundary-condition elimination
-------------------------------------------

Global assembly scatters per-element ``Kₑ`` / ``Mₑ`` blocks into a
shared CSR pattern.  Prescribed DOFs are eliminated by
partitioning — the free / fixed partitioning of ``K`` produces the
reduced system ``K_ff u_f = f_f − K_fc u_c`` that the linear
solver sees.  No Lagrange multipliers are used for simple
Dirichlet BCs.

**Reference:** Saad, *Iterative Methods for Sparse Linear Systems*,
2nd ed., 2003, §3.3–§3.4 (assembly / CSR); Cook §2.10; Bathe §3.4.

Eigenvalue problem and shift-invert
-----------------------------------

Modal analysis solves the generalised symmetric-definite problem

.. math::

    \mathbf{K}\boldsymbol{\phi} = \omega^{2} \mathbf{M}\boldsymbol{\phi}.

femorph-solver uses a shift-invert Lanczos iteration (ARPACK
``eigsh`` with ``sigma=0`` by default) to extract the lowest
``n_modes``.  Spectral clustering near :math:`\sigma` is resolved by
re-factoring ``(K - σM)`` through the registered linear backend
(Pardiso / CHOLMOD / MUMPS / SuperLU auto-chain).

**Reference:** Parlett, *The Symmetric Eigenvalue Problem*, SIAM,
1998, §13.  Lehoucq, Sorensen, Yang, *ARPACK Users' Guide*, 1998.
Ericsson & Ruhe, *Math. Comp.* 35 (1980), pp. 1251–1268
(shift-invert).

Cyclic-symmetry reduction
-------------------------

For a rotor of :math:`N` identical sectors, the global spectrum
decomposes per harmonic index :math:`k \in \{0, 1, \dots, N/2\}`.
Each harmonic yields a reduced complex-Hermitian problem
:math:`\mathbf{K}_k\boldsymbol{\phi}_k = \omega^{2}
\mathbf{M}_k\boldsymbol{\phi}_k` on a single base sector; the real
``2n``-augmentation identity lets scipy's real-symmetric eigsh
handle it.  See :doc:`/user-guide/solving/cyclic` for the driving
API.

**Reference:** Thomas, *J. Sound Vib.* 66 (1979), pp. 585–597;
Wildheim, *J. Appl. Mech.* 46 (1979), pp. 891–893; Grimes, Lewis,
Simon, *SIAM J. Matrix Anal. Appl.* 15 (1994) (real 2n
augmentation).

.. _references-by-kernel:

Where the cites live
--------------------

Each element / solver module carries a top-of-file ``References``
section with the full citations — section numbers, edition,
publication year — for every formula, shape function,
integration rule, and algorithm variant it implements.  The
reference pages here summarise; the docstrings are authoritative.

* :doc:`../elements/index` — element kernels + their cites.
* :doc:`../solvers/index` — solver algorithms + their cites.