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:

\[\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:

\[\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

\(\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

\(\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 \(\mathbf{K} \boldsymbol{\phi} = \omega^2 \mathbf{M} \boldsymbol{\phi}\) satisfy:

\[\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

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

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	\(\delta_{ij}\) / \(\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: tests/analytical/test_conservation.py — landed with TA-9c (#146) earlier this session.