API tools faq
paste
Advertisement
backstreetimrul

fld.c 5

backstreetimrul
Jul 29th, 2019
84
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 21.33 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. ///------------Density/Velocity/Force ----------------------
  13. float DensityM, DensityM1, DensityN, DensityP;
  14. float fMx, fM1x, fNx, fPx;
  15. float fMy, fM1y, fNy, fPy;
  16. int dfv = 0;
  17. //--------------------------------------------------------
  18.  
  19. //--- SIMULATION PARAMETERS ------------------------------------------------------------------------
  20. const int DIM = 50;             //size of simulation grid
  21. double dt = 0.4;                //simulation time step
  22. float visc = 0.001;             //fluid viscosity
  23. fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current moment
  24. fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment
  25. fftw_real *fx, *fy;             //(fx,fy)   = user-controlled simulation forces, steered with the mouse
  26. fftw_real *rho, *rho0;          //smoke density at the current (rho) and previous (rho0) moment
  27. rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization
  28.  
  29.  
  30. //--- VISUALIZATION PARAMETERS ---------------------------------------------------------------------
  31. int   winWidth, winHeight;      //size of the graphics window, in pixels
  32. int   color_dir = 0;            //use direction color-coding or not
  33. float vec_scale = 1000;         //scaling of hedgehogs
  34. int   draw_smoke = 0;           //draw the smoke or not
  35. int   draw_vecs = 1;            //draw the vector field or not
  36. const int COLOR_BLACKWHITE=0;   //different types of color mapping: black-and-white, rainbow, banded
  37. const int COLOR_RAINBOW=1;
  38. const int COLOR_BANDS=2;
  39. int   scalar_col = 0;           //method for scalar coloring
  40. int   frozen = 0;               //toggles on/off the animation
  41. int cone = 0; //toogle cone on/off
  42.  
  43.  
  44. //------ SIMULATION CODE STARTS HERE -----------------------------------------------------------------
  45.  
  46. //init_simulation: Initialize simulation data structures as a function of the grid size 'n'.
  47. //                 Although the simulation takes place on a 2D grid, we allocate all data structures as 1D arrays,
  48. //                 for compatibility with the FFTW numerical library.
  49. void init_simulation(int n)
  50. {
  51.     int i; size_t dim;
  52.  
  53.     dim     = n*2*(n/2+1)*sizeof(fftw_real);        //Allocate data structures
  54.     vx       = (fftw_real*) malloc(dim);
  55.     vy       = (fftw_real*) malloc(dim);
  56.     vx0      = (fftw_real*) malloc(dim);
  57.     vy0      = (fftw_real*) malloc(dim);
  58.     dim     = n * n * sizeof(fftw_real);
  59.     fx      = (fftw_real*) malloc(dim);
  60.     fy      = (fftw_real*) malloc(dim);
  61.     rho     = (fftw_real*) malloc(dim);
  62.     rho0    = (fftw_real*) malloc(dim);
  63.     plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE);
  64.     plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE);
  65.  
  66.     for (i = 0; i < n * n; i++)                      //Initialize data structures to 0
  67.     { vx[i] = vy[i] = vx0[i] = vy0[i] = fx[i] = fy[i] = rho[i] = rho0[i] = 0.0f; }
  68. }
  69.  
  70.  
  71. //FFT: Execute the Fast Fourier Transform on the dataset 'vx'.
  72. //     'dirfection' indicates if we do the direct (1) or inverse (-1) Fourier Transform
  73. void FFT(int direction,void* vx)
  74. {
  75.     if(direction==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real*)vx,(fftw_complex*)vx);
  76.     else             rfftwnd_one_complex_to_real(plan_cr,(fftw_complex*)vx,(fftw_real*)vx);
  77. }
  78.  
  79. int clamp(float x)
  80. { return ((x)>=0.0?((int)(x)):(-((int)(1-(x))))); }
  81.  
  82. //solve: Solve (compute) one step of the fluid flow simulation
  83. void solve(int n, fftw_real* vx, fftw_real* vy, fftw_real* vx0, fftw_real* vy0, fftw_real visc, fftw_real dt)
  84. {
  85.     fftw_real x, y, x0, y0, f, r, U[2], V[2], s, t;
  86.     int i, j, i0, j0, i1, j1;
  87.  
  88.     for (i=0;i<n*n;i++)
  89.     { vx[i] += dt*vx0[i]; vx0[i] = vx[i]; vy[i] += dt*vy0[i]; vy0[i] = vy[i]; }
  90.  
  91.     for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n )
  92.        for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n )
  93.        {
  94.           x0 = n*(x-dt*vx0[i+n*j])-0.5f;
  95.           y0 = n*(y-dt*vy0[i+n*j])-0.5f;
  96.           i0 = clamp(x0); s = x0-i0;
  97.           i0 = (n+(i0%n))%n;
  98.           i1 = (i0+1)%n;
  99.           j0 = clamp(y0); t = y0-j0;
  100.           j0 = (n+(j0%n))%n;
  101.           j1 = (j0+1)%n;
  102.           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]);
  103.           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]);
  104.        }
  105.  
  106.     for(i=0; i<n; i++)
  107.       for(j=0; j<n; j++)
  108.       {  vx0[i+(n+2)*j] = vx[i+n*j]; vy0[i+(n+2)*j] = vy[i+n*j]; }
  109.  
  110.     FFT(1,vx0);
  111.     FFT(1,vy0);
  112.  
  113.     for (i=0;i<=n;i+=2)
  114.     {
  115.        x = 0.5f*i;
  116.        for (j=0;j<n;j++)
  117.        {
  118.           y = j<=n/2 ? (fftw_real)j : (fftw_real)j-n;
  119.           r = x*x+y*y;
  120.           if ( r==0.0f ) continue;
  121.           f = (fftw_real)exp(-r*dt*visc);
  122.           U[0] = vx0[i  +(n+2)*j]; V[0] = vy0[i  +(n+2)*j];
  123.           U[1] = vx0[i+1+(n+2)*j]; V[1] = vy0[i+1+(n+2)*j];
  124.  
  125.           vx0[i  +(n+2)*j] = f*((1-x*x/r)*U[0]     -x*y/r *V[0]);
  126.           vx0[i+1+(n+2)*j] = f*((1-x*x/r)*U[1]     -x*y/r *V[1]);
  127.           vy0[i+  (n+2)*j] = f*(  -y*x/r *U[0] + (1-y*y/r)*V[0]);
  128.           vy0[i+1+(n+2)*j] = f*(  -y*x/r *U[1] + (1-y*y/r)*V[1]);
  129.        }
  130.     }
  131.  
  132.     FFT(-1,vx0);
  133.     FFT(-1,vy0);
  134.  
  135.     f = 1.0/(n*n);
  136.     for (i=0;i<n;i++)
  137.        for (j=0;j<n;j++)
  138.        { vx[i+n*j] = f*vx0[i+(n+2)*j]; vy[i+n*j] = f*vy0[i+(n+2)*j]; }
  139. }
  140.  
  141.  
  142. // diffuse_matter: This function diffuses matter that has been placed in the velocity field. It's almost identical to the
  143. // velocity diffusion step in the function above. The input matter densities are in rho0 and the result is written into rho.
  144. void diffuse_matter(int n, fftw_real *vx, fftw_real *vy, fftw_real *rho, fftw_real *rho0, fftw_real dt)
  145. {
  146.     fftw_real x, y, x0, y0, s, t;
  147.     int i, j, i0, j0, i1, j1;
  148.     float maximum = rho[0], minimum = rho[0];
  149.  
  150.     for (x = 0.5f / n, i = 0; i < n; i++, x += 1.0f / n)
  151.         for (y = 0.5f / n, j = 0; j < n; j++, y += 1.0f / n)
  152.         {
  153.             x0 = n * (x - dt * vx[i + n * j]) - 0.5f;
  154.             y0 = n * (y - dt * vy[i + n * j]) - 0.5f;
  155.             i0 = clamp(x0);
  156.             s = x0 - i0;
  157.             i0 = (n + (i0 % n)) % n;
  158.             i1 = (i0 + 1) % n;
  159.             j0 = clamp(y0);
  160.             t = y0 - j0;
  161.             j0 = (n + (j0 % n)) % n;
  162.             j1 = (j0 + 1) % n;
  163.             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]);
  164.             if (rho[i + n * j] > maximum) {
  165.                 maximum = rho[i + n * j];
  166.                 DensityP = maximum;
  167.             }
  168.             if (rho[i + n * j] < minimum) {
  169.                 minimum = rho[i + n * j];
  170.                 DensityN = minimum;
  171.             }
  172.             DensityM = ((DensityN + DensityP) / 2);
  173.             DensityM1 = ((DensityM + DensityN) / 2);
  174.  
  175.         }
  176. }
  177.  
  178. //set_forces: copy user-controlled forces to the force vectors that are sent to the solver.
  179. //            Also dampen forces and matter density to get a stable simulation.
  180. void set_forces(void)
  181. {
  182.     int i;
  183.     float mafx = fx[0], mifx = fx[0], mafy = fy[0], mify = fy[0];
  184.     for (i = 0; i < DIM * DIM; i++)
  185.     {
  186.  
  187.         /*if (rho[i]>maximum){
  188.             maximum=rho[i];
  189.             DensityP=maximum;
  190.         }
  191.         if (rho[i]<minimum){
  192.             minimum=rho[i];s
  193.             DensityN=minimum;
  194.         }*/
  195.         rho0[i] = 0.995 * rho[i];
  196.         fx[i] *= 0.85;
  197.         fy[i] *= 0.85;
  198.         vx0[i] = fx[i];
  199.         vy0[i] = fy[i];
  200.  
  201.         if (fx[i] > mafx) {
  202.             mafx = fx[i];
  203.             fPx = mafx;
  204.         }
  205.         if (fx[i] < mifx) {
  206.             mifx = fx[i];
  207.             fNx = mifx;
  208.         }
  209.         if (fy[i] > mafy) {
  210.             mafy = fy[i];
  211.             fPy = mafy;
  212.         }
  213.         if (fy[i] < mify) {
  214.             mify = fy[i];
  215.             fNy = mify;
  216.         }
  217.         fMx = ((fNx + fPx) / 2);
  218.         fMy = ((fNy + fPy) / 2);
  219.  
  220.         fM1x = ((fNx + fMx) / 2);
  221.         fM1y = ((fNy + fMy) / 2);
  222.         /*DensityM=((DensityN+DensityP)/2);
  223.         DensityM1=((DensityM+DensityN)/2);*/
  224.     }
  225. }
  226.  
  227.  
  228. //do_one_simulation_step: Do one complete cycle of the simulation:
  229. //      - set_forces:
  230. //      - solve:            read forces from the user
  231. //      - diffuse_matter:   compute a new set of velocities
  232. //      - gluPostRedisplay: draw a new visualization frame
  233. void do_one_simulation_step(void)
  234. {
  235.     if (!frozen)
  236.     {
  237.       set_forces();
  238.       solve(DIM, vx, vy, vx0, vy0, visc, dt);
  239.       diffuse_matter(DIM, vx, vy, rho, rho0, dt);
  240.       glutPostRedisplay();
  241.     }
  242. }
  243.  
  244.  
  245. //------ VISUALIZATION CODE STARTS HERE -----------------------------------------------------------------
  246.  
  247.  
  248. //rainbow: Implements a color palette, mapping the scalar 'value' to a rainbow color RGB
  249. void rainbow(float value,float* R,float* G,float* B)
  250. {
  251.    const float dx=0.8;
  252.    if (value<0) value=0; if (value>1) value=1;
  253.    value = (6-2*dx)*value+dx;
  254.    *R = max(0.0,(3-fabs(value-4)-fabs(value-5))/2);
  255.    *G = max(0.0,(4-fabs(value-2)-fabs(value-4))/2);
  256.    *B = max(0.0,(3-fabs(value-1)-fabs(value-2))/2);
  257. }
  258.  
  259. //set_colormap: Sets three different types of colormaps
  260. void set_colormap(float vy)
  261. {
  262.    float R,G,B;
  263.  
  264.    if (scalar_col==COLOR_BLACKWHITE)
  265.        R = G = B = vy;
  266.    else if (scalar_col==COLOR_RAINBOW)
  267.        rainbow(vy,&R,&G,&B);
  268.    else if (scalar_col==COLOR_BANDS)
  269.        {
  270.           const int NLEVELS = 7;
  271.           vy *= NLEVELS; vy = (int)(vy); vy/= NLEVELS;
  272.           rainbow(vy,&R,&G,&B);
  273.        }
  274.    glColor3f(R,G,B);
  275. }
  276.  
  277.  
  278. //direction_to_color: Set the current color by mapping a direction vector (x,y), using
  279. //                    the color mapping method 'method'. If method==1, map the vector direction
  280. //                    using a rainbow colormap. If method==0, simply use the white color
  281. void direction_to_color(float x, float y, int method)
  282. {
  283.     //float r,g,b,f;
  284.     if (method)
  285.     {
  286.       f = atan2(y,x) / 3.1415927 + 1;
  287.       r = f;
  288.       if(r > 1) r = 2 - r;
  289.       g = f + .66667;
  290.       if(g > 2) g -= 2;
  291.       if(g > 1) g = 2 - g;
  292.       b = f + 2 * .66667;
  293.       if(b > 2) b -= 2;
  294.       if(b > 1) b = 2 - b;
  295.     }
  296.     else
  297.     { r = g = b = 1; }
  298.     glColor3f(r,g,b);
  299. }
  300.  
  301. //Function to write in graphics window
  302. void drawString(void* font, float x, float y, char* str)
  303. {
  304.     /* Draws string ’str’ in font ’font’, at world (x,y,0) */
  305.     char* ch;
  306.     glRasterPos3f(x, y, 0.0);
  307.     for (ch = str; *ch; ch++)
  308.         glutBitmapCharacter(font, (int)* ch);
  309. }
  310.  
  311. void colorbar()
  312. {
  313.     //drawString(GLUT_BITMAP_HELVETICA_12, 600, 420, fx);
  314.     if (color_dir==1)
  315.     {
  316.         glBegin(GL_QUADS);
  317.         glColor3f(r, g, b);
  318.         glColor3f(0.58, 0.0, 0.82);
  319.         //glColor3f(1.0, 0.0, 0.0);
  320.         glVertex2f(600, 200);
  321.         glVertex2f(650, 200);
  322.         glColor3f(r, g, b);
  323.         glColor3f(0.3, 0.0, 0.51);
  324.         //glColor3f(0.0, 0.0, 1.0);
  325.         glVertex2f(650, 233.33);
  326.         glVertex2f(600, 233.33);
  327.         glColor3f(r, g, b);
  328.         glColor3f(0.3, 0.0, 0.51);
  329.         //glColor3f(0.0, 0.0, 1.0);
  330.         glVertex2f(650, 233.33);
  331.         glVertex2f(600, 233.33);
  332.         glColor3f(r, g, b);
  333.         glColor3f(0.0, 0.0, 1.0);
  334.         //glColor3f(0.0, 0.0, 1.0);
  335.         glVertex2f(600, 266.66);
  336.         glVertex2f(650, 266.66);
  337.         glColor3f(r, g, b);
  338.         glColor3f(0.0, 0.0, 1.0);
  339.         //glColor3f(0.0, 0.0, 1.0);
  340.         glVertex2f(600, 266.66);
  341.         glVertex2f(650, 266.66);
  342.         glColor3f(r, g, b);
  343.         glColor3f(0.0, 1.0, 0.0);
  344.         //glColor3f(0.0, 0.0, 1.0);
  345.         glVertex2f(650, 299.99);
  346.         glVertex2f(600, 299.99);
  347.         glColor3f(r, g, b);
  348.         glColor3f(0.0, 1.0, 0.0);
  349.         //glColor3f(0.0, 0.0, 1.0);
  350.         glVertex2f(650, 299.99);
  351.         glVertex2f(600, 299.99);
  352.         glColor3f(r, g, b);
  353.         glColor3f(1.0, 1.0, 0.0);
  354.         //glColor3f(0.0, 0.0, 1.0);
  355.         glVertex2f(600, 333.33);
  356.         glVertex2f(650, 333.33);
  357.         glColor3f(r, g, b);
  358.         glColor3f(1.0, 1.0, 0.0);
  359.         //glColor3f(0.0, 0.0, 1.0);
  360.         glVertex2f(600, 333.33);
  361.         glVertex2f(650, 333.33);
  362.         glColor3f(r, g, b);
  363.         glColor3f(1.0, 0.5, 0.0);
  364.         //glColor3f(0.0, 0.0, 1.0);
  365.         glVertex2f(650, 366.66);
  366.         glVertex2f(600, 366.66);
  367.         glColor3f(r, g, b);
  368.         glColor3f(1.0, 0.5, 0.0);
  369.         //glColor3f(0.0, 0.0, 1.0);
  370.         glVertex2f(650, 366.66);
  371.         glVertex2f(600, 366.66);
  372.         glColor3f(r, g, b);
  373.         glColor3f(1.0, 0.0, 0.0);
  374.         //glColor3f(0.0, 0.0, 1.0);
  375.         glVertex2f(600, 400);
  376.         glVertex2f(650, 400);
  377.  
  378.         glEnd();
  379.         glFlush();
  380.     }
  381.     else if(color_dir==0)
  382.     {
  383.         glBegin(GL_QUADS);
  384.         glColor3f(1.0, 1.0, 1.0);
  385.         glVertex2f(600, 200);
  386.         glVertex2f(650, 200);
  387.         //glColor3f(0.0, 0.0, 1.0);
  388.         glVertex2f(650, 233.33);
  389.         glVertex2f(600, 233.33);
  390.         //glColor3f(0.0, 0.0, 1.0);
  391.         glVertex2f(650, 233.33);
  392.         glVertex2f(600, 233.33);
  393.         //glColor3f(0.0, 0.0, 1.0);
  394.         glVertex2f(600, 266.66);
  395.         glVertex2f(650, 266.66);
  396.         //glColor3f(0.0, 0.0, 1.0);
  397.         glVertex2f(600, 266.66);
  398.         glVertex2f(650, 266.66);
  399.         //glColor3f(0.0, 0.0, 1.0);
  400.         glVertex2f(650, 299.99);
  401.         glVertex2f(600, 299.99);
  402.         //glColor3f(0.0, 0.0, 1.0);
  403.         glVertex2f(650, 299.99);
  404.         glVertex2f(600, 299.99);
  405.         //glColor3f(0.0, 0.0, 1.0);
  406.         glVertex2f(600, 333.33);
  407.         glVertex2f(650, 333.33);
  408.         //glColor3f(0.0, 0.0, 1.0);
  409.         glVertex2f(600, 333.33);
  410.         glVertex2f(650, 333.33);
  411.         //glColor3f(0.0, 0.0, 1.0);
  412.         glVertex2f(650, 366.66);
  413.         glVertex2f(600, 366.66);
  414.         //glColor3f(0.0, 0.0, 1.0);
  415.         glVertex2f(650, 366.66);
  416.         glVertex2f(600, 366.66);
  417.         //glColor3f(0.0, 0.0, 1.0);
  418.         glVertex2f(600, 400);
  419.         glVertex2f(650, 400);
  420.  
  421.         glEnd();
  422.         glFlush();
  423.     }
  424. }
  425.  
  426. //visualize: This is the main visualization function
  427. void visualize(void)
  428. {
  429.     colorbar();
  430.  
  431.     int i, j, idx;
  432.     fftw_real  wn = (fftw_real)(winWidth - 300) / (fftw_real)(DIM + 1);   // Grid cell width
  433.     fftw_real  hn = (fftw_real)winHeight / (fftw_real)(DIM + 1);  // Grid cell heigh    printf(winWidth);
  434.  
  435.  
  436.     if (draw_smoke)
  437.     {
  438.         //int idx0, idx1, idx2, idx3;
  439.         //double px0, py0, px1, py1, px2, py2, px3, py3;
  440.         glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  441.        
  442.  
  443.         glBegin(GL_TRIANGLES);
  444.         for (j = 0; j < DIM - 1; j++)            //draw smoke
  445.         {
  446.            
  447.             for (i = 0; i < DIM - 1; i++)
  448.             {
  449.                
  450.                 px0 = wn + (fftw_real)i * wn;
  451.                 py0 = hn + (fftw_real)j * hn;
  452.                 idx0 = (j * DIM) + i;
  453.  
  454.  
  455.                 px1 = wn + (fftw_real)i * wn;
  456.                 py1 = hn + (fftw_real)(j + 1) * hn;
  457.                 idx1 = ((j + 1) * DIM) + i;
  458.                 px2 = wn + (fftw_real)(i + 1) * wn;
  459.                 py2 = hn + (fftw_real)(j + 1) * hn;
  460.                 idx2 = ((j + 1) * DIM) + (i + 1);
  461.  
  462.  
  463.                 px3 = wn + (fftw_real)(i + 1) * wn;
  464.                 py3 = hn + (fftw_real)j * hn;
  465.                 idx3 = (j * DIM) + (i + 1);
  466.  
  467.                 set_colormap(rho[idx0]);    
  468.                                             glVertex2f(px0, py0);
  469.                 set_colormap(rho[idx1]);  
  470.                                             glVertex2f(px1, py1);
  471.                 set_colormap(rho[idx2]);    
  472.                                             glVertex2f(px2, py2);
  473.  
  474.                 set_colormap(rho[idx0]);    
  475.                                             glVertex2f(px0, py0);
  476.                 set_colormap(rho[idx2]);    
  477.                                             glVertex2f(px2, py2);
  478.                 set_colormap(rho[idx3]);    
  479.                                             glVertex2f(px3, py3);
  480.    
  481.             }
  482.         }
  483.         glEnd();
  484.     }
  485.  
  486.     if (draw_vecs)
  487.     {
  488.         if (cone == 1)
  489.         {
  490.             glBegin(GL_QUADS);              //draw velocities
  491.             for (i = 0; i < DIM; i++)
  492.                 for (j = 0; j < DIM; j++)
  493.                 {
  494.                     idx = (j * DIM) + i;
  495.                     direction_to_color(vx[idx], vy[idx], color_dir);
  496.                     glVertex2f(wn + (fftw_real)i * wn, hn + (fftw_real)j * hn);
  497.                     glVertex2f((wn + (fftw_real)i * wn) + vec_scale * vx[idx], (hn + (fftw_real)j * hn) + vec_scale * vy[idx]);
  498.                 }
  499.             glEnd();
  500.         }
  501.         else if(cone==0)
  502.         {
  503.             glBegin(GL_LINES);              //draw velocities
  504.             for (i = 0; i < DIM; i++)
  505.                 for (j = 0; j < DIM; j++)
  506.                 {
  507.                     idx = (j * DIM) + i;
  508.                     direction_to_color(vx[idx], vy[idx], color_dir);
  509.                     glVertex2f(wn + (fftw_real)i * wn, hn + (fftw_real)j * hn);
  510.                     glVertex2f((wn + (fftw_real)i * wn) + vec_scale * vx[idx], (hn + (fftw_real)j * hn) + vec_scale * vy[idx]);
  511.                 }
  512.             glEnd();
  513.         }
  514.     }
  515. }
  516.  
  517.  
  518. //------ INTERACTION CODE STARTS HERE -----------------------------------------------------------------
  519.  
  520. //display: Handle window redrawing events. Simply delegates to visualize().
  521. void range(void)
  522. {
  523.     char* str = "Density(rho)";
  524.     char* strr = "Force";
  525.     char Str1[100];
  526.     char Str2[100];
  527.     char Str3[100];
  528.     char Str4[100];
  529.  
  530.     char Strr1[100];
  531.     char Strr2[100];
  532.     char Strr3[100];
  533.     char Strr4[100];
  534.  
  535.     char Strrr1[100];
  536.     char Strrr2[100];
  537.     char Strrr3[100];
  538.     char Strrr4[100];
  539.  
  540.  
  541.     ///***
  542.     if(dfv==0)
  543.     {
  544.         glColor3f(1.0, 1.0, 1.0);
  545.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 605, 420, str);
  546.         glEnd();
  547.         glFlush();
  548.  
  549.         glColor3f(1.0, 1.0, 1.0);
  550.         gcvt(DensityP, 4, Str1);
  551.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 400, Str1);
  552.         glEnd();
  553.         glFlush();
  554.  
  555.         glColor3f(1.0, 1.0, 1.0);
  556.         gcvt(DensityM, 4, Str2);
  557.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 330, Str2);
  558.         glEnd();
  559.         glFlush();
  560.  
  561.         glColor3f(1.0, 1.0, 1.0);
  562.         gcvt(DensityM1, 4, Str3);
  563.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 260, Str3);
  564.         glEnd();
  565.         glFlush();
  566.  
  567.         glColor3f(1.0, 1.0, 1.0);
  568.         gcvt(DensityN, 4, Str4);
  569.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 200, Str4);
  570.         glEnd();
  571.         glFlush();
  572.     }
  573.     //***/
  574.  
  575.     ///***
  576.     else if(dfv==1)
  577.     {
  578.         glColor3f(1.0, 1.0, 1.0);
  579.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 605, 420, strr);
  580.         glEnd();
  581.         glFlush();
  582.  
  583.         glColor3f(1.0, 1.0, 1.0);
  584.         gcvt(fPx, 4, Strr1);
  585.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 400, Strr1);
  586.         glEnd();
  587.         glFlush();
  588.  
  589.         glColor3f(1.0, 1.0, 1.0);
  590.         gcvt(fPy, 4, Strrr1);
  591.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 360, Strrr1);
  592.         glEnd();
  593.         glFlush();
  594.  
  595.         glColor3f(1.0, 1.0, 1.0);
  596.         gcvt(fMx, 4, Strr2);
  597.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 330, Strr2);
  598.         glEnd();
  599.         glFlush();
  600.  
  601.         glColor3f(1.0, 1.0, 1.0);
  602.         gcvt(fMy, 4, Strrr2);
  603.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 300, Strrr2);
  604.         glEnd();
  605.         glFlush();
  606.  
  607.         glColor3f(1.0, 1.0, 1.0);
  608.         gcvt(fM1x, 4, Strr3);
  609.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 279, Strr3);
  610.         glEnd();
  611.         glFlush();
  612.  
  613.         glColor3f(1.0, 1.0, 1.0);
  614.         gcvt(fM1y, 4, Strrr3);
  615.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 246, Strrr3);
  616.         glEnd();
  617.         glFlush();
  618.  
  619.         glColor3f(1.0, 1.0, 1.0);
  620.         gcvt(fNx, 4, Strr4);
  621.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 226, Strr4);
  622.         glEnd();
  623.         glFlush();
  624.  
  625.         glColor3f(1.0, 1.0, 1.0);
  626.         gcvt(fNy, 4, Strrr4);
  627.         drawString(GLUT_BITMAP_TIMES_ROMAN_10, 670, 200, Strrr4);
  628.         glEnd();
  629.         glFlush();
  630.     }
  631.     //***/
  632. }
  633. ///***
  634. void display(void)
  635. {
  636.     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  637.     glMatrixMode(GL_MODELVIEW);
  638.     glLoadIdentity();
  639.     visualize();
  640.     glFlush();
  641.  
  642.     range();
  643.  
  644.     glutSwapBuffers();
  645. }
  646.  
  647. //***/
  648.  
  649. //reshape: Handle window resizing (reshaping) events
  650. void reshape(int w, int h)
  651. {
  652.     glViewport(0.0f, 0.0f, (GLfloat)w, (GLfloat)h);
  653.     glMatrixMode(GL_PROJECTION);
  654.     glLoadIdentity();
  655.     gluOrtho2D(0.0, (GLdouble)w, 0.0, (GLdouble)h);
  656.     winWidth = w; winHeight = h;
  657. }
  658.  
  659. //keyboard: Handle key presses
  660. void keyboard(unsigned char key, int x, int y)
  661. {
  662.     switch (key)
  663.     {
  664.       case 't': dt -= 0.001; break;
  665.       case 'T': dt += 0.001; break;
  666.       case 'c': color_dir = 1 - color_dir; break;
  667.       case 'S': vec_scale *= 1.2; break;
  668.       case 's': vec_scale *= 0.8; break;
  669.       case 'V': visc *= 5; break;
  670.       case 'vy': visc *= 0.2; break;
  671.       case 'x': draw_smoke = 1 - draw_smoke;
  672.             if (draw_smoke==0) draw_vecs = 1; break;
  673.       case 'y': draw_vecs = 1 - draw_vecs;
  674.             if (draw_vecs==0) draw_smoke = 1; break;
  675.       case 'm': scalar_col++; if (scalar_col>COLOR_BANDS) scalar_col=COLOR_BLACKWHITE; break;
  676.       case 'a': frozen = 1-frozen; break;
  677.       case '3': cone = 1 - cone; break;
  678.       case '4': dfv = 1 - dfv; break;
  679.       case 'q': exit(0);
  680.     }
  681. }
  682.  
  683.  
  684.  
  685. // drag: When the user drags with the mouse, add a force that corresponds to the direction of the mouse
  686. //       cursor movement. Also inject some new matter into the field at the mouse location.
  687. void drag(int mx, int my)
  688. {
  689.     int xi,yi,X,Y; double  dx, dy, len;
  690.     static int lmx=0,lmy=0;             //remembers last mouse location
  691.  
  692.     // Compute the array index that corresponds to the cursor location
  693.     xi = (int)clamp((double)(DIM + 1) * ((double)mx / (double)winWidth));
  694.     yi = (int)clamp((double)(DIM + 1) * ((double)(winHeight - my) / (double)winHeight));
  695.  
  696.     X = xi; Y = yi;
  697.  
  698.     if (X > (DIM - 1))  X = DIM - 1; if (Y > (DIM - 1))  Y = DIM - 1;
  699.     if (X < 0) X = 0; if (Y < 0) Y = 0;
  700.  
  701.     // Add force at the cursor location
  702.     my = winHeight - my;
  703.     dx = mx - lmx; dy = my - lmy;
  704.     len = sqrt(dx * dx + dy * dy);
  705.     if (len != 0.0) {  dx *= 0.1 / len; dy *= 0.1 / len; }
  706.     fx[Y * DIM + X] += dx;
  707.     fy[Y * DIM + X] += dy;
  708.     rho[Y * DIM + X] = 10.0f;
  709.     lmx = mx; lmy = my;
  710. }
  711.  
  712.  
  713. //main: The main program
  714. int main(int argc, char **argv)
  715. {
  716.     printf("Fluid Flow Simulation and Visualization\n");
  717.     printf("=======================================\n");
  718.     printf("Click and drag the mouse to steer the flow!\n");
  719.     printf("T/t:   increase/decrease simulation timestep\n");
  720.     printf
  721.     ("S/s:   increase/decrease hedgehog scaling\n");
  722.     printf("c:     toggle direction coloring on/off\n");
  723.     printf("V/vy:   increase decrease fluid viscosity\n");
  724.     printf("x:     toggle drawing matter on/off\n");
  725.     printf("y:     toggle drawing hedgehogs on/off\n");
  726.     printf("m:     toggle thru scalar coloring\n");
  727.     printf("a:     toggle the animation on/off\n");
  728.     printf("3:     toogle down cone on/off\n");
  729.     printf("4:     toggle among datasets and divergence: rho/||f||");
  730.     printf("q:     quit\n\n");
  731.  
  732.     glutInit(&argc, argv);
  733.     glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
  734.     //glutInitWindowSize(500,500);
  735.     glutInitWindowSize(800, 500);
  736.     glutCreateWindow("Real-time smoke simulation and visualization");
  737.     glutDisplayFunc(display);
  738.     glutReshapeFunc(reshape);
  739.     glutIdleFunc(do_one_simulation_step);
  740.     glutKeyboardFunc(keyboard);
  741.     glutMotionFunc(drag);
  742.  
  743.  
  744.     init_simulation(DIM);   //initialize the simulation data structures
  745.     glutMainLoop();         //calls do_one_simulation_step, keyboard, display, drag, reshape
  746.     return 0;
  747. }
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!