Nonlinear Solver (fem-nonlinear)

The fem-nonlinear crate provides an OpenSees-style iterative nonlinear static analysis framework. It is crate-independent - it defines traits that any structural element can implement, and algorithms that drive the Newton‑Raphson loop generically.

Crate design

fem-nonlinear/
├── element.rs         NonlinearElement trait (K_t, F_int, state management)
├── convergence.rs     ConvergenceCriterion trait + 6 concrete tests
├── sparse_solver.rs   SparseSolver trait + sprs-ldl wrapper (SPD + 1×1)
├── integrator/
│   └── mod.rs         Integrator trait, LoadControl, DisplacementControl
├── algorithm/
│   └── mod.rs         SolverAlgorithm trait, Linear, NewtonRaphson, ModifiedNewton
└── tests/
    └── integration_tests.rs

Traits

Trait	Purpose
NonlinearElement	Bridge from solver to element - tangent_stiffness(), resisting_force(), add_displacement(), commit(), revert()
Integrator	Assembles K and R, applies displacement updates, manages λ
SparseSolver	Factor-solve K·Δu = R (sprs-ldl for symmetric, extensible to Bunch-Kaufman)
ConvergenceCriterion	Test `
SolverAlgorithm	Runs the Newton loop: { tangent, solve, update, unbalance, test } while()

Algorithms

Algorithm	Behaviour
Linear	One-shot solve, no iteration
NewtonRaphson	Full tangent stiffness every iteration
ModifiedNewton	Factor K once, reuse for all iterations

Integrators

Integrator	What it controls
LoadControl	Fixed load increment Δλ per step; λ ramps from 0→1 in N steps
DisplacementControl	Prescribed displacement at a control DOF; solves for λ at each step

The Newton‑Raphson loop for LoadControl is 15 lines:

form unbalance (λ_target·P − F_int)
for iter in 0..max_iter:
    form tangent K_t(mode)
    factor·solve K_t·Δu = R
    update elements: u += Δu
    form unbalance
    test convergence → Converged / Continue / Failed
commit or revert elements

Key design decisions

DisplacementControl is NOT an Integrator impl. It uses an augmented system [K −p_ref; 1 0] [Δu; Δλ] = [−R; prescribed] and solves two RHS per iteration (Δu_hat and Δu_I). This requires custom loop control, so it's a standalone struct rather than going through the generic SolverAlgorithm.

NonlinearElement handles restrained DOFs. dof_indices() returns &[Option<usize>] where None means the DOF is restrained. Assembly skips these rows/columns, so the global system is non-singular without a separate constraint-handling layer.

Relative convergence is the default for real models. The SI-unit absolute tolerance problem (tol=1e-9 unreachable when forces are 1e6 N) is solved by RelativeNormUnbalance. The reference norm is latched by start_step(), which the algorithms call with the step's INITIAL unbalance before any correction - latching inside test() (post-correction) would pin the ratio at 1.0 on a near-linear step, which can then never converge.

Box delegation. impl NonlinearElement for Box<dyn NonlinearElement> allows heterogeneous element collections at integration time (beams + trusses in one solve).

Floating-point tolerance awareness. The proptest suite uses tol=1e-8 because the f64 round-trip error for k * (p/k) is ≈ ε·p - for large k,p values this exceeds 1e-12.

Test suite

28 tests total (4 unit + 24 integration):

Category	Tests	What's covered
Algorithm correctness	6	Linear, NR, Modified Newton on elastic and hardening elements
Convergence criteria	5	NormUnbalance, NormDispIncr, EnergyIncr, RelativeNormUnbalance, RelativeNormDispIncr
Integrators	8	LoadControl lambda progression, DisplacementControl drive-to-load, parallel DOFs, series springs, restrained DOF exclusion
State management	2	commit/revert called by framework on convergence and failure
Proptest invariants	4	NR = 1 iter for linear, displacement ∝ load, N steps = 1 step superposition
Framework plumbing	3	sprs-ldl n≥2 path, crate solver + criterion exercised, InitialThenCurrent behaviour

Integration with fem-core

The bridge lives in fem-core/src/nonlinear.rs: a private ElasticNlElement wrapper pairs each beam/truss element's global-coordinate stiffness with its reduced DOF map (None = restrained) and trial/committed displacement state, keeping the hot linear-solve path untouched.

solve_nonlinear_static(model, load_case, n_steps, max_iter, tol) applies the reference load in equal increments (LoadControl), equilibrating each step with full Newton-Raphson against RelativeNormUnbalance, and returns one AnalysisResults per converged step. Member and station forces are recovered against the lambda-scaled member loads, so mid-push internal forces are exact.

Deliberate v1 limits, all loud: the elastic tangent equals the initial stiffness (the smoke tests pin every step to scale exactly linearly with the load factor against the linear solver), plates are refused, capacity checks are not computed on this path (each step's warnings say so), and singular models error instead of being auto-restrained. Material and geometric nonlinearity land later as element capabilities behind the same trait - tension/compression-only members are the planned first real use.

Differences from XC / OpenSees

Aspect	XC / OpenSees	fem-nonlinear
Element interface	Virtual methods on DomainComponent	Rust trait (static dispatch for known types)
Constraint handling	Separate ConstraintHandler layer	Option<usize> in DOF indices (simpler, fewer allocations)
Solver back-end	SuperLU, UMFPACK, MUMPS, PETSc via polymorphic SOE	SparseSolver trait, currently sprs-ldl
Convergence test	NormUnbalance, NormDispIncr, NormEnergy + relative variants	Same tests + relative variants with self-latching
Tangent modes	CURRENT_TANGENT, INITIAL_TANGENT, INITIAL_THEN_CURRENT_TANGENT, etc.	Same three modes
Language	C++ with Fortran	Pure Rust (no unsafe, #![forbid(unsafe_code)])