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.