diff options
Diffstat (limited to 'Smoke/fluids.c')
| -rw-r--r-- | Smoke/fluids.c | 162 |
1 files changed, 129 insertions, 33 deletions
diff --git a/Smoke/fluids.c b/Smoke/fluids.c index 6071d36..0dc8469 100644 --- a/Smoke/fluids.c +++ b/Smoke/fluids.c @@ -19,9 +19,10 @@ #include "fluids.h" #include "smoke.h" #include "colormap.h" +#include "glyphs.h" +#include "heightplots.h" //--- SIMULATION PARAMETERS ------------------------------------------------------------------------ -const int DIM = 50; //size of simulation grid int var_dims = 25; double dt = 0.4; //simulation time step float visc = 0.001f; //fluid viscosity @@ -49,53 +50,55 @@ float threshold2 = 2.0f; int isolines_nr = 1; static int fluids_calculate = TRUE; +static int DIM; //size of simulation grid //------ SIMULATION CODE STARTS HERE ----------------------------------------------------------------- //init_simulation: Initialize simulation data structures as a function of the grid size 'n'. // Although the simulation takes place on a 2D grid, we allocate all data structures as 1D arrays, // for compatibility with the FFTW numerical library. -fftw_real *init_simulation(int n) +fftw_real *init_simulation(int n, struct vis_data_arrays *vis_data) { - int i; size_t dim; + int i; size_t dim1, dim2; fftw_real *return_value; + + DIM = n; - dim = n * 2*(n/2+1)*sizeof(fftw_real); //Allocate data structures - return_value = (fftw_real *)malloc(dim); - vx = (fftw_real*) malloc(dim); - vy = (fftw_real*) malloc(dim); - vx0 = (fftw_real*) malloc(dim); - vy0 = (fftw_real*) malloc(dim); - dim = n * n * sizeof(fftw_real); - fx = (fftw_real*) malloc(dim); - fy = (fftw_real*) malloc(dim); - rho = (fftw_real*) malloc(dim); - rho0 = (fftw_real*) malloc(dim); + dim1 = n * 2*(n/2+1)*sizeof(fftw_real); //Allocate data structures + return_value = (fftw_real *)malloc(dim1); + vx = (fftw_real*) malloc(dim1); + vy = (fftw_real*) malloc(dim1); + vx0 = (fftw_real*) malloc(dim1); + vy0 = (fftw_real*) malloc(dim1); + dim2 = n * n * sizeof(fftw_real); + fx = (fftw_real*) malloc(dim2); + fy = (fftw_real*) malloc(dim2); + rho = (fftw_real*) malloc(dim2); + rho0 = (fftw_real*) malloc(dim2); plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE); plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE); - height_array = (fftw_real*) malloc(dim); - normal_array = (struct point*) malloc(dim); - frame_hist = (fftw_real*) malloc(dim * n); + height_array = (fftw_real*) malloc(dim2); + normal_array = (struct point*) malloc(dim2); + frame_hist = (fftw_real*) malloc(dim2 * n); for (i = 0; i <= n; i++) { - frame_hist[i] = (fftw_real*) malloc(dim); + frame_hist[i] = (fftw_real*) malloc(dim1); } for (i = 0; i < n * n; i++) //Initialize data structures to 0 { vx[i] = vy[i] = vx0[i] = vy0[i] = fx[i] = fy[i] = rho[i] = rho0[i] = height_array[i] = 0.0f; } - return return_value; -} + vis_data->rho = rho; + vis_data->force = (fftw_real *)malloc(dim2); + vis_data->vel = (fftw_real *)malloc(dim1); + vis_data->div_force = (fftw_real *)malloc(dim2); + vis_data->div_vel = (fftw_real *)malloc(dim1); + vis_data->height = (fftw_real *)malloc(dim2); + vis_data->normals = (struct point *)malloc(dim2); -void fluids_init(int dim) -{ - size_t frame_size; - - frame_size = dim * 2*(dim/2+1)*sizeof(fftw_real); - - smoke_init(frame_size); + return return_value; } int rescale_to_winwidth(float value) @@ -286,7 +289,7 @@ void calculate_normal_vectors(void) } -void calculate_height_plot(void) +void calculate_height_plots(void) { int i; for (i = 0; i < DIM * DIM; i++) @@ -294,10 +297,10 @@ void calculate_height_plot(void) switch (vis_dataset) { case DATASET_FORCE: - height_array[i] = quake_root(fx[i] * fx[i] + fy[i] * fy[i]) * 100; + height_array[i] = vec_len2f(fx[i], fy[i]) * 100; break; case DATASET_VEL: - height_array[i] = quake_root(vx[i] * vx[i] + vy[i] * vy[i]) * 5000; + height_array[i] = vec_len2f(vx[i], vy[i]) * 5000; break; case DATASET_RHO: height_array[i] = rho[i] * 16; @@ -311,11 +314,88 @@ void calculate_height_plot(void) } } +void calculate_height_plot(struct vis_data_arrays *vis_data, int dataset, int index) +{ + switch(dataset) + { + default: + case DATASET_RHO: + vis_data->height[index] = rho[index]; + break; + case DATASET_VEL: + vis_data->height[index] = vec_len2f(vx[index], vy[index]); + break; + case DATASET_FORCE: + vis_data->height[index] = vec_len2f(fx[index], fy[index]); + break; + } +} + void copy_frames(fftw_real *dataset) { memcpy(smoke_get_frame(), dataset, DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real)); } +float par_der(float xy1, float xy2, fftw_real cell) +{ + return (xy1 - (xy2/ cell)); +} + +void populate_arrays(struct vis_data_arrays *vis_data) +{ + int idx, i, j; + + fftw_real wn = (fftw_real)winWidth / (fftw_real)(DIM + 1); // Grid cell width + fftw_real hn = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell height + + for (j = 0; j < DIM; j++) + { + for (i = 0; i < DIM; i++) + { + idx = (j * DIM) + i; + + vis_data->rho[idx] = rho[idx]; + vis_data->vel[idx] = vec_len2f(vx[idx], vy[idx]); + vis_data->force[idx] = vec_len2f(fx[idx], fy[idx]); + vis_data->div_vel[idx] = par_der(vx[idx +1], vx[idx], wn) + +par_der(vy[idx +1], vy[idx], winHeight); + vis_data->div_force[idx] = par_der(fx[idx +1], fx[idx], hn) + +par_der(fy[idx +1], fy[idx], winHeight); + + calculate_height_plot(vis_data, heightplots_get_dataset(), idx); + } + } +} + +fftw_real *get_frame(struct vis_data_arrays *vis_data, int dataset) +{ + fftw_real *return_value; + + switch(dataset) + { + default: + case DATASET_RHO: + return_value = vis_data->rho; + break; + case DATASET_VEL: + return_value = vis_data->vel; + break; + case DATASET_FORCE: + return_value = vis_data->force; + break; + case DATASET_DIVV: + return_value = vis_data->div_vel; + break; + case DATASET_DIVF: + return_value = vis_data->div_force; + break; + case DATASET_HIST: + return_value = vis_data->history[0]; + } + + return return_value; +} + //do_one_simulation_step: Do one complete cycle of the simulation: // - set_forces: // - solve: read forces from the user @@ -323,17 +403,28 @@ void copy_frames(fftw_real *dataset) // - gluPostRedisplay: draw a new visualization frame -void calculate_one_simulation_step(fftw_real *field) +void calculate_one_simulation_step(struct vis_data_arrays *vis_data) { + int dataset; + fftw_real *frame; + if (fluids_calculate) { set_forces(); solve(DIM, vx, vy, vx0, vy0, visc, dt); diffuse_matter(DIM, vx, vy, rho, rho0, dt); - calculate_height_plot(); + calculate_height_plots(); calculate_normal_vectors(); - copy_frames(rho); colormap_autoscale(); } + populate_arrays(vis_data); + + dataset = smoke_get_dataset(); + frame = get_frame(vis_data, dataset); + smoke_set_frame(frame); + + dataset = glyphs_get_dataset_color(); + frame = get_frame(vis_data, dataset); + glyphs_set_frame(frame); } //------ VISUALIZATION CODE STARTS HERE ----------------------------------------------------------------- @@ -491,3 +582,8 @@ int fluids_get_calculate(void) { return fluids_calculate; } + +int fluids_get_dim(void) +{ + return DIM; +} |