From 7db784bb10166490ebe7e519c85e2a7440359b33 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 12:25:17 +0200 Subject: [PATCH 01/23] docs(forge): add M1.1.6 brief --- briefs/M1.1.6-solver-sequential-impulses.md | 156 ++++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 briefs/M1.1.6-solver-sequential-impulses.md diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md new file mode 100644 index 0000000..d12aad2 --- /dev/null +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -0,0 +1,156 @@ +# M1.1.6 — Forge 3D: Sequential Impulses contact solver + +> **Status:** PLANNED +> **Phase:** 1 +> **Branch:** `phase-1/forge/solver-sequential-impulses` +> **Planned tag:** `v0.11.6-solver-sequential-impulses` +> **Dependencies:** M1.1.0–M1.1.5, hotfix `m1.1.3-hf-epa-frame-dependence` (merged, untagged per 2026-07-14 hotfix policy) +> **Opened:** 2026-07-24 +> **Closed:** — + +--- + +# FROZEN SECTION + +*Produced by Claude.ai. Not modifiable by Claude Code outside a Claude.ai round-trip (cf. § Recorded deviations).* + +## Context + +Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the velocity-level Sequential Impulses contact solver (Catto lineage, Box2D/Jolt tradition) — the FIRST real consumer of the M1.1.3/4 contact manifolds and of the frame-stable `feature_id` warm-start contract. Per `engine-physics-forge.md` §1.1/§1.2, the rigid branch is "Sequential Impulses + NGS position solver" and lives in `rigid/`, distinct from the shared `pipeline/`; NGS is M1.1.7. The M1.1.6 velocity solver therefore carries NO positional bias of any kind (no Baumgarte, no split impulse) — its only velocity bias is restitution. Accepted interim artifact until M1.1.7: a resting body keeps a small constant penetration (the solver kills approach velocity; it does not recover existing penetration). + +## Scope + +- **Body material columns (E1).** `Body` gains `friction: Real` and `restitution: Real` SoA columns, populated from `BodyDescriptor` by `addBody` (currently silently dropped). Stale-safe getters `friction(id) ?Real` / `restitution(id) ?Real` (parity with existing getters). Domain guards on the `addBody` path, debug-assert discipline (M1.1.0 `mass > 0` precedent): `friction >= 0` and finite; `restitution` in `[0, 1]` and finite. +- **Integration split (E1).** `pipeline/integration.zig` splits the monolithic pass into `integrateVelocities(bm, dt, gravity)` (gravity + force/torque accumulators + clamped damping + per-tick accumulator clear, ascending slot order) and `integratePositions(bm, dt)` (position from CURRENT velocity + first-order orientation update + normalize). `integrate(bm, dt, gravity)` remains as the exact composition of the two — existing M1.1.5 tests stay green unmodified. Rationale: SI requires solve to sit between the velocity and position passes; solving around the monolithic pass leaves the per-tick `g·dt` residual advancing positions (steady sinking ≈ `g·dt²` per tick). +- **`rigid/` package (E2).** New `src/modules/forge/forge_3d/rigid/` — flat: `contact_constraint.zig`, `contact_cache.zig`, `velocity_solver.zig`, `root.zig` (facade + lazy-analysis pins per `engine-zig-conventions.md` §13). M1.1.7 (NGS) and joints land here as additive sibling files. Import discipline: `foundation`, `weld_forge` (handle types), `../config.zig`, `../body_manager.zig`, `../pipeline/narrowphase/root.zig`. NEVER `broadphase.zig` (candidate pairs are consumed as data) and never any `pipeline/` file beyond the narrowphase facade. +- **Material combine rules (E2).** `combineFriction(a, b) = @sqrt(a * b)` (geometric mean), `combineRestitution(a, b) = @max(a, b)` — named pub functions in `contact_constraint.zig`, Box2D/Jolt convention. +- **Contact constraint setup (E2).** `ContactConstraint` (one per manifold, points `[4]` inline) built from broadphase candidate pairs: unpack `user_data` as `BodyId`, call `BodyManager.collidePair` per pair. Candidate pairs arrive canonical (`min<<32|max`, sorted, deduped — the M1.1.1 `computePairs` contract) and `collidePair` internally re-canonicalizes by `BodyId`; the solver consumes them as-is and asserts the canonical order (the feature_id frame-stability chain is guaranteed end-to-end, never re-derived). Constraints stored in a contiguous array sorted ascending by the packed pair key — the deterministic iteration order (no hash containers anywhere, M1.1.14 discipline). Per-point precompute in `prepare`: lever arms `r_a = p − x_a`, `r_b = p − x_b`, world inverse inertias, normal effective mass, two tangent effective masses (independent per axis, coupling ignored — Jolt convention), and the pre-solve relative normal velocity `v_n⁻` (restitution capture). A constraint whose total inverse mass along the normal is EXACTLY zero (e.g. kinematic×static reaching the pair filter) is skipped at build — a true-zero guard, not a geometric epsilon (M1.1.4 threshold discipline). Deep-band contacts with `penetration ≈ 0` (+ `r_sum`) are legitimate constraints (they forbid approach velocity) and are NEVER rejected on penetration; `penetration` is carried but unconsumed in M1.1.6 (NGS consumes it at M1.1.7). +- **Tangent basis (E2).** Pure function of the world normal, trig-free: branch on the dominant |component| with fixed tie-breaks, `t2 = n × t1`, both normalized. Deterministic; continuous while the normal moves continuously; DISCONTINUOUS when the dominant axis flips — the warm-start tangent reprojection (below) is the compensation, never an assumption of basis stability. +- **Warm-start cache (E3).** `ContactCache`: double-buffered flat arrays (prev/current), key `(pair_key: u64 = min(BodyId)<<32|max, subshape_id: u32 = 0, feature_id: u32)` — `subshape_id` reserved at 0 for the Compound/Mesh co-design flagged at M1.1.3, zero-cost now. Value: accumulated normal impulse `lambda_n: Real` + accumulated tangent impulse as a WORLD-SPACE `Vec3r` (not two basis scalars — basis flips would rotate a scalar pair arbitrarily in the tangent plane). Matching: binary search in the sorted prev buffer, per individual feature key — a point without a match cold-starts; no assumption of manifold count stability across frames (contact topology legitimately changes at band/topology boundaries; RD-5's fixed-order count divergence is one instance, not the contract). `beginTick`/`endTick` exactly once per tick, independent of any future per-range solve calls; `endTick` sorts the current buffer by the full key then swaps buffers — eviction is implicit (unmatched entries do not survive the swap). Hit/miss counters exposed (future `PhysicsDebugProvider.get_warm_start_data` feed). +- **Warm start application (E3).** For each matched point: re-apply `lambda_n · n_new` plus the cached world tangent impulse REPROJECTED onto the new tangent plane — `j_t' = j_t − (j_t·n_new)·n_new`, norm clamped to `mu · lambda_n`, then decomposed onto the new basis to seed the accumulated `lambda_t1/lambda_t2`. Impulses applied to both bodies' velocities before the first iteration. +- **Velocity iterations — normal + restitution (E4).** Gauss-Seidel over constraints in ascending pair-key order, points in manifold order, `velocity_iterations` passes. Normal: accumulated-impulse clamp `lambda_n >= 0` (Catto). Restitution: bias `−e · v_n⁻` applied iff `|v_n⁻| > restitution_threshold`, with `v_n⁻` captured in `prepare` — AFTER `integrateVelocities`, BEFORE `warmStart` (warm-started velocities must never suppress a bounce; pinned by a regression test). +- **Velocity iterations — friction (E5).** Per point, per iteration, AFTER the normal update of that point: independent tangential impulses on the two axes, then a CIRCULAR clamp on the accumulated pair — `‖(lambda_t1, lambda_t2)‖ <= mu · lambda_n` (against the current accumulated normal impulse), by rescaling the accumulated vector. Ratified choice over the box clamp: the box form is anisotropic (up to `√2·mu·lambda_n` on the diagonal) and makes the result depend on the arbitrary basis; the circular clamp is isotropic and basis-independent, coherent with the world-space tangent cache. +- **`SolverConfig` (E4).** `velocity_iterations: u32 = 8` (named for the M1.1.7 `position_iterations` sibling), `restitution_threshold: Real = 1.0` (m/s — a PHYSICAL velocity constant, config field, not a buried literal; not a geometric epsilon). +- **Island-additivity seam (E4).** The iteration entry solves over an explicit constraint index range (`solveRange(from, to)`-shaped); M1.1.6 always passes the full range. M1.1.8 islands become a sort + contiguous ranges — purely additive. +- **Test harness (E6).** Full-pipeline composition IN TESTS ONLY (production `step()` is M1.1.15), normative per-tick cycle frozen here: (1) `Broadphase.computePairs` on current poses → (2) `collidePair` per candidate → (3) `integrateVelocities(dt, gravity)` → (4) solver `beginTick` → build+`prepare` (captures `v_n⁻`) → `warmStart` → iterate ×`velocity_iterations` → (5) `integratePositions(dt)` → (6) solver `endTick` (sort + swap) → (7) broadphase proxy updates. + +## Out of scope + +- **NGS position solver** (M1.1.7) — and consequently ANY positional bias in this milestone: no Baumgarte, no split impulse, no pseudo-velocities. Resting penetration is not recovered here. +- Island manager, sleep detection, parallel solving (M1.1.8). SIMD/batched constraint solving. +- Joints (dedicated milestone). `PhysicsWorld` / `step()` orchestration / `PhysicsModule` instantiation / ECS Transform sync (M1.1.15). +- Raycast / shapecast / overlap / point queries (M1.1.9–10). +- CCD / speculative contacts (`Body.flags.continuous` stays stored, unconsumed — consistent with the M1.1.2/M1.1.5 scope notes). +- The M1.1.4 clip-recentring precision boundary (contact-point position noise far from A's centre). The solver consumes `r = p − x` and inherits, without structurally amplifying, the narrowphase's absolute position noise; the fix (clip-recentring) is upstream in `generateManifold` and stays a future precision hotfix. f64 build covers >10 km worlds meanwhile. +- Per-shape / per-face `PhysicsMaterial` implementation. The PRECEDENCE decision is recorded now (material overrides body properties when present — spec patch at close); the implementation lands with per-shape materials. +- Kinematic position-from-velocity (M1.1.5 additive branch, unchanged). Kinematic bodies participate correctly in the solve already: `inv_mass == 0` (infinite mass) with their velocity read. +- `forge_2d`. Shapes beyond sphere/box/capsule (M1.1.11+). + +## Specs to read first + +1. `engine-physics-forge.md` — §1 (architecture, rigid branch, `rigid/` placement), §2 (body components, materials defaults), §3 (Collision System: three-regime bands as reformulated, manifold contract, `feature_id` frame-stability, order-independence), §4 (Physical Materials). +2. `engine-phase-1-plan.md` — M1.1 table (M1.1.6 line + downstream M1.1.7/8/14/15 boundaries). +3. `engine-zig-conventions.md` — §13 (tests, lazy-analysis guard), naming/formatting. + +## Files to create or modify + +- `src/modules/forge/forge_3d/body.zig` — modify — `friction`/`restitution` fields on `Body`; domain debug asserts documented. +- `src/modules/forge/forge_3d/body_manager.zig` — modify — store both fields in `addBody`; stale-safe getters. +- `src/modules/forge/forge_3d/pipeline/integration.zig` — modify — split into `integrateVelocities`/`integratePositions`; `integrate` kept as composition. +- `src/modules/forge/forge_3d/rigid/contact_constraint.zig` — create — `ContactConstraint`, build/`prepare`, combine rules, tangent basis, true-zero inverse-mass skip. +- `src/modules/forge/forge_3d/rigid/contact_cache.zig` — create — double-buffered sorted cache, key/value types, matching, counters, `beginTick`/`endTick`. +- `src/modules/forge/forge_3d/rigid/velocity_solver.zig` — create — `SolverConfig`, `warmStart` (incl. tangent reprojection), iterations (normal + restitution + friction), range-shaped entry. +- `src/modules/forge/forge_3d/rigid/root.zig` — create — facade re-exports + comptime lazy-analysis pins. +- `src/modules/forge/forge_3d/root.zig` — modify — re-export `rigid` at `Real`. +- `src/modules/forge/forge_3d/tests/solver_test.zig` — create — harness + acceptance scenarios (below). +- `CLAUDE.md` — modify — §3.4 update in the closing PR (patch content authored by Claude.ai during gate reviews). + +## Acceptance criteria + +### Tests + +- `tests/integration_test.zig` — `test "split integration composes to the monolithic pass"` — free-flight trajectory via `integrateVelocities`+`integratePositions` is bit-identical to `integrate`; all existing M1.1.5 tests untouched and green. +- `tests/body_manager_test.zig` — `test "friction and restitution survive addBody"` — descriptor round-trip through the SoA columns and getters. +- `tests/solver_test.zig` — `test "box dropped on static ground comes to rest"` — e=0: approach velocity killed, penetration bounded by its at-impact value across 300 ticks (NO sinking; absolute non-recovery accepted, M1.1.7). +- `tests/solver_test.zig` — `test "restitution bounces above threshold and kills below"` — e=0.8 impact rebounds at ≈ e·|v| (tolerance on solver convergence); an impact below `restitution_threshold` does not bounce. +- `tests/solver_test.zig` — `test "restitution is captured before warm start"` — a warm-started resting contact followed by a real impact still bounces (regression pin for the prepare→warmStart order). +- `tests/solver_test.zig` — `test "inclined static box holds below friction angle and slides above"` — rotated static box (Plane is M1.1.11): `tan θ < μ` holds (drift below a named tolerance), `tan θ > μ` slides. +- `tests/solver_test.zig` — `test "friction clamp is isotropic"` — sliding direction rotated in the tangent plane yields the same friction magnitude (circular clamp, basis-independent). +- `tests/solver_test.zig` — `test "tangent basis flip preserves warm-started friction direction"` — normal rotated across a dominant-axis flip: reprojected world tangent impulse direction is continuous (no arbitrary rotation). +- `tests/solver_test.zig` — `test "zero total inverse mass pair is skipped without NaN"` — kinematic×static contact produces no constraint and no NaN anywhere. +- `tests/solver_test.zig` — `test "warm start hits on a resting contact and misses on BodyId generation reuse"` — hit-rate counter ≥ named threshold on a resting box after N ticks; remove+re-add reusing the slot (bumped generation) yields a cold start. +- `tests/solver_test.zig` — `test "manifold point disappearance and reappearance cold-starts only that point"` — per-feature matching; surviving points keep their impulses. +- `tests/solver_test.zig` — `test "solve is deterministic across identical runs"` — two identical 300-tick runs bit-exact (fixed constraint order, no hash containers). + +### Benchmarks + +- None mandated (M1.1.3/4/5 precedent: absolute perf targets land with the full `step`). + +### Observable behavior + +- Harness scenario (test-runnable): box dropped from 2 m on a static ground box — falls, impacts, rests (e=0) / bounces with decaying apexes (e=0.8). Logged positions demonstrate no post-rest sinking. + +### CI + +- `zig build` clean, zero warnings, configured matrix. +- `zig build test` and `zig build test-forge-3d` green — debug + ReleaseSafe, f32 AND `-Dphysics_f64=true`. +- `zig fmt --check` green; `commit-msg` hook green on every commit. + +## Gates + +- **E1** — Material columns + getters + domain asserts; integration split (RED-first: composition-equivalence test written failing against the monolithic-only API first). +- **E2** — `rigid/` package skeleton + combine rules + tangent basis + constraint build/`prepare` (effective masses, `v_n⁻` capture, true-zero skip, deterministic order). +- **E3** — Contact cache (double-buffer, full key, matching, counters, `beginTick`/`endTick`) + `warmStart` with world-space tangent reprojection. +- **E4** — Velocity iterations: normal accumulated clamp + restitution (threshold, prepare-time capture) + `SolverConfig` + range-shaped entry. +- **E5** — Friction: independent tangent effective masses + circular clamp. +- **E6** — Harness (normative cycle), acceptance scenarios, determinism run, f64 run, CLAUDE.md §3.4 patch, closing notes. + +## Conventions + +- **Branch:** `phase-1/forge/solver-sequential-impulses` +- **Final tag:** `v0.11.6-solver-sequential-impulses` +- **PR title:** `Phase 1 / Forge / Sequential Impulses contact solver` +- **Commit convention:** Conventional Commits (cf. `engine-development-workflow.md` §4.3); squash title `feat(forge): add sequential impulses contact solver (M1.1.6)` +- **Merge strategy:** squash-and-merge (Guy merges and tags; Claude Code never merges or tags) + +## Notes + +- **Inherited hotfix contracts, binding on this design:** (a) deep depth is always a realizable support hyperplane's; (b) `feature_id` frame-stability holds ONLY through `collidePair` (BodyId order) — the solver must never call bare `collide`; (c) manifold count is NOT stable across frames — matching is per individual feature key; (d) deep-band contacts at penetration ≈ 0 (+ `r_sum`) are legitimate and solved normally; (e) S2 MTV ties resolve by a fixed deterministic order — the tangent basis handles normal jumps via reprojection, never by assuming continuity. +- **No positional bias, restated:** any Baumgarte-style term found in review is a STOP. The spec's position-correction answer is NGS at M1.1.7. +- **Threshold discipline (M1.1.4):** no absolute geometric thresholds; guards fire at true zero (`inv_mass_sum == 0`, vector normalization at `floatMin`-class checks only); `restitution_threshold` is a physical velocity config field, exempt by nature. +- **Determinism by construction (anticipates M1.1.14):** sorted flat arrays only; `std.sort` with a total-order key; constraint iteration ascending by pair key; per-manifold point order as produced by `collidePair`. +- **Friction order within an iteration:** normal impulse update first, then the point's tangential update against the CURRENT accumulated `lambda_n` (Jolt order) — fixed and documented. +- **Alternatives examined and rejected:** Baumgarte velocity bias (energy injection; discarded at M1.1.7 → throwaway); box friction clamp (anisotropic, basis-dependent); per-basis scalar tangent cache (rotates arbitrarily on basis flips); hash-map cache (M1.1.14); solver in `pipeline/` (spec §1.2 places SI in `rigid/`; `pipeline/` is the branch-shared layer). +- **Spec patch at close (Claude.ai-produced, KB):** `engine-physics-forge.md` — combine rules + material domains (`friction >= 0` finite, `restitution ∈ [0,1]` finite) + body-vs-material precedence + a §-level Sequential Impulses solver paragraph (cycle, clamps, warm-start semantics, config). + +--- + +# LIVING SECTION + +*Maintained by Claude Code during the milestone.* + +## Specs read + +- [ ] `engine-physics-forge.md` (§1, §2, §3, §4) — read +- [ ] `engine-phase-1-plan.md` (M1.1 table) — read +- [ ] `engine-zig-conventions.md` — read + +## Execution log + +- — + +## Recorded deviations + +- — + +## Blockers encountered + +- — resolved by or + +## Closing notes + +- **What worked:** +- **What deviated from the original spec:** +- **What to flag explicitly in review:** +- **Final measurements:** +- **Residual risks / tech debt left intentionally:** From 6c9002fb23f8cdaad6d27489443feac1a76224ea Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 12:35:25 +0200 Subject: [PATCH 02/23] feat(forge): store body friction and restitution (M1.1.6) --- src/modules/forge/forge_3d/body.zig | 12 ++++++++- src/modules/forge/forge_3d/body_manager.zig | 22 ++++++++++++++++ .../forge_3d/tests/body_manager_test.zig | 25 +++++++++++++++++++ 3 files changed, 58 insertions(+), 1 deletion(-) diff --git a/src/modules/forge/forge_3d/body.zig b/src/modules/forge/forge_3d/body.zig index bf5c68e..350b33b 100644 --- a/src/modules/forge/forge_3d/body.zig +++ b/src/modules/forge/forge_3d/body.zig @@ -3,7 +3,9 @@ //! `MotionProperties` holds the inverse mass + inverse local inertia the solver //! 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. +//! (infinite effective mass). `Body` is the SoA row `BodyManager` stores. The +//! `friction`/`restitution` columns are the per-body material coefficients the +//! Sequential Impulses contact solver reads (M1.1.6). //! //! The `force`/`torque` columns are world-space per-tick accumulators (N, N·m): //! `BodyManager.addForce`/`addTorque` add into them and `integrate` (E2) reads @@ -66,6 +68,14 @@ pub const Body = struct { torque: Vec3r, /// Derived inverse mass/inertia + damping/gravity. motion: MotionProperties, + /// Coulomb friction coefficient (>= 0), stored from the descriptor. Consumed + /// by the Sequential Impulses contact solver (M1.1.6); the M1.1.5 `addBody` + /// dropped it. + friction: Real, + /// Restitution / bounciness in [0, 1], stored from the descriptor. Consumed + /// by the Sequential Impulses contact solver (M1.1.6); the M1.1.5 `addBody` + /// dropped it. + restitution: Real, /// Collision-shape handle (into a `ShapeStore`). shape: ShapeId, /// Simulation class. diff --git a/src/modules/forge/forge_3d/body_manager.zig b/src/modules/forge/forge_3d/body_manager.zig index 31f5cf2..cd30afe 100644 --- a/src/modules/forge/forge_3d/body_manager.zig +++ b/src/modules/forge/forge_3d/body_manager.zig @@ -69,6 +69,12 @@ pub const BodyManager = struct { /// `error.InvalidShape` on a stale/invalid `desc.shape`. pub fn addBody(self: *BodyManager, gpa: std.mem.Allocator, store: *const ShapeStore, desc: BodyDescriptor) !BodyId { const shape = store.get(desc.shape) orelse return error.InvalidShape; + // Material domain guards (debug-only; the M1.1.0 `mass > 0` precedent). + // Friction is a non-negative Coulomb coefficient; restitution is a [0, 1] + // ratio. Both must be finite. Typed-error descriptor validation is a later + // milestone; this guards the otherwise-unchecked material path. + std.debug.assert(std.math.isFinite(desc.friction) and desc.friction >= 0); + std.debug.assert(std.math.isFinite(desc.restitution) and desc.restitution >= 0 and desc.restitution <= 1); const body = Body{ .position = convVec3(desc.position), .rotation = convQuat(desc.rotation), @@ -77,6 +83,8 @@ pub const BodyManager = struct { .force = Vec3r.zero, .torque = Vec3r.zero, .motion = body_mod.computeMotion(desc, shape), + .friction = desc.friction, + .restitution = desc.restitution, .shape = desc.shape, .body_type = desc.body_type, .collision_layer = desc.collision_layer, @@ -124,6 +132,20 @@ pub const BodyManager = struct { return self.bodies.items(.motion)[idx]; } + /// Safe getter: the Coulomb friction coefficient, or null if `id` is + /// stale/invalid. Consumed by the Sequential Impulses contact solver (M1.1.6). + pub fn friction(self: *const BodyManager, id: BodyId) ?Real { + const idx = self.alloc.validate(id) orelse return null; + return self.bodies.items(.friction)[idx]; + } + + /// Safe getter: the restitution / bounciness, or null if `id` is + /// stale/invalid. Consumed by the Sequential Impulses contact solver (M1.1.6). + pub fn restitution(self: *const BodyManager, id: BodyId) ?Real { + const idx = self.alloc.validate(id) orelse return null; + return self.bodies.items(.restitution)[idx]; + } + /// Safe getter: the body's object collision-layer index, or null if `id` is /// stale/invalid. pub fn collisionLayer(self: *const BodyManager, id: BodyId) ?u8 { 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 125d85d..806e255 100644 --- a/src/modules/forge/forge_3d/tests/body_manager_test.zig +++ b/src/modules/forge/forge_3d/tests/body_manager_test.zig @@ -437,3 +437,28 @@ test "isAliveIndex tracks liveness across create, free, and reuse" { try testing.expectEqual(ib, api.PackedId.unpack(c).index); try testing.expect(bm.alloc.isAliveIndex(ib)); } + +test "friction and restitution survive addBody" { + 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 = .{} }); + + // The descriptor material coefficients are dropped by M1.1.5's `addBody`; + // M1.1.6 stores them on `Body` for the Sequential Impulses solver. 0.25 and + // 0.75 are exactly representable in f32 and f64, so the f32→Real widening + // (and the f64 build) both read back the literal exactly. + var desc = descOf(0, .dynamic, s); + desc.friction = 0.25; + desc.restitution = 0.75; + const id = try bm.addBody(gpa, &store, desc); + try testing.expectEqual(@as(Real, 0.25), bm.friction(id).?); + try testing.expectEqual(@as(Real, 0.75), bm.restitution(id).?); + + // Stale handle ⇒ null (parity with the other stale-safe getters). + bm.removeBody(id); + try testing.expectEqual(@as(?Real, null), bm.friction(id)); + try testing.expectEqual(@as(?Real, null), bm.restitution(id)); +} From 166ad2addb5cdc3220c23731baf7ef4c12bac4b3 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 12:40:01 +0200 Subject: [PATCH 03/23] feat(forge): split integration into velocity/position passes (M1.1.6) --- .../forge/forge_3d/pipeline/integration.zig | 126 ++++++++++++------ .../forge/forge_3d/tests/integration_test.zig | 71 ++++++++++ 2 files changed, 154 insertions(+), 43 deletions(-) diff --git a/src/modules/forge/forge_3d/pipeline/integration.zig b/src/modules/forge/forge_3d/pipeline/integration.zig index a81daac..c7275ca 100644 --- a/src/modules/forge/forge_3d/pipeline/integration.zig +++ b/src/modules/forge/forge_3d/pipeline/integration.zig @@ -1,12 +1,23 @@ //! `forge_3d/pipeline/integration.zig` — semi-implicit (symplectic) Euler -//! integration over the live `BodyManager` SoA store. +//! integration over the live `BodyManager` SoA store, split into a velocity pass +//! and a position pass. //! -//! 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 +//! `integrateVelocities` applies gravity + the force/torque accumulators + the +//! clamped damping (and clears the accumulators); `integratePositions` advances +//! position and orientation from the CURRENT velocity. `integrate` is their exact +//! composition. The split exists because the Sequential Impulses contact solver +//! (M1.1.6) must run BETWEEN the two — it corrects velocities after gravity but +//! before positions advance. Solving around the fused pass would leave the +//! per-tick `g·dt` residual advancing positions (a resting body would sink +//! ≈ g·dt² per tick). With no solve in between, the two halves reproduce the +//! fused pass bit-for-bit (pinned by `tests/integration_test.zig`). +//! +//! Each pass is one pure per-tick sweep in ascending slot-index order +//! (deterministic — no hash container on the path, the `BodyManager` discipline; +//! M1.1.14). The store does not compact, so each 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 contact response +//! is the solver's job (M1.1.6) — under gravity alone a dynamic body follows an //! unobstructed trajectory. //! //! Semi-implicit Euler: velocity is integrated first, position from the *new* @@ -15,18 +26,18 @@ //! 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). +//! (`engine-physics-forge.md` §2) by the velocity pass. `integrateVelocities` +//! 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). +//! integrated then clamp-damped (velocity pass), and the orientation advances by +//! the first-order (linearized) rule `q ← normalize(q + ½·dt·(ω_quat ⊗ q))` with +//! the world-space `ω` on the LEFT (position pass). 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"); @@ -37,21 +48,17 @@ 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. +/// Integrate velocities one fixed tick of `dt` seconds under world-space +/// `gravity` (m/s²), in ascending slot-index order, then clear the per-tick +/// force/torque accumulators for every live body. /// /// 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); +/// `v += a·dt`; `v *= max(0, 1 − linear_damping·dt)`. Angular: +/// `α = (R·I_local_inv·Rᵀ)·torque`; `ω += α·dt`; +/// `ω *= max(0, 1 − angular_damping·dt)`. Static and kinematic bodies keep their +/// velocity. Positions and orientations are NOT touched — `integratePositions` +/// advances them (in the full pipeline, after the contact solver has run). +pub fn integrateVelocities(bm: *BodyManager, dt: Real, gravity: Vec3r) void { const rotations = bm.bodies.items(.rotation); const linear_velocities = bm.bodies.items(.linear_velocity); const angular_velocities = bm.bodies.items(.angular_velocity); @@ -74,12 +81,9 @@ pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { // 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 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)); + linear_velocities[i] = v.scale(damp); // --- Angular --- // World-space inverse inertia I_world_inv = R · I_local_inv · Rᵀ @@ -87,16 +91,9 @@ pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { 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 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(); + angular_velocities[i] = w.scale(ang_damp); } // Reset the per-tick accumulators for every live body (§2), including @@ -105,3 +102,46 @@ pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { torques[i] = Vec3r.zero; } } + +/// Advance every live DYNAMIC body's position and orientation one fixed tick of +/// `dt` seconds from its CURRENT velocity, in ascending slot-index order. Call +/// this AFTER `integrateVelocities` (in the full pipeline, after the contact +/// solver has corrected velocities). +/// +/// Linear: `x += linear_velocity·dt`. Angular (first-order, world-space `ω` on +/// the left): `q ← normalize(q + ½·dt·(ω_quat ⊗ q))`. Static and kinematic bodies +/// are not moved (kinematic position-from-velocity is a later additive concern). +pub fn integratePositions(bm: *BodyManager, dt: Real) 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 body_types = bm.bodies.items(.body_type); + + const n: u32 = @intCast(bm.bodies.len); + var i: u32 = 0; + while (i < n) : (i += 1) { + if (!bm.alloc.isAliveIndex(i)) continue; + if (body_types[i] != .dynamic) continue; + + // Position from the current (post-solve) velocity. + positions[i] = positions[i].add(linear_velocities[i].scale(dt)); + + // First-order orientation update: q ← normalize(q + ½·dt·(ω_quat ⊗ q)), + // world-space ω on the LEFT. No |ω| division ⇒ singularity-free at ω = 0. + const wa = angular_velocities[i].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(); + } +} + +/// Advance every live body one fixed tick of `dt` seconds under world-space +/// `gravity` (m/s²): the exact composition of `integrateVelocities` then +/// `integratePositions`, with no contact solve between (free flight). Existing +/// M1.1.5 call sites and tests use this fused form unchanged; the full pipeline +/// (M1.1.15) calls the two halves directly with the solver in between. +pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { + integrateVelocities(bm, dt, gravity); + integratePositions(bm, dt); +} diff --git a/src/modules/forge/forge_3d/tests/integration_test.zig b/src/modules/forge/forge_3d/tests/integration_test.zig index fba8679..54cfe4d 100644 --- a/src/modules/forge/forge_3d/tests/integration_test.zig +++ b/src/modules/forge/forge_3d/tests/integration_test.zig @@ -3,6 +3,10 @@ //! 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. +//! +//! M1.1.6 adds the velocity/position split: `integrate` is now the composition +//! of `integrateVelocities` + `integratePositions`, and the composition is pinned +//! bit-for-bit here so the Sequential Impulses solve can sit between the two. const std = @import("std"); const config = @import("../config.zig"); @@ -541,3 +545,70 @@ test "orientation update uses the left (world-space) quaternion product" { try testing.expect(bm.rotation(id).?.approxEql(left, match_tol)); try testing.expect(!bm.rotation(id).?.approxEql(right, discrimination_tol)); } + +// --- M1.1.6 integration split ------------------------------------------------ + +test "split integration composes to the monolithic pass" { + // The Sequential Impulses solve (M1.1.6) sits BETWEEN the velocity and + // position updates, so `integrate` is split into `integrateVelocities` + + // `integratePositions`. Running the two halves back-to-back (no solve + // between) must reproduce the fused `integrate` BIT-FOR-BIT — the composition + // contract that lets the solver slot in without perturbing free flight. + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + const s = try store.createShape(gpa, .{ .box = .{ .half_extents = Vec3.fromArray(.{ 1, 0.5, 0.25 }) } }); + + // Two identical worlds: `fused` driven by `integrate`, `split` by the halves. + var bm_fused = BodyManager{}; + defer bm_fused.deinit(gpa); + var bm_split = BodyManager{}; + defer bm_split.deinit(gpa); + + // A dynamic, rotated, spinning, damped body exercises every term both passes + // touch: gravity, force·inv_mass, clamped linear + angular damping, the + // R·I·Rᵀ torque map, and the left-product orientation update. + var desc = dynDesc(0, s); + desc.mass = 2; + desc.linear_damping = 0.1; + desc.angular_damping = 0.2; + desc.rotation = Quatf.fromAxisAngle(Vec3.unit_z, 0.5); + const id_fused = try bm_fused.addBody(gpa, &store, desc); + const id_split = try bm_split.addBody(gpa, &store, desc); + bm_fused.setLinearVelocity(id_fused, vr(1, 2, -3)); + bm_split.setLinearVelocity(id_split, vr(1, 2, -3)); + bm_fused.setAngularVelocity(id_fused, vr(0.3, -0.6, 0.9)); + bm_split.setAngularVelocity(id_split, vr(0.3, -0.6, 0.9)); + + const dt: Real = 1.0 / 60.0; + const g = vr(0, -9.81, 0); + var step: u32 = 0; + while (step < 45) : (step += 1) { + // Standing force + torque re-applied each tick (both are cleared per tick + // by the velocity pass, so the two drivers see the same accumulator state). + bm_fused.addForce(id_fused, vr(5, 0, 0)); + bm_split.addForce(id_split, vr(5, 0, 0)); + bm_fused.addTorque(id_fused, vr(0, 0.4, 0)); + bm_split.addTorque(id_split, vr(0, 0.4, 0)); + + integration.integrate(&bm_fused, dt, g); + integration.integrateVelocities(&bm_split, dt, g); + integration.integratePositions(&bm_split, dt); + } + + // Bit-identical position, linear + angular velocity, and orientation. + const pf = bm_fused.position(id_fused).?.toArray(); + const ps = bm_split.position(id_split).?.toArray(); + const vf = bm_fused.linearVelocity(id_fused).?.toArray(); + const vs = bm_split.linearVelocity(id_split).?.toArray(); + const wf = bm_fused.angularVelocity(id_fused).?.toArray(); + const ws = bm_split.angularVelocity(id_split).?.toArray(); + inline for (0..3) |k| { + try testing.expectEqual(pf[k], ps[k]); + try testing.expectEqual(vf[k], vs[k]); + try testing.expectEqual(wf[k], ws[k]); + } + const qf = bm_fused.rotation(id_fused).?.toArray(); + const qs = bm_split.rotation(id_split).?.toArray(); + inline for (0..4) |k| try testing.expectEqual(qf[k], qs[k]); +} From 5eaa1061bcb79a96eba6dc7aa64a88f8063361a4 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 12:41:02 +0200 Subject: [PATCH 04/23] docs(forge): tick specs and log E1 (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index d12aad2..ace3625 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -131,13 +131,18 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel ## Specs read -- [ ] `engine-physics-forge.md` (§1, §2, §3, §4) — read -- [ ] `engine-phase-1-plan.md` (M1.1 table) — read -- [ ] `engine-zig-conventions.md` — read +- [x] `engine-physics-forge.md` (§1, §2, §3, §4) — read 2026-07-24 12:26 +- [x] `engine-phase-1-plan.md` (M1.1 table) — read 2026-07-24 12:26 +- [x] `engine-zig-conventions.md` — read 2026-07-24 12:26 ## Execution log -- — +- 2026-07-24 12:26 — Étape 0/1: read the brief + all three specs in full; branched `phase-1/forge/solver-sequential-impulses` from `main` HEAD `7e63912`; committed the brief verbatim (`7db784b`) and pushed. +- 2026-07-24 12:40 — E1 (material columns + integration split), RED-first: + - Sub-step A (`6c9002f`): observed RED (`no field or member function named 'friction' in BodyManager`) on the round-trip test, then added `friction`/`restitution` `Real` SoA columns on `Body`, stored them in `addBody` behind debug domain guards (`friction >= 0` and finite; `restitution ∈ [0, 1]` and finite — the M1.1.0 `mass > 0` precedent), and added stale-safe `friction(id) ?Real` / `restitution(id) ?Real` getters. GREEN. + - Sub-step B (`166ad2a`): observed RED (`no member named 'integrateVelocities'`) on the composition-equivalence test, then split `integration.zig` into `integrateVelocities` (gravity + force/torque accumulators + clamped damping + per-tick accumulator clear) and `integratePositions` (position from current velocity + first-order orientation + normalize), keeping `integrate` as their exact composition. GREEN. + - Verified: `test-forge-3d` green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. + - Note for review: the Acceptance criteria place the two E1 tests in `tests/integration_test.zig` and `tests/body_manager_test.zig`, which are not enumerated in the brief's "Files to create or modify"; touched per the AC (brief-mandated), flagged here for transparency. ## Recorded deviations From 1bfebf15e7aeec67cf2377d08c1537f18e1038ce Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 14:31:15 +0200 Subject: [PATCH 05/23] docs(forge): fix french token in E1 log (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index ace3625..cc0abf2 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -137,7 +137,7 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel ## Execution log -- 2026-07-24 12:26 — Étape 0/1: read the brief + all three specs in full; branched `phase-1/forge/solver-sequential-impulses` from `main` HEAD `7e63912`; committed the brief verbatim (`7db784b`) and pushed. +- 2026-07-24 12:26 — Step 0/1: read the brief + all three specs in full; branched `phase-1/forge/solver-sequential-impulses` from `main` HEAD `7e63912`; committed the brief verbatim (`7db784b`) and pushed. - 2026-07-24 12:40 — E1 (material columns + integration split), RED-first: - Sub-step A (`6c9002f`): observed RED (`no field or member function named 'friction' in BodyManager`) on the round-trip test, then added `friction`/`restitution` `Real` SoA columns on `Body`, stored them in `addBody` behind debug domain guards (`friction >= 0` and finite; `restitution ∈ [0, 1]` and finite — the M1.1.0 `mass > 0` precedent), and added stale-safe `friction(id) ?Real` / `restitution(id) ?Real` getters. GREEN. - Sub-step B (`166ad2a`): observed RED (`no member named 'integrateVelocities'`) on the composition-equivalence test, then split `integration.zig` into `integrateVelocities` (gravity + force/torque accumulators + clamped damping + per-tick accumulator clear) and `integratePositions` (position from current velocity + first-order orientation + normalize), keeping `integrate` as their exact composition. GREEN. From 9a4ca1f691344fc5aa9bead705a8f8616a5f143e Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 15:22:18 +0200 Subject: [PATCH 06/23] feat(forge): re-export ShapeStore from BodyManager (M1.1.6) --- src/modules/forge/forge_3d/body_manager.zig | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/modules/forge/forge_3d/body_manager.zig b/src/modules/forge/forge_3d/body_manager.zig index cd30afe..c1e9ec5 100644 --- a/src/modules/forge/forge_3d/body_manager.zig +++ b/src/modules/forge/forge_3d/body_manager.zig @@ -37,7 +37,10 @@ const Mat3r = config.Mat3r; const Aabbr = config.Aabbr; const BodyId = api.BodyId; const BodyDescriptor = api.BodyDescriptor; -const ShapeStore = shape_mod.ShapeStore; +/// Generational store of collision shapes. Re-exported so sibling packages (the +/// `rigid/` solver) can name the `collidePair` store parameter type without +/// importing `shape.zig` directly (import-discipline boundary). +pub const ShapeStore = shape_mod.ShapeStore; const Shape = shape_mod.Shape; const Body = body_mod.Body; const MotionProperties = body_mod.MotionProperties; From 03243b66ea90de5ffe7f68a26c926e8ca0834a94 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 15:22:24 +0200 Subject: [PATCH 07/23] feat(forge): add rigid contact constraint setup (M1.1.6) --- .../forge_3d/rigid/contact_constraint.zig | 430 ++++++++++++++++++ src/modules/forge/forge_3d/rigid/root.zig | 32 ++ src/modules/forge/forge_3d/root.zig | 11 + 3 files changed, 473 insertions(+) create mode 100644 src/modules/forge/forge_3d/rigid/contact_constraint.zig create mode 100644 src/modules/forge/forge_3d/rigid/root.zig diff --git a/src/modules/forge/forge_3d/rigid/contact_constraint.zig b/src/modules/forge/forge_3d/rigid/contact_constraint.zig new file mode 100644 index 0000000..d6fe374 --- /dev/null +++ b/src/modules/forge/forge_3d/rigid/contact_constraint.zig @@ -0,0 +1,430 @@ +//! `forge_3d/rigid/contact_constraint.zig` — the velocity-level contact +//! constraint of the rigid-body Sequential Impulses solver (Catto lineage, +//! Box2D/Jolt tradition), plus its build/precompute front end. +//! +//! The rigid branch lives in `rigid/` (distinct from the branch-shared +//! `pipeline/`; `engine-physics-forge.md` §1.2). This file owns: +//! - the material combine rules (`combineFriction` = √(a·b), `combineRestitution` +//! = max(a, b) — the Box2D/Jolt convention); +//! - `tangentBasis`, a trig-free deterministic orthonormal basis for a contact +//! normal; +//! - `ContactConstraint` (one per manifold, ≤ 4 inline points) with its +//! per-point solver data; +//! - `build`, which turns canonical broadphase candidate pairs into a +//! deterministically-ordered constraint array, running the narrowphase and +//! the per-point `prepare` precompute (lever arms, world inverse inertias, +//! normal + two tangent effective masses, pre-solve normal velocity `v_n⁻`). +//! +//! Import discipline (brief): `foundation`, `weld_forge` (handle types), +//! `../config.zig`, `../body_manager.zig`, `../pipeline/narrowphase/root.zig`. +//! NEVER `broadphase.zig` — candidate pairs are consumed as data (packed `u64` +//! keys), never re-derived — and never any `pipeline/` file beyond the +//! narrowphase facade. +//! +//! M1.1.6 is a VELOCITY solver: it carries NO positional bias (no Baumgarte, no +//! split impulse). `ContactPoint.penetration` is carried but unconsumed here; the +//! NGS position solver (M1.1.7) is the spec's position-correction answer. + +const std = @import("std"); +const config = @import("../config.zig"); +const bm_mod = @import("../body_manager.zig"); +const narrowphase = @import("../pipeline/narrowphase/root.zig"); +const api = @import("weld_forge"); +const foundation = @import("foundation"); + +const Real = config.Real; +const Vec3r = config.Vec3r; +const Quatr = config.Quatr; +const Mat3r = config.Mat3r; +const BodyManager = bm_mod.BodyManager; +const ShapeStore = bm_mod.ShapeStore; +const BodyId = api.BodyId; +const ContactManifold = narrowphase.ContactManifold(Real); + +// --- Material combine rules (Box2D/Jolt convention) -------------------------- + +/// Combined friction for a contact pair: the geometric mean √(a·b). Both inputs +/// are non-negative (the `addBody` domain guard), so the product is ≥ 0. +pub fn combineFriction(a: Real, b: Real) Real { + return @sqrt(a * b); +} + +/// Combined restitution for a contact pair: the maximum. A single bouncy surface +/// makes the contact bouncy. +pub fn combineRestitution(a: Real, b: Real) Real { + return @max(a, b); +} + +// --- Tangent basis ----------------------------------------------------------- + +/// An orthonormal tangent basis spanning the plane ⊥ a contact normal. +pub const TangentBasis = struct { + t1: Vec3r, + t2: Vec3r, +}; + +/// Orthonormal, right-handed tangent basis `(t1, t2)` for the unit world normal +/// `n` (so `t1 × t2 == n`). Trig-free and deterministic: branch on the dominant +/// |component| of `n` (fixed tie-breaks x ≥ y ≥ z) to pick a seed world axis +/// guaranteed non-parallel to `n`, Gram-Schmidt it against `n` for `t1`, then +/// `t2 = n × t1`. Continuous while `n` moves continuously, but DISCONTINUOUS +/// across a dominant-axis flip — the warm-start tangent reprojection (E3) +/// compensates, so basis stability is never assumed. +pub fn tangentBasis(n: Vec3r) TangentBasis { + const a = n.toArray(); + const ax = @abs(a[0]); + const ay = @abs(a[1]); + const az = @abs(a[2]); + // Seed = the world axis cyclically after the dominant one, which is + // guaranteed not parallel to `n`, so the Gram-Schmidt residual is nonzero. + var seed: Vec3r = undefined; + if (ax >= ay and ax >= az) { + seed = Vec3r.unit_y; // x dominant + } else if (ay >= az) { + seed = Vec3r.unit_z; // y dominant + } else { + seed = Vec3r.unit_x; // z dominant + } + const t1 = seed.sub(n.scale(seed.dot(n))).normalize(); + const t2 = n.cross(t1).normalize(); + return .{ .t1 = t1, .t2 = t2 }; +} + +// --- Contact constraint ------------------------------------------------------ + +/// One contact point's precomputed solver data. +pub const ConstraintPoint = struct { + /// Lever arm from body A's centre to the contact point (world). + r_a: Vec3r, + /// Lever arm from body B's centre to the contact point (world). + r_b: Vec3r, + /// Normal effective mass = 1 / kₙ (the impulse scale along the normal). + normal_mass: Real, + /// Tangent effective mass on axis 1 = 1 / k_t1. + tangent1_mass: Real, + /// Tangent effective mass on axis 2 = 1 / k_t2. + tangent2_mass: Real, + /// Pre-solve relative normal velocity `v_n⁻` captured in `prepare` (for the + /// restitution bias in E4). Negative when the bodies are approaching. + rel_normal_velocity: Real, + /// Surface penetration along the normal (≥ 0). Carried but UNCONSUMED in + /// M1.1.6 (no positional bias); the NGS position solver reads it at M1.1.7. + penetration: Real, + /// Per-contact feature id — the warm-start matching key (E3). Frame-stable + /// via `BodyManager.collidePair`'s BodyId order. + feature_id: u32, + /// Accumulated normal impulse (warm-started in E3, clamped ≥ 0 in E4). + normal_impulse: Real = 0, + /// Accumulated tangent impulses on the `(t1, t2)` basis (E5). + tangent1_impulse: Real = 0, + tangent2_impulse: Real = 0, +}; + +/// One manifold's worth of contact constraints between a canonical body pair. +/// Points are stored inline (≤ 4); `count` marks the valid entries. +pub const ContactConstraint = struct { + /// Canonical body A (the min BodyId of the pair). + body_a: BodyId, + /// Canonical body B (the max BodyId of the pair). + body_b: BodyId, + /// Packed canonical pair key `min(BodyId)<<32 | max` — the deterministic sort + /// key (ascending iteration order; M1.1.14). + pair_key: u64, + /// World-space contact normal A→B, shared by all points. + normal: Vec3r, + /// Orthonormal tangent basis spanning the contact plane. + tangent1: Vec3r, + tangent2: Vec3r, + /// Combined friction / restitution for the pair. + friction: Real, + restitution: Real, + /// Cached inverse mass and world-space inverse inertia per body — precomputed + /// once in `prepare`, reused by warm start (E3) and the iterations (E4/E5). + inv_mass_a: Real, + inv_mass_b: Real, + inv_inertia_a: Mat3r, + inv_inertia_b: Mat3r, + /// Per-point solver data. + points: [4]ConstraintPoint, + /// Valid entries in `points`, 1..4. + count: u8, +}; + +/// Effective mass along a unit direction `d` at a contact with lever arms `r_a`, +/// `r_b`: `1 / (invMassA + invMassB + (r_a×d)ᵀ Iₐ⁻¹ (r_a×d) + (r_b×d)ᵀ I_b⁻¹ (r_b×d))`. +/// The denominator is strictly positive for any constraint that survives the +/// true-zero skip (invMassA + invMassB > 0), so the inverse is finite. +fn effectiveMass( + inv_mass_a: Real, + inv_mass_b: Real, + inv_inertia_a: Mat3r, + inv_inertia_b: Mat3r, + r_a: Vec3r, + r_b: Vec3r, + d: Vec3r, +) Real { + const rda = r_a.cross(d); + const rdb = r_b.cross(d); + const k = inv_mass_a + inv_mass_b + + rda.dot(inv_inertia_a.mulVec(rda)) + + rdb.dot(inv_inertia_b.mulVec(rdb)); + return 1.0 / k; +} + +/// Build and precompute one contact constraint for the canonical pair `(a, b)` +/// from `manifold`: tangent basis, combined material, cached inverse mass/inertia, +/// and per-point lever arms, effective masses, and pre-solve normal velocity. +fn prepare(bm: *const BodyManager, a: BodyId, b: BodyId, pair_key: u64, manifold: ContactManifold) ?ContactConstraint { + const mp_a = bm.motionProperties(a).?; + const mp_b = bm.motionProperties(b).?; + + // True-zero skip: both bodies have infinite mass (static/kinematic ⇒ inv_mass + // and inv_inertia both zero), so the total inverse mass along the normal is + // exactly zero and there is nothing to solve. A true-zero guard, not a + // geometric epsilon (M1.1.4 threshold discipline). + if (mp_a.inv_mass + mp_b.inv_mass == 0) return null; + + const pos_a = bm.position(a).?; + const pos_b = bm.position(b).?; + const rot_a = bm.rotation(a).?; + const rot_b = bm.rotation(b).?; + const vel_a = bm.linearVelocity(a).?; + const vel_b = bm.linearVelocity(b).?; + const ang_a = bm.angularVelocity(a).?; + const ang_b = bm.angularVelocity(b).?; + + // World-space inverse inertia I_world_inv = R · I_local_inv · Rᵀ per body. + const ra_mat = Mat3r.fromQuat(rot_a); + const inv_inertia_a = ra_mat.mul(mp_a.local_inv_inertia).mul(ra_mat.transpose()); + const rb_mat = Mat3r.fromQuat(rot_b); + const inv_inertia_b = rb_mat.mul(mp_b.local_inv_inertia).mul(rb_mat.transpose()); + + const normal = manifold.normal; + const tb = tangentBasis(normal); + + var c = ContactConstraint{ + .body_a = a, + .body_b = b, + .pair_key = pair_key, + .normal = normal, + .tangent1 = tb.t1, + .tangent2 = tb.t2, + .friction = combineFriction(bm.friction(a).?, bm.friction(b).?), + .restitution = combineRestitution(bm.restitution(a).?, bm.restitution(b).?), + .inv_mass_a = mp_a.inv_mass, + .inv_mass_b = mp_b.inv_mass, + .inv_inertia_a = inv_inertia_a, + .inv_inertia_b = inv_inertia_b, + .points = undefined, + .count = manifold.count, + }; + + for (0..manifold.count) |i| { + const p = manifold.points[i]; + const r_a = p.position.sub(pos_a); + const r_b = p.position.sub(pos_b); + + // v_n⁻: n · (v_b + ω_b×r_b − v_a − ω_a×r_a) — the contact-point relative + // velocity projected onto the normal, captured BEFORE warm start. + const v_pa = vel_a.add(ang_a.cross(r_a)); + const v_pb = vel_b.add(ang_b.cross(r_b)); + + c.points[i] = .{ + .r_a = r_a, + .r_b = r_b, + .normal_mass = effectiveMass(c.inv_mass_a, c.inv_mass_b, inv_inertia_a, inv_inertia_b, r_a, r_b, normal), + .tangent1_mass = effectiveMass(c.inv_mass_a, c.inv_mass_b, inv_inertia_a, inv_inertia_b, r_a, r_b, tb.t1), + .tangent2_mass = effectiveMass(c.inv_mass_a, c.inv_mass_b, inv_inertia_a, inv_inertia_b, r_a, r_b, tb.t2), + .rel_normal_velocity = normal.dot(v_pb.sub(v_pa)), + .penetration = p.penetration, + .feature_id = p.feature_id, + }; + } + return c; +} + +/// Build the contact constraints for `pairs` (canonical packed keys +/// `min(BodyId)<<32 | max`, the M1.1.1 `computePairs` contract: sorted, deduped) +/// into `out`. `out` is cleared first. For each pair the narrowphase runs via +/// `bm.collidePair` (canonical BodyId order → frame-stable feature ids); a +/// non-null manifold becomes one `ContactConstraint` with per-point data +/// precomputed by `prepare`. +pub fn build( + gpa: std.mem.Allocator, + out: *std.ArrayListUnmanaged(ContactConstraint), + bm: *const BodyManager, + store: *const ShapeStore, + pairs: []const u64, +) !void { + out.clearRetainingCapacity(); + for (pairs) |key| { + const a: BodyId = @intCast(key >> 32); + const b: BodyId = @intCast(key & 0xFFFF_FFFF); + // Pairs arrive canonical; the solver asserts the order, never re-derives it. + std.debug.assert(a <= b); + const manifold = bm.collidePair(store, a, b) orelse continue; // separated + const constraint = prepare(bm, a, b, key, manifold) orelse continue; + try out.append(gpa, constraint); + } + // Deterministic iteration order: sort ascending by the packed pair key (a + // total order; the pairs are deduped so keys are unique). No hash containers + // anywhere on the path (M1.1.14). + std.mem.sort(ContactConstraint, out.items, {}, lessByPairKey); +} + +fn lessByPairKey(_: void, x: ContactConstraint, y: ContactConstraint) bool { + return x.pair_key < y.pair_key; +} + +// --- tests ------------------------------------------------------------------- + +const testing = std.testing; + +fn vr(x: Real, y: Real, z: Real) Vec3r { + return Vec3r.fromArray(.{ x, y, z }); +} + +fn descOf(entity_index: u32, body_type: api.BodyType, shape: api.ShapeId) api.BodyDescriptor { + return .{ + .entity = .{ .index = entity_index, .generation = 0 }, + .body_type = body_type, + .shape = shape, + }; +} + +fn pairKey(a: BodyId, b: BodyId) u64 { + return (@as(u64, @min(a, b)) << 32) | @max(a, b); +} + +test "combine rules: friction is the geometric mean, restitution is the max" { + try testing.expectApproxEqAbs(@as(Real, 0.5), combineFriction(0.25, 1.0), 1e-6); // √0.25 + try testing.expectEqual(@as(Real, 0), combineFriction(0, 0.9)); + try testing.expectEqual(@as(Real, 0.8), combineRestitution(0.2, 0.8)); + try testing.expectEqual(@as(Real, 0.8), combineRestitution(0.8, 0.2)); // symmetric +} + +test "tangent basis is orthonormal, right-handed, and deterministic" { + const normals = [_]Vec3r{ + vr(1, 0, 0), vr(0, 1, 0), vr(0, 0, 1), + vr(-1, 0, 0), vr(0, -1, 0), vr(0, 0, -1), + vr(1, 1, 1).normalize(), vr(1, -2, 0.5).normalize(), vr(-3, 0.2, 1).normalize(), + }; + for (normals) |n| { + const tb = tangentBasis(n); + try testing.expectApproxEqAbs(@as(Real, 1), tb.t1.length(), 1e-6); + try testing.expectApproxEqAbs(@as(Real, 1), tb.t2.length(), 1e-6); + try testing.expectApproxEqAbs(@as(Real, 0), tb.t1.dot(n), 1e-6); + try testing.expectApproxEqAbs(@as(Real, 0), tb.t2.dot(n), 1e-6); + try testing.expectApproxEqAbs(@as(Real, 0), tb.t1.dot(tb.t2), 1e-6); + try testing.expect(tb.t1.cross(tb.t2).approxEql(n, 1e-6)); // right-handed + const tb2 = tangentBasis(n); + try testing.expect(tb.t1.approxEql(tb2.t1, 0)); // deterministic + try testing.expect(tb.t2.approxEql(tb2.t2, 0)); + } +} + +test "prepare captures normal effective mass and pre-solve normal velocity" { + 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 } }); + + // Two unit-mass spheres overlapping along +X (centres 0.9 apart, radii sum + // 1.0). For spheres the lever arm r = radius·n ∥ n, so r×n = 0 ⇒ no rotational + // term ⇒ normal_mass = 1/(invMassA + invMassB) = 1/2 exactly. + var da = descOf(0, .dynamic, s); + da.mass = 1; + var db = descOf(1, .dynamic, s); + db.mass = 1; + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + const id_a = try bm.addBody(gpa, &store, da); + const id_b = try bm.addBody(gpa, &store, db); + bm.setLinearVelocity(id_a, vr(1, 0, 0)); + bm.setLinearVelocity(id_b, vr(-1, 0, 0)); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try build(gpa, &constraints, &bm, &store, &.{pairKey(id_a, id_b)}); + + try testing.expectEqual(@as(usize, 1), constraints.items.len); + const c = constraints.items[0]; + try testing.expectEqual(@as(u8, 1), c.count); + try testing.expect(c.normal.approxEql(vr(1, 0, 0), 1e-5)); // A→B is +X + try testing.expectApproxEqAbs(@as(Real, 0.5), c.points[0].normal_mass, 1e-5); + // v_n⁻ = n·(v_b − v_a) = (1,0,0)·(−2,0,0) = −2 (approaching). + try testing.expectApproxEqAbs(@as(Real, -2), c.points[0].rel_normal_velocity, 1e-5); +} + +test "kinematic vs static pair yields no constraint (true-zero inverse mass skip)" { + 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 } }); + + const id_a = try bm.addBody(gpa, &store, descOf(0, .static, s)); + var db = descOf(1, .kinematic, s); + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + const id_b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try build(gpa, &constraints, &bm, &store, &.{pairKey(id_a, id_b)}); + + // Both inv_mass == 0 ⇒ total inverse mass exactly zero ⇒ skipped at build. + try testing.expectEqual(@as(usize, 0), constraints.items.len); +} + +test "separated pair yields no constraint" { + 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 } }); + + const id_a = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + var db = descOf(1, .dynamic, s); + db.position = foundation.math.Vec3.fromArray(.{ 5, 0, 0 }); // far apart + const id_b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try build(gpa, &constraints, &bm, &store, &.{pairKey(id_a, id_b)}); + + try testing.expectEqual(@as(usize, 0), constraints.items.len); +} + +test "constraints come back sorted ascending by pair key" { + 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 } }); + + // Three dynamic spheres, all pairwise distances < 1.0 ⇒ all three pairs + // overlap and produce a constraint. + const d0 = descOf(0, .dynamic, s); + var d1 = descOf(1, .dynamic, s); + d1.position = foundation.math.Vec3.fromArray(.{ 0.8, 0, 0 }); + var d2 = descOf(2, .dynamic, s); + d2.position = foundation.math.Vec3.fromArray(.{ 0.4, 0.6, 0 }); + const b0 = try bm.addBody(gpa, &store, d0); + const b1 = try bm.addBody(gpa, &store, d1); + const b2 = try bm.addBody(gpa, &store, d2); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + // Feed the keys OUT of ascending order — build must sort them. + try build(gpa, &constraints, &bm, &store, &.{ + pairKey(b1, b2), pairKey(b0, b1), pairKey(b0, b2), + }); + + try testing.expectEqual(@as(usize, 3), constraints.items.len); + try testing.expect(constraints.items[0].pair_key < constraints.items[1].pair_key); + try testing.expect(constraints.items[1].pair_key < constraints.items[2].pair_key); +} diff --git a/src/modules/forge/forge_3d/rigid/root.zig b/src/modules/forge/forge_3d/rigid/root.zig new file mode 100644 index 0000000..67112b8 --- /dev/null +++ b/src/modules/forge/forge_3d/rigid/root.zig @@ -0,0 +1,32 @@ +//! `forge_3d/rigid/root.zig` — facade for the rigid-body branch (Sequential +//! Impulses + NGS position solver; `engine-physics-forge.md` §1.2). M1.1.6 lands +//! the velocity-level contact constraint setup (combine rules, tangent basis, +//! `ContactConstraint`, build/prepare); the contact cache (E3) and velocity +//! solver (E4/E5) land here as additive sibling files, and the NGS position +//! solver (M1.1.7) + joints follow the same way. +//! +//! Re-exported at `Real` by `forge_3d/root.zig`. The comptime pin analyses the +//! package's inline tests when built as a test target (engine-zig-conventions.md +//! §13 lazy-analysis guard). + +const contact_constraint = @import("contact_constraint.zig"); + +/// One manifold's velocity-solver contact constraint (≤ 4 inline points). +pub const ContactConstraint = contact_constraint.ContactConstraint; +/// One contact point's precomputed solver data. +pub const ConstraintPoint = contact_constraint.ConstraintPoint; +/// An orthonormal tangent basis for a contact normal. +pub const TangentBasis = contact_constraint.TangentBasis; + +/// Combined pair friction: geometric mean √(a·b). +pub const combineFriction = contact_constraint.combineFriction; +/// Combined pair restitution: max(a, b). +pub const combineRestitution = contact_constraint.combineRestitution; +/// Trig-free deterministic orthonormal tangent basis for a unit normal. +pub const tangentBasis = contact_constraint.tangentBasis; +/// Build a deterministically-ordered constraint array from canonical pairs. +pub const build = contact_constraint.build; + +comptime { + _ = contact_constraint; +} diff --git a/src/modules/forge/forge_3d/root.zig b/src/modules/forge/forge_3d/root.zig index 61f3585..725e77b 100644 --- a/src/modules/forge/forge_3d/root.zig +++ b/src/modules/forge/forge_3d/root.zig @@ -23,6 +23,9 @@ 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"); +// M1.1.6 — rigid-body branch (Sequential Impulses contact solver). Re-exported +// as the `rigid` namespace below; the comptime pin analyses its inline tests. +const rigid_mod = @import("rigid/root.zig"); // --- Solver scalar + math aliases --- @@ -142,6 +145,13 @@ pub fn integrate(bm: *BodyManager, dt: Real, gravity: Vec3r) void { return integration.integrate(bm, dt, gravity); } +// --- Rigid solver (Sequential Impulses contact solver) --- + +/// The rigid-body branch: contact constraint setup + material combine rules + +/// tangent basis (M1.1.6), with the contact cache + velocity solver as additive +/// siblings. Bound to `Real` through the package's `../config.zig` import. +pub const rigid = rigid_mod; + // 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 { @@ -152,6 +162,7 @@ comptime { _ = broadphase; _ = narrowphase; _ = integration; + _ = rigid_mod; _ = @import("tests/body_manager_test.zig"); _ = @import("tests/integration_test.zig"); _ = @import("tests/broadphase_test.zig"); From baa320777fc5299243fcdcd6ef0cbe62ef0f1145 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 15:23:01 +0200 Subject: [PATCH 08/23] docs(forge): log E2 (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index cc0abf2..edad8d9 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -143,6 +143,13 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - Sub-step B (`166ad2a`): observed RED (`no member named 'integrateVelocities'`) on the composition-equivalence test, then split `integration.zig` into `integrateVelocities` (gravity + force/torque accumulators + clamped damping + per-tick accumulator clear) and `integratePositions` (position from current velocity + first-order orientation + normalize), keeping `integrate` as their exact composition. GREEN. - Verified: `test-forge-3d` green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - Note for review: the Acceptance criteria place the two E1 tests in `tests/integration_test.zig` and `tests/body_manager_test.zig`, which are not enumerated in the brief's "Files to create or modify"; touched per the AC (brief-mandated), flagged here for transparency. +- 2026-07-24 15:22 — E2 (rigid package + constraint setup), RED-first: + - `9a4ca1f`: `pub const ShapeStore` re-export on `BodyManager` so the `rigid/` package can name the `collidePair` store type without importing `shape.zig` (staying inside the brief's import discipline). + - `03243b6`: new `rigid/` package — `contact_constraint.zig` (combine rules `combineFriction` = √(a·b) / `combineRestitution` = max; trig-free deterministic `tangentBasis`; `ContactConstraint` + `ConstraintPoint`; `build` from canonical `u64` pair keys via `collidePair`; `prepare` per-point precompute: lever arms, world inverse inertias, normal + two tangent effective masses, `v_n⁻` capture) + `rigid/root.zig` facade with the §13 comptime pin. `forge_3d/root.zig` re-exports the `rigid` namespace at `Real` and pins it. + - RED observed (135 total, 2 fail) BEFORE adding: the true-zero inverse-mass skip (`kinematic vs static` test expected 0, found 1) and the ascending-pair-key sort (`sorted ascending` test failed); both GREEN after adding the `inv_mass_a + inv_mass_b == 0` skip in `prepare` and `std.mem.sort` in `build`. + - Import discipline honoured: no `broadphase.zig` / no `pipeline/` beyond the narrowphase facade; candidate pairs consumed as data (packed `u64` keys), never re-derived; canonical order asserted (`a <= b`). + - Verified: 135/135 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. + - Notes for review: (1) `contact_cache.zig` (E3) and `velocity_solver.zig` (E4/E5) are NOT created yet — they land as additive sibling files at their gates (the brief's own "additive sibling files" model), avoiding empty stubs / dead public surface; (2) `forge_3d/root.zig` exposes `rigid` as a namespace (`pub const rigid`) rather than flattening each symbol, matching "re-export `rigid`" and scaling as the package grows; (3) `ContactConstraint` already carries the accumulated-impulse fields (default 0) and cached inverse mass/inertia — consumed by warm start (E3) and the iterations (E4/E5). ## Recorded deviations From 61242f929ac1d5cf76f61067da19220f5a38b89d Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 15:55:21 +0200 Subject: [PATCH 09/23] feat(forge): add warm-start cache and warm start (M1.1.6) --- .../forge/forge_3d/rigid/contact_cache.zig | 185 +++++++++++++ src/modules/forge/forge_3d/rigid/root.zig | 16 ++ .../forge/forge_3d/rigid/velocity_solver.zig | 243 ++++++++++++++++++ 3 files changed, 444 insertions(+) create mode 100644 src/modules/forge/forge_3d/rigid/contact_cache.zig create mode 100644 src/modules/forge/forge_3d/rigid/velocity_solver.zig diff --git a/src/modules/forge/forge_3d/rigid/contact_cache.zig b/src/modules/forge/forge_3d/rigid/contact_cache.zig new file mode 100644 index 0000000..3f73acc --- /dev/null +++ b/src/modules/forge/forge_3d/rigid/contact_cache.zig @@ -0,0 +1,185 @@ +//! `forge_3d/rigid/contact_cache.zig` — the Sequential Impulses warm-start cache. +//! +//! A per-tick, double-buffered flat cache of accumulated contact impulses keyed +//! by contact feature. The solver seeds each new tick's constraints from the +//! previous tick's solved impulses (warm start; `velocity_solver.warmStart`), +//! which is what makes the iterative solver converge in a handful of iterations. +//! +//! Determinism by construction (M1.1.14): sorted flat arrays only, NO hash +//! containers. The `prev` buffer is sorted ascending by the full key; matching is +//! a binary search per individual feature key. Contact topology legitimately +//! changes between frames (a manifold gaining/losing points, an MTV-tie count +//! divergence), so matching is per-key — a point without a match cold-starts; the +//! manifold count is never assumed stable. +//! +//! Lifecycle, exactly once per tick (independent of any per-range solve calls): +//! - `beginTick` — clear the current buffer, reset the hit/miss counters. +//! - `store` — record a solved contact impulse into the current buffer +//! (the solver harvests every constraint point after the +//! velocity iterations). +//! - `endTick` — sort the current buffer by the full key, then swap it with +//! `prev`. Eviction is implicit: entries in the old `prev` that +//! were not re-stored do not survive the swap. +//! +//! The value's tangent impulse is stored in WORLD space (a `Vec3r`), NOT as two +//! basis scalars: the tangent basis flips discontinuously across a dominant-axis +//! change, so basis scalars would rotate arbitrarily in the tangent plane. A +//! world vector is reprojected onto the new basis at warm start instead. + +const std = @import("std"); +const config = @import("../config.zig"); + +const Real = config.Real; +const Vec3r = config.Vec3r; + +/// A contact's warm-start identity. `subshape_id` is reserved at 0 for the +/// Compound/Mesh sub-shape co-design flagged at M1.1.3 (zero-cost now); the full +/// triple is the sort/match key so it extends without a format change. +pub const CacheKey = struct { + /// Packed canonical body pair `min(BodyId)<<32 | max`. + pair_key: u64, + /// Sub-shape index within a Compound/Mesh pair — always 0 in M1.1.6. + subshape_id: u32 = 0, + /// Per-contact feature id (from the manifold), unique within a manifold. + feature_id: u32, +}; + +/// The accumulated impulses carried between ticks for one contact feature. +pub const CacheValue = struct { + /// Accumulated normal impulse λₙ. + lambda_n: Real, + /// Accumulated tangent impulse in WORLD space (reprojected at warm start). + tangent_impulse: Vec3r, +}; + +const Entry = struct { + key: CacheKey, + value: CacheValue, +}; + +/// Double-buffered warm-start cache. `prev` (sorted) is read at warm start; +/// `current` collects this tick's solved impulses and becomes `prev` at `endTick`. +pub const ContactCache = struct { + /// Last tick's impulses, sorted ascending by the full key (read at warm start). + prev: std.ArrayListUnmanaged(Entry) = .empty, + /// This tick's impulses, filled by `store`, sorted + promoted at `endTick`. + current: std.ArrayListUnmanaged(Entry) = .empty, + /// Warm-start matches since the last `beginTick` (debug/telemetry feed). + hits: u32 = 0, + /// Warm-start cold-starts since the last `beginTick`. + misses: u32 = 0, + + /// Release both buffers. + pub fn deinit(self: *ContactCache, gpa: std.mem.Allocator) void { + self.prev.deinit(gpa); + self.current.deinit(gpa); + self.* = undefined; + } + + /// Start a tick: clear the current buffer (retaining capacity) and reset the + /// hit/miss counters. `prev` (last tick's sorted impulses) is untouched. + pub fn beginTick(self: *ContactCache) void { + self.current.clearRetainingCapacity(); + self.hits = 0; + self.misses = 0; + } + + /// Record a solved contact impulse into the current buffer. + pub fn store(self: *ContactCache, gpa: std.mem.Allocator, key: CacheKey, value: CacheValue) !void { + try self.current.append(gpa, .{ .key = key, .value = value }); + } + + /// Look `key` up in the sorted `prev` buffer (binary search). Returns the + /// cached value on a hit (and counts it) or null on a cold start (miss). + pub fn lookup(self: *ContactCache, key: CacheKey) ?CacheValue { + var lo: usize = 0; + var hi: usize = self.prev.items.len; + while (lo < hi) { + const mid = lo + (hi - lo) / 2; + switch (orderKey(self.prev.items[mid].key, key)) { + .lt => lo = mid + 1, + .gt => hi = mid, + .eq => { + self.hits += 1; + return self.prev.items[mid].value; + }, + } + } + self.misses += 1; + return null; + } + + /// End a tick: sort the current buffer by the full key, then swap it into + /// `prev`. The old `prev` becomes the scratch buffer that `beginTick` clears; + /// entries not re-stored this tick are implicitly evicted. + pub fn endTick(self: *ContactCache) void { + std.mem.sort(Entry, self.current.items, {}, lessByKey); + std.mem.swap(std.ArrayListUnmanaged(Entry), &self.prev, &self.current); + } +}; + +/// Total order over keys: `pair_key`, then `subshape_id`, then `feature_id`. +fn orderKey(a: CacheKey, b: CacheKey) std.math.Order { + if (a.pair_key != b.pair_key) return std.math.order(a.pair_key, b.pair_key); + if (a.subshape_id != b.subshape_id) return std.math.order(a.subshape_id, b.subshape_id); + return std.math.order(a.feature_id, b.feature_id); +} + +fn lessByKey(_: void, x: Entry, y: Entry) bool { + return orderKey(x.key, y.key) == .lt; +} + +// --- tests ------------------------------------------------------------------- + +const testing = std.testing; + +test "lookup finds stored keys after endTick and records hits and misses" { + const gpa = testing.allocator; + var cache = ContactCache{}; + defer cache.deinit(gpa); + + // Store three entries out of key order, then finalize (sort + promote to prev). + try cache.store(gpa, .{ .pair_key = 5, .feature_id = 2 }, .{ .lambda_n = 2, .tangent_impulse = Vec3r.zero }); + try cache.store(gpa, .{ .pair_key = 1, .feature_id = 9 }, .{ .lambda_n = 1, .tangent_impulse = Vec3r.zero }); + try cache.store(gpa, .{ .pair_key = 5, .feature_id = 1 }, .{ .lambda_n = 3, .tangent_impulse = Vec3r.zero }); + cache.endTick(); + + try testing.expectEqual(@as(Real, 1), cache.lookup(.{ .pair_key = 1, .feature_id = 9 }).?.lambda_n); + try testing.expectEqual(@as(Real, 3), cache.lookup(.{ .pair_key = 5, .feature_id = 1 }).?.lambda_n); + try testing.expectEqual(@as(Real, 2), cache.lookup(.{ .pair_key = 5, .feature_id = 2 }).?.lambda_n); + try testing.expect(cache.lookup(.{ .pair_key = 9, .feature_id = 0 }) == null); + try testing.expectEqual(@as(u32, 3), cache.hits); + try testing.expectEqual(@as(u32, 1), cache.misses); +} + +test "endTick swaps buffers and evicts entries not re-stored" { + const gpa = testing.allocator; + var cache = ContactCache{}; + defer cache.deinit(gpa); + + // Tick 1: two contacts persist. + try cache.store(gpa, .{ .pair_key = 1, .feature_id = 1 }, .{ .lambda_n = 1, .tangent_impulse = Vec3r.zero }); + try cache.store(gpa, .{ .pair_key = 2, .feature_id = 1 }, .{ .lambda_n = 2, .tangent_impulse = Vec3r.zero }); + cache.endTick(); + + // Tick 2: only the first is re-stored — the second must be evicted by the swap. + cache.beginTick(); + try cache.store(gpa, .{ .pair_key = 1, .feature_id = 1 }, .{ .lambda_n = 10, .tangent_impulse = Vec3r.zero }); + cache.endTick(); + + try testing.expectEqual(@as(Real, 10), cache.lookup(.{ .pair_key = 1, .feature_id = 1 }).?.lambda_n); + try testing.expect(cache.lookup(.{ .pair_key = 2, .feature_id = 1 }) == null); +} + +test "lookup distinguishes feature_id within one pair_key" { + const gpa = testing.allocator; + var cache = ContactCache{}; + defer cache.deinit(gpa); + + try cache.store(gpa, .{ .pair_key = 7, .feature_id = 20 }, .{ .lambda_n = 20, .tangent_impulse = Vec3r.zero }); + try cache.store(gpa, .{ .pair_key = 7, .feature_id = 10 }, .{ .lambda_n = 10, .tangent_impulse = Vec3r.zero }); + cache.endTick(); + + try testing.expectEqual(@as(Real, 10), cache.lookup(.{ .pair_key = 7, .feature_id = 10 }).?.lambda_n); + try testing.expectEqual(@as(Real, 20), cache.lookup(.{ .pair_key = 7, .feature_id = 20 }).?.lambda_n); +} diff --git a/src/modules/forge/forge_3d/rigid/root.zig b/src/modules/forge/forge_3d/rigid/root.zig index 67112b8..febef9c 100644 --- a/src/modules/forge/forge_3d/rigid/root.zig +++ b/src/modules/forge/forge_3d/rigid/root.zig @@ -10,6 +10,8 @@ //! §13 lazy-analysis guard). const contact_constraint = @import("contact_constraint.zig"); +const contact_cache = @import("contact_cache.zig"); +const velocity_solver = @import("velocity_solver.zig"); /// One manifold's velocity-solver contact constraint (≤ 4 inline points). pub const ContactConstraint = contact_constraint.ContactConstraint; @@ -27,6 +29,20 @@ pub const tangentBasis = contact_constraint.tangentBasis; /// Build a deterministically-ordered constraint array from canonical pairs. pub const build = contact_constraint.build; +/// The Sequential Impulses warm-start cache (double-buffered, sorted flat). +pub const ContactCache = contact_cache.ContactCache; +/// A warm-start cache key (pair + sub-shape + feature). +pub const CacheKey = contact_cache.CacheKey; +/// A warm-start cache value (accumulated normal + world tangent impulse). +pub const CacheValue = contact_cache.CacheValue; + +/// Warm-start constraints from the previous tick's cache (tangent reprojection). +pub const warmStart = velocity_solver.warmStart; +/// Harvest solved constraint impulses into the cache's current buffer. +pub const storeContacts = velocity_solver.storeContacts; + comptime { _ = contact_constraint; + _ = contact_cache; + _ = velocity_solver; } diff --git a/src/modules/forge/forge_3d/rigid/velocity_solver.zig b/src/modules/forge/forge_3d/rigid/velocity_solver.zig new file mode 100644 index 0000000..8278419 --- /dev/null +++ b/src/modules/forge/forge_3d/rigid/velocity_solver.zig @@ -0,0 +1,243 @@ +//! `forge_3d/rigid/velocity_solver.zig` — the Sequential Impulses velocity +//! solver (Catto lineage). M1.1.6 lands warm starting here; the velocity +//! iterations (normal + restitution, then friction) and the `SolverConfig` + +//! range-shaped entry land at E4/E5 as additive changes to this same file. +//! +//! This is a VELOCITY solver only: no positional bias (no Baumgarte, no split +//! impulse). Warm start re-applies the previous tick's accumulated impulses +//! (from `ContactCache`) before the first iteration, which is what lets the +//! iterative solve converge in a handful of passes. +//! +//! Import discipline (brief): `foundation`, `weld_forge`, `../config.zig`, +//! `../body_manager.zig`, sibling `contact_constraint.zig`/`contact_cache.zig`. +//! Velocities are read/written through `BodyManager`'s stale-safe getters/setters +//! (immediate Gauss-Seidel propagation across contacts and points). + +const std = @import("std"); +const config = @import("../config.zig"); +const bm_mod = @import("../body_manager.zig"); +const cc = @import("contact_constraint.zig"); +const cache_mod = @import("contact_cache.zig"); + +const Real = config.Real; +const Vec3r = config.Vec3r; +const BodyManager = bm_mod.BodyManager; +const ContactConstraint = cc.ContactConstraint; +const ContactCache = cache_mod.ContactCache; + +/// Apply a world-space impulse `p` at the contact levers `r_a`/`r_b` to both +/// bodies' velocities: body A receives −p, body B receives +p (the normal is +/// A→B). Naturally a no-op on an infinite-mass body (inv_mass and inv_inertia +/// are both zero for static/kinematic), so no `body_type` branch is needed. +fn applyImpulse(bm: *BodyManager, c: *const ContactConstraint, r_a: Vec3r, r_b: Vec3r, p: Vec3r) void { + bm.setLinearVelocity(c.body_a, bm.linearVelocity(c.body_a).?.sub(p.scale(c.inv_mass_a))); + bm.setAngularVelocity(c.body_a, bm.angularVelocity(c.body_a).?.sub(c.inv_inertia_a.mulVec(r_a.cross(p)))); + bm.setLinearVelocity(c.body_b, bm.linearVelocity(c.body_b).?.add(p.scale(c.inv_mass_b))); + bm.setAngularVelocity(c.body_b, bm.angularVelocity(c.body_b).?.add(c.inv_inertia_b.mulVec(r_b.cross(p)))); +} + +/// Warm-start every constraint from the previous tick's cache. For each point +/// that matches (per full feature key), re-apply the cached normal impulse along +/// the NEW normal plus the cached WORLD tangent impulse reprojected onto the new +/// tangent plane and clamped to μ·λₙ; seed the accumulated impulses and apply +/// them to both bodies' velocities before the first iteration. A per-feature miss +/// cold-starts (impulses left at 0, nothing applied). +pub fn warmStart(bm: *BodyManager, cache: *ContactCache, constraints: []ContactConstraint) void { + for (constraints) |*c| { + for (0..c.count) |i| { + const pt = &c.points[i]; + const cached = cache.lookup(.{ .pair_key = c.pair_key, .feature_id = pt.feature_id }) orelse continue; + const lambda_n = cached.lambda_n; + + // Reproject the cached WORLD tangent onto the new tangent plane + // (remove the component along the new normal — the basis may have + // rotated since it was cached), then circular-clamp its norm to μ·λₙ. + var j_t = cached.tangent_impulse.sub(c.normal.scale(cached.tangent_impulse.dot(c.normal))); + const max_t = c.friction * lambda_n; + const len_sq = j_t.lengthSq(); + if (len_sq > max_t * max_t) { + j_t = j_t.scale(max_t / @sqrt(len_sq)); + } + + // Seed the accumulated impulses (decompose the tangent onto the basis) + // and apply the total impulse to both bodies. + pt.normal_impulse = lambda_n; + pt.tangent1_impulse = j_t.dot(c.tangent1); + pt.tangent2_impulse = j_t.dot(c.tangent2); + const impulse = c.normal.scale(lambda_n).add(j_t); + applyImpulse(bm, c, pt.r_a, pt.r_b, impulse); + } + } +} + +/// Harvest each constraint point's solved impulse into the cache's current buffer +/// (key `(pair_key, 0, feature_id)`, tangent = λ_t1·t1 + λ_t2·t2 in world space). +/// Called after the velocity iterations, before `ContactCache.endTick` sorts and +/// swaps. +pub fn storeContacts(gpa: std.mem.Allocator, cache: *ContactCache, constraints: []const ContactConstraint) !void { + for (constraints) |c| { + for (0..c.count) |i| { + const pt = c.points[i]; + const tangent = c.tangent1.scale(pt.tangent1_impulse).add(c.tangent2.scale(pt.tangent2_impulse)); + try cache.store( + gpa, + .{ .pair_key = c.pair_key, .feature_id = pt.feature_id }, + .{ .lambda_n = pt.normal_impulse, .tangent_impulse = tangent }, + ); + } + } +} + +// --- tests ------------------------------------------------------------------- + +const api = @import("weld_forge"); +const foundation = @import("foundation"); +const ShapeStore = bm_mod.ShapeStore; +const testing = std.testing; + +fn vr(x: Real, y: Real, z: Real) Vec3r { + return Vec3r.fromArray(.{ x, y, z }); +} + +fn descOf(entity_index: u32, body_type: api.BodyType, shape: api.ShapeId) api.BodyDescriptor { + return .{ + .entity = .{ .index = entity_index, .generation = 0 }, + .body_type = body_type, + .shape = shape, + }; +} + +fn pairKey(a: api.BodyId, b: api.BodyId) u64 { + return (@as(u64, @min(a, b)) << 32) | @max(a, b); +} + +/// Two overlapping unit-mass spheres along +X (centres 0.9 apart, radii sum 1.0). +fn twoOverlappingSpheres(gpa: std.mem.Allocator, store: *ShapeStore, bm: *BodyManager) ![2]api.BodyId { + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + var da = descOf(0, .dynamic, s); + da.mass = 1; + var db = descOf(1, .dynamic, s); + db.mass = 1; + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + const id_a = try bm.addBody(gpa, store, da); + const id_b = try bm.addBody(gpa, store, db); + return .{ id_a, id_b }; +} + +test "warm start applies the cached normal impulse and records a hit" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try twoOverlappingSpheres(gpa, &store, &bm); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + try testing.expectEqual(@as(usize, 1), constraints.items.len); + + // Seed the cache with a known normal impulse for this contact's feature. + var cache = ContactCache{}; + defer cache.deinit(gpa); + try cache.store(gpa, .{ + .pair_key = constraints.items[0].pair_key, + .feature_id = constraints.items[0].points[0].feature_id, + }, .{ .lambda_n = 5, .tangent_impulse = Vec3r.zero }); + cache.endTick(); + + cache.beginTick(); + warmStart(&bm, &cache, constraints.items); + + // λₙ = 5 seeded; applied along n = +X ⇒ A gets −5X (inv_mass 1), B gets +5X. + try testing.expectApproxEqAbs(@as(Real, 5), constraints.items[0].points[0].normal_impulse, 1e-5); + try testing.expectApproxEqAbs(@as(Real, -5), bm.linearVelocity(ids[0]).?.toArray()[0], 1e-5); + try testing.expectApproxEqAbs(@as(Real, 5), bm.linearVelocity(ids[1]).?.toArray()[0], 1e-5); + try testing.expectEqual(@as(u32, 1), cache.hits); +} + +test "warm start cold-starts and records a miss on an unknown key" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try twoOverlappingSpheres(gpa, &store, &bm); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + + var cache = ContactCache{}; // empty ⇒ every lookup misses + defer cache.deinit(gpa); + cache.beginTick(); + warmStart(&bm, &cache, constraints.items); + + try testing.expectEqual(@as(Real, 0), constraints.items[0].points[0].normal_impulse); + try testing.expect(bm.linearVelocity(ids[0]).?.approxEql(Vec3r.zero, 0)); // unchanged + try testing.expectEqual(@as(u32, 1), cache.misses); +} + +test "warm start reprojects the cached world tangent onto the new normal plane" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try twoOverlappingSpheres(gpa, &store, &bm); // both friction default 0.5 + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + + // The contact normal is +X. Inject a cached world tangent with a LARGE + // component along the normal (X) and a small in-plane component (Y): the raw + // magnitude far exceeds μ·λₙ = 0.5, but the IN-PLANE magnitude 0.3 is below + // it. Reprojecting first clamps the in-plane part (0.3 < 0.5 ⇒ preserved); + // clamping the RAW magnitude would wrongly shrink the in-plane part. + var cache = ContactCache{}; + defer cache.deinit(gpa); + try cache.store(gpa, .{ + .pair_key = constraints.items[0].pair_key, + .feature_id = constraints.items[0].points[0].feature_id, + }, .{ .lambda_n = 1, .tangent_impulse = vr(4.0, 0.3, 0.0) }); + cache.endTick(); + + cache.beginTick(); + warmStart(&bm, &cache, constraints.items); + + // t1 = +Y, t2 = +Z for n = +X ⇒ the in-plane part decomposes to λ_t1 = 0.3. + try testing.expectApproxEqAbs(@as(Real, 0.3), constraints.items[0].points[0].tangent1_impulse, 1e-5); + try testing.expectApproxEqAbs(@as(Real, 0), constraints.items[0].points[0].tangent2_impulse, 1e-5); +} + +test "solved impulses round-trip through the cache to the next tick's warm start" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try twoOverlappingSpheres(gpa, &store, &bm); + const key = pairKey(ids[0], ids[1]); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + var cache = ContactCache{}; + defer cache.deinit(gpa); + + // Tick 1: build, simulate a solve (λₙ = 7), harvest into the cache, finalize. + cache.beginTick(); + try cc.build(gpa, &constraints, &bm, &store, &.{key}); + constraints.items[0].points[0].normal_impulse = 7; + try storeContacts(gpa, &cache, constraints.items); + cache.endTick(); + + // Tick 2: rebuild (impulses reset to 0), warm start seeds from tick 1. + cache.beginTick(); + try cc.build(gpa, &constraints, &bm, &store, &.{key}); + try testing.expectEqual(@as(Real, 0), constraints.items[0].points[0].normal_impulse); // fresh + warmStart(&bm, &cache, constraints.items); + + try testing.expectApproxEqAbs(@as(Real, 7), constraints.items[0].points[0].normal_impulse, 1e-5); + try testing.expectEqual(@as(u32, 1), cache.hits); +} From 32775fd1ba1867ce3c34a9df5b1c81236641b555 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 15:55:55 +0200 Subject: [PATCH 10/23] docs(forge): log E3 (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index edad8d9..ce8a0b6 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -150,6 +150,13 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - Import discipline honoured: no `broadphase.zig` / no `pipeline/` beyond the narrowphase facade; candidate pairs consumed as data (packed `u64` keys), never re-derived; canonical order asserted (`a <= b`). - Verified: 135/135 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - Notes for review: (1) `contact_cache.zig` (E3) and `velocity_solver.zig` (E4/E5) are NOT created yet — they land as additive sibling files at their gates (the brief's own "additive sibling files" model), avoiding empty stubs / dead public surface; (2) `forge_3d/root.zig` exposes `rigid` as a namespace (`pub const rigid`) rather than flattening each symbol, matching "re-export `rigid`" and scaling as the package grows; (3) `ContactConstraint` already carries the accumulated-impulse fields (default 0) and cached inverse mass/inertia — consumed by warm start (E3) and the iterations (E4/E5). +- 2026-07-24 15:55 — E3 (contact cache + warm start), RED-first (`61242f9`): + - `contact_cache.zig`: `ContactCache` — double-buffered flat arrays (`prev` sorted / `current`), full key `{pair_key, subshape_id=0, feature_id}`, `beginTick` (clear current + reset counters), `store`, `lookup` (binary search of `prev`, per-feature; counts hits/misses), `endTick` (`std.mem.sort` current + `std.mem.swap` into prev — implicit eviction). Value = `{lambda_n, tangent_impulse: Vec3r}` WORLD-space (not basis scalars). No hash containers. + - `velocity_solver.zig`: `warmStart` (per matched point: re-apply cached `lambda_n·n_new` + cached world tangent REPROJECTED onto the new plane `j_t − (j_t·n)·n`, circular-clamped to `μ·λₙ`, decomposed onto the new basis, applied to both bodies) + `storeContacts` (harvest solved impulses → `cache.store`) + `applyImpulse` (via `BodyManager` getters/setters; natural no-op on infinite mass). `warmStart`/iterations land in this file. + - RED observed: the tangent-reprojection test failed under a RAW-magnitude clamp (`tangent1_impulse` = 0.0374, expected 0.3) BEFORE reprojecting-then-clamping the in-plane magnitude; GREEN after the reproject-before-clamp. The world-space cache + reprojection is what preserves the tangent's world direction across a basis flip (vs the rejected per-basis-scalar cache). + - `rigid/root.zig` re-exports `ContactCache`/`CacheKey`/`CacheValue`/`warmStart`/`storeContacts` + pins both new files. + - Verified: 142/142 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. + - Note for review: `endTick` is `sort + swap` per the brief; the harvest (`storeContacts`) is a distinct step run after the iterations and before `endTick` (the current buffer must be filled before it is sorted) — the E6 normative cycle's implicit fill step. ## Recorded deviations From 56c2878b67105d7834f948c977e5d6f71b4106fb Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:11:50 +0200 Subject: [PATCH 11/23] feat(forge): add velocity iterations with restitution (M1.1.6) --- .../forge/forge_3d/rigid/velocity_solver.zig | 195 ++++++++++++++++++ 1 file changed, 195 insertions(+) diff --git a/src/modules/forge/forge_3d/rigid/velocity_solver.zig b/src/modules/forge/forge_3d/rigid/velocity_solver.zig index 8278419..bd6fec2 100644 --- a/src/modules/forge/forge_3d/rigid/velocity_solver.zig +++ b/src/modules/forge/forge_3d/rigid/velocity_solver.zig @@ -25,6 +25,18 @@ const BodyManager = bm_mod.BodyManager; const ContactConstraint = cc.ContactConstraint; const ContactCache = cache_mod.ContactCache; +/// Velocity-solver tuning. +pub const SolverConfig = struct { + /// Gauss-Seidel velocity iteration passes per tick (named for the M1.1.7 + /// `position_iterations` sibling). + velocity_iterations: u32 = 8, + /// Restitution cutoff (m/s): a bounce is applied only when the pre-solve + /// relative normal speed exceeds this — a PHYSICAL velocity constant (config + /// field), not a geometric epsilon. Below it, low-speed contacts settle + /// without jitter. + restitution_threshold: Real = 1.0, +}; + /// Apply a world-space impulse `p` at the contact levers `r_a`/`r_b` to both /// bodies' velocities: body A receives −p, body B receives +p (the normal is /// A→B). Naturally a no-op on an infinite-mass body (inv_mass and inv_inertia @@ -88,6 +100,50 @@ pub fn storeContacts(gpa: std.mem.Allocator, cache: *ContactCache, constraints: } } +/// Solve the velocity constraints over the constraint index range `[from, to)` +/// with `cfg.velocity_iterations` Gauss-Seidel passes. Constraints are visited in +/// ascending pair-key order (as `build` sorted them) and points in manifold +/// order. M1.1.6 always passes the full range; the explicit range is the +/// island-additivity seam (M1.1.8 sorts constraints into contiguous per-island +/// ranges — purely additive). +pub fn solveRange(bm: *BodyManager, constraints: []ContactConstraint, from: usize, to: usize, cfg: SolverConfig) void { + var iter: u32 = 0; + while (iter < cfg.velocity_iterations) : (iter += 1) { + for (constraints[from..to]) |*c| { + for (0..c.count) |i| { + solveNormalPoint(bm, c, &c.points[i], cfg); + } + } + } +} + +/// One Gauss-Seidel normal-impulse update for a contact point: drive the relative +/// normal velocity toward its restitution target, with the accumulated-impulse +/// clamp `λₙ ≥ 0` (Catto — the solver can only push, never pull). +fn solveNormalPoint(bm: *BodyManager, c: *const ContactConstraint, pt: *cc.ConstraintPoint, cfg: SolverConfig) void { + const v_a = bm.linearVelocity(c.body_a).?; + const w_a = bm.angularVelocity(c.body_a).?; + const v_b = bm.linearVelocity(c.body_b).?; + const w_b = bm.angularVelocity(c.body_b).?; + const v_rel = v_b.add(w_b.cross(pt.r_b)).sub(v_a.add(w_a.cross(pt.r_a))); + const v_n = c.normal.dot(v_rel); + + // Restitution bias — the solver's ONLY velocity bias (no Baumgarte, no + // positional bias): target a separating speed −e·v_n⁻ when the pre-solve + // approach speed exceeds the threshold, otherwise target rest (0). + const restitution_bias: Real = if (@abs(pt.rel_normal_velocity) > cfg.restitution_threshold) + -c.restitution * pt.rel_normal_velocity + else + 0; + + // Impulse driving v_n toward the target (Δv_n = kₙ·Δλ, kₙ = 1/normal_mass). + const dlambda = pt.normal_mass * (restitution_bias - v_n); + const new_lambda = @max(@as(Real, 0), pt.normal_impulse + dlambda); + const applied = new_lambda - pt.normal_impulse; + pt.normal_impulse = new_lambda; + applyImpulse(bm, c, pt.r_a, pt.r_b, c.normal.scale(applied)); +} + // --- tests ------------------------------------------------------------------- const api = @import("weld_forge"); @@ -241,3 +297,142 @@ test "solved impulses round-trip through the cache to the next tick's warm start try testing.expectApproxEqAbs(@as(Real, 7), constraints.items[0].points[0].normal_impulse, 1e-5); try testing.expectEqual(@as(u32, 1), cache.hits); } + +/// Dynamic sphere A at origin + static sphere B at +0.9X (overlapping along +X), +/// both with restitution `e`. The caller sets A's approach velocity, builds, and +/// solves. +fn sphereHitScene(gpa: std.mem.Allocator, store: *ShapeStore, bm: *BodyManager, e: f32) ![2]api.BodyId { + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + var da = descOf(0, .dynamic, s); + da.mass = 1; + da.restitution = e; + var db = descOf(1, .static, s); + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + db.restitution = e; + const id_a = try bm.addBody(gpa, store, da); + const id_b = try bm.addBody(gpa, store, db); + return .{ id_a, id_b }; +} + +test "solver config defaults" { + const cfg = SolverConfig{}; + try testing.expectEqual(@as(u32, 8), cfg.velocity_iterations); + try testing.expectEqual(@as(Real, 1), cfg.restitution_threshold); +} + +test "normal solve kills the approach velocity at e = 0" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try sphereHitScene(gpa, &store, &bm, 0); + bm.setLinearVelocity(ids[0], vr(3, 0, 0)); // approaching +X + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + + // Approach killed; the static body never moves; λₙ respected the Catto clamp. + try testing.expectApproxEqAbs(@as(Real, 0), bm.linearVelocity(ids[0]).?.toArray()[0], 1e-5); + try testing.expect(bm.linearVelocity(ids[1]).?.approxEql(Vec3r.zero, 0)); + try testing.expect(constraints.items[0].points[0].normal_impulse >= 0); +} + +test "restitution bounces above the threshold and is inert below it" { + const gpa = testing.allocator; + + // Above threshold: 3 m/s approach, e = 0.8 ⇒ rebound at 0.8·3 = 2.4 m/s (−X). + { + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try sphereHitScene(gpa, &store, &bm, 0.8); + bm.setLinearVelocity(ids[0], vr(3, 0, 0)); + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + try testing.expectApproxEqAbs(@as(Real, -2.4), bm.linearVelocity(ids[0]).?.toArray()[0], 1e-4); + } + + // Below threshold: 0.5 m/s approach < 1.0 ⇒ no bounce, approach killed. + { + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try sphereHitScene(gpa, &store, &bm, 0.8); + bm.setLinearVelocity(ids[0], vr(0.5, 0, 0)); + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + try testing.expectApproxEqAbs(@as(Real, 0), bm.linearVelocity(ids[0]).?.toArray()[0], 1e-5); + } +} + +test "accumulated normal impulse stays non-negative (Catto clamp)" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try sphereHitScene(gpa, &store, &bm, 0); + bm.setLinearVelocity(ids[0], vr(-3, 0, 0)); // moving AWAY (separating) + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + + // A separating contact needs no push: the clamp keeps λₙ at 0 (never pulls) + // and A keeps moving away unchanged. + try testing.expectEqual(@as(Real, 0), constraints.items[0].points[0].normal_impulse); + try testing.expectApproxEqAbs(@as(Real, -3), bm.linearVelocity(ids[0]).?.toArray()[0], 1e-5); +} + +test "solveRange solves only the given constraint index range" { + 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 } }); + + // Pair (0,1) near the origin, pair (2,3) ten metres away — both a dynamic A + // approaching a static B along +X. pair_key(0,1) = 1 < pair_key(2,3), so the + // first pair is constraint 0. + var a0 = descOf(0, .dynamic, s); + a0.mass = 1; + a0.restitution = 0; + var b0 = descOf(1, .static, s); + b0.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + b0.restitution = 0; + var a1 = descOf(2, .dynamic, s); + a1.mass = 1; + a1.restitution = 0; + a1.position = foundation.math.Vec3.fromArray(.{ 10, 0, 0 }); + var b1 = descOf(3, .static, s); + b1.position = foundation.math.Vec3.fromArray(.{ 10.9, 0, 0 }); + b1.restitution = 0; + const id_a0 = try bm.addBody(gpa, &store, a0); + const id_b0 = try bm.addBody(gpa, &store, b0); + const id_a1 = try bm.addBody(gpa, &store, a1); + const id_b1 = try bm.addBody(gpa, &store, b1); + bm.setLinearVelocity(id_a0, vr(3, 0, 0)); + bm.setLinearVelocity(id_a1, vr(3, 0, 0)); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{ pairKey(id_a0, id_b0), pairKey(id_a1, id_b1) }); + try testing.expectEqual(@as(usize, 2), constraints.items.len); + + // Solve only the first constraint's range; the second pair is untouched. + solveRange(&bm, constraints.items, 0, 1, .{}); + + try testing.expectApproxEqAbs(@as(Real, 0), bm.linearVelocity(id_a0).?.toArray()[0], 1e-5); + try testing.expectApproxEqAbs(@as(Real, 3), bm.linearVelocity(id_a1).?.toArray()[0], 1e-5); +} From b504206a8de8b57ff97f92fa998bf9443dd63c2d Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:12:32 +0200 Subject: [PATCH 12/23] docs(forge): log E4 (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index ce8a0b6..85a0381 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -157,6 +157,12 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - `rigid/root.zig` re-exports `ContactCache`/`CacheKey`/`CacheValue`/`warmStart`/`storeContacts` + pins both new files. - Verified: 142/142 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - Note for review: `endTick` is `sort + swap` per the brief; the harvest (`storeContacts`) is a distinct step run after the iterations and before `endTick` (the current buffer must be filled before it is sorted) — the E6 normative cycle's implicit fill step. +- 2026-07-24 16:12 — E4 (velocity iterations + restitution), RED-first (`56c2878`): + - `SolverConfig` (`velocity_iterations: u32 = 8`, `restitution_threshold: Real = 1.0` m/s — a physical velocity constant, config field, not a geometric epsilon). + - `solveRange(bm, constraints, from, to, cfg)`: Gauss-Seidel, `cfg.velocity_iterations` passes, constraints in ascending pair-key order (as `build` sorted them), points in manifold order. M1.1.6 always passes the full range; the explicit range is the island-additivity seam (M1.1.8). + - `solveNormalPoint`: normal impulse `Δλ = normal_mass·(restitution_bias − v_n)`, accumulated clamp `λₙ ≥ 0` (Catto). Restitution bias `−e·v_n⁻` applied iff `|v_n⁻| > restitution_threshold` (`v_n⁻` captured in `prepare`, before `warmStart`) — the solver's ONLY velocity bias (no Baumgarte, no positional bias). Friction is E5. + - RED observed: the "restitution bounces above the threshold" test failed (actual 0, expected −2.4) under a no-restitution normal solve, then GREEN after adding the restitution bias. + - Verified: 147/147 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. `velocity_solver.zig` = 438 lines (under the 500-line guideline). ## Recorded deviations From 9b0add674aa95c589c8382feb5d699ee5408c431 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:27:00 +0200 Subject: [PATCH 13/23] fix(forge): restrict restitution bias to approaching contacts (M1.1.6) --- .../forge/forge_3d/rigid/velocity_solver.zig | 41 +++++++++++++++++-- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/src/modules/forge/forge_3d/rigid/velocity_solver.zig b/src/modules/forge/forge_3d/rigid/velocity_solver.zig index bd6fec2..f973447 100644 --- a/src/modules/forge/forge_3d/rigid/velocity_solver.zig +++ b/src/modules/forge/forge_3d/rigid/velocity_solver.zig @@ -129,9 +129,12 @@ fn solveNormalPoint(bm: *BodyManager, c: *const ContactConstraint, pt: *cc.Const const v_n = c.normal.dot(v_rel); // Restitution bias — the solver's ONLY velocity bias (no Baumgarte, no - // positional bias): target a separating speed −e·v_n⁻ when the pre-solve - // approach speed exceeds the threshold, otherwise target rest (0). - const restitution_bias: Real = if (@abs(pt.rel_normal_velocity) > cfg.restitution_threshold) + // positional bias): when the pre-solve contact is APPROACHING faster than the + // threshold (v_n⁻ < −threshold), target a separating rebound speed −e·v_n⁻ + // (> 0); otherwise target rest (0). Keyed on the approach direction, never on + // |v_n⁻|, so a capture-time SEPARATING contact never yields a negative target + // that would under-enforce non-penetration. + const restitution_bias: Real = if (pt.rel_normal_velocity < -cfg.restitution_threshold) -c.restitution * pt.rel_normal_velocity else 0; @@ -436,3 +439,35 @@ test "solveRange solves only the given constraint index range" { try testing.expectApproxEqAbs(@as(Real, 0), bm.linearVelocity(id_a0).?.toArray()[0], 1e-5); try testing.expectApproxEqAbs(@as(Real, 3), bm.linearVelocity(id_a1).?.toArray()[0], 1e-5); } + +test "capture-time separating contact still enforces non-penetration" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const ids = try sphereHitScene(gpa, &store, &bm, 0.8); + + // At capture (build/prepare) the contact is SEPARATING faster than the + // threshold (v_n⁻ = +3 > 1.0), so restitution must NOT arm here — a bias + // keyed on |v_n⁻| would (wrongly) target a NEGATIVE separating speed. + bm.setLinearVelocity(ids[0], vr(-3, 0, 0)); + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(ids[0], ids[1])}); + + // The body then actually approaches at solve time (relative v_n = −2). A + // capture-time separating velocity must not grant a penetration allowance: + // the solve must still drive v_n to rest (≈ 0), never to a negative target. + bm.setLinearVelocity(ids[0], vr(2, 0, 0)); + solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + + const c = constraints.items[0]; + const va = bm.linearVelocity(ids[0]).?; + const wa = bm.angularVelocity(ids[0]).?; + const vb = bm.linearVelocity(ids[1]).?; + const wb = bm.angularVelocity(ids[1]).?; + const v_rel = vb.add(wb.cross(c.points[0].r_b)).sub(va.add(wa.cross(c.points[0].r_a))); + const v_n = c.normal.dot(v_rel); + try testing.expectApproxEqAbs(@as(Real, 0), v_n, 1e-4); +} From 917a3a22a89656282492d0c13dea5e59ce6ec73d Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:27:34 +0200 Subject: [PATCH 14/23] docs(forge): record restitution condition deviation (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index 85a0381..9f38546 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -163,10 +163,11 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - `solveNormalPoint`: normal impulse `Δλ = normal_mass·(restitution_bias − v_n)`, accumulated clamp `λₙ ≥ 0` (Catto). Restitution bias `−e·v_n⁻` applied iff `|v_n⁻| > restitution_threshold` (`v_n⁻` captured in `prepare`, before `warmStart`) — the solver's ONLY velocity bias (no Baumgarte, no positional bias). Friction is E5. - RED observed: the "restitution bounces above the threshold" test failed (actual 0, expected −2.4) under a no-restitution normal solve, then GREEN after adding the restitution bias. - Verified: 147/147 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. `velocity_solver.zig` = 438 lines (under the 500-line guideline). +- 2026-07-24 16:27 — E4 review fix (`9b0add6`): restitution condition `|v_n⁻| > threshold` → `v_n⁻ < −threshold` (approaching-only) + aligned doc comment (see Recorded deviations for the frozen-text deviation). RED-first test `"capture-time separating contact still enforces non-penetration"`: capture the contact SEPARATING (`v_n⁻ = +3`), then approach at solve time (relative `v_n = −2`), assert post-solve `v_n ≈ 0`. Observed RED `actual −2` under the old `|v_n⁻|` condition (the negative bias fully suppresses the corrective impulse via the clamp — same defect class as the ≈ −2.4 estimate, which corresponds to an approach faster than the bias magnitude), GREEN after the fix. 148/148 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; `zig build` + `zig fmt --check` clean. ## Recorded deviations -- — +- `9b0add6` — **Restitution condition deviates from the FROZEN Scope.** The frozen text says the bias `−e·v_n⁻` is applied "iff `|v_n⁻| > restitution_threshold`"; implemented instead as `v_n⁻ < −restitution_threshold` (approaching-only). Reason: the frozen `|v_n⁻|` form yields a NEGATIVE bias target on a capture-time SEPARATING contact (`v_n⁻ > +threshold`), which under-enforces non-penetration when the body actually approaches at solve time (the negative target either suppresses the corrective impulse via the `λₙ ≥ 0` clamp or, for a faster approach, lets the bodies close at the bias speed). The corrected condition keys the bias on the approach direction, so the bias is never negative. Claude.ai spec defect, ratified in the E4 review (Guy, 2026-07-24); the frozen brief text is left as-is pending a Claude.ai round-trip / KB re-upload. ## Blockers encountered From 482e3d2dd417d7e6e6a8cf4feb4b2178b16a5fe1 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:47:23 +0200 Subject: [PATCH 15/23] feat(forge): add friction with circular clamp (M1.1.6) --- .../forge/forge_3d/rigid/velocity_solver.zig | 119 ++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/src/modules/forge/forge_3d/rigid/velocity_solver.zig b/src/modules/forge/forge_3d/rigid/velocity_solver.zig index f973447..f4613a4 100644 --- a/src/modules/forge/forge_3d/rigid/velocity_solver.zig +++ b/src/modules/forge/forge_3d/rigid/velocity_solver.zig @@ -111,7 +111,10 @@ pub fn solveRange(bm: *BodyManager, constraints: []ContactConstraint, from: usiz while (iter < cfg.velocity_iterations) : (iter += 1) { for (constraints[from..to]) |*c| { for (0..c.count) |i| { + // Jolt order: the normal impulse first, then friction clamped + // against the CURRENT accumulated normal impulse of this point. solveNormalPoint(bm, c, &c.points[i], cfg); + solveFrictionPoint(bm, c, &c.points[i]); } } } @@ -147,6 +150,47 @@ fn solveNormalPoint(bm: *BodyManager, c: *const ContactConstraint, pt: *cc.Const applyImpulse(bm, c, pt.r_a, pt.r_b, c.normal.scale(applied)); } +/// One Gauss-Seidel friction update for a contact point (run after the normal +/// update): two INDEPENDENT per-axis tangential impulses driving the tangential +/// velocity toward 0 (coupling ignored — Jolt convention), then a CIRCULAR clamp +/// of the accumulated `(λ_t1, λ_t2)` pair to the friction cone `‖λ_t‖ ≤ μ·λₙ` +/// (against the current accumulated normal impulse). The circular clamp is +/// isotropic and basis-independent — coherent with the world-space tangent cache; +/// the box clamp is anisotropic (up to √2·μ·λₙ on the diagonal) and basis-biased. +fn solveFrictionPoint(bm: *BodyManager, c: *const ContactConstraint, pt: *cc.ConstraintPoint) void { + const v_a = bm.linearVelocity(c.body_a).?; + const w_a = bm.angularVelocity(c.body_a).?; + const v_b = bm.linearVelocity(c.body_b).?; + const w_b = bm.angularVelocity(c.body_b).?; + const v_rel = v_b.add(w_b.cross(pt.r_b)).sub(v_a.add(w_a.cross(pt.r_a))); + + // Independent per-axis tangential impulses driving each axis toward 0. + const v_t1 = c.tangent1.dot(v_rel); + const v_t2 = c.tangent2.dot(v_rel); + var new_t1 = pt.tangent1_impulse - pt.tangent1_mass * v_t1; + var new_t2 = pt.tangent2_impulse - pt.tangent2_mass * v_t2; + + // Circular (isotropic, basis-independent) clamp of the accumulated tangent + // pair to the friction cone μ·λₙ — rescale the vector, never per-axis (the box + // form is anisotropic, up to √2·μ·λₙ on the diagonal, and basis-biased). The + // branch guarantees `len_sq > max_friction² ≥ 0`, so the √ divisor is nonzero. + const max_friction = c.friction * pt.normal_impulse; + const len_sq = new_t1 * new_t1 + new_t2 * new_t2; + if (len_sq > max_friction * max_friction) { + const scale = max_friction / @sqrt(len_sq); + new_t1 *= scale; + new_t2 *= scale; + } + + // Apply the accumulated delta on both axes. + const applied_t1 = new_t1 - pt.tangent1_impulse; + const applied_t2 = new_t2 - pt.tangent2_impulse; + pt.tangent1_impulse = new_t1; + pt.tangent2_impulse = new_t2; + const impulse = c.tangent1.scale(applied_t1).add(c.tangent2.scale(applied_t2)); + applyImpulse(bm, c, pt.r_a, pt.r_b, impulse); +} + // --- tests ------------------------------------------------------------------- const api = @import("weld_forge"); @@ -471,3 +515,78 @@ test "capture-time separating contact still enforces non-penetration" { const v_n = c.normal.dot(v_rel); try testing.expectApproxEqAbs(@as(Real, 0), v_n, 1e-4); } + +test "friction cancels tangential sliding below the cone" { + 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 } }); + var da = descOf(0, .dynamic, s); + da.mass = 1; + da.friction = 1; + var db = descOf(1, .dynamic, s); + db.mass = 1; + db.friction = 1; + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + const id_a = try bm.addBody(gpa, &store, da); + const id_b = try bm.addBody(gpa, &store, db); + bm.setLinearVelocity(id_a, vr(0, 2, 0)); // sliding tangentially (+Y) + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(id_a, id_b)}); + const c = &constraints.items[0]; + c.points[0].normal_impulse = 10; // large λₙ ⇒ wide cone ⇒ no clamp + solveFrictionPoint(&bm, c, &c.points[0]); + + // The relative contact-point tangential velocity is driven to ≈ 0. + const va = bm.linearVelocity(id_a).?; + const wa = bm.angularVelocity(id_a).?; + const vb = bm.linearVelocity(id_b).?; + const wb = bm.angularVelocity(id_b).?; + const v_rel = vb.add(wb.cross(c.points[0].r_b)).sub(va.add(wa.cross(c.points[0].r_a))); + try testing.expectApproxEqAbs(@as(Real, 0), c.tangent1.dot(v_rel), 1e-5); + try testing.expectApproxEqAbs(@as(Real, 0), c.tangent2.dot(v_rel), 1e-5); +} + +/// Clamped friction impulse magnitude for a body sliding at `slide` against a +/// fixed normal impulse λₙ = 1 (cone μ·λₙ = 1). Used to compare slide directions. +fn frictionMagnitude(gpa: std.mem.Allocator, slide: Vec3r) !Real { + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + var da = descOf(0, .dynamic, s); + da.mass = 1; + da.friction = 1; + var db = descOf(1, .dynamic, s); + db.mass = 1; + db.friction = 1; + db.position = foundation.math.Vec3.fromArray(.{ 0.9, 0, 0 }); + const id_a = try bm.addBody(gpa, &store, da); + const id_b = try bm.addBody(gpa, &store, db); + bm.setLinearVelocity(id_a, slide); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try cc.build(gpa, &constraints, &bm, &store, &.{pairKey(id_a, id_b)}); + const c = &constraints.items[0]; + c.points[0].normal_impulse = 1; + solveFrictionPoint(&bm, c, &c.points[0]); + const pt = c.points[0]; + return @sqrt(pt.tangent1_impulse * pt.tangent1_impulse + pt.tangent2_impulse * pt.tangent2_impulse); +} + +test "friction clamp is isotropic (circular, basis-independent)" { + const gpa = testing.allocator; + // Slide fast (10 m/s) so the cone μ·λₙ = 1 clamps. Along a single tangent axis + // vs the tangent-plane diagonal must yield the SAME clamped magnitude (circular + // clamp); a box clamp would give √2 on the diagonal. + const along = try frictionMagnitude(gpa, vr(0, 10, 0)); + const diagonal = try frictionMagnitude(gpa, vr(0, 7.0710678, 7.0710678)); + try testing.expectApproxEqAbs(@as(Real, 1), along, 1e-4); // clamped to μ·λₙ = 1 + try testing.expectApproxEqAbs(along, diagonal, 1e-4); // isotropic +} From 305bd150c073dd4ef2bf6ccca0ad4da409049874 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 16:47:56 +0200 Subject: [PATCH 16/23] docs(forge): log E5 (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index 9f38546..6ecd597 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -164,6 +164,11 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - RED observed: the "restitution bounces above the threshold" test failed (actual 0, expected −2.4) under a no-restitution normal solve, then GREEN after adding the restitution bias. - Verified: 147/147 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. `velocity_solver.zig` = 438 lines (under the 500-line guideline). - 2026-07-24 16:27 — E4 review fix (`9b0add6`): restitution condition `|v_n⁻| > threshold` → `v_n⁻ < −threshold` (approaching-only) + aligned doc comment (see Recorded deviations for the frozen-text deviation). RED-first test `"capture-time separating contact still enforces non-penetration"`: capture the contact SEPARATING (`v_n⁻ = +3`), then approach at solve time (relative `v_n = −2`), assert post-solve `v_n ≈ 0`. Observed RED `actual −2` under the old `|v_n⁻|` condition (the negative bias fully suppresses the corrective impulse via the clamp — same defect class as the ≈ −2.4 estimate, which corresponds to an approach faster than the bias magnitude), GREEN after the fix. 148/148 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; `zig build` + `zig fmt --check` clean. +- 2026-07-24 16:47 — E5 (friction), RED-first (`482e3d2`): + - `solveFrictionPoint` (run after `solveNormalPoint` per point — Jolt order, so it clamps against the current accumulated λₙ): two INDEPENDENT per-axis tangential impulses (using the precomputed tangent effective masses, coupling ignored) driving the contact-point tangential velocity toward 0, then a CIRCULAR clamp of the accumulated `(λ_t1, λ_t2)` pair to `‖λ_t‖ ≤ μ·λₙ` by rescaling — isotropic and basis-independent (the box clamp is anisotropic, up to √2·μ·λₙ on the diagonal, and basis-biased). `solveRange` now calls normal then friction per point. + - RED observed: the isotropy test (`"friction clamp is isotropic"`) failed under a per-axis BOX clamp (diagonal slide clamped to √2·μ·λₙ vs μ·λₙ along an axis), GREEN after the circular rescale. Second friction test: friction cancels sub-cone tangential sliding (contact-point tangential velocity → 0). + - Verified: 150/150 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. + - Note for review: `velocity_solver.zig` = 592 lines (over the 500-line Review guideline) — the brief mandates warm start + iterations + friction + `SolverConfig` + range entry in this single file, plus inline tests (§13); conscious brief-driven overage, matching the `narrowphase.zig`/`broadphase.zig` precedent. Splitting the inline tests into a sibling `_test.zig` is available if preferred. ## Recorded deviations From 9d13e9123900dc9cf19e52d0eb80a55a6cb2fe2a Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 17:57:40 +0200 Subject: [PATCH 17/23] feat(forge): add solver harness and acceptance suite (M1.1.6) --- src/modules/forge/forge_3d/rigid/root.zig | 5 + src/modules/forge/forge_3d/root.zig | 1 + .../forge/forge_3d/tests/solver_test.zig | 490 ++++++++++++++++++ 3 files changed, 496 insertions(+) create mode 100644 src/modules/forge/forge_3d/tests/solver_test.zig diff --git a/src/modules/forge/forge_3d/rigid/root.zig b/src/modules/forge/forge_3d/rigid/root.zig index febef9c..cfc3839 100644 --- a/src/modules/forge/forge_3d/rigid/root.zig +++ b/src/modules/forge/forge_3d/rigid/root.zig @@ -36,10 +36,15 @@ pub const CacheKey = contact_cache.CacheKey; /// A warm-start cache value (accumulated normal + world tangent impulse). pub const CacheValue = contact_cache.CacheValue; +/// Velocity-solver tuning (iteration count, restitution threshold). +pub const SolverConfig = velocity_solver.SolverConfig; + /// Warm-start constraints from the previous tick's cache (tangent reprojection). pub const warmStart = velocity_solver.warmStart; /// Harvest solved constraint impulses into the cache's current buffer. pub const storeContacts = velocity_solver.storeContacts; +/// Solve the velocity constraints over an index range (M1.1.6 passes the full range). +pub const solveRange = velocity_solver.solveRange; comptime { _ = contact_constraint; diff --git a/src/modules/forge/forge_3d/root.zig b/src/modules/forge/forge_3d/root.zig index 725e77b..101d6fd 100644 --- a/src/modules/forge/forge_3d/root.zig +++ b/src/modules/forge/forge_3d/root.zig @@ -171,4 +171,5 @@ comptime { _ = @import("tests/manifold_test.zig"); _ = @import("tests/fast_paths_test.zig"); _ = @import("tests/epa_robustness_test.zig"); + _ = @import("tests/solver_test.zig"); } diff --git a/src/modules/forge/forge_3d/tests/solver_test.zig b/src/modules/forge/forge_3d/tests/solver_test.zig new file mode 100644 index 0000000..704d96b --- /dev/null +++ b/src/modules/forge/forge_3d/tests/solver_test.zig @@ -0,0 +1,490 @@ +//! M1.1.6 acceptance suite for the Sequential Impulses contact solver. +//! +//! `World` composes the full per-tick pipeline IN TESTS ONLY (the production +//! `step()` orchestration is M1.1.15). Normative per-tick cycle (brief E6): +//! (1) `Broadphase.computePairs` on the current poses (moved-driven deltas) +//! (2) narrowphase `collidePair` per candidate — done inside `rigid.build` +//! (3) `integrateVelocities(dt, gravity)` +//! (4) `cache.beginTick` → `build` + `prepare` (captures v_n⁻) → `warmStart` +//! → `solveRange` × `velocity_iterations` +//! (5) `integratePositions(dt)` +//! (6) `storeContacts` (harvest) → `cache.endTick` (sort + swap) +//! (7) broadphase proxy updates. +//! +//! `computePairs` is moved-driven (it reports only pairs touching a proxy that +//! moved since the last call), so the harness keeps a PERSISTENT candidate set +//! (`active`, a sorted-deduped key list): it merges each tick's deltas and, after +//! the solve, retains only the pairs that produced a constraint (contacting +//! ones), so a resting contact — whose proxies stop moving and thus stop being +//! reported — stays live and keeps being solved (no sinking). A real `step()` / +//! `PhysicsWorld` at M1.1.15 owns the same persistent pair set. + +const std = @import("std"); +const config = @import("../config.zig"); +const shape_mod = @import("../shape.zig"); +const bm_mod = @import("../body_manager.zig"); +const broadphase = @import("../pipeline/broadphase.zig"); +const integration = @import("../pipeline/integration.zig"); +const rigid = @import("../rigid/root.zig"); +const api = @import("weld_forge"); +const foundation = @import("foundation"); + +const Real = config.Real; +const Vec3r = config.Vec3r; +const ShapeStore = shape_mod.ShapeStore; +const BodyManager = bm_mod.BodyManager; +const BodyId = api.BodyId; +const Bp = broadphase.Broadphase(Real); +const ContactConstraint = rigid.ContactConstraint; +const ContactCache = rigid.ContactCache; +const SolverConfig = rigid.SolverConfig; +const testing = std.testing; + +fn vr(x: Real, y: Real, z: Real) Vec3r { + return Vec3r.fromArray(.{ x, y, z }); +} + +fn av3(x: f32, y: f32, z: f32) foundation.math.Vec3 { + return foundation.math.Vec3.fromArray(.{ x, y, z }); +} + +fn broadphaseLayer(bt: api.BodyType) broadphase.BroadphaseLayer { + return if (bt == .static) .static else .dynamic; +} + +fn sortDedup(list: *std.ArrayListUnmanaged(u64)) void { + std.mem.sort(u64, list.items, {}, std.sort.asc(u64)); + if (list.items.len == 0) return; + var w: usize = 1; + var i: usize = 1; + while (i < list.items.len) : (i += 1) { + if (list.items[i] != list.items[w - 1]) { + list.items[w] = list.items[i]; + w += 1; + } + } + list.shrinkRetainingCapacity(w); +} + +const BodyProxy = struct { id: BodyId, proxy: Bp.Proxy }; + +/// A minimal physics world composing the full contact-solver pipeline for tests. +const World = struct { + store: ShapeStore = .{}, + bm: BodyManager = .{}, + bp: Bp, + cache: ContactCache = .{}, + cfg: SolverConfig = .{}, + gravity: Vec3r, + dt: Real, + bodies: std.ArrayListUnmanaged(BodyProxy) = .empty, + active: std.ArrayListUnmanaged(u64) = .empty, + constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty, + scratch: std.ArrayListUnmanaged(Bp.Pair) = .empty, + + fn init(gravity: Vec3r, dt: Real) World { + return .{ .bp = Bp.init(.{}), .gravity = gravity, .dt = dt }; + } + + fn deinit(self: *World, gpa: std.mem.Allocator) void { + self.store.deinit(gpa); + self.bm.deinit(gpa); + self.bp.deinit(gpa); + self.cache.deinit(gpa); + self.bodies.deinit(gpa); + self.active.deinit(gpa); + self.constraints.deinit(gpa); + self.scratch.deinit(gpa); + self.* = undefined; + } + + fn addBody(self: *World, gpa: std.mem.Allocator, desc: api.BodyDescriptor) !BodyId { + const id = try self.bm.addBody(gpa, &self.store, desc); + const aabb = self.bm.bodyAabb(&self.store, id).?; + const proxy = try self.bp.insert(gpa, broadphaseLayer(desc.body_type), aabb, id); + try self.bodies.append(gpa, .{ .id = id, .proxy = proxy }); + return id; + } + + fn step(self: *World, gpa: std.mem.Allocator) !void { + // (1) broadphase candidate deltas → persistent active set. + try self.bp.computePairs(gpa, &self.scratch); + for (self.scratch.items) |p| try self.active.append(gpa, (@as(u64, p.a) << 32) | p.b); + sortDedup(&self.active); + + // (3) integrate velocities (gravity + accumulators + clamped damping). + integration.integrateVelocities(&self.bm, self.dt, self.gravity); + + // (4) solve: prepare captures v_n⁻ (post-gravity), then warm start + iterate. + self.cache.beginTick(); + try rigid.build(gpa, &self.constraints, &self.bm, &self.store, self.active.items); + rigid.warmStart(&self.bm, &self.cache, self.constraints.items); + rigid.solveRange(&self.bm, self.constraints.items, 0, self.constraints.items.len, self.cfg); + + // (5) integrate positions from the solved velocities. + integration.integratePositions(&self.bm, self.dt); + + // (6) harvest solved impulses into the cache, then finalize (sort + swap). + try rigid.storeContacts(gpa, &self.cache, self.constraints.items); + self.cache.endTick(); + + // Persist only contacting pairs; separated ones drop (re-added by + // computePairs when their proxies move back into overlap). + self.active.clearRetainingCapacity(); + for (self.constraints.items) |c| try self.active.append(gpa, c.pair_key); + + // (7) broadphase proxy updates to the new poses. + for (self.bodies.items) |b| { + if (self.bm.bodyAabb(&self.store, b.id)) |aabb| try self.bp.update(gpa, b.proxy, aabb); + } + } +}; + +/// Ground (static box, top at y = 0.5) + a dynamic unit box dropped from y = `drop_y`. +/// Returns the dynamic box's BodyId. Restitution `e` on both. +fn groundAndBox(gpa: std.mem.Allocator, world: *World, drop_y: Real, e: f32) !BodyId { + const ground_shape = try world.store.createShape(gpa, .{ .box = .{ .half_extents = av3(5, 0.5, 5) } }); + const box_shape = try world.store.createShape(gpa, .{ .box = .{ .half_extents = av3(0.5, 0.5, 0.5) } }); + + var ground = api.BodyDescriptor{ .entity = .{ .index = 0, .generation = 0 }, .body_type = .static, .shape = ground_shape }; + ground.restitution = e; + _ = try world.addBody(gpa, ground); // centre at origin ⇒ top face at y = 0.5 + + var box = api.BodyDescriptor{ .entity = .{ .index = 1, .generation = 0 }, .body_type = .dynamic, .shape = box_shape }; + box.mass = 1; + box.restitution = e; + box.position = av3(0, @floatCast(drop_y), 0); + return world.addBody(gpa, box); +} + +test "box dropped on a static ground comes to rest without sinking (e = 0)" { + const gpa = testing.allocator; + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const box = try groundAndBox(gpa, &world, 2.0, 0); + + // Settle (2 s). + var t: u32 = 0; + while (t < 120) : (t += 1) try world.step(gpa); + const y_settled = world.bm.position(box).?.toArray()[1]; + + // Run 3 more seconds — the box must not sink further. + while (t < 300) : (t += 1) try world.step(gpa); + const y_final = world.bm.position(box).?.toArray()[1]; + const v_final = world.bm.linearVelocity(box).?.toArray()[1]; + + // Analytical rest y = ground_top(0.5) + box_half(0.5) = 1.0; the velocity + // solver leaves a small constant penetration (no NGS recovery until M1.1.7). + try testing.expect(y_final > 0.5); // never fell through + try testing.expect(y_final < 1.02); // resting at/near the contact + try testing.expect(y_final > 0.8); // penetration bounded (< 0.2) + try testing.expect(@abs(y_final - y_settled) < 1e-3); // no sinking after settling + try testing.expect(@abs(v_final) < 0.05); // at rest +} + +test "solve is deterministic across identical runs" { + const gpa = testing.allocator; + + var w1 = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer w1.deinit(gpa); + const b1 = try groundAndBox(gpa, &w1, 2.0, 0.5); + var w2 = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer w2.deinit(gpa); + const b2 = try groundAndBox(gpa, &w2, 2.0, 0.5); + + var t: u32 = 0; + while (t < 300) : (t += 1) { + try w1.step(gpa); + try w2.step(gpa); + } + + const p1 = w1.bm.position(b1).?.toArray(); + const p2 = w2.bm.position(b2).?.toArray(); + const v1 = w1.bm.linearVelocity(b1).?.toArray(); + const v2 = w2.bm.linearVelocity(b2).?.toArray(); + const q1 = w1.bm.rotation(b1).?.toArray(); + const q2 = w2.bm.rotation(b2).?.toArray(); + inline for (0..3) |k| { + try testing.expectEqual(p1[k], p2[k]); + try testing.expectEqual(v1[k], v2[k]); + } + inline for (0..4) |k| try testing.expectEqual(q1[k], q2[k]); +} + +/// Max upward (+Y) velocity the dropped box reaches over `ticks` — the rebound +/// speed after impact (≈ 0 when it does not bounce). +fn maxReboundVy(gpa: std.mem.Allocator, drop_y: Real, e: f32, ticks: u32) !Real { + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const box = try groundAndBox(gpa, &world, drop_y, e); + var max_vy: Real = -1e30; + var t: u32 = 0; + while (t < ticks) : (t += 1) { + try world.step(gpa); + const vy = world.bm.linearVelocity(box).?.toArray()[1]; + if (vy > max_vy) max_vy = vy; + } + return max_vy; +} + +test "restitution bounces above the threshold and not below it" { + const gpa = testing.allocator; + // Drop from 2 m ⇒ impact ≈ 4.4 m/s > 1.0 threshold, e = 0.8 ⇒ strong rebound. + const high = try maxReboundVy(gpa, 2.0, 0.8, 120); + try testing.expect(high > 2.0); + // A box resting flush only ever approaches at the per-tick gravity speed + // (≈ 0.16 m/s ≪ the 1.0 threshold), so it settles without a bounce. (The + // genuine sub-threshold IMPACT case is pinned exactly by the E4 unit test.) + const low = try maxReboundVy(gpa, 1.0, 0.8, 120); + try testing.expect(low < 0.5); +} + +test "restitution is captured before warm start" { + const gpa = testing.allocator; + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const box = try groundAndBox(gpa, &world, 1.0, 0.8); // starts flush, e = 0.8 + + // Rest for 1 s: the per-tick gravity approach (≈ 0.16 m/s) is below the + // restitution threshold, so the box settles and the cache warms up. + var t: u32 = 0; + while (t < 60) : (t += 1) try world.step(gpa); + try testing.expect(world.cache.hits > 0); // warm start is matching the resting contact + try testing.expect(@abs(world.bm.linearVelocity(box).?.toArray()[1]) < 0.05); // resting + + // A real downward impact after resting: v_n⁻ is captured (in `prepare`) BEFORE + // `warmStart`, so the warm-started resting impulse does NOT suppress the bounce. + world.bm.setLinearVelocity(box, vr(0, -3, 0)); + var max_vy: Real = -1e30; + t = 0; + while (t < 60) : (t += 1) { + try world.step(gpa); + const vy = world.bm.linearVelocity(box).?.toArray()[1]; + if (vy > max_vy) max_vy = vy; + } + try testing.expect(max_vy > 1.5); // bounced (≈ e·3 = 2.4) +} + +/// Along-surface displacement of a box placed flush on a static incline (rotated +/// θ about Z) with combined friction from `mu`, after `ticks`. +fn inclineDrift(gpa: std.mem.Allocator, theta: Real, mu: f32, ticks: u32) !Real { + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const ground_shape = try world.store.createShape(gpa, .{ .box = .{ .half_extents = av3(5, 0.5, 5) } }); + const box_shape = try world.store.createShape(gpa, .{ .box = .{ .half_extents = av3(0.5, 0.5, 0.5) } }); + const rot = foundation.math.Quatf.fromAxisAngle(av3(0, 0, 1), @floatCast(theta)); + + var ground = api.BodyDescriptor{ .entity = .{ .index = 0, .generation = 0 }, .body_type = .static, .shape = ground_shape }; + ground.rotation = rot; + ground.friction = mu; + _ = try world.addBody(gpa, ground); + + // Box flush on the incline: centre at R·(0,1,0) (ground half + box half), same rotation. + var box = api.BodyDescriptor{ .entity = .{ .index = 1, .generation = 0 }, .body_type = .dynamic, .shape = box_shape }; + box.mass = 1; + box.friction = mu; + box.rotation = rot; + box.position = rot.rotateVec3(av3(0, 1, 0)); + const box_id = try world.addBody(gpa, box); + + const p0 = world.bm.position(box_id).?; + var t: u32 = 0; + while (t < ticks) : (t += 1) try world.step(gpa); + return world.bm.position(box_id).?.sub(p0).length(); +} + +test "inclined static box holds below the friction angle and slides above it" { + const gpa = testing.allocator; + const theta = std.math.pi / 6.0; // 30° ⇒ tan θ ≈ 0.577 + const holds = try inclineDrift(gpa, theta, 0.8, 180); // μ = 0.8 > tan θ ⇒ holds + const slides = try inclineDrift(gpa, theta, 0.3, 180); // μ = 0.3 < tan θ ⇒ slides + try testing.expect(holds < 0.1); + try testing.expect(slides > 0.5); +} + +/// Horizontal speed of a box sliding on flat ground after `ticks`, given the +/// initial horizontal velocity `v0` (default friction 0.5 on both). +fn slideSpeedAfter(gpa: std.mem.Allocator, v0: Vec3r, ticks: u32) !Real { + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const box = try groundAndBox(gpa, &world, 1.0, 0); // starts flush on the ground + world.bm.setLinearVelocity(box, v0); + var t: u32 = 0; + while (t < ticks) : (t += 1) try world.step(gpa); + const v = world.bm.linearVelocity(box).?.toArray(); + return @sqrt(v[0] * v[0] + v[2] * v[2]); +} + +test "friction deceleration is isotropic on flat ground" { + const gpa = testing.allocator; + // Same initial speed (3 m/s) along +X vs the XZ diagonal must decay equally — + // the circular clamp is basis-independent. + const along = try slideSpeedAfter(gpa, vr(3, 0, 0), 20); + const diagonal = try slideSpeedAfter(gpa, vr(2.1213203, 0, 2.1213203), 20); + try testing.expectApproxEqAbs(along, diagonal, 1e-2); +} + +// --- constraint/cache-level scenarios (direct build, no broadphase) ----------- + +fn descOf(idx: u32, bt: api.BodyType, shape: api.ShapeId) api.BodyDescriptor { + return .{ .entity = .{ .index = idx, .generation = 0 }, .body_type = bt, .shape = shape }; +} + +fn pairKey(a: BodyId, b: BodyId) u64 { + return (@as(u64, @min(a, b)) << 32) | @max(a, b); +} + +/// Warm-start a single sphere contact whose normal is `normalize(n_dir)` with a +/// cached WORLD tangent `t_world`, and return the reconstructed applied tangent +/// (λ_t1·t1 + λ_t2·t2) — the world tangent seeded into the new basis. +fn reconstructWarmTangent(gpa: std.mem.Allocator, n_dir: foundation.math.Vec3, t_world: Vec3r) !Vec3r { + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + var da = descOf(0, .dynamic, s); + da.friction = 1; + var db = descOf(1, .dynamic, s); + db.friction = 1; + db.position = n_dir.normalize().scale(0.9); // overlap along n_dir ⇒ normal ≈ n̂ + const a = try bm.addBody(gpa, &store, da); + const b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try rigid.build(gpa, &constraints, &bm, &store, &.{pairKey(a, b)}); + const c = &constraints.items[0]; + + var cache = ContactCache{}; + defer cache.deinit(gpa); + try cache.store(gpa, .{ .pair_key = c.pair_key, .feature_id = c.points[0].feature_id }, .{ .lambda_n = 10, .tangent_impulse = t_world }); + cache.endTick(); + cache.beginTick(); + rigid.warmStart(&bm, &cache, constraints.items); + return c.tangent1.scale(c.points[0].tangent1_impulse).add(c.tangent2.scale(c.points[0].tangent2_impulse)); +} + +test "warm-started tangent is continuous across a tangent-basis flip" { + const gpa = testing.allocator; + // Two normals straddling the x = y dominant-axis boundary (so the tangent + // basis flips discontinuously). The cached WORLD tangent must reconstruct to + // the SAME world direction for both — a per-basis-scalar cache would rotate it. + const t_world = vr(0.3, -0.3, 0.5); + const r1 = try reconstructWarmTangent(gpa, av3(1.0, 1.02, 0.0), t_world); // y-dominant + const r2 = try reconstructWarmTangent(gpa, av3(1.02, 1.0, 0.0), t_world); // x-dominant + try testing.expect(r1.approxEql(r2, 1e-2)); +} + +test "kinematic-static contact produces no constraint and no NaN" { + 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 } }); + const a = try bm.addBody(gpa, &store, descOf(0, .static, s)); + var db = descOf(1, .kinematic, s); + db.position = av3(0.9, 0, 0); + const b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try rigid.build(gpa, &constraints, &bm, &store, &.{pairKey(a, b)}); + rigid.solveRange(&bm, constraints.items, 0, constraints.items.len, .{}); + + try testing.expectEqual(@as(usize, 0), constraints.items.len); + for (bm.position(a).?.toArray()) |v| try testing.expect(!std.math.isNan(v)); + for (bm.position(b).?.toArray()) |v| try testing.expect(!std.math.isNan(v)); + for (bm.linearVelocity(b).?.toArray()) |v| try testing.expect(!std.math.isNan(v)); +} + +test "warm start hits a resting contact and misses on generation reuse" { + const gpa = testing.allocator; + + // Part A: a resting box in the full pipeline warm-starts (cache hits > 0). + { + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + _ = try groundAndBox(gpa, &world, 1.0, 0); + var t: u32 = 0; + while (t < 60) : (t += 1) try world.step(gpa); + try testing.expect(world.cache.hits > 0); + } + + // Part B: reusing a freed slot bumps the generation ⇒ a new BodyId ⇒ a new + // pair_key ⇒ the previous tick's cache entry no longer matches (cold start). + { + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .sphere = .{ .radius = 0.5 } }); + const a = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + var db = descOf(1, .dynamic, s); + db.position = av3(0.9, 0, 0); + const b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + var cache = ContactCache{}; + defer cache.deinit(gpa); + + cache.beginTick(); + try rigid.build(gpa, &constraints, &bm, &store, &.{pairKey(a, b)}); + constraints.items[0].points[0].normal_impulse = 5; + try rigid.storeContacts(gpa, &cache, constraints.items); + cache.endTick(); + + bm.removeBody(b); + const b2 = try bm.addBody(gpa, &store, db); // reuses b's slot, generation bumped + try testing.expect(b2 != b); + + cache.beginTick(); + try rigid.build(gpa, &constraints, &bm, &store, &.{pairKey(a, b2)}); + rigid.warmStart(&bm, &cache, constraints.items); + try testing.expectEqual(@as(u32, 0), cache.hits); + try testing.expect(cache.misses > 0); + } +} + +test "per-feature warm start: surviving points keep impulses, a vanished one cold-starts" { + const gpa = testing.allocator; + var store = ShapeStore{}; + defer store.deinit(gpa); + var bm = BodyManager{}; + defer bm.deinit(gpa); + const s = try store.createShape(gpa, .{ .box = .{ .half_extents = av3(0.5, 0.5, 0.5) } }); + const a = try bm.addBody(gpa, &store, descOf(0, .dynamic, s)); + var db = descOf(1, .dynamic, s); + db.position = av3(0, 0.9, 0); // stacked ⇒ a multi-point face-face manifold + const b = try bm.addBody(gpa, &store, db); + + var constraints: std.ArrayListUnmanaged(ContactConstraint) = .empty; + defer constraints.deinit(gpa); + try rigid.build(gpa, &constraints, &bm, &store, &.{pairKey(a, b)}); + const c = &constraints.items[0]; + const n: usize = c.count; + try testing.expect(n >= 2); + + // Seed the cache for every point EXCEPT the last (simulate that feature vanishing). + var cache = ContactCache{}; + defer cache.deinit(gpa); + var i: usize = 0; + while (i + 1 < n) : (i += 1) { + try cache.store(gpa, .{ .pair_key = c.pair_key, .feature_id = c.points[i].feature_id }, .{ .lambda_n = 3, .tangent_impulse = Vec3r.zero }); + } + cache.endTick(); + + cache.beginTick(); + rigid.warmStart(&bm, &cache, constraints.items); + + // Surviving points seeded from the cache; the last point cold-starts at 0. + i = 0; + while (i + 1 < n) : (i += 1) { + try testing.expectApproxEqAbs(@as(Real, 3), c.points[i].normal_impulse, 1e-5); + } + try testing.expectEqual(@as(Real, 0), c.points[n - 1].normal_impulse); + try testing.expectEqual(@as(u32, @intCast(n - 1)), cache.hits); + try testing.expectEqual(@as(u32, 1), cache.misses); +} From 3d7521fe0454005376bedc644734fc4da8bf3216 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 17:58:22 +0200 Subject: [PATCH 18/23] docs(forge): log E6 harness and acceptance suite (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index 6ecd597..7f895e2 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -169,6 +169,11 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - RED observed: the isotropy test (`"friction clamp is isotropic"`) failed under a per-axis BOX clamp (diagonal slide clamped to √2·μ·λₙ vs μ·λₙ along an axis), GREEN after the circular rescale. Second friction test: friction cancels sub-cone tangential sliding (contact-point tangential velocity → 0). - Verified: 150/150 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - Note for review: `velocity_solver.zig` = 592 lines (over the 500-line Review guideline) — the brief mandates warm start + iterations + friction + `SolverConfig` + range entry in this single file, plus inline tests (§13); conscious brief-driven overage, matching the `narrowphase.zig`/`broadphase.zig` precedent. Splitting the inline tests into a sibling `_test.zig` is available if preferred. +- 2026-07-24 17:57 — E6 (harness + acceptance suite) (`9d13e91`): + - `tests/solver_test.zig` (490 lines): `World` composes the normative per-tick cycle IN TESTS (computePairs → integrateVelocities → beginTick+build+prepare+warmStart+solveRange → integratePositions → storeContacts+endTick → broadphase proxy updates). Since `computePairs` is moved-driven, the harness keeps a PERSISTENT candidate set (sorted-deduped keys, merged each tick, pruned to contacting pairs post-solve) so resting contacts stay solved — a real `step()`/`PhysicsWorld` (M1.1.15) owns the same set. `rigid/root.zig` gained the `SolverConfig`/`solveRange` re-exports the harness consumes; `forge_3d/root.zig` pins the new test file. + - Acceptance scenarios (10): box rests without sinking (e=0, 300 ticks); bounces above / settles below the restitution threshold; restitution captured before warm start (resting box + real impact still bounces); inclined static box holds (μ>tanθ) / slides (μ Date: Fri, 24 Jul 2026 18:39:27 +0200 Subject: [PATCH 19/23] fix(forge): make combineFriction finite across the domain (M1.1.6) --- .../forge/forge_3d/rigid/contact_constraint.zig | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/modules/forge/forge_3d/rigid/contact_constraint.zig b/src/modules/forge/forge_3d/rigid/contact_constraint.zig index d6fe374..6de958f 100644 --- a/src/modules/forge/forge_3d/rigid/contact_constraint.zig +++ b/src/modules/forge/forge_3d/rigid/contact_constraint.zig @@ -43,10 +43,13 @@ const ContactManifold = narrowphase.ContactManifold(Real); // --- Material combine rules (Box2D/Jolt convention) -------------------------- -/// Combined friction for a contact pair: the geometric mean √(a·b). Both inputs -/// are non-negative (the `addBody` domain guard), so the product is ≥ 0. +/// Combined friction for a contact pair: the geometric mean √(a·b), computed as +/// √a·√b. Both inputs are non-negative (the `addBody` domain guard). The √-first +/// form is the identical geometric mean but stays finite across the whole finite +/// domain — it avoids the intermediate `a*b` overflow (→ inf, then `inf·0 = NaN` +/// in the friction clamp when λₙ = 0) and underflow of a naive `√(a·b)`. pub fn combineFriction(a: Real, b: Real) Real { - return @sqrt(a * b); + return @sqrt(a) * @sqrt(b); } /// Combined restitution for a contact pair: the maximum. A single bouncy surface @@ -303,6 +306,14 @@ test "combine rules: friction is the geometric mean, restitution is the max" { try testing.expectEqual(@as(Real, 0.8), combineRestitution(0.8, 0.2)); // symmetric } +test "combineFriction stays finite and positive across the domain" { + // The sqrt-first form avoids the intermediate a*b overflow/underflow of a + // naive √(a·b): at f32, a*b = 1e40 overflows (→ inf) and 1e-60 underflows (→ 0). + const big = combineFriction(1e20, 1e20); + try testing.expect(std.math.isFinite(big) and big > 0); + try testing.expect(combineFriction(1e-30, 1e-30) > 0); +} + test "tangent basis is orthonormal, right-handed, and deterministic" { const normals = [_]Vec3r{ vr(1, 0, 0), vr(0, 1, 0), vr(0, 0, 1), From 7460e8dca4afb2f48231a03d3cb1b2cca6c870a2 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 18:39:28 +0200 Subject: [PATCH 20/23] fix(forge): guard solver restitution threshold (M1.1.6) --- src/modules/forge/forge_3d/rigid/velocity_solver.zig | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/modules/forge/forge_3d/rigid/velocity_solver.zig b/src/modules/forge/forge_3d/rigid/velocity_solver.zig index f4613a4..36f8866 100644 --- a/src/modules/forge/forge_3d/rigid/velocity_solver.zig +++ b/src/modules/forge/forge_3d/rigid/velocity_solver.zig @@ -28,7 +28,7 @@ const ContactCache = cache_mod.ContactCache; /// Velocity-solver tuning. pub const SolverConfig = struct { /// Gauss-Seidel velocity iteration passes per tick (named for the M1.1.7 - /// `position_iterations` sibling). + /// `position_iterations` sibling). 0 is allowed — warm start only, no solve. velocity_iterations: u32 = 8, /// Restitution cutoff (m/s): a bounce is applied only when the pre-solve /// relative normal speed exceeds this — a PHYSICAL velocity constant (config @@ -107,6 +107,10 @@ pub fn storeContacts(gpa: std.mem.Allocator, cache: *ContactCache, constraints: /// island-additivity seam (M1.1.8 sorts constraints into contiguous per-island /// ranges — purely additive). pub fn solveRange(bm: *BodyManager, constraints: []ContactConstraint, from: usize, to: usize, cfg: SolverConfig) void { + // A finite, non-negative restitution threshold — a negative one would rearm + // the E4 defect (restitution bias applied to a separating contact). Zero + // `velocity_iterations` is legal (warm start only) and needs no guard. + std.debug.assert(std.math.isFinite(cfg.restitution_threshold) and cfg.restitution_threshold >= 0); var iter: u32 = 0; while (iter < cfg.velocity_iterations) : (iter += 1) { for (constraints[from..to]) |*c| { From 553cdd82f185a767da7798b7531e963798bc848a Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 18:39:30 +0200 Subject: [PATCH 21/23] fix(forge): retain contact pairs across separation in harness (M1.1.6) --- .../forge/forge_3d/tests/solver_test.zig | 88 ++++++++++++++----- 1 file changed, 64 insertions(+), 24 deletions(-) diff --git a/src/modules/forge/forge_3d/tests/solver_test.zig b/src/modules/forge/forge_3d/tests/solver_test.zig index 704d96b..dbb7aa2 100644 --- a/src/modules/forge/forge_3d/tests/solver_test.zig +++ b/src/modules/forge/forge_3d/tests/solver_test.zig @@ -11,13 +11,20 @@ //! (6) `storeContacts` (harvest) → `cache.endTick` (sort + swap) //! (7) broadphase proxy updates. //! -//! `computePairs` is moved-driven (it reports only pairs touching a proxy that -//! moved since the last call), so the harness keeps a PERSISTENT candidate set -//! (`active`, a sorted-deduped key list): it merges each tick's deltas and, after -//! the solve, retains only the pairs that produced a constraint (contacting -//! ones), so a resting contact — whose proxies stop moving and thus stop being -//! reported — stays live and keeps being solved (no sinking). A real `step()` / -//! `PhysicsWorld` at M1.1.15 owns the same persistent pair set. +//! `computePairs` is moved-driven with fat-AABB hysteresis (it reports a pair only +//! when a proxy moves enough to exit its fat AABB), so the consumer keeps a +//! PERSISTENT candidate set (`active`, a sorted-deduped key list) and merges each +//! tick's deltas into it. NORMATIVE retention rule for the M1.1.15 +//! `step()`/`PhysicsWorld` (b2ContactManager semantics): a pair is retained while +//! the two FAT broadphase AABBs overlap, dropped only when they separate — which +//! is safe, since any later re-contact first requires a proxy to exit its fat +//! AABB, which re-marks it moved and re-emits the pair. This test harness keeps +//! EVERY emitted pair (never drops) — a conservative superset of that rule, valid +//! in test: the narrowphase filters non-touching pairs (a pair with no manifold +//! produces no constraint), so a retained separated pair only costs a redundant +//! `collidePair`, never a wrong contact. Dropping a contacting pair on transient +//! separation (e.g. a small hop within the fat margin), by contrast, would lose +//! the contact until the body sank past the margin. const std = @import("std"); const config = @import("../config.zig"); @@ -128,10 +135,9 @@ const World = struct { try rigid.storeContacts(gpa, &self.cache, self.constraints.items); self.cache.endTick(); - // Persist only contacting pairs; separated ones drop (re-added by - // computePairs when their proxies move back into overlap). - self.active.clearRetainingCapacity(); - for (self.constraints.items) |c| try self.active.append(gpa, c.pair_key); + // `active` is NOT pruned — every emitted pair is retained (see the file + // header). The narrowphase filters non-touching pairs at `build`, so a + // retained separated pair costs only a redundant `collidePair`. // (7) broadphase proxy updates to the new poses. for (self.bodies.items) |b| { @@ -163,23 +169,57 @@ test "box dropped on a static ground comes to rest without sinking (e = 0)" { defer world.deinit(gpa); const box = try groundAndBox(gpa, &world, 2.0, 0); - // Settle (2 s). + // Run 5 s. Capture the height at the FIRST tick that produces a contact + // constraint (impact), then require the box never sinks below it afterwards — + // the AC's "penetration bounded by its at-impact value" (absolute non-recovery + // accepted until the NGS position solver, M1.1.7). + var y_impact: ?Real = null; + var min_after_impact: Real = 1e30; var t: u32 = 0; - while (t < 120) : (t += 1) try world.step(gpa); - const y_settled = world.bm.position(box).?.toArray()[1]; - - // Run 3 more seconds — the box must not sink further. - while (t < 300) : (t += 1) try world.step(gpa); + while (t < 300) : (t += 1) { + try world.step(gpa); + const y = world.bm.position(box).?.toArray()[1]; + if (y_impact == null and world.constraints.items.len > 0) y_impact = y; + if (y_impact != null and y < min_after_impact) min_after_impact = y; + } const y_final = world.bm.position(box).?.toArray()[1]; const v_final = world.bm.linearVelocity(box).?.toArray()[1]; - // Analytical rest y = ground_top(0.5) + box_half(0.5) = 1.0; the velocity - // solver leaves a small constant penetration (no NGS recovery until M1.1.7). - try testing.expect(y_final > 0.5); // never fell through - try testing.expect(y_final < 1.02); // resting at/near the contact - try testing.expect(y_final > 0.8); // penetration bounded (< 0.2) - try testing.expect(@abs(y_final - y_settled) < 1e-3); // no sinking after settling - try testing.expect(@abs(v_final) < 0.05); // at rest + // The box made contact and never sank below the at-impact height thereafter. + try testing.expect(y_impact != null); + try testing.expect(min_after_impact >= y_impact.? - 1e-3); + // Envelope: analytical rest y = ground_top(0.5) + box_half(0.5) = 1.0; a small + // constant penetration remains (no NGS recovery until M1.1.7). Never fell + // through, rests at/near the contact, penetration bounded, at rest. + try testing.expect(y_final > 0.8 and y_final < 1.02); + try testing.expect(@abs(v_final) < 0.05); +} + +test "small hop within the fat margin keeps the contact pair alive" { + const gpa = testing.allocator; + var world = World.init(vr(0, -9.81, 0), 1.0 / 60.0); + defer world.deinit(gpa); + const box = try groundAndBox(gpa, &world, 1.0, 0); // starts flush, e = 0 + + // Settle, then record the resting height. + var t: u32 = 0; + while (t < 120) : (t += 1) try world.step(gpa); + const y_rest = world.bm.position(box).?.toArray()[1]; + + // A small upward hop: apex ≈ v²/2g = 0.25/19.62 ≈ 1.3 cm, well within the + // broadphase fat-AABB margin (0.1 m). The box never moves far enough to be + // re-emitted by the moved-driven `computePairs`, so retaining the pair is what + // keeps the contact live and catches the box on its way down. Dropping the pair + // on separation would lose it until the box sank past the fat margin (~0.1 m). + world.bm.addImpulse(box, vr(0, 0.5, 0)); + var min_y: Real = y_rest; + t = 0; + while (t < 120) : (t += 1) { + try world.step(gpa); + const y = world.bm.position(box).?.toArray()[1]; + if (y < min_y) min_y = y; + } + try testing.expect(min_y >= y_rest - 0.02); } test "solve is deterministic across identical runs" { From 1e129d9e2503edade68cd341e0e939c00c190d9f Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 18:40:10 +0200 Subject: [PATCH 22/23] docs(forge): record E6 pre-merge fixes and deviation (M1.1.6) --- briefs/M1.1.6-solver-sequential-impulses.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index 7f895e2..00bcfe4 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -169,6 +169,12 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel - RED observed: the isotropy test (`"friction clamp is isotropic"`) failed under a per-axis BOX clamp (diagonal slide clamped to √2·μ·λₙ vs μ·λₙ along an axis), GREEN after the circular rescale. Second friction test: friction cancels sub-cone tangential sliding (contact-point tangential velocity → 0). - Verified: 150/150 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - Note for review: `velocity_solver.zig` = 592 lines (over the 500-line Review guideline) — the brief mandates warm start + iterations + friction + `SolverConfig` + range entry in this single file, plus inline tests (§13); conscious brief-driven overage, matching the `narrowphase.zig`/`broadphase.zig` precedent. Splitting the inline tests into a sibling `_test.zig` is available if preferred. +- 2026-07-24 18:39 — E6 pre-merge review fixes (`6c21786`, `7460e8d`, `553cdd8`), RED-first where applicable: + - P1 (`553cdd8`): the harness NO LONGER prunes `active` — every emitted candidate pair is retained (sort-deduped at merge); the narrowphase filters (a pair with no manifold produces no constraint). Header rewritten to the normative M1.1.15 retention rule (retain while the FAT broadphase AABBs overlap — b2ContactManager semantics; the harness keep-all is the conservative superset valid in test). RED-first test `"small hop within the fat margin keeps the contact pair alive"` (rest, +0.5 m/s hop with apex ≈ 1.3 cm < margin 0.1, 120 ticks): RED under the old prune (box sank ~0.1 m before re-emission), GREEN with retention. + - P2a (`6c21786`): `combineFriction` `@sqrt(a*b)` → `@sqrt(a)*@sqrt(b)` (identical geometric mean, finite across the domain; kills the `inf·0 = NaN` clamp path). RED-first test `"combineFriction stays finite and positive across the domain"` (RED: `inf` at f32 for `combineFriction(1e20, 1e20)`). See Recorded deviations. + - P2b (`7460e8d`): `solveRange` debug-asserts `isFinite(restitution_threshold) and >= 0` (a negative threshold rearms the E4 defect); `velocity_iterations == 0` allowed (warm start only) — documented on `SolverConfig`, no assert. + - P2c (`553cdd8`): the rest-without-sinking test strengthened to the frozen AC — capture `y_impact` at the first tick with a constraint, then assert `min(y) ≥ y_impact − 1e-3` over all 300 ticks; envelope asserts kept. + - Verified: 162/162 tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe; default `zig build` clean; `zig fmt --check` clean. - 2026-07-24 17:57 — E6 (harness + acceptance suite) (`9d13e91`): - `tests/solver_test.zig` (490 lines): `World` composes the normative per-tick cycle IN TESTS (computePairs → integrateVelocities → beginTick+build+prepare+warmStart+solveRange → integratePositions → storeContacts+endTick → broadphase proxy updates). Since `computePairs` is moved-driven, the harness keeps a PERSISTENT candidate set (sorted-deduped keys, merged each tick, pruned to contacting pairs post-solve) so resting contacts stay solved — a real `step()`/`PhysicsWorld` (M1.1.15) owns the same set. `rigid/root.zig` gained the `SolverConfig`/`solveRange` re-exports the harness consumes; `forge_3d/root.zig` pins the new test file. - Acceptance scenarios (10): box rests without sinking (e=0, 300 ticks); bounces above / settles below the restitution threshold; restitution captured before warm start (resting box + real impact still bounces); inclined static box holds (μ>tanθ) / slides (μ restitution_threshold`"; implemented instead as `v_n⁻ < −restitution_threshold` (approaching-only). Reason: the frozen `|v_n⁻|` form yields a NEGATIVE bias target on a capture-time SEPARATING contact (`v_n⁻ > +threshold`), which under-enforces non-penetration when the body actually approaches at solve time (the negative target either suppresses the corrective impulse via the `λₙ ≥ 0` clamp or, for a faster approach, lets the bodies close at the bias speed). The corrected condition keys the bias on the approach direction, so the bias is never negative. Claude.ai spec defect, ratified in the E4 review (Guy, 2026-07-24); the frozen brief text is left as-is pending a Claude.ai round-trip / KB re-upload. ## Blockers encountered From 672ec6dd130478914c9874eb1e97a79ee59ea5c7 Mon Sep 17 00:00:00 2001 From: Guy Senpai Date: Fri, 24 Jul 2026 19:40:36 +0200 Subject: [PATCH 23/23] docs(forge): close M1.1.6 and patch CLAUDE.md state (M1.1.6) --- CLAUDE.md | 10 ++++++---- briefs/M1.1.6-solver-sequential-impulses.md | 14 +++++++------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index ce128a9..3823596 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.3-hf-epa-frame-dependence` — narrowphase EPA deep-path order-dependence hotfix — code-complete, PR open. M1.1.5 (integration) CLOSED, merged + tagged `v0.11.5-integration-euler`. | +| Current milestone | M1.1.6 — Sequential Impulses contact solver — code-complete, PR open. The standalone hotfix `m1.1.3-hf-epa-frame-dependence` is CLOSED, squash-merged to `main` (untagged by hotfix policy, commit `7e63912`). | | Last released tag | `v0.11.5-integration-euler` (posted by Guy after merge) | -| Active branch | `phase-1/forge/epa-frame-dependence-hotfix` (PR open, not merged) | -| Next planned milestone | 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–M1.1.5 CLOSED (tags `v0.11.0`–`v0.11.5`); the standalone hotfix `m1.1.3-hf-epa-frame-dependence` is the CURRENT milestone (see Current milestone row), untagged by hotfix policy. 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`. | +| Active branch | `phase-1/forge/solver-sequential-impulses` (PR open, not merged) | +| Next planned milestone | M1.1.7 — NGS position solver (stacking stability; consumes the per-point `penetration` the M1.1.6 constraints carry but do not use — the velocity solver has NO positional bias by design). Then M1.1.8 (islands + sleep, additive over the M1.1.6 `solveRange` index-range seam). M1.1.0–M1.1.6 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 @@ -62,6 +62,7 @@ knowledge base — see § Quick links spec. | `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. | +| `v0.11.6-solver-sequential-impulses` | 2026-07-24 | M1.1.6 — Sequential Impulses contact solver | Seventh M1.1 sub-milestone — the FIRST consumer of the M1.1.3/4 manifolds + frame-stable `feature_id`s. New `forge_3d/rigid/` package (spec §1.2 rigid branch; `pipeline/` stays branch-shared): `contact_constraint.zig` (combine rules — friction geometric mean computed `√a·√b`, finite across the whole domain (a naive `√(a·b)` overflows/underflows), restitution `max(ea,eb)`; trig-free dominant-axis tangent basis, fixed tie-breaks; `build` over canonical broadphase pair keys → `collidePair` (never bare `collide`) → sort by pair key; `prepare` precomputes lever arms, world inverse inertias, normal + two independent tangent effective masses, captures `v_n⁻`; true-zero skip `inv_mass_a+inv_mass_b == 0`), `contact_cache.zig` (double-buffered sorted flat arrays, key `(pair u64, subshape_id u32 = 0, feature_id u32)`, value `λn` + WORLD-space tangent impulse Vec3, binary-search matching per individual feature, hit/miss counters, `beginTick`/`endTick` = sort + swap, implicit eviction, zero hash containers), `velocity_solver.zig` (`SolverConfig{velocity_iterations=8 — 0 legal, warm-start only; restitution_threshold=1.0 m/s — entry-asserted finite ≥ 0, a negative value would rearm the E4 defect}`; `warmStart` reprojects the cached world tangent onto the new tangent plane + clamps to `μ·λn` — compensates dominant-axis basis flips; Gauss-Seidel `solveRange(from,to)` — the M1.1.8 island seam — normal accumulated clamp `λn ≥ 0`, restitution bias `−e·v_n⁻` iff `v_n⁻ < −threshold` (approaching-only — recorded deviation from the frozen brief's `\|v_n⁻\|` form, which under-enforces non-penetration; KB spec patch pending), friction after normal per point with a CIRCULAR clamp `‖(λt1,λt2)‖ ≤ μ·λn` — isotropic, basis-independent). `Body` gains `friction`/`restitution` SoA columns (were silently dropped from `BodyDescriptor` since M1.1.0) + domain debug asserts (`friction ≥ 0` finite, `restitution ∈ [0,1]` finite). `integrate` split into `integrateVelocities`/`integratePositions` (`integrate` = exact composition, bit-identical, M1.1.5 tests untouched) so the solve sits between the passes. NO positional bias of any kind (no Baumgarte/split impulse — NGS is M1.1.7); resting penetration is bounded, constant, unrecovered until M1.1.7 (pinned at its at-impact value across 300 ticks). Deep-band contacts at penetration ≈ 0 are solved normally (never rejected on penetration). Test harness composes the full pipeline IN TESTS ONLY, with a persistent candidate-pair set retained while the FAT broadphase AABBs overlap (`computePairs` is moved-driven with fat-AABB hysteresis; the test harness keeps every emitted pair — a conservative superset) — retention rule NORMATIVE for the M1.1.15 `step()`/`PhysicsWorld`; a purge on transient sub-margin separation loses the contact until the body sinks past the fat margin (pinned by the small-hop regression test, pre-merge Codex P1). 162/162 tests green f32 + `-Dphysics_f64=true`, debug + ReleaseSafe. OUT (later): NGS (M1.1.7), islands/sleep/parallel (M1.1.8), joints, `step()`/`PhysicsModule` (M1.1.15), queries (M1.1.9–10), CCD/speculative, per-shape `PhysicsMaterial` (precedence decision recorded: material overrides body when present), clip-recentring precision boundary (future hotfix), `forge_2d`. | ### Hotfixes (untagged) @@ -122,6 +123,7 @@ Hotfix milestones are merged to `main` without a tag (Guy decision, - **M1.1.4 — generic EPA deep-rotated frame-dependence (M1.1.3 DEFECT) — RESOLVED (hotfix `m1.1.3-hf-epa-frame-dependence`, before M1.1.6).** `collideOrderedGeneric` (GJK/EPA in the frame of A) was order-dependent for a DEEP, ROTATED convex pair even at 1:1 aspect (violating the `engine-physics-forge.md §3` order-independence contract), on the residual capsule/box + sphere/capsule generic path (box/box was already neutralised by the M1.1.4 SAT fast path). Resolved by three fixes, all in the hotfix: (1) `epa.zig` expansion robustness — coplanar-tolerant visibility, connected visible region (flood-fill from the closest face), transactional commit-on-success expansion, skip-face non-convergence semantics, a realizable exhaustion fallback, and a flat-face terminal-face reconstruction restricted to the converged plane; (2) the E3 intrinsic point⊖segment world normal (segment-owner's world rotation + explicit ownership sign → exact bit negation across orders); (3) **RD-4** — a `gjk.zig` post-loop deep band `if (dist <= contact_margin) deep`, reusing the existing `conv_k` accumulated-rounding bound (`noise_k` untouched, so the P1b/P2 calibration is preserved), closing a deep/shallow classification stall that mis-classified a genuinely-deep overlap `.shallow` at dist≈0 (arm64 f32-unit / f64-×0.01, scalar/scale/order-dependent). Validated by a 390-config × 2-order × 3-scale order-equivalence sweep with mandatory SAT classification + EPA diagnostics (no corrupt path fires), f32 + `-Dphysics_f64=true`, debug + ReleaseSafe, arm64 + x86-64 CI. See `briefs/m1.1.3-hf-epa-frame-dependence.md`. - **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. +- **M1.1.6 KB spec patch pending (engine-physics-forge.md)**: material combine rules (`combineFriction = √(μa·μb)`, computed `√a·√b`; `combineRestitution = max(ea,eb)`) + material domains (`friction ≥ 0` finite, `restitution ∈ [0,1]` finite) + body-vs-`PhysicsMaterial` precedence (material overrides body when present; implementation lands with per-shape materials) + a Sequential Impulses solver paragraph (normative per-tick cycle incl. the fat-AABB-overlap pair-retention rule, prepare→warmStart→iterate→storeCache order, restitution approaching-only `v_n⁻ < −threshold` — corrects the M1.1.6 frozen brief's `|v_n⁻|` form —, circular friction clamp, world-space tangent cache + reprojection, cache key incl. `subshape_id = 0` reserved for Compound/Mesh). Claude.ai produces the re-uploadable file; Guy re-uploads. ## Non-negotiable rules @@ -255,4 +257,4 @@ The `briefs/` directory is the source of truth for milestone state. The brief's --- -Last updated: 2026-07-23 +Last updated: 2026-07-24 diff --git a/briefs/M1.1.6-solver-sequential-impulses.md b/briefs/M1.1.6-solver-sequential-impulses.md index 00bcfe4..ee5b5a5 100644 --- a/briefs/M1.1.6-solver-sequential-impulses.md +++ b/briefs/M1.1.6-solver-sequential-impulses.md @@ -1,12 +1,12 @@ # M1.1.6 — Forge 3D: Sequential Impulses contact solver -> **Status:** PLANNED +> **Status:** CLOSED > **Phase:** 1 > **Branch:** `phase-1/forge/solver-sequential-impulses` > **Planned tag:** `v0.11.6-solver-sequential-impulses` > **Dependencies:** M1.1.0–M1.1.5, hotfix `m1.1.3-hf-epa-frame-dependence` (merged, untagged per 2026-07-14 hotfix policy) > **Opened:** 2026-07-24 -> **Closed:** — +> **Closed:** 2026-07-24 (code-complete, PR open; merge + tag are Guy's) --- @@ -192,8 +192,8 @@ Seventh core sub-milestone of the M1.1 rigid arc (M1.1.0–15). Delivers the vel ## Closing notes -- **What worked:** -- **What deviated from the original spec:** -- **What to flag explicitly in review:** -- **Final measurements:** -- **Residual risks / tech debt left intentionally:** +- **What worked:** The gate decomposition (E1 material columns + integration split → E2 constraint setup → E3 cache + warm start → E4 normal + restitution → E5 friction → E6 harness + acceptance) was clean and each gate closed RED-first. The integration split composed bit-identically to the M1.1.5 monolithic pass (existing tests untouched). The M1.1.3/4 frame-stable `feature_id` contract fed the warm-start cache directly (per-feature binary-search matching); the world-space tangent cache + reprojection made warm start robust across tangent-basis flips. Determinism by construction held (sorted flat arrays, zero hash containers) — the two-run 300-tick bit-exact test passed first try. +- **What deviated from the original spec:** Two FROZEN-SECTION deviations, both recorded above and ratified in review: (1) restitution condition `|v_n⁻| > threshold` → `v_n⁻ < −threshold` (approaching-only; the `|v_n⁻|` form applied a negative bias to a capture-time separating contact, under-enforcing non-penetration — Claude.ai spec defect); (2) `combineFriction` `@sqrt(a*b)` → `@sqrt(a)*@sqrt(b)` (identical geometric mean, finite across the domain; kills an `inf·0 = NaN` clamp path). Both need the KB spec patch (see the new Open-decisions bullet). Two files carrying the E1 acceptance tests (`tests/integration_test.zig`, `tests/body_manager_test.zig`) are not in the brief's "Files to create or modify" list but are named by the Acceptance criteria — touched per the AC. +- **What to flag explicitly in review:** (a) The harness keeps a PERSISTENT candidate-pair set (never pruned) because `computePairs` is moved-driven with fat-AABB hysteresis; the retention rule (retain while the FAT AABBs overlap) is NORMATIVE for the M1.1.15 `step()`/`PhysicsWorld` — flagged in the harness header, the tag row, and pinned by the small-hop regression test (pre-merge Codex P1). (b) `velocity_solver.zig` = 592 lines, over the 500-line Review guideline (brief-mandated single-file solver + inline tests, matching the `narrowphase.zig`/`broadphase.zig` precedent). (c) `forge_3d/root.zig` exposes `rigid` as a namespace (`pub const rigid`) rather than flattening each symbol. +- **Final measurements:** 162/162 forge-3d tests green in f32 AND `-Dphysics_f64=true`, debug AND ReleaseSafe. Full `zig build test` (debug + ReleaseSafe) green via the pre-push gate. `zig build` clean; `zig fmt --check` clean. No benchmarks mandated (absolute perf targets land with the full `step`, per the M1.1.3/4/5 precedent). +- **Residual risks / tech debt left intentionally:** No positional recovery (velocity solver only) — resting penetration stays constant at its at-impact value until the NGS position solver (M1.1.7), which consumes the per-point `penetration` these constraints already carry. Contact-point position noise far from the origin is the inherited M1.1.4 narrowphase precision boundary (f64 build covers >10 km worlds; clip-recentring is a future upstream hotfix). The KB `engine-physics-forge.md` spec patch (combine rules, material domains, body-vs-material precedence, the SI solver paragraph) is pending (Claude.ai-produced, Guy re-uploads).