#ifdef __APPLE__ #include #else #include #endif #include #include #include #include #define WINDOW_WIDTH 70 #define WINDOW_HEIGHT 70 #define LATTICE_WIDTH 50 #define LATTICE_HEIGHT 50 #define LATTICE_DEPTH 5 #define VECTOR_SPACING 2 #define POINT_SIZE 9 #define SPACING 1 #define Q 19 #define DIM 3 double omega = 2.0; double lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][Q]; double tmp_lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][Q]; int e[Q][DIM] = {{0,0,0}, {1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}, {1, 1, 0}, {1, -1, 0}, {-1, 1, 0}, {-1, -1, 0}, {0, 1, 1}, {0, 1, -1}, {0, -1, 1}, {0, -1, -1}, {1, 0, 1}, {1, 0, -1}, {-1, 0, 1}, {-1, 0, -1}}; double weights[Q] = {1.0/3.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0}; void resolve_collisions( ) { int x,y,z, i,j; int rho; double u[DIM]; double u_squared; double u_dot_e[Q]; double f_eq; for(z = 0; z < LATTICE_DEPTH; z++) { for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { // Clear Variables rho = 0; u_squared = 0; u[0] = u[1] = u[2] = 0; for(i = 0; i < Q; i++) { u_dot_e[i] = 0; } for(i = 0; i < Q; i++) { rho += lattice[z][y][x][i]; for(j = 0; j < DIM; j++) { u[j] += e[i][j] * lattice[z][y][x][i]; } } for(j = 0; j < DIM; j++) { u[j] /= ((rho) + 0.1); u_squared += u[j] * u[j]; } for(i = 0; i < Q; i++) { for(j = 0; j < DIM; j++) { u_dot_e[i] += e[i][j] * u[j]; } } for(i = 0; i < Q; i++) { f_eq = weights[i] * rho * (1 - 1.5 * u_squared + 3 * u_dot_e[i] + 4.5 * u_dot_e[i] * u_dot_e[i]); tmp_lattice[z][y][x][i] = (1.0 - omega) * lattice[z][y][x][i] + omega * f_eq; } } } } } void propagate( ) { int x, y, z, i; int e_x, e_y, e_z; int new_x, new_y, new_z; for(z = 0; z < LATTICE_DEPTH; z++) { for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { for(i = 0; i < Q; i++) { e_x = e[i][0]; e_y = e[i][1]; e_z = e[i][2]; new_x = (x + e_x + LATTICE_WIDTH) % LATTICE_WIDTH; new_y = (y + e_y + LATTICE_HEIGHT) % LATTICE_HEIGHT; new_z = (z + e_z + LATTICE_DEPTH) % LATTICE_DEPTH; lattice[new_z][new_y][new_x][i] = tmp_lattice[z][y][x][i]; } } } } } void next( ) { resolve_collisions( ); propagate( ); } void display( ) { int x, y, z, i, j; int v_x, v_y, v_z; double slope[DIM]; glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); glColor3f(1.0, 1.0, 1.0); glBegin(GL_LINE); for(z = 0; z < LATTICE_DEPTH; z += VECTOR_SPACING) { for(x = 0; x < LATTICE_WIDTH; x += VECTOR_SPACING) { for(y = 0; y < LATTICE_HEIGHT; y += VECTOR_SPACING) { for(j = 0; j < DIM; j++) { slope[j] = 0.0; } for(v_z = z; v_z < z + VECTOR_SPACING; v_z++) { for(v_x = x; v_x < x + VECTOR_SPACING; v_x++) { for(v_y = y; v_y < y + VECTOR_SPACING; v_y++) { for(i = 0; i < Q; i++) { for(j = 0; j < DIM; j++) { slope[j] += e[i][j] * lattice[v_z][v_y][v_x][i]; } } } } } for(j = 0; j < DIM; j++) { slope[j] /= (Q * VECTOR_SPACING * VECTOR_SPACING); } glVertex3f(x, y, z); glVertex3f(x + slope[0], y + slope[1], z + slope[2]); } } } glEnd( ); /* glutSolidSphere(0.1, 100, 80);*/ glutSwapBuffers( ); } void idle( ) { next( ); glutPostRedisplay(); } void init_vars( ) { int x,y,z, i; for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { for(z = 0; z < LATTICE_DEPTH; z++) { for(i = 0; i < Q; i++) { lattice[z][y][x][i] = 0; tmp_lattice[z][y][x][i] = 0; } } } } for(x = LATTICE_HEIGHT/3; x < 3*LATTICE_WIDTH/4; x++) { for(y = LATTICE_WIDTH/3; y < 3*LATTICE_HEIGHT/4; y++) { for(z = 0; z < LATTICE_DEPTH; z++) { for(i = 0; i < Q/2; i++) { lattice[z][y][x][i] = 1; tmp_lattice[z][y][x][i] = 1; } } } } } void init_gl() { glClearColor(0.0, 0.0, 0.0, 0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); /* glFrustum(0, 2*WINDOW_WIDTH, 0, 2*WINDOW_HEIGHT, 80, -80.0);*/ glFrustum(-WINDOW_WIDTH, 2*WINDOW_WIDTH, -WINDOW_HEIGHT, 2*WINDOW_HEIGHT, LATTICE_DEPTH, -LATTICE_DEPTH); glMatrixMode(GL_MODELVIEW); glShadeModel(GL_SMOOTH); glPointSize(POINT_SIZE); glLineWidth(POINT_SIZE); } int main(int argc, char** argv) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH); glutInitWindowSize(WINDOW_WIDTH, WINDOW_HEIGHT); glutCreateWindow("D3Q19"); glEnable(GL_DEPTH_TEST); init_vars(); init_gl(); glutDisplayFunc(display); glutIdleFunc(idle); glutMainLoop(); return 0; }