Границы моделирования жидкости и адвект

Это моделирование жидкости основано на статья Стама. На странице 7 он описывает основную идею адвекции:

Начните с двух сеток: одна содержит значения плотности из предыдущего временного шага, а вторая —
это будет содержать новые значения. Для каждой ячейки сетки мы прослеживаем центр ячейки
Положение назад через поле скоростей. Затем мы линейно интерполируем из сетки
предыдущие значения плотности и присвойте это значение текущей ячейке сетки.

Advect код. Две сетки плотности 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 вектор bools со строками и столбцами 0 в N+1, Граничные объекты устанавливаются перед основным циклом путем установки координат ячейки в bound в 1,

В документе смутно говорится: «Тогда мы просто должны добавить
некоторый код set_bnd() рутина для заполнения значений для занятых ячеек из значений
их непосредственные соседи «, что-то вроде того, что я делаю. Я ищу способ более точно реализовать границы, который имеет нежидкие твердые границы и, возможно, в конечном итоге поддерживает границы для нескольких жидкостей. Качество изображения гораздо важнее чем физика правильности.