# Lodestone — Stepless Rigid Body Physics: Design Notes & Post-Mortem **Status: experimental, shelved.** This document explains what Lodestone was trying to be, how the pieces fit together, and — more importantly — *why the central goal runs into obstructions that are fundamental rather than incidental*. Read this before trying to "finish" Lodestone or the Annulus swept-collision code it stands on. The short version: the ideas are good and parts of the machinery are genuinely worth keeping, but the headline ambition, an exact, stepless, penalty-free rigid body engine, is somewhere between really, really hard and impossible. ## 1. The goal Most production physics engines are **time-stepping**: they advance the world by a fixed timestep `h`, detect overlaps, and push penetrating bodies apart with stabilization terms (Baumgarte feedback, position projection, split impulses). Those depenetration terms are *fudges*. They are not derived from the mechanics, but are controllers bolted on to fight the error that fixed-step integration and discrete collision detection introduce. Lodestone was an attempt at the opposite: a **stepless, event-driven** engine. 1. Integrate every body forward along its *exact* rigid trajectory. 2. Find the exact **time of impact (TOI)** of the next contact event, analytically. 3. Resolve that single event with an exact impulse or contact force derived from the mechanics. 4. Repeat, advancing time from one exact event to the next — never a fixed step, never an overlap to clean up, never a penalty term. The appeal is obvious: if you never take a blind step, you never penetrate, so you never need a fudge to un-penetrate. Contacts are caught precisely at the instant `gap = 0`, and resolved with physics rather than feedback gains. This document is the record of *why that appeal is a trap*. ## 2. Architecture (the intended pipeline) The stack is three layers. From the bottom: **Continuous collision detection (Annulus.SweptCollisionDetection).** For a pair of moving bodies, define the scalar **gap function** `g(t)` — the signed separation along the contact normal. A contact event is a root `g(t*) = 0`. The primitives build up from point-vs-line (`LinevsPoint.cs`), through segment-vs-point and polygon-vs-point, to polygon-vs-polygon (`PolygonvsPolygon.cs`). `LinevsPoint.SeparationCalculator` constructs `g(t)` and its exact first and second derivatives in closed form. **Root finding (Annulus + Azimuth.Approximation).** Because `g(t)` has no closed-form roots in general (see §4), it is approximated over the impact window by a **Chebyshev interpolant** (chebfun-style spectral approximation, `VerticalDistanceChebyshevBasis`), whose roots are extracted as the real eigenvalues of the **colleague matrix** (`Azimuth.SparseLinearAlgebra.ColleagueMatrix`, the Chebyshev-basis analogue of a companion matrix), then polished with **Newton's method** using the analytic derivatives, falling back to **bisection** (`RefineRoots`). Each root is classified as *crossing* or *grazing* by the **parity of its multiplicity** (`BuildHits`: odd → the function changes sign → penetration; even → it touches and returns → graze). **Contact resolution (Lodestone).** A root becomes a `ContactPoint`, which becomes a `NonPenetrationConstraint` carrying a velocity Jacobian `J` (`NonPenetrationConstraint.cs`). At each event, `Universe.SolveIslandConstraints` gathers the active constraints for a connected island of bodies and hands them to `GlobalSolver.SolveIsland`, which forms and solves a **bounded linear least-squares / LCP**: ``` minimize ‖ b − (J M⁻¹ Jᵀ) λ ‖ subject to λ ≥ 0 (contacts only push, never pull) where b = (1 + restitution) · (J · v) (J·v is the relative normal velocity = g′) ``` This is a genuinely penalty-free formulation: `λ` are the exact impulses that make the contact non-approaching in a least-squares sense, with no gains and no spring. The `SolverLevel` enum (`Solvers/Solver.cs`) names two intended regimes — `VELOCITY` (impulses, for impacts) and `ACCELERATION` (forces, for resting contact) — matching the two-level structure of nonsmooth contact dynamics. Note only `VELOCITY` was ever wired up. The integration loop (`Universe`) advances bodies along their rigid trajectories, stops at the next TOI, solves, applies the velocity delta, and continues. ## 3. The contact classification, and why it was hard The swept-collision layer tries to label each event with a `TransitionState`: `Entering`, `Exiting`, `GlancingInside`, `GlancingOutside`, `GlancingParallel`, `GlancingParallelEntering/Exiting`, `FeaturesAligning`, `Stationary`. Lodestone's `Universe` branches on these (`Entering` → build a constraint; `Exiting`/`GlancingOutside` → treat as separating). **The idea is sound; the specific taxonomy is not a partition.** Two independent things are being asked at every contact, and the enum fuses them into one flat value: - **Feature axis** — *what* touches: vertex-vertex, vertex-edge, edge-edge. (`GlancingParallel` is essentially "the pair is edge-edge" — a feature fact.) - **Kinematic axis** — *how* they touch: approaching, separating, or tangential. These are orthogonal, and the kinematic axis has a clean, well-founded definition: **the sign of the first non-vanishing derivative of the gap** at the TOI. ``` g′ < 0 → Entering (approaching) g′ > 0 → Exiting (separating) g′ = 0, g″ > 0 → GlancingOutside (touches, curves away — a graze) g′ = 0, g″ < 0 → GlancingInside (touches, curves in — penetrates) first nonzero derivative of odd order → crossing (penetration) of even order → graze g ≡ 0 on an interval → resting / sliding contact ``` That is a *total order*: the type is (order n, sign) of the first nonzero derivative. `Inside` / `Outside` / `Parallel` are not arbitrary categories — they are successive terms of the Taylor expansion of the gap. The mistake in the code was to flatten this order together with the feature axis into non-mutually-exclusive enum values and then force a single label, which requires tie-breaking between things that are simultaneously true. That is the source of the branchy, never-quite-converging classification logic (the layered `Transition` rewrite passes in `PolygonvsPolygon.FindIntersectionsInWindow`). Note, **`g′ = g″ = 0` does not mean "resting."** It can be third-order *penetration* (see §4). Resting is the much stronger condition `g ≡ 0`: *all* derivatives vanish over an interval, not merely the first two. ## 4. The fundamental obstructions These toy problems are fundamental issues with the approach above, and likely unresolvable. ### 4.1 Rotation makes the gap transcendental Pure *translation* keeps `g(t)` polynomial and cheap to solve. The moment any body **rotates, even at constant angular velocity**, the position of its vertices goes as `cos(θ₀ + ωt)`, `sin(θ₀ + ωt)`, and the gap becomes an amplitude-modulated sinusoid: ``` g(t) = [leverArm]·(sin,cos)(Δθ + Δω·t) + (p + v·t + ½a·t²)·(−sin,cos)(θ₁ + ω₁·t) − elevation ``` This is a transcendental function. It has no closed-form roots, and (critically!) **no finite Taylor truncation is exact.** All the derivative-based reasoning in §3 is correct in principle, but the object it reasons about is an infinite series. This is why the root finder must be a spectral *approximation* (Chebyshev) rather than an algebraic solver: there is nothing to solve in closed form. ### 4.2 Unbounded root multiplicity (high-order tangency) Because `g(t)` is transcendental, contact events can be roots with arbitrarily high multiplicity. The canonical example: a downward-pointing vertex on a body that spins at `ω = 1` while translating straight down at unit speed traces ``` V(t) = (cos t, sin t − t), so against the line y = 0: g(t) = sin t − t g(0) = 0, g′(0) = 0, g″(0) = 0, g‴(0) = −1 ``` Near contact, `g(t) ≈ −t³/6 < 0`: the corner touches the surface tangent to **second order** (neither velocity nor acceleration would cause it to penetrate the surface) and then screws straight through it. This is exactly the "fulcrum" scenario: a box corner acting as an instantaneous centre of rotation, but pushed one order deeper. A classifier that stops at velocity + acceleration sees a graze, when the truth is penetration. With two mutually rotating bodies (a sum of sinusoids at different frequencies) you can drive still more leading derivatives to zero. **There is no bound on the order of tangency** achievable between two rigid bodies. ### 4.3 High-order roots are numerically ill-posed — and the classifier reads the fragile quantity A root of multiplicity `m` is conditioned like `ε^(1/m)`: with a Chebyshev noise floor around `1e-14`, a triple root is located to only `~2e-5`, and the function sits at the noise floor across a whole neighbourhood. The eigenvalue solver returns *some* cluster there — three real roots, one real plus a discarded complex pair, or two — essentially at the mercy of rounding. The damning part is that the entering-vs-grazing decision (`BuildHits`) is made from the root's **multiplicity parity** — and multiplicity is the *first* thing floating point destroys near a multiple root. A true triple root (odd → should be **Entering / penetration**) can numerically resolve as a double root (even → reported as a harmless **graze**). That is a **silently missed penetration**, the worst failure a collision engine can have. So the elegant "parity of the first nonzero derivative" criterion is exactly the criterion floating point cannot preserve where it matters. ### 4.4 The solver is blind above its own order The velocity-level LCP resolves only the first derivative. Its driving term is `b = (1 + e)·(J·v) = (1 + e)·g′`. At any contact where `g′ = 0` (relative normal velocity zero), the right-hand side is zero, and the least-squares minimum-norm solution is **λ = 0**, meaning no impulse, regardless of `g″`, `g‴`, or whether the bodies are actually about to penetrate. It is not that `g′ = g″ = 0` is special to the *solver*; the solver is blind the moment `g′ = 0`. The high-order root is just the case where that zero impulse is *wrong*. An acceleration-level solve (the unbuilt `ACCELERATION` path) would catch `g′ = 0, g″ < 0` — but not `g′ = g″ = 0, g‴ < 0`. It only moves the wall one derivative up. The solver hierarchy has the *same unbounded-order structure* as the root classifier. Both bottom out at a fixed order while the hard cases live above it. ### 4.5 Impulses fix velocity, not position An impulse changes velocity; it has no mechanism to change position. So an impulse can make two bodies stop *approaching* (`g′ ≥ 0`), but it **cannot remove existing penetration** (drive a negative gap back to zero). Depenetration is a zeroth-order/positional operation, and momentum-level mechanics only speaks first order. Therefore *any* operation that removes an existing overlap (Baumgarte, projection, extra restitution, position-based correction, etc.) is necessarily **extra-physical**. Newtonian rigid dynamics simply contains no operator for "un-overlap two bodies." Once §4.2–4.4 have let a penetration slip through, getting out of it is always a fudge. ### 4.6 Painlevé: rigid bodies + Coulomb friction are inconsistent This is the deepest obstruction, and the one that converts "very hard" into "provably impossible." At the acceleration level, a single frictional contact reduces to `g̈ = a·λₙ + b`, where `a` couples normal force to normal acceleration *through friction*. Normally `a > 0` (push harder, separate faster). For large friction coefficient and certain geometries, friction feeds the tangential force `μλₙ` back into normal acceleration with the wrong sign and **`a < 0`**. Then with the contact tending to penetrate (`b < 0`): ``` λₙ = 0 branch → g̈ = b < 0 (penetrates — infeasible) g̈ = 0 branch → λₙ = −b/a < 0 (pulls — violates λₙ ≥ 0) ``` **No solution exists.** This is the Painlevé paradox, and its signature failure is *non-existence*, not mere non-uniqueness. A tempting fix is to say "when the LCP has infinitely many solutions, pick the one minimizing the sum of squares of the forces", and this is a legitimate and standard regularizer for *statically indeterminate* contact (e.g. a four-legged table), but it is powerless here: you cannot least-squares-select from an empty solution set. And even in the isolated-multiplicity form of Painlevé, minimum-norm is not physically valid, because **Coulomb friction is non-associated** (there is no maximum-dissipation potential that makes norm-minimization the correct selector), and the distinct solution branches are genuinely different physical outcomes chosen by how the system passes *dynamically* through the singularity. There is no instantaneous, finite-force, penalty-free resolution, because in these configurations the rigid + Coulomb model has no consistent behaviour. ### 4.7 Zeno An event-driven scheme stalls on infinitely many events in finite time — a bouncing body with restitution `< 1` collides infinitely often as it settles, and the stepless loop never advances past the accumulation point. ## 5. The reframing: rigidity is the fudge Every obstruction above traces to one fact: **the rigid-body + Coulomb-friction + instantaneous-contact model is internally inconsistent.** You can't have all three of these requirements in a consistent model. ## 6. Why the field time-steps The mainstream answer is **velocity-level time-stepping**. Instead of resolving each event exactly and instantaneously, these schemes advance by a step `h` and impose the contact conditions as a complementarity problem over the step, in terms of **impulses** rather than forces. The crucial property: this formulation is **provably solvable** (Painlevé configurations included) precisely *because* an impulse-over-a-step can absorb exactly what an instantaneous finite force cannot represent. The impulsive resolution that Painlevé demands is built in automatically, and Zeno is handled by integrating through the accumulation rather than chasing it. ## 7. Where the code lives - **Gap function & TOI:** `Annulus/SweptCollisionDetection/LinevsPoint.cs` (`SeparationCalculator`), `SegmentvsPoint.cs`, `PolygonvsPoint.cs`, `PolygonvsPolygon.cs`. - **Spectral root finding:** `Azimuth/Approximation/Chebyshev.cs`, `Azimuth/SparseLinearAlgebra/ColleagueMatrix.cs`. - **Contact & solver:** `Lodestone/Universe.cs`, `Lodestone/Solvers/GlobalSolver.cs`, `Lodestone/Solvers/Solver.cs` (the `SolverLevel` enum), `Lodestone/Constraints/`. - **Disabled tests:** `Annulus.UnitTests/SweptCollisionDetection/PolygonvsPolygonTests.cs` — the `GlancingHits.SquareSpinning` and `GlancingHits.SquaresSlideBack` cases are commented out; they fail because they demand a single unambiguous label for the genuinely multi-valued / high-order configurations described in §3–§4.