Code: Alles auswählen
) 
; Push right and left mouse button (left: set density, right: set velocity)
EnableExplicit
#N = 100 ; Anzahl der Gitterpunkte, die berechnet werden
#DSTEP = 5 ; Pixel per Gitterpunkt in der Darstellung
#SIZE = (#N + 2) * (#N + 2)
Global plot.l = 0 ; 0 - zeichnet Dichtefeld, sonst Geschwindigkeitsfeld
Global dense.f=0.18, grav.f=9.81 ;dense - Diche des Anfangsfeldes und am Rand, grav - Gravitationskonstante
Global visc.f = 0.0001, dt.f = 0.01, diff.f = 0.00018 ; visc - (turbulente) Vikosität
; dt - Zeitschritt
; diff - (turbulente) Diffusionskonstante
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_gravity(n.l, *x.FArray, *s.FArray, *d.FArray, dt.f)
Define.l i, j
For i = 0 To n+1
For j = 0 To n+1
*x\f[IX(i, j)] + dt* (*s\f[IX(i, j)] - grav * (dense-*d\f[IX(i, j)])/dense)
Next
Next
EndProcedure
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, *d.FArray, visc.f, dt.f)
add_source(n, *u, *u0, dt) : add_source_gravity(n, *v, *v0, *d, 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, 3, *div) : set_bnd(n, 3, *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, 3, *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
Select b
Case 0 ; Dichte Rand
For i = 0 To n+1
*x\f[IX(0, i)] = dense
*x\f[IX(n+1, i)] = dense
*x\f[IX(i, 0)] = dense
*x\f[IX(i, n+1)] = dense
Next
Case 1 ; Ost- bzw. Westrand, fester Rand u=0
For i = 0 To n+1
*x\f[IX(0, i)] = 0.0
*x\f[IX(n+1, i)] = 0.0
Next
Case 2 ; Nord- bzw. Suedrand, fester Rand v=0
For i = 0 To n+1
*x\f[IX(i, 0)] = 0.0
*x\f[IX(i, n+1)] = 0.0
Next
Case 3 ; Raender fuer die Naeherungsverfahren
For i = 0 To n+1
*x\f[IX(0, i)] = 0
*x\f[IX(n+1, i)] = 0
*x\f[IX(i, 0)] = 0
*x\f[IX(i, n+1)] = 0
Next
EndSelect
EndProcedure
;- End Simulation
Procedure.l color_gr( f1.d,f2.d)
Protected r.l, g.l, b.l,v.f
v=$500*Sqr(f1*f1+f2*f2)
g=$FF
b=$80
r = $FF -v
If r < 0 : g + r : r=0
If g < 0 : r - g : g=0
If g < $80 : b + g : EndIf
If g >= $80 : b +$FF-g : EndIf
If r < $80 : b + r : EndIf
If r >= $80
If r < $17F
b =$17F-r : r = $80
Else
r-$FF : b=0
If r>=$1FF:r=$FF:g=$FF:b=$FF :EndIf
If r>$FF :g=r-$FF:b=r-$FF:r=$FF:EndIf
EndIf
EndIf
EndIf
EndIf
ProcedureReturn = RGB(r, g, b)
EndProcedure
Procedure draw_(n.l, *u.FArray, *v.FArray, im.l, jm.l)
Define.l i, j
Define.f um, vm
Define.s info
ClearScreen(0)
ExamineKeyboard()
If KeyboardPushed(#PB_Key_Escape)
End
EndIf
StartDrawing(ScreenOutput())
For i = 1 To n
For j = 1 To n
Box(i * #DSTEP-#DSTEP, j * #DSTEP-#DSTEP, #DSTEP, #DSTEP, color_gr(*u\f[IX(i, j)],*v\f[IX(i, j)]))
Next
Next
;Mouse position
Line(im*#DSTEP-#DSTEP/2,jm*#DSTEP-4.5*#DSTEP,0,8*#DSTEP,RGB(255,0,0))
Line(im*#DSTEP-4.5*#DSTEP,jm*#DSTEP-#DSTEP/2,8*#DSTEP,0,RGB(255,0,0))
Box( n * #DSTEP,0, 200, n * #DSTEP,RGB(255,255,255))
If plot = 0
um =*u\f[IX(im, jm)]
info= "Dichte= "+StrF(um,3)+"kg/m^3"
DrawText (n * #DSTEP+10,40, info)
Else
um=*u\f[IX(im, jm)]
vm=*v\f[IX(im, jm)]
info= "u= "+StrF(um,3)+"m/s"
DrawText (n * #DSTEP+50,40, info)
info= "v= "+StrF(vm,3)+"m/s"
DrawText (n * #DSTEP+50,70, info)
info= "Sqrt(u^2+v^2)= "+StrF(Sqr(um*um+vm*vm),3)+"m/s"
DrawText (n * #DSTEP+10,100, info)
EndIf
info= "RGB= "+HexQ(color_gr(*u\f[IX(im, jm)],*v\f[IX(im, jm)]))
DrawText (n * #DSTEP+10,130, info)
StopDrawing()
WindowEvent()
FlipBuffers()
EndProcedure
Procedure set_dens (n.l, *d.FArray)
Define.l i, Size = (n + 2) * (n + 2)
For i = 0 To Size - 1
*d\f[i]= dense
Next
EndProcedure
Procedure set_datas(n.l, *dens0.FArray, *u0.FArray, *v0.FArray, im.l, jm.l)
Define.l z, Size = (n + 2) * (n + 2)
fillmemory_(*dens0, (n + 2) * (n + 2) * SizeOf(Float), 0)
If MouseButton(1)
*dens0\f[IX(im+z, jm)] = - 0.4/dt ;Dichteaenderung pro Zeitschritt
EndIf
If MouseButton(2)
*u0\f[IX(im+z, jm)] = 10/dt ;Geschwindigkeitsaenderung pro Zeitschritt
*v0\f[IX(im+z, jm)] = 8 /dt ;Geschwindigkeitsaenderung pro Zeitschritt
EndIf
*u0\f[IX(90, 90)] = -2.0 /dt ; hier ist ein Quirl eingebaut, rausnehmen
*v0\f[IX(90, 90)] = -2.0 /dt
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))
Static im.l, jm.l
set_dens(#N, *dens) ; Dichtefeld fuer den Anfang belegen
Repeat
ExamineMouse()
im = im+MouseDeltaX()
jm = jm+MouseDeltaY()
If im > #N : im = #N : EndIf
If im < 1 : im = 1 : EndIf ; 0 ist Randpunkt, darf nicht gesetzt werden
If jm > #N : jm = #N : EndIf
If jm < 1 : jm = 1 : EndIf
set_datas(#N, *dens_prev, *u_prev, *v_prev, im, jm)
vel_step(#N, *u, *v, *u_prev, *v_prev, *dens, visc, dt)
dens_step(#N, *dens, *dens_prev, *u, *v, diff, dt)
If plot = 0
draw_(#N, *dens,*dens, im, jm)
Else
draw_(#N, *u,*v, im, jm)
EndIf
ForEver
EndProcedure
InitSprite()
InitMouse()
InitKeyboard()
OpenWindow(0, 100, 100, #N * #DSTEP+200, #N * #DSTEP, "Fluid Dynamics")
CreateGadgetList(WindowID(0))
OpenWindowedScreen(WindowID(0), 0, 0, #N * #DSTEP+200, #N * #DSTEP, 0, 0, 0)
main()