@holotope/physics

Headless higher-dimensional mechanics for Holotope. The first Stage C release implements 4D convex mass properties, principal-plane inertia, dynamic RigidBody4 state, a momentum-primary ballistic PhysicsWorld4, and fixed-step interpolation into the renderer-neutral ObjectN scene graph. The current Stage D layer adds support/GJK queries, complete oriented-hyperbox contact patches, warm-started contact response with a coupled R4 tangent friction ball, a coupled four-coordinate bilateral point joint, scalar rigid Jacobian rows with generalized-force bounds, rigid distance equalities, two-guardian distance intervals, force-limited distance motors, and deterministic mixed-shape collider/body orchestration. Its rotational foundation also exposes paired-bivector coordinates, branch-aware relative SO(4) logarithms, analytic exponential/logarithm Jacobians, and the exact angular-velocity operator norm. A common one-to-six-row equality-block solver now serves point joints and three genuinely R4 rotational policies: preservation of one oriented material direction with its SO(3) stabilizer free, and preservation of an ordered two-frame with one complementary SO(2) rotation free, or preservation of a complete relative material frame with no rotational stabilizer free. Candidate generation is dimension-independent and includes exhaustive and temporally coherent sweep-and-prune providers; static and linearly swept AABBs share the same candidate contract, while infinite planes remain in an explicit exhaustive boundary lane. A capability-aware dispatcher distinguishes general distance, rounded shallow contact, bounded general R4 EPA penetration, complete vertex-polytope and exact hyperbox deep manifolds, analytic N-ball/N-ball and N-ball/hyperplane deep contact, exact R4 glome/hyperbox, hyperbox/hyperplane, and general vertex-polytope/hyperplane contact, and unsupported requests. Dimension-independent conservative advancement adds compact/compact linear casts, while compact/infinite-plane casts are analytic. R4 also has explicit constant-generator rigid trajectories and conservative compact/compact and compact/plane casts whose angular closing bound uses the exact SO(4) operator norm.

A separate XpbdConstraintSolverN supplies an auditable dimension-generic Float64 position-level kernel for compliant scalar relations. Equalities are unbounded and declared C(x) >= 0 inequalities project total multipliers onto the non-negative ray. Results expose the total XPBD multiplier, signed force estimate, raw compliant residual, and projected KKT residual, while exact RN distance, unsigned intrinsic simplex measure, and signed full-dimensional simplex measure constraints provide equality consumers and exact RN particle–hyperplane contact provides the first inequality consumer. This does not replace or silently couple to the velocity-level R4 rigid constraint solver. XpbdWorldN wraps that kernel in explicit RN point-mass prediction, velocity reconstruction, ordered post-projection velocity responses, force accumulation, substeps, read-only accepted-state guards, bounded adaptive retry, ownership checks, and atomic world-step rollback. Adaptive stepping retries only typed guard rejection and reports every attempted subdivision; it does not mask ordinary errors or provide a continuous no-inversion proof. Its first responses provide exact particle–plane Coulomb friction over the complete RN tangent ball and named timestep-invariant exponential damping.

stepXpbdIncrementalPotentialWorldN() advances one authored XpbdWorldN through that transaction. The world is authoritative for dimension, particle order, gravity, and the conservative provider registry, so those four stop being repeated at every call site where they could drift from the scene being rendered. Step filters stay explicit because the base world owns no filter registry. Registered scalar constraints, velocity responses, state guards, and non-conservative force providers cannot be represented by this path and are named as configuration errors rather than skipped. The return carries the complete lower-level step alongside its diagnosis and an immutable registration snapshot — not a success boolean.

The separate incremental-potential reference step is transactional and keeps its minimizer base explicit. Its historical default is the inertial prediction; callers may instead select the previous live positions or opt into bounded feasible-prediction chord sampling. That recovery retains every accepted or typed-domain-refused trial and makes no depenetration, nearest-point, global feasibility, or performance claim. It exists to initialize open-domain objectives without turning a solver policy into hidden behavior.

Its stop test is authored, and the unit that test bounds is a choice the library cannot make on the scene's behalf. The packed objective is F(x) = 1/2||x - xPrediction||^2_M + deltaTime^2 * U(x), so a packed gradient entry carries masslength — forcetime^2, not force. The shipped gradientTolerance (default 1e-8) is an absolute bound on the norm of that packed gradient and therefore resolves forces only down to gradientTolerance / deltaTime^2. That floor rises as the timestep falls, so refining the step makes the criterion less sensitive to force, not more. Driven end to end through the public world step, one free unit-mass particle under a constant 1000 N force over a fixed 1e-3 s horizon reaches the exact backward-Euler answer of 1 m/s at deltaTime = 5e-6 and identically 0 m/s at 2e-6. A 2.5-fold refinement takes a 1000 N force from exactly right to nothing, and all 500 steps of the 2e-6 run report applied — converging at the warm start is a legitimate outcome, not a refusable condition, and no field said otherwise.

convergence: { kind: 'packed-gradient'; tolerance } names that legacy criterion explicitly; { kind: 'maximum-acceleration-residual'; tolerance } instead bounds max_i ||gradient_i|| / (mass_i * deltaTime^2) over the free particles, in length/time^2. Fixed particles are excluded because they hold no packed coordinate and their gradient is identically zero, so they can neither raise nor lower it. The option is available on both minimizeXpbdIncrementalPotentialN and the world step's minimization policy. Authoring gradientTolerance and convergence together throws before the problem is evaluated once, since the two carry different units and no reconciliation between them is defensible. Every terminal that evaluated the base carries convergence: { kind, tolerance, initialResidual, finalResidual } whichever criterion decided it, and retains gradientNorm on its initial and final evaluations. The one exception is initial-state-refused, where no evaluation was accepted: it carries the criterion and its tolerance but no residuals and no evaluations, so narrow that status away before reading either. Diagnosis gains matching convergenceKind, convergenceTolerance, convergenceResidualInitial, and convergenceResidualFinal facts with 'lower-convergence-tolerance' and 'timestep-independent-convergence' as levers — an author who never wrote a gradientTolerance is no longer told to lower one.

The union is discriminated rather than a second scalar because measurement found no criterion to prefer outright. Six candidates were ranked on delivered physical error against a closed-form minimizer over an eight-fold refinement at fixed authored tolerance, as delivered acceleration spread / position-error spread: packed gradient norm 45.43 / 1.48, force-scaled residual 1.0153 / 62.93, maximum acceleration residual 1.0169 / 62.83, mass-weighted residual 1.0153 / 62.93, relative residual 1.0120 / 63.14, per-particle position residual 60.61 / 1.72. There are two families and no criterion spans both, because they differ by exactly deltaTime^2: the packed norm holds position error, the acceleration residual holds acceleration. The acceleration criterion is therefore not better — it trades position stability for acceleration stability. Author it when the timestep may change and a force resolution has to hold; keep the packed norm for a fixed timestep, or when a per-step position residual is the quantity that matters. Among the acceleration-stable candidates the per-particle one is the only one that is also mass-aware and count-invariant, since a global force norm grows as sqrt(N) and lets a heavy particle hide a light particle's acceleration behind it.

Its authored obstacle terms include both an oriented infinite hyperplane and one finite persistent source simplex. The latter retains the closest barycentric source coordinate. Its line-segment, triangle, and tetrahedron queries decide affine rank, zero distance, and the closest active face exactly on the supplied Float64 values, then publish one coherent Float64 witness with outward-rounded error bounds. There are no geometric tolerance knobs. The unsigned barrier pairs that query with a conservative convexity/Lipschitz segment certificate. A source-indexed family now lifts that pair over dynamic bound vertices and a separate static simplex mesh: exhaustive swept-AABB rejection keeps possible, retained, exact-active, and blocking-pair evidence separate while presenting one stable provider and one paired filter to the solver. It is not inside/outside classification, moving-simplex contact, or a claim of mesh self-collision.

P56 extends contact from constrained points to constrained features: evaluateSourceSimplexPairDistanceN is the dimension- and arity-generic minimum distance between two finite source simplices, certified by a variational inequality against every input vertex, with witnesses as source-ordered barycentric coordinates on both sides. Its result union keeps the mathematics honest — separated-unique carries the measured uniqueness margin that justifies the envelope-form gradient; separated-multiple returns every tied optimal witness (parallel edges are the canonical case) and no gradient, because none uniquely exists; zero-distance is certified with no invented normal; indeterminate refuses with its own residuals. XpbdSourceSimplexPairBarrierN lifts the clamped-log law over that distance with forces distributed through the witness weights (net internal force and the RN antisymmetric first moment cancel for two moving sides), and its paired filter certifies segment prefixes by the two-sided Hausdorff/Lipschitz bound d(t) ≥ d(0) − t·(maxDispA + maxDispB) — a certified fraction, never a collision time. compileXpbdSourceSimplexPairBarrierFamilyN sweeps one such pair per source cell of a deforming group against one static feature; the summed energy's density is discretization-defined (a shared edge carries both adjacent cells' terms — measured at exactly 2×), which the family documents instead of averaging away. This closes the P53d boundary: a sheet triangle can now be held off an obstacle that pierces its interior while every vertex is legally separated. It is not self-contact or mesh–mesh CCD.

compileXpbdSourceSimplexMeasureBarrierN offers the other weighting of the same contact. Where the pair family carries one term per source cell — so a shared edge carries both adjacent cells' terms — this law carries the cell's reference k-measure once and averages a clamped-log barrier over k + 1 fixed interior nodes, so splitting a cell does not answer twice. It compiles one conservative provider and one paired step filter, for k = 1, 2, 3 (the range over which the exact point–simplex query publishes a direction enclosure), and a successful evaluation carries exactly potentialEnergy and forces — there is no inspection surface and no Layer-2 record.

Measure consistency is not invariance under subdivision, and the two are kept apart deliberately. The integrand is a nonlinear barrier of a distance field and the rule is a fixed finite quadrature: subdivision is exactly additive only when the sampled barrier is constant over the cell, and otherwise it moves the sample locations and changes the estimate — measured at about 27% for a tilted cell split in half and about 44% for an uneven split of a curved arrangement. The refinement sequence does converge to the continuum integral, measured at second order against an independent composite Gauss–Legendre reference, with the single-cell estimate about 28% below it; that is a measurement on a named fixture, and no truncation bound is proved or claimed. No portable timing or performance multiplier is claimed either.

The quadrature rule is not authorable through the public API: the compiled terms are frozen and hold every non-authorable value in closure, so there is no rule option to pass and no rule, reference measure, obstacle snapshot or conservative scale property to overwrite.

That is a statement about the public surface, and not a concealment claim. Same-realm JavaScript metaprogramming can observe otherwise-private arrays — numeric accessors installed on Array.prototype before compilation retain the static-obstacle snapshot, the fixed rule and a private particle partition, and a replaced inherited operation receives whatever is used as its receiver. What the law guarantees is a consequence boundary rather than concealment: those retained arrays are frozen, so once the intrinsic is restored they cannot be modified to change a later evaluation, and the per-call geometry handed to released code is freshly allocated, so retaining or mutating it changes nothing later either. Persistent state is read by index with counts carried separately, precisely so that reading it does not hand it to a replaceable function.

The provider's published particles are excluded from that guarantee by design: they are the caller's own live inputs, and moving them changes later evaluations, which is the point of a contact term that reads live state.

The companion filter is required, not optional — the law measures unsigned distance and has no notion of side, so without the filter a step can leap clean through the obstacle with both endpoints admissible. This is normal contact only; friction is the separate lagged pair-friction term.

The generic pair query is still an experimental surface. A later scale audit found that its Float64 rank, zero, and optimality bands can change a result under exact similarity transforms. Do not use its 0-simplex specialization as the point–simplex authority. evaluateExactPointSimplexResult and the point–simplex barrier/family use the exact-on-supplied-Float64 path instead; the moving simplex-pair replacement remains current research.

Higher-dimensional point–simplex barriers remain available through the legacy Float64 projector so existing RN experiments can still be inspected, but they do not carry pointSimplex exact-decision evidence and inherit that projector's tolerance boundary. They are not part of the exact claim above.

P57 adds the first dissipative contact term to the objective itself. XpbdSourceSimplexPairFrictionN freezes one lag at an accepted state — the certified contact frame, the source-ordered witness weights, and the paired barrier's own normal-force magnitude — and is then conservative for exactly that snapshot, so every line-search trial sees one consistent objective. Dissipation happens between accepted states, when the lag is refreshed; calling it a globally conservative force would be wrong, and the vocabulary says so. The tangent projector I − n nᵀ is applied directly, with no authored basis, so the term is dimension-generic; the regularized Coulomb law is C¹ with a force that stays linear through zero slip (u/‖u‖ is never evaluated) and satisfies ‖f‖ ≤ μ·λ_lag by construction rather than by clamping. Only a separated-unique pair may create a lag — tied witnesses, certified zero distance, uncertified comparisons and sub-minimum distances refuse by type. compileXpbdSourceSimplexPairFrictionFamilyN lifts it over a contact family with atomic consume/rollback, and states plainly that effective friction follows mesh topology (a shared-edge contact resists exactly 2× one cell; a four-cell refinement 4×) rather than averaging that away.

slipRegularization poses the same shape of choice as the stop test above. A bare number is a world length and is never reinterpreted as anything else; it normalizes to { kind: 'slip-length'; length } carrying that exact value. A fixed length does not survive timestep refinement. Per-step slip is ||tangential velocity|| * deltaTime, so once the slip falls inside the regularized branch the force is forceLimit * slip / length and goes as deltaTime, one step's impulse as deltaTime^2, and a fixed horizon of T/deltaTime steps totals T*deltaTime — friction vanishes under refinement. Measured over an eight-fold refinement on two scenes, the tangential impulse falls to 0.133 of its coarse value with a last-halving ratio of 1.98 against the 2.00 that scaling predicts, confirmed through two independent channels — force-side impulse and velocity-side energy — agreeing to 0.19%, and cross-checked against a momentum audit to 1e-4. { kind: 'slip-velocity'; velocity } resolves the length as velocity * deltaTime, which cancels deltaTime out of slip / length exactly; under the same refinement the impulse holds to 1.06 of its coarse value, last halving 0.99. It remains a smoothing scale and nothing more: a velocity-derived scale does not establish static friction and does not give the law finite-support retention.

The resolved length is frozen into the lag, so prepare takes { deltaTime } — required under a slip velocity and refused under a slip length, because supplying it under an authored length would suggest that length responds to the timestep. Freezing is load-bearing rather than incidental: conservativeness within one lag is what lets the Armijo search evaluate the term repeatedly, and a length that moved mid-solve would leave the search minimizing a function whose own shape changed under it.

Evaluations separate two axes that read like one. regime — 'sticking', 'transition', 'sliding' — is a statement about slip alone; a lag carrying no normal force still has a slip and still reports a regime. contactActive is exactly forceLimit > 0 and is what decides whether the term can exert any tangential force at all. The two are orthogonal, and neither may be inferred from the other: in the sheet probe 144 of 192 evaluations read 'sliding' while exerting exactly zero force, so a population statistic that does not split on activity is mostly reporting about terms that are not touching anything. Friction work is likewise measured rather than inferred from total-energy decay — an integrator loses energy at mu = 0, and the measured mu = 0 control drifts 0.0395%, which is 24% of the smallest signal it certifies.

This is a different mechanism from XpbdParticleHyperplaneFrictionN, which is a post-projection Coulomb velocity response for the projected-XPBD path. The two are not interchangeable, and the incremental-potential path still refuses velocity responses outright.

When the obstacle's cells are only a decomposition of one solid rather than independently meaningful features, the per-cell sum is the wrong composition: each cell's barrier pushes away from itself, so a point over a flat support accumulates a decomposition-dependent tangential force. XpbdParticleSourceConvexHullBarrierFamilyN is the set-shaped alternative — the convex hull of the obstacle vertices its source group selects, one certified closest-point query per bound particle, one force along the separation normal, and a witness retaining which authoritative source vertices support the closest feature. The cells select vertices; they are not summed, and concavities between the selected vertices are filled, so a non-convex obstacle needs explicitly managed convex pieces. The hull is static for the family's lifetime — coordinates are snapshotted at compilation and a moved source is refused, never followed — and proximity to a lower-dimensional hull is unsigned and two-sided, because such a set has no ambient inside. A distance query that cannot certify separation or intersection within its bounded budget surfaces as a typed closest-point-indeterminate refusal rather than an answer, and the paired filter certifies conservative prefixes with the same convexity/Lipschitz proof as the point–simplex specialization.

That candidate scan stays the default and the oracle. XpbdSourceSimplexAabbHierarchyN is an opt-in immutable AABB tree over the same static obstacle, selected by passing it as candidateHierarchy — never by mesh size or a mode string, and only when it indexes the same obstacle and group objects the family does. It changes which pairs are asked, not what a retained pair means: candidate identity and order are exactly the exhaustive ones, and the exact barrier and paired prefix filter still decide contact. Because it caches bounds at compilation it requires a static obstacle, snapshots the coordinates it indexed, and refuses a moved source by naming the vertex and axis rather than rebuilding itself. Its diagnostics are operation counts; on an obstacle it cannot separate, work is linear and the counts say so.

compileXpbdSourceSimplexCosineBendingFamilyN() adds source-retained extrinsic stiffness over adjacent simplices. It is a discrete cosine-fold stiffness, not a continuum shell calibration. The coordinate is c = -uA . uB over the shared-face conormals — the orientation-neutral cosine of the fold from flat, never a signed dihedral — and the energy 0.5 k (c - cRest)^2 is therefore quartic in the fold angle where a continuum bending energy is quadratic. It is not mesh-convergent: a fixed strip refined in place has its total fall as n^-2.99, so stiffness values are discretization-dependent and do not transfer to a refinement. At a flat rest the first derivative vanishes, so small folds produce a weak restoring force.

Only unit weighting exists. The gradient is closed-form over all d+2 hinge vertices and cancels the translation and rotation modes algebraically, so netForceResidual and rotationalFirstMomentResidual are roundoff-scale evidence to compare against a tolerance rather than quantities guaranteed to be bitwise zero. The family is first-order only, so Newton-CG refuses the mixture with named unsupported-provider evidence rather than dropping bending curvature silently. Its paired filter is not optional: a search segment can begin and end with valid hinges while passing through zero conormal height in between, so the filter reuses analyzeLinearSimplexMeasureN over each distinct source simplex to certify a conservative admissible prefix. That prefix is an intrinsic-rank certificate, not an exact collapse time.

World-frame angular momentum is authoritative. Free flight therefore does not numerically integrate a gyroscopic force or silently lose momentum; angular velocity is derived through the body's principal inertia each step and the orientation remains on Spin(4) through paired-quaternion normalization.

The static source-simplex hierarchy above is the only deformable-candidate spatial index: it covers one unmoving obstacle and is opt-in. Refit for moving obstacles, moving--moving candidate trees, moving infinite-plane pose policies, distance servos, rolling resistance, and sleeping are not yet part of this package. R4 Coulomb friction is represented by one rotationally symmetric three-dimensional tangent ball, never by three independent scalar clamps.

The dimensional boundary is explicit: particle XPBD, broadphase bounds, GJK, and linear casts have RN contracts where their names say N; rigid-body state, penetration/manifold generation, and contact response currently have R4 contracts. An N-dimensional query is therefore not evidence of an N-dimensional rigid response path.

import {
  ObjectN,
  SceneN,
  createHypercube,
  tetrahedralizeCuboidCells
} from '@holotope/core';
import {
  PhysicsWorld4,
  RigidBody4,
  RigidBodyObject4Binding,
  massPropertiesFromCellComplex4,
  rebasePositionsToPrincipalFrame4
} from '@holotope/physics';

const geometry = tetrahedralizeCuboidCells(createHypercube({ dim: 4 }));
const mass = massPropertiesFromCellComplex4(geometry);
const principalPositions = rebasePositionsToPrincipalFrame4(geometry.positions, mass);
const body = RigidBody4.fromMassProperties(mass);
const scene4 = new SceneN(4);
const object4 = new ObjectN(4);
scene4.add(object4);
const binding = new RigidBodyObject4Binding(body, object4);

new PhysicsWorld4().addBody(body).step(1 / 60);
const alpha = 0.5; // normally renderAccumulator / fixedStep
binding.capture().apply(alpha);
scene4.updateWorld();

A browser render loop should keep simulation time fixed and rendering time variable. The accumulator below is the complete handoff; the first animation frame has zero elapsed time, and PhysicsWorld4.step(0) is also defined as a no-op for clocks that forward that value directly. Seed previousTime from the first animation-frame timestamp as shown—an earlier performance.now() can be slightly newer than that timestamp and produce a negative first delta.

const fixedDt = 1 / 120;
let previousTime: number | undefined;
let accumulator = 0;

function frame(timeMilliseconds: number) {
  const elapsed = previousTime === undefined
    ? 0
    : Math.min((timeMilliseconds - previousTime) / 1000, 0.25);
  previousTime = timeMilliseconds;
  accumulator += elapsed;

  while (accumulator >= fixedDt) {
    world.step(fixedDt, 2);
    binding.capture();
    accumulator -= fixedDt;
  }

  binding.apply(accumulator / fixedDt);
  scene4.updateWorld();
  renderer.render(scene, camera);
  requestAnimationFrame(frame);
}
requestAnimationFrame(frame);

The 0.25 clamp prevents a backgrounded tab from demanding an unbounded catch-up burst. Forces and torques survive a zero-time no-op and clear only after a positive completed step.

The optimization path uses the same accumulator, with two differences. Its deltaTime must be strictly positive — a zero interval is not a physical optimization step, so the while guard is what skips it rather than a no-op inside the step — and a mathematical refusal is a typed result to read rather than an exception to catch, while a configuration problem still throws. XpbdWorldN.step() and stepAdaptive() also reject a zero interval, so the guard is the policy for both of an XpbdWorldN's paths; only the rigid PhysicsWorld4.step(0) above is a defined no-op. And an XpbdWorldN has two solver paths, so running both over one interval integrates that interval twice; pick one per frame.

import {
  stepXpbdIncrementalPotentialWorldN
} from '@holotope/physics';

function optimizationFrame(timeMilliseconds: number) {
  const elapsed = previousTime === undefined
    ? 0
    : Math.min((timeMilliseconds - previousTime) / 1000, 0.25);
  previousTime = timeMilliseconds;
  accumulator += elapsed;

  // The guard is the zero-interval policy: `accumulator >= fixedDt` is never
  // true for an idle frame, so no step of length zero is ever requested.
  while (accumulator >= fixedDt) {
    const advance = stepXpbdIncrementalPotentialWorldN({
      world: particleWorld,          // never particleWorld.step() as well
      deltaTime: fixedDt,
      stepFilters: contactTerms.stepFilters,
      warmStart: 'feasible-inertial-prediction'
    });
    if (advance.step.status !== 'applied') {
      // Nothing moved and nothing threw. The diagnosis names the condition
      // and the caller-controlled levers that legitimately address it.
      reportStall(advance.diagnosis.condition, advance.diagnosis.levers);
      accumulator = 0;
      break;
    }
    binding.writeSourcePositions();
    accumulator -= fixedDt;
  }

  scene4.updateWorld();
  renderer.render(scene, camera);
  requestAnimationFrame(optimizationFrame);
}
requestAnimationFrame(optimizationFrame);

PointJoint4 binds a body-local anchor to another body or a fixed world point. Resolve it inside the world's velocity-constraint callback and pass the result to PointJointSolver4; the solver exposes the complete 4x4 point response and solves all four bilateral coordinates as one block.

ConstraintRowSolver4 solves scalar rigid-Jacobian rows with optional minForce and maxForce. It converts those generalized-force bounds to impulse bounds using the substep duration, then projects every accumulated impulse. Omitting both bounds gives an unrestricted equality row. Results separate raw equality residual from same-sign projected KKT residual, so valid saturation is distinguishable from a row that has not converged. Aggregate coordinate impulse, error, and residual values are scale-dependent solver diagnostics, not physical totals across unlike rows.

DistanceCoordinate4 is the persistent anchor binding shared by three policies. DistanceJoint4 enforces a positive rest length. DistanceIntervalJoint4.constraints(dt) always returns stable minimum and maximum unilateral guardian rows. Inside [minLength, maxLength], their speed targets bound the next first-order position; outside, signed error produces recovery bias. Keeping both guardians in the solve lets them catch unsafe radial velocity introduced by motors or other rows during iteration. interval(dt) reports the currently observed crossing state for diagnostics only. DistanceMotor4 tracks radial speed with symmetric maxForce, with positive speed lengthening the coordinate. Place its row before both guardians when solving them together. The coordinate geometry is also exposed through evaluateDistanceCoordinateN() and evaluateDistanceConstraintN() for any VecN dimension. At exact coincidence, solver rows require an authored scalar direction branch and refuse transverse or negative-branch relative motion; diagnostics may still observe one-sided distance growth without manufacturing a solve gradient.

XpbdConstraintSolverN instead projects scalar relations over mutable RN point coordinates. One solve batch has an explicit dimension and initializes one total multiplier per constraint. Compliance is physical inverse stiffness and is scaled by 1 / dt^2 inside the update; results report the corresponding signed force and C + alpha/dt^2 * lambda residual. A greater-than-or-equal relation projects the total multiplier to lambda >= 0; its projected KKT residual treats valid positive slack as zero error. Custom evaluators are pure, dimension-checked functions with one gradient per unique point. Invalid batches restore every participating position. XpbdDistanceConstraintN reuses the same exact distance coordinate and coincidence-branch rule as the rigid adapter.

evaluateSimplexSquaredMeasureN() evaluates the intrinsic k-measure of any k-simplex embedded in RN from det(E^T E) / (k!)^2, together with Float64 ambient gradients. XpbdSimplexSquaredMeasureConstraintN constrains that squared coordinate directly, so its compliance units depend on k. The coordinate is translation- and rotation-invariant and needs no dimension- specific cross product. It is deliberately unsigned: it preserves measure magnitude but is not an inversion barrier. Cofactor gradients remain finite at singular Gram matrices; a collapsed simplex whose first derivative vanishes reports no-dynamic-response rather than receiving an invented recovery normal.

evaluateOrientedSimplexMeasureN() instead evaluates det([x1 - x0, ..., xN - x0]) / N! for exactly N + 1 points in R^N. Its cofactor gradients transform covariantly under SO(N), while reflection or an odd vertex permutation reverses the scalar sign. XpbdOrientedSimplexMeasureConstraintN can therefore preserve and report material-cell orientation as well as magnitude. The full-dimensional restriction is intentional: an embedded k < N simplex needs an additional normal-frame convention before it has a scalar orientation. This equality is also not a no-tunnelling barrier; a sufficiently large discrete update may cross or land on the zero-measure set.

XpbdParticleN adds velocity, force, gravity scale, and a stable world-local id to that point coordinate. XpbdWorldN.step() performs semi-implicit prediction, XPBD projection, velocity reconstruction, then ordered XpbdVelocityResponseN policies for every substep. Responses may change only declared registered velocities and retain their evidence beside the matching solve result. Forces are held across the outer step and clear only on success. Read-only XpbdStateGuardN policies then accept or reject the completed substep. Late evaluator, response, or guard errors restore position, velocity, force, and gravity scale transactionally. Fixed particles remain outside prediction and do not acquire an inferred kinematic trajectory.

compileXpbdParticleBindingN() owns the topology-neutral one-particle-per- source-vertex correspondence and transactional source write-back. It keeps positive physical mass separate from the fixed mobility policy, so pinning a vertex does not erase its mass evidence. lumpSimplexMassesN() supplies an auditable diagonal reference mass by integrating density against intrinsic simplex rest measure and equally accumulating each element mass onto its incident vertices. It reports element and vertex totals independently.

XpbdParticleHyperplaneConstraintN declares the normalized point gap to an oriented RN hyperplane as a non-negative scalar relation. compileXpbdParticleHyperplaneFamilyN() composes one such constraint per source vertex over an existing particle binding, retaining source ordinal, compile-time gap, clearance, compliance, and exact particle identity. It is a discrete point-contact reference. The optional compileXpbdParticleHyperplaneFrictionFamilyN() consumes the same normal solves after velocity reconstruction and projects the desired stopping impulse onto the complete RN Coulomb tangent ball. In R4 that is an isotropic three-ball, not three scalar clamps. XpbdExponentialVelocityDampingN provides separate timestep-invariant decay with an inverse-seconds rate. These are not deformable surface contact, restitution, or continuous collision.

For one standalone point, construct XpbdParticleHyperplaneConstraintN directly. The compile*FamilyN form intentionally requires a real CellComplex and one bound particle per source vertex because its additional purpose is to preserve that source correspondence; it is not a more general single-particle constructor.

compileXpbdDistanceNetworkN() turns one explicitly selected two-vertex CellComplex 1-cell group into distance constraints. It can retain its compatible self-contained particle-authoring path or compose over an existing source-indexed particle binding. In the composed path it preserves exact particle identities and takes rest lengths from source geometry, so compiling constraints after deformation does not silently redefine rest. Every edge retains structural source identity. Source positions do not alias the simulation; an explicit binding or standalone-network write synchronizes them only after complete validation.

compileXpbdSimplexMeasureFamilyN() compiles one explicitly selected simplex cell group onto an existing source-indexed XpbdParticleN array. It retains a structural source id and vertex tuple per cell, derives default rest measure from source geometry rather than possibly deformed live particles, and keeps rest/compliance policies separate from topology. The family owns no particles and performs no write-back. addToWorld() requires the exact particle objects to be registered already, then preflights every lineage and constraint id before attaching the family atomically. This lets distance and local measure coordinates share one RN state without implying a complete deformable-body model.

compileXpbdOrientedCuboidFamilyN() accepts an explicitly selected full-dimensional cuboid group and applies the core's deterministic Kuhn simplexization internally. Each generated signed-measure constraint retains the structural id of its authored parent cuboid, the parent-cell ordinal, the axis-permutation ordinal and tuple, and both source vertex tuples. The raw simplex signs alternate with permutation parity; the compiler preserves that auditable ordering rather than silently rewinding cells. Rest coordinates come from source geometry, material callbacks remain separate, and the family shares an existing source-indexed particle array. World attachment preflights all parent lineage, particle ownership, and constraint ids before adding any constraint.

evaluateSimplexMetricDeformationN() compares matching rest and current k-simplices in RN through their intrinsic edge Gram metrics. Cholesky normalization expresses the current metric in an orthonormal rest-material basis, yielding the right Cauchy–Green tensor, Green–Lagrange strain, ordered principal stretches, measure ratio, rest-conditioning evidence, and spectral residual. It applies equally to embedded curves/membranes and full-dimensional solids without an ambient cross product. Only the full-dimensional case reports a signed measure ratio and preserved/inverted/collapsed state; embedded simplices require an authored normal frame before scalar orientation is meaningful. The coordinate selects no constitutive energy and produces no forces by itself.

SimplexConstitutiveEvaluationN is the shared rest-measure, energy, second Piola stress, and analytic current-gradient contract. The package supplies two Float64 laws over it: evaluateSimplexStVenantKirchhoffN() for the polynomial small-strain reference and evaluateSimplexCompressibleNeoHookeanN() for a large-strain logarithmic-volume reference. Neo-Hookean embedded elements use positive intrinsic measure; full-dimensional elements must preserve signed orientation. Collapse and inversion refuse explicitly. This evaluator is not an inversion barrier or an implicit solver.

SimplexConstitutiveLawN and compileSimplexConstitutiveFamilyN() assemble a typed law over one explicit source simplex group while retaining copied rest state, structural cell ids, live lineage, exact particle identities, and deterministic shared-vertex forces. Immutable StVK and Neo-Hookean descriptors are built in. Their named family compilers are typed convenience wrappers over the same implementation and preserve the existing StVK API/provider identity. The three shipped laws—StVK, compressible Neo-Hookean, and the smooth lower-measure barrier—also provide exact matrix-free potential Hessian-vector products. The law-level evaluations retain directional right-Cauchy–Green and second-Piola tensors; family products assemble by source vertex and plug into the complete incremental-objective analytic curvature protocol. Custom laws may remain first-order-only and are then refused explicitly by that protocol. No dense Hessian, definiteness modification, or Newton/Krylov solver is implied by the provider capability itself.

Analytic composition and the Newton APIs default to the providers' exact curvature. Authors may instead select curvaturePolicy: { kind: 'provider-local-psd' }. That explicit modified-Newton reference reconstructs each complete provider Hessian from basis HVPs, audits symmetry, diagonalizes it with the deterministic Float64 eigensolver, and clamps negative eigenvalues to zero. Results retain raw and projected spectra, clipped counts, symmetry error, eigensystem residuals, and operator cost.

Providers may expose a finer exact additive decomposition through XpbdConservativeHessianBlockProviderN. Selecting curvaturePolicy: { kind: 'provider-block-psd' } reconstructs and projects each declared block independently, then audits the raw block sum against the provider's authoritative aggregate HVP. Providers without that capability remain valid and visibly use one implicit-provider block. Constitutive families declare one source-ordered block per simplex, retaining element lineage in SimplexConstitutiveFamilyHessianBlockN.

Both modes are deterministic cubic-cost CPU golden paths. Provider-local cost is cubic in the whole provider variable count; block-local cost is the sum of the dense block costs. The latter supplies an auditable element-local reference for simplex materials, not a sparse matrix, production preconditioner, or large-mesh factorization.

compileXpbdIncrementalPotentialAnalyticHessianOperatorN() fixes one candidate coordinate and separates curvature construction from application. Exact curvature remains matrix-free. Provider-local and provider-block basis HVPs, symmetry audits, and eigendecompositions are paid once; subsequent products reuse the stored projected matrices. Block-local products still request one exact aggregate provider HVP per direction so the authored block sum remains audited rather than assumed. The compilation evidence states both the one-time and per-product provider costs.

solveXpbdIncrementalPotentialNewtonDirectionN() composes the complete analytic objective product into a bounded, non-mutating preconditioned-CG reference for H(q) p = -gradient(Phi(q)). Identity and exact inertial mass-diagonal preconditioners are available. Results retain per-iteration residual and curvature evidence, plus actual construction/application provider HVP counts, and distinguish convergence, an exact zero gradient, budget exhaustion, unsupported providers, and non-positive or numerically unresolved curvature. The solver compiles projected curvature once at each linearization coordinate and reuses it throughout CG. In exact mode it assembles no matrix and does not modify definiteness. Provider-local and provider-block PSD are the explicit exceptions described above. No mode chooses a nonlinear step, runs Armijo, or applies state.

compileSimplexConstitutiveFamilyStateGuardN() is an optional post-substep policy over that generic family. It rejects typed law-domain refusal, full-dimensional orientation change, or a configured positive minimum measure ratio. XpbdWorldN.stepAdaptive() rolls those typed rejections back and retries the same outer duration with bounded deterministic subdivision. It does not retry arbitrary failures, repair an invalid material, or prove continuous orientation preservation between accepted endpoints.

relativeOrientationCoordinates4() provides the analogous local coordinate for rotation. It chooses one lift of the paired-quaternion double cover, returns a reusable branch token for coherent timesteps, and reports the full SO(4) logarithm cut locus as a discriminated result rather than manufacturing an axis. orientationDexp4() and orientationDlog4() expose the matching 6x6 Jacobians in either world-left or body-right trivialization. The right factor uses the opposite Jacobian sign because Rotor4 composes that quaternion in reverse order. These are proof-kernel primitives; no hinge, cone, limit, or motor policy is implied yet.

ConstraintBlockSolver4 couples one to six rows through their complete J M^-1 J^T response. Equality blocks retain the original default. An explicit one-bounded projection may add exactly one force-limited coordinate: the solver eliminates the remaining equalities through a scalar Schur complement, clamps that coordinate, and re-solves the equality subspace exactly. Diagnostics distinguish raw speed error from the projected KKT residual. The default rank policy refuses lost coordinates; an explicit minimum-norm policy remains available only for unbounded equality diagnostics. Bias limiting and warm-start transport preserve orthogonal equality-basis invariance. PointJointSolver4 is a compatibility wrapper over this shared kernel.

DirectionJoint4 binds one body-local unit direction to another local or fixed-world direction. constraint() returns either a regular three-row block or a typed antipodal refusal. The three rows constrain the tangent space of the direction sphere and leave the non-abelian SO(3) stabilizer free, so the joint deliberately exposes no fictitious scalar “hinge angle.”

OrientationJoint4 binds a complete body-local frame to another material frame or a fixed world frame. Its six equality rows are the full rotational analogue of a weld; combine them with the four rows of PointJoint4 when both orientation and translation must be fixed. The error is log(inverse(frameB) * frameA), expressed in frame B, and the analytic world-left rate rows are exact negatives for the two participants. This makes the coordinate invariant under common world rotation and keeps internal angular impulses equal and opposite. The paired-quaternion lift is retained across evaluations, while the non-unique SO(4) cut locus returns a typed result with no solver block.

PlanarRotationJoint4 binds an ordered body-local orthonormal two-frame to a local or fixed-world frame. Its five-row Stiefel constraint fixes that frame and leaves only SO(2) rotation in the orthogonal plane. First-axis antipodes and degenerate second bisectors are typed refusals. This is deliberately distinct from oriented-plane preservation, which would leave a two-angle torus free.

PlanarRotationCoordinate4 attaches one phase-reference direction to each side of that joint. It reports a signed wrapped angle, a persistent unwrapped angle, the positively oriented complementary-plane bivector, and its angular speed. A sample exactly half a turn from the preceding branch is a typed unwrap-ambiguous result until the caller chooses its sign; samples must be frequent enough that an unobserved advance never reaches pi.

PlanarRotationMotor4 adds the oriented phase row to the five frame rows and tracks signed angular speed under a symmetric maxTorque bound. PlanarRotationIntervalJoint4 returns two persistent guardian blocks over the continuous unwrapped angle. Their first-order speed corridor catches unsafe motion introduced by earlier motor or constraint blocks in the same projected iteration. Singular frame charts and ambiguous half-turn lifts remain typed results rather than implicit branch choices.

For automatic mixed contact, register GlomeCollider4, PolytopeCollider4, HyperplaneContactCollider4, and/or HyperboxCollider4 instances with ContactPipeline4, then call pipeline.stepWorld(world, fixedDt). Finite colliders share conservative AABBs and temporally coherent sweep-and-prune; infinite planes are paired explicitly with every admitted compact collider. HyperboxContactPipeline4 remains the narrower homogeneous box path. The exhaustive O(n²) finite provider remains available as the CPU golden reference. For fast rigidly moving bodies, pipeline.stepWorldContinuous() is an opt-in event loop which advances to certified first impact before using the same manifold solver. Compact pairs are pruned by conservative swept AABBs, and each event scan retains broadphase diagnostics; the exhaustive provider remains the differential reference. Its result explicitly reports prescribed-motion and cast-uncertainty fallbacks; the discrete default is unchanged.

RigidTrajectory4 makes an R4 screw path an explicit reusable value rather than an assumption hidden inside a solver. convexRigidCast4() and supportShapeHyperplaneRigidCast4() conservatively advance along that exact path. Their closing-speed certificate adds each body's tight angularVelocityOperatorNorm4(generator) * boundingRadius contribution to the linear normal closure. Built-in glomes, rounded shapes, transformed shapes, and vertex-enumerable polytopes have auditable inferred radii; opaque support functions must provide a validated explicit bound. RigidBodyPosePlan4 freezes the same momentum-derived Lie-midpoint generator used by ordinary free flight. stepWorldContinuous() gives each event scan those plans and applies the exact same plans to the selected impact; response then changes momentum and causes the remainder to be replanned. No-impact rotational advancement is therefore endpoint-identical to PhysicsWorld4.integratePoses().

rigidTrajectoryFromTransforms4() constructs the principal screw segment between two coherent R4 poses. KinematicBody4 attaches a physical duration to that segment, owns its current position and rotation, and exposes the exact linear and world-left angular rates used by contact response. It accepts no impulses. Discrete world seams and stepWorldContinuous() advance registered kinematic compact colliders through the same absolute subplans used by swept broadphase and casting. Centered glomes preserve the analytic linear fast path; rotating hyperboxes, polytopes, and offset glomes use rigid casts. A legacy velocity-only RigidMotion4 still produces a typed partial fallback because no geometry path can be inferred honestly from velocity alone.

KinematicTrackDriver4 produces those segments from one position sampler and one Rotor4Track on a fixed clock. It samples each accepted boundary once, caches the shared endpoint between consecutive segments, and refuses to replace a segment before the body reaches it. CCD therefore consumes a frozen physical trajectory even when an event step is subdivided; animation is never resampled inside the collision loop. The adapter has no renderer or mixer dependency.

gjkDistance separates a stable numerical estimate from a certified result: separated is reported only with a support-gap certificate and intersecting only with an origin-enclosure proof, while iteration-limit is an explicit refusal whose accompanying distance is an estimate, never a claim. Equal and nearly tied support directions terminate — a repeated support point triggers a certificate-aware reprojection of the complete sampled support set — and a proved fixpoint refuses immediately as duplicate-support rather than burning the remaining budget on an identical cycle.

NarrowphaseDispatcherN is the common query boundary. Its best mode selects the strongest honest capability for the configured pair and margins; explicit requests never silently fall back. Stable ordered pair IDs provide coherent GJK warm starts and deterministic batch retirement. Zero-margin compact R4 pairs can return a bounded EPA minimum-translation witness. When both shapes also enumerate stable source vertices, polytope4 derives their facet halfspaces and clips a complete response-grade contact manifold whose face-pair IDs persist under coherent rigid motion. Its dimension-independent topology compiler turns the exhaustive facet search into reusable source-ID incidence; live queries reconstruct and validate only the current facet planes. PolytopeCollider4 caches this product by source identity by default. The same incidence product supplies complete point-through-polyhedron support-face contact against an infinite plane, with stable source-vertex IDs and affine-span-preserving solver-point reduction. Deep results retain an algorithm discriminator, so a smooth point patch is never mistaken for a polyhedral patch.

contactConstraintFromSmoothPointPatch4() connects either analytic smooth family to the existing coupled R4 friction solver. Coincident glome centers remain observable but non-responding because their minimum-translation normal is not unique. Mixed R4 adapters preserve either the single glome/box witness or the complete point-through-polyhedron box/plane support feature; an interior glome/box tie likewise stays observable without manufacturing a direction.

MIT © Nikolay Petrov