00001
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 #include <stdio.h>
00069 #include <stdlib.h>
00070 #include "glut.h"
00071
00072 #define IX(i,j) ((i)+(N+2)*(j))
00073 #define SWAP(x0,x) {float *tmp = x0; x0 = x; x = tmp;}
00074 #define DIM 128
00075 #define SIZE DIM*DIM
00076 #define max(x,y) x>y?x:y
00077
00078
00088 void printArray (int N, float *x);
00089
00090
00098 void init (void);
00099
00100
00107 void display (void);
00108
00109
00115 void idle (void);
00116
00117
00138 void reshape (int x, int y);
00139
00140
00151 void keyboard (unsigned char key, int x, int y);
00152
00153
00165 void click (int button, int state, int x, int y);
00166
00167
00178 void motion (int x, int y);
00179
00180
00190 void timer (int value);
00191
00192
00205 void updateFluid (int button, int mouseMove, int state);
00206
00207
00219 void add_source (int N, float *x, float *s, float dt);
00220
00221
00240 void diffuse (int N, int b, float *x, float *x0, float diff, float dt);
00241
00242
00258 void advect (int N, int b, float *d, float *d0, float *u, float *v, float dt);
00259
00260
00274 void project (int N, float *u, float *v, float *p, float *div);
00275
00276
00290 void dens_step (int N, float *x, float *x0, float *u, float *v, float diff, float dt);
00291
00292
00306 void vel_step (int N, float *u, float *v, float *u0, float *v0, float visc, float dt);
00307
00308
00323 void set_bnd (int N, int b, float *x);
00324
00325
00326 float u[SIZE];
00327 float v[SIZE];
00328 float u_prev[SIZE];
00329 float v_prev[SIZE];
00330 float dens[SIZE];
00331 float dens_prev[SIZE];
00332 int boundary[SIZE];
00333 int wWidth = max(512,DIM);
00334 int wHeight = max(512,DIM);
00335 float cellSize;
00336 float diff = 0.5f;
00337 float dt = .2f;
00338 float visc = 0.0f;
00339 int stepToggle = 0;
00340 int runToggle = 1;
00341 int consoleToggle = 0;
00342 int fluidToggle = 0;
00343 int densityToggle = 1;
00344 int velocityToggle = 0;
00345 int gridToggle = 0;
00346 int boundaryToggle = 1;
00347 int drawBoundaryToggle = 0;
00348 int id;
00349 int lastX = 0;
00350 int lastY = 0;
00351 int currentX = 0;
00352 int currentY = 0;
00353 int clicked = 0;
00354 int leftClick = 0;
00355 int N = DIM-2;
00356 int emitterLoc = DIM/2;
00357 int frameRate = 60;
00358 int updown = 0;
00359 int lastCellX;
00360 int lastCellY;
00361 float totalDensity = 0;
00362
00363
00364 int main(int argc, char *argv[]) {
00365 int i;
00366
00367
00368 for (i=0; i<SIZE; i++)
00369 {
00370 u[i] = 0.0f;
00371 u_prev[i] = 0.0f;
00372 v[i] = 0.0f;
00373 v_prev[i] = 0.0f;
00374 dens[i] = 0.0f;
00375 dens_prev[i] = 0.0f;
00376 boundary[i] = 0;
00377 }
00378
00379
00380 glutInit(&argc, argv);
00381 glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
00382 glutInitWindowSize(wWidth, wHeight);
00383 glutCreateWindow("Compute Stable Fluids");
00384 init();
00385 glutDisplayFunc(display);
00386 glutReshapeFunc(reshape);
00387 glutIdleFunc(idle);
00388 glutKeyboardFunc(keyboard);
00389 glutMouseFunc(click);
00390 glutMotionFunc(motion);
00391 glutTimerFunc(100, timer, frameRate);
00392
00393
00394 glutMainLoop();
00395
00396 return 0;
00397 }
00398
00399 void printArray (int N, float *x)
00400 {
00401 int i, j;
00402
00403
00404 for (i=N+1; i>=0; i--)
00405 {
00406
00407 for (j=0; j<N+2; j++)
00408 {
00409 printf("%f ", x[IX(i,j)]);
00410 }
00411
00412
00413 printf("\n");
00414 }
00415 }
00416
00417 void init (void)
00418 {
00419
00420 cellSize = wWidth/N;
00421
00422
00423 id = glGenLists(1);
00424
00425
00426 glNewList(id, GL_COMPILE);
00427 {
00428 glBegin(GL_QUADS);
00429 {
00430 glVertex3f(-cellSize/2.0f, cellSize/2.0f, 0);
00431 glVertex3f(-cellSize/2.0f, -cellSize/2.0f, 0);
00432 glVertex3f(cellSize/2.0f, -cellSize/2.0f, 0);
00433 glVertex3f(cellSize/2.0f, cellSize/2.0f, 0);
00434 }
00435 glEnd();
00436 }
00437 glEndList();
00438
00439
00440 glClearColor(0x00, 0x00, 0x00, 1);
00441 }
00442
00443 void display (void)
00444 {
00445 int i, j;
00446 float c, temp;
00447
00448
00449 glClear (GL_COLOR_BUFFER_BIT);
00450
00451
00452 if(runToggle || stepToggle)
00453 {
00454
00455 vel_step (N, u, v, u_prev, v_prev, visc, dt);
00456 dens_step (N, dens, dens_prev, u, v, diff, dt);
00457
00458
00459 totalDensity = 0;
00460 for (i=N; i>=1; i--)
00461 {
00462 for (j=1; j<=N; j++)
00463 {
00464 totalDensity += dens[IX(i,j)];
00465 }
00466 }
00467
00468
00469 if (consoleToggle)
00470 {
00471
00472 printf("\n\nDensity Field\n");
00473 printArray(N, dens);
00474
00475
00476 printf("\n\nu Field\n");
00477 printArray(N, u);
00478
00479
00480 printf("\n\nv Field\n");
00481 printArray(N, v);
00482
00483
00484 printf("\n\nTotal Density: %f\n", totalDensity);
00485 }
00486
00487
00488 stepToggle = 0;
00489 }
00490
00491
00492 for (i=N; i>=1; i--)
00493 {
00494 for (j=1; j<=N; j++)
00495 {
00496
00497 if(boundaryToggle)
00498 {
00499
00500 if(boundary[IX(i,j)] == 1)
00501 glColor3ub(0x00, 0x00, 0xFF);
00502
00503 else
00504 glColor3ub(0x00, 0x00, 0x00);
00505
00506 glPushMatrix();
00507 {
00508
00509 glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);
00510
00511
00512 glCallList(id);
00513 }
00514 glPopMatrix();
00515 }
00516
00517
00518 if(densityToggle)
00519 {
00520 if(dens[IX(i,j)] > 0)
00521 {
00522
00523 c = (dens[IX(i,j)]);
00524 if(c == 0)
00525 glColor3f(0,0,0);
00526 else if(c < .33)
00527 glColor3f(0,0,c);
00528 else if(c < .66)
00529 glColor3f(0,c,0);
00530 else
00531 glColor3f(c,0,0);
00532
00533 glPushMatrix();
00534 {
00535
00536 glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);
00537
00538
00539 glCallList(id);
00540 }
00541 glPopMatrix();
00542 }
00543 }
00544
00545
00546 if(velocityToggle)
00547 {
00548
00549 glColor3ub(0xFF, 0x00, 0x00);
00550
00551 glPushMatrix();
00552 {
00553 float x, y;
00554
00555
00556 x = u[IX(i,j)] * cellSize * 10;
00557
00558
00559 y = v[IX(i,j)] * cellSize * 10;
00560
00561
00562 glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);
00563
00564
00565 glBegin(GL_LINES);
00566 {
00567 glVertex3f(0, 0, 0);
00568 glVertex3f(x, y, 0);
00569 }
00570 glEnd();
00571 }
00572 glPopMatrix();
00573 }
00574
00575
00576 if(gridToggle)
00577 {
00578
00579 glColor3ub(0x00, 0xFF, 0x00);
00580
00581 glPushMatrix();
00582 {
00583 glPushAttrib(GL_POLYGON_BIT);
00584 {
00585
00586 glPolygonMode (GL_FRONT, GL_LINE);
00587
00588
00589 glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);
00590
00591
00592 glCallList(id);
00593 }
00594 glPopAttrib();
00595 }
00596 glPopMatrix();
00597 }
00598 }
00599 }
00600
00601
00602 glutSwapBuffers();
00603 }
00604
00605 void idle (void)
00606 {
00607
00608 glutPostRedisplay();
00609 }
00610
00611
00612 void reshape (int x, int y)
00613 {
00614
00615 wWidth = x;
00616
00617
00618 wHeight = y;
00619
00620
00621 glDeleteLists(id,0);
00622
00623
00624 init();
00625
00626
00627 glViewport(0, 0, x, y);
00628 glMatrixMode(GL_PROJECTION);
00629 glLoadIdentity();
00630 glOrtho(0, wWidth, 0, wHeight, 0, 1);
00631 glMatrixMode(GL_MODELVIEW);
00632 glLoadIdentity();
00633
00634
00635 glutPostRedisplay();
00636 }
00637
00638 void keyboard (unsigned char key,int x, int y)
00639 {
00640 int i;
00641
00642 switch(key)
00643 {
00644
00645 case 27:
00646 exit(0);
00647
00648
00649 case 'c':
00650 consoleToggle = !consoleToggle;
00651 break;
00652
00653
00654 case 's':
00655 runToggle = !runToggle;
00656 break;
00657
00658
00659 case 'f':
00660 dens_prev[IX(1,emitterLoc)] = N * cellSize;
00661 v_prev[IX(1,emitterLoc)] = cellSize;
00662 break;
00663
00664
00665 case 'd':
00666 densityToggle = !densityToggle;
00667 break;
00668
00669
00670 case 'v':
00671 velocityToggle = !velocityToggle;
00672 break;
00673
00674
00675 case 'g':
00676 gridToggle = !gridToggle;
00677 break;
00678
00679
00680 case 'b':
00681 boundaryToggle = !boundaryToggle;
00682 break;
00683
00684
00685 case 'B':
00686 drawBoundaryToggle = !drawBoundaryToggle;
00687 break;
00688
00689
00690 case 'r':
00691
00692 for (i=0; i<SIZE; i++)
00693 {
00694 u[i] = 0.0f;
00695 u_prev[i] = 0.0f;
00696 v[i] = 0.0f;
00697 v_prev[i] = 0.0f;
00698 dens[i] = 0.0f;
00699 dens_prev[i] = 0.0f;
00700 boundary[i] = 0;
00701 }
00702 break;
00703
00704
00705 case 32:
00706 stepToggle = 1;
00707 break;
00708 }
00709
00710
00711 glutPostRedisplay();
00712 }
00713 void click (int button, int state, int x, int y)
00714 {
00715 updown = state;
00716
00717
00718 if(state == 0)
00719 {
00720 lastCellX = -1;
00721 lastCellY = -1;
00722 }
00723
00724
00725 lastX = currentX;
00726 lastY = currentY;
00727
00728
00729 currentX = x;
00730 currentY = y;
00731
00732
00733 updateFluid(button, 0, updown);
00734 }
00735 void motion (int x, int y)
00736 {
00737
00738 if(updown == 0)
00739 {
00740 lastX = currentX;
00741 lastY = currentY;
00742
00743 currentX = x;
00744 currentY = y;
00745
00746
00747 updateFluid(0, 1, updown);
00748 }
00749 }
00750
00751 void timer(int value) {
00752
00753 glutPostRedisplay();
00754
00755
00756 glutTimerFunc(1000/frameRate, timer, value);
00757 }
00758
00759 void updateFluid (int button, int mouseMove, int state)
00760 {
00761 float fx, fy;
00762 int nx, ny;
00763
00764
00765 fx = currentX / (float)wWidth;
00766 fy = (wHeight - currentY) / (float)wHeight;
00767 nx = (int)(fx * DIM);
00768 ny = (int)(fy * DIM);
00769
00770
00771 if(nx < 1)
00772 nx = 1;
00773 else if(nx > N)
00774 nx = N;
00775 if(ny < 1)
00776 ny = 1;
00777 else if(ny > N)
00778 ny = N;
00779
00780
00781 if(button == 2)
00782 {
00783
00784 dens_prev[IX(ny,nx)] = N * cellSize;
00785
00786
00787
00788 u_prev[IX(ny,nx-1)] = -cellSize;
00789 u_prev[IX(ny,nx+1)] = cellSize;
00790 v_prev[IX(ny-1,nx)] = -cellSize;
00791 v_prev[IX(ny+1,nx)] = cellSize;
00792 }
00793
00794
00795
00796 if(button == 0 && mouseMove && drawBoundaryToggle == 0 && updown == 0)
00797 {
00798
00799
00800 u_prev[IX(ny,nx)] += (currentX-lastX) * cellSize;
00801 v_prev[IX(ny,nx)] += (lastY-currentY) * cellSize;
00802 }
00803
00804
00805 if(button == 0 && drawBoundaryToggle && updown == 0)
00806 {
00807
00808 if(boundary[IX(ny,nx)] == 0)
00809 {
00810
00811 if(nx != lastCellX || ny != lastCellY)
00812 {
00813
00814 boundary[IX(ny,nx)] = 1;
00815
00816
00817 lastCellX = nx;
00818 lastCellY = ny;
00819 }
00820 }
00821
00822
00823 else
00824 {
00825
00826 if(lastCellX != nx || lastCellY != ny)
00827 {
00828
00829 boundary[IX(ny,nx)] = 0;
00830
00831
00832 lastCellX = nx;
00833 lastCellY = ny;
00834 }
00835 }
00836 }
00837
00838
00839 glutPostRedisplay();
00840 }
00841 void add_source (int N, float *x, float *s, float dt)
00842 {
00843 int i, size = (N+2)*(N+2);
00844
00845
00846
00847 for (i=0; i<size; i++)
00848 {
00849 x[i] += dt*s[i];
00850 }
00851 }
00852
00853 void diffuse (int N, int b, float *x, float *x0, float diff, float dt)
00854 {
00855 int i, j, k;
00856 float a = dt*diff*N*N;
00857
00858
00859 for (k=0; k<20; k++)
00860 {
00861
00862 for (i=1; i<=N; i++)
00863 {
00864 for (j=1; j<=N; j++)
00865 {
00866
00867
00868
00869 if(boundary[IX(i,j)])
00870 x[IX(i,j)] = 0.0f;
00871
00872
00873
00874
00875
00876 else
00877 x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)] +
00878 x[IX(i+1,j)] +
00879 x[IX(i,j-1)] +
00880 x[IX(i,j+1)])) / (1+4*a);
00881 }
00882 }
00883
00884
00885 set_bnd (N, b, x);
00886 }
00887 }
00888
00889 void advect (int N, int b, float *d, float *d0, float *u, float *v, float dt)
00890 {
00891 int i, j, i0, j0, i1, j1;
00892 float x, y, s0, t0, s1, t1, dt0, temp;
00893
00894
00895 dt0 = dt*N;
00896
00897
00898 for (i=1; i<=N; i++)
00899 {
00900 for (j=1; j<=N; j++)
00901 {
00902 x = i-dt0*u[IX(i,j)];
00903 y = j-dt0*v[IX(i,j)];
00904
00905 if (x<0.5f)
00906 x = 0.5f;
00907 if (x>N+0.5f)
00908 x = N+0.5f;
00909 i0 = (int)x;
00910 i1 = i0+1;
00911
00912 if (y<0.5f)
00913 y = 0.5f;
00914 if (y>N+0.5f)
00915 y = N+0.5f;
00916 j0 = (int)y;
00917 j1 = j0+1;
00918
00919 s1 = x-i0;
00920 s0 = 1-s1;
00921
00922 t1 = y-j0;
00923 t0 = 1-t1;
00924
00925
00926
00927 if(boundary[IX(i,j)])
00928 d[IX(i,j)] = 0.0f;
00929
00930
00931 else
00932 {
00933
00934 temp = s0 * (t0*d0[IX(i0,j0)] + t1*d0[IX(i0,j1)])+
00935 s1 * (t0*d0[IX(i1,j0)] + t1*d0[IX(i1,j1)]);
00936
00937
00938
00939 if(b == 1 && temp < 0 && boundary[IX(i,j-1)])
00940
00941 d[IX(i,j)] = -temp;
00942
00943
00944
00945 else if(b == 1 && temp > 0 && boundary[IX(i,j+1)])
00946
00947 d[IX(i,j)] = -temp;
00948
00949
00950
00951 else if(b == 2 && temp < 0 && boundary[IX(i-1,j)])
00952
00953 d[IX(i,j)] = -temp;
00954
00955
00956
00957 else if(b == 2 && temp > 0 && boundary[IX(i+1,j)])
00958
00959 d[IX(i,j)] = -temp;
00960
00961
00962 else
00963 d[IX(i,j)] = temp;
00964 }
00965 }
00966 }
00967
00968
00969 set_bnd (N, b, d);
00970 }
00971
00972 void project (int N, float *u, float *v, float *p, float *div)
00973 {
00974 int i, j, k;
00975 float h, temp;
00976
00977
00978 h = 1.0f/N;
00979
00980
00981 for (i=1; i<=N; i++)
00982 {
00983 for (j=1; j<=N; j++)
00984 {
00985
00986 temp = -0.5f*h*(u[IX(i+1,j)] - u[IX(i-1,j)]+
00987 v[IX(i,j+1)] - v[IX(i,j-1)]);
00988
00989
00990
00991 if(boundary[IX(i,j)])
00992 div[IX(i,j)] = 0.0f;
00993
00994
00995
00996 else if(temp < 0 && boundary[IX(i-1,j)])
00997
00998 div[IX(i,j)] = -temp;
00999
01000
01001
01002 else if(temp > 0 && boundary[IX(i+1,j)])
01003
01004 div[IX(i,j)] = -temp;
01005
01006
01007 else
01008 div[IX(i,j)] = temp;
01009
01010
01011 p[IX(i,j)] = 0.0f;
01012 }
01013 }
01014
01015
01016 set_bnd (N, 0, div);
01017 set_bnd (N, 0, p);
01018
01019
01020 for (k=0; k<20; k++)
01021 {
01022
01023 for (i=1; i<=N; i++)
01024 {
01025 for (j=1; j<=N; j++)
01026 {
01027
01028 temp = (div[IX(i,j)] + p[IX(i-1,j)] +
01029 p[IX(i+1,j)] +
01030 p[IX(i,j-1)] +
01031 p[IX(i,j+1)]) / 4;
01032
01033
01034
01035 if(boundary[IX(i,j)])
01036 p[IX(i,j)] = 0.0f;
01037
01038
01039
01040 else if(temp < 0 && boundary[IX(i,j-1)])
01041
01042 p[IX(i,j)] = -temp;
01043
01044
01045
01046 else if(temp > 0 && boundary[IX(i,j+1)])
01047
01048 p[IX(i,j)] = -temp;
01049
01050
01051 else
01052 p[IX(i,j)] = temp;
01053 }
01054 }
01055
01056
01057 set_bnd (N, 0, p);
01058 }
01059
01060
01061 for (i=1; i<=N; i++)
01062 {
01063 for (j=1; j<=N; j++)
01064 {
01065
01066
01067 if(boundary[IX(i,j)])
01068 {
01069 u[IX(i,j)] = 0.0f;
01070 v[IX(i,j)] = 0.0f;
01071 }
01072
01073
01074 else
01075 {
01076
01077 temp = 0.5f * (p[IX(i+1,j)] - p[IX(i-1,j)]) / h;
01078
01079
01080
01081 if(temp < 0 && boundary[IX(i,j-1)])
01082
01083 u[IX(i,j)] -= -temp;
01084
01085
01086
01087 else if(temp > 0 && boundary[IX(i,j+1)])
01088
01089 u[IX(i,j)] -= -temp;
01090
01091
01092 else
01093 u[IX(i,j)] -= temp;
01094
01095
01096 temp = 0.5f * (p[IX(i,j+1)] - p[IX(i,j-1)]) / h;
01097
01098
01099
01100 if(temp < 0 && boundary[IX(i-1,j)])
01101
01102 v[IX(i,j)] -= -temp;
01103
01104
01105
01106 else if(temp > 0 && boundary[IX(i+1,j)])
01107
01108 v[IX(i,j)] -= -temp;
01109
01110
01111 else
01112 v[IX(i,j)] -= temp;
01113 }
01114 }
01115 }
01116
01117
01118 set_bnd (N, 0, u);
01119 set_bnd (N, 0, v);
01120 }
01121
01122 void dens_step (int N, float *x, float *x0, float *u, float *v, float diff, float dt)
01123 {
01124
01125 add_source (N, x, x0, dt);
01126
01127
01128 SWAP (x0, x);
01129
01130
01131 diffuse (N, 1, x, x0, diff, dt);
01132
01133
01134 SWAP (x0, x);
01135
01136
01137 advect (N, 1, x, x0, u, v, dt);
01138 }
01139
01140 void vel_step (int N, float *u, float *v, float *u0, float *v0, float visc, float dt)
01141 {
01142
01143 add_source (N, u, u0, dt);
01144
01145
01146 add_source (N, v, v0, dt);
01147
01148
01149 SWAP (u0, u);
01150
01151
01152 diffuse (N, 0, u, u0, visc, dt);
01153
01154
01155 SWAP (v0, v);
01156
01157
01158 diffuse (N, 0, v, v0, visc, dt);
01159
01160
01161 project (N, u, v, u0, v0);
01162
01163
01164 SWAP (u0, u);
01165
01166
01167 SWAP (v0, v);
01168
01169
01170 advect (N, 0, u, u0, u0, v0, dt);
01171
01172
01173 advect (N, 0, v, v0, u0, v0, dt);
01174
01175
01176 project (N, u, v, u0, v0);
01177 }
01178
01179 void set_bnd (int N, int b, float *x)
01180 {
01181 int i;
01182
01183
01184 if(b == 0)
01185 {
01186 for (i=1; i<=N; i++)
01187 {
01188 x[IX(0 ,i)] = x[IX(1,i)];
01189 x[IX(N+1,i)] = x[IX(N,i)];
01190 x[IX(i,0 )] = x[IX(i,1)];
01191 x[IX(i,N+1)] = x[IX(i,N)];
01192 }
01193
01194 x[IX(0 ,0 )] = 0.5f*(x[IX(1,0 )] + x[IX(0 ,1)]);
01195 x[IX(0 ,N+1)] = 0.5f*(x[IX(1,N+1)] + x[IX(0 ,N)]);
01196 x[IX(N+1,0 )] = 0.5f*(x[IX(N,0 )] + x[IX(N+1,1)]);
01197 x[IX(N+1,N+1)] = 0.5f*(x[IX(N,N+1)] + x[IX(N+1,N)]);
01198 }
01199
01200 else
01201 {
01202 for (i=1; i<=N; i++)
01203 {
01204 x[IX(0 ,i)] = 0;
01205 x[IX(N+1,i)] = 0;
01206 x[IX(i,0 )] = 0;
01207 x[IX(i,N+1)] = 0;
01208 }
01209
01210 x[IX(0 ,0 )] = 0;
01211 x[IX(0 ,N+1)] = 0;
01212 x[IX(N+1,0 )] = 0;
01213 x[IX(N+1,N+1)] = 0;
01214 }
01215 }