fea-core

Imported from Teased-oChroid-orrA/powershell_tool (fea-core/AGENTS.md). Install upstream with npx skills add Teased-oChroid-orrA/powershell_tool --skill fea-core. Copyright stays with the author.

fea-core

TL;DR: General finite-element kernel: 2D (plane stress / plane strain / axisymmetric) and 3D solid elements (Tri3/6, Quad4/8/9, Tet4/10, Hex8/20/27), colour-parallel symmetric block assembly, Dirichlet elimination, faer sparse Cholesky (AMD or geometric nested dissection) or PCG + smoothed-aggregation AMG for large 3D, multi-point constraints and reference-point rigid coupling, isotropic and anisotropic elasticity, point / surface (constant or field) / body / thermal loads, nodal stress (averaged or SPR) with a ZZ error estimate, boundary-face extraction, .vtu writer, structured mesh generators. Pure function of inputs; no file I/O, no UI. Plan, validation and measured speed: docs/fea-core.md, decision: docs/adr/ADR-012-general-fea-kernel.md.

Purpose

Owns: element library, element kernels, assembly, linear solve, loads, stress recovery, mesh generators. Does not own: any UI, material tables (mechanics-core), the lug physics (lug-solver: FeaLug in lug-solver/src/fea*.rs runs the lug on this crate and is the toolbox's main solver; its condensed / finite-strain solvers are the legacy comparison, see docs/fea-core.md Phase 7). lug-solver depends on this crate, never the reverse. The app-tui FEA Workbench (problem definition, mesh preview, contours) is fea-problem + app-tui/src/toolboxes/fea_workbench/ (docs/fea-core.md Phase 8).

Code Map

Derivative engine (src/derivative/, epic #39)

Generic Scalar kernels (dual / reverse tape / hyper-dual / complex), a policy engine that picks the method with a reason, exact linear, plastic-path, contact, transient and second-order sensitivities. Read docs/derivative-engine.md (architecture, support matrix, limits) before changing it. Rules: derivatives are of the computed response, not the continuous one; never loosen a finite-difference tolerance to pass (use Richardson steps); refuse or flag nonsmooth points instead of returning a classical derivative; any new kernel registers its capabilities in kernels.rs; keep ordered (deterministic) reductions; CostCalibration defaults come from tests/derivative_bench.rs.

Contracts

• Long contact solves are steered through NlOptions: interrupt (Interrupt: wall-clock deadline and/or shared AtomicBool, checked at every Newton iteration and step; ends with Stop::Interrupted, not complete()), observer (StepObserver, called after every converged load step with a StepEvent: live progress). NlOptions is Clone, not Copy.
• stick_slip_guard (on in Tuning::friction(), off elsewhere) is for frictional contact whose Newton iteration cycles. Root cause found (eccentric bushing, 2026-10-08): points at the edge of a contact patch, at ~zero pressure, enter and leave the contact from one iteration to the next (NL_TRACE shows the active set of the first contact changing between iterations); each flip jumps the residual by more than the tolerance, so Newton alternates forever, and which way it falls depends on round-off (the case passed on macOS and failed on Linux CI). NlSolution::held_solves counts the solves that needed it (65 % of the random eccentric-bushing cases: it is the workhorse, not a rare net). When 8 iterations bring no 3 % gain of the best residual the solve holds the active set (ContactSet::hold_active_set: a point is in contact iff it was at that iterate, its pressure following the smooth branch; released when the solve ends) and gets 15 more iterations: every case that had failed all attempts (fine mesh, concentric thick wall, soft parts, the reported 8856 lbf case) then converged on the first attempt, 2-10 x faster. Only if that also stagnates: with stall_tol the lowest iterate is accepted when its residual is below stall_tol x force scale (NlSolution::stalled_solves; the force scale includes large constant forces, so check what the answer must still satisfy), else the solve fails (STAGNATED, 3 step cuts). step_memory (each Newton step started from the last accepted backtrack fraction) was removed: it rescued one case, ruined another, and no case needed it once the hold existed.
• Node order is VTK (corners, edge midsides, face centres, volume centre); a mesh writes to .vtu without renumbering. Boundary faces/edges are listed counter-clockwise seen from outside (2D: body on the left): Pressure(p) pushes into the body.
• Units are free (the kernel is unit-free); conventionally inch / psi / lbf. Plane problems carry a thickness in Physics; axisymmetric loads and stiffness are for the full 360 degrees (x = radius).
• Isotropic law is sigma = lambda tr(eps) I + 2 mu eps in every case; plane stress uses lambda* = E nu / (1 - nu^2). Thermal: sigma_th = -(3 lambda + 2 mu) alpha dT (plane stress: 2 lambda* + 2 mu).
• Matrices are the lower block triangle of d x d node blocks (BlockMatrix, Pattern); a node no element touches still gets a diagonal block and must be constrained.
• Assembly writes disjoint blocks per colour class through a raw pointer (assembly.rs, one unsafe block): do not break the invariant "elements of one colour share no node".
• Constrained dofs are eliminated, never penalised (Dirichlet in linear.rs, multi-point constraints by T^T K T in mpc.rs; extra dofs >= mesh.n_dofs() carry reference-point translations/rotations and have no stiffness of their own).
• Anisotropic Elastic (aniso: Some) uses the B-matrix kernel and Layout::NodeBlock; e/nu are then only a Voigt-average display value, so never read mat.e/mat.nu/lambda()/mu() outside the isotropic branches (stiffness, stiffness_fast, thermal_load, stress_from_strain, complementary_energy all branch on aniso). Plane problems need a material symmetric about z; axisymmetric axes are (r, z, theta).
• solve_static / solve_static_many refuse a model that leaves a rigid-body motion of any connected piece free (Model::check_constrained: the rigid modes against the constrained dofs): a singular system can factorise anyway and return the answer plus an arbitrary rigid motion. A node no element touches is ignored by the check and still needs its own constraint. Reference-point / MPC solves (mpc.rs) and the nonlinear driver do their own constraining and do not call it.
• SolveMethod::Auto factors (faer) except 3D systems whose factor holds AUTO_ITERATIVE_FACTOR_NNZ (60M) entries or more (looked at from AUTO_ANALYSE_DOFS = 20k free dofs) or exceeding AUTO_ITERATIVE_DOFS (150k), where it runs PCG + AMG; an unconverged iteration is an Err, never a wrong answer, and Auto falls back to the direct solve (Solution::fallback says why). The numeric factorisation is parallel from 1M (3D) / 8M (2D) factor entries (linear.rs).

Pitfalls

• Eigen residuals (dynamics.rs): the reported residual is a backward error less a 100 eps max(K_ii) floating-point allowance, so a rigid-body mode of a free structure (K x = round-off) reads ~0 rather than ~1; it is a verification, not a convergence flag, and a pair above the tolerance is an Err. Buckling shifts below the smallest positive factor (the proof is a successful Cholesky of K - sigma B); a pencil with a single sign needs no shift trouble, a bending stress state (both signs) does.

• Temperature fields: Loads::temperature is a nodal field on top of the uniform delta_t; stresses must be recovered with the same field (gauss_stresses_field / nodal_stresses_field), otherwise the thermal strain is missing from the stress (a test shows the error is the full E alpha dT). The nonlinear driver refuses it with plasticity.

• Tiny strains in finite-strain laws: write the law in the displacement gradient (H), not in F = I + H: F - F^-T loses eps absolute digits and Newton cannot reach a relative tolerance below G eps / force (neo-Hookean fixed this way; the Hencky J2 path has the same floor in principle).

• stagnated_cuts (frictional stick_slip_guard) counts per load step and resets after a converged step.

• Serendipity and quadratic faces: consistent nodal forces of a uniform traction are not equal (corners carry negative share on Quad8/Hex20). Apply end loads as SurfaceLoad::Traction on a face, never as equal nodal forces (a test did this wrong and showed +4.6 % tip deflection).

• Plane strain free thermal expansion is (1 + nu) alpha dT, not alpha dT (eps_zz = 0).

• Hex8 under self weight is only exact to ~5e-7: the exact field needs x^2 terms trilinear elements lack.

• Tri3 (CST) and Tet4 are poor on pressure-driven problems (Lame at 16x16: 2 % displacement error); they are in the library for meshers, not recommended for results.

• Ordering: AMD wins in 2D; geometric nested dissection wins 2-3x on trilinear hexahedra, ties on quadratic hex/tets. Ordering::Auto analyses both in 3D and keeps the sparser factor. faer solve is parallel by default (global rayon parallelism).

• Dev profile is optimised ([profile.dev.package.fea-core] in the root Cargo.toml).

• Speed: element stiffness is not the bottleneck (< 1 % of a 3D linear solve); the factorization is. A register-tiled GEMM-style stiffness kernel and pulp runtime SIMD dispatch (x86 Windows builds are SSE2-only by default) are deferred until the Phase 3 nonlinear loop re-assembles every Newton iteration and a profile shows assembly matters.

• SPR needs nt + 2 samples per patch (a two-triangle Tri6 patch has 6 samples for 6 unknowns and interpolates noise); boundary nodes take an interior neighbour's polynomial (one-sided boundary patches are inaccurate). The ZZ estimate is only meaningful (effectivity ~1) on smooth solutions; near a stress concentration the coarse-mesh estimate overshoots.

• Mesh grading toward a hole: an exponent above ~1.2 makes the Jacobian vanish at the hole edge and ruins the first element ring (found validating Kirsch).

• AMG: 2D and strongly anisotropic-aspect meshes converge slowly (~100 iterations); a strength-of-connection filter cuts iterations but raises setup (coarsening stalls into a big coarse factorization), so default is no filter. AmgOptions is the tuning surface.

• Test-design lesson: a quadratic displacement field is only an FE-exact solution if its -div(sigma) is applied as body force.

• Meshing: boundary vertices come from the geometry's curves (not the chord); Region outer loop is CCW and holes CW whatever order is given. The mesher conforms segments by splitting (never by flipping), so every boundary subsegment is a triangle edge; a line piece splits at the midpoint of its actual chord, an arc piece at the curve point. Ray casting must use one closed polyline (resampling each segment independently leaves a rounding gap that swallows a crossing).

• Quad meshes come from splitting each triangle into three quads (Catmull-Clark style): quad edge ~ half the triangle edge, so adapt::adapted_size_field(.., quad_split = true) triples the element area to recover the triangle size.

• Importers verify quadratic node order geometrically (midside nodes within 0.3 element sizes of the position their corners imply); keep that check when adding element types, and build test files from the format's own edge/face definitions, not from the permutation table under test.

• Hex sweeps cannot touch the revolution axis (x = 0 nodes would create zero-length edges); use axisymmetric elements there.

• Nonlinear: Block::plasticity (J2) switches a block to the plastic path; the linear solve_static ignores it. Elastic blocks in a nonlinear model use the linear kernels (small strain), so a finite-strain analysis needs J2::elastic_finite() on its elastic blocks. Stress/tangent conventions: S = Cauchy (small) / first Piola-Kirchhoff (finite), A = dS/dH as a 9x9 [3i+J][3k+L]; the 4th-order tangent contracts with the shape-function gradients (plus the hoop entry u_r / r for axisymmetry). Plane-stress convergence of sigma_zz = 0 must be scaled by the elastic stress scale, never by the yield stress (an elastic law has sigma_y = 1e30; that bug shipped to the first test run).

• Displacement control needs the tangent predictor (done in Ctx::newton); arc-length uses a secant predictor (the sign of K^-1 f is meaningless where the tangent is singular) and falls back to LDL^T when Cholesky fails; load control must NOT fall back (an indefinite tangent there is an unstable equilibrium and Newton would converge to it).

• The nonlinear driver factorizes with direct faer only; do not point it at the AMG path without a plan for indefinite tangents.

• Contact: a slave point's gap/normal come from the exact master (rigid analytic or Newton projection onto the quadratic face), weights from the reference slave surface. The collocation rule matters: full Gauss on quadratic edges oscillates (over-constrained), nodal is smooth against rigid masters in 2D and 3D, reduced Gauss is right for deformable masters (nodal fails the patch test there). Mean pressure, transmitted force and contact extent are reliable; pointwise Gauss-point pressures on non-matching deformable meshes scatter.

• Friction: keep eps_n, eps_t near E / h; the symmetrised sliding tangent is indefinite, so contact analyses factor with LDL^T fallback and freeze the sliding friction in the tangent. The slip history must be seeded for points touching at the start (initial_state) and committed for points within the activation gap, not only active ones, or the first step transmits no friction.

• Contact needs load control; rigid masters move by shift + lambda * travel; the matrix pattern is extended once with ContactSet::pattern_pairs from the initial proximity (margin): a contact reaching a node pair outside it is an Err.

• Free rigid translations (RigidMaster::with_free_translation): the translation components are extra unknowns after the mesh dofs; ContactSet::eval fills ExtraOut { rq, kuq, kqq } and Ctx::newton_step condenses them with a Schur complement over the existing factorization of K_uu. Residual sign: r_y = -(force on slave), force on master = sum r_y, R_q = sum r_q - lambda load with r_q = -r_y. A free body with no contact is singular (error), so start it touching (shift). The first tangent of a step uses the near-contact stabiliser and is deliberately not consistent; every later one is exact.

• Never take a central difference of max(0, .) (found cross-checking the lug): the contact tangent differentiates the smooth active branch (PointCtx::smooth). Use NL_FDCHECK=1 (NL_FDH, NL_TRACE, NL_DEBUG) to compare any analysis' tangent with finite differences; expect ~1e-11 from the second iteration on.

• Loads::block_body loads one block only (a pin inside an assembly of bodies); Loads::body loads all blocks.

• Editing hazard (cost one lost test module): never replace a file tail with s[:i] + new when the original had content after i; check git diff --stat / test counts after scripted edits.

• Contact additions (lug integration): ContactSpec::with_overlap (an interference fit: gap = projected gap - overlap, deformable masters); solve_nonlinear_contact_from + Start (continue from a converged state, e.g. the fit, with a shorter displacement vector for bodies added later); Loads::ground (weak springs (dof, k) that hold a body only contact holds; leak k u); NlOptions::{stop_on_plateau, stop_on_fall, stall_tol, first_step}. stop_on_plateau reads travel not steps (load_flat: gained < 1.2 % of the top over the last 5 % of the travel; also a failed step at < 4 % gain is a limit state), because tiny post-cut steps made any load look flat.

• Friction tangent: a smooth-branch evaluation can see p < 0 at a point on the active-set boundary: the capacity is max(mu p, 0) (a negative one flipped the traction and 0/0 made the tangent NaN: frictional runs near full load died with ZeroPivot). factor_ldlt retries an exactly zero pivot with a 1e-10..1e-6 diagonal shift and refuses a non-finite tangent. Adding the mu dp coupling to the symmetrised sliding tangent stalls Newton, and a load-ratio predictor made friction worse (both tried, reverted).

• Speed facts (10k-dof 2D lug): faer is sequential below 40k free dofs (parallelism_for; its threads cost more than they save), the assembly keeps rayon; the 2D plastic tangent uses accumulate_tangent_2d (5x fewer operations); AMD beats geometric nested dissection ~2x in 2D; chord Newton did not help once evaluation was cheap; NL_PROFILE=1 prints evaluation / factorisation / solve time.

• Kernel upgrades of the eccentric-bushing effort (record: docs/fea-core.md Phase 9, evidence in docs/eccentric-bushing.md section 9):

  • Arc-length works with deformable contact (not with free rigid masters): after every converged arc step the AL passes run at the fixed load, every step starts from a load-controlled Newton polish, and a step that overshoots lambda_max is replaced by a load-controlled step that lands on it exactly. Start::lambda continues a run that already applied part of its load; NlOptions::stop_at_force presses a body in by prescribed displacement until the first master carries a target force (an overshooting step is halved).
  • The frictional tangent is made consistent by defect correction: the factorised matrix stays symmetric with sliding friction frozen; the exact tangent minus it (the mu dp coupling and the unsymmetric part, per frictional point, ContactStats::defect; likewise the follower-load stiffness, Eval::load_defect) is applied by refine_consistent (iterative refinement over the existing factorisation, <= 6 back-solves). Restricting the defect to sliding points was tried: Newton stalls again. Not available with free rigid translations / rotation (the Schur complement assumes symmetry: near the friction limit a moment-controlled run slows).
  • Follower pressure (Loads::followers): the load follows the deformed surface; verified against the exact Hencky compression (to 1e-9, 2D and 3D) and finite differences of the load stiffness. LDL^T is used for any run with followers (the symmetric part of the load stiffness is indefinite well before the body is).
  • Anisotropic elasticity with J2 plasticity (small strain; finite strain stays isotropic) solves the return mapping for stress and multiplier together (the flow direction changes during the return because D n is not parallel to n); with an isotropic stiffness it reproduces the isotropic law to 1e-9. Thermal expansion tensor for anisotropic materials; plasticity with a temperature change is still refused.
  • Convergence is measured against the elastic internal force (Eval::f_elem) and an absolute floor 1e-10 |f_ext|; a run whose bodies are held by nothing is refused (check_not_floating); NlSolution::ground_leak is the net force the weak grounding springs carry; a node pair outside the contact matrix pattern retries with doubled margins (4 times).
  • The closest-point projection's convergence tolerance is round-off aware (eps |x| / |tau|); a fixed 1e-14 made contact points silently leave contact once coordinates were of order one (a hole at (1, 1) failed, at the origin it solved): a_contact_solution_does_not_depend_on_where_the_model_sits.
  • The tangent predictor of prescribed displacements is evaluated without the near-contact stabiliser (ContactSet::set_stabilise), which would drag an approaching body along.
  • Contact robustness found by random corpora (2026-10-08, fea-problem/tests/sweep.rs, lug-solver/tests/kernel_sweep.rs, fea-core/tests/contact_sweep.rs; each failure was traced before it was fixed, and each fix was tried on the other corpora):
    1. Master-face choice is sticky. A slave point exactly on an element edge of the master surface (matching meshes, nodal collocation: every vertex of a concentric fit) lies on two faces whose normals differ (the surface is only C0). The minimum-gap rule picked either on round-off, so the residual jumped ~5 lbf at those points between iterations and Newton stagnated at every step size (the active-set hold cannot help: the set was constant). A point now stays on its committed face while its projection is within STICKY_TOL of that face's domain, and hold_active_set also freezes the master face of every point in contact (a point crossing an edge for real, in a stagnating solve, flipped faces again: tests/data/face_hold_eccentric.json). Diagnosis recipe: NL_TRACE prints per iteration the residual, |dx|, min gap, per contact the points in contact / slipping / off their committed face, and the two largest residual dofs; equal and opposite residuals at interface nodes that do not move with dx, plus "off their face" flickering, is this. Regressions: a_slave_point_on_a_master_element_edge_does_not_flip_faces, every_corpus_case_solves (--ignored).
    2. The state a load step returns is an equilibrated pair. refine_contact ends the augmented-Lagrangian passes with a multiplier update that is not re-solved; the state then paired the displacements with the updated multipliers, so a force read off it (NlSolution::contact.master_force, support reactions) missed the pin load by the update (4.5 % of a small pin load next to a large fit pressure; Verification::force_balance failed 8 of 24 random lugs). The normal multiplier is now the one the displacements were solved for; the friction multiplier is the friction of the iterate (committing the slip history zeroes the slip, so a re-evaluation must reproduce it). All 24 pass.
    3. One tangential penalty ratio (fit::EPS_T_RATIO = 0.03, Tuning::eps_t): stacked 3D blocks sheared near the friction limit stagnated at 0.1, converged at 0.03; 0.01 and 0.003 fail other cases. The two hand-written 0.03 in lug-solver and fit.rs now use the constant.
    4. Tried and rejected: holding the active set in every stagnating contact solve (not only frictional ones) breaks a_lightly_loaded_stiff_penalty_contact_converges_in_few_factorisations (a free pin that has not touched yet is held out of contact: singular rigid-body dofs); a stricter convergence test on the free body's own force scale breaks the_kernel_solves_every_load_angle_symmetrically (singular first step); history of a point that crossed to another master face evaluated at the current position of its committed material point (dropped: the tangent ignores that point's motion and Newton diverged on 2 of 14 stacked-block cases that pass without it).
    5. 3D partial slip: fixed (was: 1/3 of random stacked-block cases above 0.4 of the friction limit did not converge). Root cause: refine_consistent solved (K + D) x = b by plain Richardson refinement over the symmetric factorisation, which diverges when most points are at their friction limit (the defect D is then as large as K) and left the inconsistent symmetric step. It is now right-preconditioned GMRES (30 vectors), kept only if the true residual falls. 0 of 24 failures, the sweep is 6x faster; contact_sweep gates 3D to 0.8. Fixing it exposed a latent wrong answer: a hold taken at a stagnating iterate that had lost the pin's contact (see hold_plausible in nonlinear.rs).
  • Debug aids: NL_TOPRES=1 prints the four largest initial-residual dofs (node, component, position, displacement) of every step; NL_PROFILE=1 splits evaluation into element assembly and contact.

• Stress recovery (gauss_stresses_field, recover_spr, zz_error_field) runs its per-element / per-node stages with rayon above 1000 elements / 2000 nodes. Each item is computed independently and the sums are taken in element order, so the result is bit-identical for any thread count (tests/recovery.rs::recovery_is_bit_identical_on_one_thread_and_many); keep it that way (no parallel sum()/reduce of floats). tests/bench.rs::stress_recovery_throughput prints timings plus an SPR checksum to compare an optimisation against (Quad9 100x100: 66 -> 18 ms, Hex20 14^3: 101 -> 30 ms, Tet10 10^3: 29 -> 11 ms for gauss + SPR + ZZ on 8 cores). Mesh::boundary_faces sorts (sorted-corner-key, owner) pairs instead of hashing (it was 7 of 12 ms of SPR on a Tet10 mesh); the test the_sorted_face_search_equals_counting_every_face keeps the old counting method as the reference.

• Epic #25 lessons (details and numbers: docs/hybrid-meshing-amr-evidence.md):

  • Simplex consistent mass: the stiffness rule of Tri3 / Tet4 (one point) and Tri6 / Tet10 (degree 2) is too poor for int rho N_a N_b; kernel::mass uses ElementKind::mass_table (collapsed Gauss, exact to the product degree). Before, a Tri3 element mass was rank one and Tri6 frequencies carried a quadrature error larger than the discretisation error.
  • Error integrals are not stiffness integrals: the Gauss points of the stiffness rule are (nearly) the superconvergent points of the finite-element stress, so sigma* - sigma_h sampled there misses the error between them. zz_error_field and stress_error_energy integrate with recover::error_rule; the ZZ estimate then equals the true energy error on a smooth field (tests/recovery.rs).
  • SPR across materials: recover_spr patches span all elements around a node; for several materials use recover_spr_by_material (what ZZ and DWR call). A layered bar with an exact solution gives an estimate of zero only that way.
  • Which goal-oriented estimate: the recovered-stress DWR is dominated by singular corners and does not fall with the error on adapted meshes; a residual DWR with a polynomial-fit dual had the wrong sign (removed); the h-enriched dual (adaptivity::enrich, 2D triangles and quadrilaterals) is the controller's estimate (effectivity 0.7-1.0). None of them sees boundary-approximation error (a circle by quadratic edges).
  • Controller: refinement is size-field remeshing (NodalSize::Min so one marked element is honoured at its nodes); without max_growth an unreachable physics tolerance makes every element mandatory and the mesh grows fourfold per cycle.
  • Nonlinear: NlSolution::prev_state is the state at the start of the last accepted step (the consistent tangent of a return map needs it); plastic_update_at is the single constitutive path of the element kernel and of every stress read-out. Plastic history interpolation between meshes is approximate and, for monotone proportional loading, worse than no history: use warm_start_displacement. The nonlinear loop judges nonlinear and contact quantities by their change between meshes, not by an error estimate.
  • Hybrid mesher, outlines: the cells a non-axis-aligned boundary crosses (or comes within BOUNDARY_PAD sizes of) are unstructured and trace_regions follows the true boundary pieces cut at the grid lines (nodes keyed by lattice index or rounded position; a failure to close is an explicit fallback). The grid subdivision must count only structured and boundary cells of a strip: the smallest size over the whole other axis sliced regular blocks into slivers (the template plate with a hole fell back at size jump 8.1). balanced_divisions hands cells that would stretch past 3x to the unstructured mesher; remove_pinches gives up the smaller cell of a checkerboard corner (the unstructured loop would touch itself).
  • Hybrid mesher, multi-material: the unstructured mesher splits an encroached boundary segment where the interior asks for it, so two regions do not necessarily place the same interface nodes; split_shared cuts the shared segments at the union of both sides' corner nodes and meshes again.
  • Hybrid mesher: isotropic structured blocks are not cheaper than a graded Delaunay mesh; only an anisotropic block on a slender member with a local feature is (38 % fewer unknowns at equal energy error). Selection is by predicted cost (meshing::cost), never by default.

Boundaries

Always

• Add or extend a closed-form / convergence-rate test for every numerical change (tests/linear_static.rs); speed claims need a tests/bench.rs measurement.

Never

• Add UI, material tables, or file I/O here (the .vtu writer must return a String).
• Weaken a test tolerance to make a failing case pass; find the cause.

Navigation

Plan and measured numbers: docs/fea-core.md. Sibling solvers: lug-solver/AGENTS.md, edge-check/AGENTS.md.

• Native-meshing foundations: meshing::{tet_topology,tet_delaunay,tet_flip,tet_smooth,metric,metric_refine,curved,compatibility}. tetrahedralize meshes a point-set hull; Delaunay adds checked incremental insertion; PLC recovery lives in solid-kernel/src/plc_conform.rs, not here. tet_flip (exact 2-3 flips / edge removal, optional caller quality via improve_by, boundary and protected faces locked unless the opt-in planar_boundary flips edges of exactly coplanar convex boundary quadrilaterals: Phase 86, face names are not tracked) never moves vertices. transition::{tetrahedralize,tetrahedralize_wedges} convert selected Wedge6 (and Hex8) into Tet4 with a global-node diagonal rule (interface_diagonal; conforming interface; refuses quad-separating selections, first order only; no Pyramid5: converted regions only, an unconverted hex next to tets needs a pyramid element or the node tie of weld::merge, see Phase 94). check_conformity cannot see a crack between blocks that triangulate a shared quad differently: count faces against the exterior. adaptivity::MeshStore keeps the converged mesh per WorkloadKey as the starting-mesh recommendation (0.56x time, Phase 98); adaptivity::feedback summarises finished AMR runs (regional demand) and advises initial size multipliers per exact WorkloadKey from converged runs only (never coarsening; run_amr untouched). learned_cost is the optional, per-machine learned stage-time model (log-linear ridge, 90 % leave-one-out interval, explicit Fallback(reason) to cost; Phase 88); cost prices a mesh (CAD/generation/adaptation/transfer once per mesh, solver terms per analysis, memory/time limits, break_even_analyses). Its rates are single-mesh calibrations (off by 0.2-7x across sizes: calibrate near the target size) and model the direct solve only: MeshCost::iterative flags meshes Auto solves with PCG+AMG, which choose_within refuses to check against limits. motion::move_mesh / move_mesh_sliding (axis-aligned Slides) relocate nodes (any element mix) by size-stiffened pseudo-elastic sub-steps with inversion/quality gates (Gauss-point Jacobians included: a warped Hex8 can be negative between corners), rollback and a NeedsRemesh hand-over; adaptivity::moving::advance loops motion, the caller's remesher and NodalInterpolation of a nodal field; it is not ALE; ale::transport_step is the conservative upwind frame on a moving mesh (swept-volume fluxes, GCL exact; Phase 92). tet_collapse removes interior vertices by edge collapse (exact orientation, no duplicate cells, length and quality guards, returns the old-to-new map; boundary and protected vertices stay). tet_smooth moves only interior unprotected vertices (orientation-exact, monotone worst radius ratio; optimize alternates it with tet_flip); slide_planar_boundary (off by default) also slides boundary vertices inside an exactly planar patch; edges, corners and curved boundaries stay locked. metric is SPD fields (validation, gradation, intersection, complexity), metric_refine splits by Rivara terminal edges then flips in the metric then smooths interior and planar-boundary vertices in the metric (refine_projected projects inserted boundary vertices onto a curved surface with an orientation guard: Phase 90; strong anisotropy needs boundary_first_passes with FlipOptions::planar_boundary: Phase 89 converges 6:1 to 40:1, the defaults do not above 6:1) (size-penalised so grading survives; RefineOptions::smooth_passes); then collapses short edges (coarsen_passes, never below the current worst cell or the length target); curved boundaries stay polyhedral. curved::simplex_validity proves Tet10/Tri6 Jacobian positivity by Bernstein bound; tetmesh::tet10_ok requires it (snap_to_features must keep every violated face within reach active, not only those within the edge reach: Phase 85; some vertices are still refused by the quality guard on a thick cylinder, tests/high_order_convergence.rs), because its 11 sample points alone accept ~0.5 % of folded elements. compatibility registry v2 gates NativeConstrained on the shell's smallest angle (NATIVE_MIN_SHELL_ANGLE_DEG); none of this is selected by production strategy (OCCT stays production). Exact-sign degeneracies and quality gates must remain intact. Evidence: docs/issue-46-phase-1.md, -58..-79.