http://www.dgp.toronto.edu/people/stam/ ... /GDC03.pdf
daraufhin habe ich den C-Code dort übersetzt nach PB. Das Resultat ist
hier (Maustaste drücken und bewegen):
Code: Alles auswählen
EnableExplicit
Structure FArray
f.f[0]
EndStructure
Macro IX(i, j)
((i) + (n + 2) * (j))
EndMacro
Declare set_bnd(n.l, b.l, *x.FArray)
Declare project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)
Procedure add_source(n.l, *x.FArray, *s.FArray, dt.f)
Define.l i, Size = (n + 2) * (n + 2)
For i = 0 To Size - 1
*x\f[i] + dt * *s\f[i]
Next
EndProcedure
Procedure diffuse(n.l, b.l, *x.FArray, *x0.FArray, diff.f, dt.f)
Define.l i, j, k
Define.f a
a = dt * diff * n * n
For k = 0 To 19
For i = 1 To n
For j = 1 To n
*x\f[IX(i, j)] = (*x0\f[IX(i, j)] + a * (*x\f[IX(i-1, j)] + *x\f[IX(i+1, j)] + *x\f[IX(i, j-1)] + *x\f[IX(i, j+1)])) / (1.0 + 4.0 * a)
;*x\f[IX(i, j)] = *x0\f[IX(i, j)] + a * (*x0\f[IX(i-1, j)] + *x0\f[IX(i+1, j)] + *x0\f[IX(i, j-1)] + *x0\f[IX(i, j+1)] - 4 * *x0\f[IX(i, j)])
Next
Next
set_bnd(n, b, *x)
Next
EndProcedure
Procedure advect (n.l, b.l, *d.FArray, *d0.FArray, *u.FArray, *v.FArray, dt.f)
Define.l i, j, i0, j0, i1, j1
Define.f x, y, s0, t0, s1, t1, dt0
dt0 = dt * n
For i = 1 To n
For j = 1 To n
x = i - dt0 * *u\f[IX(i, j)]
y = j - dt0 * *v\f[IX(i, j)]
If (x < 0.5)
x = 0.5
EndIf
If (x > n + 0.5)
x = n + 0.5
EndIf
i0 = Int(x)
i1 = i0 + 1
If (y < 0.5)
y = 0.5
EndIf
If (y > n + 0.5)
y = n + 0.5
EndIf
j0 = Int(y)
j1 = j0 + 1
s1 = x - i0 : s0 = 1 - s1 : t1 = y - j0 : t0 = 1 - t1
*d\f[IX(i, j)] = s0 * (t0 * *d0\f[IX(i0, j0)] + t1 * *d0\f[IX(i0, j1)]) + s1 * (t0 * *d0\f[IX(i1, j0)] + t1 * *d0\f[IX(i1, j1)])
Next
Next
set_bnd(n, b, *d)
EndProcedure
Procedure dens_step(n.l, *x.FArray, *x0.FArray, *u.FArray, *v.FArray, diff.f, dt.f)
add_source(n, *x, *x0, dt)
diffuse(n, 0, *x0, *x, diff, dt )
advect(n, 0, *x, *x0, *u, *v, dt )
EndProcedure
Procedure vel_step(n.l, *u.FArray, *v.FArray, *u0.FArray, *v0.FArray, visc.f, dt.f)
add_source(n, *u, *u0, dt) : add_source(n, *v, *v0, dt)
diffuse(n, 1, *u0, *u, visc, dt)
diffuse(n, 2, *v0, *v, visc, dt)
project(n, *u0, *v0, *u, *v)
advect(n, 1, *u, *u0, *u0, *v0, dt) : advect(n, 2, *v, *v0, *u0, *v0, dt)
project(n, *u, *v, *u0, *v0)
EndProcedure
Procedure project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)
Define.l i, j, k
Define.f h = 1.0 / n
For i = 1 To n
For j = 1 To n
*div\f[IX(i, j)] = -0.5 * h * (*u\f[IX(i+1, j)] - *u\f[IX(i-1, j)] + *v\f[IX(i, j+1)] - *v\f[IX(i, j-1)])
*p\f[IX(i, j)] = 0
Next
Next
set_bnd(n, 0, *div) : set_bnd(n, 0, *p)
For k = 0 To 19
For i = 1 To n
For j = 1 To n
*p\f[IX(i, j)] = (*div\f[IX(i, j)] + *p\f[IX(i-1, j)] + *p\f[IX(i+1, j)] + *p\f[IX(i, j-1)] + *p\f[IX(i, j+1)]) / 4
Next
Next
set_bnd(n, 0, *p)
Next
For i = 1 To n
For j = 1 To n
*u\f[IX(i, j)] - (0.5 * (*p\f[IX(i+1, j)] - *p\f[IX(i-1, j)]) / h)
*v\f[IX(i, j)] - (0.5 * (*p\f[IX(i, j+1)] - *p\f[IX(i, j-1)]) / h)
Next
Next
set_bnd(n, 1, *u) : set_bnd(n, 2, *v)
EndProcedure
Procedure set_bnd(n.l, b.l, *x.FArray)
Define.l i
For i = 1 To n
If b = 1
*x\f[IX(0, i)] = -*x\f[IX(1, i)]
*x\f[IX(n+1, i)] = -*x\f[IX(n, i)]
Else
*x\f[IX(0, i)] = *x\f[IX(1, i)]
*x\f[IX(n+1, i)] = *x\f[IX(n, i)]
EndIf
If b = 2
*x\f[IX(i, 0)] = -*x\f[IX(i, 1)]
*x\f[IX(i, n+1)] = -*x\f[IX(i, n)]
Else
*x\f[IX(i, 0)] = *x\f[IX(i, 1)]
*x\f[IX(i, n+1)] = *x\f[IX(i, n)]
EndIf
Next
*x\f[IX(0, 0)] = 0.5 * (*x\f[IX(1, 0)] + *x\f[IX(0, 1)])
*x\f[IX(0, n+1)] = 0.5 * (*x\f[IX(1, n+1)] + *x\f[IX(0, n)])
*x\f[IX(n+1, 0)] = 0.5 * (*x\f[IX(n, 0)] + *x\f[IX(n+1, 1)])
*x\f[IX(n+1, n+1)] = 0.5 * (*x\f[IX(n, n+1)] + *x\f[IX(n+1, n)])
EndProcedure
;- End Simulation
#N = 60
#SIZE = (#N + 2) * (#N + 2)
Procedure array_normalize(n.l, *a.FArray)
Define.l i, j, Size = (n + 2) * (n + 2)
Define.f max
For i = 0 To n - 1
For j = 0 To n - 1
If Abs(*a\f[IX(i, j)]) > max
max = Abs(*a\f[IX(i, j)])
EndIf
Next
Next
For i = 0 To Size - 1
*a\f[i] / (max + 0.001)
Next
EndProcedure
Procedure draw_dens(n.l, *dens.FArray)
Define.l i, j
Define.f max
ClearScreen(0)
ExamineKeyboard()
If KeyboardPushed(#PB_Key_Escape)
End
EndIf
StartDrawing(ScreenOutput())
array_normalize(n, *dens)
For i = 0 To n - 1
For j = 0 To n - 1
Box(i * 10, j * 10, 10, 10, Abs(*dens\f[IX(i, j)]) * 255)
Next
Next
StopDrawing()
WindowEvent()
FlipBuffers()
EndProcedure
Procedure set_datas(n.l, *dens0.FArray, *u0.FArray, *v0.FArray)
Define.l i, j
Define.l z, Size = (n + 2) * (n + 2)
ExamineMouse()
i = MouseX() / 10
j = MouseY() / 10
If MouseButton(1)
fillmemory_(*dens0, (n + 2) * (n + 2) * SizeOf(Float), 0)
fillmemory_(*u0, (n + 2) * (n + 2) * SizeOf(Float), 0)
fillmemory_(*v0, (n + 2) * (n + 2) * SizeOf(Float), 0)
*u0\f[IX(i, j)] = 0
*v0\f[IX(i, j)] = -99
;Beep_(1000, 50)
For z = 0 To Size - 1
*dens0\f[z] = 0.0005
Next
EndIf
EndProcedure
Procedure main()
Define *u.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define *v.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define *u_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define *v_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define *dens.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define *dens_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
Define.f visc = 8.0e-4, dt = 0.5, diff = 0.01
Repeat
set_datas(#N, *dens_prev, *u_prev, *v_prev)
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)
Delay(30)
ForEver
EndProcedure
InitSprite()
InitMouse()
InitKeyboard()
OpenWindow(0, 100, 100, #N * 10, #N * 10, "Fluid Dynamics")
OpenWindowedScreen(WindowID(0), 0, 0, #N * 10, #N * 10, 0, 0, 0)
main()
die Funktion array_normalize() ist unnötig, da sie das Resultat verfälscht.
Ev. hat es also noch Fehler oder die Parameter sind nicht im richtigen
Bereich. Wenns jemand gleich sieht, danke für eine Rückmeldung

greetz
Remi