diff --git a/CLAUDE.md b/CLAUDE.md index 98cbd6a..1983eaa 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -10,10 +10,10 @@ knowledge base — see § Quick links spec. | Field | Value | |---|---| | Phase | 1 (Etch ↔ ECS) | -| Current milestone | M1.1.4 — Forge 3D narrowphase fast paths — code-complete, PR #57 open | -| Last released tag | `v0.11.4-narrowphase-fast-paths` | -| Active branch | `phase-1/forge/narrowphase-fast-paths` (PR #57 open, not merged) | -| Next planned milestone | M1.1.5 — Integration: semi-implicit Euler + gravity + damping — sixth core sub-milestone of the M1.1 rigid arc (M1.1.0–15, `PhysicsModule` freeze at M1.1.15). M1.1.0 (foundations), M1.1.1 (broadphase BVH), M1.1.2 (narrowphase GJK), M1.1.3 (narrowphase EPA + manifold), and M1.1.4 (narrowphase analytic fast paths ss/sb/cc/bb) are CLOSED. The M1.0 core-language series (M1.0.0–M1.0.18) is CLOSED — C1.6 tag `v0.10.18-extension-additive-warning`. Reserved (Tier-1-dependent, NOT core gaps): `future`, `override`, `quantize`. | +| Current milestone | M1.1.5 — Forge 3D integration (semi-implicit Euler) — code-complete, PR #58 open | +| Last released tag | `v0.11.5-integration-euler` (posted by Guy after merge) | +| Active branch | `phase-1/forge/integration-euler` (PR #58 open, not merged) | +| Next planned milestone | The standalone hotfix `m1.1.3-hf-epa-frame-dependence` (scheduled BEFORE M1.1.6 — its first real consumer; see Open decisions), then M1.1.6 — contact solver (Sequential Impulses, warm-start consuming the M1.1.3/4 `feature_id`s), seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15, `PhysicsModule` freeze at M1.1.15). M1.1.0 (foundations), M1.1.1 (broadphase BVH), M1.1.2 (narrowphase GJK), M1.1.3 (narrowphase EPA + manifold), M1.1.4 (narrowphase analytic fast paths ss/sb/cc/bb), and M1.1.5 (integration) are CLOSED. The M1.0 core-language series (M1.0.0–M1.0.18) is CLOSED — C1.6 tag `v0.10.18-extension-additive-warning`. Reserved (Tier-1-dependent, NOT core gaps): `future`, `override`, `quantize`. | ## Tags @@ -61,6 +61,7 @@ knowledge base — see § Quick links spec. | `v0.11.2-narrowphase-gjk` | 2026-07-18 | M1.1.2 — Forge 3D narrowphase: GJK convex intersection | Third M1.1 sub-milestone. `pipeline/narrowphase.zig` (imports `foundation` math ONLY; scalar comptime `T`, instantiated at `config.Real`): distance-based GJK on convex **cores** + inflation radius (Jolt convex-radius architecture). `SupportShape(T)` — `Core` union (point / segment:half_height / box:half_extents) + `radius` — with fixed-tie-break `support` (`dir.y == 0` → +Y endpoint; box component `== 0` → +half_extent). `RelativePose(T)` — frame-of-A precompute `rot_rel = conj(rot_a)·rot_b`, `pos_rel = conj(rot_a)·(pos_b − pos_a)` (conjugate, never `inverse()`) + `supportB`. `Simplex(T)` — triplet `Vertex (w, support_a, support_b)` + Ericson RTCD §5.1 Voronoi closest-origin solver over point/segment/triangle/tetrahedron (fixed region order; every division guarded → NaN-safe on duplicated/collinear/coplanar degeneracies). `gjk(T, shape_a, pose_a, shape_b, pose_b) GjkResult(T)` — bounded descent (`max_gjk_iterations = 32`, named relative progress tolerance, NEVER normalizes the search direction — squared distances, one `sqrt`), three regimes, classification hardened over six post-delivery correction cycles — NO classification decision depends on a geometric quantity (radius, shape size, `w`-vertex magnitude); every threshold is `k·floatEps(T)·coordScale` (coordScale = absolute support magnitude), absorbing only float noise: `separated` (`dist − r_sum` beyond `conv_k·floatEps·coordScale`, an ABSOLUTE additive margin, not `r_sum`-proportional — P1; `conv_k=16` bounds the pipeline's ACCUMULATED rounding so exact tangencies stay `.shallow` — P1b; `coord_scale = |Δpos| + coreExtent_a + coreExtent_b` is SYMMETRIC, so `.deep`/`.shallow`/`.separated` is invariant under an A/B swap — P1c) / `shallow` (else, cores disjoint; exact touch counts) both carry core distance + world closest points; `deep` (geometric ENCLOSURE — a non-degenerate tetra `count==4`, NO distance threshold on that path; plus a machine-epsilon numerical-NOISE floor `mach_eps=noise_k·floatEps` for degenerate coincident/collinear Minkowski configs scaled by the absolute support magnitude — P2, and an anti-cycling duplicate-support guard) carries the terminal simplex (EPA seed M1.1.3), no depth/normal. forge_3d integration: `shape.supportShape`, `BodyManager.rotation(id)` getter + `gjkPair(store,a,b)` (broadphase→narrowphase adapter, null on stale). Broadphase fat-AABB false positives classified `.separated`. `narrowphase.zig` 726 lines (brief-mandated single-file GJK stack, cf. `broadphase.zig` 686). Full suite green at f32 AND `-Dphysics_f64=true` (local). Out (later): EPA/manifold/depth/normal (M1.1.3), fast paths (M1.1.4). | | `v0.11.3-narrowphase-epa-manifold` | 2026-07-19 | M1.1.3 — Forge 3D narrowphase: EPA + contact manifold | Fourth M1.1 sub-milestone. narrowphase.zig → `pipeline/narrowphase/` package (pure relocation). EPA (`epa.zig`) on the cores in A's frame from the `.deep` simplex → `EpaResult{normal world A→B, core depth, closest_a/b}`. Manifold (`manifold.zig`): `collide()`/`collideOrdered()` → one generator both regimes (shallow from GJK closest points, deep from EPA); face-face by A-frame supporting-face clipping, segment by clipSegment, edge/vertex/point-core by a single witness contact along the EPA axis; per-point `penetration = r_sum − s`. `feature_id`: five class-tagged producers in disjoint ranges, each encoding its real sub-feature — kept-vertex (face,vertex), edge×plane (side-plane,edge), reference corner (ref-vertex,incident-face), clipSegment (face,capsule-endpoint), point-core/single-contact (real vertex-or-face,real vertex-or-face) — unique per simultaneous contact (131,751-manifold sweep + face_face_min + capsule-pair sweeps) and frame-stable via `collidePair`'s BodyId order; bare `collide` is pose-canonical (no inter-frame guarantee). `ContactManifold`/`ContactPoint` FROZEN. Order-independence via `collidePair` BodyId order + normal negation (RD-4 exception: bit-identical coincident shapes). Integration: `BodyManager.collidePair`, `Real` re-exports. feature_id contract hardened across an extended Codex review cycle (FIX-1..12). Note: P1d extreme-aspect box `.deep` = GJK classification limit upstream of EPA → M1.1.4. Fast paths M1.1.4. | | `v0.11.4-narrowphase-fast-paths` | 2026-07-20 | M1.1.4 — Forge 3D narrowphase fast paths | Fifth M1.1 sub-milestone. Four analytic per-pair fast paths dispatched from `collideOrdered` before GJK — sphere/sphere, sphere/box, capsule/capsule, box/box (radius-0 box only; rounded box → generic). Each kernel computes only a `ContactSeed {normal A→B, core closest points, base_penetration}` fed to the UNCHANGED `generateManifold`, so feature_ids/clipping/≤4-reduction/order-independence/frame-stability are inherited by construction; the fast path replaces only the GJK descent + EPA expansion. `fastSeed` 3-state; `separated` mirrors gjk.zig (conv_k=16, coord_scale=|Δpos|+coreExtent). sphere/box deep = closed-form least-penetration face (point/box P1d fix, >200:1). box/box = 15-axis SAT, order-independent; separation on every non-zero raw edge×edge axis; edge-edge witness via full-degenerate Ericson closestSegSeg; fixes box/box P1d and bypasses the M1.1.3 deep-rotated EPA frame-dependence for box/box. capsule/capsule = closestSegSeg → 3 regimes via generateManifold (crossed via E1 fix; collinear/parallel radial, never axial, fallback normal). E1 fix: segment×segment non-parallel → single witness. Hardened through a 5-round external review to a threshold-class discipline, all RED-first in-milestone: normalization guards fire ONLY at true zero (`floatMin`, never a geometric scale); point-set dedup is relative to the candidates' PER-AXIS EXTENT (translation+scale+anisotropy invariant); separation margins are `k·floatEps·coordScale`. Defects closed: degenerate segment closest-points, SAT separating-axis completeness, collinear-capsule normal, and the length-threshold class (fast_paths + shared manifold.zig generator). Validation: differential vs `collideOrderedGeneric` (gated on oracle self-consistency), closed-form P1d/extreme-aspect, oracle-free deep+rotated box/box, feature_id producer×pair×order matrix, scale×translation×anisotropy×decentre×order invariance matrix. Bench (ReleaseFast, dispatched/generic): 0.55 ss / 0.27 sb / 0.27 bb / 0.12 cc. Full suite green f32 + `-Dphysics_f64=true`, debug + ReleaseSafe. No gjk.zig/epa.zig touched. | +| `v0.11.5-integration-euler` | 2026-07-23 | M1.1.5 — Forge 3D integration (semi-implicit Euler) | Sixth M1.1 sub-milestone — the FIRST where bodies move. `pipeline/integration.zig`: `integrate(bm, dt, gravity)`, one pure per-tick pass in ascending slot-index order over the live SoA store (`IdAllocator.isAliveIndex` filters the non-compacting store's dead slots), FREE-FLIGHT only (broadphase/narrowphase exist but are NOT invoked; contact response is M1.1.6). Semi-implicit (symplectic) Euler: `v` first, `x` from the NEW `v`. Gravity = ACCELERATION × `gravity_factor` (mass-independent), never a force. Damping = Jolt clamped-linear `v *= max(0, 1 − d·dt)` (the clamp is physical, not a geometric epsilon); applied exactly once per call — `integrate` has no substep opinion (orchestrator's concern, M1.1.15; the form is substep-count-sensitive, `(1−d·dt/N)^N → e^{−d·dt}`). Angular: `α = (R·I_local_inv·Rᵀ)·τ`, `ω` clamp-damped, orientation FIRST-ORDER `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with world-space ω on the LEFT — no `|ω|` division ⇒ singularity-free at ω=0 (M1.1.4 threshold discipline: zero guards added; `ω_quat⊗q ⊥ q` in 4D so `‖q+dq‖ ≥ 1`, normalize never near-zero). Gyroscopic term DROPPED (Jolt/PhysX/Bevy default). `Body` gains `force`/`torque` world-space per-tick accumulator columns (layout-clean, position/rotation bulk-sync invariant untouched); `BodyManager` gains stale-safe `linearVelocity`/`angularVelocity` getters, `setLinearVelocity`/`setAngularVelocity`, `addForce`/`addTorque` (accumulate), `addImpulse` (immediate `Δv = impulse·inv_mass`, naturally no-op on static/kinematic); `addTorque`/`setAngularVelocity` INTERNAL (the day-1 `PhysicsModule` 3D surface has no angular mutators — an M1.1.15 freeze decision). §2 accumulator reset each tick is UNIFORM over LIVE bodies only; static/kinematic cleared but never moved. Discrete free-fall oracle `x_n = x₀ + n·v₀·dt + g·dt²·n(n+1)/2` (−4.98675 m at 60 Hz/1 s, NOT the continuous −4.905). Codex review (2 P2, test-only): dead-slot skip now OBSERVED bit-exact (stale pos/vel unchanged AND stale force NOT cleared, checked before slot reuse), and the LEFT quaternion product LOCKED by a one-step non-commutative oracle (left vs right explicitly discriminated, discrimination guard protects the test's power). 112/112 green f32 + `-Dphysics_f64=true`, debug + ReleaseSafe. Production untouched by the review fixes. | ### Hotfixes (untagged) @@ -118,8 +119,9 @@ Hotfix milestones are merged to `main` without a tag (Guy decision, - **M1.1.2 scope boundary (narrowphase GJK convex detection)**: `pipeline/narrowphase.zig` imports `foundation` (math) ONLY — never `weld_forge`/`body*`/`config`/`broadphase`; the scalar is the comptime `T` (forge_3d instantiates at `config.Real` via `root.zig` re-exports). GJK runs on the convex **cores** (point/segment/box) + an inflation radius (Jolt convex-radius architecture); touch = `dist(cores) ≤ r_a + r_b`, `separated` iff the core distance exceeds `r_a + r_b` beyond an ABSOLUTE noise-scaled contact margin `conv_k·floatEps(T)·coordScale` (exact touch → `shallow`, face-inclusive convention). Distance-based (NOT boolean) GJK: the simplex stores `(w, support_a, support_b)` triplets so closest points on A/B reconstruct from the barycentrics; the Voronoi simplex solver is Ericson (not Johnson, not signed-volumes). Determinism by construction (anticipates M1.1.14): no hash containers, no trig (dot/cross only), `max_gjk_iterations = 32`, named relative progress tolerance, fixed support tie-breaks; the search direction is NEVER normalized (squared distances, one `sqrt`); computation frozen in the frame of A (`rot_rel`/`pos_rel` precompute, conjugate-rotate never `inverse()`) — changing the frame after M1.1.14 would break bit-exactness. `GjkResult.deep` carries the terminal origin-enclosing simplex as the EPA seed — no depth/normal/manifold here. The forge_3d-side adapters (`shape.supportShape`, `BodyManager.gjkPair`) own the `Shape → SupportShape` conversion + the `BodyId`-level entry, mirroring the broadphase `user_data` discipline. **Classification robustness (six post-delivery correction cycles, Guy-directed — RESOLVES the E3 absolute-`deep_eps`, the Fix 3 vertex-relative residual, the P1/P2 remaining geometric dependencies, the P1b accumulated-rounding under-dimensioning, the P1c A/B-order asymmetry, AND the P1d anisotropic-tetra false-degeneracy; the extreme-aspect radius-0 box `.deep` residual is DOCUMENTED and deferred to M1.1.4/M1.1.3):** GUIDING PRINCIPLE — no `.deep`/`.shallow`/`.separated` decision depends on a GEOMETRIC quantity (radius, shape size, `w`-vertex magnitude); `.deep` = geometric enclosure; every remaining threshold is `k·floatEps(T)·coordScale` (coordScale = absolute support magnitude, `@sqrt(maxSupportMagSq)`), absorbing only float rounding noise. (a) `.deep` = a non-degenerate tetra enclosing the origin (`count==4`, degeneracy tested DIMENSIONLESSLY `|bary_det|² ≤ deg_rel²·(|b−a|²·|c−a|²·|d−a|²)` = squared normalized volume, aspect-ratio-invariant — P1d replaced the `maxEdgeSq³` normalization that rejected valid anisotropic tetra and read a sharp box's interior `.separated`), NO distance threshold on that path; (a2, Fix 3b→P2) the sole distance early-out (`degenerateOriginReached`) is a numerical-NOISE floor `closest·closest ≤ mach_eps²·maxSupportMagSq` (scaled by the ABSOLUTE support magnitude), with `mach_eps = noise_k·floatEps(T)` (P2 recalibrated from the 8×-eps `1e-6` literal — which was a visible geometric threshold, ~8.7e-5 at half-extent 50 — down to `noise_k`≈2 ULPs); the earlier vertex-relative `rel_deep²·maxVertexMag²` (Fix 3) reproduced the bug for large shapes and was removed; (b) an anti-cycling duplicate-support guard kills degenerate Minkowski tetrahedra at the source; (c, Fix 4→P1) `.separated iff dist − r_sum > conv_k·floatEps(T)·coordScale` — an ABSOLUTE additive margin (P1 removed the `r_sum·(1+contact_rel)` form, whose `contact_rel·r_sum` slack — 0.1 m at r_sum 1000 — was a collision margin beyond the core radius, out of scope; it mis-classified two r=500 spheres 5 cm apart as `.shallow`), honoring the frozen touch=shallow rule via the convergence-noise margin; (c2, P1b) `conv_k` = 16, not 2 — 2 ULP under-dimensioned the ACCUMULATED pipeline rounding (quaternion + Voronoi tetra + sqrt), so an exact tangency (Codex case: `dist − r_sum` = 7.15e-7 ≈ 2.1 ULP relative, zero convergence residue) tipped `.separated`; 16 ULP bounds it and stays at the noise level (the deep noise floor keeps `noise_k = 2`, a distinct point-noise floor — aligning them recreates a P1b/P2 false positive). Scale-robust, verified at O(1), half-extents 50–500, and at the f32 noise floor (~9 orders tighter in f64 — an f32 fundamental limit, not a geometric margin); a measure-zero EXACT tangency can still, in rare pathological configs (~1/20000 in sweep), tip `.separated` by one ULP — irreducible float limit, corrected next frame by penetration, no absolute guarantee claimed. (c3, P1c) the contact margin's `coord_scale` is SYMMETRIC by construction — `|pos_b − pos_a| + coreExtent(a) + coreExtent(b)` (`coreExtent` = the core's max local support magnitude, radius excluded: point→0, segment→half_height, box→`|half_extents|`) — replacing the A-frame terminal-support magnitude `@sqrt(maxSupportMagSq)`, which was frame-dependent and made one tangency classify `.separated` in one A/B order and `.shallow` in the other (a structural defect, not calibration). `.deep`/`.shallow`/`.separated` is therefore invariant under an A/B swap by construction, asserted by `gjk classification is order-independent`; `maxSupportMagSq` stays (anti-cycling guard + deep noise floor — `w`-magnitude, a distinct quantity). **File-length note:** `narrowphase.zig` is 726 lines (over the 500 Review guideline) — the brief mandates the whole GJK stack (support shapes + simplex solver + loop + result) in this single file (mirror of `broadphase.zig` at 686); conscious brief-driven overage. **OUT (not gaps; additive/later):** EPA + penetration depth + contact normal + manifold (M1.1.3 — in particular do NOT compute `depth = r_a+r_b−dist` or a normal for `.shallow`); pair-type fast paths ss/sb/cc/bb (M1.1.4); speculative-contact thresholds / margins beyond the core radius / box convex radius (box radius 0 here); `collision_layer`/`CollisionConfig` matrix wiring; stepping/`PhysicsWorld`/`PhysicsModule` instantiation (M1.1.15); raycast/shapecast/point queries (M1.1.9–10); shapes beyond sphere/box/capsule (`SupportShape.Core` stays additively extensible); batched/SIMD pair processing; `forge_2d` (M1.8); robust `.deep` for radius-0 box cores of EXTREME aspect ratio (>~50:1) — a GJK-f32 limit on sharp cores (tetra degeneracy at extreme anisotropy + premature anti-cycling termination before the enclosure tetra forms), reliable to ~30:1, deferred to M1.1.4 analytic box fast paths / M1.1.3 EPA (P1d). - **Extension-conflict semantics (last-wins vs reject) — RESOLVED (M1.1.1-HF4)**: the dedicated post-HF3 spec session ratified **reject** as the normative, permanent policy for additive extension conflicts, superseding the former "last-extension-wins" prose. The invariant is **{base components} ∪ {active extensions' components} conflict-free**; three rejected forms: (a) two extensions declare the same component; (b) an extension declares a component already carried by the base (or an earlier extension); (c) the same extension is listed twice. Static conflict → **fatal cook error** `E1797 ExtensionAdditiveConflict` (all three forms, distinct static message; no output; `weld check` fails in CI), guaranteeing `cooked ⇒ loadable`; dynamic activation → runtime `error.ExtensionComponentConflict` (a/b) / `error.ExtensionAlreadyActive` (c), atomic (entity unchanged). Specs corrected (`engine-scene-serialization.md` additive-conflict note §695–706, `etch-reference-part2.md §30.4–30.5`, `etch-diagnostics.md §18.3` E1797 row + the W1791→E1797 reclassification note); code comments cite `engine-scene-serialization.md` (the authority), not the HF3 brief. Deactivation needs no provenance — the reject invariant means no two active declarants ever share a component. Cook detection: form (c) is a resolver-independent lexical pass (so the public `cook` null-resolver path rejects duplicates); forms (a)/(b) need the prefab resolver and are `validate.structure`-guarded before any accessor getter. Alternative examined and rejected in the session: last-wins-with-provenance (masking + shadowed-byte side table + out-of-order deactivate + `.sav` shadow-stack serialization — ≈5–10× the code, touches the ECS reserve-then-mutate invariant and the serialization format, undoes the HF3 deactivate simplification — to support a design error). `@exclusive_with` remains a future, unimplemented escape hatch. - **Scope boundary M1.1.3** — Deferred to M1.1.4: extreme-aspect (>~50:1) radius-0 box `.deep` classification (GJK-f32 limit upstream of EPA), and analytic fast paths superseding the generic clipper for common pairs. Deferred until the shapes ship (convex hull/compound/mesh, ~M1.1.19-20, co-designed with M1.1.6 warm-start): generic-clipper feature ids for non-fast-path convex shapes. In M1.1.3 the feature_id contract is fully met for box/sphere/capsule (all producers class-tagged, unique, frame-stable via collidePair). **Closed for M1.1.3 scope.** -- **M1.1.4 — generic EPA deep-rotated frame-dependence (M1.1.3 DEFECT, tracked).** `collideOrderedGeneric` (GJK/EPA in the frame of A) is order-dependent for a DEEP, ROTATED convex pair even at 1:1 aspect (sweep config: pen 1.047 forward, pen 0/count 1 reversed) — violates the frozen order-independence contract (`engine-physics-forge.md §303`), independent of the P1d aspect limit. The box/box SAT fast path is order-independent by construction and neutralizes the exposure FOR box/box; **residual: capsule/box + sphere/capsule stay on the generic EPA path** (demo pairs). Deferred to a dedicated hotfix `m1.1.3-hf-epa-frame-dependence` (out of M1.1.4 scope — gjk.zig/epa.zig untouched). Detected + recorded at M1.1.4 E3. +- **M1.1.4 — generic EPA deep-rotated frame-dependence (M1.1.3 DEFECT, tracked).** `collideOrderedGeneric` (GJK/EPA in the frame of A) is order-dependent for a DEEP, ROTATED convex pair even at 1:1 aspect (sweep config: pen 1.047 forward, pen 0/count 1 reversed) — violates the frozen order-independence contract (`engine-physics-forge.md §303`), independent of the P1d aspect limit. The box/box SAT fast path is order-independent by construction and neutralizes the exposure FOR box/box; **residual: capsule/box + sphere/capsule stay on the generic EPA path** (demo pairs). Deferred to a dedicated hotfix `m1.1.3-hf-epa-frame-dependence` (out of M1.1.4 scope — gjk.zig/epa.zig untouched). Detected + recorded at M1.1.4 E3. SCHEDULED: before M1.1.6 (Sequential Impulses is the first real consumer of deep manifolds under sustained contact — the frame-dependence would corrupt warm-starting; M1.1.5 was free-flight and never invoked EPA). - **M1.1.4 scope boundary — fast paths.** Four frozen pairs only (ss/sb/cc/bb), radius-0 box. OUT (generic, correct): capsule/box, sphere/capsule (capsule/box a post-demo additive candidate after profiling). OUT (later): warm-start (M1.1.6), stepping/solver (M1.1.5+), Plane/Mesh (M1.1.11). `collideOrdered` count>1 point-set is fixed-order; full order-independence via `collide`/`collidePair` (frozen M1.1.3 contract, inherited). **Precision boundary (out of scope):** contact-point POSITION accuracy far from body A's centre is f32-limited (absolute noise ~floatEps·|pos|); order-independence + count are correct (extent-based dedup), but a feature below f32 resolution at its position has an inherently ambiguous count — a `generateManifold` precision characteristic (would need clip-recentring), noted for a future precision hotfix. **Closed for M1.1.4 scope.** +- **M1.1.5 scope boundary (integration, semi-implicit Euler).** Integration is INTERNAL to forge_3d (`pipeline/integration.zig`, re-exported at `Real` by `root.zig`); NOT wired into `PhysicsModule` — the `step` orchestration (broadphase → narrowphase → solve → integrate) and the ECS Transform sync land at M1.1.15. `integrate(dt)` is substep-neutral: damping applies once per call; the substep-count sensitivity of the clamped-linear form (`(1−d·dt/N)^N → e^{−d·dt}`) is the orchestrator's decision at M1.1.15. OUT (later/additive, NOT debt): contact response / Sequential Impulses / friction / restitution / warm-start (M1.1.6); sleeping (M1.1.8); CCD; kinematic position-from-velocity (additive branch, no refactor); the gyroscopic term `ω × (I·ω)` (additive implicit-gyroscopic option); exact non-linearized orientation integration (additive strategy flag — introduces an `|ω|` division singular at zero, deliberately avoided in the first-order default); ECS ↔ BodyManager sync + `PhysicsModule` instantiation (M1.1.15); SIMD / island-parallel / batched integration (M1.1.8+/Phase 4+); absolute integration perf targets / bench (land with the full `step`, per the M1.1.3/M1.1.4 precedent). The day-1 `PhysicsModule` 3D interface asymmetry (no `addTorque`/`setAngularVelocity`, unlike 2D's `addTorque2D`/`setAngularVelocity2D`) is FLAGGED for the M1.1.15 freeze review — the internal `BodyManager` mutators exist either way. ## Non-negotiable rules @@ -253,4 +255,4 @@ The `briefs/` directory is the source of truth for milestone state. The brief's --- -Last updated: 2026-07-20 +Last updated: 2026-07-23 diff --git a/briefs/M1.1.5-integration-euler.md b/briefs/M1.1.5-integration-euler.md new file mode 100644 index 0000000..f9e3b37 --- /dev/null +++ b/briefs/M1.1.5-integration-euler.md @@ -0,0 +1,293 @@ +# M1.1.5 — Forge 3D integration (semi-implicit Euler) + +> **Status:** CLOSED +> **Phase:** 1 +> **Branch:** `phase-1/forge/integration-euler` +> **Planned tag:** `v0.11.5-integration-euler` +> **Dependencies:** M1.1.0 (SoA `BodyManager` store, `MotionProperties` with analytic inverse mass/inertia + `linear_damping`/`angular_damping`/`gravity_factor`, FROZEN `RigidBody`/`Velocity`/`PhysicsForces` ECS components, `IdAllocator` generational slots). Sequence predecessor / branch base: M1.1.4 (`v0.11.4-narrowphase-fast-paths`). +> **Opened:** 2026-07-23 +> **Closed:** 2026-07-23 + +--- + +# FROZEN SECTION + +*Produced by Claude.ai. Not modifiable by Claude Code outside a Claude.ai round-trip (cf. § Recorded deviations).* + +## Context + +This is the FIRST sub-milestone where bodies move. M1.1.0–4 laid the static machinery: the SoA `BodyManager` store (position, rotation, linear/angular velocity, `MotionProperties`), the broadphase BVH, and the narrowphase (GJK/EPA/manifold + analytic fast paths). No stepping exists yet — M1.1.0 explicitly deferred it (`PhysicsModule` freezes at M1.1.15). This milestone adds the integration stage the plan freezes at the M1.1.5 row of `engine-phase-1-plan.md` — "Integration : semi-implicit Euler + gravité + damping", demonstrator "Box en chute libre, trajectoire correcte" — landing it in `pipeline/integration.zig`, the file the spec directory tree (`engine-directory-structure.md`; `engine-spec.md §3.5`) reserves for "semi-implicit Euler, sleep, CCD" (only the Euler + gravity + damping part lands here; sleep is M1.1.8, CCD is later). + +It is FREE-FLIGHT integration only. Broadphase/narrowphase exist but are NOT invoked — there is no contact response (Sequential Impulses is M1.1.6). A dynamic body under gravity follows an unobstructed trajectory. Integration lives INSIDE `forge_3d` (a free function over the `BodyManager` SoA store); it is NOT wired into the `PhysicsModule` interface, whose `step` orchestration composes broadphase → narrowphase → solve → integrate at M1.1.15. + +## Scope + +Integration is a pure per-tick function over the live SoA store, in ascending slot-index order (deterministic; no hash container on the path — the `BodyManager` discipline), written against `Real` (the comptime f32/f64 flip). Semi-implicit (symplectic) Euler: velocity is integrated first, position from the *new* velocity. Gravity is applied as an ACCELERATION scaled by `gravity_factor` (mass-independent — Galileo; the demonstrator's whole point), never as a force. Damping is the Jolt clamped-linear form `v *= max(0, 1 − d·dt)`. Orientation is integrated first-order (linearized), `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with the world-space `ω` on the LEFT — no `|ω|` division, hence no singularity at `ω = 0` (respects the M1.1.4-frozen threshold discipline: normalization guards fire only at true zero). The gyroscopic term `ω × (I·ω)` is DROPPED (Jolt/PhysX/Bevy default; explicit integration of it injects energy). + +- **E1 — Body force/torque accumulators + `BodyManager` mutators/getters + `IdAllocator.isAliveIndex`.** + - **(a) Accumulator columns.** Add `force: Vec3r` and `torque: Vec3r` (world-space, N and N·m) to the SoA `Body` row. `addBody` initializes both to `Vec3r.zero`. The position/rotation layout-clean invariant for the M1.1.15 `Transform` bulk-sync is UNAFFECTED (a `MultiArrayList` stores each field as an independent column); update the `body.zig`/`body_manager.zig` doc comments to mention the two columns. + - **(b) Slot liveness by index.** Add `pub fn isAliveIndex(self: *const IdAllocator, index: u32) bool` to `slot_alloc.zig` — `index < slots.items.len and slots.items[index].alive`. Required because `removeBody`/`free` does NOT compact: `bodies.len` is the high-water mark including dead slots, so an index-ascending pass must filter liveness (`validate` takes a packed id, not a bare index, so it cannot serve here). + - **(c) Mutators + getters on `BodyManager`** — all validate the `BodyId` and are a no-op / return `null` on a stale handle (parity with the existing `position`/`rotation` getters): + - `linearVelocity(id) ?Vec3r`, `angularVelocity(id) ?Vec3r` (getters); + - `setLinearVelocity(id, Vec3r)`, `setAngularVelocity(id, Vec3r)`; + - `addForce(id, Vec3r)`, `addTorque(id, Vec3r)` — ACCUMULATE into the `force`/`torque` columns (any live body; the per-tick reset in E2 clears them uniformly); + - `addImpulse(id, Vec3r)` — IMMEDIATE linear velocity change `Δv = impulse · inv_mass` (naturally a no-op on static/kinematic, `inv_mass == 0`; no `body_type` branch needed). + - `addTorque` and `setAngularVelocity` are INTERNAL to `BodyManager`. Do NOT extend the public `PhysicsModule` 3D interface (it is frozen only at M1.1.15, and its day-1 surface — `engine-tier-interfaces.md §1` — carries `addForce`/`addImpulse`/`setLinearVelocity` but no `addTorque`/`setAngularVelocity`; the asymmetry is a M1.1.15 freeze decision, not this milestone's). + - No integration yet. Non-regression: the full existing `test-forge-3d` suite is green (the new columns break no size/AABB/inertia pin). Tests: each mutator accumulates/sets, each getter reads back; a stale handle is a no-op / `null`; `addImpulse` gives `Δv = impulse·inv_mass` and is a no-op on a static body; `isAliveIndex` is correct across create / free / reuse. + +- **E2 — `integrate` (linear + gravity + damping + position) + the free-fall oracle.** Create `pipeline/integration.zig` with `pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void`. Single ascending pass over `0..bodies.len` filtered by `isAliveIndex` (accessing the SoA columns directly — DOD). Per live body: clear its `force`/`torque` accumulators (the `engine-physics-forge.md §2` "reset each fixed tick" invariant, holds for ALL live bodies). For DYNAMIC bodies only, before the clear: `a = gravity·gravity_factor + force·inv_mass`; `v += a·dt`; `v *= max(0, 1 − linear_damping·dt)`; `x += v·dt`. Static/kinematic: accumulators cleared, no motion (kinematic position-from-velocity is a later additive concern; no kinematic demonstrator here). Tests: **discrete semi-implicit free-fall oracle** — `linear_damping = 0`, `x_n = x₀ + n·v₀·dt + g·dt²·n(n+1)/2`, so after N=60 @ dt=1/60 with v₀=0: `x = −4.98675 m` (NOT −4.905; a fixed-step explicit integrator does not reproduce the continuous parabola), `v = −9.81 m/s` (the velocity DOES match the continuous value, `v_n = n·g·dt`); mass-independence (two masses, same trajectory); `gravity_factor` scaling (0.5 → half acceleration; 0 → no gravity contribution); damping bends the trajectory (`|v|` strictly below undamped after one step) and clamps (`d·dt > 1` → `v = 0`, no sign flip); forces consumed once (a second `integrate` with no re-apply shows no residual force); `addImpulse` immediate `Δv`; static/kinematic bodies unmoved but their accumulators cleared; a freed dynamic slot's stale data is NOT integrated; determinism (two runs bit-identical). + +- **E3 — Angular integration (torque → world inertia → ω → damping → first-order quaternion) + close.** Extend `integrate`'s dynamic branch with the angular half: `I_world_inv = R · I_local_inv · Rᵀ` (`R = Mat3r.fromQuat(rotation)`, `I_local_inv` the diagonal `MotionProperties.local_inv_inertia`); `α = I_world_inv · torque`; `ω += α·dt`; `ω *= max(0, 1 − angular_damping·dt)`; `q ← normalize(q + (½·dt)·(ω_quat ⊗ q))` where `ω_quat = Quatr{ ω.x, ω.y, ω.z, 0 }` and `⊗` is `Quatr.mul` with `ω_quat` on the LEFT (`ω_quat.mul(q)` — world-space `ω`). Gyroscopic term dropped. Tests: torque on a rotated anisotropic box → `ω` after one step equals `(R·I_local_inv·Rᵀ)·torque·dt` (closed-form for a known `R` + diagonal `I`); angular damping reduces `ω`, `d·dt > 1` clamps to 0; the orientation quaternion stays unit (`|q|` within a named tolerance of 1 after N steps); `ω = 0` leaves `q` unchanged with no NaN (the first-order path needs no zero guard); determinism. Then: (1) `CLAUDE.md §3.4` patch (content produced by Claude.ai at the E3 review, applied by CC on-branch). (2) Close. Full suite green at f32 AND `-Dphysics_f64=true`, debug + ReleaseSafe. + +## Out of scope + +- **Contact response / solver.** No Sequential Impulses, no friction, no restitution, no warm-start — M1.1.6. `integrate` does NOT call broadphase or narrowphase; a body under gravity is unobstructed. +- **Constraints / joints, sleeping (M1.1.8), CCD, sub-stepping.** `PhysicsConfig.max_substeps` exists but is not consulted here — `integrate(dt)` applies damping exactly once with the `dt` it is given and has NO opinion on substeps; call cadence is the orchestrator's concern (M1.1.15). +- **Gyroscopic term `ω × (I·ω)`.** Dropped (Jolt/PhysX/Bevy default). An implicit-gyroscopic option is a future additive flag. +- **Exact (non-linearized) orientation integration.** First-order only. `fromAxisAngle(ω/|ω|, |ω|·dt) ⊗ q` is a future additive strategy flag if a fast-spinner demonstrator ever needs it (it introduces an `|ω|` division singular at zero, requiring a true-zero guard — deliberately avoided here). +- **Kinematic position-from-velocity integration.** Additive branch, no refactor; no kinematic demonstrator in M1.1.5. +- **ECS ↔ `BodyManager` sync**, **`PhysicsModule` instantiation / `step` orchestration**, **double-precision sync** — all M1.1.15. +- **Extending the public `PhysicsModule` 3D interface.** `addTorque`/`setAngularVelocity` stay internal to `BodyManager`; the interface is frozen at M1.1.15. +- **SIMD / island-parallel / batched integration.** Scalar, single-threaded, ascending-index. Island-parallel integration is M1.1.8+; SIMD/GPU is Phase 4+. +- **Absolute performance targets / a bench.** Integration perf lands with the full `step` (later), per the M1.1.3/M1.1.4 precedent. No `build.zig` change, no bench file. +- **Spec file edits on-branch.** Per `engine-development-workflow.md §3.5` the spec lives only in the Claude.ai KB. A new `engine-physics-forge.md` "Integration" section (capturing the Euler order, gravity-as-acceleration, the damping form, the first-order quaternion, the dropped gyroscopic term) is produced by Claude.ai as a complete re-uploadable file at the E3 review; CC edits NO spec file. +- **The EPA hotfix `m1.1.3-hf-epa-frame-dependence`.** A separate hotfix milestone, scheduled before M1.1.6 (its first real consumer). NOT this milestone — M1.1.5 does not touch EPA. + +## Specs to read first + +1. `engine-physics-forge.md` — §1.2 (pipeline architecture: "Integration : semi-implicit Euler, sleep detection, CCD" — the file placement), §2 (Physics Bodies: `RigidBody`/`Velocity`/`PhysicsForces` component fields, and the `PhysicsForces` "reset chaque fixed tick" contract), §1.5–1.6 (the `Real` scalar / f32 default / `-Dphysics_f64` flag). +2. `engine-project-settings.md` — §3.2 Physique: `PhysicsConfig.gravity = [0, -9.81, 0]` (the global-gravity source the integrator's `gravity` parameter will carry at M1.1.15), `fixed_rate = 60`, `max_substeps = 4`. +3. `engine-coordinate-system.md` — §1 (right-handed, Y-up, `Vec3.down = (0,-1,0)`; SI units — m, m/s, m/s², kg, N), §3 (quaternion storage + right-hand-rule sign convention). +4. `engine-tier-interfaces.md` — §1 `PhysicsModule` (the day-1 contract the internal mutators anticipate: `step`/`addForce`/`addImpulse`/`setLinearVelocity`; note NO `addTorque`/`setAngularVelocity` in 3D). +5. `engine-phase-1-plan.md` — the M1.1.5 row + the M1.1 rigid-arc context (comptime core, `PhysicsModule` freeze at M1.1.15, M1.1.6 solver consumes contacts). +6. `engine-zig-conventions.md` — Zig 0.16.x conventions (no `@cImport` outside `*_c` modules, no `usingnamespace`, `unreachable` only on proven invariants, English doc comments on public API, naming). + +## Files to create or modify + +- `src/modules/forge/forge_3d/pipeline/integration.zig` — **create** — `pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void` (linear E2 + angular E3). Imports `foundation` (math), `config.zig`, `body_manager.zig`. Operates on the SoA columns directly, ascending index, `isAliveIndex`-filtered. +- `src/modules/forge/forge_3d/slot_alloc.zig` — **modify** — add `pub fn isAliveIndex(self: *const IdAllocator, index: u32) bool`. +- `src/modules/forge/forge_3d/body.zig` — **modify** — add `force: Vec3r` + `torque: Vec3r` columns to `Body`; update the doc comment. +- `src/modules/forge/forge_3d/body_manager.zig` — **modify** — initialize `force`/`torque` to zero in `addBody`; add `linearVelocity`/`angularVelocity` getters, `setLinearVelocity`/`setAngularVelocity`, `addForce`/`addTorque`/`addImpulse` (all validate + no-op/`null` on stale); doc-comment touch-up. +- `src/modules/forge/forge_3d/root.zig` — **modify** — pin `integration.zig` in the comptime block; pin `tests/integration_test.zig`; re-export `integrate` at `Real` (`pub fn integrate(bm, dt, gravity)` forwarding to `integration.integrate`). +- `src/modules/forge/forge_3d/tests/integration_test.zig` — **create** — the full E2 + E3 test matrix. +- `CLAUDE.md` — **modify** — §3.4 update (E3; patch content produced by Claude.ai at review). + +## Acceptance criteria + +### Tests + +Named tolerances: assert the discrete free-fall closed form to a named fp tolerance (documented at the test site; NOT `1e-6` for the f32 path — sized to accumulated fixed-step rounding over N steps). The quaternion-unit tolerance and the angular closed-form tolerance are likewise named at their sites. + +- `tests/integration_test.zig` — `test "free fall matches the discrete semi-implicit oracle"` — `linear_damping = 0`, N=60 @ dt=1/60, v₀=0: `x ≈ −4.98675 m`, `v ≈ −9.81 m/s`. +- `tests/integration_test.zig` — `test "free fall is mass-independent"` — two dynamic bodies, `mass₁ ≠ mass₂`, identical trajectory (gravity as acceleration). +- `tests/integration_test.zig` — `test "gravity_factor scales the fall"` — `gravity_factor = 0.5` → half the gravitational acceleration; `gravity_factor = 0` → no gravity contribution. +- `tests/integration_test.zig` — `test "linear damping bends the trajectory and clamps"` — `d > 0` → `|v|` strictly below the undamped value after one step; `d·dt > 1` → `v = 0`, no sign flip. +- `tests/integration_test.zig` — `test "forces are consumed once per tick"` — `addForce` then `integrate` applies it; a second `integrate` with no re-apply shows no residual force contribution (accumulator cleared). +- `tests/integration_test.zig` — `test "impulse is an immediate velocity change"` — `addImpulse` → `Δv = impulse·inv_mass` immediately (before any `integrate`); no-op on a static body. +- `tests/integration_test.zig` — `test "torque on a rotated anisotropic box"` — a rotated box with a diagonal inertia; after one step `ω = (R·I_local_inv·Rᵀ)·torque·dt` (closed-form for the known `R` + `I`). +- `tests/integration_test.zig` — `test "angular damping and clamp"` — `angular_damping > 0` reduces `ω`; `d·dt > 1` clamps `ω` to 0. +- `tests/integration_test.zig` — `test "orientation quaternion stays unit"` — after N steps with nonzero `ω`, `|q|` within the named tolerance of 1. +- `tests/integration_test.zig` — `test "zero angular velocity leaves orientation unchanged"` — `ω = 0` → `q` unchanged, no NaN. +- `tests/integration_test.zig` — `test "static and kinematic bodies are not integrated"` — position/velocity/orientation unchanged after `integrate`; their `force`/`torque` accumulators ARE cleared. +- `tests/integration_test.zig` — `test "freed slots are skipped"` — a freed dynamic body's stale data is not integrated; a body created in a reused slot integrates correctly. +- `tests/integration_test.zig` — `test "integration is deterministic"` — two identical N-step runs over several bodies, bit-identical. +- All tests green in `debug` AND `ReleaseSafe`, at f32 AND `-Dphysics_f64=true`. + +### Benchmarks + +None (out of scope — absolute integration perf lands with the full `step`, later). + +### Observable behavior + +- `zig build test-forge-3d` green (f32 + f64, debug + ReleaseSafe). + +### CI + +- `zig build` clean, zero warnings, on the configured matrix +- `zig build test` green (debug + ReleaseSafe), including the `-Dphysics_f64=true` local sweep +- `zig fmt --check` green +- `zig build lint` green +- `commit-msg` hook green on every commit of the branch + +## Conventions + +- **Branch:** `phase-1/forge/integration-euler` +- **Final tag:** `v0.11.5-integration-euler` +- **PR title:** `Phase 1 / Forge / Forge 3D integration (semi-implicit Euler)` +- **Commit convention:** Conventional Commits (cf. `engine-development-workflow.md §4.3`) +- **Merge strategy:** squash-and-merge (cf. `engine-development-workflow.md §4.6`) + +## Notes + +- **Discrete oracle, not the continuous parabola.** Semi-implicit Euler integrates `v` then `x` from the new `v`, giving `x_n = x₀ + n·v₀·dt + g·dt²·n(n+1)/2` — after 1 s at 60 Hz, `−4.98675 m`, not the continuous `−4.905 m`. A fixed-step explicit integrator is *not supposed* to match the analytic parabola; "trajectoire correcte" means it matches its own discrete scheme. `v_n = n·g·dt` DOES match the continuous velocity — the test asserts both. The default `linear_damping` is 0.05, so the oracle test forces `linear_damping = 0`. +- **Gravity as acceleration.** `a += gravity·gravity_factor`, mass-independent (Galileo) — the demonstrator checks exactly this. `F = m·g` then `a = F·inv_mass = g` wastes a multiply and loses fp precision at large mass. +- **Damping form.** Jolt clamped-linear `v *= max(0, 1 − d·dt)`. At d=0.05, dt=1/60 the deviation from `exp(−d·dt)` is ~3.5e-7 — negligible; the clamp is physical (prevents sign flip at `d·dt > 1`), not a geometric epsilon. Under sub-stepping (M1.1.15) this form is substep-count-sensitive (`(1 − d·dt/N)^N → e^{−d·dt}`); that is the orchestrator's call — `integrate(dt)` stays a pure per-call function. +- **First-order orientation.** `q ← normalize(q + ½·dt·(ω_quat ⊗ q))`, `ω` world-space on the LEFT. The math layer was provisioned for this (`Quat.add` doc: "for angular-velocity integration; renormalize"). No `|ω|` division → no singularity at `ω = 0`, so the M1.1.4-frozen threshold discipline (guards fire only at true zero) is honored for free. Error is imperceptible at game `ω·dt`. +- **Gyroscopic term dropped** (Jolt/PhysX/Bevy default); explicit integration of `ω × (I·ω)` injects energy. Additive implicit option later. +- **Determinism.** Ascending slot-index order, fixed `dt`, no hash container on the path, `Real`-generic (f32/f64 comptime flip). The only guard is the damping clamp (physical). No absolute geometric threshold introduced — consistent with the M1.1.4-frozen discipline. +- **Codex second-opinion (verified against source) folded in:** the discrete oracle (#1), `isAliveIndex` for the non-compacting store (#2, verified in `slot_alloc.zig`), clear accumulators for all live bodies (#3), `addImpulse` in the mutator set (#4, verified in `engine-tier-interfaces.md §1`), and the removal of the contradictory substep phrasing (#5). + +--- + +# LIVING SECTION + +*Maintained by Claude Code during the milestone. The log is not a marketing report: it serves review and post-mortem debugging.* + +## Specs read + +- [x] `engine-physics-forge.md` (§1.2, §2, §1.5–1.6) — read 2026-07-23 11:21 +- [x] `engine-project-settings.md` (§3.2) — read 2026-07-23 11:21 +- [x] `engine-coordinate-system.md` (§1, §3) — read 2026-07-23 11:21 +- [x] `engine-tier-interfaces.md` (§1) — read 2026-07-23 11:21 +- [x] `engine-phase-1-plan.md` (M1.1.5 row, M1.1 arc) — read 2026-07-23 11:21 +- [x] `engine-zig-conventions.md` — read 2026-07-23 11:21 + +## Execution log + +### E1 — accumulators + mutators/getters + `isAliveIndex` + +- 2026-07-23 — (a) Added `force`/`torque` `Vec3r` columns to the SoA `Body` + row (`body.zig`), initialised to `Vec3r.zero` in `addBody` (`body_manager.zig`); + updated both file-header doc comments. Independent `MultiArrayList` columns, so + the position/rotation bulk-sync layout invariant (M1.1.15) is untouched. +- 2026-07-23 — (b) Added `IdAllocator.isAliveIndex(index) bool` (`slot_alloc.zig`) + — bare-index liveness for the non-compacting store's index-ascending pass + (`validate` takes a packed id, so it cannot serve). +- 2026-07-23 — (c) Added `BodyManager` `linearVelocity`/`angularVelocity` getters, + `setLinearVelocity`/`setAngularVelocity`, `addForce`/`addTorque` (accumulate), + `addImpulse` (immediate `Δv = impulse·inv_mass`). All validate the `BodyId` and + no-op / return null on a stale handle. `addTorque`/`setAngularVelocity` kept + INTERNAL (not on the frozen `PhysicsModule` 3D interface). +- 2026-07-23 — Tests placed in the existing `tests/body_manager_test.zig` (the + established `BodyManager` acceptance-test home) rather than a new file — this + file is outside the brief's "Files to create or modify" list; justified as the + consistent home and it avoids touching `root.zig` (an E2 file, per the brief) in + E1. `isAliveIndex` is exercised via the public `bm.alloc` field across + create/free/reuse, so `slot_alloc.zig` needs no inline test (and no `root.zig` + pin) at E1; its body is covered by that empirical consumer (§13 surface guard). +- 2026-07-23 — Green: `zig build test-forge-3d` 98/98 at f32 AND `-Dphysics_f64=true`, + Debug AND ReleaseSafe; `zig fmt --check`, `zig build lint`, `zig build` all clean. + +### E2 — `integrate` linear + gravity + damping + free-fall oracle + +- 2026-07-23 — Created `pipeline/integration.zig` with + `pub fn integrate(bm, dt, gravity)`: one ascending `0..bodies.len` pass, + `isAliveIndex`-filtered, reading the SoA columns directly. Dynamic branch: + `a = gravity·gravity_factor + force·inv_mass`; `v += a·dt`; + `v *= max(0, 1 − linear_damping·dt)`; `x += v·dt` (position from the new + velocity). Every live body then has `force`/`torque` cleared (§2, uniform). + Static/kinematic are not moved. The only guard is the damping clamp (physical, + not a geometric epsilon) — M1.1.4 threshold discipline honoured. +- 2026-07-23 — `root.zig`: imported `integration`, re-exported `integrate` at + `Real`, pinned `integration` + `tests/integration_test.zig` in the comptime + block (§13 lazy-analysis guard). +- 2026-07-23 — Created `tests/integration_test.zig` (E2 subset of the acceptance + matrix): discrete free-fall oracle (N=60 @ dt=1/60, v₀=0 → x ≈ −4.98675 m, + v ≈ −9.81 m/s; named tolerance 1e-3, documented at the site — excludes the + continuous −4.905), mass-independence, `gravity_factor` scaling (0.5 / 0), + damping bends + clamps (`d·dt>1` → v=0, no sign flip), forces consumed once, + impulse immediate Δv, static/kinematic unmoved with accumulators cleared, + freed-slot skipping + reused-slot correctness, determinism (bit-identical). +- 2026-07-23 — f64 compile fix (caught by the `-Dphysics_f64=true` sweep): the + determinism test assigned `Real`-typed exprs to the descriptor's `mass`/ + `gravity_factor` fields, which are always `f32` (matched only at `Real = f32`); + computed them as `f32`. +- 2026-07-23 — Green: `zig build test-forge-3d` 107/107 at f32 AND + `-Dphysics_f64=true`, Debug AND ReleaseSafe; `zig fmt --check`, + `zig build lint`, `zig build` all clean. + +### E3 — angular integration + close + +- 2026-07-23 — Extended `integrate`'s dynamic branch with the angular half: + `I_world_inv = R · I_local_inv · Rᵀ` (`R = Mat3r.fromQuat(rotation)`, + `I_local_inv` the diagonal `MotionProperties.local_inv_inertia`); + `α = I_world_inv · torque`; `ω += α·dt`; `ω *= max(0, 1 − angular_damping·dt)`; + `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with `ω_quat = {ω.x, ω.y, ω.z, 0}` and + `ω_quat.mul(q)` (world-space ω on the LEFT). Gyroscopic term dropped. The + first-order path divides by no `|ω|` ⇒ no zero guard (M1.1.4 discipline). + `integration.zig` gains the `Mat3r`/`Quatr` aliases; header + `integrate` doc + updated. +- 2026-07-23 — `tests/integration_test.zig` angular tests: torque on a rotated + anisotropic box (`ω = (R·I_local_inv·Rᵀ)·τ·dt`, recomputed independently from + the pre-step pose + inertia; tol 1e-5), angular damping + clamp (`d·dt>1` → 0), + orientation stays unit (|q| within 1e-5 of 1 after 120 spinning steps), + `ω = 0` leaves `q` unchanged with no NaN. Extended the E2 static/kinematic test + with orientation-unchanged + torque-cleared checks, and the determinism test to + cover rotation + angular velocity + torque (bit-identical). +- 2026-07-23 — Green: `zig build test-forge-3d` 111/111 at f32 AND + `-Dphysics_f64=true`, Debug AND ReleaseSafe; `zig fmt --check`, + `zig build lint`, `zig build` all clean. +- 2026-07-23 — Codex review findings P2×2 (test-only; production code untouched), + both in `tests/integration_test.zig`: (1) `test "freed slots are skipped"` now + seeds the freed slot with a distinctive velocity + force and snapshots its raw + columns after the free, asserting bit-exact equality after `integrate` — the + dead slot's position is not integrated, its velocity is untouched by gravity, + and its stale force is NOT cleared (locking that the §2 reset is live-only). + (2) New `test "orientation update uses the left (world-space) quaternion + product"` — a non-commutative (q0, ω) config with an in-site left-vs-right + oracle, a discrimination guard (`!left.approxEql(right)`), and matches the + engine to the LEFT product only. Re-validated 112/112 across the full matrix + (f32 + `-Dphysics_f64=true`, Debug + ReleaseSafe; fmt/lint/build clean). +- 2026-07-23 — Closure steps (CLAUDE.md §3.4 patch, §3.6.1 language + drift + audits, Closing notes, `Status: CLOSED`) pending the E3 GO and the Claude.ai + `CLAUDE.md` patch. + +## Recorded deviations + +## Blockers encountered + +## Closing notes + +**Delivered.** The first sub-milestone where bodies move. `pipeline/integration.zig` +`integrate(bm, dt, gravity)` — one pure per-tick pass, ascending slot-index order, +`isAliveIndex`-filtered over the non-compacting SoA store, free-flight only +(broadphase/narrowphase exist but are never invoked; contact response is M1.1.6). +Linear: gravity as an acceleration × `gravity_factor` (mass-independent), Jolt +clamped-linear damping, position from the new velocity. Angular: +`α = (R·I_local_inv·Rᵀ)·τ`, clamp-damped `ω`, first-order orientation +`q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with world-space ω on the LEFT; gyroscopic +term dropped. `Body` gained `force`/`torque` accumulator columns; `BodyManager` +gained the stale-safe velocity getters/setters, `addForce`/`addTorque`, +`addImpulse` (`addTorque`/`setAngularVelocity` internal — no public interface +change). No EPA/GJK/manifold file touched; `PhysicsModule` not instantiated. + +**Threshold discipline (M1.1.4 class).** The only guard is the damping clamp +`max(0, 1 − d·dt)` (physical, not a geometric epsilon). The first-order quaternion +path divides by no `|ω|` and needs no zero guard — `ω_quat ⊗ q ⊥ q` in 4D, so +`‖q + dq‖ ≥ 1` and the renormalisation is never near-zero. + +**Tests.** 112/112 green at f32 AND `-Dphysics_f64=true`, Debug AND ReleaseSafe; +`zig fmt --check`, `zig build lint`, `zig build` clean. Discrete semi-implicit +free-fall oracle (−4.98675 m / −9.81 m/s, tolerance justified at the site), +mass-independence, `gravity_factor` scaling, damping bend + clamp, force +consume-once, impulse, static/kinematic invariance (position + orientation + +accumulators-cleared), freed-slot skip observed bit-exact, angular closed form on +a rotated anisotropic box, angular damping + clamp, unit-quaternion preservation, +ω=0 no-op/no-NaN, left-product lock, and full-integrator determinism (linear + +angular, bit-identical). + +**Recorded deviations / decisions** (detail in Execution log): (1) E1 mutator/ +getter/`isAliveIndex` tests placed in the existing `tests/body_manager_test.zig` +(BodyManager's established acceptance-test home) — outside the brief's file list, +justified. (2) f32 fix: descriptor `mass`/`gravity_factor` are always `f32` +(caught by the f64 sweep). (3) Codex review P2×2, TEST-ONLY (production untouched): +dead-slot skip now observed bit-exact incl. stale-force-not-cleared, and the LEFT +quaternion product locked by a non-commutative one-step oracle with a +discrimination guard. + +**§3.6.1 audits.** Language: production code, `CLAUDE.md`, and the brief's LIVING +section are French-free (one leaked French spec quote in `body.zig`'s doc comment +was removed — `docs(forge): drop french spec quote…`). The brief's FROZEN SECTION +retains verbatim French spec citations authored by Claude.ai — the Context +demonstrator label ("Box en chute libre, trajectoire correcte"), the plan row +("Integration : … + gravité + damping"), and the Specs `§2 "reset chaque fixed +tick"` contract quote — left untouched per the FROZEN-section rule (not CC prose; +they are quotations of the French KB spec). Drift: the introduced symbols +(`integrate`, `isAliveIndex`, `force`/`torque`, the velocity/force mutators) are +referenced consistently across `integration.zig`/`root.zig`/`body_manager.zig`/ +tests — no orphaned or half-renamed references. + +**Out of scope — confirmed untouched.** No solver / Sequential Impulses / friction +/ restitution / warm-start, no sleeping, no CCD, no substep handling, no kinematic +position-from-velocity, no gyroscopic term, no ECS↔BodyManager sync, no +`PhysicsModule` wiring, no bench, no spec-file edit, no EPA hotfix. diff --git a/src/modules/forge/forge_3d/body.zig b/src/modules/forge/forge_3d/body.zig index 518a05a..bf5c68e 100644 --- a/src/modules/forge/forge_3d/body.zig +++ b/src/modules/forge/forge_3d/body.zig @@ -4,6 +4,12 @@ //! integrates against; `computeMotion` derives them from a `BodyDescriptor` and //! its `Shape`. Static and kinematic bodies have zero inverse mass and inertia //! (infinite effective mass). `Body` is the SoA row `BodyManager` stores. +//! +//! The `force`/`torque` columns are world-space per-tick accumulators (N, N·m): +//! `BodyManager.addForce`/`addTorque` add into them and `integrate` (E2) reads +//! then clears them each fixed tick (the `engine-physics-forge.md` §2 per-tick +//! reset contract). Being independent SoA columns, they leave the +//! position/rotation bulk-sync layout invariant (M1.1.15) untouched. const std = @import("std"); const api = @import("weld_forge"); @@ -52,6 +58,12 @@ pub const Body = struct { linear_velocity: Vec3r, /// World-space angular velocity. angular_velocity: Vec3r, + /// Accumulated world-space force (N) for the current tick; `addForce` + /// accumulates, `integrate` reads then clears it each fixed tick. + force: Vec3r, + /// Accumulated world-space torque (N·m) for the current tick; `addTorque` + /// accumulates, `integrate` reads then clears it each fixed tick. + torque: Vec3r, /// Derived inverse mass/inertia + damping/gravity. motion: MotionProperties, /// Collision-shape handle (into a `ShapeStore`). diff --git a/src/modules/forge/forge_3d/body_manager.zig b/src/modules/forge/forge_3d/body_manager.zig index cb4a50d..31f5cf2 100644 --- a/src/modules/forge/forge_3d/body_manager.zig +++ b/src/modules/forge/forge_3d/body_manager.zig @@ -12,6 +12,15 @@ //! Execution-log note flagging this against decision 7. Id allocation is //! deterministic (no hash-map on the path — M1.1.14). World AABBs are computed //! exactly per primitive on demand. +//! +//! Velocity/force/torque mutators (`setLinearVelocity`/`setAngularVelocity`, +//! `addForce`/`addTorque`/`addImpulse`) and their getters validate the handle +//! and no-op / return null on a stale one (parity with `position`/`rotation`). +//! `addForce`/`addTorque` accumulate into the per-tick `force`/`torque` columns +//! (cleared by `integrate`); `addImpulse` is an immediate `Δv = impulse·inv_mass` +//! (a natural no-op on static/kinematic bodies, `inv_mass == 0`). `addTorque` +//! and `setAngularVelocity` are INTERNAL — the public `PhysicsModule` 3D +//! interface (frozen M1.1.15) carries no angular mutators. const std = @import("std"); const api = @import("weld_forge"); @@ -65,6 +74,8 @@ pub const BodyManager = struct { .rotation = convQuat(desc.rotation), .linear_velocity = Vec3r.zero, .angular_velocity = Vec3r.zero, + .force = Vec3r.zero, + .torque = Vec3r.zero, .motion = body_mod.computeMotion(desc, shape), .shape = desc.shape, .body_type = desc.body_type, @@ -120,6 +131,62 @@ pub const BodyManager = struct { return self.bodies.items(.collision_layer)[idx]; } + /// Safe getter: world-space linear velocity, or null if `id` is stale/invalid. + pub fn linearVelocity(self: *const BodyManager, id: BodyId) ?Vec3r { + const idx = self.alloc.validate(id) orelse return null; + return self.bodies.items(.linear_velocity)[idx]; + } + + /// Safe getter: world-space angular velocity, or null if `id` is stale/invalid. + pub fn angularVelocity(self: *const BodyManager, id: BodyId) ?Vec3r { + const idx = self.alloc.validate(id) orelse return null; + return self.bodies.items(.angular_velocity)[idx]; + } + + /// Set the world-space linear velocity. No-op on a stale/invalid handle. + pub fn setLinearVelocity(self: *BodyManager, id: BodyId, velocity: Vec3r) void { + const idx = self.alloc.validate(id) orelse return; + self.bodies.items(.linear_velocity)[idx] = velocity; + } + + /// Set the world-space angular velocity. No-op on a stale/invalid handle. + /// Internal to `BodyManager`: the public `PhysicsModule` 3D interface has no + /// angular-velocity mutator (frozen decision at M1.1.15). + pub fn setAngularVelocity(self: *BodyManager, id: BodyId, velocity: Vec3r) void { + const idx = self.alloc.validate(id) orelse return; + self.bodies.items(.angular_velocity)[idx] = velocity; + } + + /// Accumulate a world-space force (N) into the body's per-tick force + /// accumulator. No-op on a stale/invalid handle. Any live body accumulates + /// (the per-tick clear in `integrate` is uniform); a force must be + /// re-applied every tick it should act. + pub fn addForce(self: *BodyManager, id: BodyId, force: Vec3r) void { + const idx = self.alloc.validate(id) orelse return; + const forces = self.bodies.items(.force); + forces[idx] = forces[idx].add(force); + } + + /// Accumulate a world-space torque (N·m) into the body's per-tick torque + /// accumulator. No-op on a stale/invalid handle. Internal to `BodyManager` + /// (the interface has no torque mutator — see `setAngularVelocity`). + pub fn addTorque(self: *BodyManager, id: BodyId, torque: Vec3r) void { + const idx = self.alloc.validate(id) orelse return; + const torques = self.bodies.items(.torque); + torques[idx] = torques[idx].add(torque); + } + + /// Apply an instantaneous linear impulse (N·s) as an IMMEDIATE velocity + /// change `Δv = impulse · inv_mass`. No-op on a stale/invalid handle, and + /// naturally a no-op on a static/kinematic body (`inv_mass == 0`) — no + /// `body_type` branch needed. + pub fn addImpulse(self: *BodyManager, id: BodyId, impulse: Vec3r) void { + const idx = self.alloc.validate(id) orelse return; + const inv_mass = self.bodies.items(.motion)[idx].inv_mass; + const vel = self.bodies.items(.linear_velocity); + vel[idx] = vel[idx].add(impulse.scale(inv_mass)); + } + /// Safe getter: the exact world-space AABB of the body's shape, or null if /// `id` (or its shape) is stale/invalid. pub fn bodyAabb(self: *const BodyManager, store: *const ShapeStore, id: BodyId) ?Aabbr { diff --git a/src/modules/forge/forge_3d/pipeline/integration.zig b/src/modules/forge/forge_3d/pipeline/integration.zig new file mode 100644 index 0000000..a81daac --- /dev/null +++ b/src/modules/forge/forge_3d/pipeline/integration.zig @@ -0,0 +1,107 @@ +//! `forge_3d/pipeline/integration.zig` — semi-implicit (symplectic) Euler +//! integration over the live `BodyManager` SoA store. +//! +//! One pure per-tick pass in ascending slot-index order (deterministic — no hash +//! container on the path, the `BodyManager` discipline; M1.1.14). The store does +//! not compact, so the pass walks `0..bodies.len` and filters liveness with +//! `IdAllocator.isAliveIndex`. This is FREE-FLIGHT integration: broadphase and +//! narrowphase exist but are NOT invoked here, and there is no contact response +//! (Sequential Impulses is M1.1.6) — a dynamic body under gravity follows an +//! unobstructed trajectory. +//! +//! Semi-implicit Euler: velocity is integrated first, position from the *new* +//! velocity. Gravity is applied as an ACCELERATION scaled by `gravity_factor` +//! (mass-independent — Galileo), never as a force. Damping is the Jolt +//! clamped-linear form `v *= max(0, 1 − d·dt)`; the clamp is physical (it +//! prevents a sign flip when `d·dt > 1`), not a geometric epsilon. The +//! `force`/`torque` accumulators are reset every fixed tick for ALL live bodies +//! (`engine-physics-forge.md` §2). `integrate` applies damping exactly once with +//! the `dt` it is given and has no opinion on substeps — call cadence is the +//! orchestrator's concern (M1.1.15). +//! +//! Angular integration mirrors the linear half: the torque is mapped through the +//! world-space inverse inertia `I_world_inv = R · I_local_inv · Rᵀ`, `ω` is +//! integrated then clamp-damped, and the orientation advances by the first-order +//! (linearized) rule `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with the world-space +//! `ω` on the LEFT. That form divides by no `|ω|`, so it is singularity-free at +//! `ω = 0` (no zero guard needed — the M1.1.4 threshold discipline is honoured). +//! The gyroscopic term `ω × (I·ω)` is dropped (Jolt/PhysX/Bevy default; its +//! explicit integration injects energy). + +const config = @import("../config.zig"); +const body_manager = @import("../body_manager.zig"); + +const Real = config.Real; +const Vec3r = config.Vec3r; +const Quatr = config.Quatr; +const Mat3r = config.Mat3r; +const BodyManager = body_manager.BodyManager; + +/// Advance every live body one fixed tick of `dt` seconds under world-space +/// `gravity` (m/s²), in ascending slot-index order. +/// +/// For each DYNAMIC body — linear: `a = gravity·gravity_factor + force·inv_mass`; +/// `v += a·dt`; `v *= max(0, 1 − linear_damping·dt)`; `x += v·dt` (position from +/// the new velocity). Angular: `α = (R·I_local_inv·Rᵀ)·torque`; `ω += α·dt`; +/// `ω *= max(0, 1 − angular_damping·dt)`; `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` +/// with world-space `ω` on the left. Static and kinematic bodies are not moved +/// (kinematic position-from-velocity is a later additive concern). Every live +/// body then has its `force`/`torque` accumulators cleared for the next tick. +/// +/// This is a pure function of the store and `(dt, gravity)`: no broadphase, +/// narrowphase, contacts, or sleeping. The gyroscopic term is dropped. +pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { + const positions = bm.bodies.items(.position); + const rotations = bm.bodies.items(.rotation); + const linear_velocities = bm.bodies.items(.linear_velocity); + const angular_velocities = bm.bodies.items(.angular_velocity); + const forces = bm.bodies.items(.force); + const torques = bm.bodies.items(.torque); + const motions = bm.bodies.items(.motion); + const body_types = bm.bodies.items(.body_type); + + const n: u32 = @intCast(bm.bodies.len); + var i: u32 = 0; + while (i < n) : (i += 1) { + // The store never compacts, so dead slots sit between live ones — skip + // them (their stale data must not be integrated). + if (!bm.alloc.isAliveIndex(i)) continue; + + if (body_types[i] == .dynamic) { + const mp = motions[i]; + + // --- Linear --- + // Gravity as an acceleration (mass-independent), plus the accumulated + // force divided by mass. Read `force` before the clear below. + const accel = gravity.scale(mp.gravity_factor).add(forces[i].scale(mp.inv_mass)); + // Integrate velocity first (symplectic), then damp, then position. + var v = linear_velocities[i].add(accel.scale(dt)); + const damp = @max(@as(Real, 0), 1 - mp.linear_damping * dt); + v = v.scale(damp); + linear_velocities[i] = v; + positions[i] = positions[i].add(v.scale(dt)); + + // --- Angular --- + // World-space inverse inertia I_world_inv = R · I_local_inv · Rᵀ + // (Rᵀ == R⁻¹ for an orthonormal rotation), then α = I_world_inv·τ. + const r = Mat3r.fromQuat(rotations[i]); + const i_world_inv = r.mul(mp.local_inv_inertia).mul(r.transpose()); + const alpha = i_world_inv.mulVec(torques[i]); + var w = angular_velocities[i].add(alpha.scale(dt)); + const ang_damp = @max(@as(Real, 0), 1 - mp.angular_damping * dt); + w = w.scale(ang_damp); + angular_velocities[i] = w; + // First-order orientation update: q ← normalize(q + ½·dt·(ω_quat ⊗ q)), + // world-space ω on the LEFT. No |ω| division ⇒ singularity-free at ω = 0. + const wa = w.toArray(); + const w_quat = Quatr{ .x = wa[0], .y = wa[1], .z = wa[2], .w = 0 }; + const dq = w_quat.mul(rotations[i]).scale(0.5 * dt); + rotations[i] = rotations[i].add(dq).normalize(); + } + + // Reset the per-tick accumulators for every live body (§2), including + // static/kinematic — the clear is uniform. + forces[i] = Vec3r.zero; + torques[i] = Vec3r.zero; + } +} diff --git a/src/modules/forge/forge_3d/root.zig b/src/modules/forge/forge_3d/root.zig index f488b8b..c8dbb2e 100644 --- a/src/modules/forge/forge_3d/root.zig +++ b/src/modules/forge/forge_3d/root.zig @@ -20,6 +20,9 @@ const broadphase = @import("pipeline/broadphase.zig"); // (engine-zig-conventions.md §13 lazy-analysis guard — an unreferenced module's // tests are silently skipped). const narrowphase = @import("pipeline/narrowphase/root.zig"); +// M1.1.5 — semi-implicit Euler integration over the `BodyManager` SoA store. +// Re-exported at `Real` below; the comptime pin analyses its acceptance tests. +const integration = @import("pipeline/integration.zig"); // --- Solver scalar + math aliases --- @@ -129,6 +132,16 @@ pub fn collideOrderedGeneric(shape_a: SupportShape, pos_a: Vec3r, rot_a: Quatr, /// The `(normal, closest points, base penetration)` fast-path seed at solver precision. pub const ContactSeed = narrowphase.ContactSeed(Real); +// --- Integration (semi-implicit Euler) --- + +/// Advance every live body one fixed tick of `dt` under world-space `gravity` +/// (m/s²) with semi-implicit Euler + gravity·gravity_factor + clamped-linear +/// damping — the `Real`-bound entry. Free-flight only (no contacts); called by +/// the `step` orchestrator at M1.1.15. +pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { + return integration.integrate(bm, dt, gravity); +} + // Pins so the inline tests + the acceptance suite are analysed when this module // is built as a test target (engine-zig-conventions.md §13). comptime { @@ -138,7 +151,9 @@ comptime { _ = body_manager; _ = broadphase; _ = narrowphase; + _ = integration; _ = @import("tests/body_manager_test.zig"); + _ = @import("tests/integration_test.zig"); _ = @import("tests/broadphase_test.zig"); _ = @import("tests/gjk_test.zig"); _ = @import("tests/epa_test.zig"); diff --git a/src/modules/forge/forge_3d/slot_alloc.zig b/src/modules/forge/forge_3d/slot_alloc.zig index 22eedf2..5017e22 100644 --- a/src/modules/forge/forge_3d/slot_alloc.zig +++ b/src/modules/forge/forge_3d/slot_alloc.zig @@ -94,4 +94,14 @@ pub const IdAllocator = struct { if (!meta.alive or meta.generation != p.generation) return null; return p.index; } + + /// Whether the slot at bare column `index` is currently live. Distinct from + /// `validate` (which takes a packed id and checks the generation): `free` + /// does NOT compact, so `slots.items.len` is the high-water mark including + /// dead slots. An index-ascending pass (the integrator) walks that range and + /// must filter liveness per slot — a bare index carries no generation, so + /// `validate` cannot serve here. + pub fn isAliveIndex(self: *const IdAllocator, index: u32) bool { + return index < self.slots.items.len and self.slots.items[index].alive; + } }; diff --git a/src/modules/forge/forge_3d/tests/body_manager_test.zig b/src/modules/forge/forge_3d/tests/body_manager_test.zig index 842693b..125d85d 100644 --- a/src/modules/forge/forge_3d/tests/body_manager_test.zig +++ b/src/modules/forge/forge_3d/tests/body_manager_test.zig @@ -302,3 +302,138 @@ test "rotation getter round-trips and rejects stale handles" { bm.removeBody(id); try testing.expect(bm.rotation(id) == null); } + +// --- E1 velocity / force / torque / impulse mutators & getters --------------- + +test "linear and angular velocity set/get round-trip and reject stale handles" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const id = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + + // Velocities start at zero. + try testing.expect(bm.linearVelocity(id).?.approxEql(Vec3r.zero, 0)); + try testing.expect(bm.angularVelocity(id).?.approxEql(Vec3r.zero, 0)); + + const lin = vr(1, -2, 3); + const ang = vr(-0.5, 4, 0.25); + bm.setLinearVelocity(id, lin); + bm.setAngularVelocity(id, ang); + try testing.expect(bm.linearVelocity(id).?.approxEql(lin, 0)); + try testing.expect(bm.angularVelocity(id).?.approxEql(ang, 0)); + + // Stale handle ⇒ getters null, setters no-op (no crash). + bm.removeBody(id); + try testing.expect(bm.linearVelocity(id) == null); + try testing.expect(bm.angularVelocity(id) == null); + bm.setLinearVelocity(id, vr(9, 9, 9)); // no-op + bm.setAngularVelocity(id, vr(9, 9, 9)); // no-op + try testing.expect(!bm.isValid(id)); +} + +test "addForce and addTorque accumulate into the per-tick columns" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const id = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + const idx = bm.alloc.validate(id).?; + + // Accumulators start at zero. + try testing.expect(bm.bodies.items(.force)[idx].approxEql(Vec3r.zero, 0)); + try testing.expect(bm.bodies.items(.torque)[idx].approxEql(Vec3r.zero, 0)); + + bm.addForce(id, vr(1, 0, 0)); + bm.addForce(id, vr(0, 2, -1)); + bm.addTorque(id, vr(0, 0, 3)); + bm.addTorque(id, vr(1, 0, 0)); + try testing.expect(bm.bodies.items(.force)[idx].approxEql(vr(1, 2, -1), 1e-6)); + try testing.expect(bm.bodies.items(.torque)[idx].approxEql(vr(1, 0, 3), 1e-6)); + + // Stale handle ⇒ no-op. + bm.removeBody(id); + bm.addForce(id, vr(5, 5, 5)); // no-op, no crash + bm.addTorque(id, vr(5, 5, 5)); // no-op, no crash + try testing.expect(!bm.isValid(id)); +} + +test "addForce accumulates on a static body too" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const id = try bm.addBody(gpa, &store, descOf(0, .static, s)); + const idx = bm.alloc.validate(id).?; + + // Any live body accumulates; the uniform per-tick clear is `integrate`'s job. + bm.addForce(id, vr(7, 0, 0)); + try testing.expect(bm.bodies.items(.force)[idx].approxEql(vr(7, 0, 0), 1e-6)); +} + +test "impulse is an immediate velocity change and is a no-op on a static body" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + + // Dynamic: Δv = impulse · inv_mass, applied immediately (no integrate). + var dyn = descOf(0, .dynamic, s); + dyn.mass = 4; + const dyn_id = try bm.addBody(gpa, &store, dyn); + bm.addImpulse(dyn_id, vr(8, 0, -4)); + try testing.expect(bm.linearVelocity(dyn_id).?.approxEql(vr(2, 0, -1), 1e-6)); // /4 + + // A second impulse accumulates onto the current velocity. + bm.addImpulse(dyn_id, vr(0, 4, 0)); + try testing.expect(bm.linearVelocity(dyn_id).?.approxEql(vr(2, 1, -1), 1e-6)); + + // Static: inv_mass == 0 ⇒ no velocity change. + const stat_id = try bm.addBody(gpa, &store, descOf(1, .static, s)); + bm.addImpulse(stat_id, vr(100, 100, 100)); + try testing.expect(bm.linearVelocity(stat_id).?.approxEql(Vec3r.zero, 0)); + + // Stale handle ⇒ no-op. + bm.removeBody(dyn_id); + bm.addImpulse(dyn_id, vr(1, 1, 1)); // no-op, no crash + try testing.expect(!bm.isValid(dyn_id)); +} + +test "isAliveIndex tracks liveness across create, free, and reuse" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + // Bare index 0 is not alive before any slot exists. + try testing.expect(!bm.alloc.isAliveIndex(0)); + + const a = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + const b = try bm.addBody(gpa, &store, descOf(1, .dynamic, s)); + const ia = api.PackedId.unpack(a).index; + const ib = api.PackedId.unpack(b).index; + try testing.expect(bm.alloc.isAliveIndex(ia)); + try testing.expect(bm.alloc.isAliveIndex(ib)); + try testing.expect(!bm.alloc.isAliveIndex(ib + 1)); // out of range + + // Free b: `removeBody` does NOT compact, so the slot count is unchanged but + // that index is now dead — exactly what an index-ascending pass must skip. + bm.removeBody(b); + try testing.expect(bm.alloc.isAliveIndex(ia)); + try testing.expect(!bm.alloc.isAliveIndex(ib)); + + // Reuse b's slot (LIFO): the same bare index goes live again. + const c = try bm.addBody(gpa, &store, descOf(2, .dynamic, s)); + try testing.expectEqual(ib, api.PackedId.unpack(c).index); + try testing.expect(bm.alloc.isAliveIndex(ib)); +} diff --git a/src/modules/forge/forge_3d/tests/integration_test.zig b/src/modules/forge/forge_3d/tests/integration_test.zig new file mode 100644 index 0000000..fba8679 --- /dev/null +++ b/src/modules/forge/forge_3d/tests/integration_test.zig @@ -0,0 +1,543 @@ +//! M1.1.5 acceptance suite for semi-implicit Euler integration. E2 covers the +//! linear half (gravity as acceleration, clamped-linear damping, position from +//! the new velocity) plus the discrete free-fall oracle, force consumption, +//! impulse, static/kinematic invariance, freed-slot skipping, and determinism. +//! E3 adds the angular tests to this same file. + +const std = @import("std"); +const config = @import("../config.zig"); +const shape_mod = @import("../shape.zig"); +const bm_mod = @import("../body_manager.zig"); +const integration = @import("../pipeline/integration.zig"); +const api = @import("weld_forge"); +const math = @import("foundation").math; + +const Real = config.Real; +const Vec3r = config.Vec3r; +const Quatr = config.Quatr; +const Mat3r = config.Mat3r; +const ShapeStore = shape_mod.ShapeStore; +const BodyManager = bm_mod.BodyManager; +const Vec3 = math.Vec3; // f32 descriptor vector +const Quatf = math.Quatf; // f32 descriptor quaternion +const testing = std.testing; + +// --- helpers ----------------------------------------------------------------- + +fn vr(x: Real, y: Real, z: Real) Vec3r { + return Vec3r.fromArray(.{ x, y, z }); +} + +/// Euclidean norm of a quaternion (Quat.length is not public). +fn quatNorm(q: Quatr) Real { + const a = q.toArray(); + return @sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3]); +} + +/// A dynamic descriptor with damping forced to 0 (the descriptor default is +/// 0.05; the integration tests set damping explicitly where they exercise it). +fn dynDesc(entity_index: u32, shape: api.ShapeId) api.BodyDescriptor { + return .{ + .entity = .{ .index = entity_index, .generation = 0 }, + .body_type = .dynamic, + .shape = shape, + .linear_damping = 0, + }; +} + +// --- E2 tests ---------------------------------------------------------------- + +test "free fall matches the discrete semi-implicit oracle" { + // Semi-implicit Euler integrates v then x from the *new* v, so after N steps + // x = x0 + N·v0·dt + g·dt²·N(N+1)/2 — NOT the continuous parabola. At 60 Hz + // over 1 s with v0 = 0 that is −4.98675 m (the continuous value is −4.905, + // 0.08 m away), while v = N·g·dt = −9.81 m/s matches the continuous velocity. + // + // Named tolerance: 60 f32 accumulations at |x| ≈ 5 carry ~60 ULP ≈ 1e-4 of + // rounding; 1e-3 bounds that and still excludes the continuous −4.905. f64 is + // ~9 orders tighter, so the same bound holds there too. + const oracle_tol: Real = 1e-3; + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const id = try bm.addBody(gpa, &store, dynDesc(0, s)); + + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + var step: u32 = 0; + while (step < 60) : (step += 1) integration.integrate(&bm, dt, g); + + const pos = bm.position(id).?.toArray(); + const vel = bm.linearVelocity(id).?.toArray(); + try testing.expectApproxEqAbs(@as(Real, -4.98675), pos[1], oracle_tol); + try testing.expectApproxEqAbs(@as(Real, -9.81), vel[1], oracle_tol); + // X and Z stay put — gravity is Y-only. + try testing.expectApproxEqAbs(@as(Real, 0), pos[0], oracle_tol); + try testing.expectApproxEqAbs(@as(Real, 0), pos[2], oracle_tol); +} + +test "free fall is mass-independent" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + var light = dynDesc(0, s); + light.mass = 1; + var heavy = dynDesc(1, s); + heavy.mass = 1000; + const id_light = try bm.addBody(gpa, &store, light); + const id_heavy = try bm.addBody(gpa, &store, heavy); + + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + var step: u32 = 0; + while (step < 30) : (step += 1) integration.integrate(&bm, dt, g); + + // Gravity is an acceleration, so the two trajectories are bit-identical. + const pl = bm.position(id_light).?.toArray(); + const ph = bm.position(id_heavy).?.toArray(); + const vl = bm.linearVelocity(id_light).?.toArray(); + const vh = bm.linearVelocity(id_heavy).?.toArray(); + inline for (0..3) |k| { + try testing.expectEqual(pl[k], ph[k]); + try testing.expectEqual(vl[k], vh[k]); + } +} + +test "gravity_factor scales the fall" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + var full = dynDesc(0, s); + full.gravity_factor = 1.0; + var half = dynDesc(1, s); + half.gravity_factor = 0.5; + var none = dynDesc(2, s); + none.gravity_factor = 0.0; + const id_full = try bm.addBody(gpa, &store, full); + const id_half = try bm.addBody(gpa, &store, half); + const id_none = try bm.addBody(gpa, &store, none); + + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + integration.integrate(&bm, dt, g); + + const v_full = bm.linearVelocity(id_full).?.toArray()[1]; + const v_half = bm.linearVelocity(id_half).?.toArray()[1]; + const v_none = bm.linearVelocity(id_none).?.toArray()[1]; + try testing.expectApproxEqAbs(@as(Real, -9.81) * dt, v_full, 1e-6); + try testing.expectApproxEqAbs(v_full * 0.5, v_half, 1e-6); // half the acceleration + try testing.expectEqual(@as(Real, 0), v_none); // no gravity contribution +} + +test "linear damping bends the trajectory and clamps" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + + // Bend: after one step, |v_damped| < |v_undamped| (both fall, damped less). + var undamped = dynDesc(0, s); + undamped.linear_damping = 0; + var damped = dynDesc(1, s); + damped.linear_damping = 0.5; + const id_undamped = try bm.addBody(gpa, &store, undamped); + const id_damped = try bm.addBody(gpa, &store, damped); + integration.integrate(&bm, dt, g); + const v_undamped = bm.linearVelocity(id_undamped).?.toArray()[1]; + const v_damped = bm.linearVelocity(id_damped).?.toArray()[1]; + try testing.expect(v_damped < 0); // still moving down + try testing.expect(v_damped > v_undamped); // less negative ⇒ smaller magnitude + + // Clamp: d·dt > 1 ⇒ factor max(0, 1−d·dt) = 0 ⇒ velocity zeroed, no sign flip. + var clamp = dynDesc(2, s); + clamp.linear_damping = 100; // 100/60 = 1.667 > 1 + const id_clamp = try bm.addBody(gpa, &store, clamp); + bm.setLinearVelocity(id_clamp, vr(3, -5, 2)); + integration.integrate(&bm, dt, vr(0, 0, 0)); // no gravity, isolate the clamp + const v_clamp = bm.linearVelocity(id_clamp).?; + try testing.expect(v_clamp.approxEql(Vec3r.zero, 0)); // exactly zero, no flip +} + +test "forces are consumed once per tick" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + var d = dynDesc(0, s); + d.mass = 2; + const id = try bm.addBody(gpa, &store, d); + const idx = bm.alloc.validate(id).?; + + const dt: Real = 1.0 / 60.0; + const zero_gravity = vr(0, 0, 0); + bm.addForce(id, vr(10, 0, 0)); + integration.integrate(&bm, dt, zero_gravity); // a = F/m = 5 ⇒ v = 5·dt + const v_after_1 = bm.linearVelocity(id).?.toArray()[0]; + try testing.expectApproxEqAbs(@as(Real, 5) * dt, v_after_1, 1e-6); + // Accumulator cleared by the first tick. + try testing.expect(bm.bodies.items(.force)[idx].approxEql(Vec3r.zero, 0)); + + // Second tick with no re-apply: no force contribution ⇒ velocity unchanged. + integration.integrate(&bm, dt, zero_gravity); + const v_after_2 = bm.linearVelocity(id).?.toArray()[0]; + try testing.expectEqual(v_after_1, v_after_2); +} + +test "impulse is an immediate velocity change" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + + // Dynamic: Δv = impulse · inv_mass, applied immediately (before any integrate). + var d = dynDesc(0, s); + d.mass = 4; + const id = try bm.addBody(gpa, &store, d); + bm.addImpulse(id, vr(8, 0, -4)); + try testing.expect(bm.linearVelocity(id).?.approxEql(vr(2, 0, -1), 1e-6)); + + // Static: inv_mass == 0 ⇒ no change. + const stat = try bm.addBody(gpa, &store, .{ + .entity = .{ .index = 1, .generation = 0 }, + .body_type = .static, + .shape = s, + }); + bm.addImpulse(stat, vr(100, 100, 100)); + try testing.expect(bm.linearVelocity(stat).?.approxEql(Vec3r.zero, 0)); +} + +test "static and kinematic bodies are not integrated" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + + const rot = Quatf.fromAxisAngle(Vec3.unit_y, 0.7); + const rot_r = Quatr.fromArray(.{ rot.x, rot.y, rot.z, rot.w }); + + // Static at a fixed pose + orientation, with a force and torque applied. + var stat = api.BodyDescriptor{ + .entity = .{ .index = 0, .generation = 0 }, + .body_type = .static, + .shape = s, + }; + stat.position = Vec3.fromArray(.{ 1, 2, 3 }); + stat.rotation = rot; + const id_stat = try bm.addBody(gpa, &store, stat); + bm.addForce(id_stat, vr(50, 50, 50)); + bm.addTorque(id_stat, vr(9, 9, 9)); + + // Kinematic with a velocity set (no position-from-velocity in M1.1.5), + // an orientation, and a force + torque. + var kin = api.BodyDescriptor{ + .entity = .{ .index = 1, .generation = 0 }, + .body_type = .kinematic, + .shape = s, + }; + kin.position = Vec3.fromArray(.{ 4, 5, 6 }); + kin.rotation = rot; + const id_kin = try bm.addBody(gpa, &store, kin); + bm.setLinearVelocity(id_kin, vr(1, 0, 0)); + bm.setAngularVelocity(id_kin, vr(0, 2, 0)); + bm.addForce(id_kin, vr(50, 50, 50)); + bm.addTorque(id_kin, vr(9, 9, 9)); + + integration.integrate(&bm, dt, g); + + // Neither moved or rotated. + try testing.expect(bm.position(id_stat).?.approxEql(vr(1, 2, 3), 0)); + try testing.expect(bm.position(id_kin).?.approxEql(vr(4, 5, 6), 0)); + try testing.expect(bm.rotation(id_stat).?.approxEql(rot_r, 0)); + try testing.expect(bm.rotation(id_kin).?.approxEql(rot_r, 0)); + // Kinematic velocities are untouched (not consumed in M1.1.5). + try testing.expect(bm.linearVelocity(id_kin).?.approxEql(vr(1, 0, 0), 0)); + try testing.expect(bm.angularVelocity(id_kin).?.approxEql(vr(0, 2, 0), 0)); + // But their accumulators ARE cleared (§2 uniform reset). + const i_stat = bm.alloc.validate(id_stat).?; + const i_kin = bm.alloc.validate(id_kin).?; + try testing.expect(bm.bodies.items(.force)[i_stat].approxEql(Vec3r.zero, 0)); + try testing.expect(bm.bodies.items(.force)[i_kin].approxEql(Vec3r.zero, 0)); + try testing.expect(bm.bodies.items(.torque)[i_stat].approxEql(Vec3r.zero, 0)); + try testing.expect(bm.bodies.items(.torque)[i_kin].approxEql(Vec3r.zero, 0)); +} + +test "freed slots are skipped" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + const one_step: Real = -9.81 * dt * dt; // x after 1 semi-implicit step from rest + + const a = try bm.addBody(gpa, &store, dynDesc(0, s)); + const b = try bm.addBody(gpa, &store, dynDesc(1, s)); + + // Seed a's slot with a distinctive velocity + force before freeing it, so the + // snapshot below can prove the dead slot was left untouched by `integrate`. + bm.setLinearVelocity(a, vr(7, 3, -2)); + bm.addForce(a, vr(11, 0, 5)); + const ia = api.PackedId.unpack(a).index; + bm.removeBody(a); // dead slot between live ones — must be skipped, no crash + + // Snapshot the dead slot's raw columns after the free, before integrating. + const pos_dead = bm.bodies.items(.position)[ia].toArray(); + const vel_dead = bm.bodies.items(.linear_velocity)[ia].toArray(); + const force_dead = bm.bodies.items(.force)[ia].toArray(); + + integration.integrate(&bm, dt, g); // step 1: b falls one step, a's slot skipped + + // Double lock on the `isAliveIndex` filter: the dead slot is bit-unchanged — + // position not integrated, velocity not touched by gravity, AND the stale + // force NOT cleared (the §2 uniform reset runs only for LIVE slots, so a dead + // slot keeps its stale accumulator). Checked before `c` overwrites the slot. + try testing.expectEqual(pos_dead, bm.bodies.items(.position)[ia].toArray()); + try testing.expectEqual(vel_dead, bm.bodies.items(.linear_velocity)[ia].toArray()); + try testing.expectEqual(force_dead, bm.bodies.items(.force)[ia].toArray()); + + // A body created in the reused slot integrates correctly from its own start. + const c = try bm.addBody(gpa, &store, dynDesc(2, s)); + try testing.expectEqual( + api.PackedId.unpack(a).index, + api.PackedId.unpack(c).index, + ); // c reused a's freed slot + + integration.integrate(&bm, dt, g); // step 2: b falls again, c falls once + + try testing.expectApproxEqAbs(one_step * 3.0, bm.position(b).?.toArray()[1], 1e-6); // 2 steps: g·dt²·(1+2) + try testing.expectApproxEqAbs(one_step, bm.position(c).?.toArray()[1], 1e-6); // 1 step +} + +test "integration is deterministic" { + const gpa = testing.allocator; + const Out = struct { + pos: [4][3]Real = undefined, + vel: [4][3]Real = undefined, + rot: [4][4]Real = undefined, + ang: [4][3]Real = undefined, + }; + const Runner = struct { + fn run(g_alloc: std.mem.Allocator, out: *Out) !void { + var store = ShapeStore{}; + defer store.deinit(g_alloc); + var bm = BodyManager{}; + defer bm.deinit(g_alloc); + const box = try store.createShape(g_alloc, .{ .box = .{ .half_extents = Vec3.fromArray(.{ 1, 2, 3 }) } }); + + var ids: [4]api.BodyId = undefined; + inline for (0..4) |k| { + const kf = @as(f32, @floatFromInt(k)); // f32 descriptor fields + const kr = @as(Real, @floatFromInt(k)); // Real velocities/torques + var d = dynDesc(@intCast(k), box); + d.mass = kf + 1.0; + d.gravity_factor = 1.0 - kf * 0.2; + d.rotation = Quatf.fromAxisAngle(Vec3.fromArray(.{ 1, kf + 1, 0.5 }).normalize(), 0.3 + kf * 0.2); + ids[k] = try bm.addBody(g_alloc, &store, d); + bm.setLinearVelocity(ids[k], vr(kr, 0, -kr)); + bm.setAngularVelocity(ids[k], vr(0.5 * kr, 1.0, -0.3 * kr)); + bm.addTorque(ids[k], vr(kr, 0.5, -kr)); // consumed on the first tick + } + + const dt: Real = 1.0 / 60.0; + const g = vr(0.5, -9.81, 0.25); + var step: u32 = 0; + while (step < 45) : (step += 1) integration.integrate(&bm, dt, g); + + inline for (0..4) |k| { + out.pos[k] = bm.position(ids[k]).?.toArray(); + out.vel[k] = bm.linearVelocity(ids[k]).?.toArray(); + out.rot[k] = bm.rotation(ids[k]).?.toArray(); + out.ang[k] = bm.angularVelocity(ids[k]).?.toArray(); + } + } + }; + var a: Out = .{}; + var b: Out = .{}; + try Runner.run(gpa, &a); + try Runner.run(gpa, &b); + inline for (0..4) |k| { + try testing.expectEqual(a.pos[k], b.pos[k]); // bit-identical + try testing.expectEqual(a.vel[k], b.vel[k]); + try testing.expectEqual(a.rot[k], b.rot[k]); + try testing.expectEqual(a.ang[k], b.ang[k]); + } +} + +// --- E3 angular tests -------------------------------------------------------- + +test "torque on a rotated anisotropic box" { + // With angular_damping = 0 and ω₀ = 0, one step gives + // ω = (R·I_local_inv·Rᵀ)·τ·dt for the pre-step rotation R = fromQuat(q₀) and + // the diagonal local inverse inertia — recomputed independently here. + // Named tolerance: O(1) values over one step, f32 ~1e-5. + const angular_tol: Real = 1e-5; + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const box = try store.createShape(gpa, .{ .box = .{ .half_extents = Vec3.fromArray(.{ 1, 2, 3 }) } }); + + var d = dynDesc(0, box); + d.mass = 6; + d.angular_damping = 0; // isolate the closed form + d.rotation = Quatf.fromAxisAngle(Vec3.fromArray(.{ 0.3, 1, 0.2 }).normalize(), 0.9); + const id = try bm.addBody(gpa, &store, d); + + const torque = vr(5, -3, 2); + bm.addTorque(id, torque); + + // Capture the pre-step pose + inertia the integrator will use. + const rot_before = bm.rotation(id).?; + const local_inv_inertia = bm.motionProperties(id).?.local_inv_inertia; + + const dt: Real = 1.0 / 60.0; + integration.integrate(&bm, dt, vr(0, 0, 0)); // no gravity + + const r = Mat3r.fromQuat(rot_before); + const i_world_inv = r.mul(local_inv_inertia).mul(r.transpose()); + const expected = i_world_inv.mulVec(torque).scale(dt); + try testing.expect(bm.angularVelocity(id).?.approxEql(expected, angular_tol)); +} + +test "angular damping and clamp" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + const dt: Real = 1.0 / 60.0; + const no_gravity = vr(0, 0, 0); + + // Reduce: d > 0, no torque, initial ω ⇒ |ω| strictly smaller after one step. + var damped = dynDesc(0, s); + damped.angular_damping = 0.5; + const id_damped = try bm.addBody(gpa, &store, damped); + bm.setAngularVelocity(id_damped, vr(0, 4, 0)); + integration.integrate(&bm, dt, no_gravity); + const w = bm.angularVelocity(id_damped).?.toArray()[1]; + try testing.expect(w > 0 and w < 4); // damped toward zero, not past it + + // Clamp: d·dt > 1 ⇒ factor max(0, 1−d·dt) = 0 ⇒ ω zeroed, no sign flip. + var clamp = dynDesc(1, s); + clamp.angular_damping = 100; // 100/60 = 1.667 > 1 + const id_clamp = try bm.addBody(gpa, &store, clamp); + bm.setAngularVelocity(id_clamp, vr(1, -2, 3)); + integration.integrate(&bm, dt, no_gravity); + try testing.expect(bm.angularVelocity(id_clamp).?.approxEql(Vec3r.zero, 0)); +} + +test "orientation quaternion stays unit" { + // The first-order update renormalises each step, so |q| stays ≈ 1 even after + // many steps of nonzero spin. Named tolerance: accumulated f32 renormalise + // noise over 120 steps is well under 1e-5. + const unit_tol: Real = 1e-5; + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + var d = dynDesc(0, s); + d.angular_damping = 0; + d.rotation = Quatf.fromAxisAngle(Vec3.unit_z, 0.3); + const id = try bm.addBody(gpa, &store, d); + bm.setAngularVelocity(id, vr(1.5, -2.0, 0.7)); + + const dt: Real = 1.0 / 60.0; + var step: u32 = 0; + while (step < 120) : (step += 1) integration.integrate(&bm, dt, vr(0, 0, 0)); + + try testing.expectApproxEqAbs(@as(Real, 1), quatNorm(bm.rotation(id).?), unit_tol); +} + +test "zero angular velocity leaves orientation unchanged" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + const q0 = Quatf.fromAxisAngle(Vec3.fromArray(.{ 1, 1, 0 }).normalize(), 0.6); + var d = dynDesc(0, s); + d.angular_damping = 0; + d.rotation = q0; + const id = try bm.addBody(gpa, &store, d); + // ω = 0 (default); no torque. The first-order path never divides by |ω|. + + const dt: Real = 1.0 / 60.0; + var step: u32 = 0; + while (step < 60) : (step += 1) integration.integrate(&bm, dt, vr(0, 0, 0)); + + const q = bm.rotation(id).?; + const q0r = Quatr.fromArray(.{ q0.x, q0.y, q0.z, q0.w }); + try testing.expect(q.approxEql(q0r, 1e-6)); // unchanged + // No NaN produced by the ω = 0 path. + for (q.toArray()) |c| try testing.expect(!std.math.isNan(c)); +} + +test "orientation update uses the left (world-space) quaternion product" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{} }); + + // Non-commutative config: an X-axis initial orientation with a Y-axis spin, + // so the left and right quaternion products give distinct results. + var d = dynDesc(0, s); + d.angular_damping = 0; + const q0 = Quatf.fromAxisAngle(Vec3.unit_x, 0.8); + d.rotation = q0; + const id = try bm.addBody(gpa, &store, d); + bm.setAngularVelocity(id, vr(0, 3, 0)); + + const dt: Real = 1.0 / 60.0; + integration.integrate(&bm, dt, vr(0, 0, 0)); // no gravity, no torque + + // Oracle: replicate the first-order formula both ways. The engine uses the + // LEFT product (world-space ω): q ← normalize(q + ½·dt·(ω_quat ⊗ q)). + const w_quat = Quatr{ .x = 0, .y = 3, .z = 0, .w = 0 }; + const q0r = Quatr.fromArray(.{ q0.x, q0.y, q0.z, q0.w }); + const left = q0r.add(w_quat.mul(q0r).scale(0.5 * dt)).normalize(); + const right = q0r.add(q0r.mul(w_quat).scale(0.5 * dt)).normalize(); + + // Discrimination guard: this (q0, ω) is non-commutative, so left ≠ right — + // protects the test's discriminating power if q0/ω are ever changed. + const discrimination_tol: Real = 1e-4; + try testing.expect(!left.approxEql(right, discrimination_tol)); + + // The engine matches the LEFT product, and NOT the right one. + const match_tol: Real = 1e-6; + try testing.expect(bm.rotation(id).?.approxEql(left, match_tol)); + try testing.expect(!bm.rotation(id).?.approxEql(right, discrimination_tol)); +}