Conservation + orthogonality — four-in-one sanity check
========================================================

Four fundamental correctness properties every structural FE code
must satisfy — grouped here because they share the same fixture
(a 2×2×2 HEX8 unit cube) and the same "identity, not
approximation" posture of the patch test.  If any of the four
fails, the formulation has a bug upstream of any convergence
analysis.

1. Rigid-body modes
-------------------

An unconstrained 3D elastic solid has a six-dimensional null
space: three translations + three rotations.  The modal solve
must recover:

.. math::

    \omega_1^2 = \omega_2^2 = \cdots = \omega_6^2 = 0

to within the eigensolver's zero-frequency floor (set by the
``K``-normalisation used internally).  Modes 7+ must be strictly
positive — first elastic mode sits orders of magnitude above the
floor.

Reference: ``ker(K) = RBM`` for any symmetric PSD 3D-elasticity
operator.  This follows directly from translation and rotation
invariance of the strain-energy functional.

2. Energy balance (Clapeyron identity)
--------------------------------------

For a linear-elastic static solve the strain-energy stored
equals the external-work done:

.. math::

    \tfrac{1}{2} \mathbf{u}^T \mathbf{K} \mathbf{u}
        = \tfrac{1}{2} \mathbf{f}^T \mathbf{u}

This is a consequence of Clapeyron's theorem applied to linear-
elastic systems.  Any deviation beyond floating-point residuals
flags a stiffness-matrix assembly bug or a load-vector bug —
you can't tell from the residual alone which side is wrong, so
this test doubles as a combined sanity check on both.

Reference: Cook et al. *Concepts and Applications of FEA*, §2.2;
Zienkiewicz / Taylor §1.3.

3. Mass conservation
--------------------

Two conditions, one per mass matrix formulation:

**Lumped mass**
  :math:`\text{trace}(\mathbf{M}_{\text{lumped}}) = 3 \rho V`.
  Each diagonal entry carries one translational DOF's share;
  the three-dimensionality multiplier is exact.

**Consistent mass**
  :math:`\mathbf{1}^T \mathbf{M} \mathbf{1} = 3 \rho V`.
  Follows from direction-separability of the consistent-mass
  block at a node — each Cartesian direction carries ρV
  independently, sharing no stiffness coupling.

Reference: Zienkiewicz/Taylor Vol. 1 §20; Hinton/Rock/Zienkiewicz
*EESD* 1976 (HRZ row-sum lumping derivation).

4. Modal orthogonality
----------------------

Eigenvectors of the generalised eigenproblem
:math:`\mathbf{K} \boldsymbol{\phi} = \omega^2 \mathbf{M} \boldsymbol{\phi}`
satisfy:

.. math::

    \boldsymbol{\Phi}^T \mathbf{M} \boldsymbol{\Phi} = \mathbf{I}
    \quad \text{and} \quad
    \boldsymbol{\Phi}^T \mathbf{K} \boldsymbol{\Phi} = \text{diag}(\omega^2)

when ARPACK or LAPACK returns M-normalised modes — which all
femorph-solver eigen backends (ARPACK via ``eigsh``, LAPACK via
``eigh``, PRIMME, LOBPCG) do by default.

Reference: Hughes, *The Finite Element Method*, Dover, §9.3;
Bathe, *Finite Element Procedures*, §10.2.

Fixture
-------

One 2 × 2 × 2 HEX8 unit cube.  Aluminium-ish properties
(``E = 70 GPa``, ``ν = 0.3``, ``ρ = 2700 kg/m³``) so ``ρV = 2700
kg`` is a clean round figure in every conservation assertion.
81 DOFs total → the dense LAPACK modal path handles the
singular unconstrained pencil cleanly (the default shift-invert
ARPACK can't factor an unconstrained PSD ``K``).

femorph-solver result
---------------------

.. list-table::
    :header-rows: 1
    :widths: 45 30 25

    * - Check
      - Expected
      - femorph-solver
    * - 6 rigid-body modes at ω² ≈ 0
      - Identity
      - ≤ 1e-6 × max(\|diag K\|)
    * - 7th mode strictly > 0
      - :math:`\omega_7^2 > 10^3 \omega_{RBM}`
      - ✓
    * - :math:`\tfrac{1}{2} u^T K u = \tfrac{1}{2} f^T u`
      - Identity
      - ≤ 1e-10 relative
    * - :math:`\text{trace}(M_{\text{lumped}}) = 3\rho V`
      - Identity
      - ≤ 1e-10 relative
    * - :math:`\mathbf{1}^T M \mathbf{1} = 3\rho V`
      - Identity
      - ≤ 1e-10 relative
    * - :math:`\Phi^T M \Phi = I`
      - Identity
      - ≤ 1e-8 absolute
    * - :math:`\Phi^T K \Phi = \text{diag}(\omega^2)`
      - Identity
      - ≤ 1e-6 relative

All five assertions pass at the tolerances above.  These are
identity tests — a tolerance above 1e-4 would flag a bug in K
assembly, M assembly, or eigensolve normalisation; in every case
the tolerance is tight enough to catch a bug but loose enough
for machine arithmetic noise on the Gauss-integration path.

Cross-references
----------------

.. list-table::
    :header-rows: 1
    :widths: 50 15 35

    * - Source
      - Expected
      - Problem ID / location
    * - Rigid-body count: Hughes *The FEM*
      - 6
      - Dover (2000), §4.4
    * - Clapeyron theorem: Cook *CAFEA*
      - Identity
      - §2.2
    * - Mass-trace identity: Zienkiewicz/Taylor
      - 3ρV
      - Vol. 1 §20
    * - Modal orthogonality: Hughes *FEM*
      - :math:`\delta_{ij}` / :math:`\omega^2 \delta_{ij}`
      - §9.3
    * - NAFEMS *Background to Benchmarks*
      - All four identities
      - BtB-5.1 (rigid body) + BtB-5.3 (orthogonality)
    * - Abaqus Verification Manual
      - All four identities
      - AVM 1.3.5 (rigid-body recovery), AVM 1.3.7 (modal orthogonality)
    * - MAPDL Verification Manual
      - All four identities
      - VM-67 (rigid-body modes), VM-45 (modal orthogonality)

Every reference treats these as identities that hold for *any*
correct FE implementation.  They're the minimum set below which
a solver is provably wrong; femorph-solver clears all four at
machine-precision-ish tolerances, confirming the assembly +
eigensolve stack is consistent.

Source
------

Backing regression test:
:file:`tests/analytical/test_conservation.py`
— landed with TA-9c (#146) earlier this session.