Censer
A bronze censer hangs on a long chain and breathes incense from its vents. On the CPU, its heavy body swings rigid on three short chains and a long one, in small position-based steps: lift it over the pivot and the chain falls slack, drop it and the chain snaps taut. Drag it to throw it, and a gentle swing returns once you leave it alone. The smoke is a 3D grid of air stepped on the GPU after Crane, Llamas and Tariq's chapter in GPU Gems 3, with MacCormack advection, vorticity confinement and a multigrid pressure solve. The moving censer pushes the air, which carries the smoke on a grid twice as fine. At 120 Hz, each frame records half a step. Each step lights the smoke from one key light through itself and the censer, and each frame marches it at half resolution in steps that shorten at its edges, skipping empty bricks of air, before laying it over the censer by depth. The embers inside light the smoke near the lid. The smoke setting scales what the vents breathe out, from a faint thread to a billowing cloud.
Files, in reading order
- scene.tsCreates the scene: its resources, its settings and what each frame records.
- censer.tsA thurible hanging by its chain: its pose, the solid and smoke it gives the air, and a hand that may take hold of it and throw it.
- march.tsThe smoke lit by one key light from the scene's sun and by the embers in the censer, marched at half the output's size through each pixel's view up to the scene's depth and laid over the picture by depth.
- model.tsThe censer's bronze body and its chains, swinging with
censer. - smoke.tsA box of smoke on a grid of cells: the air's velocity, advected by MacCormack's method, lifted by heat, kept swirling by vorticity confinement and made incompressible by a multigrid pressure solve, carrying smoke and heat on a grid twice as fine.
- chain.tsA heavy body hanging by three legs of chain from a crown, and the crown by a chain from a fixed pivot, stepped by extended position-based dynamics in small steps.
- cells.ts
- pressure.tsSolves the pressure's Poisson equation on a box of cells by geometric multigrid: a V-cycle of red–black Gauss–Seidel sweeps over grids halved down to a coarse one, so a push reaches across the whole box in a step.
Kit modules cameras, clock, finish, lighting, materials, memo, owner, passes, pipelines, shading, targets, tone. Also @typegpu/noise, typegpu.
scene.ts
import { cameras, clock, finish, lighting, owner, passes, pipelines, targets, tone, type Frame, type Surface,} from '@typegpu/kit';import type { TgpuRoot } from 'typegpu';import { createCenser, STEP } from './censer.ts';import { createMarch } from './march.ts';import { createModel } from './model.ts';import type { Settings } from '../settings.ts';import { createSmoke } from './smoke.ts';/** The smoke's box around the swing: cells, their size and lowest corner. */const CELLS = [112, 112, 88] as const;const CELL = 0.04;const CORNER = [-2.24, 0.24, -1.76] as const;/** Radians the camera turns for a drag across the whole surface. */const ORBIT = 3;export interface SceneOptions { /** The smoke the vents breathe out, 1 by default. */ readonly emission?: number;}export async function createScene(root: TgpuRoot, surface: Surface, options: SceneOptions = {}) { const own = owner(); const screen = own( targets.create(root, surface, { hdr: { format: 'rgba16float', sampled: true }, depth: { format: 'depth24plus', sampled: true }, }), ); const view = own( cameras.orbit(root, { size: surface, target: [0, 2.3, 0], distance: 4.6, yaw: 0, pitch: 0.05, fov: Math.PI / 4, pitchRange: [-0.3, 0.9], distanceRange: [2, 9], near: 0.05, far: 30, }), ); const light = own( lighting.create(root, { settings: { sunDirection: [-1, 1.2, 0.35], sunColor: [1, 0.97, 0.92], sunIntensity: 2.4, skyTop: [0.3, 0.34, 0.4], // The embers and the lit smoke warm the body's underside skyBottom: [0.9, 0.6, 0.38], skyIntensity: 0.08, skyReflection: 0.5, }, }), ); const scene = own( passes .create( root, screen, { hdr: passes.color, depth: passes.depth }, { clear: [0.003, 0.003, 0.0035, 1] }, ) .pipe(view.inject()) .beneath(light.inject()), ); const censer = createCenser(root); const model = own(await createModel(root, censer)); model.update(); const smoke = own( createSmoke(root, { cells: CELLS, origin: CORNER, cell: CELL, solid: censer.solid, source: censer.source, force: censer.force, carry: censer.carry, sponge: 3, settings: { dt: STEP, buoyancy: 2.5, weight: 0.02, vorticity: 16, smokeDecay: 0.15, heatDecay: 0.4, drag: 0.05, }, }), ); const air = own( createMarch(root, smoke, censer, { view, light, size: surface, picture: () => screen.textures.hdr, depth: () => screen.textures.depth, }), ); const picture = own( finish.create(root, { format: surface.format, source: () => air.output, tone: tone.aces, settings: { exposure: 0.6, vignette: 0.25, dither: 0 }, }), ); const drawn = [scene.add(model.drawing)]; await pipelines.compile([smoke, scene, air, picture]); // Half steps, so a frame at 120 Hz records one half of a smoke step and one at 60 Hz both. A // frame that misses the display records no catch-up on top of its own half const sim = clock.fixedStep({ step: STEP / 2, maxSteps: 4, steady: true }); let advecting = true; let time = 0; // A frame that records only the second half moves the chain half a step on, so both frames at // 120 Hz carry half its cost, but only the first half shows the censer where it is let ahead = 0; // A stroke that starts on the censer takes hold of it; any other orbits the camera const ray = { origin: [0, 0, 0], direction: [0, 0, 0] }; let holding = false; const last = [0, 0]; const values = { emission: options.emission ?? 1 }; const settings: Settings<typeof values> = { values, knobs: { emission: { label: 'smoke', min: 0.25, max: 4, step: 0.05 } }, write() { censer.emit(values.emission); }, }; return { view, settings, pointing: { drag(x: number, y: number, start: boolean) { view.rayAt(x, y, ray); if (start) holding = censer.grab(ray.origin, ray.direction); else if (holding) censer.pull(ray.origin, ray.direction); else cameras.turn(view.goal, x - last[0]!, y - last[1]!, { sensitivity: ORBIT }); last[0] = x; last[1] = y; }, end() { if (holding) censer.release(); holding = false; }, }, frame(frame: Frame) { const halves = sim.advance(frame.dt); if (halves > 0) { const pass = frame.encoder.beginComputePass(); let moved = false; for (let i = 0; i < halves; i++) { if (advecting) { time += STEP; censer.swing(STEP - ahead, time); model.update(); ahead = 0; smoke.advect(pass); moved = true; } else { if (halves === 1) { censer.chain.advance(STEP / 2); ahead = STEP / 2; } smoke.project(pass); } advecting = !advecting; } // The light changes only with the smoke and the censer, which move in the first half if (moved) air.light(pass); pass.end(); } screen.fit(); view.update(frame.dt); scene.record(frame, drawn); air.record(frame); picture.record(frame, surface.view(frame)); }, destroy: own.destroy, };}censer.ts
import { perlin2d } from '@typegpu/noise';import { d, std, type TgpuRoot } from 'typegpu';import { createChain, LENGTH, SIZE } from './chain.ts';type Vec3 = [number, number, number];type Mat3 = [number, number, number, number, number, number, number, number, number];/** The body as a capsule along the chain: half its length and its radius. */const HALF = 0.32 * SIZE;const RADIUS = 0.5 * SIZE;/** How near the body's middle a pointer's ray must pass to take hold of it, in metres. */const REACH = 0.32;/** * The vents breathe out in puffs from seven places around the top of the lid, each a soft blob * this wide, this far out from the chain and this high above the body's middle, just clear of its * solid. */const VENTS = 7;const PUFF = 0.085;const VENT_OUT = 0.13;const VENT_UP = 0.37;/** * Smoke and heat at the middle of a puff at the default emission, and its jet's acceleration, * per second. */const SMOKE = 1.3;const HEAT = 17;const JET = 30;/** * Air within this many metres of the body is dragged along with it, at up to this rate per second, * and air in a puff takes on its vent's motion at this rate, so the censer stirs what it passes. */const SHELL = 0.1;const SHELL_RATE = 40;const INHERIT = 30;/** How fast each vent's breath comes and goes, per second. */const BREATHING = 2.3;/** The smoke's step, in seconds, over which each puff is swept along the body's path. */export const STEP = 1 / 60;const Pose = d.struct({ /** Turns the model of the hanging censer to the chain's direction. */ rotation: d.mat3x3f, /** The middle of the body, and its velocity. */ middle: d.vec3f, velocity: d.vec3f, /** The body's angular velocity, in radians per second. */ spin: d.vec3f, /** The smoke's time in seconds, which the vents breathe to. */ time: d.f32, /** The smoke and heat the vents add per second at the middle of a puff. */ vents: d.vec2f,});/** How far a point of the body's frame lies outside the body. */const distanceTo = (body: d.v3f) => { 'use gpu'; const along = d.vec3f(0, std.clamp(body.y, -HALF, HALF), 0); return std.length(body - along) - RADIUS;};const dot = (a: ArrayLike<number>, b: ArrayLike<number>) => a[0]! * b[0]! + a[1]! * b[1]! + a[2]! * b[2]!;/** * A thurible hanging by its chain: its pose, the solid and smoke it gives the air, and a hand that * may take hold of it and throw it. */export function createCenser(root: TgpuRoot) { const chain = createChain(); // Plain arrays, which TypeGPU writes as they are, where setting a vector's components would make // garbage that Hermes collects in frames const state = { rotation: [1, 0, 0, 0, 1, 0, 0, 0, 1] as Mat3, middle: chain.middle, velocity: chain.velocity, spin: chain.spin, time: 0, vents: [SMOKE, HEAT] as [number, number], }; const pose = root.createUniform(Pose).$name('censer pose'); // The plane the hand moves in, through where it took hold and facing along the ray, and where it // pulls to const anchor: Vec3 = [0, 0, 0]; const facing: Vec3 = [0, 0, 1]; const pulled: Vec3 = [0, 0, 0]; /** Writes the chain's state as the body's pose. */ const write = (time: number) => { const x = chain.turn[0]; const y = chain.turn[1]; const z = chain.turn[2]; const w = chain.turn[3]; // Column by column const r = state.rotation; r[0] = 1 - 2 * (y * y + z * z); r[1] = 2 * (x * y + z * w); r[2] = 2 * (x * z - y * w); r[3] = 2 * (x * y - z * w); r[4] = 1 - 2 * (x * x + z * z); r[5] = 2 * (y * z + x * w); r[6] = 2 * (x * z + y * w); r[7] = 2 * (y * z - x * w); r[8] = 1 - 2 * (x * x + y * y); state.time = time; pose.write(state); }; /** A point of the body, measured from its middle, along and across the chain. */ const bodyOf = (position: d.v3f) => { 'use gpu'; return std.transpose(pose.$.rotation) * (position - pose.$.middle); }; /** * The puff of the vent nearest a point: its jet's direction in the body's frame, and how much of * it reaches the point. Each vent breathes on its own rhythm and sweeps its puff along the way * it moved over the last step, so a fast swing leaves a trail rather than a string of beads. */ const puffAt = (position: d.v3f) => { 'use gpu'; const body = bodyOf(position); const swept = std.transpose(pose.$.rotation) * pose.$.velocity * STEP; const margin = PUFF * 3 + std.length(swept); const rim = std.length(body.xz) - VENT_OUT; if (std.abs(body.y - VENT_UP) > margin || std.abs(rim) > margin) return d.vec4f(); const around = std.atan2(body.z, body.x) * (VENTS / (2 * Math.PI)); const vent = std.round(around); const angle = vent * ((2 * Math.PI) / VENTS); const out = d.vec3f(std.cos(angle), 0, std.sin(angle)); const middle = out * VENT_OUT + d.vec3f(0, VENT_UP, 0); const length2 = std.max(std.dot(swept, swept), 1e-8); const back = std.saturate(std.dot(middle - body, swept) / length2); const offset = body - (middle - swept * back); const spread = PUFF / (PUFF + std.sqrt(length2) * 0.56); const reach = std.exp(-std.dot(offset, offset) / (PUFF * PUFF)) * spread; const breath = perlin2d.sample(d.vec2f(pose.$.time * BREATHING, vent * 7.3 + 11)); const pulse = std.saturate(breath * 2 + 0.5); const jet = std.normalize(out * 0.3 + d.vec3f(0, 1.5, 0)); return d.vec4f(jet, reach * pulse * pulse); }; return { chain, /** Places a point of the model, given from the pivot as the censer hangs at rest. */ place: (position: d.v3f) => { 'use gpu'; return pose.$.rotation * (position + d.vec3f(0, LENGTH, 0)) + pose.$.middle; }, /** Turns a direction of the model as the censer hangs. */ turn: (direction: d.v3f) => { 'use gpu'; return pose.$.rotation * direction; }, /** The middle of the body, where the embers burn. */ middle: () => { 'use gpu'; return d.vec3f(pose.$.middle); }, /** How far a point lies outside the censer's body, in metres. */ distance: (position: d.v3f) => { 'use gpu'; return distanceTo(bodyOf(position)); }, /** The censer's body as a solid, moving as it does. */ solid: (position: d.v3f) => { 'use gpu'; if (distanceTo(bodyOf(position)) > 0) return d.vec4f(); const moving = pose.$.velocity + std.cross(pose.$.spin, position - pose.$.middle); return d.vec4f(moving, 1); }, /** * The motion air near the body is dragged toward: its own, in a thin shell around it and in * each puff. */ carry: (position: d.v3f) => { 'use gpu'; const outside = distanceTo(bodyOf(position)); const shell = std.saturate(1 - outside / SHELL) * std.step(0, outside); const rate = SHELL_RATE * shell * shell + INHERIT * puffAt(position).w; const moving = pose.$.velocity + std.cross(pose.$.spin, position - pose.$.middle); return d.vec4f(moving, rate); }, /** Smoke and heat from the vents in the lid, in puffs. */ source: (position: d.v3f) => { 'use gpu'; const puff = puffAt(position); return pose.$.vents * puff.w; }, /** The jet each puff leaves its vent with, up and out of the lid. */ force: (position: d.v3f) => { 'use gpu'; const puff = puffAt(position); return pose.$.rotation * puff.xyz * (JET * puff.w); }, /** * Scales the smoke and heat the vents breathe out, 1 by default. Hotter puffs rise faster, so * more smoke billows higher rather than piling up. */ emit(rate: number) { state.vents[0] = SMOKE * rate; state.vents[1] = HEAT * rate; }, /** Swings the censer `dt` seconds on, to `time`. */ swing(dt: number, time: number) { chain.advance(dt); write(time); }, /** * Takes hold of the body if a pointer's ray passes near its middle, and returns whether it did. * The hand then moves in the plane through the grip that faces along the ray. */ grab(origin: ArrayLike<number>, direction: ArrayLike<number>) { const middle = chain.middle; const to = [middle[0] - origin[0]!, middle[1] - origin[1]!, middle[2] - origin[2]!]; const along = dot(to, direction); if (along <= 0 || dot(to, to) - along * along > REACH * REACH) return false; for (let i = 0; i < 3; i++) { facing[i] = direction[i]!; anchor[i] = origin[i]! + direction[i]! * along; } chain.hold(chain.localOf(anchor)); return true; }, /** Moves the hand to where a pointer's ray meets its plane. */ pull(origin: ArrayLike<number>, direction: ArrayLike<number>) { const square = dot(direction, facing); if (!chain.held || square < 1e-3) return; for (let i = 0; i < 3; i++) pulled[i] = anchor[i]! - origin[i]!; const t = dot(pulled, facing) / square; for (let i = 0; i < 3; i++) pulled[i] = origin[i]! + direction[i]! * t; chain.pull(pulled); }, /** Lets go, leaving the censer swinging as the hand threw it. */ release() { chain.release(); }, };}export type Censer = ReturnType<typeof createCenser>;march.ts
import { cameras, memo, shading, targets, type CameraBase, type Frame, type Lighting, type SizeSource, type finish,} from '@typegpu/kit';import { randf } from '@typegpu/noise';import { common, d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';import type { Censer } from './censer.ts';import { dispatch, groups, kernel, volume } from './cells.ts';import type { Smoke } from './smoke.ts';/** Cells along each edge of a brick, the unit the march skips empty air in. */const BRICK = 8;/** * Smoke stops light more than in proportion to how much there is, so a dense core shades into * wisps that let the dark through: `EXTINCTION` per metre at `REFERENCE` smoke, rising as its * `CURVE` power, and easing off toward `DENSEST` per metre so even the core shows the shadows * inside it. Smoke under `EMPTY` stops nothing, and bricks holding no more are skipped. */const EXTINCTION = 2;const REFERENCE = 0.05;const CURVE = 1.5;const DENSEST = 6;const EMPTY = 0.002;/** Toward the key, smoke shades as if this many times denser, so its folds show light and shadow. */const SHADOW = 3;/** How much light reaches a point through smoke: a share that scatters on through thick smoke. */const DEEP = 0.2;const DEEP_SHARE = 0.2;/** The censer's body stops light this much per metre, its edge softened over `SOFT` metres. */const BODY = 60;const SOFT = 0.06;/** A light cell looks toward the key in steps of this many cells, at most `LIGHT_STEPS` of them. */const LIGHT_STRIDE = 2;const LIGHT_STEPS = 128;/** * The march's steps in metres: short where the density changes fast, so edges stay sharp, and up * to `FAR` through even smoke, each changing it by about `EDGE` per metre. Never more than * `MOST_STEPS` a pixel. */const NEAR = 0.01;const FAR = 0.04;const EDGE = 0.4;const MOST_STEPS = 400;/** Smoke fades out this far inside the box's open faces, so it leaves without an edge. */const FADE = 0.24;/** A broad lobe and a forward one, with its share, so smoke before the light glows at its rim. */const BROAD = 0.2;const FORWARD = 0.75;const FORWARD_SHARE = 0.15;/** The faint cool light smoke takes in its shade, as a share of the key's. */const SHADE = d.vec3f(0.035, 0.04, 0.05);/** The embers' light on the smoke near the vents: colour, strength and reach in metres. */const EMBER = d.vec3f(1, 0.36, 0.09);const EMBER_POWER = 0.5;const EMBER_REACH = 0.2;/** * The march runs at this fraction of the output's size, each of its pixels marching up to the * nearest surface among the output pixels it covers, and an output pixel takes the marched ones * around it whose depths match its own within `MATCH` of its distance. */const REDUCE = 2;const MATCH = 0.03;const marchLayout = tgpu.bindGroupLayout({ depth: { texture: d.textureDepth2d() },});const compositeLayout = tgpu.bindGroupLayout({ picture: { texture: d.texture2d(d.f32) }, depth: { texture: d.textureDepth2d() }, marched: { texture: d.texture2d(d.f32) }, nearest: { texture: d.texture2d(d.f32), sampleType: 'unfilterable-float' },});/** Henyey and Greenstein's phase function, times 4π so light scattered evenly is 1. */const lobe = (cosine: number, g: number) => { 'use gpu'; return (1 - g * g) / std.pow(1 + g * g - 2 * g * cosine, 1.5);};/** The density that stops light at a point, per metre, from the smoke it holds. */const densityOf = (smoke: number) => { 'use gpu'; const thick = std.max(smoke - EMPTY, 0) / REFERENCE; const linear = std.pow(thick, CURVE) * EXTINCTION; return DENSEST * (1 - std.exp(-linear / DENSEST));};/** * The smoke lit by one key light from the scene's sun and by the embers in the censer, marched at * half the output's size through each pixel's view up to the scene's depth and laid over the * picture by depth. */export function createMarch( root: TgpuRoot, smoke: Smoke, censer: Censer, options: { view: CameraBase; light: Lighting; size: SizeSource; picture: () => finish.FloatTexture; depth: () => targets.Depth; },) { const { view, light, size, picture, depth } = options; const [nx, ny, nz] = smoke.cells; const bricks = [nx / BRICK, ny / BRICK, nz / BRICK] as const; const output = targets.create(root, size, { hdr: { format: 'rgba16float', sampled: true } }); const reduced = targets.create(root, size, { hdr: { format: 'rgba16float', sampled: true, scale: 1 / REDUCE }, nearest: { format: 'r32float', sampled: true, scale: 1 / REDUCE }, }); const smooth = root.createSampler({ minFilter: 'linear', magFilter: 'linear' }); // Per cell of the air's grid the key light let through to it and the embers' light; per brick // the most smoke near it const lit = volume(root, smoke.cells, 'rgba16float', 'march light'); const occupied = volume(root, bricks, 'rgba16float', 'march bricks'); const low = d.vec3f(...smoke.low); const high = d.vec3f(...smoke.high); const { cell } = smoke; const extent = d.vec3f(nx * cell, ny * cell, nz * cell); const brickSize = cell * BRICK; const brickCount = d.vec3i(...bricks); /** How much of the smoke at a point shows, fading out toward the box's open faces. */ const fadeAt = (position: d.v3f) => { 'use gpu'; const inside = std.min(position - low, high - position); const edge = std.min(std.min(inside.x, inside.y), inside.z); return std.smoothstep(0, FADE, edge); }; /** The smoke at a point, faded out toward the box's open faces. */ const smokeAt = (position: d.v3f) => { 'use gpu'; return smoke.smokeAt(position) * fadeAt(position); }; /** * The smoke averaged over the air's cells at a point, sampled between them, so light and bricks * need not read every finer cell. */ const coarseAt = (position: d.v3f) => { 'use gpu'; const uvw = (position - low) / extent; return std.textureSampleLevel(smoke.averaged.$, smooth.$, uvw, 0).x * fadeAt(position); }; /** * For each cell, the key light that reaches it through the smoke and the censer, and the embers'. * Only cells in bricks holding smoke look toward the key: the march reads no others. */ const lightCells = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; const at = low + (d.vec3f(id) + 0.5) * cell; const ember = censer.middle() + censer.turn(d.vec3f(0, 0.18, 0)); const near = std.length(at - ember) / EMBER_REACH; const glow = EMBER_POWER / (1 + near * near); if (std.textureLoad(occupied.read.$, d.vec3i(id) / BRICK, 0).x < EMPTY) { std.textureStore(lit.write.$, id, d.vec4f(1, glow, 0, 0)); return; } const toward = shading.towardsSun(); let depthOf = d.f32(0); for (let i = d.i32(0); i < LIGHT_STEPS; i++) { const q = at + toward * ((d.f32(i) + 0.5) * cell * LIGHT_STRIDE); if (std.any(std.gt(q, high)) || std.any(std.lt(q, low))) break; // The censer's body casts a shadow, softened over a few centimetres at its edge const body = std.saturate(0.5 - censer.distance(q) / SOFT); depthOf += (densityOf(coarseAt(q)) * SHADOW + body * BODY) * cell * LIGHT_STRIDE; } const through = std.mix(std.exp(-depthOf), std.exp(-depthOf * DEEP), DEEP_SHARE); std.textureStore(lit.write.$, id, d.vec4f(through, glow, 0, 0)); }); /** For each brick, the most smoke in it and a cell around it, where the march must look. */ const occupancy = kernel(({ id }) => { 'use gpu'; if (id.x >= bricks[0] || id.y >= bricks[1] || id.z >= bricks[2]) return; const first = d.vec3i(id) * BRICK - 1; const last = d.vec3i(nx - 1, ny - 1, nz - 1); let most = d.f32(0); for (let z = d.i32(0); z < BRICK + 2; z++) { for (let y = d.i32(0); y < BRICK + 2; y++) { for (let x = d.i32(0); x < BRICK + 2; x++) { const at = std.clamp(first + d.vec3i(x, y, z), d.vec3i(), last); const shown = fadeAt(low + (d.vec3f(at) + 0.5) * cell); most = std.max(most, std.textureLoad(smoke.averaged.$, at, 0).x * shown); } } } std.textureStore(occupied.write.$, id, d.vec4f(most, 0, 0, 0)); }); /** The view distance of the nearest surface among the output pixels a reduced one covers. */ const nearestOf = (reducedPixel: d.v2i) => { 'use gpu'; const last = d.vec2i(std.textureDimensions(marchLayout.$.depth)) - 1; let nearest = d.f32(3.4e38); for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const pixel = std.min(reducedPixel * REDUCE + d.vec2i(x, y), last); nearest = std.min( nearest, cameras.distanceOf(cameras.current.$, std.textureLoad(marchLayout.$.depth, pixel, 0)), ); } } return nearest; }; const march = tgpu.fragmentFn({ in: { uv: d.vec2f, position: d.builtin.position }, // The smoke a reduced pixel marched, and the view distance it marched up to in `x` out: { smoke: d.vec4f, nearest: d.vec4f }, })(({ uv, position }) => { 'use gpu'; const camera = cameras.current.$; const ray = cameras.rayOf(camera, uv * d.vec2f(2, -2) + d.vec2f(-1, 1)); const nearest = nearestOf(d.vec2i(position.xy)); const end = nearest / std.max(0.1, std.dot(ray.direction, cameras.forwardOf(camera))); // Where the view crosses the smoke's box const inverse = d.vec3f(1) / ray.direction; const near = (low - ray.origin) * inverse; const far = (high - ray.origin) * inverse; const enter = std.max( std.max(std.min(near.x, far.x), std.min(near.y, far.y)), std.max(std.min(near.z, far.z), 0), ); const exit = std.min( std.min(std.max(near.x, far.x), std.max(near.y, far.y)), std.min(std.max(near.z, far.z), end), ); if (exit <= enter) return { smoke: d.vec4f(0, 0, 0, 1), nearest: d.vec4f(nearest, 0, 0, 0) }; randf.seed2(std.floor(position.xy)); const offset = randf.sample(); const env = shading.settings(); const cosine = std.dot(ray.direction, shading.towardsSun()); const phase = std.mix(lobe(cosine, BROAD), lobe(cosine, FORWARD), FORWARD_SHARE); const key = env.sunColor * (env.sunIntensity * phase); const shade = SHADE * env.sunIntensity; let gathered = d.vec3f(); let through = d.f32(1); let t = enter + offset * NEAR; let step = d.f32(NEAR); let before = d.f32(0); for (let i = d.i32(0); i < MOST_STEPS; i++) { if (t >= exit || through < 0.004) break; const world = ray.origin + ray.direction * t; const brick = std.clamp( d.vec3i(std.floor((world - low) / brickSize)), d.vec3i(), brickCount - 1, ); if (std.textureLoad(occupied.read.$, brick, 0).x < EMPTY) { // Skip to where the view leaves the brick const corner = low + d.vec3f(brick) * brickSize; const a = (corner - ray.origin) * inverse; const b = (corner + brickSize - ray.origin) * inverse; const out = std.min(std.min(std.max(a.x, b.x), std.max(a.y, b.y)), std.max(a.z, b.z)); t = std.max(out, t) + offset * NEAR + 1e-4; step = NEAR; before = 0; continue; } const density = densityOf(smokeAt(world).x); step = std.clamp((EDGE * step) / (std.abs(density - before) + 1e-4), NEAR, FAR); before = density; if (density > 1e-3) { // The last step is cut where the view meets the censer or leaves the box const length = std.min(step, exit - t); const here = std.textureSampleLevel(lit.read.$, smooth.$, (world - low) / extent, 0); const fade = std.exp(-density * length); const scattered = key * here.x + shade + EMBER * here.y; gathered += scattered * (through * (1 - fade)); through *= fade; } t += step; } return { smoke: d.vec4f(gathered, through), nearest: d.vec4f(nearest, 0, 0, 0) }; }); /** The smoke marched at the reduced size over the picture, at the output's. */ const composite = tgpu.fragmentFn({ in: { position: d.builtin.position }, out: d.vec4f, })(({ position }) => { 'use gpu'; const pixel = d.vec2i(position.xy); const seen = std.textureLoad(compositeLayout.$.picture, pixel, 0).xyz; const own = cameras.distanceOf( cameras.current.$, std.textureLoad(compositeLayout.$.depth, pixel, 0), ); const last = d.vec2i(std.textureDimensions(compositeLayout.$.marched)) - 1; const at = position.xy / REDUCE - 0.5; const base = d.vec2i(std.floor(at)); const f = at - std.floor(at); let sum = d.vec4f(); let weight = d.f32(0); let best = d.vec4f(0, 0, 0, 1); let closest = d.f32(3.4e38); for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const marched = std.clamp(base + d.vec2i(x, y), d.vec2i(), last); const gap = std.abs(std.textureLoad(compositeLayout.$.nearest, marched, 0).x - own); const value = std.textureLoad(compositeLayout.$.marched, marched, 0); const near = std.mix(1 - f, f, d.vec2f(x, y)); const w = near.x * near.y * std.exp(-gap / (MATCH * own)); sum += value * w; weight += w; if (gap < closest) { closest = gap; best = d.vec4f(value); } } } const laid = std.select(best, sum / weight, weight > 1e-4); return d.vec4f(seen * laid.w + laid.xyz, 1); }); const lit3d = root.pipe(view.inject()).pipe(light.inject()); const pipelines = { light: lit3d.createComputePipeline({ compute: lightCells }), occupancy: root.createComputePipeline({ compute: occupancy }), march: lit3d.createRenderPipeline({ vertex: common.fullScreenTriangle, fragment: march, targets: { smoke: { format: 'rgba16float' }, nearest: { format: 'r32float' } }, }), composite: root.pipe(view.inject()).createRenderPipeline({ vertex: common.fullScreenTriangle, fragment: composite, targets: { format: 'rgba16float' }, }), }; const marching = memo((deep: targets.Depth, matter: Smoke['matter']) => pipelines.march.with(root.createBindGroup(marchLayout, { depth: deep })).with(matter), ); const compositing = memo( (color: finish.FloatTexture, deep: targets.Depth, marched: typeof reduced.textures) => pipelines.composite.with( root.createBindGroup(compositeLayout, { picture: color, depth: deep, marched: marched.hdr, nearest: marched.nearest, }), ), ); const attachments = { colorAttachments: { view: output.views.hdr } }; const reducedAttachments = { colorAttachments: [{ view: reduced.views.hdr }, { view: reduced.views.nearest }], }; const cellGroups = groups(smoke.cells); const brickGroups = groups(bricks); return { get output() { return output.textures.hdr; }, initAsync: () => Promise.all(Object.values(pipelines).map((each) => each.initAsync())), /** Records finding the bricks the smoke fills and lighting it, after the smoke's `advect`. */ light(pass: TgpuComputePass) { pass.setPipeline(pipelines.occupancy); dispatch(pass, brickGroups); pass.setPipeline(pipelines.light); dispatch(pass, cellGroups); }, /** Marches the smoke at the reduced size, then lays it over the picture into `output`. */ record(frame: Frame) { reduced.fit(); output.fit(); reducedAttachments.colorAttachments[0]!.view = reduced.views.hdr; reducedAttachments.colorAttachments[1]!.view = reduced.views.nearest; const marchPass = frame.encoder.beginRenderPass(reducedAttachments); marchPass.setPipeline(marching(depth(), smoke.matter)); marchPass.draw(3); marchPass.end(); attachments.colorAttachments.view = output.views.hdr; const pass = frame.encoder.beginRenderPass(attachments); pass.setPipeline(compositing(picture(), depth(), reduced.textures)); pass.draw(3); pass.end(); }, destroy() { output.destroy(); reduced.destroy(); lit.texture.destroy(); occupied.texture.destroy(); }, };}model.ts
import { cameras, materials, passes, shading } from '@typegpu/kit';import { perlin3d } from '@typegpu/noise';import { meshes } from '@typegpu/kit/geometry';import { d, std, tgpu, type TgpuRoot } from 'typegpu';import type { Censer } from './censer.ts';import { CROWN, LENGTH, PIVOT, PROFILE, SIZE } from './chain.ts';const Corner = d.struct({ position: d.vec4f, world: d.vec3f, normal: d.vec3f, local: d.vec3f, iron: d.u32,});const Seen = d.struct({ color: d.vec4f });/** A chain link's length, width and wire radius, in metres, and the length each adds to a chain. */const LINK = [0.042, 0.025, 0.0045] as const;const PITCH = LINK[0] - LINK[2] * 2;/** How far the chain runs on up from the pivot, past the top of the view. */const ABOVE = 1.2;/** Points around the body, and around each half of a link and across its wire. */const SIDES = 72;const BEND = 8;const WIRE = 6;/** The profile's points with their normals, radius and height in `xy` and the normal in `zw`. */const OUTLINE = tgpu.const( d.arrayOf(d.vec4f, PROFILE.length), PROFILE.map(([r, y], i) => { const before = PROFILE[Math.max(0, i - 1)]!; const after = PROFILE[Math.min(PROFILE.length - 1, i + 1)]!; const dr = after[0] - before[0]; const dy = after[1] - before[1]; const length = Math.hypot(dr, dy) || 1; return d.vec4f(r, y, -dy / length, dr / length); }),);const outlineAt = (v: number) => { 'use gpu'; return d.vec4f(OUTLINE.$[d.u32(std.round(v * (PROFILE.length - 1)))]!);};const aroundAt = (u: number) => { 'use gpu'; const a = u * 2 * Math.PI; return d.vec3f(std.cos(a), 1, std.sin(a));};/** The censer's body, the profile turned about the chain, its middle `LENGTH` below the pivot. */const body = meshes.parametric( { at: (u, v) => { 'use gpu'; return aroundAt(u) * outlineAt(v).xyx * SIZE - d.vec3f(0, LENGTH, 0); }, normalAt: (u, v) => { 'use gpu'; return aroundAt(u) * outlineAt(v).zwz; }, }, { cols: SIDES, rows: PROFILE.length - 1 },);/** * The way out from a link's centre line in its plane: each end bends through half a turn, and the * straight sides join them. */const outAt = (u: number) => { 'use gpu'; const column = std.round(u * (BEND * 2 + 2)); const a = std.min(std.select(column - 1, column, column <= BEND), BEND * 2) * (Math.PI / BEND); return d.vec3f(std.cos(a), std.sin(a), 0);};const centreAt = (u: number) => { 'use gpu'; const column = std.round(u * (BEND * 2 + 2)); const half = (LINK[0] - LINK[1]) / 2; const end = std.select(d.f32(-half), half, column <= BEND || column > BEND * 2 + 1); return d.vec3f(0, end, 0) + outAt(u) * (LINK[1] / 2 - LINK[2]);};const acrossAt = (u: number, v: number) => { 'use gpu'; const b = v * 2 * Math.PI; return outAt(u) * std.cos(b) - d.vec3f(0, 0, std.sin(b));};/** One chain link about its own middle: a stadium of wire, long along y and wide along x. */const link = meshes.parametric( { at: (u, v) => { 'use gpu'; return centreAt(u) + acrossAt(u, v) * LINK[2]; }, normalAt: acrossAt, }, { cols: BEND * 2 + 2, rows: WIRE },);/** * The links along one piece of chain: the first in `links`, how many, the particles they run * through, -1 standing for the point above the pivot, and the foot they end at, if any. */interface Run { first: number; count: number; particles: readonly number[]; foot: number;}/** * The censer's bronze body and its chains, swinging with `censer`. Each chain link is an instance, * laid along a smooth curve through the chain's particles and turned a quarter from the last. */export async function createModel(root: TgpuRoot, censer: Censer) { const mesh = await meshes.bakeAsync(root, meshes.concat(body, link)); const { chain } = censer; // The curve's control points, gathered for each run in turn const points = new Float64Array(3 * 24); const footAt: [number, number, number] = [0, 0, 0]; const put = (n: number, p: ArrayLike<number>) => { points[n * 3] = p[0]!; points[n * 3 + 1] = p[1]!; points[n * 3 + 2] = p[2]!; }; const runs: Run[] = []; let total = 0; const run = (length: number, particles: readonly number[], foot = -1) => { const count = Math.max(1, Math.round(length / PITCH)); runs.push({ first: total, count, particles, foot }); total += count; }; run(ABOVE, [-1, 0]); run(LENGTH - CROWN, chain.chain); chain.legs.forEach((leg, f) => { const [a, b] = [chain.positions[leg[0]!]!, chain.positions[leg[1]!]!]; run(Math.hypot(b[0] - a[0], b[1] - a[1], b[2] - a[2]) * leg.length, leg, f); }); const above = [PIVOT[0], PIVOT[1] + ABOVE, PIVOT[2]]; /** Puts a run's control points, and returns how many. */ const gather = ({ particles, foot }: Run) => { for (let n = 0; n < particles.length; n++) { put(n, particles[n]! < 0 ? above : chain.positions[particles[n]!]!); } if (foot < 0) return particles.length; put(particles.length, chain.footOf(foot, footAt)); return particles.length + 1; }; // Per link: its middle, the way it runs and the way across it, each a vec4 const links = root.createReadonly(d.arrayOf(d.vec4f, total * 3)).$name('censer links'); const laid = new Float32Array(total * 12); const at = [0, 0, 0]; const slope = [0, 0, 0]; /** * Writes into `at` and `slope` the point and slope of a Catmull–Rom curve through the first `n` * points, at `t` from 0 to `n - 1`. */ const curve = (n: number, t: number) => { const i = Math.min(n - 2, Math.floor(t)); const f = t - i; for (let k = 0; k < 3; k++) { const p1 = points[i * 3 + k]!; const p2 = points[(i + 1) * 3 + k]!; const p0 = i > 0 ? points[(i - 1) * 3 + k]! : 2 * p1 - p2; const p3 = i + 2 < n ? points[(i + 2) * 3 + k]! : 2 * p2 - p1; const a = -p0 + p2; const b = 2 * p0 - 5 * p1 + 4 * p2 - p3; const c = -p0 + 3 * p1 - 3 * p2 + p3; at[k] = 0.5 * (2 * p1 + a * f + b * f * f + c * f * f * f); slope[k] = 0.5 * (a + 2 * b * f + 3 * c * f * f); } }; /** Lays every link along its run, carrying the way across each from the last without twist. */ const lay = () => { for (let r = 0; r < runs.length; r++) { const { first, count } = runs[r]!; const n = gather(runs[r]!); let ux = 1; let uy = 0; let uz = 0; for (let k = 0; k < count; k++) { curve(n, ((k + 0.5) / count) * (n - 1)); const length = Math.hypot(slope[0]!, slope[1]!, slope[2]!) || 1; const wx = slope[0]! / length; const wy = slope[1]! / length; const wz = slope[2]! / length; // The last way across, less its part along this link, or any way across if none is left const along = ux * wx + uy * wy + uz * wz; ux -= wx * along; uy -= wy * along; uz -= wz * along; let across = Math.hypot(ux, uy, uz); if (across < 1e-4) { const useX = Math.abs(wx) < 0.9; ux = useX ? 1 - wx * wx : -wy * wx; uy = useX ? -wx * wy : 1 - wy * wy; uz = useX ? -wx * wz : -wy * wz; across = Math.hypot(ux, uy, uz); } ux /= across; uy /= across; uz /= across; // Every other link turns a quarter about the chain const odd = (first + k) % 2 === 1; const o = (first + k) * 12; laid[o] = at[0]!; laid[o + 1] = at[1]!; laid[o + 2] = at[2]!; laid[o + 4] = wx; laid[o + 5] = wy; laid[o + 6] = wz; laid[o + 8] = odd ? wy * uz - wz * uy : ux; laid[o + 9] = odd ? wz * ux - wx * uz : uy; laid[o + 10] = odd ? wx * uy - wy * ux : uz; } } links.buffer.write(laid.buffer); }; const drawing: passes.Drawing<typeof Corner, d.Infer<typeof Seen>> = { label: 'censer', Vertex: Corner, at: (vertexIndex: number, instanceIndex: number) => { 'use gpu'; const vertex = mesh.vertexAt(mesh.indexAt(vertexIndex)); const iron = vertexIndex >= body.indexCount; let world = censer.place(vertex.position); let normal = censer.turn(vertex.normal); if (iron) { const middle = links.$[instanceIndex * 3]!.xyz; const along = links.$[instanceIndex * 3 + 1]!.xyz; const across = links.$[instanceIndex * 3 + 2]!.xyz; const turned = d.mat3x3f(across, along, std.cross(across, along)); world = middle + turned * vertex.position; normal = turned * vertex.normal; } return Corner({ position: cameras.current.$.viewProjection * d.vec4f(world, 1), world, normal, local: vertex.position, iron: std.select(d.u32(0), d.u32(1), iron), }); }, surfaceOf: (corner) => { 'use gpu'; const normal = std.normalize(corner.normal); const shape = corner.local / SIZE; const grain = perlin3d.sample(shape * 9); const stain = perlin3d.sample(shape * 1.3); const patina = std.smoothstep(0.02, 0.43, stain + grain * 0.18); let color = std.mix(d.vec3f(0.68, 0.36, 0.12), d.vec3f(0.12, 0.2, 0.15), patina * 0.5); let roughness = 0.31 + patina * 0.18; let metallic = d.f32(0.9); let emission = d.vec3f(); // A polished rim, and above it the pierced lid const height = shape.y + LENGTH / SIZE; if (height > 0.1 && height < 0.22) { color = d.vec3f(0.85, 0.57, 0.22); roughness = 0.2; } if (height > 0.26 && height < 0.74) { // Lancet vents in the lid, lit from inside by the embers const around = std.fract(std.atan2(shape.z, shape.x) * (14 / (2 * Math.PI))); const taper = std.saturate(1 - std.abs(height - 0.48) / 0.2); const vent = 1 - std.smoothstep(0.12, 0.16, std.abs(around - 0.5) / std.max(taper, 0.01)); color = std.mix(color, d.vec3f(0.02, 0.012, 0.008), vent); emission = d.vec3f(1, 0.32, 0.07) * (vent * (1.4 + stain)); } if (corner.iron === 1) { color = d.vec3f(0.34, 0.32, 0.3); metallic = 0.85; roughness = 0.3; } const material = shading.Material({ color, roughness, metallic, occlusion: 1, emission }); const lit = shading.shade( materials.antialiased(material, normal), shading.pointAt(corner.world, normal), ); return Seen({ color: d.vec4f(lit, 1) }); }, color: passes.shaded, draw: (pass) => { pass.draw(body.indexCount); pass.draw(link.indexCount, total, body.indexCount); }, }; return { drawing, /** Lays the chain's links where its particles are now. */ update: lay, destroy() { mesh.destroy(); links.buffer.destroy(); }, };}smoke.ts
import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';import { dispatch, groups, inside, kernel, volume, WORKGROUP } from './cells.ts';import { ALL_SOLID, createPressure, STEPS } from './pressure.ts';type Vec3 = readonly [number, number, number];/** How fast the sponge at the box's faces takes smoke and motion away, per second, at the face. */const SPONGE_RATE = 6;/** How the smoke moves. Accelerations are in metres per second squared, decays per second. */const SmokeSettings = d.struct({ /** The step's length in seconds. */ dt: d.f32, /** Upward acceleration per unit of heat. */ buoyancy: d.f32, /** Downward acceleration per unit of smoke. */ weight: d.f32, /** How strongly vorticity confinement feeds back the swirls the grid smooths away. */ vorticity: d.f32, smokeDecay: d.f32, heatDecay: d.f32, drag: d.f32,});/** Smoke and heat fainter than this are dropped as they move: too faint to show or lift the air. */const FAINT = 1e-4;/** The most steps a path is traced back in, one per cell it crosses. */const MOST_TRACE = 8;/** * How many times finer than the air's cells the smoke and heat are carried, along each axis, so * they keep edges sharper than the motion that carries them. Each workgroup covers whole cells of * the air's grid, which it averages its smoke into. */const DETAIL = 2;/** A thread's place in its workgroup's shared array. */const slotOf = (local: d.v3u) => { 'use gpu'; return local.x + WORKGROUP[0] * (local.y + WORKGROUP[1] * local.z);};/** A cell of a field, or the nearest one inside it. */const load = (texture: d.texture3d<d.F32>, cell: d.v3i) => { 'use gpu'; const edge = d.vec3i(std.textureDimensions(texture)) - 1; return std.textureLoad(texture, std.clamp(cell, d.vec3i(), edge), 0);};/** Whether the neighbour a bit of the mask names is solid. */const isSolid = (mask: number, bit: number) => { 'use gpu'; return (mask & d.u32(bit)) !== 0;};interface SmokeOptions { /** Cells along x, y and z, each a multiple of 4. */ readonly cells: Vec3; /** The box's lowest corner, in metres. */ readonly origin: Vec3; /** A cell's edge in metres. */ readonly cell: number; /** * What fills a point, a `'use gpu'` function of its position: the velocity of the solid there in * `xyz`, and 1 in `w` inside a solid, 0 in the air. */ readonly solid: (position: d.v3f) => d.v4f; /** Smoke in `x` and heat in `y` added per second at a point, a `'use gpu'` function of it. */ readonly source: (position: d.v3f) => d.v2f; /** The acceleration a point's air takes, in metres per second squared, such as a jet's. */ readonly force: (position: d.v3f) => d.v3f; /** * The velocity a point's air is carried toward in `xyz`, and how fast in `w`, per second, such as * air a moving solid drags along beside it. */ readonly carry: (position: d.v3f) => d.v4f; /** * Cells inside each face over which a sponge takes smoke, heat and motion away, more strongly * toward the face, so nothing reaches the walls. */ readonly sponge: number; readonly settings: d.InferInput<typeof SmokeSettings>;}/** * A box of smoke on a grid of cells: the air's velocity, advected by MacCormack's method, lifted by * heat, kept swirling by vorticity confinement and made incompressible by a multigrid pressure * solve, carrying smoke and heat on a grid twice as fine. Solids slip freely and push the air * with their velocity. The box lets smoke out of every face, through a sponge that takes it away * before it reaches one. */export function createSmoke(root: TgpuRoot, options: SmokeOptions) { const [nx, ny, nz] = options.cells; const { cell: h, solid, source, force, carry, sponge } = options; const [fx, fy, fz] = [nx * DETAIL, ny * DETAIL, nz * DETAIL]; const origin = d.vec3f(...options.origin); const size = d.vec3f(nx, ny, nz); const box = d.vec3i(nx, ny, nz); const last = d.vec3i(nx - 1, ny - 1, nz - 1); const settings = root.createUniform(SmokeSettings, options.settings).$name('smoke settings'); // Smoke and heat are a pair, which `texture-formats-tier1` stores in half the bytes const pair = root.device.features.has('texture-formats-tier1') ? 'rg16float' : 'rgba16float'; const field = ( name: string, format: 'rgba16float' | 'rg16float' = 'rgba16float', cells: Vec3 = options.cells, ) => volume(root, cells, format, `smoke ${name}`); // The velocity holds each step's result and the advected one holds it before the projection, with // 1 in `w` in solid cells const velocity = field('velocity'); const advected = field('advected'); const solids = field('solids'); const curls = field('curls'); const averaged = field('averaged', pair); const solver = createPressure(root, options.cells); const written = d.textureStorage3d('r32float', 'write-only'); const spread = solver.rhs.createView(written); const masks = { write: solver.mask.createView(written), read: solver.mask.createView(d.texture3d(d.f32), { sampleType: 'unfilterable-float' }), }; const pressure = solver.pressure.createView(d.textureStorage3d('r32float', 'read-write')); const smooth = root.createSampler({ magFilter: 'linear', minFilter: 'linear' }); // Smoke in x and heat in y, on the finer grid: each step carries one of the pair into the other const matter = [ field('matter a', pair, [fx, fy, fz]), field('matter b', pair, [fx, fy, fz]), ] as const; const matterLayout = tgpu.bindGroupLayout({ matter: { texture: d.texture3d(d.f32) } }); const stepLayout = tgpu.bindGroupLayout({ from: { texture: d.texture3d(d.f32) }, to: { storageTexture: d.textureStorage3d(pair, 'write-only') }, }); const reading = matter.map(({ texture }) => root.createBindGroup(matterLayout, { matter: texture }), ); const stepping = [0, 1].map((i) => root.createBindGroup(stepLayout, { from: matter[i]!.texture, to: matter[1 - i]!.texture }), ); let current = 0; /** Samples a field at a point in cells, with cell centres at half steps. */ const sampled = (texture: d.texture3d<d.F32>, at: d.v3f) => { 'use gpu'; return std.textureSampleLevel(texture, smooth.$, at / size, 0); }; const positionOf = (cell: d.v3i) => { 'use gpu'; return origin + (d.vec3f(cell) + 0.5) * h; }; /** * Where the air at a point in cells was `dt` seconds before, or will be with a negative `dt`, * traced in steps of at most a cell so fast air follows its curve rather than cutting across it. */ const trace = (at: d.v3f, dt: number) => { 'use gpu'; const scale = dt / h; const first = sampled(velocity.read.$, at).xyz * scale; const steps = std.clamp(std.ceil(std.length(first)), 1, MOST_TRACE); let point = at - first / steps; for (let i = d.i32(1); i < d.i32(steps); i++) { point = point - (sampled(velocity.read.$, point).xyz * scale) / steps; } return point; }; /** * Advects a field from where the air came from, corrects it by tracing that forwards and back * again, and clamps it to the cells the path started between, so the correction adds no new * extremes. `predicted` is the field where the air came from. Points are in the air's cells; the * field's cells are `scale` times finer. */ const macCormackOf = (scale: number) => { return ( texture: d.texture3d<d.F32>, predicted: d.v4f, at: d.v3f, back: d.v3f, restored: d.v3f, ) => { 'use gpu'; const corrected = predicted + (sampled(texture, at) - sampled(texture, restored)) * 0.5; const base = d.vec3i(std.floor(back * scale - 0.5)); let low = load(texture, base); let high = d.vec4f(low); for (const z of tgpu.unroll([0, 1])) { for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const corner = load(texture, base + d.vec3i(x, y, z)); low = std.min(low, corner); high = std.max(high, corner); } } } return std.clamp(corrected, low, high); }; }; const macCormack = macCormackOf(1); const macCormackFine = macCormackOf(DETAIL); /** What the sponge leaves of anything at a point in cells over a step. */ const spongeAt = (at: d.v3f) => { 'use gpu'; const wall = std.min( std.min(at.x, at.y), std.min(at.z, std.min(std.min(size.x - at.x, size.y - at.y), size.z - at.z)), ); const deep = 1 - std.saturate(wall / sponge); return std.exp(-SPONGE_RATE * settings.$.dt * deep * deep); }; const findSolids = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; std.textureStore(solids.write.$, id, solid(positionOf(d.vec3i(id)))); }); const findCurl = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; const cell = d.vec3i(id); const dx = load(velocity.read.$, cell + d.vec3i(1, 0, 0)) - load(velocity.read.$, cell - d.vec3i(1, 0, 0)); const dy = load(velocity.read.$, cell + d.vec3i(0, 1, 0)) - load(velocity.read.$, cell - d.vec3i(0, 1, 0)); const dz = load(velocity.read.$, cell + d.vec3i(0, 0, 1)) - load(velocity.read.$, cell - d.vec3i(0, 0, 1)); const curl = d.vec3f(dy.z - dz.y, dz.x - dx.z, dx.y - dy.x) * (0.5 / h); std.textureStore(curls.write.$, id, d.vec4f(curl, std.length(curl))); }); const advect = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; const cell = d.vec3i(id); const here = std.textureLoad(solids.read.$, cell, 0); if (here.w > 0.5) { std.textureStore(advected.write.$, id, here); return; } const dt = settings.$.dt; const at = d.vec3f(cell) + 0.5; const back = trace(at, dt); const restored = trace(trace(at, -dt), dt); let moved = macCormack(velocity.read.$, sampled(velocity.read.$, back), at, back, restored).xyz; const held = sampled(matterLayout.$.matter, at).xy; // Vorticity confinement turns the air around the strongest swirls nearby const swirl = std.textureLoad(curls.read.$, cell, 0).xyz; const rise = d.vec3f( load(curls.read.$, cell + d.vec3i(1, 0, 0)).w - load(curls.read.$, cell - d.vec3i(1, 0, 0)).w, load(curls.read.$, cell + d.vec3i(0, 1, 0)).w - load(curls.read.$, cell - d.vec3i(0, 1, 0)).w, load(curls.read.$, cell + d.vec3i(0, 0, 1)).w - load(curls.read.$, cell - d.vec3i(0, 0, 1)).w, ); const toward = rise / (std.length(rise) + 1e-5); const confinement = std.cross(toward, swirl) * (settings.$.vorticity * h); const position = positionOf(cell); const lift = d.vec3f(0, settings.$.buoyancy * held.y - settings.$.weight * held.x, 0); const pushed = confinement + lift + force(position); moved = (moved + pushed * dt) * std.exp(-settings.$.drag * dt); const dragged = carry(position); moved = std.mix(moved, dragged.xyz, 1 - std.exp(-dragged.w * dt)); std.textureStore(advected.write.$, id, d.vec4f(moved * spongeAt(at), 0)); }); const carried = tgpu.workgroupVar(d.arrayOf(d.vec2f, WORKGROUP[0] * WORKGROUP[1] * WORKGROUP[2])); /** * Smoke and heat, carried on the finer grid by the air's velocity at the start of the step, and * averaged over each of the air's cells. */ const advectMatter = tgpu.computeFn({ workgroupSize: [...WORKGROUP], in: { id: d.builtin.globalInvocationId, local: d.builtin.localInvocationId }, })(({ id, local }) => { 'use gpu'; const inBox = id.x < fx && id.y < fy && id.z < fz; const at = (d.vec3f(id) + 0.5) / DETAIL; let stored = d.vec4f(); if (inBox && std.textureLoad(solids.read.$, d.vec3i(at), 0).w <= 0.5) { const dt = settings.$.dt; const back = trace(at, dt); const predicted = sampled(stepLayout.$.from, back); const added = source(origin + at * h) * dt; // Where the air brings in only traces, they are dropped rather than corrected if (std.max(predicted.x, predicted.y) < FAINT) { stored = d.vec4f(added * spongeAt(at), 0, 0); } else { const restored = trace(trace(at, -dt), dt); let moved = macCormackFine(stepLayout.$.from, predicted, at, back, restored).xy; moved = moved + added; moved = moved * d.vec2f(std.exp(-settings.$.smokeDecay * dt), std.exp(-settings.$.heatDecay * dt)); stored = d.vec4f(std.max(moved * spongeAt(at), d.vec2f()), 0, 0); } } if (inBox) std.textureStore(stepLayout.$.to, id, stored); // As the texture holds it, in half precision carried.$[slotOf(local)] = std.quantizeToF16(stored.xy); std.workgroupBarrier(); const corner = d.vec3u( std.intdiv(id.x, DETAIL), std.intdiv(id.y, DETAIL), std.intdiv(id.z, DETAIL), ); const first = local.x % DETAIL === 0 && local.y % DETAIL === 0 && local.z % DETAIL === 0; if (first && corner.x < nx && corner.y < ny && corner.z < nz) { let sum = d.vec2f(); for (const z of tgpu.unroll([0, 1])) { for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) sum += carried.$[slotOf(local + d.vec3u(x, y, z))]!; } } std.textureStore(averaged.write.$, corner, d.vec4f(sum / DETAIL ** 3, 0, 0)); } }); /** A neighbour's velocity, and whether it is solid. Outside the box is open. */ const neighbour = (cell: d.v3i, offset: d.v3i, own: d.v3f) => { 'use gpu'; const next = cell + offset; if (!inside(next, box)) return d.vec4f(own, 0); return std.textureLoad(advected.read.$, next, 0); }; const divergence = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; const cell = d.vec3i(id); const here = std.textureLoad(advected.read.$, cell, 0); if (here.w > 0.5) { std.textureStore(spread.$, id, d.vec4f()); std.textureStore(masks.write.$, id, d.vec4f(ALL_SOLID, 0, 0, 0)); return; } const own = here.xyz; const left = neighbour(cell, d.vec3i(-1, 0, 0), own); const right = neighbour(cell, d.vec3i(1, 0, 0), own); const below = neighbour(cell, d.vec3i(0, -1, 0), own); const above = neighbour(cell, d.vec3i(0, 1, 0), own); const back = neighbour(cell, d.vec3i(0, 0, -1), own); const front = neighbour(cell, d.vec3i(0, 0, 1), own); const outflow = (right.x - left.x + above.y - below.y + front.z - back.z) * 0.5; const mask = left.w + right.w * 2 + below.w * 4 + above.w * 8 + back.w * 16 + front.w * 32; std.textureStore(spread.$, id, d.vec4f(outflow, 0, 0, 0)); std.textureStore(masks.write.$, id, d.vec4f(mask, 0, 0, 0)); }); /** Pressure at a neighbour, 0 outside the box, or the cell's own where the neighbour is solid. */ const pressureNear = (cell: d.v3i, offset: d.v3i, solidNeighbour: boolean, own: number) => { 'use gpu'; if (solidNeighbour) return own; const next = cell + offset; if (!inside(next, box)) return d.f32(0); return std.textureLoad(pressure.$, next).x; }; /** The velocity of the solid at a neighbour. */ const wall = (cell: d.v3i, offset: d.v3i) => { 'use gpu'; return std.textureLoad(solids.read.$, std.clamp(cell + offset, d.vec3i(), last), 0).xyz; }; const project = kernel(({ id }) => { 'use gpu'; if (id.x >= nx || id.y >= ny || id.z >= nz) return; const cell = d.vec3i(id); const mask = d.u32(std.textureLoad(masks.read.$, cell, 0).x); if (mask === ALL_SOLID) { std.textureStore( velocity.write.$, id, d.vec4f(std.textureLoad(advected.read.$, cell, 0).xyz, 0), ); return; } const own = std.textureLoad(pressure.$, cell).x; const moved = std.textureLoad(advected.read.$, cell, 0).xyz; // Air beside a solid moves with it across the face between them and slips along it for (const axis of tgpu.unroll([0, 1, 2])) { const before = isSolid(mask, 1 << (axis * 2)); const after = isSolid(mask, 2 << (axis * 2)); const down = STEPS[axis * 2]!; const up = STEPS[axis * 2 + 1]!; const near = pressureNear(cell, up, after, own) - pressureNear(cell, down, before, own); moved[axis] = moved[axis]! - near * 0.5; if (before) moved[axis] = wall(cell, down)[axis]!; if (after) moved[axis] = wall(cell, up)[axis]!; } std.textureStore(velocity.write.$, id, d.vec4f(moved, 0)); }); const pipelines = { solids: root.createComputePipeline({ compute: findSolids }), curl: root.createComputePipeline({ compute: findCurl }), advect: root.createComputePipeline({ compute: advect }), matter: root.createComputePipeline({ compute: advectMatter }), divergence: root.createComputePipeline({ compute: divergence }), project: root.createComputePipeline({ compute: project }), }; const cellGroups = groups(options.cells); const fineGroups = groups([fx, fy, fz]); const advecting = reading.map((group) => pipelines.advect.with(group)); const carrying = stepping.map((group) => pipelines.matter.with(group)); const low: Vec3 = options.origin; const high: Vec3 = [low[0] + nx * h, low[1] + ny * h, low[2] + nz * h]; const extent = d.vec3f(nx * h, ny * h, nz * h); return { /** Cells along x, y and z, and a cell's edge in metres. */ cells: options.cells, cell: h, /** The box's lowest and highest corners, in metres. */ low, high, /** Binds the smoke as it stands for pipelines that read `smokeAt`: changes with each step. */ get matter() { return reading[current]!; }, /** The smoke and heat at a point in metres, 0 outside the box, read where `matter` is bound. */ smokeAt: (position: d.v3f) => { 'use gpu'; const at = (position - origin) / extent; if (std.any(std.lt(at, d.vec3f())) || std.any(std.gt(at, d.vec3f(1)))) return d.vec2f(); return std.textureSampleLevel(matterLayout.$.matter, smooth.$, at, 0).xy; }, /** The smoke and heat averaged over each of the air's cells, as the last step left them. */ averaged: averaged.read, initAsync: () => Promise.all([ ...Object.values(pipelines).map((pipeline) => pipeline.initAsync()), solver.initAsync(), ]), /** * Records the first half of a step into a compute pass: finds the solids, then moves the air * and the smoke and measures how much the moved air spreads. `matter` changes here. */ advect(pass: TgpuComputePass) { pass.setPipeline(pipelines.solids); dispatch(pass, cellGroups); pass.setPipeline(pipelines.curl); dispatch(pass, cellGroups); pass.setPipeline(advecting[current]!); dispatch(pass, cellGroups); pass.setPipeline(carrying[current]!); dispatch(pass, fineGroups); current = 1 - current; pass.setPipeline(pipelines.divergence); dispatch(pass, cellGroups); }, /** * Records the second half of a step, after `advect`: solves the pressure and makes the air * incompressible. The two halves cost about the same, so a program may record them in * different frames. */ project(pass: TgpuComputePass) { solver.solve(pass); pass.setPipeline(pipelines.project); dispatch(pass, cellGroups); }, destroy() { for (const each of [velocity, advected, solids, curls, averaged, ...matter]) each.texture.destroy(); solver.destroy(); }, };}export type Smoke = ReturnType<typeof createSmoke>;chain.ts
type Vec3 = [number, number, number];type Quat = [number, number, number, number];const GRAVITY = 9.81;/** The step the chain is solved in, in seconds: small steps let light links hold a heavy body. */const SUBSTEP = 1 / 2400;/** Air slowing the links, the body and the body's turning, per second. */const LINK_AIR = 0.4;const BODY_AIR = 0.05;const TURN_AIR = 0.6;/** A hand pulls like a critically damped spring of this frequency, in hertz, at most this hard. */const HAND = 2.6;const HAND_DAMPING = 0.8;const HAND_MOST = 80;/** The hand also holds the body's turn, loosely, like a wrist: hertz and damping ratio. */const WRIST = 1.2;const WRIST_DAMPING = 0.7;/** * After a hand lets go, the idle swing waits this many seconds and then eases in over as many * again, pumping or braking the swing toward its amplitude at this rate per second. */const IDLE_WAIT = 3;const IDLE_EASE = 3;const IDLE_RATE = 0.9;/** Where the censer's chain hangs from, in metres. */export const PIVOT: Vec3 = [0, 3.8, 0];/** From the pivot to the middle of the censer's body. */export const LENGTH = 2.6;/** The body's size: its model is drawn in units of this many metres. */export const SIZE = 0.4;/** * The censer's outline from its finial to its foot: radius and height from the body's middle, in * units of `SIZE`. */export const PROFILE = [ 0, 1, 0.05, 0.98, 0.075, 0.91, 0.04, 0.84, 0.1, 0.8, 0.2, 0.74, 0.3, 0.64, 0.38, 0.52, 0.44, 0.38, 0.47, 0.26, 0.5, 0.22, 0.55, 0.2, 0.55, 0.12, 0.51, 0.1, 0.53, 0, 0.55, -0.12, 0.52, -0.28, 0.44, -0.42, 0.3, -0.53, 0.16, -0.58, 0.11, -0.64, 0.1, -0.72, 0.2, -0.78, 0.25, -0.84, 0.24, -0.88, 0, -0.89,].flatMap((v, i, all) => (i % 2 ? [] : [[v, all[i + 1]!] as const]));/** Where the legs meet above the body's middle, and where they end on its rim, in metres. */export const CROWN = 1.85 * SIZE;const FEET = [0, 1, 2].map((i) => { const a = (i / 3) * Math.PI * 2 + 0.5; return [Math.cos(a) * 0.53 * SIZE, 0.16 * SIZE, Math.sin(a) * 0.53 * SIZE] as Vec3;});/** * A bronze body of a few kilograms on a light iron chain: its mass, its moments of inertia about * its axis and across it, and the chain's mass per metre. */const MASS = 3.5;const INERTIA: Vec3 = [0.12, 0.07, 0.12];const DENSITY = 0.3;/** Links in the chain to the crown and in each leg. */const LINKS = 16;const LEG_LINKS = 4;/** The amplitude the idle swing settles to, in radians. */const IDLE = 0.5;// The chain steps thousands of times a second, so its vector maths writes into vectors made once// rather than returning new ones; `out` may be an inputconst dot = (a: Vec3, b: Vec3) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];/** Writes `a × b` into `out`. */const cross = (out: Vec3, a: Readonly<Vec3>, b: Readonly<Vec3>) => { const x = a[1] * b[2] - a[2] * b[1]; const y = a[2] * b[0] - a[0] * b[2]; const z = a[0] * b[1] - a[1] * b[0]; out[0] = x; out[1] = y; out[2] = z; return out;};/** Writes `v` turned by the unit quaternion (x, y, z, w) into `out`. */const turnBy = (out: Vec3, x: number, y: number, z: number, w: number, v: Readonly<Vec3>) => { const v0 = v[0]; const v1 = v[1]; const v2 = v[2]; const t0 = (y * v2 - z * v1) * 2; const t1 = (z * v0 - x * v2) * 2; const t2 = (x * v1 - y * v0) * 2; out[0] = v0 + w * t0 + (y * t2 - z * t1); out[1] = v1 + w * t1 + (z * t0 - x * t2); out[2] = v2 + w * t2 + (x * t1 - y * t0); return out;};/** Writes `v` turned by the unit quaternion `q`, or by its inverse, into `out`. */const rotate = (out: Vec3, q: Readonly<Quat>, v: Readonly<Vec3>) => turnBy(out, q[0], q[1], q[2], q[3], v);const unrotate = (out: Vec3, q: Readonly<Quat>, v: Readonly<Vec3>) => turnBy(out, -q[0], -q[1], -q[2], q[3], v);/** Turns `q` on by the small rotation vector `w`, and keeps it unit. */const spin = (q: Quat, w: Readonly<Vec3>, scale: number) => { const x = w[0] * scale * 0.5; const y = w[1] * scale * 0.5; const z = w[2] * scale * 0.5; const a = q[0]; const b = q[1]; const c = q[2]; const s = q[3]; q[0] += x * s + y * c - z * b; q[1] += y * s + z * a - x * c; q[2] += z * s + x * b - y * a; q[3] -= x * a + y * b + z * c; const length = Math.hypot(q[0], q[1], q[2], q[3]); for (let i = 0; i < 4; i++) q[i]! /= length;};/** Writes the small rotation vector that turns `from` to `to`, the shorter way, into `out`. */const turnFrom = (out: Vec3, from: Readonly<Quat>, to: Readonly<Quat>) => { const a = from[0]; const b = from[1]; const c = from[2]; const s = from[3]; const x = to[0]; const y = to[1]; const z = to[2]; const w = to[3]; const sign = w * s + x * a + y * b + z * c < 0 ? -2 : 2; out[0] = (w * -a + x * s + y * -c - z * -b) * sign; out[1] = (w * -b + y * s + z * -a - x * -c) * sign; out[2] = (w * -c + z * s + x * -b - y * -a) * sign; return out;};/** * A heavy body hanging by three legs of chain from a crown, and the crown by a chain from a fixed * pivot, stepped by extended position-based dynamics in small steps. Links pull but never push, so * a chain lifted toward its pivot falls slack and drapes until it snaps taut again. The body is * rigid: it tilts and turns with its own inertia as the legs pull on its feet. A hand pulls it as a * damped spring toward a point, so it keeps its weight and momentum and flies on when let go. Left * alone for a while, a gentle drive keeps it swinging. It starts swung `start` radians out in the * x–y plane. Units are metres, kilograms and seconds. */export function createChain(start = 0.4) { const length = LENGTH - CROWN; const outline = PROFILE.map(([r, y]) => [r * SIZE, y * SIZE] as const); const top = outline[0]![1]; const bottom = outline[outline.length - 1]![1]; /** The body's radius at a height from its middle. */ const widthAt = (y: number) => { for (let i = 1; i < outline.length; i++) { const r1 = outline[i]![0]; const y1 = outline[i]![1]; if (y >= y1) { const r0 = outline[i - 1]![0]; const y0 = outline[i - 1]![1]; return r1 + ((r0 - r1) * (y - y1)) / Math.max(1e-6, y0 - y1); } } return 0; }; const turned = (v: Vec3): Vec3 => [ v[0] * Math.cos(start) - v[1] * Math.sin(start), v[0] * Math.sin(start) + v[1] * Math.cos(start), v[2], ]; const hung = (v: Vec3): Vec3 => { const t = turned(v); return [PIVOT[0] + t[0], PIVOT[1] + t[1], PIVOT[2] + t[2]]; }; // The body: its middle, which is its centre of mass, its turn and their rates const below = -length - CROWN; const middle = hung([0, below, 0]); const halfStart = start / 2; const turn: Quat = [0, 0, Math.sin(halfStart), Math.cos(halfStart)]; const velocity: Vec3 = [0, 0, 0]; const spinning: Vec3 = [0, 0, 0]; const lastMiddle: Vec3 = [...middle]; const lastTurn: Quat = [...turn]; // Particles: the pivot, the chain down to the crown, then each leg's links between crown and foot const positions: Vec3[] = []; const inverse: number[] = []; const linkMass = (DENSITY * length) / LINKS; for (let i = 0; i <= LINKS; i++) { positions.push(hung([0, (-length * i) / LINKS, 0])); inverse.push(i === 0 ? 0 : 1 / (i === LINKS ? linkMass * 2 : linkMass)); } const crownAt = LINKS; type Rod = { a: number; b: number; foot: number; rest: number }; const rods: Rod[] = []; for (let i = 0; i < LINKS; i++) rods.push({ a: i, b: i + 1, foot: -1, rest: length / LINKS }); const legs: number[][] = []; FEET.forEach((foot, f) => { const span = Math.hypot(foot[0], foot[1] - CROWN, foot[2]); const legMass = (DENSITY * span) / LEG_LINKS; const leg = [crownAt]; for (let i = 1; i < LEG_LINKS; i++) { const t = i / LEG_LINKS; positions.push(hung([foot[0] * t, below + CROWN + (foot[1] - CROWN) * t, foot[2] * t])); inverse.push(1 / legMass); leg.push(positions.length - 1); } for (let i = 0; i + 1 < leg.length; i++) { rods.push({ a: leg[i]!, b: leg[i + 1]!, foot: -1, rest: span / LEG_LINKS }); } rods.push({ a: leg[leg.length - 1]!, b: -1, foot: f, rest: span / LEG_LINKS }); legs.push(leg); }); const count = positions.length; const previous = positions.map((p) => [...p] as Vec3); const speeds = positions.map(() => [0, 0, 0] as Vec3); const hand = { held: false, grip: [0, 0, 0] as Vec3, target: [0, 0, 0] as Vec3, turn: [0, 0, 0, 1] as Quat, // Where the target was at the last advance, where it is within this one and how fast it moves from: [0, 0, 0] as Vec3, now: [0, 0, 0] as Vec3, moving: [0, 0, 0] as Vec3, }; let free = IDLE_WAIT + IDLE_EASE; let left = 0; const reach = length + CROWN; /** A point of the body given from its middle in its own frame, in the world. */ const pointOf = (local: Vec3): Vec3 => { const r = rotate([0, 0, 0], turn, local); return [middle[0] + r[0], middle[1] + r[1], middle[2] + r[2]]; }; const given: Vec3 = [0, 0, 0]; /** How freely the body gives along `n` at offset `r` from its middle. */ const giveOf = (r: Vec3, n: Vec3) => { const k = unrotate(given, turn, cross(given, r, n)); return 1 / MASS + k[0] ** 2 / INERTIA[0] + k[1] ** 2 / INERTIA[1] + k[2] ** 2 / INERTIA[2]; }; const shoved: Vec3 = [0, 0, 0]; /** Pushes the body by `impulse` at offset `r` from its middle, in position. */ const shove = (r: Vec3, impulse: Vec3) => { for (let i = 0; i < 3; i++) middle[i]! += impulse[i]! / MASS; const local = unrotate(shoved, turn, cross(shoved, r, impulse)); for (let i = 0; i < 3; i++) local[i] = local[i]! / INERTIA[i]!; spin(turn, rotate(local, turn, local), 1); }; /** The angle from straight down the body's energy would carry it to. */ const amplitudeOf = () => { const height = middle[1] - (PIVOT[1] - reach); const energy = dot(velocity, velocity) / 2 + GRAVITY * height; return Math.acos(Math.max(-1, Math.min(1, 1 - energy / (GRAVITY * reach)))); }; const force: Vec3 = [0, 0, 0]; const torque: Vec3 = [0, 0, 0]; // Scratch vectors for one substep, each for one use const gripped: Vec3 = [0, 0, 0]; const moving: Vec3 = [0, 0, 0]; const pull: Vec3 = [0, 0, 0]; const twist: Vec3 = [0, 0, 0]; const error: Vec3 = [0, 0, 0]; const turning: Vec3 = [0, 0, 0]; const wrist: Vec3 = [0, 0, 0]; const local: Vec3 = [0, 0, 0]; const pushed: Vec3 = [0, 0, 0]; const delta: Vec3 = [0, 0, 0]; const normal: Vec3 = [0, 0, 0]; const lever: Vec3 = [0, 0, 0]; const impulse: Vec3 = [0, 0, 0]; const inBody: Vec3 = [0, 0, 0]; const outward: Vec3 = [0, 0, 0]; const substep = (h: number) => { // Forces on the body: its weight, a hand or the idle drive force[0] = 0; force[1] = -GRAVITY * MASS; force[2] = 0; torque[0] = torque[1] = torque[2] = 0; if (hand.held) { const r = rotate(gripped, turn, hand.grip); cross(moving, spinning, r); for (let i = 0; i < 3; i++) moving[i] = moving[i]! + velocity[i]! - hand.moving[i]!; const w = 2 * Math.PI * HAND; for (let i = 0; i < 3; i++) { const v = hand.now[i]! - (middle[i]! + r[i]!); pull[i] = MASS * (w * w * v - 2 * HAND_DAMPING * w * moving[i]!); } const strength = Math.hypot(pull[0], pull[1], pull[2]); const most = MASS * HAND_MOST; if (strength > most) for (let i = 0; i < 3; i++) pull[i]! *= most / strength; for (let i = 0; i < 3; i++) force[i]! += pull[i]!; cross(twist, r, pull); // The wrist turns the body back toward how it was taken hold of, about its own axes unrotate(error, turn, turnFrom(error, turn, hand.turn)); unrotate(turning, turn, spinning); const v = 2 * Math.PI * WRIST; for (let k = 0; k < 3; k++) { wrist[k] = INERTIA[k]! * (v * v * error[k]! - 2 * WRIST_DAMPING * v * turning[k]!); } rotate(wrist, turn, wrist); for (let i = 0; i < 3; i++) torque[i]! += twist[i]! + wrist[i]!; } else { const share = Math.min(1, Math.max(0, (free - IDLE_WAIT) / IDLE_EASE)); if (share > 0) { // Pump or brake along the way the body moves, from straight down along x when it is still const speed = Math.hypot(velocity[0], velocity[1], velocity[2]); const still = speed <= 0.05; const push = share * IDLE_RATE * (IDLE - amplitudeOf()) * Math.sqrt(GRAVITY * reach) * MASS; force[0] += (still ? (middle[0] >= PIVOT[0] ? -1 : 1) : velocity[0] / speed) * push; force[1] += (still ? 0 : velocity[1] / speed) * push; force[2] += (still ? 0 : velocity[2] / speed) * push; } } // Move everything freely for (let i = 0; i < count; i++) { if (inverse[i] === 0) continue; const p = positions[i]!; const v = speeds[i]!; v[1] -= GRAVITY * h; previous[i]![0] = p[0]; previous[i]![1] = p[1]; previous[i]![2] = p[2]; for (let k = 0; k < 3; k++) p[k]! += v[k]! * h; } for (let k = 0; k < 3; k++) { velocity[k]! += (force[k]! / MASS) * h; lastMiddle[k] = middle[k]!; middle[k]! += velocity[k]! * h; } // Torque turns the body about its own axes, which hold its inertia unrotate(local, turn, spinning); const h0 = local[0] * INERTIA[0]; const h1 = local[1] * INERTIA[1]; const h2 = local[2] * INERTIA[2]; unrotate(pushed, turn, torque); pushed[0] = (pushed[0] - (local[1] * h2 - local[2] * h1)) / INERTIA[0]; pushed[1] = (pushed[1] - (local[2] * h0 - local[0] * h2)) / INERTIA[1]; pushed[2] = (pushed[2] - (local[0] * h1 - local[1] * h0)) / INERTIA[2]; const change = rotate(pushed, turn, pushed); for (let k = 0; k < 3; k++) spinning[k]! += change[k]! * h; for (let k = 0; k < 4; k++) lastTurn[k] = turn[k]!; spin(turn, spinning, h); // Each link pulls its ends back to its length, and never pushes for (let j = 0; j < rods.length; j++) { const rod = rods[j]!; const a = positions[rod.a]!; const wa = inverse[rod.a]!; if (rod.foot < 0) { const b = positions[rod.b]!; const wb = inverse[rod.b]!; for (let k = 0; k < 3; k++) delta[k] = a[k]! - b[k]!; const apart = Math.hypot(delta[0], delta[1], delta[2]); const stretch = apart - rod.rest; if (stretch <= 0 || wa + wb === 0) continue; const move = stretch / (wa + wb) / apart; for (let k = 0; k < 3; k++) { a[k]! -= delta[k]! * move * wa; b[k]! += delta[k]! * move * wb; } } else { const r = rotate(lever, turn, FEET[rod.foot]!); for (let k = 0; k < 3; k++) delta[k] = a[k]! - (middle[k]! + r[k]!); const apart = Math.hypot(delta[0], delta[1], delta[2]); const stretch = apart - rod.rest; if (stretch <= 0) continue; for (let k = 0; k < 3; k++) normal[k] = delta[k]! / apart; const lambda = stretch / (wa + giveOf(r, normal)); for (let k = 0; k < 3; k++) a[k]! -= normal[k]! * lambda * wa; for (let k = 0; k < 3; k++) impulse[k] = normal[k]! * lambda; shove(r, impulse); } } // Links rest on the body rather than passing through it for (let i = 1; i < count; i++) { const p = positions[i]!; for (let k = 0; k < 3; k++) inBody[k] = p[k]! - middle[k]!; unrotate(inBody, turn, inBody); if (inBody[1] > top || inBody[1] < bottom) continue; const apart = Math.hypot(inBody[0], inBody[2]); const radius = widthAt(inBody[1]); if (apart >= radius) continue; // Out the shortest way: through the side, the foot or the finial let depth = radius - apart; outward[0] = apart > 1e-6 ? inBody[0] / apart : 1; outward[1] = 0; outward[2] = apart > 1e-6 ? inBody[2] / apart : 0; if (inBody[1] - bottom < depth) { depth = inBody[1] - bottom; outward[0] = 0; outward[1] = -1; outward[2] = 0; } if (top - inBody[1] < depth) { depth = top - inBody[1]; outward[0] = 0; outward[1] = 1; outward[2] = 0; } const n = rotate(normal, turn, outward); for (let k = 0; k < 3; k++) lever[k] = inBody[k]! + outward[k]! * depth; const r = rotate(lever, turn, lever); const lambda = depth / (inverse[i]! + giveOf(r, n)); for (let k = 0; k < 3; k++) p[k]! += n[k]! * lambda * inverse[i]!; for (let k = 0; k < 3; k++) impulse[k] = n[k]! * -lambda; shove(r, impulse); } // Velocities from the moves, slowed by the air const linkKeep = Math.exp(-LINK_AIR * h); for (let i = 0; i < count; i++) { if (inverse[i] === 0) continue; for (let k = 0; k < 3; k++) { speeds[i]![k] = ((positions[i]![k]! - previous[i]![k]!) / h) * linkKeep; } } const bodyKeep = Math.exp(-BODY_AIR * h); for (let k = 0; k < 3; k++) velocity[k] = ((middle[k]! - lastMiddle[k]!) / h) * bodyKeep; // The turn from the last to this one, as a rate const turnKeep = Math.exp(-TURN_AIR * h); turnFrom(spinning, lastTurn, turn); for (let k = 0; k < 3; k++) spinning[k] = (spinning[k]! / h) * turnKeep; }; return { /** The body's middle, its turn as a quaternion, and their rates. Read only. */ middle: middle as Readonly<Vec3>, turn: turn as Readonly<Quat>, velocity: velocity as Readonly<Vec3>, spin: spinning as Readonly<Vec3>, /** Every particle's position: the pivot, the chain to the crown, then each leg. Read only. */ positions: positions as readonly Readonly<Vec3>[], /** The particles of the chain from the pivot to the crown. */ chain: Array.from({ length: LINKS + 1 }, (_, i) => i), /** The particles of each leg from the crown, ending before its foot. */ legs: legs as readonly (readonly number[])[], /** Writes where a foot is in the world into `out`. */ footOf(index: number, out: Vec3) { rotate(out, turn, FEET[index]!); for (let k = 0; k < 3; k++) out[k] = middle[k]! + out[k]!; return out; }, /** A point in the world in the body's own frame, from its middle. */ localOf: (point: ArrayLike<number>): Vec3 => unrotate([0, 0, 0], turn, [ point[0]! - middle[0], point[1]! - middle[1], point[2]! - middle[2], ]), /** Whether a hand holds the body. */ get held() { return hand.held; }, /** Moves the censer on by `dt` seconds, in fixed steps. */ advance(dt: number) { left += dt; // The hand's target moves evenly from where it was to where it is over the advance const steps = Math.floor(left / SUBSTEP); for (let k = 0; k < 3; k++) hand.moving[k] = dt > 0 ? (hand.target[k]! - hand.from[k]!) / dt : 0; let step = 0; while (left >= SUBSTEP) { const t = ++step / steps; for (let k = 0; k < 3; k++) { hand.now[k] = hand.from[k]! + (hand.target[k]! - hand.from[k]!) * t; } substep(SUBSTEP); left -= SUBSTEP; if (!hand.held) free += SUBSTEP; } for (let k = 0; k < 3; k++) hand.from[k] = hand.target[k]!; }, /** Takes hold of the body at a point of it, given from its middle in its own frame. */ hold(grip: ArrayLike<number>) { hand.held = true; for (let k = 0; k < 3; k++) hand.grip[k] = grip[k]!; const at = pointOf(hand.grip); for (let k = 0; k < 3; k++) hand.target[k] = hand.from[k] = hand.now[k] = at[k]!; for (let k = 0; k < 4; k++) hand.turn[k] = turn[k]!; }, /** Pulls the held point toward a point in the world. */ pull(target: ArrayLike<number>) { for (let k = 0; k < 3; k++) hand.target[k] = target[k]!; }, /** Lets go of the body, which keeps its velocity. */ release() { hand.held = false; free = 0; }, /** The angle from straight down the body would swing out to, in radians. */ amplitude: amplitudeOf, };}export type Chain = ReturnType<typeof createChain>;cells.ts
import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';type Vec3 = readonly [number, number, number];export const WORKGROUP = [8, 4, 2] as const;/** A compute entry point over the cells of a box. */export const kernel = tgpu.computeFn({ workgroupSize: [...WORKGROUP], in: { id: d.builtin.globalInvocationId },});/** The workgroups over a box of cells, or over one colour's half of its rows. */export const groups = (size: Vec3, halved = false) => [ Math.ceil((halved ? Math.ceil(size[0] / 2) : size[0]) / WORKGROUP[0]), Math.ceil(size[1] / WORKGROUP[1]), Math.ceil(size[2] / WORKGROUP[2]),];/** Whether a cell lies in a box of `size` cells. */export const inside = (cell: d.v3i, size: d.v3i) => { 'use gpu'; return std.all(std.ge(cell, d.vec3i())) && std.all(std.lt(cell, size));};// Spreading would allocate on every dispatch, which Hermes collects in framesexport const dispatch = (pass: TgpuComputePass, counts: number[]) => pass.dispatchWorkgroups(counts[0]!, counts[1], counts[2]);/** A 3D texture that kernels store to and sample, with a view for each. */export const volume = ( root: TgpuRoot, size: Vec3, format: 'rgba16float' | 'rg16float', name: string,) => { const texture = root .createTexture({ size: [...size], format, dimension: '3d' }) .$usage('storage', 'sampled') .$name(name); return { texture, read: texture.createView(d.texture3d(d.f32)), write: texture.createView(d.textureStorage3d(format, 'write-only')), };};pressure.ts
import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';import { dispatch, groups, inside, kernel } from './cells.ts';type Vec3 = readonly [number, number, number];/** A cell's solid neighbours as bits: -x, +x, -y, +y, -z, +z. All six mark a solid cell. */export const ALL_SOLID = 63;/** The six neighbours of a cell, in the order of their bits in a mask. */export const STEPS = [ d.vec3i(-1, 0, 0), d.vec3i(1, 0, 0), d.vec3i(0, -1, 0), d.vec3i(0, 1, 0), d.vec3i(0, 0, -1), d.vec3i(0, 0, 1),];const exact = { texture: d.texture3d(d.f32), sampleType: 'unfilterable-float' } as const;const solved = d.textureStorage3d('r32float', 'read-write');const written = d.textureStorage3d('r32float', 'write-only');/** One grid of the hierarchy: its unknowns, their right-hand side and each cell's solid neighbours. */const Grid = tgpu.bindGroupLayout({ value: { storageTexture: solved }, rhs: exact, mask: exact,});/** The next coarser grid, as the finer one hands it the residual and takes back its correction. */const Coarser = tgpu.bindGroupLayout({ value: { storageTexture: solved }, rhs: { storageTexture: written },});/** A coarser grid's solid cells, which the coarsening marks in its right-hand side, and its mask. */const Marking = tgpu.bindGroupLayout({ solid: exact, mask: { storageTexture: written },});const maskAt = (cell: d.v3i) => { 'use gpu'; return d.u32(std.textureLoad(Grid.$.mask, cell, 0).x);};/** * One V-cycle over four grids, each half the last along every axis. Red and black sweep once on * the finest grid each way, the dearest, twice on the others and 8 times on the coarsest. Each * sweep moves a value 1.2 times as far as Gauss–Seidel's. */const LEVELS = 4;const FINEST = 1;const SMOOTHING = 2;const COARSEST = 8;const RELAXATION = 1.2;/** * Solves the pressure's Poisson equation on a box of cells by geometric multigrid: a V-cycle of * red–black Gauss–Seidel sweeps over grids halved down to a coarse one, so a push reaches across * the whole box in a step. The equation is the 7-point Laplacian in cells: the sum of a cell's open * neighbours less their count times its own value equals the right-hand side. A solid neighbour * drops out of both, holding the pressure's gradient across it at zero, and a neighbour outside the * box holds 0, so air leaves freely. Each solve starts from the last. */export function createPressure(root: TgpuRoot, cells: Vec3) { const sizes: Vec3[] = [cells]; for (let k = 1; k < LEVELS; k++) { const [x, y, z] = sizes[k - 1]!; sizes.push([Math.ceil(x / 2), Math.ceil(y / 2), Math.ceil(z / 2)]); } const scalar = (size: Vec3, name: string) => root .createTexture({ size: [...size], format: 'r32float', dimension: '3d' }) .$usage('storage', 'sampled') .$name(name); const grids = sizes.map((size, k) => ({ size, value: scalar(size, `pressure ${k}`), rhs: scalar(size, `pressure rhs ${k}`), mask: scalar(size, `pressure mask ${k}`), })); /** * The pressure is 0 half a finest cell past the box on every grid, so a cell just beyond a grid * holds its neighbour inside's value times this, on the line through that 0. */ const beyondOf = (size: d.v3i) => { 'use gpu'; // A cell of this grid spans `span` finest ones const span = d.f32(cells[0]) / d.f32(size.x); return (1 - span) / (1 + span); }; /** * The sum of a cell's open neighbours' values in the box, and the weight of its own value: one * per open neighbour, less what each beyond the box holds of it. */ const neighbours = (cell: d.v3i, mask: number, size: d.v3i) => { 'use gpu'; let sum = d.f32(0); let beyond = d.f32(0); for (const k of tgpu.unroll([0, 1, 2, 3, 4, 5])) { const next = cell + STEPS[k]!; if ((mask & d.u32(1 << k)) === 0) { if (inside(next, size)) sum += std.textureLoad(Grid.$.value, next).x; else beyond += 1; } } return d.vec2f(sum, d.f32(6 - std.countOneBits(mask)) - beyond * beyondOf(size)); }; /** What the cell's value leaves of its right-hand side. */ const residualAt = (cell: d.v3i, size: d.v3i) => { 'use gpu'; const mask = maskAt(cell); if (mask === ALL_SOLID) return d.f32(0); const around = neighbours(cell, mask, size); const own = std.textureLoad(Grid.$.value, cell).x; return std.textureLoad(Grid.$.rhs, cell, 0).x - (around.x - around.y * own); }; // Red and black cells alternate, so each half reads only the other half's latest values const sweep = (parity: number) => kernel(({ id }) => { 'use gpu'; const size = d.vec3i(std.textureDimensions(Grid.$.value)); const cell = d.vec3i( d.i32(id.x * 2 + ((id.y + id.z + parity) & 1)), d.i32(id.y), d.i32(id.z), ); if (!inside(cell, size)) return; const mask = maskAt(cell); if (mask === ALL_SOLID) { std.textureStore(Grid.$.value, cell, d.vec4f()); return; } const around = neighbours(cell, mask, size); const rhs = std.textureLoad(Grid.$.rhs, cell, 0).x; const was = std.textureLoad(Grid.$.value, cell).x; const next = std.mix(was, (around.x - rhs) / around.y, RELAXATION); std.textureStore(Grid.$.value, cell, d.vec4f(next, 0, 0, 0)); }); /** * Hands the coarser grid the residual its eight cells leave, scaled to its doubled spacing, and * clears its correction. */ const restrict = kernel(({ id }) => { 'use gpu'; const coarse = d.vec3i(id); if (!inside(coarse, d.vec3i(std.textureDimensions(Coarser.$.value)))) return; const size = d.vec3i(std.textureDimensions(Grid.$.value)); let total = d.f32(0); for (const z of tgpu.unroll([0, 1])) { for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const fine = coarse * 2 + d.vec3i(x, y, z); if (inside(fine, size)) total += residualAt(fine, size); } } } // The average of eight, times four for a Laplacian over twice the spacing std.textureStore(Coarser.$.rhs, coarse, d.vec4f(total * 0.5, 0, 0, 0)); std.textureStore(Coarser.$.value, coarse, d.vec4f()); }); /** * Adds the coarser grid's correction to each cell, interpolated between its eight nearest, with * those beyond the box continuing the line to its 0. */ const prolong = kernel(({ id }) => { 'use gpu'; const fine = d.vec3i(id); if (!inside(fine, d.vec3i(std.textureDimensions(Grid.$.value)))) return; if (maskAt(fine) === ALL_SOLID) return; const coarse = d.vec3i(std.textureDimensions(Coarser.$.value)); const last = coarse - 1; const beyond = beyondOf(coarse); const at = d.vec3f(fine) * 0.5 - 0.25; const base = d.vec3i(std.floor(at)); const t = at - d.vec3f(base); let correction = d.f32(0); for (const z of tgpu.unroll([0, 1])) { for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const wanted = base + d.vec3i(x, y, z); const corner = std.clamp(wanted, d.vec3i(), last); const reach = std.select(d.vec3f(1), d.vec3f(beyond), std.ne(wanted, corner)); const weight = std.mix(1 - t, t, d.vec3f(x, y, z)) * reach; correction += std.textureLoad(Coarser.$.value, corner).x * weight.x * weight.y * weight.z; } } } const own = std.textureLoad(Grid.$.value, fine).x; std.textureStore(Grid.$.value, fine, d.vec4f(own + correction, 0, 0, 0)); }); /** Marks a coarser cell solid when all its finer cells are, in its right-hand side for now. */ const coarsen = kernel(({ id }) => { 'use gpu'; const coarse = d.vec3i(id); if (!inside(coarse, d.vec3i(std.textureDimensions(Coarser.$.value)))) return; const size = d.vec3i(std.textureDimensions(Grid.$.value)); let solid = true; for (const z of tgpu.unroll([0, 1])) { for (const y of tgpu.unroll([0, 1])) { for (const x of tgpu.unroll([0, 1])) { const fine = coarse * 2 + d.vec3i(x, y, z); if (inside(fine, size) && maskAt(fine) !== ALL_SOLID) solid = false; } } } std.textureStore(Coarser.$.rhs, coarse, d.vec4f(std.select(0, 1, solid), 0, 0, 0)); }); /** Gathers each coarser cell's solid neighbours from the marks into its mask. */ const mark = kernel(({ id }) => { 'use gpu'; const cell = d.vec3i(id); const size = d.vec3i(std.textureDimensions(Marking.$.solid)); if (!inside(cell, size)) return; let mask = d.u32(ALL_SOLID); if (std.textureLoad(Marking.$.solid, cell, 0).x < 0.5) { mask = 0; for (const k of tgpu.unroll([0, 1, 2, 3, 4, 5])) { const next = cell + STEPS[k]!; if (inside(next, size) && std.textureLoad(Marking.$.solid, next, 0).x > 0.5) { mask |= d.u32(1 << k); } } } std.textureStore(Marking.$.mask, cell, d.vec4f(d.f32(mask), 0, 0, 0)); }); const pipelines = { red: root.createComputePipeline({ compute: sweep(0) }), black: root.createComputePipeline({ compute: sweep(1) }), restrict: root.createComputePipeline({ compute: restrict }), prolong: root.createComputePipeline({ compute: prolong }), coarsen: root.createComputePipeline({ compute: coarsen }), mark: root.createComputePipeline({ compute: mark }), }; const read = d.texture3d(d.f32); const bound = grids.map(({ value, rhs, mask }, k) => { const coarser = grids[k + 1]; return { grid: root.createBindGroup(Grid, { value: value.createView(solved), rhs: rhs.createView(read, { sampleType: 'unfilterable-float' }), mask: mask.createView(read, { sampleType: 'unfilterable-float' }), }), coarser: coarser && root.createBindGroup(Coarser, { value: coarser.value.createView(solved), rhs: coarser.rhs.createView(written), }), marking: coarser && root.createBindGroup(Marking, { solid: coarser.rhs.createView(read, { sampleType: 'unfilterable-float' }), mask: coarser.mask.createView(written), }), }; }); const whole = sizes.map((size) => groups(size)); const half = sizes.map((size) => groups(size, true)); const smooth = (pass: TgpuComputePass, k: number, count: number) => { pass.setBindGroup(Grid, bound[k]!.grid); for (let i = 0; i < count; i++) { pass.setPipeline(pipelines.red); dispatch(pass, half[k]!); pass.setPipeline(pipelines.black); dispatch(pass, half[k]!); } }; /** Runs a pipeline that reads grid `k` and writes the next coarser one, over grid `over`. */ const between = ( pass: TgpuComputePass, name: keyof typeof pipelines, k: number, over: number, ) => { pass.setPipeline(pipelines[name]); pass.setBindGroup(Grid, bound[k]!.grid); pass.setBindGroup(Coarser, bound[k]!.coarser!); dispatch(pass, whole[over]!); }; return { /** The finest grid's pressure, which each solve refines, as `read-write` r32float. */ pressure: grids[0]!.value, /** The finest grid's right-hand side, which the caller writes before each solve. */ rhs: grids[0]!.rhs, /** The finest grid's cells' solid neighbours as bits, which the caller writes before each solve. */ mask: grids[0]!.mask, initAsync: () => Promise.all(Object.values(pipelines).map((each) => each.initAsync())), /** Records a solve into a compute pass, after the right-hand side and mask are written. */ solve(pass: TgpuComputePass) { for (let k = 0; k + 1 < LEVELS; k++) { between(pass, 'coarsen', k, k + 1); pass.setPipeline(pipelines.mark); pass.setBindGroup(Marking, bound[k]!.marking!); dispatch(pass, whole[k + 1]!); } for (let k = 0; k + 1 < LEVELS; k++) { smooth(pass, k, k === 0 ? FINEST : SMOOTHING); between(pass, 'restrict', k, k + 1); } smooth(pass, LEVELS - 1, COARSEST); for (let k = LEVELS - 2; k >= 0; k--) { between(pass, 'prolong', k, k); smooth(pass, k, k === 0 ? FINEST : SMOOTHING); } }, destroy() { for (const each of grids) { each.value.destroy(); each.rhs.destroy(); each.mask.destroy(); } }, };}