API tools faq
paste
Advertisement
backstreetimrul

fluid.c 4

backstreetimrul
Jul 29th, 2019
117
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 17.64 KB | None | 0 0
raw download clone embed print report
  1. // Usage: Drag with the mouse to add smoke to the fluid. This will also move a "rotor" that disturbs
  2. //        the velocity field at the mouse location. Press the indicated keys to change options
  3. //--------------------------------------------------------------------------------------------------
  4. #include <rfftw.h>              //the numerical simulation FFTW library
  5. #include <GL/glut.h>            //the GLUT graphics library
  6. #include <stdio.h>              //for printing the help text
  7. #include<math.h>
  8. int idx0, idx1, idx2, idx3;
  9. double px0, py0, px1, py1, px2, py2, px3, py3;
  10. float r, g, b, f;
  11.  
  12. //--- SIMULATION PARAMETERS ------------------------------------------------------------------------
  13. const int DIM = 50;             //size of simulation grid
  14. double dt = 0.4;                //simulation time step
  15. float visc = 0.001;             //fluid viscosity
  16. fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current moment
  17. fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment
  18. fftw_real *fx, *fy;             //(fx,fy)   = user-controlled simulation forces, steered with the mouse
  19. fftw_real *rho, *rho0;          //smoke density at the current (rho) and previous (rho0) moment
  20. rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization
  21.  
  22.  
  23. //--- VISUALIZATION PARAMETERS ---------------------------------------------------------------------
  24. int   winWidth, winHeight;      //size of the graphics window, in pixels
  25. int   color_dir = 0;            //use direction color-coding or not
  26. float vec_scale = 1000;         //scaling of hedgehogs
  27. int   draw_smoke = 0;           //draw the smoke or not
  28. int   draw_vecs = 1;            //draw the vector field or not
  29. const int COLOR_BLACKWHITE=0;   //different types of color mapping: black-and-white, rainbow, banded
  30. const int COLOR_RAINBOW=1;
  31. const int COLOR_BANDS=2;
  32. int   scalar_col = 0;           //method for scalar coloring
  33. int   frozen = 0;               //toggles on/off the animation
  34. int cone = 0; //toogle cone on/off
  35.  
  36.  
  37. //------ SIMULATION CODE STARTS HERE -----------------------------------------------------------------
  38.  
  39. //init_simulation: Initialize simulation data structures as a function of the grid size 'n'.
  40. //                 Although the simulation takes place on a 2D grid, we allocate all data structures as 1D arrays,
  41. //                 for compatibility with the FFTW numerical library.
  42. void init_simulation(int n)
  43. {
  44.     int i; size_t dim;
  45.  
  46.     dim     = n*2*(n/2+1)*sizeof(fftw_real);        //Allocate data structures
  47.     vx       = (fftw_real*) malloc(dim);
  48.     vy       = (fftw_real*) malloc(dim);
  49.     vx0      = (fftw_real*) malloc(dim);
  50.     vy0      = (fftw_real*) malloc(dim);
  51.     dim     = n * n * sizeof(fftw_real);
  52.     fx      = (fftw_real*) malloc(dim);
  53.     fy      = (fftw_real*) malloc(dim);
  54.     rho     = (fftw_real*) malloc(dim);
  55.     rho0    = (fftw_real*) malloc(dim);
  56.     plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE);
  57.     plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE);
  58.  
  59.     for (i = 0; i < n * n; i++)                      //Initialize data structures to 0
  60.     { vx[i] = vy[i] = vx0[i] = vy0[i] = fx[i] = fy[i] = rho[i] = rho0[i] = 0.0f; }
  61. }
  62.  
  63.  
  64. //FFT: Execute the Fast Fourier Transform on the dataset 'vx'.
  65. //     'dirfection' indicates if we do the direct (1) or inverse (-1) Fourier Transform
  66. void FFT(int direction,void* vx)
  67. {
  68.     if(direction==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real*)vx,(fftw_complex*)vx);
  69.     else             rfftwnd_one_complex_to_real(plan_cr,(fftw_complex*)vx,(fftw_real*)vx);
  70. }
  71.  
  72. int clamp(float x)
  73. { return ((x)>=0.0?((int)(x)):(-((int)(1-(x))))); }
  74.  
  75. //solve: Solve (compute) one step of the fluid flow simulation
  76. void solve(int n, fftw_real* vx, fftw_real* vy, fftw_real* vx0, fftw_real* vy0, fftw_real visc, fftw_real dt)
  77. {
  78.     fftw_real x, y, x0, y0, f, r, U[2], V[2], s, t;
  79.     int i, j, i0, j0, i1, j1;
  80.  
  81.     for (i=0;i<n*n;i++)
  82.     { vx[i] += dt*vx0[i]; vx0[i] = vx[i]; vy[i] += dt*vy0[i]; vy0[i] = vy[i]; }
  83.  
  84.     for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n )
  85.        for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n )
  86.        {
  87.           x0 = n*(x-dt*vx0[i+n*j])-0.5f;
  88.           y0 = n*(y-dt*vy0[i+n*j])-0.5f;
  89.           i0 = clamp(x0); s = x0-i0;
  90.           i0 = (n+(i0%n))%n;
  91.           i1 = (i0+1)%n;
  92.           j0 = clamp(y0); t = y0-j0;
  93.           j0 = (n+(j0%n))%n;
  94.           j1 = (j0+1)%n;
  95.           vx[i+n*j] = (1-s)*((1-t)*vx0[i0+n*j0]+t*vx0[i0+n*j1])+s*((1-t)*vx0[i1+n*j0]+t*vx0[i1+n*j1]);
  96.           vy[i+n*j] = (1-s)*((1-t)*vy0[i0+n*j0]+t*vy0[i0+n*j1])+s*((1-t)*vy0[i1+n*j0]+t*vy0[i1+n*j1]);
  97.        }
  98.  
  99.     for(i=0; i<n; i++)
  100.       for(j=0; j<n; j++)
  101.       {  vx0[i+(n+2)*j] = vx[i+n*j]; vy0[i+(n+2)*j] = vy[i+n*j]; }
  102.  
  103.     FFT(1,vx0);
  104.     FFT(1,vy0);
  105.  
  106.     for (i=0;i<=n;i+=2)
  107.     {
  108.        x = 0.5f*i;
  109.        for (j=0;j<n;j++)
  110.        {
  111.           y = j<=n/2 ? (fftw_real)j : (fftw_real)j-n;
  112.           r = x*x+y*y;
  113.           if ( r==0.0f ) continue;
  114.           f = (fftw_real)exp(-r*dt*visc);
  115.           U[0] = vx0[i  +(n+2)*j]; V[0] = vy0[i  +(n+2)*j];
  116.           U[1] = vx0[i+1+(n+2)*j]; V[1] = vy0[i+1+(n+2)*j];
  117.  
  118.           vx0[i  +(n+2)*j] = f*((1-x*x/r)*U[0]     -x*y/r *V[0]);
  119.           vx0[i+1+(n+2)*j] = f*((1-x*x/r)*U[1]     -x*y/r *V[1]);
  120.           vy0[i+  (n+2)*j] = f*(  -y*x/r *U[0] + (1-y*y/r)*V[0]);
  121.           vy0[i+1+(n+2)*j] = f*(  -y*x/r *U[1] + (1-y*y/r)*V[1]);
  122.        }
  123.     }
  124.  
  125.     FFT(-1,vx0);
  126.     FFT(-1,vy0);
  127.  
  128.     f = 1.0/(n*n);
  129.     for (i=0;i<n;i++)
  130.        for (j=0;j<n;j++)
  131.        { vx[i+n*j] = f*vx0[i+(n+2)*j]; vy[i+n*j] = f*vy0[i+(n+2)*j]; }
  132. }
  133.  
  134.  
  135. // diffuse_matter: This function diffuses matter that has been placed in the velocity field. It's almost identical to the
  136. // velocity diffusion step in the function above. The input matter densities are in rho0 and the result is written into rho.
  137. void diffuse_matter(int n, fftw_real *vx, fftw_real *vy, fftw_real *rho, fftw_real *rho0, fftw_real dt)
  138. {
  139.     fftw_real x, y, x0, y0, s, t;
  140.     int i, j, i0, j0, i1, j1;
  141.  
  142.     for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n )
  143.         for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n )
  144.         {
  145.             x0 = n*(x-dt*vx[i+n*j])-0.5f;
  146.             y0 = n*(y-dt*vy[i+n*j])-0.5f;
  147.             i0 = clamp(x0);
  148.             s = x0-i0;
  149.             i0 = (n+(i0%n))%n;
  150.             i1 = (i0+1)%n;
  151.             j0 = clamp(y0);
  152.             t = y0-j0;
  153.             j0 = (n+(j0%n))%n;
  154.             j1 = (j0+1)%n;
  155.             rho[i+n*j] = (1-s)*((1-t)*rho0[i0+n*j0]+t*rho0[i0+n*j1])+s*((1-t)*rho0[i1+n*j0]+t*rho0[i1+n*j1]);
  156.         }
  157. }
  158.  
  159. //set_forces: copy user-controlled forces to the force vectors that are sent to the solver.
  160. //            Also dampen forces and matter density to get a stable simulation.
  161. void set_forces(void)
  162. {
  163.     int i;
  164.     for (i = 0; i < DIM * DIM; i++)
  165.     {
  166.         rho0[i]  = 0.995 * rho[i];
  167.         fx[i] *= 0.85;
  168.         fy[i] *= 0.85;
  169.         vx0[i]    = fx[i];
  170.         vy0[i]    = fy[i];
  171.     }
  172. }
  173.  
  174.  
  175. //do_one_simulation_step: Do one complete cycle of the simulation:
  176. //      - set_forces:
  177. //      - solve:            read forces from the user
  178. //      - diffuse_matter:   compute a new set of velocities
  179. //      - gluPostRedisplay: draw a new visualization frame
  180. void do_one_simulation_step(void)
  181. {
  182.     if (!frozen)
  183.     {
  184.       set_forces();
  185.       solve(DIM, vx, vy, vx0, vy0, visc, dt);
  186.       diffuse_matter(DIM, vx, vy, rho, rho0, dt);
  187.       glutPostRedisplay();
  188.     }
  189. }
  190.  
  191.  
  192. //------ VISUALIZATION CODE STARTS HERE -----------------------------------------------------------------
  193.  
  194.  
  195. //rainbow: Implements a color palette, mapping the scalar 'value' to a rainbow color RGB
  196. void rainbow(float value,float* R,float* G,float* B)
  197. {
  198.    const float dx=0.8;
  199.    if (value<0) value=0; if (value>1) value=1;
  200.    value = (6-2*dx)*value+dx;
  201.    *R = max(0.0,(3-fabs(value-4)-fabs(value-5))/2);
  202.    *G = max(0.0,(4-fabs(value-2)-fabs(value-4))/2);
  203.    *B = max(0.0,(3-fabs(value-1)-fabs(value-2))/2);
  204. }
  205.  
  206. //set_colormap: Sets three different types of colormaps
  207. void set_colormap(float vy)
  208. {
  209.    float R,G,B;
  210.  
  211.    if (scalar_col==COLOR_BLACKWHITE)
  212.        R = G = B = vy;
  213.    else if (scalar_col==COLOR_RAINBOW)
  214.        rainbow(vy,&R,&G,&B);
  215.    else if (scalar_col==COLOR_BANDS)
  216.        {
  217.           const int NLEVELS = 7;
  218.           vy *= NLEVELS; vy = (int)(vy); vy/= NLEVELS;
  219.           rainbow(vy,&R,&G,&B);
  220.        }
  221.    glColor3f(R,G,B);
  222. }
  223.  
  224.  
  225. //direction_to_color: Set the current color by mapping a direction vector (x,y), using
  226. //                    the color mapping method 'method'. If method==1, map the vector direction
  227. //                    using a rainbow colormap. If method==0, simply use the white color
  228. void direction_to_color(float x, float y, int method)
  229. {
  230.     //float r,g,b,f;
  231.     if (method)
  232.     {
  233.       f = atan2(y,x) / 3.1415927 + 1;
  234.       r = f;
  235.       if(r > 1) r = 2 - r;
  236.       g = f + .66667;
  237.       if(g > 2) g -= 2;
  238.       if(g > 1) g = 2 - g;
  239.       b = f + 2 * .66667;
  240.       if(b > 2) b -= 2;
  241.       if(b > 1) b = 2 - b;
  242.     }
  243.     else
  244.     { r = g = b = 1; }
  245.     glColor3f(r,g,b);
  246. }
  247.  
  248. //Function to write in graphics window
  249. void drawString(void* font, float x, float y, char* str)
  250. {
  251.     /* Draws string ’str’ in font ’font’, at world (x,y,0) */
  252.     char* ch;
  253.     glRasterPos3f(x, y, 0.0);
  254.     for (ch = str; *ch; ch++)
  255.         glutBitmapCharacter(font, (int)* ch);
  256. }
  257.  
  258. void colorbar()
  259. {
  260.     //drawString(GLUT_BITMAP_HELVETICA_12, 600, 420, fx);
  261.     if (color_dir==1)
  262.     {
  263.         glBegin(GL_QUADS);
  264.         glColor3f(r, g, b);
  265.         glColor3f(0.58, 0.0, 0.82);
  266.         //glColor3f(1.0, 0.0, 0.0);
  267.         glVertex2f(600, 200);
  268.         glVertex2f(650, 200);
  269.         glColor3f(r, g, b);
  270.         glColor3f(0.3, 0.0, 0.51);
  271.         //glColor3f(0.0, 0.0, 1.0);
  272.         glVertex2f(650, 233.33);
  273.         glVertex2f(600, 233.33);
  274.         glColor3f(r, g, b);
  275.         glColor3f(0.3, 0.0, 0.51);
  276.         //glColor3f(0.0, 0.0, 1.0);
  277.         glVertex2f(650, 233.33);
  278.         glVertex2f(600, 233.33);
  279.         glColor3f(r, g, b);
  280.         glColor3f(0.0, 0.0, 1.0);
  281.         //glColor3f(0.0, 0.0, 1.0);
  282.         glVertex2f(600, 266.66);
  283.         glVertex2f(650, 266.66);
  284.         glColor3f(r, g, b);
  285.         glColor3f(0.0, 0.0, 1.0);
  286.         //glColor3f(0.0, 0.0, 1.0);
  287.         glVertex2f(600, 266.66);
  288.         glVertex2f(650, 266.66);
  289.         glColor3f(r, g, b);
  290.         glColor3f(0.0, 1.0, 0.0);
  291.         //glColor3f(0.0, 0.0, 1.0);
  292.         glVertex2f(650, 299.99);
  293.         glVertex2f(600, 299.99);
  294.         glColor3f(r, g, b);
  295.         glColor3f(0.0, 1.0, 0.0);
  296.         //glColor3f(0.0, 0.0, 1.0);
  297.         glVertex2f(650, 299.99);
  298.         glVertex2f(600, 299.99);
  299.         glColor3f(r, g, b);
  300.         glColor3f(1.0, 1.0, 0.0);
  301.         //glColor3f(0.0, 0.0, 1.0);
  302.         glVertex2f(600, 333.33);
  303.         glVertex2f(650, 333.33);
  304.         glColor3f(r, g, b);
  305.         glColor3f(1.0, 1.0, 0.0);
  306.         //glColor3f(0.0, 0.0, 1.0);
  307.         glVertex2f(600, 333.33);
  308.         glVertex2f(650, 333.33);
  309.         glColor3f(r, g, b);
  310.         glColor3f(1.0, 0.5, 0.0);
  311.         //glColor3f(0.0, 0.0, 1.0);
  312.         glVertex2f(650, 366.66);
  313.         glVertex2f(600, 366.66);
  314.         glColor3f(r, g, b);
  315.         glColor3f(1.0, 0.5, 0.0);
  316.         //glColor3f(0.0, 0.0, 1.0);
  317.         glVertex2f(650, 366.66);
  318.         glVertex2f(600, 366.66);
  319.         glColor3f(r, g, b);
  320.         glColor3f(1.0, 0.0, 0.0);
  321.         //glColor3f(0.0, 0.0, 1.0);
  322.         glVertex2f(600, 400);
  323.         glVertex2f(650, 400);
  324.  
  325.         glEnd();
  326.         glFlush();
  327.     }
  328.     else if(color_dir==0)
  329.     {
  330.         glBegin(GL_QUADS);
  331.         glColor3f(1.0, 1.0, 1.0);
  332.         glVertex2f(600, 200);
  333.         glVertex2f(650, 200);
  334.         //glColor3f(0.0, 0.0, 1.0);
  335.         glVertex2f(650, 233.33);
  336.         glVertex2f(600, 233.33);
  337.         //glColor3f(0.0, 0.0, 1.0);
  338.         glVertex2f(650, 233.33);
  339.         glVertex2f(600, 233.33);
  340.         //glColor3f(0.0, 0.0, 1.0);
  341.         glVertex2f(600, 266.66);
  342.         glVertex2f(650, 266.66);
  343.         //glColor3f(0.0, 0.0, 1.0);
  344.         glVertex2f(600, 266.66);
  345.         glVertex2f(650, 266.66);
  346.         //glColor3f(0.0, 0.0, 1.0);
  347.         glVertex2f(650, 299.99);
  348.         glVertex2f(600, 299.99);
  349.         //glColor3f(0.0, 0.0, 1.0);
  350.         glVertex2f(650, 299.99);
  351.         glVertex2f(600, 299.99);
  352.         //glColor3f(0.0, 0.0, 1.0);
  353.         glVertex2f(600, 333.33);
  354.         glVertex2f(650, 333.33);
  355.         //glColor3f(0.0, 0.0, 1.0);
  356.         glVertex2f(600, 333.33);
  357.         glVertex2f(650, 333.33);
  358.         //glColor3f(0.0, 0.0, 1.0);
  359.         glVertex2f(650, 366.66);
  360.         glVertex2f(600, 366.66);
  361.         //glColor3f(0.0, 0.0, 1.0);
  362.         glVertex2f(650, 366.66);
  363.         glVertex2f(600, 366.66);
  364.         //glColor3f(0.0, 0.0, 1.0);
  365.         glVertex2f(600, 400);
  366.         glVertex2f(650, 400);
  367.  
  368.         glEnd();
  369.         glFlush();
  370.     }
  371. }
  372.  
  373. //visualize: This is the main visualization function
  374. void visualize(void)
  375. {
  376.     colorbar();
  377.  
  378.     int i, j, idx;
  379.     fftw_real  wn = (fftw_real)(winWidth - 300) / (fftw_real)(DIM + 1);   // Grid cell width
  380.     fftw_real  hn = (fftw_real)winHeight / (fftw_real)(DIM + 1);  // Grid cell heigh    printf(winWidth);
  381.  
  382.  
  383.     if (draw_smoke)
  384.     {
  385.         //int idx0, idx1, idx2, idx3;
  386.         //double px0, py0, px1, py1, px2, py2, px3, py3;
  387.         glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  388.        
  389.  
  390.         glBegin(GL_TRIANGLES);
  391.         for (j = 0; j < DIM - 1; j++)            //draw smoke
  392.         {
  393.            
  394.             for (i = 0; i < DIM - 1; i++)
  395.             {
  396.                
  397.                 px0 = wn + (fftw_real)i * wn;
  398.                 py0 = hn + (fftw_real)j * hn;
  399.                 idx0 = (j * DIM) + i;
  400.  
  401.  
  402.                 px1 = wn + (fftw_real)i * wn;
  403.                 py1 = hn + (fftw_real)(j + 1) * hn;
  404.                 idx1 = ((j + 1) * DIM) + i;
  405.                 px2 = wn + (fftw_real)(i + 1) * wn;
  406.                 py2 = hn + (fftw_real)(j + 1) * hn;
  407.                 idx2 = ((j + 1) * DIM) + (i + 1);
  408.  
  409.  
  410.                 px3 = wn + (fftw_real)(i + 1) * wn;
  411.                 py3 = hn + (fftw_real)j * hn;
  412.                 idx3 = (j * DIM) + (i + 1);
  413.  
  414.                 set_colormap(rho[idx0]);    
  415.                                             glVertex2f(px0, py0);
  416.                 set_colormap(rho[idx1]);  
  417.                                             glVertex2f(px1, py1);
  418.                 set_colormap(rho[idx2]);    
  419.                                             glVertex2f(px2, py2);
  420.  
  421.                 set_colormap(rho[idx0]);    
  422.                                             glVertex2f(px0, py0);
  423.                 set_colormap(rho[idx2]);    
  424.                                             glVertex2f(px2, py2);
  425.                 set_colormap(rho[idx3]);    
  426.                                             glVertex2f(px3, py3);
  427.    
  428.             }
  429.         }
  430.         glEnd();
  431.     }
  432.  
  433.     if (draw_vecs)
  434.     {
  435.         if (cone == 1)
  436.         {
  437.             glBegin(GL_QUADS);              //draw velocities
  438.             for (i = 0; i < DIM; i++)
  439.                 for (j = 0; j < DIM; j++)
  440.                 {
  441.                     idx = (j * DIM) + i;
  442.                     direction_to_color(vx[idx], vy[idx], color_dir);
  443.                     glVertex2f(wn + (fftw_real)i * wn, hn + (fftw_real)j * hn);
  444.                     glVertex2f((wn + (fftw_real)i * wn) + vec_scale * vx[idx], (hn + (fftw_real)j * hn) + vec_scale * vy[idx]);
  445.                 }
  446.             glEnd();
  447.         }
  448.         else if(cone==0)
  449.         {
  450.             glBegin(GL_LINES);              //draw velocities
  451.             for (i = 0; i < DIM; i++)
  452.                 for (j = 0; j < DIM; j++)
  453.                 {
  454.                     idx = (j * DIM) + i;
  455.                     direction_to_color(vx[idx], vy[idx], color_dir);
  456.                     glVertex2f(wn + (fftw_real)i * wn, hn + (fftw_real)j * hn);
  457.                     glVertex2f((wn + (fftw_real)i * wn) + vec_scale * vx[idx], (hn + (fftw_real)j * hn) + vec_scale * vy[idx]);
  458.                 }
  459.             glEnd();
  460.         }
  461.     }
  462. }
  463.  
  464.  
  465. //------ INTERACTION CODE STARTS HERE -----------------------------------------------------------------
  466.  
  467. //display: Handle window redrawing events. Simply delegates to visualize().
  468. void display(void)
  469. {
  470.     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  471.     glMatrixMode(GL_MODELVIEW);
  472.     glLoadIdentity();
  473.     visualize();
  474.  
  475.     glFlush();
  476.     glutSwapBuffers();
  477. }
  478.  
  479. //reshape: Handle window resizing (reshaping) events
  480. void reshape(int w, int h)
  481. {
  482.     glViewport(0.0f, 0.0f, (GLfloat)w, (GLfloat)h);
  483.     glMatrixMode(GL_PROJECTION);
  484.     glLoadIdentity();
  485.     gluOrtho2D(0.0, (GLdouble)w, 0.0, (GLdouble)h);
  486.     winWidth = w; winHeight = h;
  487. }
  488.  
  489. //keyboard: Handle key presses
  490. void keyboard(unsigned char key, int x, int y)
  491. {
  492.     switch (key)
  493.     {
  494.       case 't': dt -= 0.001; break;
  495.       case 'T': dt += 0.001; break;
  496.       case 'c': color_dir = 1 - color_dir; break;
  497.       case 'S': vec_scale *= 1.2; break;
  498.       case 's': vec_scale *= 0.8; break;
  499.       case 'V': visc *= 5; break;
  500.       case 'vy': visc *= 0.2; break;
  501.       case 'x': draw_smoke = 1 - draw_smoke;
  502.             if (draw_smoke==0) draw_vecs = 1; break;
  503.       case 'y': draw_vecs = 1 - draw_vecs;
  504.             if (draw_vecs==0) draw_smoke = 1; break;
  505.       case 'm': scalar_col++; if (scalar_col>COLOR_BANDS) scalar_col=COLOR_BLACKWHITE; break;
  506.       case 'a': frozen = 1-frozen; break;
  507.       case '3': cone = 1 - cone; break;
  508.       case 'q': exit(0);
  509.     }
  510. }
  511.  
  512.  
  513.  
  514. // drag: When the user drags with the mouse, add a force that corresponds to the direction of the mouse
  515. //       cursor movement. Also inject some new matter into the field at the mouse location.
  516. void drag(int mx, int my)
  517. {
  518.     int xi,yi,X,Y; double  dx, dy, len;
  519.     static int lmx=0,lmy=0;             //remembers last mouse location
  520.  
  521.     // Compute the array index that corresponds to the cursor location
  522.     xi = (int)clamp((double)(DIM + 1) * ((double)mx / (double)winWidth));
  523.     yi = (int)clamp((double)(DIM + 1) * ((double)(winHeight - my) / (double)winHeight));
  524.  
  525.     X = xi; Y = yi;
  526.  
  527.     if (X > (DIM - 1))  X = DIM - 1; if (Y > (DIM - 1))  Y = DIM - 1;
  528.     if (X < 0) X = 0; if (Y < 0) Y = 0;
  529.  
  530.     // Add force at the cursor location
  531.     my = winHeight - my;
  532.     dx = mx - lmx; dy = my - lmy;
  533.     len = sqrt(dx * dx + dy * dy);
  534.     if (len != 0.0) {  dx *= 0.1 / len; dy *= 0.1 / len; }
  535.     fx[Y * DIM + X] += dx;
  536.     fy[Y * DIM + X] += dy;
  537.     rho[Y * DIM + X] = 10.0f;
  538.     lmx = mx; lmy = my;
  539. }
  540.  
  541.  
  542. //main: The main program
  543. int main(int argc, char **argv)
  544. {
  545.     printf("Fluid Flow Simulation and Visualization\n");
  546.     printf("=======================================\n");
  547.     printf("Click and drag the mouse to steer the flow!\n");
  548.     printf("T/t:   increase/decrease simulation timestep\n");
  549.     printf
  550.     ("S/s:   increase/decrease hedgehog scaling\n");
  551.     printf("c:     toggle direction coloring on/off\n");
  552.     printf("V/vy:   increase decrease fluid viscosity\n");
  553.     printf("x:     toggle drawing matter on/off\n");
  554.     printf("y:     toggle drawing hedgehogs on/off\n");
  555.     printf("m:     toggle thru scalar coloring\n");
  556.     printf("a:     toggle the animation on/off\n");
  557.     printf("3:     toogle down cone on/off\n");
  558.     printf("q:     quit\n\n");
  559.  
  560.     glutInit(&argc, argv);
  561.     glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
  562.     //glutInitWindowSize(500,500);
  563.     glutInitWindowSize(800, 500);
  564.     glutCreateWindow("Real-time smoke simulation and visualization");
  565.     glutDisplayFunc(display);
  566.     glutReshapeFunc(reshape);
  567.     glutIdleFunc(do_one_simulation_step);
  568.     glutKeyboardFunc(keyboard);
  569.     glutMotionFunc(drag);
  570.  
  571.  
  572.     init_simulation(DIM);   //initialize the simulation data structures
  573.     glutMainLoop();         //calls do_one_simulation_step, keyboard, display, drag, reshape
  574.     return 0;
  575. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!