Simulation de feu, question ?

threedslider · Message par **threedslider** » ven. 24/mars/2023 16:22

Hello à tous,

Est ce que quelqu'un a réussi de faire la simulation de feu en OpenGL avec Purebasic ? J'ai pas réussi à convertir de C++ en Purebasic car à cause de la programmation orienté objet ça complique un peu en procédural de Purebasic ... :/

Vous avez des astuces ou des conseils sur comment le créer la simulation de feu ? Car j'aimerais le faire par expérience

Sinon merci encore de votre aide.

threedslider · Message par **threedslider** » jeu. 04/mai/2023 12:49

Ok bon ben ya personne alors

Je vais lire un bouquin de fluid simulation pour tenter à construire moi même

J'espère vous montrer plus tard quand j'aurais compris, merci quand même si ya personne xD

Josue · Message par **Josue** » mar. 26/sept./2023 12:42

Salut,

Pour créer une simulation de feu en OpenGL avec Purebasic, tu peux suivre ces étapes :

Utilise les shaders : Utilise les shaders GLSL pour manipuler les effets visuels, car ils te permettront de simuler le mouvement du feu de manière plus réaliste.

Texture de flamme : Crée une texture animée qui ressemble à une flamme et applique-la sur un quad (un rectangle) en utilisant OpenGL. Cette texture donnera l'illusion d'une flamme en mouvement.

Variation de l'intensité : Anime la texture de flamme en modifiant son intensité et en jouant avec les couleurs pour simuler le mouvement et les fluctuations du feu.

Déplacement aléatoire : Applique un mouvement aléatoire aux particules de la flamme pour donner un aspect plus naturel au feu.

Intégration du temps : Utilise le temps pour mettre à jour la simulation de manière fluide, en faisant évoluer la texture de flamme au fil du temps.

N'oublie pas de consulter la documentation de Purebasic pour savoir comment utiliser OpenGL avec Purebasic de manière efficace. Bonne chance avec ton projet !

Mindphazer · Message par **Mindphazer** » mer. 27/sept./2023 14:46

Tiens, ChatGPT est parmi nous ?

threedslider · Message par **threedslider** » mer. 27/sept./2023 15:06

Merci beaucoup à tous

Pour ma part je voulais tester celui là : https://github.com/983/Fluid-Simulation mais sans succès au procédural sur Purebasic :/ Et pas besoin de perlin noise !

Et vous vous avez le code sur Purebasic et qui simule bien le feu ?

ChatGPT c'est bof pour moi, je préfères faire moi même ^^

Happy coding !

threedslider · Message par **threedslider** » lun. 16/oct./2023 14:17

Re bonjour à toutes et tous

Voici le code en wip que je voulais vous partager, peut etre pourriez vous m'aider à faire fonctionner ? Merci à l'avance

Code : Tout sélectionner

;;; 
;;; Fluid Simulation v. 0.1 
;;;

w.i = 512
h.i = 512

nx.i = 256
ny.i = 256

dt.f = 0.02
iterations.i = 5
vorticity.f = 10.0

texture = 0

Structure vec2f
  x.f
  y.f
EndStructure

Structure grid
  x.i
  y.i
  velocity.vec2f
  xf.f
  yf.f
  Array Data.i(65536)
EndStructure

MyVar.vec2f

; Dim old_velocity.grid(65536)
; Dim new_velocity.grid(65536)
; 
; Dim old_density.grid(65536)
; Dim new_density.grid(65536)
; 
; Dim pixels.grid(65536)
; 
; 
; Procedure draw(Array var.vec2f(1), n.i, mode.i)
;   glEnableClientState_(#GL_VERTEX_ARRAY)
;   glVertexPointer_(2, #GL_FLOAT, 0, var())
;   glDrawArrays_(mode, 0, n)
; EndProcedure
; 
; 
; 
; Procedure draw_circle(x.f, y.f, r.f, n.i = 100)
;   Dim pos.vec2f(n)
;   
;   Shared MyVar
;   
;   For i = 0 To n
;     angle.f = 2.0*3.14159*i/n
;     MyVar\x = x
;     MyVar\y = y
;     pos(i)\x = MyVar\x + r*Cos(angle)
;     pos(i)\y = MyVar\y + r*Sin(angle)
;   Next
;   
;   draw(pos(), n, #GL_LINE_LOOP)
;   
;   FreeArray(pos())
;   
; EndProcedure
;   
; 
Procedure init()
  
  Shared nx
  Shared ny
  
  Shared texture
  
;   Shared old_density()
;   Shared old_velocity()
;   
;   For y = 0 To ny
;     For n = 0 To nx
;       
;       old_density(n)\xf = 0.0
;       old_density(y)\yf = 0.0
;       old_velocity(n)\velocity\x = 0.0
;       old_velocity(y)\velocity\y = 0.0
;         
;     Next
;   Next
;   

    glEnable_(#GL_TEXTURE_2D)
    glEnable_(#GL_BLEND)
    glBlendFunc_(#GL_SRC_ALPHA, #GL_ONE_MINUS_SRC_ALPHA)

    glGenTextures_(1, @texture)
    glBindTexture_(#GL_TEXTURE_2D, texture)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)
    glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, nx, ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, 0)

EndProcedure
; 
; Procedure.f lerp(a.f, b.f, c.f)
;   ProcedureReturn (1.0 - c)*a + c*b
; EndProcedure
; 
; Procedure.f interpolationX_d(Array g.grid(1), *var.vec2f, ax.i)
;   ix.i = *var\x 
;   iy.i = *var\y
;   
;   ux.f = *var\x - ix
;   uy.f = *var\y - iy
;   
;   ;ProcedureReturn lerp( lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 0, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 0, ux), uy),  lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 1, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 1, ux), uy), uy)
;   
;   ProcedureReturn lerp(lerp(g(ax)\xf + 0, g(ax)\xf + 1, ux), lerp(g(ax)\xf + 1, g(ax)\xf + 0, ux), uy)
;   
;  EndProcedure
;  
;  Procedure.f interpolationY_d(Array g.grid(1), *var.vec2f, ay.i)
;   ix.i = *var\x 
;   iy.i = *var\y
;   
;   ux.f = *var\x - ix
;   uy.f = *var\y - iy
;   
;    ;ProcedureReturn lerp( lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 0, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 0, ux), uy),  lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 1, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 1, ux), uy), uy)
;   
;   ProcedureReturn lerp(lerp(g(ay)\yf + 0, g(ay)\yf + 1, ux), lerp(g(ay)\yf + 1, g(ay)\yf + 1, ux), uy)
;   
;   
; EndProcedure
; 
; 
; Procedure.f interpolationX_v(Array g.grid(1), *var.vec2f, ax.i)
;   ix.i = *var\x 
;   iy.i = *var\y
;   
;   ux.f = *var\x - ix
;   uy.f = *var\y - iy
;   
;   ;ProcedureReturn lerp( lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 0, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 0, ux), uy),  lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 1, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 1, ux), uy), uy)
;   
;   ProcedureReturn lerp(lerp(g(ax)\velocity\x + 0, g(ax)\velocity\x + 1, ux), lerp(g(ax)\velocity\x + 1, g(ax)\velocity\x + 0, ux), uy)
;   
;  EndProcedure
;  
;  Procedure.f interpolationY_v(Array g.grid(1), *var.vec2f, ay.i)
;   ix.i = *var\x 
;   iy.i = *var\y
;   
;   ux.f = *var\x - ix
;   uy.f = *var\y - iy
;   
;    ;ProcedureReturn lerp( lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 0, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 0, ux), uy),  lerp( lerp(g(ax)\xf + 0, g(ay)\yf + 1, ux), lerp(g(ax)\xf + 1, g(ay)\yf + 1, ux), uy), uy)
;   
;   ProcedureReturn lerp(lerp(g(ay)\velocity\y + 0, g(ay)\velocity\y + 1, ux), lerp(g(ay)\velocity\y + 1, g(ay)\velocity\y + 1, ux), uy)
;   
;   
; EndProcedure
; 
; 
; Procedure swap_od_2_nd_x(Array od.grid(1), Array nd.grid(1), ax.i)
;   Swap od(ax)\xf, nd(ax)\xf
; EndProcedure
; 
; Procedure swap_od_2_nd_y(Array od.grid(1), Array nd.grid(1), ay.i)
;   Swap od(ay)\yf, nd(ay)\yf
; EndProcedure
; 
; Procedure swap_ov_2_nv_x(Array ov.grid(1), Array nv.grid(1), ax.i)
;   Swap ov(ax)\velocity\x, nv(ax)\velocity\x
; EndProcedure
; 
; Procedure swap_ov_2_nv_y(Array ov.grid(1), Array nv.grid(1), ay.i)
;   Swap ov(ay)\velocity\y, nv(ay)\velocity\y
; EndProcedure
; 
;  
;  
;  Procedure advect_density()
;    Shared nx
;    Shared ny
;    
;    Shared old_density()
;    Shared new_density()
;    Shared old_velocity()
;   
;    pos.vec2f
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        pos\x = dt* old_velocity(x)\velocity\x
;        pos\y = dt* old_velocity(y)\velocity\y
;        
;        new_density(x)\xf = interpolationX_d(old_density(), @pos, x)
;        new_density(y)\yf = interpolationY_d(old_density(), @pos, y)
;      Next
;    Next
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        swap_od_2_nd_x(old_density(), new_density(), x)
;        swap_od_2_nd_y(old_density(), new_density(), y)
;      Next 
;    Next   
;  EndProcedure
;  
;  
;   Procedure advect_velocity()
;    Shared nx
;    Shared ny
;    
;    Shared old_velocity()
;    Shared old_density()
;    
;    Shared new_density()
;    Shared new_velocity()
;    
;    pos.vec2f
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        pos\x = dt* old_velocity(x)\velocity\x
;        pos\y = dt* old_velocity(y)\velocity\y
;        
;        new_velocity(x)\velocity\x = interpolationX_v(old_velocity(), @pos, x)
;        new_velocity(y)\velocity\y = interpolationY_v(old_velocity(), @pos, y)
;      Next
;    Next
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        swap_ov_2_nv_x(old_velocity(), new_velocity(), x)
;        swap_ov_2_nv_y(old_velocity(), new_velocity(), y)
;      Next 
;    Next   
;  EndProcedure
;  
;  
;  Procedure diffuse_density()
;    Shared dt
;    Shared nx
;    Shared ny
;    
;    Shared old_density()
;    Shared new_density()
;    
;    sum.vec2f
;    
;    viscosity.f = dt * 100.01
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        
;        sum\x = viscosity * ( old_density(x-1)\xf + old_density(x+1)\xf + old_density(x+0)\xf + old_density(x+0)\xf ) + old_density(x+0)\xf
;        sum\y = viscosity * ( old_density(y+0)\yf + old_density(y+0)\yf + old_density(y-1)\yf + old_density(y+1)\yf ) + old_density(y+0)\yf
;        
;        new_density(x)\xf = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
;        new_density(y)\yf = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
;      Next
;    Next
;    
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        swap_od_2_nd_x(old_density(), new_density(), x)
;        swap_od_2_nd_y(old_density(), new_density(), y)
;      Next
;    Next
;    
;  EndProcedure
;  
;  
;   Procedure diffuse_velocity()
;    Shared dt
;    Shared nx
;    Shared ny
;    
;    Shared old_velocity()
;    Shared new_velocity()
;    
;    sum.vec2f
;    
;    viscosity.f = dt * 0.000001
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        
;        sum\x = viscosity * ( old_velocity(x-1)\velocity\x + old_velocity(x+1)\velocity\x + old_velocity(x+0)\velocity\x + old_velocity(x+0)\velocity\x ) + old_velocity(x+0)\velocity\x
;        sum\y = viscosity * ( old_velocity(y+0)\velocity\y + old_velocity(y+0)\velocity\y + old_velocity(y-1)\velocity\y + old_velocity(y+1)\velocity\y ) + old_velocity(y+0)\velocity\y
;        
;        new_velocity(x)\velocity\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
;        new_velocity(y)\velocity\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
;      Next
;    Next
;    
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        swap_ov_2_nv_x(old_velocity(), new_velocity(), x)
;        swap_ov_2_nv_y(old_velocity(), new_velocity(), y)
;        
;      Next
;    Next
;    
;  EndProcedure    
;  
;  
;  Procedure project_velocity()
;    
;    Shared nx
;    Shared ny
;    
;    Shared iterations
;    
;    Shared old_velocity()
;    
;    Dim p.grid(65536)
;    Dim p2.grid(65536)
;    Dim div.grid(65536)
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        dx.f = old_velocity(x+1)\velocity\x - old_velocity(x-1)\velocity\x
;        dy.f = old_velocity(y+1)\velocity\y - old_velocity(y-1)\velocity\y
;        
;        div(x)\velocity\x = dx + dy
;        div(y)\velocity\y = dx + dy
;        
;        p(x)\velocity\x = 0.0
;        p(y)\velocity\y = 0.0
;      Next
;    Next
;    
;    
;    For k = 0 To iterations
;      
;      For y = 0 To ny
;        For x = 0 To nx
;          sumX.f = -div(x)\velocity\x + p(x+1)\velocity\x + p(x-1)\velocity\x + p(x+0)\velocity\x + p(x+0)\velocity\x
;          sumY.f = -div(y)\velocity\y + p(y+0)\velocity\y + p(y+0)\velocity\y + p(y+1)\velocity\y + p(y+1)\velocity\y
;          
;          p2(x)\velocity\x = 0.25 * sumX
;          p2(y)\velocity\y = 0.25 * sumY
;        Next
;      Next
;      
;      For y = 0 To ny
;        For x = 0 To nx
;          swap_ov_2_nv_x(p(), p2(), x)
;          swap_ov_2_nv_y(p(), p2(), y)
;        Next
;      Next
;      
;    Next
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        old_velocity(x)\velocity\x = old_velocity(x)\velocity\x - 0.5 * p(x+1)\velocity\x - p(x-1)\velocity\x
;        old_velocity(y)\velocity\y = old_velocity(y)\velocity\y - 0.5 * p(y+1)\velocity\y - p(y-1)\velocity\y
;      Next
;    Next
;    
;    
;    FreeArray(p())
;    FreeArray(p2())
;    FreeArray(div())
;  EndProcedure
;  
;  
;  Procedure.f curl(ax.i, ay.i)
;    
;    Shared old_velocity()
;    
;    ProcedureReturn old_velocity(ax)\velocity\x + old_velocity(ay+1)\y - old_velocity(ax)\velocity\x - old_velocity(ay-1)\velocity\y + old_velocity(ax-1)\velocity\x + old_velocity(ay)\velocity\y - old_velocity(ax+1)\velocity\x + old_velocity(ay)\velocity\y 
;    
;  EndProcedure
;  
;  Procedure.f dot(*A_v.vec2f, *B_v.vec2f)
;    ProcedureReturn *A_v\x * *B_v\x + *B_v\y * *B_v\y
;  EndProcedure
;   
;  Procedure.f length(*var.vec2f)
;    ProcedureReturn Sqr(dot(var, var))
;  EndProcedure
;  
;  Procedure vorticity_confinement()
;    Dim abs_curl.grid(65536)
;    
;    Shared nx
;    Shared ny
;    
;    Shared new_velocity()
;    Shared old_velocity()
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        abs_curl(x)\velocity\x = Abs(curl(x,y))
;        abs_curl(y)\velocity\y = Abs(curl(x,y))
;      Next
;    Next
;    
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        direction.vec2f
;        
;        direction\x = abs_curl(x+0)\velocity\x + abs_curl(y-1)\velocity\y - abs_curl(x+0)\velocity\x + abs_curl(y+1)\velocity\y
;        direction\y = abs_curl(x+1)\velocity\x + abs_curl(y+0)\velocity\y - abs_curl(x+1)\velocity\x + abs_curl(y+0)\velocity\y
;        
;        direction\x = vorticity / (direction\x +1e-5) * direction\x
;        direction\y = vorticity / (direction\y +1e-5) * direction\y
;  
;        If x < nx/2
;          direction\x = direction\x * 0.0
;          direction\y = direction\y * 0.0
;        EndIf
;        
;        new_velocity(x)\velocity\x = old_velocity(x)\velocity\x + dt * curl(x, y) * direction\x
;        new_velocity(y)\velocity\y = old_velocity(y)\velocity\y + dt * curl(x, y) * direction\x
;      Next
;    Next
;    
;    
;    For y = 0 To ny
;      For x = 0 To nx
;        swap_ov_2_nv_x(old_velocity(), new_velocity(), x)
;        swap_ov_2_nv_y(old_velocity(), old_velocity(), y)
;      Next
;    Next
;    
; FreeArray(abs_curl())   
;    
; EndProcedure
; 
; 
; Procedure.f clamp(x.f, a.f, b.f)
;   If x < a
;     ProcedureReturn a
;   ElseIf x > b
;     ProcedureReturn b
;   Else
;     ProcedureReturn x
;   EndIf
; EndProcedure
;   
;  
;  Procedure.f smoothstep(a.f, b.f, u.f)
;    t.f = clamp((u-a) / (b-a), 0.0, 1.0)
;    
;    ProcedureReturn t*t * (3.0 - 2.0 * t)
;  EndProcedure
;  
;    
;  
;  
;  Procedure add_density(px.i, py.i, r.i = 10, value.f = 0.5)
;    
;    Shared old_density()
;    
;    For y = -r To r
;      For x = -r To r
;        d.f = Sqr(x*x + y*y)
;        u.f = smoothstep((r), 0.0, d)
;        old_density(x)\xf = old_density(x)\xf + u*value
;        old_density(y)\yf = old_density(y)\yf + u*value
;      Next
;    Next
;  EndProcedure
 
 
 Procedure on_frame()
   
   Shared nx
   Shared ny
   
   ;Shared old_density()
   
   ;Shared pixels()
   
   Dim Mydata.f(65536)
   
    glClearColor_(0.1, 0.2, 0.3, 0.0)
    glClear_(#GL_COLOR_BUFFER_BIT)

    ;fluid_simulation_step()

    ;double t = sec()
    
    ; density field To pixels
    For y = 0 To ny
      For x = 0 To nx
        ;fd.f = old_density(x)\xf
        fd = fd*0.25 + 1.0
        f3.f = fd*fd*fd
        r.f = 1.5*fd
        g.f = 1.5*f3
        b.f = f3*f3
        ;pixels(x)\values(x) = RGBA(r, g, b, 1.0)
        ;pixels(y)\values(y) = RGBA(r, g, b, 1.0)
        MyData(x+y) = RGBA(r, g, b, 1.0)
      Next 
    Next
    
        
    
    ;double dt = sec() - t
    

    ; upload pixels To texture
    glBindTexture_(#GL_TEXTURE_2D, texture)
    glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, nx, ny, #GL_RGBA, #GL_UNSIGNED_BYTE, Mydata())

    ; draw texture
    glBegin_(#GL_QUADS)
    glTexCoord2f_(0.0, 0.0) 
    glVertex2f_(-1.0, -1.0)
    glTexCoord2f_(1.0, 0.0) 
    glVertex2f_(1.0, -1.0)
    glTexCoord2f_(1.0, 1.0) 
    glVertex2f_(1.0, 1.0)
    glTexCoord2f_(0.0, 1.0) 
    glVertex2f_(-1.0, 1.0)
    glEnd_()
    
    
    FreeArray(Mydata())

  EndProcedure
  
  
  
InitKeyboard()
InitSprite()

OpenWindow(0, 0, 0, 800, 600, "OpenGL : Cube for test", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, 800, 600)


glMatrixMode_(#GL_PROJECTION)
glLoadIdentity_() 
gluPerspective_(45.0, 800/600, 1.0, 60.0)
glMatrixMode_(#GL_MODELVIEW)
;glTranslatef_(-1/4, 0, -4)
glShadeModel_(#GL_SMOOTH)
glEnable_(#GL_LIGHT0)
glEnable_(#GL_LIGHTING)
glEnable_(#GL_COLOR_MATERIAL)
glEnable_(#GL_DEPTH_TEST)
glEnable_(#GL_CULL_FACE)   
glViewport_(0, 0, 800, 600)


glTranslatef_(0.0, 0.0, -4.0)

Repeat
;glClearColor_(0.0, 0.0, 0.0, 0.0) ; background color
 ; glClear_(#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)
  Repeat
    event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow
        quit = 1
     EndSelect
   Until event = 0
   
   init()
   
  on_frame()
   
 
   
 FlipBuffers()
 
 ExamineKeyboard()
 
Until KeyboardPushed(#PB_Key_Escape) Or quit = 1

  
       
  
   
   
 
  
; FreeArray(old_density())
; FreeArray(old_velocity())
; FreeArray(new_density())
; FreeArray(new_velocity())
; FreeArray(pixels())

Je serais ravi grandement si vous y arriver de faire marcher ça

Voilà happy coding !

PS : Veuillez voir le lien que j'ai indiqué precedement le post. (Pour ce qui est C++ pour convertir au Purebasic)

threedslider · Message par **threedslider** » mar. 17/oct./2023 1:00

Comme vous pouvez le constater mon code de simulation de fluide basé sur C++ au Purebasic est un peu trop brouillon, donc je vous demandes des conseils sur l’implémentation en Purebasic, est ce que les mémoires sont bien codés ? ou bien il y a d'autres façon de coder ? Vous verrez aussi qu'en orienté objet j'ai décomposé en chaque procédural séparé pour faciliter la programmation. Mais encore une fois sont ils tous bien codé ou pas ?

Vous avez des idées sur comment bien implémenter de cette inspiration de Fluid simulation en C++ pour ma version en purebasic ?

Avec votre aide je pourrais mieux comprendre comment le faire

, donc svp aidez moi à fixer cela !

C'est vrai le C++ en Purebasic c'est une vraie casse tête xD

Je suis ouvert à tout de votre commentaire, conseil, amélioration ou proposition ou même critique, car je suis un peu bloqué tout de même

Merci encore et vous souhaite de bonne soirée

Le toujours débutant qui veut progresser ^^

see you soon and happy coding !

threedslider · Message par **threedslider** » ven. 03/nov./2023 13:13

Bon ok je vois quand même mes erreurs au niveau d'algorithme comme mémoire et implémentation

Sinon j'espère je vais corriger cela en attendant

Merci quand même

threedslider · Message par **threedslider** » mer. 08/nov./2023 18:04

Hello à tous,

Bon j’étais tout proche mais toujours pas fonctionnel correctement

Donc l'utilisateur "pjay" dans le forum anglais de Purebasic a réécrit mes codes à sa manière maintenant il marche ! Ouf Nous sommes arrivés

Voici ses codes :

Code : Tout sélectionner

; fluid dynamics (flame) - PB conversion of https://github.com/983/Fluid-Simulation  - Phil James 2023

EnableExplicit

CompilerIf #PB_Compiler_OS=#PB_OS_Linux
  Structure point
    x.i : y.i
  EndStructure
CompilerEndIf

Structure vec2f
  x.f : y.f
EndStructure

;#nx = 256 : #ny = 256 : #Upscale = 4 ; ~7ms on 5800h
;#nx = 160 : #ny = 120 : #Upscale = 8 ; ~2ms on 5800h
;#nx = 320 : #ny = 160 : #Upscale = 6 ; ~5ms on 5800h
;#nx = 512 : #ny = 256 : #Upscale = 4 ; ~12ms on 5800h
#nx = 320 : #ny = 400 : #Upscale = 2 ; ~12ms on 5800h

#nx1 = #nx - 1 : #ny1 = #ny - 1 : #nx2 = #nx - 2 : #ny2 = #ny - 2
Global dt.f = 0.02, iterations = 5, vorticity.f = 8.0, mouse.point, lmb, w.i = #nx * #Upscale, h = #ny * #Upscale, texture
Global v1size = #nx * #ny * SizeOf(float), v2size = v1size * 2, onframe

Global Dim old_velocity.vec2f(#ny1,#nx1), Dim new_velocity.vec2f(#ny1,#nx1)
Global Dim old_density.f(#ny1, #nx1), Dim new_density.f(#ny1, #nx1)
Global Dim pixels.l(#ny1,#nx1)

Macro LengthM(x,y) : Sqr((x*x) + (y*y)) : EndMacro
Macro LerpM(a, b, t) : ( a + t * (b - a) ) : EndMacro
Macro Log2M(v) : (Log(v) / 0.6931471805599453094) : EndMacro ; not sure if this is correct Log2...
Macro RandomFloatM(a, b) : LerpM(a, b, Random(1000000) / 1000000.0) : EndMacro
 
Procedure Init()
  Protected x, y
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) = 0.0 : old_velocity(y,x)\x = 0 : old_velocity(y,x)\y = 0
  Next : Next

  OpenWindow(0,0,0,w,h,"Fluid-Simulation (Flame)",#PB_Window_ScreenCentered|#PB_Window_SystemMenu)
  OpenGLGadget(0,0,0,w,h,#PB_OpenGL_NoFlipSynchronization)
  glOrtho_(0,1,1,0,0,1)
  glEnable_(#GL_TEXTURE_2D);
  
  glGenTextures_(1, @texture);
  glBindTexture_(#GL_TEXTURE_2D, texture);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)    ;
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)    ;
  glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, #nx, #ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #Null);
EndProcedure
Procedure.f interpolate(Array grid.f(2), *p.vec2f, *d.float)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\f = LerpM(LerpM(grid(iy,ix),grid(iy,ix + 1), ux),LerpM(grid(iy+1,ix),grid(iy+1,ix+1), ux), uy)
EndProcedure
Procedure.f interpolateVec2(Array grid.vec2f(2), *p.vec2f, *d.vec2f)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\x = LerpM(LerpM(grid(iy,ix)\x,grid(iy,ix+1)\x, ux), LerpM(grid(iy+1,ix)\x,grid(iy+1,ix+1)\x, ux), uy)
  *d\y = LerpM(LerpM(grid(iy,ix)\y,grid(iy,ix+1)\y, ux), LerpM(grid(iy+1,ix)\y,grid(iy+1,ix+1)\y, ux), uy)
EndProcedure

Procedure advect_density()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x
      pos\y = y - dt * old_velocity(y,x)\y
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolate(old_density(), @pos, @new_density(y,x));
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure advect_velocity()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x;
      pos\y = y - dt * old_velocity(y,x)\y;
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolateVec2(old_velocity(), @pos,@new_velocity(y,x))
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure diffuse_density()
  Protected pos.vec2f, x, y, diffusion.f = dt * 100.01, sum.f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum = diffusion * (old_density(y, x - 1) + old_density(y, x + 1) + old_density(y - 1, x) + old_density(y + 1, x)) + old_density(y,x);
      new_density(y, x) = 1.0 / (1.0 + 4.0 * diffusion) * sum 
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure diffuse_velocity()
  Protected viscosity = dt*0.000001, x, y, sum.vec2f, gotnan.a  = #False
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum\x = viscosity * (old_velocity(y, x - 1)\x + old_velocity(y, x + 1)\x + old_velocity(y - 1, x)\x + old_velocity(y + 1, x)\x) + old_velocity(y, x)\x
      sum\y = viscosity * (old_velocity(y, x - 1)\y + old_velocity(y, x + 1)\y + old_velocity(y - 1, x)\y + old_velocity(y + 1, x)\y) + old_velocity(y, x)\y
      new_velocity(y,x)\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
      new_velocity(y,x)\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure project_velocity()
  Protected x, y, dx.f, dy.f, k, sum.f
  Static Dim p.f(#ny,#nx), Dim p2.f(#ny,#nx), Dim div.f(#ny,#nx)

  For y = 1 To #ny2 : For x = 1 To #nx2
      dx = old_velocity(y    , x + 1)\x - old_velocity(y    ,x - 1)\x
      dy = old_velocity(y + 1, x    )\y - old_velocity(y - 1,x    )\y
      div(y,x) = dx + dy : p(y,x) = 0.0
  Next : Next
  For k = 1 To iterations
    For y = 1 To #ny2 : For x = 1 To #nx2
        sum = -div(y,x) + p(y ,x + 1) + p(y ,x - 1) + p(y + 1,x) + p(y - 1,x)
        p2(y,x) = 0.25 * sum
    Next : Next
    CopyMemory(@p2(),@p(),v1size)
  Next
  For y = 1 To #ny2 : For x = 1 To #nx2
      old_velocity(y,x)\x - (0.5 * (p(y    ,x + 1) - p(y    , x - 1)))
      old_velocity(y,x)\y - (0.5 * (p(y + 1,x  ) - p(y - 1, x )))
  Next : Next
EndProcedure
Procedure.f curl(y,x)
  ProcedureReturn old_velocity(y + 1,x)\x - old_velocity(y - 1,x)\x + old_velocity(y,x - 1)\y - old_velocity(y, x + 1)\y;
EndProcedure
Procedure vorticity_confinement()
  Protected direction.vec2f, length.f, x, y;
  Static Dim abs_curl.f(#ny1,#nx1)
  For y = 1 To #ny2 : For x = 1 To #nx2
      abs_curl(y,x) = Abs(curl(y,x));
  Next : Next
  
  For y = 1 To #ny2 : For x = 1 To #nx2
      direction\x = abs_curl(y - 1, x    ) - abs_curl(y + 1, x    );
      direction\y = abs_curl(y    , x + 1) - abs_curl(y    , x - 1);
      length = lengthM(direction\x,direction\y)
      direction\x = vorticity / (length + 1e-5) * direction\x;
      direction\y = vorticity / (length + 1e-5) * direction\y;
      new_velocity(y,x)\x = old_velocity(y,x)\x + dt * curl(y,x) * direction\x;
      new_velocity(y,x)\y = old_velocity(y,x)\y + dt * curl(y,x) * direction\y;
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure.f Clamp(v.f,Low.f = 0.0, High.f = 1.0)
  If v < Low : ProcedureReturn Low : ElseIf v > High : ProcedureReturn High : EndIf
  ProcedureReturn v
EndProcedure
Procedure.f SmoothStep(Low.f, High.f, v.f)
  V = Clamp((v - low) / (high - low));
  ProcedureReturn v * v * (3.0 - 2.0 * v);
EndProcedure

Procedure add_density(px, py, r = 10, value.f = 0.5)
  Protected x, y, u.f, rf.f = r, x2, y2
  For y = -r To r : For x = -r To  r
      y2 = py + y : x2 = px + x
      If x2 < 0 : x2 = 0 : ElseIf x2 > #nx1 : x2 = #nx1 : EndIf
      If y2 < 0 : y2 = 0 : ElseIf y2 > #ny1 : y2 = #ny1 : EndIf
      u = SmoothStep(rf, 0.0, Sqr(x*x + y*y));
      old_density(y2, x2) + u * value
  Next : Next
  
EndProcedure

Procedure fluid_simulation_step()
  Protected x,y, r.f = 10.0
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x + RandomFloatM(-r, r) : old_velocity(y,x)\y + RandomFloatM(-r, r);
  Next : Next
  
  ;// dense regions rise up
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\y + (old_density(y,x) * 20.0 - 5.0) * dt;
  Next : Next
 
  ;// fast movement is dampened
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x * 0.999 : old_velocity(y,x)\y * 0.999;
  Next : Next
  
  ;// fade away
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) * 0.99;
  Next : Next

  add_density(#nx*0.25, #ny*0.1, 10,0.5);
  add_density(#nx*0.75, #ny*0.1, 10,0.5);
  
  If mouse\x > 0 And mouse\y > 0 And mouse\x < #nx And mouse\y < #ny1 : add_density(Mouse\x, #ny - Mouse\y, #nx * 0.025,0.5) : EndIf
  
  vorticity_confinement()
  advect_velocity()
  project_velocity()
  diffuse_velocity()
  diffuse_density();
  advect_density()
  
  ;/ wipe out borders as prevents crashes
  For y = 0 To #ny1 
      old_density(y,0) = 0.0 : old_velocity(y,0)\x = 0 : old_velocity(y,0)\y = 0      
      old_density(y,1) = 0.0 : old_velocity(y,1)\x = 0 : old_velocity(y,1)\y = 0    
      old_density(y,#nx1) = 0.0 : old_velocity(y,#nx1)\x = 0 : old_velocity(y,#nx1)\y = 0      
      old_density(y,#nx2) = 0.0 : old_velocity(y,#nx2)\x = 0 : old_velocity(y,#nx2)\y = 0      
  Next
  For x = 0 To #nx1 
      old_density(0,x) = 0.0 : old_velocity(0,x)\x = 0 : old_velocity(0,x)\y = 0      
      old_density(#ny1,x) = 0.0 : old_velocity(#ny1,x)\x = 0 : old_velocity(#ny1,x)\y = 0      
      old_density(1,x) = 0.0 : old_velocity(1,x)\x = 0 : old_velocity(1,x)\y = 0      
      old_density(#ny2,x) = 0.0 : old_velocity(#ny2,x)\x = 0 : old_velocity(#ny2,x)\y = 0      
  Next
EndProcedure

Procedure.l RGBAgl(r.f, g.f, b.f, a.f)
  r = clamp(r * 255, 0, 255) : g = clamp(g * 255, 0, 255) : b = clamp(b * 255, 0, 255) : a = clamp(a * 255, 0, 255)
  ProcedureReturn (Int(a) << 24) | (Int(b) << 16) | (Int(g) << 8) | Int(r);
EndProcedure
Procedure Render()
  ;// density field To pixels
  Protected x, y, f.f, f3.f, r.f, g.f, b.f
  For y = 0 To #ny1 : For x = 0 To #nx1
      f = old_density(y,x);
      f = log2m( f * 0.25 + 1.0);
      f3 = f*f*f
      r = 1.5 * f : g = 1.5 * f3 : b = (f3 *f3)
      pixels(y,x) = RGBAgl(r, g, b, 1.0);
  Next : Next
  
  ;// upload pixels() array To texture
  glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, #nx, #ny, #GL_RGBA, #GL_UNSIGNED_BYTE, @pixels());
  
  ;// draw texture
  glBegin_(#GL_QUADS);
  glTexCoord2f_(0.0, 1.0): glVertex2f_(0  , 0)
  glTexCoord2f_(1.0, 1.0): glVertex2f_(1.0, 0)
  glTexCoord2f_(1.0, 0.0): glVertex2f_(1.0, 1.0)
  glTexCoord2f_(0.0, 0.0): glVertex2f_(0  , 1.0)
  glEnd_()
  
  SetGadgetAttribute(0,#PB_OpenGL_FlipBuffers,#True)
EndProcedure
Init()

Define Event, Exit, time.i

Repeat
  Repeat
    Event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow : Exit = #True
      Case #PB_Event_Gadget
        Select EventType()
          Case #PB_EventType_MouseMove
            Mouse\x = (GetGadgetAttribute(0,#PB_OpenGL_MouseX) / (w / #nx)) / DesktopResolutionX()
            Mouse\y = (GetGadgetAttribute(0,#PB_OpenGL_MouseY) / (h / #ny)) / DesktopResolutionY()
          Case #PB_EventType_LeftButtonDown
            LMB = #True
            add_density(Mouse\x, #ny - Mouse\y, #nx * 0.2,8.0); burst
          Case #PB_EventType_LeftButtonUp
            LMB = #False
        EndSelect
    EndSelect
  Until Event = 0

  time = ElapsedMilliseconds()
  fluid_simulation_step();
  Render() : onframe + 1
  time = ElapsedMilliseconds() - time
  If onframe = 60 : onframe = 1
    SetWindowTitle(0,"Fluid dynamics w/vorticity - Phil James 2023 - Process time: "+ Str(time)+"ms.")
  EndIf
  
Until Exit = #True
glDeleteTextures_(1,@texture)

Happy coding !

threedslider · Message par **threedslider** » mer. 08/nov./2023 21:17

Et encore une amélioration par pjay

Code : Tout sélectionner

; fluid dynamics (flame) - PB conversion of https://github.com/983/Fluid-Simulation  - Phil James 2023

EnableExplicit

CompilerIf #PB_Compiler_OS=#PB_OS_Linux
  Structure point
    x.i : y.i
  EndStructure
CompilerEndIf

Structure vec2f
  x.f : y.f
EndStructure

;#nx = 256 : #ny = 256 : #Upscale = 4 ; ~7ms on 5800h
;#nx = 160 : #ny = 120 : #Upscale = 8 ; ~2ms on 5800h
;#nx = 320 : #ny = 160 : #Upscale = 6 ; ~5ms on 5800h
;#nx = 512 : #ny = 256 : #Upscale = 4 ; ~12ms on 5800h
#nx = 320 : #ny = 400 : #Upscale = 2 ; ~12ms on 5800h

#nx1 = #nx - 1 : #ny1 = #ny - 1 : #nx2 = #nx - 2 : #ny2 = #ny - 2
Global dt.f = 0.02, iterations = 5, vorticity.f = 10.0, mouse.point, lmb, w.i = #nx * #Upscale, h = #ny * #Upscale, texture
Global v1size = #nx * #ny * SizeOf(float), v2size = v1size * 2, onframe

Global Dim old_velocity.vec2f(#ny1,#nx1), Dim new_velocity.vec2f(#ny1,#nx1)
Global Dim old_density.f(#ny1, #nx1), Dim new_density.f(#ny1, #nx1)
Global Dim pixels.l(#ny1,#nx1)

Macro LengthM(x,y) : Sqr((x*x) + (y*y)) : EndMacro
Macro LerpM(a, b, t) : ( a + t * (b - a) ) : EndMacro
Macro Log2M(v) : (Log(v) / 0.6931471805599453094) : EndMacro ; not sure if this is correct Log2...
Macro RandomFloatM(a, b) : LerpM(a, b, Random(1000000) / 1000000.0) : EndMacro
 
Procedure Init()
  Protected x, y
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) = 0.0 : old_velocity(y,x)\x = 0 : old_velocity(y,x)\y = 0
  Next : Next

  OpenWindow(0,0,0,w,h,"Fluid-Simulation (Flame)",#PB_Window_ScreenCentered|#PB_Window_SystemMenu)
  OpenGLGadget(0,0,0,w,h);,#PB_OpenGL_NoFlipSynchronization)
  glOrtho_(0,1,1,0,0,1)
  glEnable_(#GL_TEXTURE_2D);
  
  glGenTextures_(1, @texture);
  glBindTexture_(#GL_TEXTURE_2D, texture);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)    ;
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)    ;
  glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, #nx, #ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #Null);
EndProcedure
Procedure.f interpolate(Array grid.f(2), *p.vec2f, *d.float)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\f = LerpM(LerpM(grid(iy,ix),grid(iy,ix + 1), ux),LerpM(grid(iy+1,ix),grid(iy+1,ix+1), ux), uy)
EndProcedure
Procedure.f interpolateVec2(Array grid.vec2f(2), *p.vec2f, *d.vec2f)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\x = LerpM(LerpM(grid(iy,ix)\x,grid(iy,ix+1)\x, ux), LerpM(grid(iy+1,ix)\x,grid(iy+1,ix+1)\x, ux), uy)
  *d\y = LerpM(LerpM(grid(iy,ix)\y,grid(iy,ix+1)\y, ux), LerpM(grid(iy+1,ix)\y,grid(iy+1,ix+1)\y, ux), uy)
EndProcedure

Procedure advect_density()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x
      pos\y = y - dt * old_velocity(y,x)\y
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolate(old_density(), @pos, @new_density(y,x));
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure advect_velocity()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x;
      pos\y = y - dt * old_velocity(y,x)\y;
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolateVec2(old_velocity(), @pos,@new_velocity(y,x))
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure diffuse_density()
  Protected pos.vec2f, x, y, diffusion.f = dt * 100.01, sum.f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum = diffusion * (old_density(y, x - 1) + old_density(y, x + 1) + old_density(y - 1, x) + old_density(y + 1, x)) + old_density(y,x);
      new_density(y, x) = 1.0 / (1.0 + 4.0 * diffusion) * sum 
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure diffuse_velocity()
  Protected viscosity = dt*0.000001, x, y, sum.vec2f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum\x = viscosity * (old_velocity(y, x - 1)\x + old_velocity(y, x + 1)\x + old_velocity(y - 1, x)\x + old_velocity(y + 1, x)\x) + old_velocity(y, x)\x
      sum\y = viscosity * (old_velocity(y, x - 1)\y + old_velocity(y, x + 1)\y + old_velocity(y - 1, x)\y + old_velocity(y + 1, x)\y) + old_velocity(y, x)\y
      new_velocity(y,x)\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
      new_velocity(y,x)\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure project_velocity()
  Protected x, y, dx.f, dy.f, k, sum.f
  Static Dim p.f(#ny,#nx), Dim p2.f(#ny,#nx), Dim div.f(#ny,#nx)

  For y = 1 To #ny2 : For x = 1 To #nx2
      dx = old_velocity(y    , x + 1)\x - old_velocity(y    ,x - 1)\x
      dy = old_velocity(y + 1, x    )\y - old_velocity(y - 1,x    )\y
      div(y,x) = dx + dy : p(y,x) = 0.0
  Next : Next
  For k = 1 To iterations
    For y = 1 To #ny2 : For x = 1 To #nx2
        sum = -div(y,x) + p(y ,x + 1) + p(y ,x - 1) + p(y + 1,x) + p(y - 1,x)
        p2(y,x) = 0.25 * sum
    Next : Next
    CopyMemory(@p2(),@p(),v1size)
  Next
  For y = 1 To #ny2 : For x = 1 To #nx2
      old_velocity(y,x)\x - (0.5 * (p(y    ,x + 1) - p(y    , x - 1)))
      old_velocity(y,x)\y - (0.5 * (p(y + 1,x  ) - p(y - 1, x )))
  Next : Next
EndProcedure
Procedure.f curl(y,x)
  ProcedureReturn old_velocity(y + 1,x)\x - old_velocity(y - 1,x)\x + old_velocity(y,x - 1)\y - old_velocity(y, x + 1)\y;
EndProcedure
Procedure vorticity_confinement()
  Protected direction.vec2f, length.f, x, y;
  Static Dim abs_curl.f(#ny1,#nx1)
  For y = 1 To #ny2 : For x = 1 To #nx2
      abs_curl(y,x) = Abs(curl(y,x));
  Next : Next
  
  For y = 1 To #ny2 : For x = 1 To #nx2
      direction\x = abs_curl(y - 1, x    ) - abs_curl(y + 1, x    );
      direction\y = abs_curl(y    , x + 1) - abs_curl(y    , x - 1);
      length = lengthM(direction\x,direction\y)
      direction\x = vorticity / (length + 1e-5) * direction\x;
      direction\y = vorticity / (length + 1e-5) * direction\y;
      
      If x < #nx * 0.5 : direction\x = 0.0 : direction\y = 0.0 : EndIf
      
      new_velocity(y,x)\x = old_velocity(y,x)\x + dt * curl(y,x) * direction\x;
      new_velocity(y,x)\y = old_velocity(y,x)\y + dt * curl(y,x) * direction\y;
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure.f Clamp(v.f,Low.f = 0.0, High.f = 1.0)
  If v < Low : ProcedureReturn Low : ElseIf v > High : ProcedureReturn High : EndIf
  ProcedureReturn v
EndProcedure
Procedure.f SmoothStep(Low.f, High.f, v.f)
  V = Clamp((v - low) / (high - low));
  ProcedureReturn v * v * (3.0 - 2.0 * v);
EndProcedure

Procedure add_density(px, py, r = 10, value.f = 0.5)
  Protected x, y, u.f, rf.f = r, x2, y2
  For y = -r To r : For x = -r To  r
      y2 = py + y : x2 = px + x
      If x2 < 0 : x2 = 0 : ElseIf x2 > #nx1 : x2 = #nx1 : EndIf
      If y2 < 0 : y2 = 0 : ElseIf y2 > #ny1 : y2 = #ny1 : EndIf
      u = SmoothStep(rf, 0.0, Sqr(x*x + y*y));
      old_density(y2, x2) + u * value
  Next : Next
  
EndProcedure

Procedure fluid_simulation_step()
  Protected x,y, r.f = 10.0 ; r is resolution dependant
  
  For y = 0 To #ny1 : For x = 0 To #nx1
      If x > #nx * 0.5
        old_velocity(y,x)\x + RandomFloatM(-r, r)
        old_velocity(y,x)\y + RandomFloatM(-r, r)
      EndIf
  Next : Next
  
  ;// dense regions rise up
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\y + (old_density(y,x) * 20.0 - 5.0) * dt;
  Next : Next
 
  ;// fast movement is dampened
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x * 0.999 : old_velocity(y,x)\y * 0.999;
  Next : Next
  
  ;// fade away
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) * 0.99;
  Next : Next

  add_density(#nx*0.25, #ny*0.1, 10.0, 0.5);
  add_density(#nx*0.75, #ny*0.1, 10.0, 0.5);
  
  If mouse\x > 0 And mouse\y > 0 And mouse\x < #nx And mouse\y < #ny1 : add_density(Mouse\x, #ny - Mouse\y, #nx * 0.025,0.5) : EndIf
  
  vorticity_confinement()
  ;diffuse_velocity()
  ;project_velocity()
  advect_velocity()
  project_velocity()
  ;diffuse_density();
  advect_density()

  ;/ wipe out borders as prevents crashes
  For y = 0 To #ny1 
      old_density(y,0) = 0.0 : old_velocity(y,0)\x = 0 : old_velocity(y,0)\y = 0      
      old_density(y,1) = 0.0 : old_velocity(y,1)\x = 0 : old_velocity(y,1)\y = 0    
      old_density(y,#nx1) = 0.0 : old_velocity(y,#nx1)\x = 0 : old_velocity(y,#nx1)\y = 0      
      old_density(y,#nx2) = 0.0 : old_velocity(y,#nx2)\x = 0 : old_velocity(y,#nx2)\y = 0      
  Next
  For x = 0 To #nx1 
      old_density(0,x) = 0.0 : old_velocity(0,x)\x = 0 : old_velocity(0,x)\y = 0      
      old_density(#ny1,x) = 0.0 : old_velocity(#ny1,x)\x = 0 : old_velocity(#ny1,x)\y = 0      
      old_density(1,x) = 0.0 : old_velocity(1,x)\x = 0 : old_velocity(1,x)\y = 0      
      old_density(#ny2,x) = 0.0 : old_velocity(#ny2,x)\x = 0 : old_velocity(#ny2,x)\y = 0      
  Next
EndProcedure

Procedure.l RGBAgl(r.f, g.f, b.f, a.f)
  r = clamp(r * 255, 0, 255) : g = clamp(g * 255, 0, 255) : b = clamp(b * 255, 0, 255) : a = clamp(a * 255, 0, 255)
  ProcedureReturn (Int(a) << 24) | (Int(b) << 16) | (Int(g) << 8) | Int(r);
EndProcedure
Procedure Render()
  ;// density field To pixels
  Protected x, y, f.f, f3.f, r.f, g.f, b.f
  For y = 0 To #ny1 : For x = 0 To #nx1
      f = old_density(y,x);
      f = log2m( f * 0.25 + 1.0);
      f3 = f*f*f
      r = 1.5 * f : g = 1.5 * f3 : b = (f3 *f3)
      pixels(y,x) = RGBAgl(r, g, b, 1.0);
  Next : Next
  
  ;// upload pixels() array To texture
  glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, #nx, #ny, #GL_RGBA, #GL_UNSIGNED_BYTE, @pixels());
  
  ;// draw texture
  glBegin_(#GL_QUADS);
  glTexCoord2f_(0.0, 1.0): glVertex2f_(0  , 0)
  glTexCoord2f_(1.0, 1.0): glVertex2f_(1.0, 0)
  glTexCoord2f_(1.0, 0.0): glVertex2f_(1.0, 1.0)
  glTexCoord2f_(0.0, 0.0): glVertex2f_(0  , 1.0)
  glEnd_()
  
  SetGadgetAttribute(0,#PB_OpenGL_FlipBuffers,#True)
EndProcedure
Init()

Define Event, Exit, time.i

Repeat
  Repeat
    Event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow : Exit = #True
      Case #PB_Event_Gadget
        Select EventType()
          Case #PB_EventType_MouseMove
            Mouse\x = (GetGadgetAttribute(0,#PB_OpenGL_MouseX) / (w / #nx)) / DesktopResolutionX()
            Mouse\y = (GetGadgetAttribute(0,#PB_OpenGL_MouseY) / (h / #ny)) / DesktopResolutionY()
          Case #PB_EventType_LeftButtonDown
            LMB = #True
            add_density(Mouse\x, #ny - Mouse\y, #nx * 0.2,8.0); burst
          Case #PB_EventType_LeftButtonUp
            LMB = #False
        EndSelect
    EndSelect
  Until Event = 0

  time = ElapsedMilliseconds()
  fluid_simulation_step();
  Render() : onframe + 1
  time = ElapsedMilliseconds() - time
  If onframe = 60 : onframe = 1
    SetWindowTitle(0,"Fluid dynamics w/vorticity - Phil James 2023 - Process time: "+ Str(time)+"ms.")
  EndIf
  
Until Exit = #True
glDeleteTextures_(1,@texture)

PureBasic

Simulation de feu, question ?

Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?

Re: Simulation de feu, question ?