API reference #

femorph_solver.solvers.static.solve_static(K: csr_array, F: ndarray, prescribed: dict[int, float], *, linear_solver: str = 'auto', thread_limit: int | None = None) → StaticResult[source]#

Solve K u = F with Dirichlet BCs at the indices in prescribed.

Parameters:

K ((N, N) scipy.sparse.csr_array)
F ((N,) float64 load vector)
prescribed (mapping {global_dof_index: prescribed_value})
linear_solver (str) – Name of a registered linear backend — see femorph_solver.solvers.linear.list_linear_solvers().
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the linear solve. None defers to the global limit (see femorph_solver.set_thread_limit()).

Solves K φ = ω² M φ for the lowest n_modes free-vibration modes with Dirichlet BCs partitioned out. DOFs with zero diagonal stiffness (for example the transverse DOFs of a pure axial truss) are dropped as rigid-body modes — the canonical free-vibration default.

For n_free <= DENSE_CUTOFF we always use the dense LAPACK path (scipy.linalg.eigh()). For larger systems the eigen backend is configurable via the eigen_solver / linear_solver kwargs — see femorph_solver.solvers.eigen and femorph_solver.solvers.linear for the registered options.

References

Bathe, K.-J. Finite Element Procedures, 2nd ed., Prentice-Hall (2014), §11. [structural eigenproblem K φ = ω² M φ and discussion of subspace / Lanczos methods]
Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, §9. [free-vibration modal analysis, M-orthonormalisation, rigid-body mode treatment]
Wilkinson, J. H. The Algebraic Eigenvalue Problem, Oxford (1965). [reference text for the symmetric generalized eigenvalue problem solved by scipy.linalg.eigh on the dense branch (LAPACK’s DSYGVD)]

class femorph_solver.solvers.modal.ModalResult(omega_sq: ndarray, frequency: ndarray, mode_shapes: ndarray, free_mask: ndarray)[source]#

Bases: object

Result of a modal solve.

omega_sq#

(n_modes,) eigenvalues ω² = (2π f)².

Type:: numpy.ndarray

frequency#

(n_modes,) cyclic frequencies f [Hz].

Type:: numpy.ndarray

mode_shapes#

(N, n_modes) DOF-indexed mode shapes — each column is one mode in the global DOF ordering produced by femorph_solver.model.Model.dof_map().

Type:: numpy.ndarray

free_mask#

(N,) bool — DOFs kept in the eigenproblem after Dirichlet reduction.

Type:: numpy.ndarray

save(path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise this modal result to a disk-backed .pv file.

Re-projects the DOF-indexed mode_shapes onto per-node (n_points, n_dof_per_node) arrays using the model’s dof_map(), then hands off to write_modal_result().

Parameters:

path (str or pathlib.Path) – Destination file (.pv extension conventional).
model (femorph_solver.model.Model) – The model whose K/M matrices this result was computed from — supplies the mesh geometry and DOF layout needed to reshape mode_shapes.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

Notes

After saving, load the file back via femorph_solver.result.read(); the resulting object is a ModalResult (disk-backed) with the same frequencies and per-mode displacements plus the full plotting / animation surface.

femorph_solver.solvers.modal.solve_modal(K: csr_array, M: csr_array, prescribed: dict[int, float], n_modes: int = 10, *, eigen_solver: str = 'arpack', linear_solver: str = 'auto', sigma: float = 0.0, tol: float = 1e-12, thread_limit: int | None = None, mixed_precision: bool | None = False, v0: ndarray | None = None) → ModalResult[source]#

Solve the generalized eigenproblem K φ = ω² M φ.

Parameters:

K (scipy.sparse.csr_array) – Global stiffness and mass matrices, same shape.
M (scipy.sparse.csr_array) – Global stiffness and mass matrices, same shape.
prescribed (mapping {global_dof_index: 0}) – Indices to drop from the eigenproblem (clamped / zero-stiffness DOFs).
n_modes (int) – Number of lowest modes to return.
eigen_solver (str) – Name of a registered eigen backend ("arpack", "lobpcg", "primme", "dense"), or "auto" to let the solver pick based on problem size. The autotune prefers PRIMME for n_modes > 20 when the package is installed (faster on the clustered eigenvalue spectra typical of rotor modal analyses); it falls back to ARPACK otherwise. See femorph_solver.solvers.eigen.list_eigen_solvers().
linear_solver (str) – Name of the registered linear backend used by ARPACK’s shift-invert step. Pass "auto" (the default) to pick the fastest available SPD direct solver (pardiso → cholmod → umfpack → superlu). See femorph_solver.solvers.linear.list_linear_solvers() for the full registry.
tol (float) – Convergence tolerance passed to the sparse eigensolver. Default 1e-12 — tight enough that the returned frequencies agree with a full machine-epsilon (tol=0) solve to at least 13 significant digits on typical modal problems, while saving ~25 % of ARPACK iterations. Pass 0.0 to force ARPACK to converge to full machine precision; pass a looser value (1e-10, 1e-8) when wall-time matters more than the last couple of digits.
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the dense and sparse eigenpaths. None (default) defers to the global limit set via femorph_solver.set_thread_limit() or the FEMORPH_SOLVER_NUM_THREADS environment variable; pass an explicit integer to override for a single call.
v0 (numpy.ndarray or None, optional) – Starting vector for ARPACK’s Lanczos iteration. Default is None — ARPACK seeds from a random vector. Pass a previous mode shape (full-system length N or reduced-system length n_free) to warm-start a parametric sweep — ARPACK’s convergence iterations typically drop 1.5–3× when the new eigenvectors are close to the old ones (e.g. material sensitivity studies, mistuning identification). Ignored by non-ARPACK eigen backends.

femorph_solver.solvers.modal.solve_modal_reduced(K_ff: csr_array, M_ff: csr_array, free_mask: ndarray, n_modes: int = 10, *, eigen_solver: str = 'arpack', linear_solver: str = 'auto', sigma: float = 0.0, tol: float = 1e-12, thread_limit: int | None = None, mixed_precision: bool | None = False, v0: ndarray | None = None) → ModalResult[source]#

Solve the eigenproblem directly on pre-reduced K_ff / M_ff.

Skips the full-matrix BC slicing that solve_modal() does, so the caller can build K_ff / M_ff via _bc_reduce(), release the original K / M (Model cache + their own refs), then call this to run the eigensolve with a peak RSS roughly K_ff + M_ff + triu + MKL factor instead of the default path’s K + M + K_ff + M_ff + triu + MKL factor.

Cyclic-symmetry modal analysis for a single base sector.

A rotor with N identical sectors spanning 360° satisfies, at each harmonic index (nodal diameter) k ∈ {0, 1, …, N/2},

u_sector_{j+1}(global) = e^{i k α} · R(α) · u_sector_j(global)

with α = 2π/N and R(α) the spatial rotation about the symmetry axis. Applied at the base-sector’s high-angle face (which is identified with sector 1’s low-angle face), this gives the constraint

u_high(global) = e^{i k α} · R(α) · u_low(global)

relating paired low-/high-face DOFs. That constraint lets a modal solve work on one sector and produce the full-rotor spectrum one k at a time — an N-times speedup relative to meshing the whole rotor.

The reduction here mirrors the standard cyclic-symmetry reduction in its complex form. Base-sector DOFs are partitioned into interior I and low-face L / high-face H sets, with a one-to-one (low-node, high-node) correspondence given by the user. High DOFs are eliminated via the constraint above; the remaining [u_I; u_L] vector satisfies a Hermitian generalised eigenproblem

K_red φ = ω² M_red φ

with K_red = P^H K P and M_red = P^H M P (K, M real SPD; P the complex constraint matrix that embeds the reduced vector back into the full one). For k=0 and k=N/2 (even N) the phase is ±1 and the reduced problem can be real; for intermediate k it is complex Hermitian.

Per-pair rotation#

Stiffness / mass are usually assembled in a global XYZ frame, which means the cyclic phase relation picks up the same rotation that maps sector j to sector j+1. Pass pair_rotation as the (d, d) matrix R(α) acting on each node pair’s d translational DOFs (3 for a 3D solid sweeping about z). Omit it (None) when face DOFs are already in a cylindrical / per-sector local frame, in which case R(α) reduces to the identity.

Harmonic-index counting#

k=0 — single real eigenmodes (in-phase family)
k=N/2 (even N) — single real eigenmodes (anti-phase family)
0 < k < N/2 — each eigenvalue corresponds to a degenerate pair in the full rotor (travelling-wave directions). The complex reduced solve returns them once per k; the full-rotor multiset is recovered by counting them twice.

References

Thomas, D. L. “Dynamics of rotationally periodic structures.” International Journal for Numerical Methods in Engineering 14 (1), 81-102 (1979). [the foundational complex-valued cyclic reduction applied here — one base sector factored at each harmonic index k recovers the full-rotor spectrum]
Wildheim, S. J. “Vibrations of rotationally periodic structures.” Ph.D. dissertation, Linköping University (1979). [companion treatment of the constraint matrix P and its Hermitian reduction P^H K P]
Bathe, K.-J. Finite Element Procedures, Prentice Hall (1996), §10.3.4 (cyclic symmetry). [textbook derivation matching the partitioned [u_I; u_L] variable convention used below]

class femorph_solver.solvers.cyclic.CyclicModalResult(harmonic_index: int, n_sectors: int, omega_sq: ndarray, frequency: ndarray, mode_shapes: ndarray, axis_point: ndarray | None = None, axis_dir: ndarray | None = None)[source]#

Bases: object

Modes for one harmonic index of a cyclic-symmetry modal solve.

axis_dir: ndarray | None = None#: (3,) — unit vector along the rotor axis. None means “not recorded”; writers default to +Z. Populated by solve_cyclic_modal() from the supplied pair_rotation.

axis_point: ndarray | None = None#: (3,) — point on the rotor rotation axis. None means “not recorded by the caller”; writers default to the origin in that case. Populated by solve_cyclic_modal() when a pair_rotation is supplied.

frequency: ndarray#: (n_modes,) — cyclic frequencies f [Hz].

mode_shapes: ndarray#: (N, n_modes) — complex mode shapes on the base sector, indexed by the original DOF layout (not the reduced layout). For k=0 and k=N/2 the imaginary part is identically zero.

n_sectors: int#: Number of sectors in the full rotor (same for every result in a sweep).

omega_sq: ndarray#: (n_modes,) — ω² = (2πf)²; numerically negative values clipped to 0.

femorph_solver.solvers.cyclic.save_cyclic_modal_results(results: Sequence[CyclicModalResult], path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise a cyclic-modal sweep to a single disk-backed .pv file.

Parameters:

results (sequence of CyclicModalResult) – Per-harmonic-index results as returned by solve_cyclic_modal(). All entries must share the same n_sectors.
path (str or pathlib.Path) – Destination file.
model (femorph_solver.model.Model) – Model whose K / M matrices this sweep was computed on; supplies the mesh and DOF layout.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

Notes

The DOF-indexed (N, n_modes) complex mode shapes on each per-k result get re-projected onto per-(k, mode), per-node (n_modes, n_points, 6) complex arrays — the layout write_cyclic_modal_result() expects.

femorph_solver.solvers.cyclic.solve_cyclic_modal(K: spmatrix | sparray, M: spmatrix | sparray, low_node_dofs: ndarray, high_node_dofs: ndarray, n_sectors: int, n_modes: int = 10, *, pair_rotation: ndarray | None = None, harmonic_indices: Sequence[int] | None = None, prescribed: dict[int, float] | None = None, thread_limit: int | None = None, linear_solver: str = 'auto') → list[CyclicModalResult][source]#

Solve the cyclic-symmetry modal problem on a single base sector.

Parameters:

K (scipy.sparse) – Global stiffness / mass of the base sector, before cyclic or Dirichlet reduction. Assumed real SPD.
M (scipy.sparse) – Global stiffness / mass of the base sector, before cyclic or Dirichlet reduction. Assumed real SPD.
low_node_dofs ((P, d) int arrays) – Global DOF indices at P paired boundary nodes, each with d DOFs (typically 3 for a 3D solid). low_node_dofs[i, :] and high_node_dofs[i, :] must address the same physical DOFs (same labels, same order) at the node pair that’s identified by the cyclic BC. A 1-D array is accepted for scalar (d=1) problems.
high_node_dofs ((P, d) int arrays) – Global DOF indices at P paired boundary nodes, each with d DOFs (typically 3 for a 3D solid). low_node_dofs[i, :] and high_node_dofs[i, :] must address the same physical DOFs (same labels, same order) at the node pair that’s identified by the cyclic BC. A 1-D array is accepted for scalar (d=1) problems.
n_sectors (int) – Number of repetitions that close the rotor (N).
n_modes (int, default 10) – Number of lowest modes per harmonic index.
pair_rotation ((d, d) array, optional) – Spatial rotation R(α) that maps sector j’s local frame to sector j+1’s, applied per node pair. Pass None when face DOFs are already in a per-sector local frame (then R=I implied).
harmonic_indices (sequence of int, optional) – Harmonic indices to solve. Defaults to 0, 1, …, N//2.
prescribed (mapping, optional) – Dirichlet BCs on the sector (global DOF index → 0). Currently only zero-value Dirichlet is supported; if a cyclic-face DOF is prescribed its partner must also be prescribed (so the cyclic constraint is trivially satisfied at that DOF).
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the per-harmonic dense Hermitian eigensolve. None defers to the global limit (see femorph_solver.set_thread_limit()).
linear_solver (str, default "auto") – Name of the registered linear backend used for the inner shift-invert factorisation in the sparse path. "auto" picks the fastest available SPD direct solver in the priority chain pardiso → cholmod → umfpack → superlu and emits a one-shot warning if pardiso is missing on a problem large enough to benefit from it. See femorph_solver.solvers.linear.list_linear_solvers().

Returns:

One CyclicModalResult per requested harmonic index, in the order requested.

Return type:

list[CyclicModalResult]

Linear transient solver (Newmark-β time integration).

Integrates M ü + C u̇ + K u = F(t) over a uniform timestep grid. Dirichlet constraints fix both the displacement and (implicitly) the velocity / acceleration at those DOFs to zero.

Newmark-β recurrence (standard form with parameters β, γ):

u_{n+1} = u_n + dt · u̇_n + dt² · [(1/2 − β) ü_n + β ü_{n+1}]
u̇_{n+1} = u̇_n + dt · [(1 − γ) ü_n + γ ü_{n+1}]

Rearranging M ü_{n+1} + C u̇_{n+1} + K u_{n+1} = F_{n+1} in terms of u_{n+1} gives the effective system:

A u_{n+1} = b_{n+1}
A = a0·M + a1·C + K
b = F_{n+1} + M(a0 u_n + a2 u̇_n + a3 ü_n) + C(a1 u_n + a4 u̇_n + a5 ü_n)

with the usual coefficient shorthands a0 = 1/(β dt²), a1 = γ/(β dt), a2 = 1/(β dt), a3 = 1/(2β) − 1, a4 = γ/β − 1, a5 = dt/2 · (γ/β − 2).

Defaults (β, γ) = (1/4, 1/2) are the unconditionally-stable average-acceleration method. This implementation is implicit and requires β > 0 — explicit (β=0) central-difference would need a different solver path that factors the mass matrix directly.

References

Newmark, N. M. “A method of computation for structural dynamics.” J. Eng. Mech. Div. (ASCE) 85 (EM3), 67-94 (1959). [original Newmark-β integration scheme]
Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (2000), §9.1-9.2. [stability + accuracy analysis of the (β, γ) family; conditions for second-order accuracy and unconditional stability]
Bathe, K.-J. Finite Element Procedures, 2nd ed., §9.2. [implementation form of the effective-stiffness recurrence used above and the reduction to a per-step linear solve]

class femorph_solver.solvers.transient.TransientResult(time: ndarray, displacement: ndarray, velocity: ndarray, acceleration: ndarray)[source]#

Bases: object

Result of a linear transient solve, sampled on the time grid.

time#

(n_steps + 1,) sample times [s].

Type:: numpy.ndarray

displacement#

(n_steps + 1, N) DOF-indexed displacement; row per step.

Type:: numpy.ndarray

velocity#

(n_steps + 1, N) DOF-indexed velocity.

Type:: numpy.ndarray

acceleration#

(n_steps + 1, N) DOF-indexed acceleration.

Type:: numpy.ndarray

save(path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise this transient result to a disk-backed .pv file.

Re-projects the DOF-indexed displacement / velocity / acceleration onto per-step, per-node (n_steps, n_points, n_dof_per_node) arrays via dof_map(), then hands off to write_transient_result().

Parameters:

path (str or pathlib.Path) – Destination file.
model (femorph_solver.model.Model) – Model whose K / M matrices this result was integrated on; supplies the mesh + DOF layout.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

femorph_solver.solvers.transient.solve_transient(K: csr_array, M: csr_array, F: ndarray | Callable[[float], ndarray], *, dt: float, n_steps: int, prescribed: dict[int, float] | None = None, C: csr_array | None = None, u0: ndarray | None = None, v0: ndarray | None = None, beta: float = 0.25, gamma: float = 0.5, thread_limit: int | None = None) → TransientResult[source]#

Integrate M ü + C u̇ + K u = F(t) with Newmark-β.

Parameters:

K – Global stiffness and mass matrices ((N, N) sparse).
M – Global stiffness and mass matrices ((N, N) sparse).
F – Load vector. Either a fixed (N,) array (step-constant load) or a callable t -> (N,).
dt – Uniform timestep.
n_steps – Number of steps to take. Returns n_steps + 1 samples including the initial state at t = 0.
prescribed – Dirichlet constraints {global_dof_index: prescribed_value}. Constrained DOFs are held at u_c with zero velocity and acceleration; their rows/cols are dropped from A.
C – Damping matrix. None means no damping.
u0 – Initial displacement and velocity. Default to zero.
v0 – Initial displacement and velocity. Default to zero.
beta – Newmark parameters. Defaults (1/4, 1/2) give the unconditionally-stable average-acceleration method. (0, 1/2) reduces to central-difference (explicit, conditionally stable).
gamma – Newmark parameters. Defaults (1/4, 1/2) give the unconditionally-stable average-acceleration method. (0, 1/2) reduces to central-difference (explicit, conditionally stable).
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the factorisation and per-step back-solve. None defers to the global limit (see femorph_solver.set_thread_limit()).

Linear solver backends#

Pluggable linear solvers.

femorph-solver provides several backends for solving A x = b with a sparse SPD or general sparse A. The default is SciPy’s built-in SuperLU; the rest are optional and require extra dependencies listed in the class’s install_hint.

Usage:

from femorph_solver.solvers.linear import get_linear_solver

Solver = get_linear_solver("cholmod")  # or "superlu", "umfpack", "pardiso", ...
s = Solver(K_csc)
x = s.solve(b)

list_linear_solvers() reports every registered backend with an available flag so UIs / benchmarks can enumerate them.

class femorph_solver.solvers.linear.CGSolver(A, *, assume_spd: bool = True, rtol: float = 1e-10, atol: float = 0.0, maxiter: int | None = None)[source]#

Preconditioned Conjugate Gradient for SPD A.

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'cg'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.CholmodSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

CHOLMOD supernodal sparse Cholesky (SPD only).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = 'install SuiteSparse system library + `pip install scikit-sparse`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'cholmod'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.GMRESSolver(A, *, assume_spd: bool = False, rtol: float = 1e-10, atol: float = 0.0, restart: int = 50, maxiter: int | None = None)[source]#

Preconditioned GMRES (general sparse).

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'gmres'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.LinearSolverBase(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

Bases: object

Protocol every linear-solver backend implements.

Concrete subclasses set the class attributes (name, kind, spd_only, install_hint), override available(), and implement _factor() and solve().

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = ''#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = ''#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.MUMPSSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False, ooc: bool = False, ordering: str = 'auto')[source]#

MUMPS sparse direct solver via python-mumps (multifrontal).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install python-mumps` (requires libmumps — e.g. apt install libmumps-scotch-dev, brew install mumps, or conda install -c conda-forge mumps)'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'mumps'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.PardisoSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

Intel MKL Pardiso via pypardiso (multi-threaded).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install pypardiso` (pulls in MKL)'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'pardiso'#: Short identifier used by get_linear_solver("...").

release() → None[source]#

Release MKL’s internal factor memory.

MKL Pardiso holds the numerical factor on a 64-slot pt handle array that survives Python-side garbage collection — dropping our PardisoSolver reference does not free it. pypardiso.free_memory(everything=True) runs phase -1, which tears down the MKL side. On a 200×200×2 HEX8 plate that’s ~390 MB of permanent factor that would otherwise live until process exit; for workflows that run multiple solves (modal + static + modal, notebooks iterating over meshes) the leak accumulates linearly.

Idempotent — safe to call multiple times. Safe to call mid-interpreter-shutdown because pypardiso.free_memory constructs a tiny dummy A / b locally and the MKL library handle is already resolved.

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.PyAMGSolver(A, *, assume_spd: bool = True, rtol: float = 1e-10, maxiter: int = 200)[source]#

Smoothed-aggregation AMG preconditioned CG (SPD).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install pyamg`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'pyamg'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

exception femorph_solver.solvers.linear.SolverUnavailableError[source]#

Bases: RuntimeError

Raised when a registered solver’s optional dependency is missing.

class femorph_solver.solvers.linear.SuperLUSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

SciPy’s bundled SuperLU (always available).

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'superlu'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

class femorph_solver.solvers.linear.UMFPACKSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

UMFPACK LU factorization via scikit-umfpack.

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = 'install SuiteSparse system library + `pip install scikit-umfpack`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'umfpack'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

femorph_solver.solvers.linear.get_linear_solver(name: str, *, n_dof: int | None = None) → type[LinearSolverBase][source]#

Look up a registered linear-solver class by name.

Passing "auto" resolves via select_default_linear_solver(); pass n_dof so the auto picker can prefer Pardiso on large problems where its parallel factorisation dominates, and skip it on small ones where its fixed overhead erases the gain.

Raises KeyError if the name is unknown, or SolverUnavailableError if the backend is registered but its optional dependency is not installed.

femorph_solver.solvers.linear.list_linear_solvers() → list[dict[str, object]][source]#: Return a list of {name, available, kind, spd_only, install_hint} rows.

femorph_solver.solvers.linear.select_default_linear_solver(*, spd: bool = True, n_dof: int | None = None) → str[source]#

Return the best available backend name for the requested workload.

Size-aware: when n_dof >= 50_000 the preference is pardiso → cholmod → umfpack → superlu because Pardiso’s threaded factor dominates. Below the threshold (or when n_dof is not known) Pardiso is demoted because its fixed setup cost erases its factor-time lead on small systems. superlu is always available as the final fallback.

If Pardiso would have been chosen by size but is not installed, a one-shot UserWarning points at the install hint.

Common base class for pluggable linear solvers.

class femorph_solver.solvers.linear._base.LinearSolverBase(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

Bases: object

Protocol every linear-solver backend implements.

Concrete subclasses set the class attributes (name, kind, spd_only, install_hint), override available(), and implement _factor() and solve().

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = ''#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = ''#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

exception femorph_solver.solvers.linear._base.SolverUnavailableError[source]#

Bases: RuntimeError

Raised when a registered solver’s optional dependency is missing.

Intel MKL Pardiso via pypardiso.

Pardiso is multi-threaded and wins big on problems ≳30k DOF. Requires Intel MKL, which pypardiso pulls in as a dependency.

The import pypardiso call eagerly loads MKL (~200 ms on first touch), so this module defers the import until a PardisoSolver is actually constructed. That keeps import femorph_solver snappy for users who don’t touch Pardiso, and avoids MKL init penalising the CHOLMOD/SuperLU paths on small models where the auto-selector never reaches Pardiso.

When the caller asserts assume_spd=True we switch to MKL’s real-SPD mode (mtype=2) which runs a sparse Cholesky factor. On a 57 k-DOF ptr.cdb K_ff: factor 0.39 s / solve 43 ms vs the default unsymmetric-LU path’s 0.51 s / 94 ms — roughly 2× faster solves with identical output (max abs diff 5 × 10⁻¹⁶).

References

Schenk, O. & Gärtner, K. “Solving unsymmetric sparse systems of linear equations with PARDISO.” Future Generation Computer Systems 20 (3), 475-487 (2004). [the PARDISO algorithm Intel MKL ships]
Schenk, O. & Gärtner, K. “On fast factorization pivoting methods for sparse symmetric indefinite systems.” ETNA 23, 158-179 (2006). [supernodal + Bunch-Kaufman pivoting paths used for indefinite systems (mtype=-2)]
Intel. Intel oneAPI Math Kernel Library Developer Reference, “Intel oneMKL PARDISO — Parallel Direct Sparse Solver Interface.” [iparm flag semantics the DirectMklPardisoSolver sibling wrapper drives directly via ctypes]

class femorph_solver.solvers.linear._pardiso.PardisoSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

Intel MKL Pardiso via pypardiso (multi-threaded).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install pypardiso` (pulls in MKL)'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'pardiso'#: Short identifier used by get_linear_solver("...").

release() → None[source]#

Release MKL’s internal factor memory.

MKL Pardiso holds the numerical factor on a 64-slot pt handle array that survives Python-side garbage collection — dropping our PardisoSolver reference does not free it. pypardiso.free_memory(everything=True) runs phase -1, which tears down the MKL side. On a 200×200×2 HEX8 plate that’s ~390 MB of permanent factor that would otherwise live until process exit; for workflows that run multiple solves (modal + static + modal, notebooks iterating over meshes) the leak accumulates linearly.

Idempotent — safe to call multiple times. Safe to call mid-interpreter-shutdown because pypardiso.free_memory constructs a tiny dummy A / b locally and the MKL library handle is already resolved.

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

Direct ctypes wrapper around MKL Pardiso.

Why not keep using pypardiso?

iparm(28)=1 silently no-ops on this MKL build. pypardiso sets the mixed-precision flag but after factorize the solver reports iparm[27]=0 — MKL chose to keep a double factor anyway. We want a direct path so we can verify (and, where the build allows, actually engage) single-precision factorisation that halves factor memory.

Per-solve ``A.indptr.astype(int32) + 1`` + ``A.indices.astype + 1`` churn. PR #112 patched this with a cached _call_pardiso override; doing it in-house from the start avoids the monkey-patch entirely and lets us tell MKL we ship 0-based indices (iparm[34]=1) so no +1 copy ever happens.

Cleanup. pypardiso does not call phase=-1 on GC, so the MKL factor leaks across repeated solves (see agent-1’s #107). We own __del__ here and release explicitly.

Phase control. pypardiso bundles analyse + factorise into factorize(); for user code that wants to re-solve with the same pattern but updated values (modal sweeps over perturbed matrices), we can split phase=11 (analyse) from phase=22 (factor only) and amortise the analyse cost.

The MKL library load path piggybacks on pypardiso’s discovery — we defer to its already-loaded handle when available so we don’t duplicate the mkl_rt search heuristic. Users who set the PYPARDISO_MKL_RT environment variable get the same resolution on this backend.

References

Schenk, O. & Gärtner, K. “Solving unsymmetric sparse systems of linear equations with PARDISO.” Future Generation Computer Systems 20 (3), 475-487 (2004). [the Pardiso algorithm implemented by MKL]
Schenk, O. & Gärtner, K. “On fast factorization pivoting methods for sparse symmetric indefinite systems.” ETNA 23, 158-179 (2006).
Intel. Intel oneAPI Math Kernel Library Developer Reference, “Intel oneMKL PARDISO — Parallel Direct Sparse Solver Interface.” [authoritative iparm / phase / mtype spec — the set of flags this wrapper drives directly]
Intel. Intel MKL PARDISO out-of-core user guide. [iparm(60) OOC mode + MKL_PARDISO_OOC_MAX_CORE_SIZE heuristic exercised by the ooc=True path]

class femorph_solver.solvers.linear._mkl_pardiso.DirectMklPardisoSolver(A, *, assume_spd: bool = False, mixed_precision: bool | None = None, ooc: bool = False, msglvl: int = 0)[source]#

ctypes-level MKL Pardiso direct solver.

Parameters:

A (scipy.sparse.csr_array) – Real-valued SPD (assume_spd=True) or general square sparse matrix. When assume_spd=True and A.femorph_triu is set, we use the supplied triu storage directly. Otherwise we extract sp.triu(A) ourselves before handing to MKL.
assume_spd (bool, default False) – Flip to MKL’s mtype=2 real-SPD Cholesky path.
mixed_precision (bool or None, default None) – True asks MKL for a single-precision factor with double- precision iterative refinement (iparm(28)=1). None (auto) stays double-precision — the MKL builds we test against silently decline this flag, and even when honoured the iterative-refinement buffers often undo the factor-memory win. False forces double precision regardless.
ooc (bool, default False) – True enables out-of-core factorisation (iparm(60)=2 — spill factor to disk). Trades wall-time for the ability to factor K that would otherwise exceed RAM.
msglvl (int, default 0) – Pass 1 to dump MKL’s factor statistics to stdout — useful when diagnosing whether a flag (MP / OOC / parallel) actually engaged.

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = 'same MKL as ``pardiso`` (``pip install pypardiso``)'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'mkl_direct'#: Short identifier used by get_linear_solver("...").

release() → None[source]#: Release MKL’s internal factor memory (phase = -1).

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

CHOLMOD sparse Cholesky via scikit-sparse.

For SPD A (stiffness matrices after BCs, for example), CHOLMOD’s supernodal Cholesky is ~2-3× faster than SuperLU’s LU and uses about half the memory.

References

Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. “Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate.” ACM Trans. Math. Softw. 35 (3), 22:1-22:14 (2008). [the shipped supernodal algorithm + its update paths]
Davis, T. A. & Hager, W. W. “Dynamic supernodes in sparse Cholesky update/downdate and triangular solves.” ACM Trans. Math. Softw. 35 (4), 27:1-27:23 (2009). [the row-and-column modification path used when updating factorisations across parametric sweeps]

class femorph_solver.solvers.linear._cholmod.CholmodSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

CHOLMOD supernodal sparse Cholesky (SPD only).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = 'install SuiteSparse system library + `pip install scikit-sparse`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'cholmod'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

MUMPS multifrontal sparse direct solver via python-mumps.

MUMPS is the go-to direct solver when Pardiso is unavailable or when the redistributable-license / out-of-core story matters:

Licence — CeCILL-C (LGPL-like). Unlike MKL Pardiso, MUMPS can be shipped with a commercial redistribution without the MKL EULA.
Out-of-core — mature OOC mode via factor(ooc=True) / MUMPS’s ICNTL(22); spills the factor to disk so RAM stops being the ceiling on single-node solves. Pardiso has an OOC mode too but it is less battle-tested.
Distributed memory — MUMPS is the on-ramp to MPI-distributed solves across nodes, which neither Pardiso nor SuperLU offer.

On a single shared-memory node at ptr.cdb scale (≤1 M DOFs, SPD) MUMPS is usually 1.5–3× slower than Pardiso and uses ~10–30 % more factor memory. Keep it available, but auto still prefers Pardiso > CHOLMOD > MUMPS on SPD workloads.

Implementation notes for the python-mumps binding (pip install python-mumps, import mumps):

The Context constructor takes no sym= argument; symmetry is declared on Context.set_matrix() as symmetric=True, which internally maps to MUMPS’s sym=2 (general-symmetric LDLᵀ). The wrapper does not expose MUMPS’s sym=1 “strictly SPD” mode because python-mumps 0.0.6 doesn’t — sym=2 handles SPD input correctly, just without the small memory win of skipping pivoting.
The Context docstring warns “only complex numbers supported”, but the underlying code dispatches on the input dtype (dmumps for float64, zmumps for complex128). Float64 works fine.

References

Amestoy, P. R., Duff, I. S., L’Excellent, J.-Y., & Koster, J. “A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling.” SIAM J. Matrix Anal. Appl. 23 (1), 15-41 (2001). [the MUMPS multifrontal + distributed-scheduling algorithm]
Amestoy, P. R., Guermouche, A., L’Excellent, J.-Y., & Pralet, S. “Hybrid scheduling for the parallel solution of linear systems.” Parallel Computing 32 (2), 136-156 (2006). [hybrid static/dynamic scheduler shipped in MUMPS 4.x+]
Agullo, E., Guermouche, A., & L’Excellent, J.-Y. “Reducing the I/O volume in sparse out-of-core multifrontal methods.” SIAM J. Scientific Computing 31 (6), 4774-4794 (2010). [the OOC path exposed here via ICNTL(22)]

class femorph_solver.solvers.linear._mumps.MUMPSSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False, ooc: bool = False, ordering: str = 'auto')[source]#

MUMPS sparse direct solver via python-mumps (multifrontal).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install python-mumps` (requires libmumps — e.g. apt install libmumps-scotch-dev, brew install mumps, or conda install -c conda-forge mumps)'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'mumps'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

UMFPACK sparse LU via scikit-umfpack.

UMFPACK tends to be faster than SuperLU on unsymmetric sparse problems with moderate fill-in. Ships as a separate pip package because it links against the SuiteSparse system library.

References

Davis, T. A. “Algorithm 832: UMFPACK V4.3 — An Unsymmetric-Pattern Multifrontal Method.” ACM Trans. Math. Softw. 30 (2), 196-199 (2004). [the shipped multifrontal algorithm]
Davis, T. A. & Duff, I. S. “An Unsymmetric-Pattern Multifrontal Method for Sparse LU Factorization.” SIAM J. Matrix Anal. Appl. 18 (1), 140-158 (1997). [foundational method paper]

class femorph_solver.solvers.linear._umfpack.UMFPACKSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

UMFPACK LU factorization via scikit-umfpack.

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = 'install SuiteSparse system library + `pip install scikit-umfpack`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'umfpack'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

SuperLU direct sparse solver via scipy.sparse.linalg.splu.

Default backend — always available because SciPy ships SuperLU. Handles general (nonsymmetric) sparse A; uses the Sherman-family permutation MMD_AT_PLUS_A by default — the standard heuristic for SPD problems.

References

Demmel, J. W., Eisenstat, S. C., Gilbert, J. R., Li, X. S., & Liu, J. W. H. “A Supernodal Approach to Sparse Partial Pivoting.” SIAM J. Matrix Anal. Appl. 20 (3), 720-755 (1999). [right-looking supernodal LU — the core SuperLU algorithm]
Li, X. S. “An Overview of SuperLU: Algorithms, Implementation, and User Interface.” ACM Trans. Math. Softw. 31 (3), 302-325 (2005). [shipped implementation wrapped by SciPy]

class femorph_solver.solvers.linear._superlu.SuperLUSolver(A: spmatrix | sparray, *, assume_spd: bool = False, mixed_precision: bool | None = False)[source]#

SciPy’s bundled SuperLU (always available).

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'direct'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'superlu'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

Conjugate Gradient iterative solver (SPD).

Uses scipy.sparse.linalg.cg() with a diagonal (Jacobi) preconditioner. Always available; best for very large sparse SPD systems where factorization memory would be prohibitive.

References

Hestenes, M. R. & Stiefel, E. “Methods of Conjugate Gradients for Solving Linear Systems.” J. Res. Natl. Bur. Stand. 49 (6), 409-436 (1952). [the original CG algorithm]
Saad, Y. Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM (2003), §6. [modern treatment + convergence theory]

class femorph_solver.solvers.linear._cg.CGSolver(A, *, assume_spd: bool = True, rtol: float = 1e-10, atol: float = 0.0, maxiter: int | None = None)[source]#

Preconditioned Conjugate Gradient for SPD A.

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'cg'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

GMRES iterative solver (general).

Nonsymmetric Krylov solver — use when A is not SPD. Jacobi preconditioner by default.

References

Saad, Y. & Schultz, M. H. “GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.” SIAM J. Sci. Stat. Comput. 7 (3), 856-869 (1986). [the original GMRES algorithm implemented by scipy.sparse.linalg.gmres]
Saad, Y. Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM (2003), §6.5. [modern treatment including restart + the restart-length trade-off scipy exposes via its restart kwarg]

class femorph_solver.solvers.linear._gmres.GMRESSolver(A, *, assume_spd: bool = False, rtol: float = 1e-10, atol: float = 0.0, restart: int = 50, maxiter: int | None = None)[source]#

Preconditioned GMRES (general sparse).

install_hint: ClassVar[str] = 'bundled with SciPy — no install needed'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'gmres'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = False#: True if the solver requires a symmetric positive-definite matrix.

Algebraic multigrid (AMG) preconditioned CG via pyamg.

AMG is the standard industrial preconditioner for elliptic problems like structural stiffness — often 10-100× fewer iterations than Jacobi-CG on large unstructured SPD systems.

References

Ruge, J. W. & Stüben, K. “Algebraic multigrid.” In Multigrid Methods (S. F. McCormick, ed.), SIAM Frontiers in Applied Mathematics, Chapter 4, pp. 73-130 (1987). [the foundational classical-AMG algorithm]
Vaněk, P., Mandel, J., & Brezina, M. “Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems.” Computing 56 (3), 179-196 (1996). [smoothed-aggregation variant, pyamg’s default smoothed_aggregation_solver]
Bell, N., Olson, L. N., & Schroder, J. “PyAMG: Algebraic Multigrid Solvers in Python v5.0.” Journal of Open Source Software 8 (81), 5495 (2023). [the shipped implementation]

class femorph_solver.solvers.linear._pyamg.PyAMGSolver(A, *, assume_spd: bool = True, rtol: float = 1e-10, maxiter: int = 200)[source]#

Smoothed-aggregation AMG preconditioned CG (SPD).

static available() → bool[source]#: Return True if the solver can be constructed on this install.

install_hint: ClassVar[str] = '`pip install pyamg`'#: User-facing hint printed when the backend is unavailable.

kind: ClassVar[str] = 'iterative'#: "direct" (factor + forward/back solve) or "iterative".

name: ClassVar[str] = 'pyamg'#: Short identifier used by get_linear_solver("...").

solve(b: ndarray) → ndarray[source]#: Solve A x = b for a single RHS.

spd_only: ClassVar[bool] = True#: True if the solver requires a symmetric positive-definite matrix.

Eigen solver backends#

Pluggable generalized-eigenvalue solvers.

All backends solve K φ = λ M φ for the n_modes smallest-λ eigenpairs. The default backend ("arpack") uses scipy.sparse.linalg.eigsh() in shift-invert mode with σ = 0 — identical to femorph-solver’s behaviour up to this release. Alternative backends let callers trade accuracy / convergence / memory against speed and dependency footprint.

Usage:

from femorph_solver.solvers.eigen import get_eigen_solver

Solver = get_eigen_solver("arpack")
w, v = Solver.solve(K, M, n_modes=10, linear_solver="cholmod")

class femorph_solver.solvers.eigen.ArpackSolver[source]#

scipy.sparse.linalg.eigsh with shift-invert, pluggable inner solver.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0, mixed_precision: bool | None = False, v0: ndarray | None = None) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

class femorph_solver.solvers.eigen.EigenSolverBase[source]#

Bases: object

Every eigen solver implements solve(K, M, n_modes, ...).

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

class femorph_solver.solvers.eigen.LapackDenseSolver[source]#

Dense LAPACK eigh — correct at any size, but O(N^3).

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

class femorph_solver.solvers.eigen.LobpcgSolver[source]#

scipy.sparse.linalg.lobpcg with pluggable preconditioner.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 1e-08, maxiter: int = 400, preconditioner: str | LinearOperator = 'factor') → tuple[ndarray, ndarray][source]#

Solve the generalised eigenproblem K φ = ω² M φ with LOBPCG.

Parameters:: preconditioner ("factor" (default), "jacobi", "none",) – or an explicit scipy.sparse.linalg.LinearOperator. "factor" reproduces the historical behaviour (full direct factor of K). "jacobi" is the memory-friendly path — diagonal inverse of K, no factor, O(n) storage. "none" disables preconditioning. When an explicit LinearOperator is supplied linear_solver is ignored.

class femorph_solver.solvers.eigen.PrimmeSolver[source]#

primme.eigsh for generalized SPD eigenproblems.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 1e-10) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

femorph_solver.solvers.eigen.get_eigen_solver(name: str) → type[EigenSolverBase][source]#: Look up a registered eigen-solver class by name.

femorph_solver.solvers.eigen.list_eigen_solvers() → list[dict[str, object]][source]#: Return a list of {name, available, kind, spd_only, install_hint} rows.

Common base class for eigen-solver backends.

class femorph_solver.solvers.eigen._base.EigenSolverBase[source]#

Bases: object

Every eigen solver implements solve(K, M, n_modes, ...).

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

ARPACK Implicitly Restarted Lanczos via scipy.sparse.linalg.eigsh.

Shift-invert mode at σ = 0 finds the n_modes eigenvalues closest to zero — exactly the lowest natural frequencies we want for modal analysis. The inner (K - σM)^-1 apply goes through the pluggable femorph-solver linear-solver registry so callers can pick SuperLU (default), CHOLMOD, Pardiso, or anything else registered.

K·v and M·v inside ARPACK’s generalized-eigenvalue iteration are wrapped in a custom OpenMP parallel CSR matvec (_core.csr_matvec) when the matrix is sufficiently large — scipy’s built-in csr_matvec is single-threaded and costs ~13 % of wall time on a 600 k-DOF ptr-sized system.

References

Lehoucq, R. B. & Sorensen, D. C. “Deflation Techniques for an Implicitly Restarted Arnoldi Iteration.” SIAM J. Matrix Anal. Appl. 17 (4), 789-821 (1996). [implicit-restart Arnoldi algorithm at the core of ARPACK’s eigsh]
Lehoucq, R. B., Sorensen, D. C., & Yang, C. ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM (1998). [shipped reference for the ARPACK library scipy wraps, including the shift-invert mode used here]
Ericsson, T. & Ruhe, A. “The Spectral Transformation Lanczos Method for the Numerical Solution of Large Sparse Generalized Symmetric Eigenvalue Problems.” Math. Comp. 35 (152), 1251-1268 (1980). [spectral-transformation theory justifying the (K - σM)^-1 shift-invert applied on top of Lanczos]

class femorph_solver.solvers.eigen._arpack.ArpackSolver[source]#

scipy.sparse.linalg.eigsh with shift-invert, pluggable inner solver.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0, mixed_precision: bool | None = False, v0: ndarray | None = None) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

LOBPCG (Locally Optimal Block Preconditioned Conjugate Gradient).

Iterative eigensolver — works well when an efficient preconditioner is available (e.g. AMG). Without one it generally loses to shift-invert Lanczos on structural problems; the preconditioner is the whole game.

This wrapper offers three preconditioner modes:

"factor" (default) — run a full direct solver on K and use the factor as a preconditioner. Convergent but has the same memory footprint as the shift-invert path; useful when the problem is too ill-conditioned for a cheaper preconditioner but you need LOBPCG’s convergence pattern on clustered eigenvalues.
"jacobi" — diagonal Jacobi preconditioner (M = diag(K)). Free. Memory-friendly: storage is O(n) floats, no factor. On well-conditioned structural problems it typically converges in 2–5× more iterations than "factor" but uses ~10× less memory at the factor stage. The right default when K is too big to factor.
"none" — no preconditioner. Occasionally wins on tiny problems where the preconditioner setup dominates; mostly for diagnostics.

Pass an explicit scipy.sparse.linalg.LinearOperator via the preconditioner= argument to override entirely.

References

Knyazev, A. V. “Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method.” SIAM J. Sci. Comput. 23 (2), 517-541 (2001). [the LOBPCG algorithm implemented by scipy.sparse.linalg.lobpcg]
Knyazev, A. V., Argentati, M. E., Lashuk, I., & Ovtchinnikov, E. E. “Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc.” SIAM J. Sci. Comput. 29 (5), 2224-2239 (2007). [block-variant stability fixes carried forward into scipy’s implementation]

class femorph_solver.solvers.eigen._lobpcg.LobpcgSolver[source]#

scipy.sparse.linalg.lobpcg with pluggable preconditioner.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 1e-08, maxiter: int = 400, preconditioner: str | LinearOperator = 'factor') → tuple[ndarray, ndarray][source]#

Solve the generalised eigenproblem K φ = ω² M φ with LOBPCG.

Parameters:: preconditioner ("factor" (default), "jacobi", "none",) – or an explicit scipy.sparse.linalg.LinearOperator. "factor" reproduces the historical behaviour (full direct factor of K). "jacobi" is the memory-friendly path — diagonal inverse of K, no factor, O(n) storage. "none" disables preconditioning. When an explicit LinearOperator is supplied linear_solver is ignored.

Dense LAPACK fallback via scipy.linalg.eigh.

For very small problems (N ≲ 500) a dense generalized eigenproblem is faster than any sparse solver and numerically bullet-proof. Serves as our reference path in tests.

References

Anderson, E. et al. LAPACK Users’ Guide, 3rd ed., SIAM (1999), Chapter 2. [generalized symmetric-definite eigenproblem — the DSYGVD and DSYGVX drivers scipy.linalg.eigh calls]
Parlett, B. N. The Symmetric Eigenvalue Problem, SIAM Classics in Applied Mathematics (1998). [the bisection + Rayleigh-quotient + divide-and-conquer algorithms underlying LAPACK’s symmetric routines]
Golub, G. H. & Van Loan, C. F. Matrix Computations, 4th ed., JHU Press (2013), §8. [reference treatment of the symmetric eigenproblem used for the dense path]

class femorph_solver.solvers.eigen._lapack_dense.LapackDenseSolver[source]#

Dense LAPACK eigh — correct at any size, but O(N^3).

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 0.0) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

PRIMME via primme.eigsh — a state-of-the-art block-Davidson solver.

Often beats ARPACK on large SPD generalized eigenproblems because of its restart + soft-locking scheme. Pip-installable and self-contained (no extra system libs).

References

Stathopoulos, A. & McCombs, J. R. “PRIMME: PReconditioned Iterative MultiMethod Eigensolver — Methods and software description.” ACM Trans. Math. Softw. 37 (2), 21:1-21:30 (2010). [the shipped PRIMME library: block-JD + GD+k restart + the soft-locking convergence machinery]
Stathopoulos, A. “Nearly Optimal Preconditioned Methods for Hermitian Eigenproblems Under Limited Memory. Part I: Seeking One Eigenvalue.” SIAM J. Sci. Comput. 29 (2), 481-514 (2007). [the GD+k restart scheme inside PRIMME]
Davidson, E. R. “The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices.” J. Comp. Phys. 17 (1), 87-94 (1975). [the original Davidson method PRIMME descends from]

class femorph_solver.solvers.eigen._primme.PrimmeSolver[source]#

primme.eigsh for generalized SPD eigenproblems.

static solve(K: spmatrix | sparray, M: spmatrix | sparray, n_modes: int, *, sigma: float = 0.0, linear_solver: str = 'auto', tol: float = 1e-10) → tuple[ndarray, ndarray][source]#: Return (eigenvalues, eigenvectors) sorted ascending by λ.

Elements#

femorph-solver element library.

Each element class implements ElementBase and registers itself via register(). The registry is how the assembler (and anything else walking a mesh) dispatches per-element math from the per-cell element-type identifier on a model’s grid.

Kernel modules carry topology-first names (hex8, hex20, tet10, wedge15 / pyr13 from _wedge15_pyr13, beam2, truss2, quad4_shell, quad4_plane, spring, point_mass). Foreign-deck nomenclature (e.g. MAPDL SOLID185 / NASTRAN CHEXA / Abaqus C3D8) is translated to neutral kernel names at each interop boundary and never reaches the public API of this package.

class femorph_solver.elements.ElementBase[source]#

Bases: ABC

abstractmethod static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

abstractmethod static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

femorph_solver.elements.get(name: str) → type[ElementBase][source]#

Return the kernel class registered under its neutral name.

Only canonical kernel names ("HEX8", "TET10", …) are accepted. Foreign-deck names must be translated at the interop boundary (e.g. femorph_solver.interop.mapdl.APDL) before being passed here.

femorph_solver.elements.is_element_spec(obj: object) → bool[source]#

Return True if obj is an element-spec class or instance.

Used by Model.assign() to detect the spec form vs. the legacy string / ElementType form.

Abstract element interface.

class femorph_solver.elements._base.ElementBase[source]#

Bases: ABC

abstractmethod static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

abstractmethod static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Hex8 — 8-node linear hexahedral solid.

Tri-linear shape functions on the reference cube ξ, η, ζ ∈ [-1, 1] with 2×2×2 Gauss integration. For each node i = 0…7 with natural coordinates (ξᵢ, ηᵢ, ζᵢ) ∈ {±1}³:

Nᵢ(ξ, η, ζ) = (1/8)(1 + ξᵢ ξ)(1 + ηᵢ η)(1 + ζᵢ ζ)

Node ordering matches VTK_HEXAHEDRON — the eight cube corners traversed (-,-,-), (+,-,-), (+,+,-), (-,+,-), (-,-,+), (+,-,+), (+,+,+), (-,+,+).

Stiffness / mass#

K_e = ∫_V Bᵀ C B dV ≈ Σ_g Bᵀ(ξ_g) C B(ξ_g) det J(ξ_g) M_e = ρ ∫_V Nᵀ N dV ≈ Σ_g ρ Nᵀ(ξ_g) N(ξ_g) det J(ξ_g)

with isotropic linear elasticity in Voigt form (engineering shears):

C = λ·[111,000] ⊕ 2G·[diag(1,1,1,½,½,½)]

where λ = Eν/((1+ν)(1-2ν)) and G = E/(2(1+ν)).

Real constants#

None required — geometry comes from node coordinates.

Material#

EX, PRXY : elastic moduli (both mandatory) DENS : density (required for M_e)

Lumped mass#

Row-sum (HRZ) lumping distributes each row of the consistent mass onto its diagonal. For 8-node hex this gives m_total / 8 per corner node — the standard and physically correct lumped form.

Integration flag#

material["_HEX8_INTEGRATION"] selects the stiffness formulation:

"full" (default) — 2×2×2 Gauss full integration with Hughes volumetric B-bar. The dilatation part of B at each Gauss point is replaced by its element-volume-weighted average so near-incompressible regimes do not lock.
"plain_gauss" — plain 2×2×2 Gauss without the volumetric correction. Useful for cross-checks against codes that don’t ship B-bar (e.g. the scikit-fem cross-FEM verification fixture); not recommended for production runs because it locks under near-incompressibility.
"enhanced_strain" — Simo-Rifai 9-parameter enhanced assumed strain (EAS) with Pian-Tong Jacobian scaling. Suppresses shear locking in bending-dominated regular hex meshes.

References

Shape functions — trilinear Lagrange on the 8-node hexahedron: Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, Ch. 6. Equivalent form in Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §6.2.
2×2×2 Gauss quadrature rule (product rule, 3 integration points → 8 product points): Zienkiewicz & Taylor §5.10 Table 5.3.
B-bar formulation ("full", default — volumetric dilatation replaced by its element-volume-weighted mean to avoid near-incompressible locking): Hughes, T.J.R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, 2000 (reprint of Prentice-Hall 1987), §4.4 (“B-bar method”). Original formulation in Hughes, T.J.R., “Generalization of selective integration procedures to anisotropic and nonlinear media,” IJNME 15 (1980), pp. 1413–1418.
Enhanced assumed strain ("enhanced_strain" — 9 internal parameters, Pian-Tong Jacobian scaling): Simo, J.C. and Rifai, M.S., “A class of mixed assumed strain methods and the method of incompatible modes,” IJNME 29 (1990), pp. 1595–1638.
Row-sum / HRZ lumped mass: Hinton, E., Rock, T. and Zienkiewicz, O.C., “A note on mass lumping and related processes in the finite element method,” Earthquake Engng. Struct. Dyn. 4 (1976), pp. 245–249. Zienkiewicz & Taylor §16.2.2.

class femorph_solver.elements.hex8.Hex8[source]#

Bases: ElementBase

static eel_batch(coords: ndarray, u_e: ndarray, material: dict[str, float] | None = None) → ndarray | None[source]#

Per-element elastic strain at every element node.

coords is (n_elem, 8, 3) and u_e is (n_elem, 24) — the element DOF vector laid out as [ux0, uy0, uz0, ux1, …]. Returns (n_elem, 8, 6) with Voigt components [εxx, εyy, εzz, γxy, γyz, γxz] (engineering shears), or None if the C++ extension is unavailable.

material is accepted for signature uniformity with plane kernels but is unused by the 3D-solid path (full 6-component Voigt strain is recovered directly from u_e).

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static ke_batch(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray | None[source]#

Vectorized Ke over (n_elem, 8, 3) coords.

Routes to the C++ B-bar kernel on the default path, to plain Gauss when the plain_gauss flag is set, and returns None for EAS (Python-only for now) so the assembler falls back to per-element ke() dispatch.

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Hex20 — 20-node quadratic hexahedral solid (serendipity).

Node ordering matches VTK_QUADRATIC_HEXAHEDRON: the 8 corner nodes come first in Hex8 order, then the 4 bottom-face mid-edge nodes, the 4 top-face mid-edge nodes, and finally the 4 vertical mid-edge nodes connecting the bottom and top faces.

Shape functions on the reference cube ξ, η, ζ ∈ [-1, 1]:

Corner node with (ξᵢ, ηᵢ, ζᵢ) ∈ {±1}³:

Nᵢ = (1/8)(1 + ξᵢ ξ)(1 + ηᵢ η)(1 + ζᵢ ζ)
     · (ξᵢ ξ + ηᵢ η + ζᵢ ζ - 2)

Mid-edge node with one natural coordinate zero (e.g. ξᵢ = 0):
```
Nᵢ = (1/4)(1 - ξ²)(1 + ηᵢ η)(1 + ζᵢ ζ)
```

Integration rule#

material["_HEX20_INTEGRATION"] selects the stiffness integration rule:

"reduced" (default) — 2×2×2 Gauss-Legendre (8 points). Single- element Ke has 6 extra zero-energy hourglass modes; on a connected mesh those modes almost always die under the assembly because a hex whose midside is shared with a neighbour has no free hourglass direction to excite. Make sure your mesh has no dangling midside nodes before using the default; if in doubt switch to "full".
"full" — 3×3×3 Gauss-Legendre (27 points). Ke is fully rank (54 non-zero eigenvalues + 6 rigid-body zeros) on a single hex.

Mass integration uses the 14-point Irons rule by default; pass material["_HEX20_MASS"] = "consistent" for the textbook 3×3×3 consistent rule. The wedge / pyramid degenerate forms also use the Irons 14-point rule (see _wedge15_pyr13.py); the rule is documented in the References section below.

Real constants#

None required — geometry comes from node coordinates.

Material#

EX, PRXY : elastic moduli (both mandatory) DENS : density (required for M_e)

References

Shape functions — 20-node serendipity hexahedron (corner nodes weighted by (ξᵢξ + ηᵢη + ζᵢζ − 2) correction, mid-edge nodes on zero-coordinate axes): Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §8.7.3. Original derivation in Ergatoudis, I., Irons, B.M. and Zienkiewicz, O.C., “Curved, isoparametric, ‘quadrilateral’ elements for finite element analysis,” Int. J. Solids Struct. 4 (1968), pp. 31–42, extended to 20-node hex in Zienkiewicz, O.C., The Finite Element Method in Engineering Science, 2nd ed., McGraw-Hill, 1971.
2×2×2 reduced Gauss integration ("reduced", default): Bedrosian, G., “Shape functions and integration formulas for three-dimensional finite element analysis,” IJNME 35 (1992), pp. 95–108 — the canonical reference for the uniform-reduced rule on the 20-node hex and its degenerate wedge / pyramid forms. Hourglass-mode analysis of the reduced rule on a single hex: Belytschko, T., Liu, W.K., Moran, B. and Elkhodary, K., Nonlinear Finite Elements for Continua and Structures, 2nd ed., Wiley, 2014, §8.6.
3×3×3 full Gauss-Legendre rule ("full", 27 points): Zienkiewicz & Taylor §5.10 Table 5.3.
14-point Irons integration rule (used for the consistent mass and in the degenerate wedge / pyramid forms): Irons, B.M., “Quadrature rules for brick based finite elements,” IJNME 3 (1971), pp. 293–294. Originally derived for the 20-node hex with degenerate wedge / pyramid as special cases, integrating polynomials through degree 5 exactly.
Voigt elastic matrix (isotropic linear elasticity, engineering shear convention): Cook, Malkus, Plesha, Witt Concepts and Applications of FEA, 4th ed., Wiley, 2002, §3.7.

class femorph_solver.elements.hex20.Hex20[source]#

Bases: ElementBase

static eel_batch(coords: ndarray, u_e: ndarray, material: dict[str, float] | None = None) → ndarray | None[source]#

Per-element elastic strain at every element node.

coords: (n_elem, 20, 3); u_e: (n_elem, 60). Returns (n_elem, 20, 6) — engineering-shear Voigt strain. material is accepted for signature uniformity with plane kernels but is unused.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Degenerate Hex20 shapes — 15-node wedge and 13-node pyramid.

These are the wedge and pyramid topologies expressed as VTK QUADRATIC_WEDGE (15 pts) and QUADRATIC_PYRAMID (13 pts). Foreign-deck readers (CDB, BDF, INP) translate their respective catalogue spellings (collapsed-corner Hex20, NASTRAN CPENTA15 / CPYRAM, Abaqus C3D15 / C3D13) at the boundary.

Two implementations live here:

Wedge15 — dedicated 15-node quadratic serendipity wedge kernel with its own shape functions and Gauss rules (9-pt for stiffness, 21-pt for consistent mass). Numerically equivalent to the 20-node-hex collapse on the meshes we have tested (the wedge formulation is mathematically the same hex with K=L and O=P collapsed — the dedicated kernel is just the compact path).
Pyr13 — currently a hex-fold wrapper K = Tᵀ K_hex T pending a dedicated 13-node pyramid implementation. Pyramids are the long pole on cross-solver parity studies (~5e-4 rel modal error, all from this path); closing the gap requires proper apex-singular shape functions (Bedrosian / Zgainski-Bossavit / Wachspress) plus a pyramid-domain Gauss rule — tracked as a dedicated follow-up.

References

Degenerate hex → wedge / pyramid topology collapse: Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §8.8.1 (discussion of collapsed hex nodes producing wedge / pyramid shapes).
15-node serendipity wedge shape functions (base corners on barycentric (ξ₁, ξ₂, ξ₃) = L₁, L₂, L₃ with ξ₃ = 1 − ξ₁ − ξ₂, quadratic through thickness in ζ): Bedrosian, G., “Shape functions and integration formulas for three-dimensional finite element analysis,” IJNME 35 (1992), pp. 95–108. Equivalent derivation in Zienkiewicz & Taylor §8.8.3 Table 8.8-1.
9-point wedge Gauss rule (3-point triangular rule × 3-point Gauss-Legendre through thickness): triangular rule from Dunavant, D.A., “High degree efficient symmetrical Gaussian quadrature rules for the triangle,” IJNME 21 (1985), pp. 1129–1148 (degree-4 triangle rule on 3 symmetric points); tensor-product with the 1-D Gauss-Legendre rule follows Zienkiewicz & Taylor §5.10.
21-point wedge mass rule (6-point triangular rule × 3-point through-thickness, integrates the degree-4 Bᵀ B polynomial on the wedge): Dunavant degree-4 triangle rule × Gauss-Legendre 3.
13-node pyramid shape functions (Pyr13, currently a hex-fold wrapper pending a dedicated kernel): Bedrosian 1992 (as above) for the apex-singular polynomial basis, alternative formulations in Zgainski, F.X., Coulomb, J.L., and Marechal, Y., “A new family of finite elements: the pyramidal elements,” IEEE Trans. Magn. 32 (1996), pp. 1393–1396 and Wachspress, E.L., A Rational Finite Element Basis, Academic Press, 1975.

class femorph_solver.elements._wedge15_pyr13.Pyr13[source]#

Bases: ElementBase

13-node degenerate Hex20 pyramid.

Default: dedicated Bedrosian apex-singular serendipity kernel integrated with a 2×2×2 Duffy collapsed-hex Gauss rule. The 2×2×2 rule (not 3×3×3) is the under-integrated quadrature that matches the foreign-deck convention; the higher-order "consistent" rule is mathematically more accurate but ~13 % apart in Frobenius norm on a single element.

The hex-fold wrapper (T^T · K_20-node · T) is retained as material["_PYR13_INTEGRATION"] = "hex_fold" for parity tests against foreign decks that emit hex-folded pyramids.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

class femorph_solver.elements._wedge15_pyr13.Wedge15[source]#

Bases: ElementBase

15-node degenerate Hex20 wedge (K=L, O=P collapse).

Dedicated quadratic serendipity kernel with shape functions in area coords (ξ₁, ξ₂, ζ), stiffness at 9-pt Gauss (3-pt triangle × 3-pt line), consistent mass at 21-pt Gauss (7-pt triangle × 3-pt line).

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Tet10 — 10-node quadratic tetrahedral solid.

Node ordering matches VTK_QUADRATIC_TETRA:

0-3 : corner nodes I, J, K, L 4 : mid-edge I-J 5 : mid-edge J-K 6 : mid-edge K-I 7 : mid-edge I-L 8 : mid-edge J-L 9 : mid-edge K-L

Reference tetrahedron with natural coordinates (ξ, η, ζ), volume coordinates (L₁, L₂, L₃, L₄) = (1-ξ-η-ζ, ξ, η, ζ):

Corner node i (Lᵢ = 1 at corner i): Nᵢ = Lᵢ(2Lᵢ - 1)
Mid-edge (i, j): N_{ij} = 4 Lᵢ Lⱼ

Both stiffness and mass use the 4-point Gauss rule on the reference tet (degree-2 exact, weights 1/24 at the 4 symmetric Keast points). The 4-pt rule integrates K exactly (degree-2 integrand from BᵀCB on quadratic shape functions); the per-element M is degree-4 and therefore rank 12 of 30 (PSD, not PD) under the same rule. That’s fine: once two elements share nodes the global M becomes PD, and lumped M is diagonal regardless.

Real constants#

None required — geometry comes from node coordinates.

Material#

EX, PRXY : elastic moduli (both mandatory) DENS : density (required for M_e)

References

Shape functions — 10-node quadratic Lagrange tetrahedron in volume coordinates (L₁, L₂, L₃, L₄): Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §8.8.2. Also Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §6.5 Table 6.5-1.
4-point Keast Gauss rule (degree-2 exact on the unit tetrahedron, weights 1/24 at the four symmetric points): Keast, P., “Moderate-degree tetrahedral quadrature formulas,” Comput. Methods Appl. Mech. Engrg. 55 (1986), pp. 339–348. The 4-point rule integrates K exactly (degree-2 integrand from BᵀCB on quadratic shape functions) and M with rank deficiency 12 of 30 that the global assembly absorbs.
Voigt elastic matrix: Cook et al. §3.7.

class femorph_solver.elements.tet10.Tet10[source]#

Bases: ElementBase

static eel_batch(coords: ndarray, u_e: ndarray, material: dict[str, float] | None = None) → ndarray | None[source]#

Per-element elastic strain at every element node.

coords: (n_elem, 10, 3); u_e: (n_elem, 30). Returns (n_elem, 10, 6) — engineering-shear Voigt strain. material is accepted for signature uniformity with plane kernels but is unused.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Beam2 — 3D 2-node linear beam.

femorph-solver implements the slender-beam (Euler–Bernoulli) form of Beam2: axial, torsional, and two independent bending planes combine into a 12×12 local stiffness matrix rotated into global coordinates. This is exactly the limit Beam2 converges to when the Timoshenko shear parameter Φ → 0, i.e. for slender sections where bending dominates over shear deformation. A future revision can introduce Φ from an explicit shear area or section integration.

Local DOF ordering at each node:

(u, v, w, θx, θy, θz)

so global assembly fills rows/cols 0..5 (node i) and 6..11 (node j) in the canonical (UX, UY, UZ, ROTX, ROTY, ROTZ) DOF order.

Real constants#

real[0]: AREA — cross-sectional area A real[1]: IZZ — bending inertia about local z (resists y-deflection) real[2]: IYY — bending inertia about local y (resists z-deflection) real[3]: J — Saint-Venant torsion constant

Material#

EX : Young’s modulus PRXY : Poisson’s ratio (→ G = E/(2(1+ν)) for torsion) DENS : mass density (required for M_e)

Orientation#

The local x-axis points from node I → node J. Without an explicit orientation node, the local y-axis is chosen so that local z stays “up” (world +Z) whenever the beam is not parallel to the global Z-axis, falling back to world +Y for vertical beams. For axisymmetric sections (IYY == IZZ) the orientation is immaterial.

Stiffness (slender-beam limit)#

K_local = [ K_axial · · · ]
[ · K_bend_xy · · ] [ · · K_bend_xz · ] [ · · · K_torsion ]

with the familiar Hermite-cubic 4×4 bending blocks. A block-diagonal direction-cosine matrix Λ rotates it: K_global = Λᵀ K_local Λ.

Mass (consistent, Hermite cubic)#

Axial / torsional: ρ A L / 6 · [[2,1],[1,2]] (resp. ρ I_p L / 6 with I_p ≈ I_y + I_z). Bending uses the standard 4×4 block ρ A L / 420 · (…) in both planes.

Lumped mass puts ρ A L / 2 on each translational DOF and zeros the rotational ones — a simple diagonal approximation useful for explicit dynamics but less accurate for modal work than the consistent form.

References

Euler–Bernoulli beam (slender-beam limit, shear neglected): Timoshenko, S.P., Vibration Problems in Engineering, 4th ed., Wiley, 1974, §5.3. Beam kinematics and derivation of the 12-DOF stiffness / consistent-mass blocks: Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §2.4 (axial), §2.5 (bending Hermite cubics), §2.6 (torsion).
Hermite cubic shape functions for transverse displacement and slope continuity at the nodes: Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §2.5.1 eqs. (2.26)–(2.27).
Consistent mass matrix from Hermite cubics (coefficients ρAL/420 on the 4×4 bending block and ρAL/6 on the axial / torsional block): Cook et al. Table 16.3-1.
Timoshenko-beam form (non-zero shear, Φ = 12EI/(κAGL²) — future revision) and comparison with the Euler-Bernoulli limit: Bathe, K.J., Finite Element Procedures, 2nd ed., Prentice Hall, 2014, §5.4.

class femorph_solver.elements.beam2.Beam2[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Quad4Shell — 4-node flat Mindlin–Reissner shell.

Kinematics (first-order shear-deformation theory):

u_point = u + z · θ_y v_point = v − z · θ_x w_point = w

giving a split into membrane / bending / transverse-shear energy after through-thickness integration. Node ordering matches VTK_QUAD:

(-,-), (+,-), (+,+), (-,+) in natural (ξ, η) ∈ [-1, 1].

The element is built in a local 2-D frame (ex along edge 1→2, ez normal to the best-fit plane, ey = ez × ex). Coordinates are projected into this frame, the plane stiffness is assembled, then a block-diagonal 3×3 direction-cosine matrix rotates every node’s six-DOF block back to global coordinates.

Integration#

Membrane + bending: 2×2 Gauss (full)
Transverse shear: 1×1 Gauss (selective reduced) — suppresses shear locking for thin shells without needing MITC4 tying.

Drilling DOF (θ_z, local) has no true stiffness for a flat shell; femorph-solver adds a small stabilization α · G · t · A_elem per corner so the 24×24 local matrix stays non-singular when every DOF is free.

Real constants#

real[0]: THICKNESS — shell thickness (mandatory)

Material#

EX, PRXY : plane-stress elastic moduli (mandatory) DENS : mass density (required for M_e)

References

Mindlin–Reissner first-order shear-deformation theory (transverse shear included so rotations are independent of in-plane slope): Mindlin, R.D., “Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates,” J. Appl. Mech. 18 (1951), pp. 31–38. Reissner, E., “The effect of transverse shear deformation on the bending of elastic plates,” J. Appl. Mech. 12 (1945), pp. A69–A77.
Bilinear isoparametric 4-node shell element with 6 DOF / node: Zienkiewicz, O.C. and Taylor, R.L., The Finite Element Method, Vol. 2: Solid Mechanics, 5th ed., Butterworth-Heinemann, 2000, Ch. 8 (plate & shell elements).
Selective reduced integration for transverse shear (1×1 on shear, 2×2 on membrane + bending — suppresses shear locking for thin shells): Malkus, D.S. and Hughes, T.J.R., “Mixed finite element methods — reduced and selective integration techniques: a unification of concepts,” Comput. Methods Appl. Mech. Engrg. 15 (1978), pp. 63–81.
Alternative MITC4 (Mixed-Interpolated Tensorial Components) shear-locking treatment — roadmap for a future variant: Bathe, K.J. and Dvorkin, E.N., “A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation,” IJNME 21 (1985), pp. 367–383.
Drilling-DOF stabilisation (α·G·t·A per corner to keep the local 24×24 non-singular when θ_z is free): Allman, D.J., “A compatible triangular element including vertex rotations for plane elasticity analysis,” Comput. Struct. 19 (1984), pp. 1–8, and Hughes, T.J.R. and Brezzi, F., “On drilling degrees of freedom,” Comput. Methods Appl. Mech. Engrg. 72 (1989), pp. 105–121.

class femorph_solver.elements.quad4_shell.Quad4Shell[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Quad4Plane — 2D 4-node structural solid.

Quad element with bilinear shape functions on the reference square ξ, η ∈ [-1, 1]:

N_i(ξ, η) = (1/4)(1 + ξ_i ξ)(1 + η_i η)

Node ordering matches VTK_QUAD — the four corners traversed counter-clockwise: (-,-), (+,-), (+,+), (-,+). Two translational DOFs per node (UX, UY); element lives in the global (X, Y) plane.

material["_QUAD4_PLANE_MODE"] selects the constitutive regime:

"stress" (default) — plane stress (σ_zz = 0)
"strain" — plane strain (ε_zz = 0)
axisymmetric mode is roadmap and not yet implemented.

Stiffness / mass#

2×2 Gauss integration for both K and M:

K_e = ∫_A Bᵀ C B t dA   ≈ Σ_g Bᵀ(ξ_g) C B(ξ_g) det J(ξ_g) t_g
M_e = ρ t ∫_A Nᵀ N dA   ≈ Σ_g ρ t Nᵀ(ξ_g) N(ξ_g) det J(ξ_g)

where t is the out-of-plane thickness (real constant, default 1.0 for plane strain) and C is the 3×3 plane elastic matrix in engineering shears.

Real constants#

real[0] : THK — out-of-plane thickness (default 1.0)

Material#

EX, PRXY : elastic moduli (both mandatory) DENS : density (required for M_e)

Lumped mass#

Row-sum (HRZ) lumping distributes each row of the consistent mass onto its diagonal.

References

Shape functions — bilinear Lagrange on the 4-node quadrilateral (Q4 / tensor-product Lagrange on the reference square): Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §6.3.2. Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §6.3.
2×2 Gauss-Legendre quadrature (product rule): Zienkiewicz & Taylor §5.10 Table 5.3.
Plane-stress / plane-strain elastic matrix: Cook et al. §3.5 (plane stress) and §3.6 (plane strain).
Row-sum / HRZ lumped mass: Hinton, Rock & Zienkiewicz, “A note on mass lumping and related processes in the finite element method,” Earthquake Engng. Struct. Dyn. 4 (1976), pp. 245–249.

class femorph_solver.elements.quad4_plane.Quad4Plane[source]#

Bases: ElementBase

static eel_batch(coords: ndarray, u_e: ndarray, material: dict[str, float] | None = None) → ndarray | None[source]#

Per-element engineering-shear Voigt strain at element nodes.

coords: (n_elem, 4, 3) or (n_elem, 4, 2). u_e: (n_elem, 8) — DOF layout [ux0, uy0, ux1, …].

Returns (n_elem, 4, 6) Voigt strain ordered [εxx, εyy, εzz, γxy, γyz, γxz] — matches the 3D-solid Voigt convention so compute_nodal_stress() can drive the same σ = C·ε contraction with the standard isotropic (6, 6) C.

For plane stress (default) we fill εzz = -ν/(1-ν)·(εxx+εyy) so σ_zz recovered via the 6×6 C is identically zero (the defining condition of plane stress). For plane strain we fill εzz = 0 (the defining condition of plane strain); σ_zz is then ν·(σ_xx + σ_yy) per Hooke’s law.

material must be provided when the kernel is operating in plane stress (the default) so ν is available; for plane strain the material argument is unused.

Strain is evaluated at each corner node (ξ_i, η_i) via the bilinear B-matrix — the same basis used in ke. This is the standard Q4 nodal recovery; the recovered strain is exact at the centroid and accurate to O(h²) at corner nodes for smooth solutions.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Truss2 — 3D 2-node spar (truss).

Linear shape functions in the natural coordinate s ∈ [-1, 1] for all three translations. Only axial strain contributes to internal work:

B = (1/L) [-l, -m, -n,  l, m, n]      (1 × 6)
K_e = (E·A/L) · [[ C, -C],
                 [-C,  C]]            where C_ij = d_i·d_j  (3 × 3)

with (l, m, n) the unit vector along the element axis.

Consistent mass:

M_e = ρ·A·L / 6 · [[2 I₃,   I₃],
                   [  I₃, 2 I₃]]

Lumped mass: M_e = ρ·A·L/2 · I₆.

Real constants#

real[0]: AREA — cross-sectional area (mandatory) real[1]: ADDMAS — added mass per unit length (optional; unused here) real[2]: ISTRN — initial strain (optional; unused here — stress path only)

Material properties#

EX : Young’s modulus (mandatory) DENS : mass density (required for M_e)

References

Shape functions — linear Lagrange on a 1-D 2-node element (N₁ = ½(1 − s), N₂ = ½(1 + s) on s ∈ [−1, 1]): Zienkiewicz, O.C. and Taylor, R.L. The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013, §2.3. Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §2.4.
Axial-only strain (spar / truss element, no bending / shear contribution, only longitudinal stiffness): Cook et al. §2.3.
Consistent mass ρAL/6 · [[2, 1], [1, 2]] on the two axial DOFs (and its 6×6 extension to 3-D for the three translational pairs): Cook et al. Table 16.3-1; Zienkiewicz & Taylor §16.2.1.
Row-sum lumped mass ρAL/2 per translational DOF: Hinton, Rock & Zienkiewicz, Earthquake Engng. Struct. Dyn. 4 (1976), pp. 245–249.

class femorph_solver.elements.truss2.Truss2[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Spring — 3D 2-node spring-damper (longitudinal mode only).

Default mode: longitudinal spring oriented along the I→J axis. Real constants:

real[0] : K — spring stiffness (force / length) real[1] : CV1 — linear damping (force · time / length, unused for K_e) real[2] : CV2 — non-linear (cubic) damping (unused here) real[3] : IL — initial length (unused; reference length is the nodal distance)

Stiffness:

K_e = K · [[ C, -C],
           [-C,  C]]        with C_ij = d_i · d_j  (3 × 3)

where (d_x, d_y, d_z) is the unit vector along I→J.

Spring is massless by convention (modal codes treat it as a stiffness-only contribution); me returns a 6×6 zero matrix so that assembly loops stay uniform.

Only the longitudinal mode is implemented; torsional and 2-D planar variants need separate code paths and will be added alongside their verification datasets.

References

Discrete spring / lumped connector element (3-D longitudinal mode): Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §2.3 and §2.10 (direction-cosine rotation of a scalar spring between two nodes in 3-D). Stiffness K · (d ⊗ d) on the I→J axis is the standard result.
Direction-cosine transformation for a bar / spring aligned with an arbitrary 3-D axis: Zienkiewicz, O.C. and Taylor, R.L., The Finite Element Method, 7th ed., 2013, §2.4.4.

class femorph_solver.elements.spring.Spring[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

PointMass — concentrated point mass.

femorph-solver implements the translational-only variant: a single-node element carrying three translational DOFs UX, UY, UZ. The rotary-inertia variant requires 6 DOFs/node and will be added alongside Beam2 once mixed-DOF assembly is in place.

Real constants#

real[0]: MASSX — translational mass along x (mandatory) real[1]: MASSY — translational mass along y (defaults to MASSX) real[2]: MASSZ — translational mass along z (defaults to MASSX)

Matrices#

K_e = 0₃×₃ (a point mass has no stiffness contribution) M_e = diag(MASSX, MASSY, MASSZ)

Lumped and consistent formulations coincide for a one-node element.

References

Concentrated / point-mass element (lumped nodal mass, no stiffness contribution): Cook, Malkus, Plesha, Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, 2002, §16.3 (discussion of lumped nodal masses). Bathe, K.J., Finite Element Procedures, 2nd ed., Prentice Hall, 2014, §4.2.4.
Anisotropic diagonal mass tensor (diag(MASSX, MASSY, MASSZ), supports different masses on each translational axis): standard generalisation of the isotropic m·I₃ lumped mass.

class femorph_solver.elements.point_mass.PointMass[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with property names as keys ("EX", "PRXY", "DENS", …))
real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

Post-processing recovery#

Public post-processing helpers — strain, stress, invariants.

A stable, non-underscore home for the per-node recovery helpers a gateway / downstream consumer needs to expose results without reaching into internal modules. Every name re-exported here is the canonical public binding; the underscore paths (femorph_solver.result._stress_recovery, femorph_solver.elements._stress) remain for backwards compatibility but should not be imported by new code.

Examples

>>> import femorph_solver as fs
>>> from femorph_solver.recover import (
...     compute_nodal_strain,
...     compute_nodal_stress,
...     stress_invariants,
... )
>>> # m: fs.Model, u: (n_dof,) DOF-indexed displacement
>>> eps = compute_nodal_strain(m, u)        # (n_points, 6) Voigt strain
>>> sig = compute_nodal_stress(m, u)        # (n_points, 6) Voigt stress
>>> inv = stress_invariants(sig)            # von_mises, s1/s2/s3, ...

femorph_solver.recover.compute_nodal_strain(model: Model, displacement: np.ndarray) → np.ndarray[source]#

Compute (n_points, 6) Voigt strain at every mesh node.

Per-element strain comes from Model.strain(), which evaluates \(\varepsilon = B(\xi_\text{node})\,u_e\) directly at each element’s own nodes. Where multiple elements touch a node the strains are averaged with equal weight (the same convention as the stress recovery in compute_nodal_stress() and as MAPDL’s PRNSOL).

Parameters:

model (femorph_solver.model.Model) – Source of the strain kernel (model.strain(u)). Element cells / nodes are read off model.grid.
displacement (numpy.ndarray) – (N,) DOF-indexed displacement in Model.dof_map() ordering. Typically a static-solve displacement or a modal mode shape.

Returns:

(n_points, 6) engineering-shear Voigt strain [ε_xx, ε_yy, ε_zz, γ_xy, γ_yz, γ_xz], averaged across every element that touches a node. Nodes with no incident elements are zero.

Return type:

numpy.ndarray

Validation framework#

The public surface is re-exported from femorph_solver.validation; the private _benchmark, _convergence, _report modules are not separately autodoc’d to avoid duplicate cross-reference targets in the rendered docs.

Validation framework — canonical FEA problems with published references.

femorph-solver’s validation pillar is a catalogue of benchmark problems whose expected result is a published value — a closed-form solution from a textbook, a handbook table, or a peer-reviewed paper — not another solver’s output. Each benchmark carries its source citation in-line so a reviewer can trace every tolerance back to the page it was derived on.

The framework provides three things:

A declarative BenchmarkProblem interface that couples a mesh-building recipe to a catalogue of PublishedValue references.
A ConvergenceStudy runner that sweeps a problem across mesh refinements, fits a convergence rate, and emits a ConvergenceRecord suitable for regression testing, verification-gallery rendering, or JSON output.
A _report module that writes JSON + Markdown reports consumable by the docs pipeline and by CI regression guards.

Usage:

from femorph_solver.validation import ConvergenceStudy
from femorph_solver.validation.problems import CantileverEulerBernoulli

study = ConvergenceStudy(
    problem=CantileverEulerBernoulli(),
    refinements=[(20, 3, 3), (40, 3, 3), (80, 3, 3)],
)
records = study.run()
for r in records:
    print(r.quantity_name, r.rel_error, r.convergence_rate)

The problems submodule is the catalogue; new problems are registered there. Every problem mirrors a regression test under tests/analytical/ so the verification pillar stays green when the docs-facing catalogue does.

class femorph_solver.validation.BenchmarkProblem[source]#

Bases: ABC

Declarative canonical FEA problem with published reference(s).

Subclasses declare:

name — short identifier, snake_case.
description — one-sentence summary.
published_values — tuple of PublishedValue.
build_model() — construct a ready-to-solve model.
extract() — pull a named quantity from the solve output.

The base class provides validate() which composes build_model → solve/modal_solve → extract and returns a list of ValidationResult.

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

abstractmethod build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

abstractmethod extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

validate(**mesh_params: Any) → list[ValidationResult][source]#

Build, solve, compare against every published value.

Returns one ValidationResult per published_values entry.

class femorph_solver.validation.ConvergenceRecord(quantity_name: str, results: list[ValidationResult], convergence_rate: float | None = None)[source]#

Bases: object

Measured convergence of one quantity across refinements.

results is the per-refinement list. convergence_rate is the fitted \(p\) in \(|e| = C \cdot n_\text{dof}^{-p}\) over the last two refinements (coarsest pair drops out so the rate reflects the fine-mesh regime). None if fewer than two refinements completed or the error did not monotonically decrease.

class femorph_solver.validation.ConvergenceStudy(problem: BenchmarkProblem, refinements: list[dict[str, Any]] = <factory>)[source]#

Bases: object

Run a benchmark at multiple refinements, capture + fit.

Parameters:

problem (BenchmarkProblem) – A BenchmarkProblem instance.
refinements (list[dict[str, Any]]) – Sequence of mesh-parameter dicts. Each dict is forwarded verbatim to problem.build_model. The order is coarse → fine; the final entry is the reference mesh for acceptance.

run() → list[ConvergenceRecord][source]#: Run every refinement; return one record per published quantity.

class femorph_solver.validation.PublishedValue(name: str, value: float, unit: str = '', source: str = '', formula: str = '', tolerance: float = 0.05)[source]#

Bases: object

One published reference value associated with a benchmark.

The source string must uniquely identify the publication: "Timoshenko 1955 §5.4", "Rao 2019 Table 8.1", "Irons & Razzaque 1972" — enough that a reviewer can find the citation in the reference list without further lookup.

tolerance is the acceptance tolerance for this quantity at the recommended mesh — looser on coarse meshes, tighter on finer ones. Convergence studies may use a stricter tolerance_fn(n_dof) later; the scalar is the floor.

class femorph_solver.validation.ValidationResult(problem_name: str, published: ~femorph_solver.validation._benchmark.PublishedValue, computed: float, mesh_params: dict[str, ~typing.Any] = <factory>, n_dof: int | None = None, wall_s: float | None = None, peak_rss_mb: float | None = None)[source]#

Bases: object

One validation attempt at a specific mesh refinement.

rel_error is the signed relative error ((computed - published) / published). passed is set when the absolute relative error is within the PublishedValue’s tolerance. All other fields are the raw inputs so the caller can re-render the comparison without re-solving.

wall_s and peak_rss_mb capture the wall-clock time and peak resident-set size of the solve that produced this row, so validation outputs can feed the femorph_solver.estimators training loader — the same way benchmark rows do.

Public alias for the problem catalogue.

Importers should use:

from femorph_solver.validation.problems import CantileverEulerBernoulli

not the private _problems path. The catalogue re-exports everything at this module level so the validation framework has a single canonical public surface.

class femorph_solver.validation.problems.AxialRodNaturalFrequency(name: str = 'axial_rod_nf', description: str = 'Fixed-free slender rod, fundamental axial natural frequency f_1 = (1/4L) sqrt(E/rho) (Rao 2017 §8.2).', analysis: str = 'modal', default_n_modes: int = 1, L: float = 1.0, area: float = 0.0001, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Fixed-free axial-rod fundamental natural frequency.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Rod length [m].

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

area: float = 0.0001#: Cross-sectional area [m²] (any positive value works — the axial natural frequency is independent of A).

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

nu: float = 0.3#: Poisson’s ratio (unused by the truss kernel but required by the material stamp).

rho: float = 7850.0#: Density [kg/m³].

class femorph_solver.validation.problems.CantileverEulerBernoulli(name: str = 'cantilever_eb', description: str = 'Slender clamped cantilever under transverse tip load — Euler-Bernoulli closed-form (Timoshenko 1955 §5.4).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P: float = 1000.0)[source]#

Cantilever tip deflection / root stress vs Euler-Bernoulli.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Cantilever length [m].

P: float = 1000.0#: Total transverse tip load [N].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

nu: float = 0.3#: Poisson’s ratio.

rho: float = 7850.0#: Density [kg/m^3] — irrelevant for the static solution but required by the material stamp.

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.CantileverHigherModes(name: str = 'cantilever_higher_modes', description: str = 'Clamped-free prismatic beam — 2nd and 3rd transverse bending modes. Euler-Bernoulli closed-form with β_n L from the characteristic equation (Rao 2017 §8.5).', analysis: str = 'modal', default_n_modes: int = 12, L: float = 4.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Cantilever second + third bending natural frequencies.

L: float = 4.0#: Cantilever length [m]. Set to 4.0 (h/L = 1/80) so the Euler-Bernoulli closed form is asymptotically exact for the first four bending modes — the FE result tracks the EB prediction to within ~0.2 % across all four modes at the default mesh, instead of drifting up to ~4 % from the Timoshenko shear correction that a stockier beam picks up.

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.CantileverNaturalFrequency(name: str = 'cantilever_nf', description: str = 'First transverse bending frequency of a clamped-free prismatic beam — Rao 2017 §8.5 Table 8.1.', analysis: str = 'modal', default_n_modes: int = 10, L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Cantilever first bending natural frequency (Rao 2017 Table 8.1).

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

default_n_modes: int = 10#: Request 10 modes so the extract can filter for the first bending mode even if axial / torsional modes appear below it.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.CantileverTipMoment(name: str = 'cantilever_tip_moment', description: str = 'Slender clamped cantilever under a pure tip moment — Euler-Bernoulli closed-form (Timoshenko 1955 §5.4).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, M: float = 50.0)[source]#

Cantilever tip deflection + rotation under a pure tip moment.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Cantilever length [m].

M: float = 50.0#: Total tip moment [N·m] (about the horizontal neutral axis).

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

nu: float = 0.3#: Poisson’s ratio.

rho: float = 7850.0#: Density [kg/m^3] — irrelevant for the static solution but required by the material stamp.

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.CantileverTorsion(name: str = 'cantilever_torsion', description: str = 'Cantilever beam with a pure tip torque — Saint-Venant closed-form tip-twist φ = T L / (G J).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, T: float = 10.0)[source]#

Cantilever tip twist under a pure tip torque.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Cantilever length [m].

T: float = 10.0#: Tip torque about the beam axis [N·m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m] (t in the Saint-Venant β table).

width: float = 0.05#: Cross-section width [m] (b in the Saint-Venant β table).

class femorph_solver.validation.problems.CantileverUDL(name: str = 'cantilever_udl', description: str = 'Slender clamped cantilever under uniformly-distributed transverse load — Euler-Bernoulli closed-form (Timoshenko 1955 §5.4).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, q: float = 1000.0)[source]#

Cantilever tip deflection under uniform distributed load.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Cantilever length [m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

nu: float = 0.3#: Poisson’s ratio.

q: float = 1000.0#: Uniform distributed load along the span [N/m] (positive value, acts in -z).

rho: float = 7850.0#: Density [kg/m^3] — irrelevant for the static solution but required by the material stamp.

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.ClampedClampedBeamCentralLoad(name: str = 'cc_beam_central_load', description: str = 'Prismatic clamped-clamped beam with central transverse point load — Euler-Bernoulli closed-form (Timoshenko 1955 §5.7).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P: float = 1000.0)[source]#

Clamped-clamped beam mid-span deflection under central point load.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Beam length [m].

P: float = 1000.0#: Total central transverse load [N] (acts in -z).

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

nu: float = 0.3#: Poisson’s ratio.

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.ClampedClampedBeamModes(name: str = 'cc_beam_modes', description: str = 'Clamped-clamped prismatic beam, first three transverse bending natural frequencies — Rao 2017 §8.5 Table 8.1.', analysis: str = 'modal', default_n_modes: int = 8, L: float = 4.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

CC beam first three bending natural frequencies.

L: float = 4.0#: Beam span [m]. Set to 4.0 (h/L = 1/80) so the Euler- Bernoulli closed form is asymptotically exact for the first three bending modes — at the original L = 1 (h/L = 1/20) the 3D solid mesh truthfully captures a Timoshenko shear / rotary-inertia correction that grows up to ~4 % at mode 3.

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.ClampedPlateStatic(name: str = 'clamped_plate_static', description: str = 'Clamped square Kirchhoff plate under uniform pressure — NAFEMS LE5 / Timoshenko & Woinowsky-Krieger 1959 §32.', analysis: str = 'static', a: float = 1.0, h: float = 0.02, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, q: float = 100000.0)[source]#

Clamped square plate under uniform pressure (NAFEMS LE5).

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.FreeFreeRodModes(name: str = 'free_free_rod_modes', description: str = 'Free-free uniform rod axial modes — f_n = n / (2L) sqrt(E/ρ) (Rao 2017 §8.2 Table 8.1).', analysis: str = 'modal', default_n_modes: int = 4, L: float = 1.0, area: float = 0.0001, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Free-free axial rod first elastic + second elastic frequencies.

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.IronsPatchTest(name: str = 'patch_test', description: str = 'Irons constant-strain patch test: prescribed linear boundary displacement must reproduce a uniform strain to machine precision.', analysis: str = 'static', L: float = 1.0, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, eps0: float = 0.0001)[source]#

2×2×2 hex patch under imposed uniform strain (Irons 1972).

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

eps0: float = 0.0001#: Target uniaxial strain imposed via boundary displacements. The published expectation is that every interior DOF reproduces u = eps0 * x at machine precision.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.KirchhoffSSPlateModes(name: str = 'ss_plate_modes', description: str = 'Simply-supported rectangular Kirchhoff plate fundamental bending frequency — Timoshenko & Woinowsky-Krieger 1959 §63.', analysis: str = 'modal', default_n_modes: int = 10, a: float = 1.0, b: float = 1.0, h: float = 0.01, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Simply-supported Kirchhoff plate fundamental frequency.

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

default_n_modes: int = 10#: Request 10 modes so the extract can filter for the (1,1) plate-bending mode even if spurious in-plane modes appear below it on coarse meshes.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.LameCylinder(name: str = 'lame_cylinder', description: str = 'Thick-walled cylinder under internal pressure — Lamé closed-form radial displacement + hoop stress at r = a (Timoshenko & Goodier 1970 §33).', analysis: str = 'static', a: float = 0.01, b: float = 0.02, t_axial: float = 0.001, E: float = 210000000000.0, nu: float = 0.3, rho: float = 7850.0, p_i: float = 100000000.0)[source]#

Thick-walled cylinder under internal pressure (plane strain).

E: float = 210000000000.0#: Young’s modulus [Pa].

a: float = 0.01#: Inner radius [m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

b: float = 0.02#: Outer radius [m].

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

p_i: float = 100000000.0#: Internal pressure [Pa].

t_axial: float = 0.001#: Axial thickness of the 3D slab [m] (plane-strain model).

class femorph_solver.validation.problems.NafemsLE1(name: str = 'nafems_le1', description: str = 'NAFEMS R0015 §2.1 LE1 — elliptic membrane under outward pressure, σ_yy at point D (stub of the inner ellipse).', analysis: str = 'static', a_in: float = 2.0, b_in: float = 1.0, a_out: float = 3.25, b_out: float = 2.75, thickness: float = 0.1, E: float = 210000000000.0, nu: float = 0.3, rho: float = 7800.0, pressure: float = 10000000.0)[source]#

NAFEMS LE1 — elliptic membrane, plane stress.

Reports σ_yy at point D = (2, 0), the “stub” of the inner ellipse. NAFEMS converged reference: 92.7 MPa.

E: float = 210000000000.0#: Young’s modulus [Pa].

a_in: float = 2.0#: Inner ellipse semi-axes [m].

a_out: float = 3.25#: Outer ellipse semi-axes [m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

pressure: float = 10000000.0#: Outward-normal pressure on the outer elliptic edge [Pa].

thickness: float = 0.1#: Out-of-plane thickness [m].

class femorph_solver.validation.problems.PinchedRing(name: str = 'pinched_ring', description: str = 'Thin circular ring under two opposed point loads — Castigliano closed-form diametrical deflection (Timoshenko & Young 1968 §79).', analysis: str = 'static', R: float = 0.1, width: float = 0.01, height: float = 0.005, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P: float = 10.0)[source]#

Pinched-ring diametrical deflection (quarter-symmetry).

E: float = 200000000000.0#: Young’s modulus [Pa].

P: float = 10.0#: Magnitude of each applied point load [N].

R: float = 0.1#: Ring mean radius [m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.005#: Cross-section height (radial-thickness) [m].

width: float = 0.01#: Cross-section width [m].

class femorph_solver.validation.problems.PlateWithHole(name: str = 'plate_with_hole', description: str = 'Quarter-symmetry plane-stress plate with a circular hole, remote uniaxial tension σ_∞. Kirsch (1898) stress concentration K_t = 3 at the hole edge, perpendicular to the remote tension axis.', analysis: str = 'static', a: float = 0.1, W: float = 1.0, thickness: float = 0.01, E: float = 210000000000.0, nu: float = 0.3, rho: float = 7850.0, sigma_infty: float = 10000000.0)[source]#

Plate with a circular hole under uniaxial tension — Kirsch K_t = 3.

E: float = 210000000000.0#: Young’s modulus [Pa].

W: float = 1.0#: Plate half-width [m] (W >> a for the infinite-plate limit).

a: float = 0.1#: Hole radius [m].

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

sigma_infty: float = 10000000.0#: Remote uniaxial tension [Pa].

thickness: float = 0.01#: Out-of-plane thickness [m] (plane stress).

class femorph_solver.validation.problems.ProppedCantileverCentralLoad(name: str = 'propped_cantilever_central_load', description: str = 'Fixed-pinned beam with central transverse point load — Euler-Bernoulli closed-form (Gere & Goodno 2018 §10.3 Table 10-1 Case 6).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P: float = 1000.0)[source]#

Propped cantilever mid-span deflection under central point load.

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.SSPlateStatic(name: str = 'ss_plate_static', description: str = 'Simply-supported square Kirchhoff plate under uniform pressure; checks the centre deflection against the Navier series — Timoshenko & Woinowsky-Krieger 1959 §30-§31.', analysis: str = 'static', a: float = 1.0, h: float = 0.02, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, q: float = 100000.0)[source]#

SS square plate under uniform pressure — centre-deflection check.

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

class femorph_solver.validation.problems.SimplySupportedBeamCentralLoad(name: str = 'ss_beam_central_load', description: str = 'Prismatic simply-supported beam with central transverse point load — Euler-Bernoulli closed-form (Timoshenko 1955 §5.6).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P: float = 1000.0)[source]#

Simply-supported beam mid-span deflection under central point load.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 1.0#: Beam length [m].

P: float = 1000.0#: Total transverse mid-span load [N] (acts in -z).

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

nu: float = 0.3#: Poisson’s ratio.

rho: float = 7850.0#: Density [kg/m^3] — irrelevant for the static solution.

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.SimplySupportedBeamModes(name: str = 'ss_beam_modes', description: str = 'Prismatic simply-supported beam, first three transverse bending natural frequencies — Euler-Bernoulli closed-form (Rao 2017 §8.5).', analysis: str = 'modal', default_n_modes: int = 6, L: float = 4.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0)[source]#

Simply-supported beam first three bending natural frequencies.

E: float = 200000000000.0#: Young’s modulus [Pa].

L: float = 4.0#: Beam length [m]. Set to 4.0 (h/L = 1/80) so the Euler- Bernoulli closed form is asymptotically exact across the first three bending modes — at L = 1 (h/L = 1/20) the 3D solid mesh truthfully captured a Timoshenko shear correction up to ~2.6 % at mode 3.

analysis: str = 'modal'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

height: float = 0.05#: Cross-section height [m].

width: float = 0.05#: Cross-section width [m].

class femorph_solver.validation.problems.SimplySupportedBeamUDL(name: str = 'ss_beam_udl', description: str = 'Prismatic simply-supported beam with uniformly-distributed transverse load — Euler-Bernoulli closed-form (Timoshenko 1955 §5.6).', analysis: str = 'static', L: float = 1.0, width: float = 0.05, height: float = 0.05, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, q: float = 1000.0)[source]#

SS beam mid-span deflection under uniformly-distributed load.

analysis: str = 'static'#: Analysis type — "static" or "modal". Subclasses override.

build_model(**mesh_params: Any) → Model[source]#: Construct a Model ready to solve.

extract(model: Model, result: Any, name: str) → float[source]#

Pull the quantity named name from the solve output.

result is whatever model.solve() / model.modal_solve() returned, unboxed. name is the PublishedValue name.

q: float = 1000.0#: Uniformly distributed load along the span [N/m] (acts in -z).

class femorph_solver.validation.problems.SingleHexUniaxial(name: str = 'single_hex_uniaxial', description: str = "Single 8-node hex under uniaxial tension — verifies 1-D Hooke's law and Poisson contraction to machine precision.", analysis: str = 'static', L: float = 1.0, E: float = 200000000000.0, nu: float = 0.3, rho: float = 7850.0, P_total: float = 1000000.0)[source]#