做物理不懂流体，就好比打 dota 不玩卡尔，缺少一种品味的美。于是花了几天时间补了一下，虽然玩的烂但是好歹能放出民工三连的水平。
读完以后感觉粒子这个东西是真的好用啊。凡是解析方法要搞一堆积分的东西，换成粒子瞬间就简单好多，梯度也好算了。
已经有非常多介绍 sph 算法推导过程的文章了，最原始的论文链接也放在这，就只简单讲下算法的思路了，网页上敲公式也是累死人。
Matthias Muller, SIGGRAPH 2003
简而言之，原先你要模拟流体，就需要求解大名鼎鼎的 Navier-Stokes 方程，长这样：
这是一个 force-based 模型，右侧三项分别是 ,,
,
。左侧是动量对 t, x, y, z 的全导数拆解，等价于动量的变化。简单来说就是
的流体版。
现在我们换粒子解决这个方程的求解，简单方便。
模拟上存在两个视角，Euler 和 Lagrangian。Euler methods 求解连续的场，例如速度场 v(x,y,z,t)。而 Lagrangian 方法求解离散的粒子。这篇文章比较清楚地总结了两种方法的优缺点。SPH_2017 对于一片规则的水域，可以采用网格法求解场，但是对于经常需要产生碰撞、形变的水体，粒子方法比网格法有着明显优势。
SPH 全称是 Smoothed Particle Hydrodynamics，Smoothed 是一个重要的特性，这使得每个粒子的特性都可以由临近的粒子计算得出。从而产生了邻域核函数的应用。N-S 方程右侧的两个梯度都可以用核函数的方式进行重写，从而得到简单的求梯度公式。论文中列的很清晰，这里不再赘述。
既然要 in a weekend 那么我们就撸起袖子开干，发挥面向 google 编程的优秀传统。
在 github 上看了几个 sph 相关的代码库。SplishSplash 是工业化做的比较好的，但是不能满足我们 in a weekend 快速阅读实现的需求。所以又找了一个 [SPH-3D-Fluid-Simulation] (https://github.com/saeedmahani/SPH-3D-Fluid-Simulation)。这个项目写的很 easy，课程作业的难度。ok 就是它了，拿来改改。Visual studio， 启动！
能运行，效果还凑合，那我们直接进行一个代码的读。
首先是基础流程，所有物理模拟都是类似这个。这里建了一个网格主要用来做邻域的圈定，每个粒子只需要计算邻域粒子的影响。
PARTICLE_SYSTEM::stepVerlet(double dt)
calculateAcceleration() // 更新所有粒子的加速度
//用加速度更新速度 位置 (Verlet)
VEC3D newPosition = particle.position() + particle.velocity()*dt + particle.acceleration()*dt*dt;
VEC3D newVelocity = (newPosition - particle.position()) / dt;
updateGrid() //粒子移动后更新所有网格中保存的粒子
求 density , pressure，参考 Matthias 论文 式(3) ， （12）
for particle 的所有邻域粒子 neighborParticle:
VEC3D diffPosition = particle.position() - neighborParticle.position();
double radiusSquared = diffPosition.dot(diffPosition);
particle.density() += Wpoly6(radiusSquared);
particle.density() *= PARTICLE_MASS;
particle.pressure() = GAS_STIFFNESS * (particle.density() - REST_DENSITY);
各种核函数及其梯度：
Wpoly6, Wpoly6Gradient, Wpoly6Laplacian : Matthias 论文 式（20）
WspikyGradient : Matthias 论文 式（21）
WviscosityLaplacian: Matthias 论文 式（22）
然后计算 N-S 右侧。对于 p0 粒子有：
for p0 的所有邻域粒子 p1:
f_pressure += p0.pressure/pow(p0.density,2) +
p1.pressure/pow(p1.density,2) * spikyGradient;
f_viscosity += (p1.v - p0.v) * WviscosityLaplacian(rsqd) / p1.density;
对于表面张力要特殊处理。计算表面粒子与表面张力如下。
for p0 的所有邻域粒子 p1:
colorFieldNormal += poly6Gradient / p1.density;
colorFieldLaplacian += Wpoly6Laplacian(rsqd) / p1.density;
colorFieldNormalMagnitude = colorFieldNormal.magnitude();
if (colorFieldNormalMagnitude > SURFACE_THRESHOLD) {
particle.flag() = true;
f_surface = -SURFACE_TENSION * colorFieldNormal / colorFieldNormalMagnitude * colorFieldLaplacian;
汇总计算加速度
f_gravity(0.0, particle.density() * -9.80665, 0.0);
particle.acceleration() =
(f_pressure + f_viscosity + f_surface + f_gravity) / particle.density();
上面跑完就是一个 SPH 的完整流程了。公式很复杂，但流程其实很简单。
至于渲染，其实粒子可以视为一种点云，而点云在点的数量比较少的情况下，可以通过 marching cubes 方法通过查找等值面来将点云重建为 mesh。然后就可以走正常的渲染了。这也是很经典的方法了，至于有没有其他的特技，就交给更专业的人吧。
SPH in a weekend，跑通收工。
