Detect global instability via the solution residual (fixes #255, #275)

[gaoflow]

Summary

Fixes #255 and #275. A globally unstable model — a mechanism (e.g. a beam with an internal hinge and no support at the hinge) or a structure with insufficient supports (rigid-body motion) — has a singular global stiffness matrix even though every diagonal term is non-zero. analyze() currently solves it anyway and silently returns erroneous reactions, shears, moments and displacements instead of raising.

The existing nodal stability check (Analysis._check_stability) only inspects the diagonal for unsupported zero-stiffness DOFs, so it cannot see this kind of global rank deficiency. The underlying solvers don't reliably report it either:

• scipy.sparse.linalg.spsolve returns NaN (with only a MatrixRankWarning), or in some cases finite-but-meaningless displacements with no warning at all (No error is thrown on unstable model in some scenarios #255 scenario 2).
• numpy.linalg.solve can miss the singularity due to round-off (No error is thrown on unstable model in some scenarios #255 scenario 3), and neither solver flags the hinge mechanism of Model produces results when unstable with a hinge #275.

Approach

I add Analysis._solve_unknown_disp(K11, rhs, sparse, check_stability), which solves for the unknown displacements and, when check_stability is enabled, verifies the solution against equilibrium using the relative residual ‖K11·x − rhs‖ / ‖rhs‖. A stable structure — even a numerically ill-conditioned one — solves to within round-off because the direct solvers are backward-stable, whereas a singular system yields a non-finite or O(1)-residual result.

I first tried the natural pivot-magnitude / rank approach (factorize and compare the smallest LU pivot to the largest). It is not sound for this codebase: several stable models in the existing test suite (e.g. the Kassimali 2D frame, and the mat-foundation / plate meshes) have a smallest-to-largest pivot ratio already down at ≈1e-17, indistinguishable from a genuine mechanism, so any pivot threshold misclassifies them. Measured across the suite, the worst stable relative residual is ≈2.4e-13, while the unstable models from #255/#275 give residuals of order 1 — a clean ~13-order-of-magnitude separation. The default tolerance is 1e-6.

The check is cheap (one matrix–vector product plus two norms, no extra factorization) and gated on the existing check_stability flag, so the historical fast path is preserved when it is disabled.

Scope

The first-order elastic / tension-compression path (Analysis._first_order, used by analyze()) and the linear path (FEModel3D.analyze_linear) now route through the helper. The P-Delta and pushover paths are intentionally left unchanged — they deliberately drive the structure toward buckling (where the tangent stiffness approaches singular by design), so a hard instability raise there would need separate, buckling-aware handling.

Reproduction & verification

All four reported scenarios now raise The stiffness matrix is singular ... instead of returning results, for both the sparse and dense solvers:

• Model produces results when unstable with a hinge #275: simply-supported two-span beam joined by an internal hinge with no support at the hinge.
• No error is thrown on unstable model in some scenarios #255 scenario 1: a single unsupported member.
• No error is thrown on unstable model in some scenarios #255 scenario 2: an unsupported frame loaded in-plane (previously returned finite garbage with no warning).
• No error is thrown on unstable model in some scenarios #255 scenario 3: an unsupported closed frame solved with sparse=False.

Added Testing/test_instability.py covering the hinge mechanism and the unsupported-model scenarios (sparse and dense), that a stable propped cantilever still solves to the correct reaction (5wL/8), and that check_stability=False bypasses the check.

The full core test suite passes (the optional-dependency test files for the visualization/VTK backends fail to collect locally for unrelated reasons and were excluded): 121 passed (115 existing + 6 new), with no regressions.

[JWock82]

Thanks for this excellent contribution. I made a few improvements to this with the help of GitHub Copilot and merged it all in PR #323.

[gaoflow]

Thanks for folding it into #323 — glad the global-instability detection was useful, and the Copilot polish on top sounds great. Happy to see it land.