Симуляція рідини SPH з нуля у 200 рядках

// Радіус згладжування const H = 30; // px — усі сусіди мають бути в межах цієї відстані const H2 = H * H; // Ядро Poly6 — для оцінки щільності (гладке навіть при r=0) function Wpoly6(r2) { if (r2 >= H2) return 0; const x = H2 - r2; return 4 / (Math.PI * H2 * H2 * H2 * H2) * x * x * x; // точний коефіцієнт: 4/(pi * h^8), але ми вкладаємо його у вибір REST_DENSITY } // Градієнт ядра Spiky — для сил ТИСКУ (ненульовий при r=0, протидіє злипанню) function gradWspiky(rx, ry, r) { if (r <= 0 || r >= H) return { x: 0, y: 0 }; const coeff = -30 / (Math.PI * H2 * H2 * H2) * (H - r) * (H - r) / r; return { x: coeff * rx, y: coeff * ry }; } // Лапласіан ядра Viscosity — для сил В'ЯЗКОСТІ (додатний, розсіює швидкість) function lapWviscosity(r) { if (r >= H) return 0; return 40 / (Math.PI * H2 * H2 * H2) * (H - r); }

const REST_DENSITY = 300; // цільова щільність (приблизно кг/м², у одиницях частинок) const MASS = 1; // маса частинки function computeDensities(particles) { for (const pi of particles) { let density = 0; for (const pj of particles) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r2 = dx*dx + dy*dy; density += MASS * Wpoly6(r2); } pi.density = Math.max(density, REST_DENSITY * 0.01); // уникаємо нуля } }

const STIFFNESS = 200; // тиск = жорсткість * (щільність - щільність_спокою) function computePressure(p) { p.pressure = STIFFNESS * (p.density - REST_DENSITY); } function pressureForce(pi, pj) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r = Math.sqrt(dx*dx + dy*dy); if (r === 0 || r >= H) return { x: 0, y: 0 }; // Симетрична сила тиску (3-й закон Ньютона) const pressure = (pi.pressure + pj.pressure) / 2; const grad = gradWspiky(dx, dy, r); return { x: -MASS * pressure / pj.density * grad.x, y: -MASS * pressure / pj.density * grad.y }; }

const VISCOSITY = 200; function viscosityForce(pi, pj) { const dx = pi.x - pj.x, dy = pi.y - pj.y; const r = Math.sqrt(dx*dx + dy*dy); if (r >= H) return { x: 0, y: 0 }; const lap = lapWviscosity(r); const dvx = pj.vx - pi.vx; const dvy = pj.vy - pi.vy; return { x: VISCOSITY * MASS / pj.density * dvx * lap, y: VISCOSITY * MASS / pj.density * dvy * lap }; } // Поєднуємо всі сили function computeForces(particles, gravity = 500) { for (const pi of particles) { pi.fx = 0; pi.fy = 0; computePressure(pi); } for (let i = 0; i < particles.length; i++) { const pi = particles[i]; for (let j = i + 1; j < particles.length; j++) { const pj = particles[j]; const pf = pressureForce(pi, pj); const vf = viscosityForce(pi, pj); pi.fx -= pf.x + vf.x; pi.fy -= pf.y + vf.y; pj.fx += pf.x + vf.x; pj.fy += pf.y + vf.y; } pi.fy += gravity * pi.density; // гравітація пропорційна щільності } }

const DT = 0.003, DAMPING = -0.5; const W = 800, HEIGHT = 600; // розмір полотна function integrate(particles) { for (const p of particles) { // Напівнеявний (симплектичний) метод Ейлера p.vx += DT * p.fx / p.density; p.vy += DT * p.fy / p.density; p.x += DT * p.vx; p.y += DT * p.vy; // Зіткнення з межами if (p.x < H) { p.x = H; p.vx *= DAMPING; } if (p.x > W - H) { p.x = W - H; p.vx *= DAMPING; } if (p.y < H) { p.y = H; p.vy *= DAMPING; } if (p.y > HEIGHT - H) { p.y = HEIGHT - H; p.vy *= DAMPING; } } } // Рендеринг function render(ctx, particles) { ctx.clearRect(0, 0, W, HEIGHT); ctx.fillStyle = 'rgba(0,10,30,1)'; ctx.fillRect(0, 0, W, HEIGHT); for (const p of particles) { const speed = Math.sqrt(p.vx*p.vx + p.vy*p.vy); const t = Math.min(speed / 300, 1); ctx.fillStyle = `hsl(${200 - t*140}, 80%, ${40 + t*30}%)`; ctx.beginPath(); ctx.arc(p.x, p.y, H * 0.3, 0, Math.PI*2); ctx.fill(); } }

class SpatialHash { constructor(cellSize) { this.cs = cellSize; this.table = new Map(); } key(x, y) { return `${Math.floor(x/this.cs)},${Math.floor(y/this.cs)}`; } clear() { this.table.clear(); } insert(p) { const k = this.key(p.x, p.y); if (!this.table.has(k)) this.table.set(k, []); this.table.get(k).push(p); } neighbours(p, radius) { const result = []; const cells = Math.ceil(radius / this.cs); const cx = Math.floor(p.x / this.cs), cy = Math.floor(p.y / this.cs); for (let dx = -cells; dx <= cells; dx++) for (let dy = -cells; dy <= cells; dy++) { const bucket = this.table.get(`${cx+dx},${cy+dy}`); if (bucket) result.push(...bucket); } return result; } } // Використання у computeForces: замініть внутрішній цикл на // const neighbours = spatialHash.neighbours(pi, H); // for (const pj of neighbours) { ... } // Часова складність: O(N * k), де k = сер. кількість сусідів (~20–50)