Probleme bei Real Time Fluids

[F34r0fTh3D4rk]

hi, ich hänge schon lange an diesem Problem und komme einfach nicht weiter.

ich habe einen C Code nach Delphi portiert, dabei ging es um Echtzeit Fluid Berechnungen.

Der Solver:
jgt.akpeters.com/pap...Stam01/solver03.html

Die Demo:
jgt.akpeters.com/papers/Stam01/demo03.html

Die Doku:
www.dgp.toronto.edu/...search/pdf/GDC03.pdf

Das Problem war, dass mein Velocity Grid, bei einer Viskosität von 0 nach einiger Zeit völlig verrückt spielte und die Dichte nicht korrekt dargestellt wurde (es bildeten sich merkwürdige Streifen). Ich denke, dass diese beiden Probleme auf die gleiche Ursache zurückzuführen ist.

Ich habe den Source also komplett neu übersetzt, dieses mal schick als Klasse damit ich diesen gleich verwenden konnte, jedoch traten die gleichen Probleme auf. Nun bin ich am Ende mit meinem Latein und brauche Hilfe, hier mein Source:

Delphi-Quellcode:

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.

crossposting DF: http://www.delphi-forum.de/viewtopic.php?t=65342

Vielen Dank schonmal,

mfg