/**
 * PBDSolver.ts
 *
 * GPU Position-Based Dynamics solver using WebGPU compute shaders.
 * Implements soft-body physics with distance, volume, bending, collision,
 * and attachment constraints solved in parallel via graph coloring.
 *
 * Pipeline per frame:
 *   [1] Apply external forces (gravity, wind)
 *   [2] Predict new positions (Euler integration)
 *   [3] Solve constraints (parallel Gauss-Seidel, N iterations)
 *       ├─ Distance constraints (edges maintain rest length)
 *       ├─ Volume constraints (tetrahedra maintain rest volume)
 *       ├─ Bending constraints (adjacent triangles resist folding)
 *       ├─ Collision constraints (SDF-based penetration resolution)
 *       └─ Attachment constraints (pinned vertices stay fixed)
 *   [4] Update velocities from position delta
 *   [5] Apply damping
 *   [6] Recompute normals
 *
 * @module physics
 */
import type { IVector3, ISoftBodyConfig, ISoftBodyState, ISDFCollider, IConstraintColoring } from './PhysicsTypes';
/**
 * WGSL: Apply external forces and predict positions via semi-implicit Euler
 */
export declare const PBD_PREDICT_SHADER = "\nstruct SimParams {\n  gravity: vec3f,\n  dt: f32,\n  wind: vec3f,\n  damping: f32,\n  numVertices: u32,\n  padding: vec3u,\n}\n\n@group(0) @binding(0) var<storage, read_write> positions: array<f32>;\n@group(0) @binding(1) var<storage, read_write> velocities: array<f32>;\n@group(0) @binding(2) var<storage, read_write> predicted: array<f32>;\n@group(0) @binding(3) var<storage, read> masses: array<f32>;\n@group(0) @binding(4) var<uniform> params: SimParams;\n\n@compute @workgroup_size(256)\nfn cs_predict(@builtin(global_invocation_id) gid: vec3u) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let mass = masses[idx];\n  let i3 = idx * 3u;\n\n  // Load current position and velocity\n  var pos = vec3f(positions[i3], positions[i3 + 1u], positions[i3 + 2u]);\n  var vel = vec3f(velocities[i3], velocities[i3 + 1u], velocities[i3 + 2u]);\n\n  // Skip pinned vertices (mass == 0)\n  if (mass > 0.0) {\n    // Apply gravity\n    vel += params.gravity * params.dt;\n\n    // Apply wind (simple drag model)\n    vel += params.wind * params.dt / mass;\n\n    // Apply damping\n    vel *= params.damping;\n\n    // Predict position\n    pos += vel * params.dt;\n  }\n\n  // Store predicted position\n  predicted[i3]      = pos[0];\n  predicted[i3 + 1u] = pos[1];\n  predicted[i3 + 2u] = pos[2];\n}\n";
/**
 * WGSL: Solve distance constraints (one color group at a time)
 *
 * XPBD formulation (Macklin 2016):
 *   alphaTilde = compliance / dt²
 *   dLambda = (−C − alphaTilde · lambdaAcc) / (wSum + alphaTilde)
 *   lambdaAcc += dLambda          (reset to 0 at substep start, not per-iteration)
 *
 * The lambda accumulator buffer is indexed by constraint index and must be
 * zeroed by the host at the start of each substep (before the first iteration).
 */
export declare const PBD_DISTANCE_SHADER = "\nstruct DistanceConstraint {\n  vertexA: u32,\n  vertexB: u32,\n  restLength: f32,\n  compliance: f32,\n}\n\nstruct SolveParams {\n  dt: f32,\n  numConstraints: u32,\n  iteration: u32,\n  padding: u32,\n}\n\n@group(0) @binding(0) var<storage, read_write> predicted: array<f32>;\n@group(0) @binding(1) var<storage, read> masses: array<f32>;\n@group(0) @binding(2) var<storage, read> constraints: array<DistanceConstraint>;\n@group(0) @binding(3) var<uniform> params: SolveParams;\n@group(0) @binding(4) var<storage, read_write> lambdaAcc: array<f32>;\n\n@compute @workgroup_size(256)\nfn cs_solve_distance(@builtin(global_invocation_id) gid: vec3u) {\n  let cIdx = gid[0];\n  if (cIdx >= params.numConstraints) { return; }\n\n  // Reset lambda accumulator at iteration 0 (start of substep)\n  if (params.iteration == 0u) {\n    lambdaAcc[cIdx] = 0.0;\n  }\n\n  let c = constraints[cIdx];\n  let iA = c.vertexA * 3u;\n  let iB = c.vertexB * 3u;\n\n  let pA = vec3f(predicted[iA], predicted[iA + 1u], predicted[iA + 2u]);\n  let pB = vec3f(predicted[iB], predicted[iB + 1u], predicted[iB + 2u]);\n\n  let diff = pB - pA;\n  let dist = length(diff);\n\n  if (dist < 1e-7) { return; }\n\n  let wA = masses[c.vertexA];\n  let wB = masses[c.vertexB];\n  let invMassA = select(1.0 / wA, 0.0, wA <= 0.0);\n  let invMassB = select(1.0 / wB, 0.0, wB <= 0.0);\n  let wSum = invMassA + invMassB;\n\n  if (wSum < 1e-7) { return; }\n\n  // XPBD: compliance scaled by dt\u00B2; accumulate lambda across iterations\n  let alphaTilde = c.compliance / (params.dt * params.dt);\n  let C = dist - c.restLength;\n  let dLambda = (-C - alphaTilde * lambdaAcc[cIdx]) / (wSum + alphaTilde);\n  lambdaAcc[cIdx] += dLambda;\n\n  let correction = (diff / dist) * dLambda;\n\n  // Apply corrections\n  if (invMassA > 0.0) {\n    predicted[iA]      -= correction[0] * invMassA;\n    predicted[iA + 1u] -= correction[1] * invMassA;\n    predicted[iA + 2u] -= correction[2] * invMassA;\n  }\n  if (invMassB > 0.0) {\n    predicted[iB]      += correction[0] * invMassB;\n    predicted[iB + 1u] += correction[1] * invMassB;\n    predicted[iB + 2u] += correction[2] * invMassB;\n  }\n}\n";
/**
 * WGSL: Solve volume constraints (one tetrahedron at a time)
 *
 * XPBD: accumulates lambda per constraint per substep (binding 4).
 * Host must zero lambdaAcc before the first iteration of each substep.
 */
export declare const PBD_VOLUME_SHADER = "\nstruct VolumeConstraint {\n  v0: u32,\n  v1: u32,\n  v2: u32,\n  v3: u32,\n  restVolume: f32,\n  compliance: f32,\n  padding: vec2u,\n}\n\nstruct SolveParams {\n  dt: f32,\n  numConstraints: u32,\n  iteration: u32,\n  padding: u32,\n}\n\n@group(0) @binding(0) var<storage, read_write> predicted: array<f32>;\n@group(0) @binding(1) var<storage, read> masses: array<f32>;\n@group(0) @binding(2) var<storage, read> constraints: array<VolumeConstraint>;\n@group(0) @binding(3) var<uniform> params: SolveParams;\n@group(0) @binding(4) var<storage, read_write> lambdaAcc: array<f32>;\n\nfn loadPos(idx: u32) -> vec3f {\n  let i = idx * 3u;\n  return vec3f(predicted[i], predicted[i + 1u], predicted[i + 2u]);\n}\n\nfn storePos(idx: u32, p: vec3f) {\n  let i = idx * 3u;\n  predicted[i]      = p[0];\n  predicted[i + 1u] = p[1];\n  predicted[i + 2u] = p[2];\n}\n\n@compute @workgroup_size(64)\nfn cs_solve_volume(@builtin(global_invocation_id) gid: vec3u) {\n  let cIdx = gid[0];\n  if (cIdx >= params.numConstraints) { return; }\n\n  // Reset lambda accumulator at substep start\n  if (params.iteration == 0u) {\n    lambdaAcc[cIdx] = 0.0;\n  }\n\n  let c = constraints[cIdx];\n  let p0 = loadPos(c.v0);\n  let p1 = loadPos(c.v1);\n  let p2 = loadPos(c.v2);\n  let p3 = loadPos(c.v3);\n\n  // Signed volume of tetrahedron = dot(p1-p0, cross(p2-p0, p3-p0)) / 6\n  let d1 = p1 - p0;\n  let d2 = p2 - p0;\n  let d3 = p3 - p0;\n  let volume = dot(d1, cross(d2, d3)) / 6.0;\n\n  let C = volume - c.restVolume;\n  if (abs(C) < 1e-7) { return; }\n\n  // Gradients of volume w.r.t. each vertex\n  let g1 = cross(d2, d3) / 6.0;\n  let g2 = cross(d3, d1) / 6.0;\n  let g3 = cross(d1, d2) / 6.0;\n  let g0 = -(g1 + g2 + g3);\n\n  let w0 = select(1.0 / masses[c.v0], 0.0, masses[c.v0] <= 0.0);\n  let w1 = select(1.0 / masses[c.v1], 0.0, masses[c.v1] <= 0.0);\n  let w2 = select(1.0 / masses[c.v2], 0.0, masses[c.v2] <= 0.0);\n  let w3 = select(1.0 / masses[c.v3], 0.0, masses[c.v3] <= 0.0);\n\n  let denom = w0 * dot(g0, g0) + w1 * dot(g1, g1) + w2 * dot(g2, g2) + w3 * dot(g3, g3);\n\n  if (denom < 1e-7) { return; }\n\n  // XPBD: accumulate lambda across iterations within a substep\n  let alphaTilde = c.compliance / (params.dt * params.dt);\n  let dLambda = (-C - alphaTilde * lambdaAcc[cIdx]) / (denom + alphaTilde);\n  lambdaAcc[cIdx] += dLambda;\n\n  if (w0 > 0.0) { storePos(c.v0, p0 + g0 * dLambda * w0); }\n  if (w1 > 0.0) { storePos(c.v1, p1 + g1 * dLambda * w1); }\n  if (w2 > 0.0) { storePos(c.v2, p2 + g2 * dLambda * w2); }\n  if (w3 > 0.0) { storePos(c.v3, p3 + g3 * dLambda * w3); }\n}\n";
/**
 * WGSL: SDF collision constraint solve
 */
export declare const PBD_COLLISION_SHADER = "\nstruct SDFParams {\n  boundsMin: vec3f,\n  cellSize: f32,\n  boundsMax: vec3f,\n  friction: f32,\n  gridSizeX: u32,\n  gridSizeY: u32,\n  gridSizeZ: u32,\n  numVertices: u32,\n  collisionMargin: f32,\n  padding: vec3f,\n}\n\n@group(0) @binding(0) var<storage, read_write> predicted: array<f32>;\n@group(0) @binding(1) var<storage, read> masses: array<f32>;\n@group(0) @binding(2) var<storage, read> sdfData: array<f32>;\n@group(0) @binding(3) var<uniform> params: SDFParams;\n\nfn sampleSDF(pos: vec3f) -> f32 {\n  let local = (pos - params.boundsMin) / params.cellSize;\n  let gi = vec3u(clamp(vec3i(floor(local)), vec3i(0), vec3i(\n    i32(params.gridSizeX) - 2,\n    i32(params.gridSizeY) - 2,\n    i32(params.gridSizeZ) - 2\n  )));\n  let f = fract(local);\n\n  // Trilinear interpolation\n  let idx000 = gi[0] + gi[1] * params.gridSizeX + gi[2] * params.gridSizeX * params.gridSizeY;\n  let idx100 = idx000 + 1u;\n  let idx010 = idx000 + params.gridSizeX;\n  let idx110 = idx010 + 1u;\n  let idx001 = idx000 + params.gridSizeX * params.gridSizeY;\n  let idx101 = idx001 + 1u;\n  let idx011 = idx001 + params.gridSizeX;\n  let idx111 = idx011 + 1u;\n\n  let c00 = mix(sdfData[idx000], sdfData[idx100], f[0]);\n  let c10 = mix(sdfData[idx010], sdfData[idx110], f[0]);\n  let c01 = mix(sdfData[idx001], sdfData[idx101], f[0]);\n  let c11 = mix(sdfData[idx011], sdfData[idx111], f[0]);\n  let c0 = mix(c00, c10, f[1]);\n  let c1 = mix(c01, c11, f[1]);\n  return mix(c0, c1, f[2]);\n}\n\nfn sdfGradient(pos: vec3f) -> vec3f {\n  let eps = params.cellSize * 0.5;\n  return normalize(vec3f(\n    sampleSDF(pos + vec3f(eps, 0.0, 0.0)) - sampleSDF(pos - vec3f(eps, 0.0, 0.0)),\n    sampleSDF(pos + vec3f(0.0, eps, 0.0)) - sampleSDF(pos - vec3f(0.0, eps, 0.0)),\n    sampleSDF(pos + vec3f(0.0, 0.0, eps)) - sampleSDF(pos - vec3f(0.0, 0.0, eps))\n  ));\n}\n\n@compute @workgroup_size(256)\nfn cs_solve_collision(@builtin(global_invocation_id) gid: vec3u) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n  if (masses[idx] <= 0.0) { return; }\n\n  let i3 = idx * 3u;\n  let pos = vec3f(predicted[i3], predicted[i3 + 1u], predicted[i3 + 2u]);\n\n  // Check if inside SDF bounds\n  if (any(pos < params.boundsMin) || any(pos > params.boundsMax)) { return; }\n\n  let dist = sampleSDF(pos);\n\n  if (dist < params.collisionMargin) {\n    let normal = sdfGradient(pos);\n    let penetration = params.collisionMargin - dist;\n\n    // Push vertex out along SDF gradient\n    let correction = normal * penetration;\n    predicted[i3]      += correction[0];\n    predicted[i3 + 1u] += correction[1];\n    predicted[i3 + 2u] += correction[2];\n  }\n}\n";
/**
 * WGSL: Update velocities from position delta and apply damping
 */
export declare const PBD_VELOCITY_SHADER = "\nstruct VelParams {\n  dt: f32,\n  damping: f32,\n  numVertices: u32,\n  padding: u32,\n}\n\n@group(0) @binding(0) var<storage, read> positions: array<f32>;\n@group(0) @binding(1) var<storage, read_write> velocities: array<f32>;\n@group(0) @binding(2) var<storage, read> predicted: array<f32>;\n@group(0) @binding(3) var<uniform> params: VelParams;\n\n@compute @workgroup_size(256)\nfn cs_update_velocity(@builtin(global_invocation_id) gid: vec3u) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let i3 = idx * 3u;\n  let invDt = 1.0 / params.dt;\n\n  velocities[i3]      = (predicted[i3]      - positions[i3])      * invDt;\n  velocities[i3 + 1u] = (predicted[i3 + 1u] - positions[i3 + 1u]) * invDt;\n  velocities[i3 + 2u] = (predicted[i3 + 2u] - positions[i3 + 2u]) * invDt;\n}\n";
/**
 * WGSL: Copy predicted back to positions (finalize step)
 */
export declare const PBD_FINALIZE_SHADER = "\nstruct FinalizeParams {\n  numVertices: u32,\n  padding: vec3u,\n}\n\n@group(0) @binding(0) var<storage, read_write> positions: array<f32>;\n@group(0) @binding(1) var<storage, read> predicted: array<f32>;\n@group(0) @binding(2) var<uniform> params: FinalizeParams;\n\n@compute @workgroup_size(256)\nfn cs_finalize(@builtin(global_invocation_id) gid: vec3u) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let i3 = idx * 3u;\n  positions[i3]      = predicted[i3];\n  positions[i3 + 1u] = predicted[i3 + 1u];\n  positions[i3 + 2u] = predicted[i3 + 2u];\n}\n";
/**
 * WGSL: Recompute vertex normals from triangle indices
 */
export declare const PBD_NORMALS_SHADER = "\nstruct NormalParams {\n  numTriangles: u32,\n  numVertices: u32,\n  padding: vec2u,\n}\n\n@group(0) @binding(0) var<storage, read> positions: array<f32>;\n@group(0) @binding(1) var<storage, read_write> normals: array<atomic<i32>>;\n@group(0) @binding(2) var<storage, read> indices: array<u32>;\n@group(0) @binding(3) var<uniform> params: NormalParams;\n\nfn loadPos(idx: u32) -> vec3f {\n  let i = idx * 3u;\n  return vec3f(positions[i], positions[i + 1u], positions[i + 2u]);\n}\n\n// Atomic accumulate normal via fixed-point (multiply by 1e6, cast to i32)\nfn atomicAddNormal(vertIdx: u32, n: vec3f) {\n  let i = vertIdx * 3u;\n  let scale = 1000000.0;\n  atomicAdd(&normals[i],      i32(n[0] * scale));\n  atomicAdd(&normals[i + 1u], i32(n[1] * scale));\n  atomicAdd(&normals[i + 2u], i32(n[2] * scale));\n}\n\n@compute @workgroup_size(256)\nfn cs_compute_normals(@builtin(global_invocation_id) gid: vec3u) {\n  let triIdx = gid[0];\n  if (triIdx >= params.numTriangles) { return; }\n\n  let i0 = indices[triIdx * 3u];\n  let i1 = indices[triIdx * 3u + 1u];\n  let i2 = indices[triIdx * 3u + 2u];\n\n  let p0 = loadPos(i0);\n  let p1 = loadPos(i1);\n  let p2 = loadPos(i2);\n\n  let faceNormal = cross(p1 - p0, p2 - p0);\n\n  atomicAddNormal(i0, faceNormal);\n  atomicAddNormal(i1, faceNormal);\n  atomicAddNormal(i2, faceNormal);\n}\n";
/**
 * WGSL: Normalize accumulated vertex normals
 */
export declare const PBD_NORMALIZE_SHADER = "\nstruct NormalizeParams {\n  numVertices: u32,\n  padding: vec3u,\n}\n\n@group(0) @binding(0) var<storage, read_write> normalsI32: array<i32>;\n@group(0) @binding(1) var<storage, read_write> normalsF32: array<f32>;\n@group(0) @binding(2) var<uniform> params: NormalizeParams;\n\n@compute @workgroup_size(256)\nfn cs_normalize_normals(@builtin(global_invocation_id) gid: vec3u) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let i3 = idx * 3u;\n  let scale = 1.0 / 1000000.0;\n  var n = vec3f(\n    f32(normalsI32[i3])      * scale,\n    f32(normalsI32[i3 + 1u]) * scale,\n    f32(normalsI32[i3 + 2u]) * scale\n  );\n\n  let len = length(n);\n  if (len > 1e-7) {\n    n /= len;\n  } else {\n    n = vec3f(0.0, 1.0, 0.0);\n  }\n\n  normalsF32[i3]      = n[0];\n  normalsF32[i3 + 1u] = n[1];\n  normalsF32[i3 + 2u] = n[2];\n\n  // Reset atomic accumulator for next frame\n  normalsI32[i3]      = 0;\n  normalsI32[i3 + 1u] = 0;\n  normalsI32[i3 + 2u] = 0;\n}\n";
/**
 * WGSL compute shader for bending constraints.
 * Maintains rest dihedral angles between adjacent triangle pairs.
 * Uses XPBD formulation with compliance and lambda accumulation.
 *
 * Buffers:
 *   binding(0) params      — uniform: dt, numConstraints, compliance, iteration
 *   binding(1) positions   — read_write predicted positions (flat xyz)
 *   binding(2) masses      — read per-vertex masses (0 = pinned)
 *   binding(3) constraints — read 4 u32 per constraint [v0, v1, v2, v3]
 *   binding(4) restAngles  — read rest dihedral angle per constraint
 *   binding(5) lambdaAcc   — read_write lambda accumulator (reset at iteration 0)
 */
export declare const PBD_BENDING_SHADER = "\nstruct Params {\n  dt: f32,\n  numConstraints: u32,\n  compliance: f32,\n  iteration: u32,\n};\n\n@group(0) @binding(0) var<uniform> params: Params;\n@group(0) @binding(1) var<storage, read_write> positions: array<f32>;\n@group(0) @binding(2) var<storage, read> masses: array<f32>;\n@group(0) @binding(3) var<storage, read> constraints: array<u32>;\n@group(0) @binding(4) var<storage, read> restAngles: array<f32>;\n@group(0) @binding(5) var<storage, read_write> lambdaAcc: array<f32>;\n\nfn loadPos(i: u32) -> vec3<f32> {\n  return vec3<f32>(positions[i * 3u], positions[i * 3u + 1u], positions[i * 3u + 2u]);\n}\n\nfn storePos(i: u32, p: vec3<f32>) {\n  positions[i * 3u] = p[0];\n  positions[i * 3u + 1u] = p[1];\n  positions[i * 3u + 2u] = p[2];\n}\n\n@compute @workgroup_size(64)\nfn main(@builtin(global_invocation_id) gid: vec3<u32>) {\n  let idx = gid[0];\n  if (idx >= params.numConstraints) { return; }\n\n  // Reset lambda accumulator at substep start\n  if (params.iteration == 0u) {\n    lambdaAcc[idx] = 0.0;\n  }\n\n  // v0-v1 = shared edge, v2 = opposite on face 1, v3 = opposite on face 2\n  let v0 = constraints[idx * 4u];\n  let v1 = constraints[idx * 4u + 1u];\n  let v2 = constraints[idx * 4u + 2u];\n  let v3 = constraints[idx * 4u + 3u];\n\n  let p0 = loadPos(v0);\n  let p1 = loadPos(v1);\n  let p2 = loadPos(v2);\n  let p3 = loadPos(v3);\n\n  // Shared edge\n  let e = p1 - p0;\n  let eLen = length(e);\n  if (eLen < 1e-7) { return; }\n  let eNorm = e / eLen;\n\n  // Vectors from v0 to opposite vertices\n  let d2 = p2 - p0;\n  let d3 = p3 - p0;\n\n  // Face normals (unnormalized, magnitude = 2 * triangle area)\n  let n1 = cross(e, d2);\n  let n2 = cross(e, d3);\n  let n1Len = length(n1);\n  let n2Len = length(n2);\n  if (n1Len < 1e-7 || n2Len < 1e-7) { return; }\n\n  let n1n = n1 / n1Len;\n  let n2n = n2 / n2Len;\n\n  // Dihedral angle via atan2 for full [-pi, pi] range\n  let cosTheta = clamp(dot(n1n, n2n), -1.0, 1.0);\n  let sinTheta = dot(cross(n1n, n2n), eNorm);\n  let theta = atan2(sinTheta, cosTheta);\n\n  // Constraint: C = theta - restAngle\n  let restAngle = restAngles[idx];\n  let C = theta - restAngle;\n  if (abs(C) < 1e-6) { return; }\n\n  // Perpendicular distances from opposite vertices to shared edge\n  let h2 = n1Len / eLen;\n  let h3 = n2Len / eLen;\n  if (h2 < 1e-7 || h3 < 1e-7) { return; }\n\n  // Gradients of theta w.r.t. vertex positions (Bridson formulation)\n  let g2 = n1n / h2;\n  let g3 = -n2n / h3;\n\n  // Parametric projections along shared edge for gradient distribution\n  let eDotE = eLen * eLen;\n  let s2 = dot(d2, e) / eDotE;\n  let s3 = dot(d3, e) / eDotE;\n\n  // Edge vertex gradients (sum of all gradients = 0, translation invariance)\n  let g0 = -(1.0 - s2) * g2 - (1.0 - s3) * g3;\n  let g1 = -s2 * g2 - s3 * g3;\n\n  // Inverse masses\n  let w0 = select(0.0, 1.0 / masses[v0], masses[v0] > 0.0);\n  let w1 = select(0.0, 1.0 / masses[v1], masses[v1] > 0.0);\n  let w2 = select(0.0, 1.0 / masses[v2], masses[v2] > 0.0);\n  let w3 = select(0.0, 1.0 / masses[v3], masses[v3] > 0.0);\n\n  // XPBD: alphaTilde = compliance / dt\u00B2; accumulate lambda across iterations\n  let alphaTilde = params.compliance / (params.dt * params.dt);\n  let denom = w0 * dot(g0, g0) + w1 * dot(g1, g1)\n            + w2 * dot(g2, g2) + w3 * dot(g3, g3) + alphaTilde;\n  if (denom < 1e-10) { return; }\n\n  let dLambda = (-C - alphaTilde * lambdaAcc[idx]) / denom;\n  lambdaAcc[idx] += dLambda;\n\n  // Apply position corrections\n  if (w0 > 0.0) { storePos(v0, p0 + g0 * dLambda * w0); }\n  if (w1 > 0.0) { storePos(v1, p1 + g1 * dLambda * w1); }\n  if (w2 > 0.0) { storePos(v2, p2 + g2 * dLambda * w2); }\n  if (w3 > 0.0) { storePos(v3, p3 + g3 * dLambda * w3); }\n}\n";
/**
 * WGSL compute shader for attachment constraints.
 * Pins vertices to target positions with XPBD compliance.
 * compliance=0 → hard pin, compliance>0 → soft spring.
 *
 * Buffers:
 *   binding(1) positions    — read_write predicted positions (flat xyz)
 *   binding(2) masses       — read per-vertex masses
 *   binding(3) vertexIndices — read vertex index per constraint
 *   binding(4) targets      — read target positions (flat xyz, 3 per constraint)
 *   binding(5) compliances  — read compliance per constraint
 */
export declare const PBD_ATTACHMENT_SHADER = "\nstruct Params {\n  dt: f32,\n  numConstraints: u32,\n  _pad0: f32,\n  _pad1: f32,\n};\n\n@group(0) @binding(0) var<uniform> params: Params;\n@group(0) @binding(1) var<storage, read_write> positions: array<f32>;\n@group(0) @binding(2) var<storage, read> masses: array<f32>;\n@group(0) @binding(3) var<storage, read> vertexIndices: array<u32>;\n@group(0) @binding(4) var<storage, read> targets: array<f32>;\n@group(0) @binding(5) var<storage, read> compliances: array<f32>;\n\n@compute @workgroup_size(64)\nfn main(@builtin(global_invocation_id) gid: vec3<u32>) {\n  let idx = gid[0];\n  if (idx >= params.numConstraints) { return; }\n\n  let vi = vertexIndices[idx];\n  let mass = masses[vi];\n  if (mass <= 0.0) { return; }  // Pinned vertex, already at target\n\n  let px = positions[vi * 3u];\n  let py = positions[vi * 3u + 1u];\n  let pz = positions[vi * 3u + 2u];\n\n  let tx = targets[idx * 3u];\n  let ty = targets[idx * 3u + 1u];\n  let tz = targets[idx * 3u + 2u];\n\n  let dx = tx - px;\n  let dy = ty - py;\n  let dz = tz - pz;\n\n  let dist = sqrt(dx * dx + dy * dy + dz * dz);\n  if (dist < 1e-7) { return; }\n\n  let compliance = compliances[idx];\n  let w = 1.0 / mass;\n  let alpha = compliance / (params.dt * params.dt);\n\n  // XPBD: C = dist, gradient = normalize(d), single vertex\n  let lambda = dist / (w + alpha);\n\n  let nx = dx / dist;\n  let ny = dy / dist;\n  let nz = dz / dist;\n\n  positions[vi * 3u]      = px + nx * lambda * w;\n  positions[vi * 3u + 1u] = py + ny * lambda * w;\n  positions[vi * 3u + 2u] = pz + nz * lambda * w;\n}\n";
/**
 * WGSL: Solve density constraints (SPH-style).
 * Each fluid particle computes local density from neighbors and projects
 * towards the rest density using XPBD compliance.
 *
 * Buffer layout:
 *   binding(0) params       — uniform: dt, numConstraints, restDensity, kernelRadius
 *   binding(1) predicted    — read/write predicted positions (flat xyz)
 *   binding(2) masses       — read per-vertex masses
 *   binding(3) particleIdx  — read: which vertices are fluid particles
 *   binding(4) neighborList — read: flattened neighbor indices (maxNeighbors per particle)
 *   binding(5) neighborCnt  — read: neighbor count per particle
 */
export declare const PBD_DENSITY_SHADER = "\nstruct DensityParams {\n  dt: f32,\n  numParticles: u32,\n  restDensity: f32,\n  kernelRadius: f32,\n  compliance: f32,\n  maxNeighbors: u32,\n  _pad0: u32,\n  _pad1: u32,\n}\n\n@group(0) @binding(0) var<uniform> params: DensityParams;\n@group(0) @binding(1) var<storage, read_write> predicted: array<f32>;\n@group(0) @binding(2) var<storage, read> masses: array<f32>;\n@group(0) @binding(3) var<storage, read> particleIdx: array<u32>;\n@group(0) @binding(4) var<storage, read> neighborList: array<u32>;\n@group(0) @binding(5) var<storage, read> neighborCnt: array<u32>;\n@group(0) @binding(6) var<storage, read_write> lambdas: array<f32>;\n\n// Poly6 kernel\nfn poly6(r2: f32, h: f32) -> f32 {\n  let h2 = h * h;\n  if (r2 >= h2) { return 0.0; }\n  let diff = h2 - r2;\n  return 315.0 / (64.0 * 3.14159265 * pow(h, 9.0)) * diff * diff * diff;\n}\n\n// Spiky gradient magnitude / r (to multiply with direction)\nfn spikyGrad(r: f32, h: f32) -> f32 {\n  if (r >= h || r < 1e-7) { return 0.0; }\n  let diff = h - r;\n  return -45.0 / (3.14159265 * pow(h, 6.0)) * diff * diff;\n}\n\nfn loadPos(idx: u32) -> vec3f {\n  let i = idx * 3u;\n  return vec3f(predicted[i], predicted[i + 1u], predicted[i + 2u]);\n}\n\n// neighborList stores PARTICLE indices (into particleIdx/lambdas arrays),\n// NOT raw vertex indices.  Using particle indices allows direct lambda lookup.\n@compute @workgroup_size(256)\nfn cs_compute_density_lambda(@builtin(global_invocation_id) gid: vec3u) {\n  let pIdx = gid[0];\n  if (pIdx >= params.numParticles) { return; }\n\n  let vi = particleIdx[pIdx];\n  let pi = loadPos(vi);\n  let h = params.kernelRadius;\n  let nCount = min(neighborCnt[pIdx], params.maxNeighbors);\n  let nBase = pIdx * params.maxNeighbors;\n\n  // Compute density (neighborList entries are PARTICLE indices)\n  var density = poly6(0.0, h) * masses[vi]; // self-contribution\n  for (var n = 0u; n < nCount; n++) {\n    let nPIdx = neighborList[nBase + n];   // neighbor particle index\n    let nj = particleIdx[nPIdx];           // neighbor vertex index\n    let pj = loadPos(nj);\n    let diff = pi - pj;\n    let r2 = dot(diff, diff);\n    density += poly6(r2, h) * masses[nj];\n  }\n\n  // Constraint: C = density / restDensity - 1\n  let C = density / params.restDensity - 1.0;\n  if (C <= 0.0) { lambdas[pIdx] = 0.0; return; }\n\n  // Compute denominator: sum of squared gradient magnitudes\n  var gradSumSq: f32 = 0.0;\n  var gradI = vec3f(0.0);\n\n  for (var n = 0u; n < nCount; n++) {\n    let nPIdx = neighborList[nBase + n];   // neighbor particle index\n    let nj = particleIdx[nPIdx];           // neighbor vertex index\n    let pj = loadPos(nj);\n    let diff = pi - pj;\n    let r = length(diff);\n    let s = spikyGrad(r, h) / params.restDensity;\n    if (r > 1e-7) {\n      let gradJ = diff * (s / r);\n      gradSumSq += dot(gradJ, gradJ);\n      gradI += gradJ;\n    }\n  }\n\n  gradSumSq += dot(gradI, gradI);\n\n  let alpha = params.compliance / (params.dt * params.dt);\n  lambdas[pIdx] = -C / (gradSumSq + alpha);\n}\n\n@compute @workgroup_size(256)\nfn cs_apply_density(@builtin(global_invocation_id) gid: vec3u) {\n  let pIdx = gid[0];\n  if (pIdx >= params.numParticles) { return; }\n\n  let vi = particleIdx[pIdx];\n  let pi = loadPos(vi);\n  let h = params.kernelRadius;\n  let nCount = min(neighborCnt[pIdx], params.maxNeighbors);\n  let nBase = pIdx * params.maxNeighbors;\n  let lambdaI = lambdas[pIdx];\n\n  var correction = vec3f(0.0);\n  for (var n = 0u; n < nCount; n++) {\n    let nPIdx = neighborList[nBase + n];   // neighbor particle index\n    let nj = particleIdx[nPIdx];           // neighbor vertex index\n    let lambdaJ = lambdas[nPIdx];          // neighbor's lambda\n    let pj = loadPos(nj);\n    let diff = pi - pj;\n    let r = length(diff);\n    let s = spikyGrad(r, h) / params.restDensity;\n    if (r > 1e-7) {\n      // Symmetric pressure gradient: (\u03BB_i + \u03BB_j) per PBD-fluid (M\u00FCller 2003)\n      correction += diff * ((lambdaI + lambdaJ) * s / r);\n    }\n  }\n\n  let i3 = vi * 3u;\n  predicted[i3]      += correction[0];\n  predicted[i3 + 1u] += correction[1];\n  predicted[i3 + 2u] += correction[2];\n}\n";
/**
 * WGSL compute shader: build spatial hash grid.
 * Hashes each vertex position into a 3D grid cell and atomically appends
 * the vertex index to that cell's list via a prefix-sum-style counter array.
 *
 * Two-pass approach:
 *   Pass 1 (this shader): count vertices per cell + write vertex→cell mapping
 *   CPU prefix sum: compute cell offsets from counts
 *   Pass 2 (resolve shader): iterate 27 neighbors, resolve penetrations
 */
export declare const PBD_HASH_BUILD_SHADER = "\nstruct Params {\n  numVertices: u32,\n  cellSize: f32,\n  gridDimX: u32,\n  gridDimY: u32,\n  gridDimZ: u32,\n  originX: f32,\n  originY: f32,\n  originZ: f32,\n};\n\n@group(0) @binding(0) var<uniform> params: Params;\n@group(0) @binding(1) var<storage, read> positions: array<f32>;\n@group(0) @binding(2) var<storage, read_write> cellCounts: array<atomic<u32>>;\n@group(0) @binding(3) var<storage, read_write> vertexCells: array<u32>;\n\nfn hashCell(ix: u32, iy: u32, iz: u32) -> u32 {\n  return ix + iy * params.gridDimX + iz * params.gridDimX * params.gridDimY;\n}\n\n@compute @workgroup_size(256)\nfn main(@builtin(global_invocation_id) gid: vec3<u32>) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let px = positions[idx * 3u];\n  let py = positions[idx * 3u + 1u];\n  let pz = positions[idx * 3u + 2u];\n\n  let ix = clamp(u32(floor((px - params.originX) / params.cellSize)), 0u, params.gridDimX - 1u);\n  let iy = clamp(u32(floor((py - params.originY) / params.cellSize)), 0u, params.gridDimY - 1u);\n  let iz = clamp(u32(floor((pz - params.originZ) / params.cellSize)), 0u, params.gridDimZ - 1u);\n\n  let cell = hashCell(ix, iy, iz);\n  vertexCells[idx] = cell;\n  atomicAdd(&cellCounts[cell], 1u);\n}\n";
/**
 * WGSL compute shader: resolve self-collisions.
 * For each vertex, iterate all vertices in the 27 neighboring cells
 * and push apart any that are closer than the collision radius.
 *
 * Buffers:
 *   binding(1) positions    — read_write predicted positions
 *   binding(2) masses       — read per-vertex masses
 *   binding(3) vertexCells  — read cell index per vertex
 *   binding(4) cellOffsets  — read prefix-sum offsets per cell (from CPU)
 *   binding(5) cellCounts   — read vertex count per cell
 *   binding(6) sortedVerts  — read vertex indices sorted by cell
 */
export declare const PBD_SELF_COLLISION_SHADER = "\nstruct Params {\n  numVertices: u32,\n  collisionRadius: f32,\n  gridDimX: u32,\n  gridDimY: u32,\n  gridDimZ: u32,\n  originX: f32,\n  originY: f32,\n  originZ: f32,\n  cellSize: f32,\n  _pad: f32,\n  _pad2: f32,\n  _pad3: f32,\n};\n\n@group(0) @binding(0) var<uniform> params: Params;\n@group(0) @binding(1) var<storage, read_write> positions: array<f32>;\n@group(0) @binding(2) var<storage, read> masses: array<f32>;\n@group(0) @binding(3) var<storage, read> vertexCells: array<u32>;\n@group(0) @binding(4) var<storage, read> cellOffsets: array<u32>;\n@group(0) @binding(5) var<storage, read> cellCounts2: array<u32>;\n@group(0) @binding(6) var<storage, read> sortedVerts: array<u32>;\n\nfn hashCell(ix: u32, iy: u32, iz: u32) -> u32 {\n  return ix + iy * params.gridDimX + iz * params.gridDimX * params.gridDimY;\n}\n\n@compute @workgroup_size(64)\nfn main(@builtin(global_invocation_id) gid: vec3<u32>) {\n  let idx = gid[0];\n  if (idx >= params.numVertices) { return; }\n\n  let mi = masses[idx];\n  if (mi <= 0.0) { return; }\n  let wi = 1.0 / mi;\n\n  let px = positions[idx * 3u];\n  let py = positions[idx * 3u + 1u];\n  let pz = positions[idx * 3u + 2u];\n\n  let ix = clamp(u32(floor((px - params.originX) / params.cellSize)), 0u, params.gridDimX - 1u);\n  let iy = clamp(u32(floor((py - params.originY) / params.cellSize)), 0u, params.gridDimY - 1u);\n  let iz = clamp(u32(floor((pz - params.originZ) / params.cellSize)), 0u, params.gridDimZ - 1u);\n\n  let radius = params.collisionRadius;\n  let radiusSq = radius * radius;\n\n  var corrX = 0.0;\n  var corrY = 0.0;\n  var corrZ = 0.0;\n\n  // Iterate 27 neighbor cells\n  for (var dz: i32 = -1; dz <= 1; dz++) {\n    for (var dy: i32 = -1; dy <= 1; dy++) {\n      for (var dx: i32 = -1; dx <= 1; dx++) {\n        let nx = i32(ix) + dx;\n        let ny = i32(iy) + dy;\n        let nz = i32(iz) + dz;\n\n        if (nx < 0 || nx >= i32(params.gridDimX) ||\n            ny < 0 || ny >= i32(params.gridDimY) ||\n            nz < 0 || nz >= i32(params.gridDimZ)) { continue; }\n\n        let cell = hashCell(u32(nx), u32(ny), u32(nz));\n        let offset = cellOffsets[cell];\n        let count = cellCounts2[cell];\n\n        for (var k: u32 = 0u; k < count; k++) {\n          let j = sortedVerts[offset + k];\n          if (j == idx) { continue; }\n\n          let mj = masses[j];\n          if (mj <= 0.0) { continue; }\n\n          let qx = positions[j * 3u] - px;\n          let qy = positions[j * 3u + 1u] - py;\n          let qz = positions[j * 3u + 2u] - pz;\n          let distSq = qx * qx + qy * qy + qz * qz;\n\n          if (distSq < radiusSq && distSq > 1e-10) {\n            let dist = sqrt(distSq);\n            let penetration = radius - dist;\n            let wj = 1.0 / mj;\n            let wSum = wi + wj;\n            let scale = penetration / (wSum * dist);\n\n            // Push this vertex away (negative direction)\n            corrX -= qx * scale * wi;\n            corrY -= qy * scale * wi;\n            corrZ -= qz * scale * wi;\n          }\n        }\n      }\n    }\n  }\n\n  positions[idx * 3u]      = px + corrX;\n  positions[idx * 3u + 1u] = py + corrY;\n  positions[idx * 3u + 2u] = pz + corrZ;\n}\n";
/**
 * Greedy edge coloring for distance/bending constraints.
 *
 * Two constraints that share a vertex MUST receive different colors so that
 * same-color constraints can execute in parallel on the GPU without data races.
 *
 * Correct algorithm: each vertex tracks the SET of all colors already
 * assigned to constraints incident to it.  For each new constraint, the
 * forbidden set is the UNION of both endpoint's incident-color sets.  The
 * constraint receives the lowest color not in that union.
 *
 * The previous implementation stored only a single color per vertex
 * (vertexColor[v] = color), overwriting earlier assignments.  For a vertex
 * shared by k constraints, only the last assignment was remembered, so two
 * earlier constraints on the same vertex could receive the same color.
 */
export declare function colorConstraints(constraints: Array<{
    vertexA: number;
    vertexB: number;
}>, numVertices: number): IConstraintColoring;
/**
 * Generate tetrahedral mesh from a surface triangle mesh.
 * Uses simple fan tetrahedralization from centroid — suitable for convex
 * and mildly concave meshes. For complex shapes, use Delaunay tetrahedralization.
 */
export declare function generateTetrahedra(positions: Float32Array, indices: Uint32Array): {
    tetIndices: Uint32Array;
    restVolumes: Float32Array;
};
/**
 * Extract unique edges from triangle indices for distance constraints.
 */
export declare function extractEdges(indices: Uint32Array, _numVertices: number): {
    edges: Uint32Array;
    restLengths: Float32Array;
    positions: Float32Array;
} | {
    edges: Uint32Array;
};
/**
 * Compute rest lengths for distance constraints.
 */
export declare function computeRestLengths(positions: Float32Array, edges: Uint32Array): Float32Array;
/**
 * Extract bending constraint pairs from a triangle mesh.
 * Finds pairs of triangles sharing an edge and computes rest dihedral angles.
 *
 * @returns constraints — flat Uint32Array, 4 indices per constraint [v0, v1, v2, v3]
 *          where (v0,v1) is the shared edge and v2, v3 are the opposite vertices.
 * @returns restAngles — Float32Array of rest dihedral angles (one per constraint).
 */
export declare function extractBendingPairs(indices: Uint32Array, positions: Float32Array): {
    constraints: Uint32Array;
    restAngles: Float32Array;
};
/**
 * Generate a signed distance field from a triangle mesh.
 *
 * Distance: brute-force closest-triangle for each grid cell.
 * Sign: ray-crossing parity test using 3-axis majority vote for robustness
 * on degenerate geometry (grazing rays, edge/vertex hits).  A point is
 * inside the mesh when the majority of the three axis-aligned ray-crossing
 * counts are odd.
 *
 * Sign convention: negative = inside, positive = outside.
 * For compliance=0 (no-compliance XPBD) the sign does not affect constraint
 * solve direction — this is purely for query correctness.
 *
 * For production use, replace the O(N·M) brute force with a BVH.
 */
export declare function generateSDF(vertices: Float32Array, indices: Uint32Array, gridSize: number, padding?: number): ISDFCollider;
/**
 * CPU-side PBD solver. Mirrors the GPU pipeline for environments without WebGPU.
 */
export declare class PBDSolverCPU {
    private config;
    private state;
    private distanceConstraints;
    private volumeConstraints;
    private bendingConstraints;
    private attachmentConstraints;
    private densityParticles;
    private densityNeighbors;
    private densityRestDensity;
    private densityKernelRadius;
    private densityCompliance;
    private coloring;
    private sdfColliders;
    /**
     * Per-constraint lambda accumulators for XPBD (Macklin 2016).
     * Reset to 0 at the start of each substep; accumulated across iterations
     * within the same substep.  Indexed in the same order as the corresponding
     * constraint arrays (distance, volume, bending).
     */
    private distanceLambdaAcc;
    private volumeLambdaAcc;
    private bendingLambdaAcc;
    constructor(config: ISoftBodyConfig);
    private buildConstraints;
    private computeTetVolume;
    /**
     * Add an SDF collider for collision detection
     */
    addCollider(collider: ISDFCollider): void;
    /**
     * Pin a vertex to a world position.
     * @param compliance 0 = hard pin, >0 = soft spring (XPBD compliant)
     */
    pinVertex(vertexIndex: number, position: IVector3, compliance?: number): void;
    /**
     * Unpin a vertex
     */
    unpinVertex(vertexIndex: number): void;
    /**
     * Apply an impulse at a position (for grab interaction)
     */
    applyImpulse(position: IVector3, force: IVector3, radius: number): void;
    /**
     * Step the simulation forward by dt seconds
     */
    step(dt: number): ISoftBodyState;
    /**
     * Solve a single distance constraint using XPBD (Macklin 2016).
     *
     * XPBD update rule (per iteration within a substep):
     *   alphaTilde = compliance / dt²
     *   dLambda = (−C − alphaTilde · lambdaAcc) / (wSum + alphaTilde)
     *   lambdaAcc += dLambda
     *
     * lambdaAcc is reset to 0 once per substep (before the first iteration).
     * For compliance=0, alphaTilde=0, the term vanishes and the update reduces
     * to standard PBD — fully backward-compatible.
     *
     * @param accIdx  Index into distanceLambdaAcc for this constraint.
     */
    private solveDistanceConstraint;
    /**
     * Solve a single volume (tetrahedral) constraint using XPBD.
     * @param accIdx  Index into volumeLambdaAcc for this constraint.
     */
    private solveVolumeConstraint;
    /**
     * Solve a single bending constraint using XPBD.
     * @param accIdx  Index into bendingLambdaAcc for this constraint.
     */
    private solveBendingConstraint;
    private solveAttachmentConstraint;
    /**
     * CPU density constraint solver (SPH-style for fluid-PBD coupling).
     * Each fluid particle corrects toward rest density using neighbor kernel.
     *
     * Position correction (Müller PBD-fluid):
     *   Δp_i = (1/ρ₀) Σ_j (λ_i + λ_j) ∇_i W(p_i − p_j)
     *
     * densityNeighbors[p] stores PARTICLE indices (indices into densityParticles
     * and the lambdas array), not raw vertex indices.  This allows O(1) lookup
     * of lambdaJ.
     */
    private solveDensityConstraints;
    /**
     * Configure density constraints for fluid particles in the unified buffer.
     * Call this to enable fluid-PBD coupling.
     */
    setDensityParticles(particleIndices: number[], neighbors: number[][], restDensity?: number, kernelRadius?: number, compliance?: number): void;
    private solveSelfCollision;
    /**
     * Update the target position of an existing attachment constraint.
     * Used for dynamic targets like grab interaction or rigid body following.
     */
    updateAttachmentTarget(vertexIndex: number, newTarget: IVector3): void;
    private solveSdfCollision;
    private recomputeNormals;
    /** Get current state (positions, normals, etc.) */
    getState(): ISoftBodyState;
    /** Pause/resume */
    setActive(active: boolean): void;
    /** Reset to rest shape */
    reset(): void;
}
/**
 * GPU-accelerated PBD solver.
 * Uses WebGPU compute shaders for massively parallel constraint solving.
 */
export declare class PBDSolverGPU {
    private config;
    private state;
    private device;
    private buffers;
    private pipelines;
    private bindGroups;
    private coloring;
    constructor(config: ISoftBodyConfig, device: GPUDevice);
    /**
     * Initialize GPU resources (buffers, pipelines, bind groups).
     */
    initialize(): Promise<void>;
    private createBuffer;
    private createPipeline;
    /**
     * Step the simulation.
     */
    step(dt: number): ISoftBodyState;
    getState(): ISoftBodyState;
}
/**
 * Create a soft-body solver with the given configuration.
 * Applies preset parameters if a preset is specified.
 */
export declare function createPBDSolver(config: ISoftBodyConfig): Promise<PBDSolverCPU | PBDSolverGPU>;
//# sourceMappingURL=PBDSolver.d.ts.map