From 06ec4cbea25430c7abe754b0d69cdc804d831def Mon Sep 17 00:00:00 2001 From: m3ultra Date: Thu, 16 Jul 2026 21:28:59 +1000 Subject: [PATCH] Add 3D sail cloth sim with per-face wind and XPBD corner loads Verlet cloth on a bilinear patch between 4 anchors, N=10 grid, structural/shear/bend constraints, 5 relaxation iterations at a fixed 1/60 substep. Wind is applied per FACE so hypar twist genuinely sheds load rather than being cosmetic. Two deviations from PLAN3D worth flagging: - Load is read from each constraint's XPBD Lagrange multiplier, not from FABRIC_K * leftover-stretch. After a fixed iteration count the leftover stretch is solver error, not fabric strain, so the naive reading came out ~50x hot (60 kN peaks on a 5x5 m sail). The multiplier is the real constraint impulse, which the statics assert confirms by balancing the corner reactions against the applied wind to 8%. - Wind uses a signed square (d*|d|) rather than clamp(d)^2, so the leeward face is pushed too. A sail is double-sided. The sim core deliberately does not import three.js: it runs headless under node today, stays allocation-free in the hot loop, and replays bit-for-bit. createSailView() pulls three in lazily for rendering. Loads land in real newtons (~1-4 kN on a 5x5 m sail in a 34 m/s storm), so hardware ratings are real working load limits. Co-Authored-By: Claude Opus 4.8 --- web/world/js/sail.js | 617 ++++++++++++++++++++++++++++++++++ web/world/js/sail.selftest.js | 352 +++++++++++++++++++ 2 files changed, 969 insertions(+) create mode 100644 web/world/js/sail.js create mode 100644 web/world/js/sail.selftest.js diff --git a/web/world/js/sail.js b/web/world/js/sail.js new file mode 100644 index 0000000..49cd8f9 --- /dev/null +++ b/web/world/js/sail.js @@ -0,0 +1,617 @@ +/** + * sail.js — shade sail cloth simulation, corner loads, hardware failure. [Lane B] + * + * A 3D verlet cloth on a bilinear patch between 4 anchors. Wind pressure is + * applied per FACE, not per node, which is the whole point: a twisted (hypar) + * sail turns most of its faces edge-on to the wind and sheds load, while a flat + * one presents every face square-on and catches everything. That difference is + * the game's thesis and it is asserted in sail.selftest.js. + * + * Units are SI throughout: metres, kilograms, seconds, newtons. Corner loads + * come out in real newtons and hardware ratings are real working load limits, + * so a 5x5 m sail in a 30 m/s gust genuinely puts ~5 kN on a corner — which is + * genuinely why real shade sails use 3 kN+ shackles. + * + * This module deliberately does NOT import three.js. The sim core is plain + * typed arrays so it runs headless under node (see sail.selftest.js) before any + * renderer exists, stays allocation-free in the hot loop, and can be replayed + * bit-for-bit. three.js is pulled in lazily by createSailView() only. + */ + +// ---------- sim tunables ---------- +const SIM_DT = 1 / 60; // sim always steps at a fixed rate; step() accumulates +const MAX_SUBSTEPS = 5; // spiral-of-death guard when the frame hitches +const RELAX_ITERS = 5; // FABRIC_K is calibrated against this; changing it rescales loads +const GRAVITY = -9.81; + +// ---------- aerodynamics ---------- +// 0.5 * air density (1.225) * flat-plate drag coefficient (~1.4). +// Newtons per m^2 of face area per (m/s)^2 of normal-on airflow. +const PRESSURE_COEFF = 0.86; +const TANGENT_COEFF = 0.02; // skin friction dragging along the face +const MAX_NORMAL_SPEED = 45; // clamp on the normal-on component, m/s — stability in extreme gusts + +// ---------- fabric ---------- +const FABRIC_DENSITY = 0.32; // kg/m^2, typical knitted shade cloth +// Axial stiffness of one grid spring, N/m — roughly E*t*width/length for +// knitted HDPE mesh. Fed to the solver as a compliance (1/k), not used to +// convert stretch into force: see _measureLoads for why that distinction is +// the whole ballgame. +const FABRIC_K = 100000; +const K_COMPRESS = 0.08; // cloth resists stretch hard, compression barely (from prototype) +const K_BEND = 0.04; +const COMP_STRETCH = 1 / FABRIC_K; +const COMP_COMPRESS = 1 / (FABRIC_K * K_COMPRESS); +const COMP_BEND = 1 / (FABRIC_K * K_BEND); +const VEL_DAMP = 0.995; // light; relative-wind drag supplies the real damping + +// ---------- failure ---------- +const OVERLOAD_SECS = 0.4; // prototype: 0.4 s sustained overload before it lets go +const OVERLOAD_RECOVER = 2.0; // prototype: overload timer bleeds off at 2x +const LOAD_TAU = 0.11; // load meter smoothing time constant, s + +/** + * Hardware tiers. Costs are the prototype's economy verbatim; ratings are + * retuned from the prototype's arbitrary 9/19/40 into real newtons, preserving + * the same relative spread. `rating` is a working load limit in N. + */ +export const HARDWARE = [ + { name: 'carabiner', cost: 5, rating: 1200, color: '#e2b04a' }, + { name: 'shackle', cost: 15, rating: 3200, color: '#c8d2d8' }, + { name: 'rated shackle', cost: 30, rating: 6500, color: '#7ee0ff' }, +]; + +export const TENSION_MIN = 0.6; +export const TENSION_MAX = 1.4; +const TRIM_MIN = 0.85; +const TRIM_MAX = 1.15; + +const clamp = (v, lo, hi) => (v < lo ? lo : v > hi ? hi : v); + +/** + * Order 4 anchors into a non-self-intersecting ring by angle around their + * centroid, projected onto the ground plane. Ported from the prototype's + * orderRing; without it, picking corners in a silly order knots the sail. + */ +export function orderRing(anchors) { + const n = anchors.length; + let cx = 0, cz = 0; + for (const a of anchors) { cx += a.pos.x; cz += a.pos.z; } + cx /= n; cz /= n; + return [...anchors].sort( + (a, b) => Math.atan2(a.pos.z - cz, a.pos.x - cx) - Math.atan2(b.pos.z - cz, b.pos.x - cx) + ); +} + +export class SailRig { + /** + * @param {object} opts + * @param {Array} opts.anchors world.anchors — [{id, pos:{x,y,z}, type, sway(t)->{x,y,z}}] + * sway(t) returns an OFFSET to add to pos, not an absolute position. + * @param {number} opts.gridN nodes per side (default 10) + * @param {number} opts.porosity 0 = solid membrane, ~0.3 = knitted shade cloth (blows through, less load) + */ + constructor({ anchors = [], gridN = 10, porosity = 0 } = {}) { + this.anchors = anchors; + this.N = gridN; + this.porosity = porosity; + this.corners = []; + this.events = []; + this.tension = 1.0; + this.t = 0; + this.rigged = false; + this._acc = 0; + // scratch, reused every face to keep the hot loop allocation-free + this._probe = { x: 0, y: 0, z: 0 }; + } + + /** + * Rig the sail across 4 anchors. + * @param {string[]} anchorIds 4 anchor ids; reordered into a ring internally + * @param {object[]} hwChoices hardware per anchor id, same order as anchorIds + * @param {number} tension 0.6 (loose, flogs) .. 1.4 (drum tight, shock-loads) + */ + attach(anchorIds, hwChoices, tension = 1.0) { + if (anchorIds.length !== 4) throw new Error(`sail needs exactly 4 corners, got ${anchorIds.length}`); + + const picked = anchorIds.map((id) => { + const a = this.anchors.find((x) => x.id === id); + if (!a) throw new Error(`unknown anchor "${id}"`); + return a; + }); + const hwById = new Map(anchorIds.map((id, i) => [id, hwChoices[i] || HARDWARE[0]])); + + const ring = orderRing(picked); + this.tension = clamp(tension, TENSION_MIN, TENSION_MAX); + this.corners = ring.map((a) => ({ + anchorId: a.id, + anchor: a, + hw: hwById.get(a.id), + load: 0, + peakLoad: 0, + overload: 0, + broken: false, + trim: 1.0, + loadVec: { x: 0, y: 0, z: 0 }, // reaction direction, not just magnitude — see _measureLoads + })); + + this._build(ring); + this.rigged = true; + return this; + } + + _build(ring) { + const N = this.N; + const nodeCount = N * N; + this.pos = new Float64Array(nodeCount * 3); + this.prev = new Float64Array(nodeCount * 3); + this.force = new Float64Array(nodeCount * 3); + this.invMass = new Float64Array(nodeCount); + + // Bilinear patch across the 4 corners. Because the anchors sit at different + // heights, this initial surface is already a hypar — the sim just relaxes it. + const [c0, c1, c2, c3] = ring.map((a) => a.pos); + for (let v = 0; v < N; v++) { + for (let u = 0; u < N; u++) { + const fu = u / (N - 1), fv = v / (N - 1); + const i = (v * N + u) * 3; + for (let k = 0; k < 3; k++) { + const ax = ['x', 'y', 'z'][k]; + const top = (1 - fu) * c0[ax] + fu * c1[ax]; + const bot = (1 - fu) * c3[ax] + fu * c2[ax]; + this.pos[i + k] = this.prev[i + k] = (1 - fv) * top + fv * bot; + } + } + } + + const idx = (u, v) => v * N + u; + this.cornerIdx = [idx(0, 0), idx(N - 1, 0), idx(N - 1, N - 1), idx(0, N - 1)]; + + // springs: structural + shear carry load; bend only resists folding + this.springs = []; + const link = (a, b, kind) => { + const ax = a * 3, bx = b * 3; + const dx = this.pos[bx] - this.pos[ax]; + const dy = this.pos[bx + 1] - this.pos[ax + 1]; + const dz = this.pos[bx + 2] - this.pos[ax + 2]; + this.springs.push({ a, b, restBase: Math.hypot(dx, dy, dz), rest: 0, kind }); + }; + for (let v = 0; v < N; v++) { + for (let u = 0; u < N; u++) { + if (u < N - 1) link(idx(u, v), idx(u + 1, v), 'struct'); + if (v < N - 1) link(idx(u, v), idx(u, v + 1), 'struct'); + if (u < N - 1 && v < N - 1) { + link(idx(u, v), idx(u + 1, v + 1), 'shear'); + link(idx(u + 1, v), idx(u, v + 1), 'shear'); + } + if (u < N - 2) link(idx(u, v), idx(u + 2, v), 'bend'); + if (v < N - 2) link(idx(u, v), idx(u, v + 2), 'bend'); + } + } + + // XPBD Lagrange multipliers, one per spring, reset every substep + this.lambda = new Float64Array(this.springs.length); + + // Springs meeting each corner, kept as {spring, index} so the load meter can + // look up each one's multiplier. Bend springs are included: the hardware + // physically carries every element that touches it, and leaving them out + // under-reports the reaction and breaks the statics balance. + this.cornerSprings = this.cornerIdx.map((ci) => + this.springs + .map((s, si) => ({ s, si })) + .filter(({ s }) => s.a === ci || s.b === ci) + ); + + // triangles: wind acts per face, and coverage raycasts against these + this.tris = new Uint16Array((N - 1) * (N - 1) * 6); + let ti = 0; + for (let v = 0; v < N - 1; v++) { + for (let u = 0; u < N - 1; u++) { + const a = idx(u, v), b = idx(u + 1, v), c = idx(u + 1, v + 1), d = idx(u, v + 1); + this.tris[ti++] = a; this.tris[ti++] = b; this.tris[ti++] = c; + this.tris[ti++] = a; this.tris[ti++] = c; this.tris[ti++] = d; + } + } + + // grid-space proximity of every node to each corner, for per-corner trim + this._cornerWeight = []; + for (let k = 0; k < 4; k++) { + const cu = [0, N - 1, N - 1, 0][k], cv = [0, 0, N - 1, N - 1][k]; + const w = new Float64Array(nodeCount); + for (let v = 0; v < N; v++) { + for (let u = 0; u < N; u++) { + const dist = Math.hypot(u - cu, v - cv) / (N - 1); + w[idx(u, v)] = Math.max(0, 1 - dist); + } + } + this._cornerWeight.push(w); + } + + this.area = this._surfaceArea(); + const mass = (FABRIC_DENSITY * this.area) / nodeCount; + this.nodeMass = mass; + this.invMass.fill(1 / mass); + + this._applyRestLengths(); + this._repin(0); + } + + /** Rest lengths shrink as tension rises (1/tension, ported), modulated per corner by trim. */ + _applyRestLengths() { + for (const s of this.springs) { + let wsum = 0, tsum = 0; + for (let k = 0; k < 4; k++) { + const w = this._cornerWeight[k][s.a] + this._cornerWeight[k][s.b]; + wsum += w; + tsum += w * this.corners[k].trim; + } + const trim = wsum > 1e-9 ? tsum / wsum : 1; + s.rest = s.restBase / (this.tension * trim); + } + } + + /** Pinned corners are infinite-mass so springs stretch honestly against them. */ + _repin(t) { + for (let k = 0; k < 4; k++) { + const c = this.corners[k]; + const ci = this.cornerIdx[k]; + if (c.broken) { + this.invMass[ci] = 1 / this.nodeMass; // freed node — flogging falls out of this + continue; + } + this.invMass[ci] = 0; + const p = this._anchorPos(c.anchor, t); + this.pos[ci * 3] = p.x; this.pos[ci * 3 + 1] = p.y; this.pos[ci * 3 + 2] = p.z; + this.prev[ci * 3] = p.x; this.prev[ci * 3 + 1] = p.y; this.prev[ci * 3 + 2] = p.z; + } + } + + _anchorPos(a, t) { + const p = a.pos; + if (!a.sway) return p; + const s = a.sway(t); + return { x: p.x + s.x, y: p.y + s.y, z: p.z + s.z }; + } + + _surfaceArea() { + let total = 0; + for (let i = 0; i < this.tris.length; i += 3) { + const a = this.tris[i] * 3, b = this.tris[i + 1] * 3, c = this.tris[i + 2] * 3; + const e1x = this.pos[b] - this.pos[a], e1y = this.pos[b + 1] - this.pos[a + 1], e1z = this.pos[b + 2] - this.pos[a + 2]; + const e2x = this.pos[c] - this.pos[a], e2y = this.pos[c + 1] - this.pos[a + 1], e2z = this.pos[c + 2] - this.pos[a + 2]; + const nx = e1y * e2z - e1z * e2y, ny = e1z * e2x - e1x * e2z, nz = e1x * e2y - e1y * e2x; + total += Math.hypot(nx, ny, nz) * 0.5; + } + return total; + } + + /** + * Advance the sim. Accumulates real time and burns it in fixed SIM_DT chunks, + * so a variable-rate render loop and a fast-forwarded selftest produce + * identical traces. Never reads a clock. + * + * @param {number} dt seconds elapsed since last call + * @param {object} wind { sample(pos, t) -> {x,y,z} } + * @param {number} t world time, seconds + */ + step(dt, wind, t) { + if (!this.rigged) return; + this._acc += dt; + let n = 0; + while (this._acc >= SIM_DT && n < MAX_SUBSTEPS) { + this._substep(SIM_DT, wind, this.t); + this._acc -= SIM_DT; + this.t += SIM_DT; + n++; + } + if (n === MAX_SUBSTEPS) this._acc = 0; // dropped frames: don't try to catch up + } + + _substep(dt, wind, t) { + this._accumulateWind(wind, t, dt); + this._integrate(dt); + this.lambda.fill(0); // XPBD multipliers are per-substep + for (let i = 0; i < RELAX_ITERS; i++) this._relax(dt * dt); + this._pinCorners(t); + this._measureLoads(dt); + this._checkFailure(dt); + } + + /** Wind force per FACE — the hypar mechanic lives here. */ + _accumulateWind(wind, t, dt) { + const pos = this.pos, prev = this.prev, F = this.force; + F.fill(0); + const coeff = PRESSURE_COEFF * (1 - this.porosity); + const tanCoeff = TANGENT_COEFF * (1 - this.porosity); + const invDt = 1 / dt; + const probe = this._probe; + + for (let i = 0; i < this.tris.length; i += 3) { + const ia = this.tris[i] * 3, ib = this.tris[i + 1] * 3, ic = this.tris[i + 2] * 3; + + const e1x = pos[ib] - pos[ia], e1y = pos[ib + 1] - pos[ia + 1], e1z = pos[ib + 2] - pos[ia + 2]; + const e2x = pos[ic] - pos[ia], e2y = pos[ic + 1] - pos[ia + 1], e2z = pos[ic + 2] - pos[ia + 2]; + // |cross| is twice the area and its direction is the face normal + let nx = e1y * e2z - e1z * e2y, ny = e1z * e2x - e1x * e2z, nz = e1x * e2y - e1y * e2x; + const len = Math.hypot(nx, ny, nz); + if (len < 1e-9) continue; + const area = len * 0.5; + nx /= len; ny /= len; nz /= len; + + probe.x = (pos[ia] + pos[ib] + pos[ic]) / 3; + probe.y = (pos[ia + 1] + pos[ib + 1] + pos[ic + 1]) / 3; + probe.z = (pos[ia + 2] + pos[ib + 2] + pos[ic + 2]) / 3; + const w = wind.sample(probe, t); + + // Relative wind, not absolute: as the cloth accelerates downwind the load + // bleeds off by itself. This is what stops flogging from exploding. + const vx = (pos[ia] - prev[ia] + pos[ib] - prev[ib] + pos[ic] - prev[ic]) / 3 * invDt; + const vy = (pos[ia + 1] - prev[ia + 1] + pos[ib + 1] - prev[ib + 1] + pos[ic + 1] - prev[ic + 1]) / 3 * invDt; + const vz = (pos[ia + 2] - prev[ia + 2] + pos[ib + 2] - prev[ib + 2] + pos[ic + 2] - prev[ic + 2]) / 3 * invDt; + const rx = w.x - vx, ry = w.y - vy, rz = w.z - vz; + + const d = clamp(rx * nx + ry * ny + rz * nz, -MAX_NORMAL_SPEED, MAX_NORMAL_SPEED); + // d*|d| rather than d^2: keeps the v^2 magnitude but points the force the + // way the wind is actually blowing. A sail is double-sided. + const p = coeff * area * d * Math.abs(d); + + const tx = (rx - nx * d) * tanCoeff * area; + const ty = (ry - ny * d) * tanCoeff * area; + const tz = (rz - nz * d) * tanCoeff * area; + + const fx = (nx * p + tx) / 3, fy = (ny * p + ty) / 3, fz = (nz * p + tz) / 3; + F[ia] += fx; F[ia + 1] += fy; F[ia + 2] += fz; + F[ib] += fx; F[ib + 1] += fy; F[ib + 2] += fz; + F[ic] += fx; F[ic + 1] += fy; F[ic + 2] += fz; + } + } + + _integrate(dt) { + const pos = this.pos, prev = this.prev, F = this.force, im = this.invMass; + const dt2 = dt * dt; + for (let n = 0; n < im.length; n++) { + if (im[n] === 0) continue; // pinned + const i = n * 3; + for (let k = 0; k < 3; k++) { + const a = F[i + k] * im[n] + (k === 1 ? GRAVITY : 0); + const x = pos[i + k]; + const nx = x + (x - prev[i + k]) * VEL_DAMP + a * dt2; + prev[i + k] = x; + pos[i + k] = nx; + } + } + } + + /** + * XPBD constraint solve. The plain-PBD version of this is simpler, but its + * position corrections carry no force information — the leftover stretch + * after a fixed iteration count is solver error, not fabric strain, so + * reading load off it measures the solver. XPBD's Lagrange multiplier lambda + * is the real constraint impulse, so lambda/dt^2 is a genuine newton value + * and the corner reactions balance the applied wind by construction. + */ + _relax(dt2) { + const pos = this.pos, im = this.invMass, lam = this.lambda; + for (let si = 0; si < this.springs.length; si++) { + const s = this.springs[si]; + const wa = im[s.a], wb = im[s.b]; + const w = wa + wb; + if (w === 0) continue; // both ends pinned + const ia = s.a * 3, ib = s.b * 3; + const dx = pos[ib] - pos[ia], dy = pos[ib + 1] - pos[ia + 1], dz = pos[ib + 2] - pos[ia + 2]; + const d = Math.hypot(dx, dy, dz); + if (d < 1e-9) continue; + const C = d - s.rest; + const compliance = s.kind === 'bend' ? COMP_BEND : C > 0 ? COMP_STRETCH : COMP_COMPRESS; + const at = compliance / dt2; + const dLambda = (-C - at * lam[si]) / (w + at); + lam[si] += dLambda; + // grad C is -n for node a and +n for node b, with n = (b - a)/d + const nx = dx / d, ny = dy / d, nz = dz / d; + pos[ia] -= nx * dLambda * wa; pos[ia + 1] -= ny * dLambda * wa; pos[ia + 2] -= nz * dLambda * wa; + pos[ib] += nx * dLambda * wb; pos[ib + 1] += ny * dLambda * wb; pos[ib + 2] += nz * dLambda * wb; + } + } + + _pinCorners(t) { + for (let k = 0; k < 4; k++) { + const c = this.corners[k]; + if (c.broken) continue; + const ci = this.cornerIdx[k] * 3; + const p = this._anchorPos(c.anchor, t); + this.pos[ci] = p.x; this.pos[ci + 1] = p.y; this.pos[ci + 2] = p.z; + } + } + + /** + * Corner load = magnitude of the VECTOR sum of the tensions in the springs + * meeting that corner, each read off its XPBD multiplier as |lambda|/dt^2. + * + * Reading tension as FABRIC_K * leftover-stretch instead looks equivalent and + * is not: after a fixed 5 iterations the leftover stretch is solver error, so + * that number measures the solver rather than the fabric and comes out ~50x + * hot. The multiplier is the actual constraint impulse, which is why the + * statics assert balances. + * + * The vector sum (rather than a scalar total) is what + * makes DESIGN.md's anchor-angle mechanic fall out for free: edges pulling in + * nearly the same direction add up, edges pulling apart partly cancel — so a + * pinched corner really does multiply its own load. + */ + _measureLoads(dt) { + const pos = this.pos, lam = this.lambda; + const invDt2 = 1 / (dt * dt); + const alpha = 1 - Math.exp(-dt / LOAD_TAU); + for (let k = 0; k < 4; k++) { + const c = this.corners[k]; + if (c.broken) { c.load = 0; c.loadVec.x = c.loadVec.y = c.loadVec.z = 0; continue; } + const ci = this.cornerIdx[k], cix = ci * 3; + let sx = 0, sy = 0, sz = 0; + for (const { s, si } of this.cornerSprings[k]) { + if (lam[si] >= 0) continue; // slack or compressed fabric pulls on nothing + const tension = -lam[si] * invDt2; // the multiplier IS the impulse; /dt^2 makes it newtons + const o = (s.a === ci ? s.b : s.a) * 3; + const dx = pos[o] - pos[cix], dy = pos[o + 1] - pos[cix + 1], dz = pos[o + 2] - pos[cix + 2]; + const d = Math.hypot(dx, dy, dz); + if (d < 1e-9) continue; + sx += (dx / d) * tension; sy += (dy / d) * tension; sz += (dz / d) * tension; + } + c.loadVec.x += (sx - c.loadVec.x) * alpha; + c.loadVec.y += (sy - c.loadVec.y) * alpha; + c.loadVec.z += (sz - c.loadVec.z) * alpha; + const raw = Math.hypot(sx, sy, sz); + c.load += (raw - c.load) * alpha; + if (c.load > c.peakLoad) c.peakLoad = c.load; + } + } + + /** Ported from the prototype: 0.4 s sustained over the rating and it lets go. */ + _checkFailure(dt) { + for (const c of this.corners) { + if (c.broken) continue; + if (c.load > c.hw.rating) c.overload += dt; + else c.overload = Math.max(0, c.overload - dt * OVERLOAD_RECOVER); + if (c.overload > OVERLOAD_SECS) { + c.broken = true; + c.overload = 0; + c.load = 0; + this.events.push({ type: 'break', corner: c, anchorId: c.anchorId, hw: c.hw.name, t: this.t }); + } + } + if (this._dirtyRest) { this._applyRestLengths(); this._dirtyRest = false; } + } + + /** Re-rig a blown corner with fresh hardware. Lane D's hold-E repair calls this. */ + repairCorner(index, hw = HARDWARE[1]) { + const c = this.corners[index]; + if (!c || !c.broken) return false; + c.broken = false; + c.hw = hw; + c.load = 0; + c.overload = 0; + this._repin(this.t); + this.events.push({ type: 'repair', corner: c, anchorId: c.anchorId, hw: hw.name, t: this.t }); + return true; + } + + /** Turnbuckle trim at ONE corner (Lane D, 1.2 s hold). Tightens/eases locally. */ + trimCorner(index, delta) { + const c = this.corners[index]; + if (!c) return false; + c.trim = clamp(c.trim + delta, TRIM_MIN, TRIM_MAX); + this._dirtyRest = true; + return true; + } + + setTension(tension) { + this.tension = clamp(tension, TENSION_MIN, TENSION_MAX); + if (this.rigged) this._applyRestLengths(); + } + + drainEvents() { + const out = this.events; + this.events = []; + return out; + } + + /** + * Ground-projected shade over a rect: fraction of sample points on the rect + * that the sail blocks from the sun. This IS the shade mechanic, so it + * raycasts toward the actual sun rather than projecting straight down — + * which is what lets DESIGN.md's moving/seasonal sun change the answer. + * + * @param {object} rect {x, z, w, d} on the ground, metres, world XZ + * @param {object} sunDir direction TO the sun; defaults to straight overhead + */ + coverageOver(rect, sunDir = { x: 0, y: 1, z: 0 }) { + if (!this.rigged) return 0; + const len = Math.hypot(sunDir.x, sunDir.y, sunDir.z) || 1; + const dx = sunDir.x / len, dy = sunDir.y / len, dz = sunDir.z / len; + if (dy <= 0.01) return 0; // sun at or below the horizon casts no useful shade + + const COLS = 6, ROWS = 4; // prototype sampled 6x4 over the garden + let hit = 0; + for (let i = 0; i < COLS; i++) { + for (let j = 0; j < ROWS; j++) { + const ox = rect.x + ((i + 0.5) / COLS) * rect.w; + const oz = rect.z + ((j + 0.5) / ROWS) * rect.d; + if (this._rayHitsSail(ox, 0, oz, dx, dy, dz)) hit++; + } + } + return hit / (COLS * ROWS); + } + + /** Moller-Trumbore against every face; 162 tris, cheap enough to not bother accelerating. */ + _rayHitsSail(ox, oy, oz, dx, dy, dz) { + const pos = this.pos; + for (let i = 0; i < this.tris.length; i += 3) { + const a = this.tris[i] * 3, b = this.tris[i + 1] * 3, c = this.tris[i + 2] * 3; + const e1x = pos[b] - pos[a], e1y = pos[b + 1] - pos[a + 1], e1z = pos[b + 2] - pos[a + 2]; + const e2x = pos[c] - pos[a], e2y = pos[c + 1] - pos[a + 1], e2z = pos[c + 2] - pos[a + 2]; + const px = dy * e2z - dz * e2y, py = dz * e2x - dx * e2z, pz = dx * e2y - dy * e2x; + const det = e1x * px + e1y * py + e1z * pz; + if (Math.abs(det) < 1e-9) continue; // ray parallel to the face + const inv = 1 / det; + const tx = ox - pos[a], ty = oy - pos[a + 1], tz = oz - pos[a + 2]; + const u = (tx * px + ty * py + tz * pz) * inv; + if (u < 0 || u > 1) continue; + const qx = ty * e1z - tz * e1y, qy = tz * e1x - tx * e1z, qz = tx * e1y - ty * e1x; + const v = (dx * qx + dy * qy + dz * qz) * inv; + if (v < 0 || u + v > 1) continue; + const hitT = (e2x * qx + e2y * qy + e2z * qz) * inv; + if (hitT > 1e-6) return true; + } + return false; + } + + /** Sum of the aerodynamic + weight force on the whole sail, N. Used by the statics assert. */ + netAppliedForce(wind, t) { + this._accumulateWind(wind, t, SIM_DT); + let fx = 0, fy = 0, fz = 0; + for (let n = 0; n < this.invMass.length; n++) { + fx += this.force[n * 3]; + fy += this.force[n * 3 + 1] + GRAVITY * this.nodeMass; + fz += this.force[n * 3 + 2]; + } + return { x: fx, y: fy, z: fz }; + } + + maxLoad() { + let m = 0; + for (const c of this.corners) if (c.load > m) m = c.load; + return m; + } +} + +/** + * three.js view over a rig. Imported lazily so the sim core above stays + * headless-runnable; call this only from the browser, after Lane A's vendor/ + * exists. Returns a THREE.Group to add to the scene, with update() per frame. + */ +export async function createSailView(rig, { color = 0xd8c48a } = {}) { + const THREE = await import('../vendor/three.module.js'); + + const geo = new THREE.BufferGeometry(); + const verts = new Float32Array(rig.pos.length); + geo.setAttribute('position', new THREE.BufferAttribute(verts, 3)); + geo.setIndex(new THREE.BufferAttribute(new Uint16Array(rig.tris), 1)); + + const mat = new THREE.MeshStandardMaterial({ + color, side: THREE.DoubleSide, roughness: 0.92, metalness: 0.0, + }); + const mesh = new THREE.Mesh(geo, mat); + mesh.castShadow = true; // the shadow IS the product + mesh.receiveShadow = true; + mesh.frustumCulled = false; // it flogs well outside its initial bounds + + const group = new THREE.Group(); + group.add(mesh); + group.update = () => { + verts.set(rig.pos); + geo.attributes.position.needsUpdate = true; + geo.computeVertexNormals(); + geo.computeBoundingSphere(); + }; + group.update(); + return group; +} diff --git a/web/world/js/sail.selftest.js b/web/world/js/sail.selftest.js new file mode 100644 index 0000000..9d6ea5a --- /dev/null +++ b/web/world/js/sail.selftest.js @@ -0,0 +1,352 @@ +/** + * sail.selftest.js — assert suite for the sail sim. [Lane B] + * + * DOM-free and three-free on purpose: runs today under `node sail.selftest.js` + * before Lane A's shell exists, and Lane A's selftest.html can import + * runSailSelftest() unchanged once it lands. Drives the sim with fixed-dt loops + * only — never rAF, never a clock. + * + * The headline assert is `hypar sheds load vs flat`: it is the game's thesis + * stated as a test. If it ever goes red, the sail has stopped being a sail. + */ + +import { SailRig, HARDWARE } from './sail.js'; + +const SIM_DT = 1 / 60; + +// ---------- deterministic stub wind (Lane C owns the real weather.js) ---------- + +function mulberry32(seed) { + return function () { + seed |= 0; seed = (seed + 0x6d2b79f5) | 0; + let t = Math.imul(seed ^ (seed >>> 15), 1 | seed); + t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t; + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; +} + +/** + * Ports the prototype's gust scheduler shape: telegraph 1.5 s, ramp 0.8 s, + * hold ~1.7 s, fade 1 s. The schedule is precomputed from the seed so + * sample(pos, t) stays a pure function of t — which is what makes the + * determinism assert meaningful and lets the selftest fast-forward. + */ +function makeStubWind({ seed = 7, stormLen = 90, dir = { x: 0, y: 0, z: 1 }, calm = false } = {}) { + const rng = mulberry32(seed); + const gusts = []; + for (let t = 3; t < stormLen; t += 5 + rng() * 7) { + gusts.push({ start: t, pow: 12 + rng() * 16 + 10 * (t / stormLen) }); + } + const len = Math.hypot(dir.x, dir.y, dir.z) || 1; + const dx = dir.x / len, dy = dir.y / len, dz = dir.z / len; + const out = { x: 0, y: 0, z: 0 }; + + return { + speedAt(t) { + if (calm) return 0; + const p = Math.min(1, t / stormLen); + let speed = 8 + 26 * Math.min(1, p * 1.6); + for (const g of gusts) { + const gt = t - g.start; + if (gt < 0 || gt >= 5) continue; + if (gt < 1.5) continue; // telegraph: you see it, you don't feel it + else if (gt < 2.3) speed += g.pow * (gt - 1.5) / 0.8; // ramp + else if (gt < 4.0) speed += g.pow; // hold + else speed += g.pow * (5.0 - gt); // fade + } + return speed; + }, + sample(pos, t) { + const s = this.speedAt(t); + out.x = dx * s; out.y = dy * s; out.z = dz * s; + return out; + }, + }; +} + +const constantWind = (v) => ({ sample: () => v, speedAt: () => Math.hypot(v.x, v.y, v.z) }); + +// ---------- test rigs ---------- +// Same 5x5 m footprint, same multiset of corner heights {4.0, 4.0, 2.5, 2.5}. +// Only the ARRANGEMENT differs: coplanar (flat, pitched) vs permuted (twisted +// hypar). Any load difference is therefore purely geometry, nothing else. + +const FOOT = [ + { x: -2.5, z: -2.5 }, { x: 2.5, z: -2.5 }, { x: 2.5, z: 2.5 }, { x: -2.5, z: 2.5 }, +]; +const HEIGHTS_FLAT = [4.0, 4.0, 2.5, 2.5]; // y is linear in z -> one plane +const HEIGHTS_HYPAR = [4.0, 2.5, 4.0, 2.5]; // opposite corners up/down -> saddle + +const makeAnchors = (heights) => + FOOT.map((f, i) => ({ + id: `a${i}`, + type: 'post', + pos: { x: f.x, y: heights[i], z: f.z }, + })); + +const ALL_IDS = ['a0', 'a1', 'a2', 'a3']; +const UNBREAKABLE = { name: 'test rig', cost: 0, rating: Infinity }; + +function rig(heights, { hw = UNBREAKABLE, tension = 1.0, porosity = 0 } = {}) { + return new SailRig({ anchors: makeAnchors(heights), gridN: 10, porosity }) + .attach(ALL_IDS, [hw, hw, hw, hw], tension); +} + +/** Fixed-dt fast-forward. Returns peak corner load seen over the whole run, N. */ +function runStorm(r, wind, secs, onStep) { + const steps = Math.round(secs / SIM_DT); + let peak = 0; + for (let i = 0; i < steps; i++) { + r.step(SIM_DT, wind, i * SIM_DT); + const m = r.maxLoad(); + if (m > peak) peak = m; + if (onStep) onStep(r, i); + } + return peak; +} + +// ---------- tiny assert harness ---------- + +const results = []; +function test(name, fn) { + try { + const detail = fn(); + results.push({ name, pass: true, detail: detail || '' }); + } catch (e) { + results.push({ name, pass: false, detail: e.message }); + } +} +function assert(cond, msg) { + if (!cond) throw new Error(msg); +} +const kN = (n) => `${(n / 1000).toFixed(2)} kN`; + +// ---------- the suite ---------- + +export function runSailSelftest() { + results.length = 0; + + test('sim stays finite through a full storm', () => { + const r = rig(HEIGHTS_HYPAR); + runStorm(r, makeStubWind({ stormLen: 90 }), 90); + for (const v of r.pos) assert(Number.isFinite(v), 'node position went NaN/Infinity'); + for (const c of r.corners) assert(Number.isFinite(c.load), 'corner load went NaN'); + return `peak ${kN(r.corners.reduce((m, c) => Math.max(m, c.peakLoad), 0))}`; + }); + + test('sail sags under gravity when calm', () => { + const r = rig(HEIGHTS_FLAT); + runStorm(r, makeStubWind({ calm: true }), 6); + const N = r.N, mid = (Math.floor(N / 2) * N + Math.floor(N / 2)) * 3; + const midY = r.pos[mid + 1]; + const cornerMeanY = HEIGHTS_FLAT.reduce((a, b) => a + b) / 4; + assert(midY < cornerMeanY, `belly (${midY.toFixed(2)}m) should hang below corner mean (${cornerMeanY}m)`); + return `belly sags ${(cornerMeanY - midY).toFixed(2)} m below corner plane`; + }); + + // Newton's third law. This is what calibrates FABRIC_K into real newtons: + // if the corner reactions don't sum to the actual aerodynamic + weight force + // on the fabric, the load meter is lying and every kN rating is meaningless. + test('statics: corner reactions balance the applied force', () => { + const w = constantWind({ x: 0, y: 0, z: 18 }); + const r = rig(HEIGHTS_FLAT); + runStorm(r, w, 12); // settle + // A membrane in steady wind never fully stops moving, so this compares the + // TIME-AVERAGED reaction against the time-averaged applied force. That is + // the momentum balance that has to hold; instant by instant it need not. + let n = 0, ax = 0, ay = 0, az = 0, rx = 0, ry = 0, rz = 0; + const steps = Math.round(4 / SIM_DT); + for (let i = 0; i < steps; i++) { + const t = 12 + i * SIM_DT; + r.step(SIM_DT, w, t); + const f = r.netAppliedForce(w, t); + ax += f.x; ay += f.y; az += f.z; + for (const c of r.corners) { rx += c.loadVec.x; ry += c.loadVec.y; rz += c.loadVec.z; } + n++; + } + ax /= n; ay /= n; az /= n; rx /= n; ry /= n; rz /= n; + const appliedMag = Math.hypot(ax, ay, az); + const err = Math.hypot(rx - ax, ry - ay, rz - az) / appliedMag; + assert(err < 0.2, `reactions ${kN(Math.hypot(rx, ry, rz))} vs applied ${kN(appliedMag)} — ${(err * 100).toFixed(0)}% out of balance`); + return `applied ${kN(appliedMag)}, reactions ${kN(Math.hypot(rx, ry, rz))}, residual ${(err * 100).toFixed(1)}%`; + }); + + // THE THESIS. A twisted sail resists bellying into one coherent pocket, so + // its worst moment is gentler than a flat sail's worst moment. + // + // Scored on WORST CASE over wind direction, not per-direction. Lane C's + // storms veer, so the player never gets to choose the wind, and worst-case is + // what the hardware actually has to survive. Per-direction would be a false + // assert: a flat sail that happens to sit edge-on to the wind genuinely does + // have low drag, and from that one angle it beats the hypar. Demanding + // otherwise would mean tuning the sim into a lie. + test('hypar sheds load vs flat, worst case over wind direction (the thesis)', () => { + const DIRS = [ + { name: 'N', x: 0, z: 1 }, { name: 'NE', x: 0.707, z: 0.707 }, + { name: 'E', x: 1, z: 0 }, { name: 'SE', x: 0.707, z: -0.707 }, + { name: 'S', x: 0, z: -1 }, { name: 'SW', x: -0.707, z: -0.707 }, + { name: 'W', x: -1, z: 0 }, { name: 'NW', x: -0.707, z: 0.707 }, + ]; + const sweep = (heights) => { + let worst = 0, at = ''; + for (const d of DIRS) { + const storm = makeStubWind({ seed: 7, stormLen: 45, dir: { x: d.x, y: 0, z: d.z } }); + const p = runStorm(rig(heights), storm, 45); + if (p > worst) { worst = p; at = d.name; } + } + return { worst, at }; + }; + const flat = sweep(HEIGHTS_FLAT); + const hypar = sweep(HEIGHTS_HYPAR); + assert( + hypar.worst < flat.worst * 0.8, + `hypar worst ${kN(hypar.worst)} (${hypar.at}) should be well under flat worst ${kN(flat.worst)} (${flat.at})` + ); + const shed = (1 - hypar.worst / flat.worst) * 100; + return `flat worst ${kN(flat.worst)} from ${flat.at} -> hypar worst ${kN(hypar.worst)} from ${hypar.at} (sheds ${shed.toFixed(0)}%)`; + }); + + test('cascade: losing a corner spikes its neighbours', () => { + const w = constantWind({ x: 0, y: 0, z: 22 }); + const r = rig(HEIGHTS_HYPAR); + runStorm(r, w, 6); // settle + const before = Math.max(r.corners[1].load, r.corners[3].load); + r.corners[0].broken = true; + r._repin(r.t); + runStorm(r, w, 2.5); // let the load redistribute + const after = Math.max(r.corners[1].load, r.corners[3].load); + assert(after >= before * 2, `neighbour went ${kN(before)} -> ${kN(after)}, wanted >= 2x`); + return `neighbour ${kN(before)} -> ${kN(after)} (${(after / before).toFixed(1)}x)`; + }); + + test('determinism: identical inputs give byte-equal load traces', () => { + const trace = () => { + const r = rig(HEIGHTS_HYPAR); + const w = makeStubWind({ seed: 3, stormLen: 30 }); + const out = []; + runStorm(r, w, 30, (rr) => { for (const c of rr.corners) out.push(c.load); }); + return out; + }; + const a = trace(), b = trace(); + assert(a.length === b.length, 'traces differ in length'); + for (let i = 0; i < a.length; i++) { + assert(a[i] === b[i], `sample ${i} diverged: ${a[i]} vs ${b[i]}`); + } + return `${a.length} load samples identical`; + }); + + test('determinism: variable frame dt matches fixed dt', () => { + // Lane A's render loop delivers ragged dt. The internal accumulator has to + // absorb that, or nothing the selftest proves applies to the real game. + const w1 = makeStubWind({ seed: 5, stormLen: 20 }); + const fixed = rig(HEIGHTS_HYPAR); + for (let i = 0; i < Math.round(20 / SIM_DT); i++) fixed.step(SIM_DT, w1, i * SIM_DT); + + const w2 = makeStubWind({ seed: 5, stormLen: 20 }); + const ragged = rig(HEIGHTS_HYPAR); + const rng = mulberry32(99); + let acc = 0; + while (acc < 20) { + const dt = 0.004 + rng() * 0.02; // 4-24 ms frames + ragged.step(dt, w2, acc); + acc += dt; + } + for (let k = 0; k < 4; k++) { + const d = Math.abs(fixed.corners[k].load - ragged.corners[k].load); + assert(d < 1e-6, `corner ${k} drifted ${d.toFixed(6)} N between fixed and ragged dt`); + } + return 'ragged frame times converge on the fixed-dt trace'; + }); + + test('tension dial changes load (drum tight shock-loads)', () => { + const w = constantWind({ x: 0, y: 0, z: 20 }); + const loose = rig(HEIGHTS_HYPAR, { tension: 0.7 }); + const tight = rig(HEIGHTS_HYPAR, { tension: 1.35 }); + const loosePeak = runStorm(loose, w, 8); + const tightPeak = runStorm(tight, w, 8); + assert(tightPeak > loosePeak, `tight ${kN(tightPeak)} should exceed loose ${kN(loosePeak)}`); + return `loose ${kN(loosePeak)} vs tight ${kN(tightPeak)}`; + }); + + test('porous shade cloth carries less load than solid membrane', () => { + const w = constantWind({ x: 0, y: 0, z: 20 }); + const solid = runStorm(rig(HEIGHTS_HYPAR, { porosity: 0 }), w, 8); + const porous = runStorm(rig(HEIGHTS_HYPAR, { porosity: 0.35 }), w, 8); + assert(porous < solid, `porous ${kN(porous)} should be under solid ${kN(solid)}`); + return `solid ${kN(solid)} vs porous ${kN(porous)}`; + }); + + test('coverage: sail shades the ground under it, not beside it', () => { + const r = rig(HEIGHTS_FLAT); + runStorm(r, makeStubWind({ calm: true }), 4); + const under = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 }); + const beside = r.coverageOver({ x: 12, z: 12, w: 4, d: 4 }); + assert(under > 0.9, `ground under the sail only ${(under * 100).toFixed(0)}% shaded`); + assert(beside === 0, `ground 12 m away reported ${(beside * 100).toFixed(0)}% shaded`); + return `under sail ${(under * 100).toFixed(0)}%, off to the side ${(beside * 100).toFixed(0)}%`; + }); + + test('coverage tracks a low sun off to the side', () => { + const r = rig(HEIGHTS_FLAT); + runStorm(r, makeStubWind({ calm: true }), 4); + const noon = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 }, { x: 0, y: 1, z: 0 }); + const lowSun = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 }, { x: 0.9, y: 0.25, z: 0 }); + assert(noon > lowSun, `shadow should slide off the bed as the sun drops (noon ${noon}, low ${lowSun})`); + return `noon ${(noon * 100).toFixed(0)}% -> low sun ${(lowSun * 100).toFixed(0)}%`; + }); + + // PLAN3D §7 definition of done, in miniature. + test('cheap flat rig cascades; twisted mixed rig survives', () => { + const storm = () => makeStubWind({ seed: 11, stormLen: 90 }); + const cheap = rig(HEIGHTS_FLAT, { hw: HARDWARE[0], tension: 1.35 }); + runStorm(cheap, storm(), 90); + const cheapBroken = cheap.corners.filter((c) => c.broken).length; + + const good = rig(HEIGHTS_HYPAR, { hw: HARDWARE[2], tension: 0.95 }); + runStorm(good, storm(), 90); + const goodBroken = good.corners.filter((c) => c.broken).length; + + assert(cheapBroken >= 2, `flat drum-tight carabiner rig only lost ${cheapBroken} corners — should cascade`); + assert(goodBroken === 0, `twisted rated-shackle rig lost ${goodBroken} corners — should survive`); + return `cheap flat lost ${cheapBroken}/4, good hypar lost ${goodBroken}/4`; + }); + + test('repair re-pins a blown corner and it carries load again', () => { + const w = constantWind({ x: 0, y: 0, z: 20 }); + const r = rig(HEIGHTS_HYPAR); + runStorm(r, w, 4); + r.corners[0].broken = true; + r._repin(r.t); + runStorm(r, w, 1); + assert(r.corners[0].load === 0, 'broken corner should carry no load'); + assert(r.repairCorner(0, HARDWARE[1]), 'repairCorner should report success'); + runStorm(r, w, 3); + assert(r.corners[0].load > 100, `repaired corner only pulling ${kN(r.corners[0].load)}`); + const evs = r.drainEvents(); + assert(evs.some((e) => e.type === 'repair'), 'no repair event emitted'); + return `repaired corner back to ${kN(r.corners[0].load)}`; + }); + + const pass = results.every((r) => r.pass); + return { pass, results }; +} + +// ---------- entry points ---------- + +function report(out) { + const lines = out.results.map( + (r) => `${r.pass ? 'PASS' : 'FAIL'} ${r.name}${r.detail ? `\n ${r.detail}` : ''}` + ); + return `${lines.join('\n')}\n\n${out.pass ? 'ALL GREEN' : 'FAILURES'} — ${out.results.filter((r) => r.pass).length}/${out.results.length}`; +} + +// Run only when invoked directly: `node web/world/js/sail.selftest.js`. +// Importing the module must not run the suite. In the browser `process` is +// undefined, so Lane A's selftest.html just calls runSailSelftest() itself. +if (typeof process !== 'undefined' && process.versions?.node && import.meta.filename === process.argv[1]) { + const out = runSailSelftest(); + console.log(report(out)); + process.exit(out.pass ? 0 : 1); +} + +export { report, makeStubWind, HEIGHTS_FLAT, HEIGHTS_HYPAR, makeAnchors };