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 Variational form and the discretised equations and Isoparametric mapping first, then visit the per-topic chapters from there.

Theory chapters

Variational form and the discretised equations#

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

\[\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 \(\boldsymbol{\varepsilon}=\tfrac{1}{2}(\nabla\mathbf{u}+\nabla\mathbf{u}^{\!\top})\) the symmetric gradient and \(\mathbf{C}\) the fourth-order elasticity tensor. After discretisation \(\mathbf{u}\approx\mathbf{N u}_e\) and \(\boldsymbol{\varepsilon}=\mathbf{B u}_e\) the element stiffness and mass take the standard form

\[\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 Variational form and the discretised equations 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 Isoparametric mapping 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:

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

\[\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 \(\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 \(N\) identical sectors, the global spectrum decomposes per harmonic index \(k \in \{0, 1, \dots, N/2\}\). Each harmonic yields a reduced complex-Hermitian problem \(\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 Cyclic-symmetry modal analysis 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).

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.