Troubleshooting flowchart#

Decision tree for the most-frequently-encountered failure modes. Walk it top-to-bottom; the first match usually isolates the root cause.

Symptom: solver fails with “matrix is singular” / “factorization failed”#

  1. Did you fix at least three translations + three rotations’ worth of DOFs?

    • No — the model has surviving rigid-body modes. Run m.solve_modal(n_modes=6): any frequency near zero is an unconstrained DOF. Add the missing Model.fix() calls. See Boundary-condition pitfalls.

    • Yes — continue.

  2. Are any element Jacobians near zero?

    • Visualise det(J) per element (most CAD meshers export this directly; pyvista’s cell_quality() filter computes it). Anything below 0.1 of the mesh maximum is a candidate for refinement / remesh.

    • Refine or remesh the offending region. See Mesh quality.

  3. Are there massive material-property ratios in the model?

    • A steel insert in a rubber matrix with E ratio 1e4 pushes condition number above 1e10 routinely. Solve the parts separately if possible, or add a rigidity scaling factor on the soft material that matches the stiff one within an order of magnitude.

Symptom: solve succeeds but displacements are absurdly large#

This is a soft form of the singular-matrix case — the matrix isn’t quite singular but is close enough that the solver returns a near-null-space vector. The diagnostics:

  1. Run the rigid-body-mode check. m.solve_modal(n_modes=6) on the un-loaded model. Any near-zero frequency is the unconstrained DOF.

  2. Check the support-reaction sum. After m.solve_static(), sum StaticResult.reaction along each global axis. The sum must equal -ΣF_applied (Newton’s third law). See Reaction forces — global equilibrium and root moment from a tip load. If it doesn’t, you have a load-application bug — usually a missing accumulate=True on Model.apply_force() overwriting a previous load on the same node.

Symptom: stress field has spurious peaks at the clamp#

This is Saint-Venant edge effects, not a bug. The recovered nodal stress at a clamped boundary is artificially elevated because the clamp imposes a non-physical perfect constraint that the structure resists with concentrated local stress.

The closed-form analytical reference assumes a smooth distributed reaction; the FE clamp is a pointwise constraint on every face node, so the local stress is higher.

The fix:

  • For engineering-margin purposes, ignore the clamp- edge peak — read the stress one cell-row in from the clamp where Saint-Venant decay has cleaned the field.

  • For detailed-design purposes, model a small fillet / pad / pin at the clamp interface so the constraint imposes itself smoothly.

Symptom: modal solve returns first frequency near 0 Hz#

You have a surviving rigid-body mode. Either:

  • The model is intentionally free-free (a pressure vessel flying through space, a tuning-fork in vacuum). In that case the leading rigid-body modes are physical, and you’d pass them through any downstream response calculation unchanged — see Free-free axial rod — natural frequencies for the canonical free-free pattern.

  • The model should be constrained but isn’t — see the Boundary-condition pitfalls rigid-body diagnostic.

Symptom: solve is much slower than expected#

  1. Check thread pinning. echo $OMP_NUM_THREADS. If unset, the BLAS pool sizes itself to the host core count (often 32+) which thrashes 5-10× from scheduler contention. Pin to 4 on workstations, 8 on bigger nodes.

  2. Check backend selection. Result.diagnostics from the static solve reports which linear backend was chosen. If it’s SuperLU, install Pardiso or CHOLMOD — 2-3× speedup for free. See Conditioning and performance.

  3. Check matrix size. Model.stiffness_matrix().shape. Above 500 k DOFs, direct factorisation memory becomes the bottleneck — switch to MUMPS out-of-core or to an iterative solver (CG + AMG via linear_solver="cg+pyamg").

Symptom: modal frequencies don’t match the closed form#

  1. Check element type. HEX8 plain-Gauss locks badly on slender bending — the recovered mode is artificially stiff. Switch to ELEMENTS.HEX8(integration="enhanced_strain") (Simo-Rifai EAS).

  2. Check slenderness. At \(L/h < 20\), the Timoshenko shear correction the 3-D solid truthfully captures pushes frequencies a few percent below the Bernoulli closed form — that’s real physics, not a solver bug. See the cantilever-higher-modes verification: Cantilever beam — first four bending modes.

  3. Check the mesh-density rule of thumb — Mesh quality. At least 8 elements per wavelength of the highest mode you care about.

Where the documentation can’t help#

If none of the above matches, the issue is more likely a genuine bug in femorph-solver or a use case the library doesn’t yet support. Either:

  • File a GitHub issue with a minimal reproducer (a .py script that builds the model from scratch). The issue tracker is at femorph/solver#issues.

  • Ask the question with code on hand — the Result.diagnostics dict + the rigid-body-mode check + the support-reaction sum are usually enough information for a maintainer to spot the cause without re-running the model.