unit fluid_solver;
interface
uses
dglOpengl;
type
TSingleArray =
array of Single;
TFluidSimulation =
class
private
fN,
fsize: integer;
fu, fv, fu_prev, fv_prev,
fdens, fdens_prev: TSingleArray;
fdiff, fvisc: single;
procedure SetN(
const Value: integer);
function IX(i, j: integer): integer;
procedure SWAP(
var x0, x: TSingleArray);
procedure add_source(
var x: TSingleArray; s: TSingleArray; dt: single);
procedure set_bnd(b: integer;
var x: TSingleArray);
procedure lin_solve(b: integer;
var x: TSingleArray; x0: TSingleArray; a, c: single);
procedure diffuse(b: integer;
var x: TSingleArray; x0: TSingleArray; diff, dt: single);
procedure advect(b: integer;
var d: TSingleArray; d0, u, v: TSingleArray; dt: single);
procedure project(
var u, v, p, d: TSingleArray);
procedure dens_step(
var x, x0: TSingleArray; u, v: TSingleArray; diff, dt: single);
procedure vel_step(
var u, v, u0, v0: TSingleArray; visc, dt: single);
public
procedure Initialize(pN: integer; pvisc, pdiff: single);
procedure Reset;
procedure NextStep(pTimestep: single);
procedure AddForce(px, py: integer; pForceX, pForceY: single);
procedure AddSource(px, py: integer; pSource: single);
procedure Draw_Velocity;
procedure Draw_Density;
end;
implementation
procedure TFluidSimulation.Initialize(pN: integer; pvisc, pdiff: single);
begin
SetN(pN);
setlength(fu, fsize);
setlength(fv, fsize);
setlength(fu_prev, fsize);
setlength(fv_prev, fsize);
setlength(fdens, fsize);
setlength(fdens_prev, fsize);
fvisc := pvisc;
fdiff := pdiff;
end;
procedure TFluidSimulation.SetN(
const Value: integer);
begin
fN := value;
fsize := (fN + 2) * (fN + 2);
end;
function TFluidSimulation.IX(i, j: integer): integer;
begin
result := i + (fN + 2) * j;
end;
procedure TFluidSimulation.SWAP(
var x0, x: TSingleArray);
var
tmp: TSingleArray;
begin
tmp := x0;
x0 := x;
x := tmp;
end;
procedure TFluidSimulation.add_source(
var x: TSingleArray; s: TSingleArray; dt: single);
var
i: integer;
begin
for i := 0
to fsize - 1
do
x[i] := x[i] + dt * s[i];
end;
procedure TFluidSimulation.set_bnd(b: integer;
var x: TSingleArray);
var
i: integer;
begin
for i := 1
to fN - 1
do
begin
if b = 1
then
begin
x[IX(0 , i)] := -x[IX(1 , i)];
x[IX(fN + 1, i)] := -x[IX(fN, i)];
end else
begin
x[IX(0 , i)] := x[IX(1 , i)];
x[IX(fN + 1, i)] := x[IX(fN, i)];
end;
if b = 2
then
begin
x[IX(i, 0 )] := -x[IX(i, 1 )];
x[IX(i, fN + 1)] := -x[IX(i, fN)];
end else
begin
x[IX(i, 0 )] := x[IX(i, 1 )];
x[IX(i, fN + 1)] := x[IX(i, fN)];
end;
end;
x[IX(0 , 0 )] := 0.5 * (x[IX(1 , 0 )] + x[IX(0 , 1 )]);
x[IX(0 , fN + 1)] := 0.5 * (x[IX(1 , fN + 1)] + x[IX(0 , fN)]);
x[IX(fN + 1, 0 )] := 0.5 * (x[IX(fN, 0 )] + x[IX(fN + 1, 1 )]);
x[IX(fN + 1, fN + 1)] := 0.5 * (x[IX(fN, fN + 1)] + x[IX(fN + 1, fN)]);
end;
procedure TFluidSimulation.lin_solve(b: integer;
var x: TSingleArray; x0: TSingleArray; a, c: single);
var
i, j, k: integer;
begin
for k := 0
to 19
do
begin
for i := 1
to fN
do
for j := 1
to fN
do
x[IX(i, j)] := (x0[IX(i, j)] + a * (x[IX(i - 1, j)] + x[IX(i + 1, j)] + x[IX(i, j - 1)] + x[IX(i, j + 1)])) / c;
set_bnd(b, x);
end;
end;
procedure TFluidSimulation.diffuse(b: integer;
var x: TSingleArray; x0: TSingleArray; diff, dt: single);
var
a: single;
begin
a := dt * diff * fN * fN;
lin_solve(b, x, x0, a, 1 + 4 * a);
end;
procedure TFluidSimulation.advect(b: integer;
var d: TSingleArray; d0, u, v: TSingleArray; dt: single);
var
i, j, i0, j0, i1, j1: integer;
x, y, s0, t0, s1, t1, dt0: single;
begin
dt0 := dt * fN;
for i := 1
to fN
do
for j := 1
to fN
do
begin
x := i - dt0 * u[IX(i, j)];
y := j - dt0 * v[IX(i, j)];
if (x < 0.5)
then
x := 0.5;
if (x > fN + 0.5)
then
x := fN + 0.5;
i0 := round(x);
i1 := i0 + 1;
if (y < 0.5)
then
y := 0.5;
if (y > fN + 0.5)
then
y := fN + 0.5;
j0 := round(y);
j1 := j0 + 1;
s1 := x - i0;
s0 := 1 - s1;
t1 := y - j0;
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)]);
end;
set_bnd(b, d);
end;
procedure TFluidSimulation.project(
var u, v, p, d: TSingleArray);
var
i, j: integer;
begin
for i := 1
to fN
do
for j := 1
to fN
do
begin
d[IX(i, j)] := -0.5 * (u[IX(i + 1, j)] - u[IX(i - 1, j)] + v[IX(i, j + 1)] - v[IX(i, j - 1)]) / fN;
p[IX(i, j)] := 0;
end;
set_bnd(0, d);
set_bnd(0, p);
lin_solve(0, p, d, 1, 4);
for i := 1
to fN
do
for j := 1
to fN
do
begin
u[IX(i,j)] := u[IX(i,j)] - (0.5 * fN * (p[IX(i + 1, j)] - p[IX(i - 1, j)]));
v[IX(i,j)] := v[IX(i,j)] - (0.5 * fN * (p[IX(i, j + 1)] - p[IX(i, j - 1)]));
end;
set_bnd(1, u);
set_bnd(2, v);
end;
procedure TFluidSimulation.dens_step(
var x, x0: TSingleArray; u, v: TSingleArray; diff, dt: single);
begin
add_source(x, x0, dt);
SWAP(x0, x);
diffuse(0, x, x0, diff, dt);
SWAP(x0, x);
advect(0, x, x0, u, v, dt);
end;
procedure TFluidSimulation.vel_step(
var u, v, u0, v0: TSingleArray; visc, dt: single);
begin
add_source(u, u0, dt);
add_source(v, v0, dt);
SWAP(u0, u);
diffuse(1, u, u0, visc, dt);
SWAP(v0, v);
diffuse(2, v, v0, visc, dt);
project(u, v, u0, v0);
SWAP(u0, u);
SWAP(v0, v);
advect(1, u, u0, u0, v0, dt);
advect(2, v, v0, u0, v0, dt);
project(u, v, u0, v0);
end;
procedure TFluidSimulation.NextStep(pTimestep: single);
begin
vel_step(fu, fv, fu_prev, fv_prev, fvisc, pTimestep);
dens_step(fdens, fdens_prev, fu, fv, fdiff, pTimestep);
end;
procedure TFluidSimulation.Reset;
var
i: integer;
begin
for i := 0
to fsize - 1
do
begin
fu[i] := 0;
fv[i] := 0;
fu_prev[i] := 0;
fv_prev[i] := 0;
fdens[i] := 0;
fdens_prev[i] := 0;
end;
end;
procedure TFluidSimulation.AddForce(px, py: integer; pForceX, pForceY: single);
begin
fu_prev[IX(px, py)] := pForceX;
fv_prev[IX(px, py)] := pForceY;
end;
procedure TFluidSimulation.AddSource(px, py: integer; pSource: single);
begin
fdens_prev[IX(px, py)] := pSource;
end;
procedure TFluidSimulation.Draw_Velocity;
var
i, j: integer;
x, y, h: single;
begin
h := 1 / fN;
glColor3f(1, 1, 1);
glLineWidth(1);
glBegin(GL_LINES);
for i := 1
to fN
do
begin
x := (i - 0.5) * h;
for j := 1
to fN
do
begin
y := (j - 0.5) * h;
glColor3f(1, 1, 1);
glVertex2f(x, y);
glVertex2f(x + fu[IX(i, j)], y + fv[IX(i, j)]);
end
end;
glEnd;
end;
procedure TFluidSimulation.Draw_Density;
var
i, j: integer;
x, y, h, d00, d01, d10, d11: single;
begin
h := 1 / fN;
glBegin(GL_QUADS);
for i := 0
to fN
do
begin
x := (i - 0.5) * h;
for j := 0
to fN
do
begin
y := (j - 0.5) * h;
d00 := fdens[IX(i , j )];
d01 := fdens[IX(i , j + 1)];
d10 := fdens[IX(i + 1, j )];
d11 := fdens[IX(i + 1, j + 1)];
glColor3f(d00, 0, 0); glVertex2f(x , y );
glColor3f(d10, 0, 0); glVertex2f(x + h, y );
glColor3f(d11, 0, 0); glVertex2f(x + h, y + h);
glColor3f(d01, 0, 0); glVertex2f(x , y + h);
end;
end;
glEnd;
end;
end.