1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
242:
243:
| unit UMain;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, ExtCtrls;
const
N = 10;
size = (N+2)*(N+2);
type
TSingleArray = array [1..size-1] of Single;
var
u, v, u_prev, v_prev,
dens, dens_prev: TSingleArray;
visc, diff,
dt: single;
type
TMainForm = class(TForm)
IMG: TImage;
procedure FormCreate(Sender: TObject);
procedure ApplicationEventsIdle(Sender: TObject; var Done: Boolean);
procedure draw_dens(N: integer; dens: TSingleArray);
private
CalcTime,
Frequency: int64;
public
end;
var
MainForm: TMainForm;
implementation
{$R *.dfm}
function IX(i, j: integer): integer;
begin
result := ((i)+(N+2)*(j));
end;
procedure SWAP(x0, x: TSingleArray);
var
tmp: TSingleArray;
begin
tmp := x0;
x0 := x;
x := tmp;
end;
procedure set_bnd(N, b: integer; x: TSingleArray);
var
i: integer;
begin
for i := 1 to N do
begin
if b = 1 then
begin
x[IX(0, i)] := -x[IX(1, i)];
x[IX(N+1,i)] := -x[IX(N, i)];
end else
begin
x[IX(0, i)] := x[IX(1, i)];
x[IX(N+1,i)] := x[IX(N, i)];
end;
if b = 2 then
begin
x[IX(i, 0)] := -x[IX(i, 1)];
x[IX(i, N+1)] := -x[IX(i, N)];
end else
begin
x[IX(i, 0)] := x[IX(i,1 )];
x[IX(i, N+1)] := x[IX(i, N)];
end;
end;
x[IX(0, 0)] := 0.5*(x[IX(1, 0)]+x[IX(0, 1)]);
x[IX(0, N+1)] := 0.5*(x[IX(1,N+1)]+x[IX(0 ,N)]);
x[IX(N+1, 0)] := 0.5*(x[IX(N, 0)]+x[IX(N+1, 1)]);
x[IX(N+1, N+1)] := 0.5*(x[IX(N, N+1)]+x[IX(N+1, N)]);
end;
procedure project(N: integer; u, v, p, divi: TSingleArray);
var
i, j, k: integer;
h: single;
begin
h := 1/N;
for i := 1 to N do
for j := 1 to N do
begin
divi[IX(i, j)] := -0.5*h*(u[IX(i+1, j)]-u[IX(i-1, j)]+v[IX(i, j+1)]-v[IX(i, j-1)]);
p[IX(i, j)] := 0;
end;
set_bnd(N, 0, divi);
set_bnd(N, 0, p);
for k := 0 to 19 do
for i := 1 to N do
for j := 1 to N do
p[IX(i, j)] := (divi[IX(i, j)]+p[IX(i-1, j)]+p[IX(i+1, j)]+p[IX(i, j-1)]+p[IX(i, j+1)])/4;
set_bnd(N, 0, p);
for i := 1 to N do
for j := 1 to N do
begin
u[IX(i, j)] := u[IX(i, j)] - (0.5*(p[IX(i+1, j)]-p[IX(i-1, j)])/h);
v[IX(i, j)] := v[IX(i, j)] - (0.5*(p[IX(i, j+1)]-p[IX(i, j-1)])/h)
end;
set_bnd(N, 1, u);
set_bnd(N, 2, v);
end;
procedure add_source(N: integer; x, s: TSingleArray; dt: single);
var
i: integer;
begin
for i := 1 to size-1 do
x[i] := x[i]+dt*s[i];
end;
procedure diffuse(N, b: integer; x, x0: TSingleArray; diff, dt: single);
var
i, j, k: integer;
a: single;
begin
a := dt*diff*N*N;
for k := 0 to 19 do
for i := 1 to N do
for j := 1 to N 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)])) / (1+4*a);
set_bnd(N, b, x);
end;
procedure advect(N, b: integer; d, 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*N;
for i := 1 to N do
for j := 1 to N 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 > N + 0.5) then
x := N + 0.5;
i0 := round(x);
i1 := i0 + 1;
if (y < 0.5) then
y := 0.5;
if (y > N + 0.5) then
y := N + 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(N, b, d);
end;
procedure dens_step(N: integer; x, x0, u, v: TSingleArray; diff, dt: single);
begin
add_source(N, x, x0, dt);
SWAP(x0, x);
diffuse(N, 0, x, x0, diff, dt);
SWAP(x0, x);
advect(N, 0, x, x0, u, v, dt);
end;
procedure vel_step(N: integer; u, v, u0, v0: TSingleArray; visc, dt: single);
begin
add_source(N, u, u0, dt);
add_source(N, v, v0, dt);
SWAP(u0, u);
diffuse(N, 1, u, u0, visc, dt);
SWAP(v0, v);
diffuse(N, 2, v, v0, visc, dt);
project(N, u, v, u0, v0);
SWAP(u0, u);
SWAP(v0, v);
advect(N, 1, u, u0, u0, v0, dt);
advect(N, 2, v, v0, u0, v0, dt);
project(N, u, v, u0, v0);
end;
procedure TMainForm.draw_dens(N: integer; dens: TSingleArray);
var
i, j,
val: integer;
begin
for i := 1 to N do
for j := 1 to N do
begin
val := round(dens[IX(i, j)]);
IMG.Canvas.Pixels[i, j] := RGB(val, val, val);
end;
end;
procedure TMainForm.ApplicationEventsIdle(Sender: TObject; var Done: Boolean);
var
CurrTime: int64;
begin
QueryPerformanceCounter(Currtime);
dt := (CurrTime - CalcTime) / Frequency;
CalcTime := CurrTime;
vel_step(N, u, v, u_prev, v_prev, visc, dt);
dens_step(N, dens, dens_prev, u, v, diff, dt);
draw_dens(N, dens);
done := false;
end;
procedure TMainForm.FormCreate(Sender: TObject);
var
i, j: integer;
begin
QueryPerformanceFrequency(Frequency);
QueryPerformanceCounter(CalcTime);
Application.OnIdle := ApplicationEventsIdle;
for i := 1 to N do
for j := 1 to N do
begin
dens_prev[IX(i, j)] := 0;
u_prev[IX(i, j)] := random(5);
v_prev[IX(i, j)] := random(5);
end;
end;
Initialization
randomize;
end. |
