这个流体模拟基于Stam的论文。他在第7页介绍了对流背后的基本思想：
这种方法简洁有效，但是由于值是向后追踪和插值的，实现对象边界对我来说很棘手。我目前的解决方案是，如果旁边有空白区域（或空格），则将密度推出边界，但这只是在视觉上准确，而且会导致密度在角落和具有对角线速度的区域积累。我现在正在寻求“正确性”。
边界代码的相关部分：
下面是代码。两个密度网格为从两个网格开始：一个包含上一时间步长的密度值，另一个将包含新值。对于后者的每个网格单元，我们通过速度场向后跟踪单元的中心位置。然后我们从以前的密度值网格进行线性插值，并将该值分配给当前网格单元。
d和d0，u和v是速度分量，dt是时间步长，N（全局）是网格大小，b可以忽略：void advect(int b, vfloat &d, const vfloat &d0, const vfloat &u, const vfloat &v, float dt, std::vector<bool> &bound)
{
float dt0 = dt*N;
for (int i=1; i<=N; i++)
{
for (int j=1; j<=N; j++)
{
float x = i - dt0*u[IX(i,j)];
float y = j - dt0*v[IX(i,j)];
if (x<0.5) x=0.5; if (x>N+0.5) x=N+0.5;
int i0=(int)x; int i1=i0+1;
if (y<0.5) y=0.5; if (y>N+0.5) y=N+0.5;
int j0=(int)y; int j1=j0+1;
float s1 = x-i0; float s0 = 1-s1; float t1 = y-j0; float t0 = 1-t1;
d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)] + t1*d0[IX(i0,j1)]) +
s1*(t0*d0[IX(i1,j0)] + t1*d0[IX(i1,j1)]);
}
}
set_bnd(b, d, bound);
}
这种方法简洁有效，但是由于值是向后追踪和插值的，实现对象边界对我来说很棘手。我目前的解决方案是，如果旁边有空白区域（或空格），则将密度推出边界，但这只是在视觉上准确，而且会导致密度在角落和具有对角线速度的区域积累。我现在正在寻求“正确性”。
边界代码的相关部分：
void set_bnd(const int b, vfloat &x, std::vector<bool> &bound)
{
//...
for (int i=1; i<=N; i++)
{
for (int j=1; j<=N; j++)
{
if (bound[IX(i,j)])
{
//...
else if (b==0)
{
// Distribute density from bound to surrounding cells
int nearby_count = !bound[IX(i+1,j)] + !bound[IX(i-1,j)] + !bound[IX(i,j+1)] + !bound[IX(i,j-1)];
if (!nearby_count) x[IX(i,j)] = 0;
else
x[IX(i,j)] = ((bound[IX(i+1,j)] ? 0 : x[IX(i+1,j)]) +
(bound[IX(i-1,j)] ? 0 : x[IX(i-1,j)]) +
(bound[IX(i,j+1)] ? 0 : x[IX(i,j+1)]) +
(bound[IX(i,j-1)] ? 0 : x[IX(i,j-1)])) / surround;
}
}
}
}
}
bound 是一个布尔型向量，其行和列范围为 0 到 N+1。在主循环之前通过将单元格坐标设置为 1 来设置边界对象。
论文模糊地说明“然后我们只需添加一些代码到 set_bnd() 例程中，从它们直接相邻的值来填写占用单元格的值”，这或多或少是我所做的。 我正在寻找更精确实现边界的方法，即具有非流体固体边界，并最终支持多流体边界。视觉质量比物理正确性更重要。
set_bnd的注释部分处理的是不需要回答的边缘情况。我的完整代码在这里。 - qwrforall bpt in bound: bpt.i,bpt.j,bpt.xxx(例如预计算更多的东西，所以循环中就少了if)。为了看得更清楚，我在渲染循环的底部做了网格:SDL_SetRenderDrawColor(renderer,100,100,0,0); SDL_RenderDrawRect(renderer,&r);但由于四舍五入，有些线条是双倍宽的。 - Craig Estey