776 lines
26 KiB
C
776 lines
26 KiB
C
/*********************************************************************************/
|
|
/* */
|
|
/* Animation of heat equation in a planar domain */
|
|
/* */
|
|
/* N. Berglund, May 2021 */
|
|
/* */
|
|
/* Feel free to reuse, but if doing so it would be nice to drop a */
|
|
/* line to nils.berglund@univ-orleans.fr - Thanks! */
|
|
/* */
|
|
/* compile with */
|
|
/* gcc -o heat heat.c */
|
|
/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */
|
|
/* */
|
|
/* To make a video, set MOVIE to 1 and create subfolder tif_heat */
|
|
/* It may be possible to increase parameter PAUSE */
|
|
/* */
|
|
/* create movie using */
|
|
/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */
|
|
/* */
|
|
/*********************************************************************************/
|
|
|
|
/*********************************************************************************/
|
|
/* */
|
|
/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
|
|
/* One could make it much faster by using a GPU */
|
|
/* */
|
|
/*********************************************************************************/
|
|
|
|
#include <math.h>
|
|
#include <string.h>
|
|
#include <GL/glut.h>
|
|
#include <GL/glu.h>
|
|
#include <unistd.h>
|
|
#include <sys/types.h>
|
|
#include <tiffio.h> /* Sam Leffler's libtiff library. */
|
|
#include <omp.h>
|
|
|
|
#define MOVIE 1 /* set to 1 to generate movie */
|
|
|
|
/* General geometrical parameters */
|
|
|
|
// #define WINWIDTH 1280 /* window width */
|
|
#define WINWIDTH 720 /* window width */
|
|
#define WINHEIGHT 720 /* window height */
|
|
|
|
// #define NX 1280 /* number of grid points on x axis */
|
|
#define NX 720 /* number of grid points on x axis */
|
|
#define NY 720 /* number of grid points on y axis */
|
|
// #define NX 640 /* number of grid points on x axis */
|
|
// #define NY 360 /* number of grid points on y axis */
|
|
|
|
/* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
|
|
/* but will multiply run time by 4 */
|
|
|
|
// #define XMIN -2.0
|
|
// #define XMAX 2.0 /* x interval */
|
|
#define XMIN -1.125
|
|
#define XMAX 1.125 /* x interval */
|
|
#define YMIN -1.125
|
|
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
|
|
|
|
#define JULIA_SCALE 0.5 /* scaling for Julia sets */
|
|
|
|
/* Choice of the billiard table */
|
|
|
|
#define B_DOMAIN 27 /* choice of domain shape, see list in global_pdes.c */
|
|
|
|
#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */
|
|
|
|
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
|
|
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
|
|
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
|
|
|
|
#define LAMBDA 1.1 /* parameter controlling the dimensions of domain */
|
|
#define MU 0.8 /* parameter controlling the dimensions of domain */
|
|
#define NPOLY 6 /* number of sides of polygon */
|
|
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
|
|
#define MDEPTH 6 /* depth of computation of Menger gasket */
|
|
#define MRATIO 5 /* ratio defining Menger gasket */
|
|
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
|
|
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
|
|
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
|
|
#define NGRIDX 15 /* number of grid point for grid of disks */
|
|
#define NGRIDY 20 /* number of grid point for grid of disks */
|
|
|
|
#define X_SHOOTER -0.2
|
|
#define Y_SHOOTER -0.6
|
|
#define X_TARGET 0.4
|
|
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
|
|
|
|
#define ISO_XSHIFT_LEFT -1.65
|
|
#define ISO_XSHIFT_RIGHT 0.4
|
|
#define ISO_YSHIFT_LEFT -0.05
|
|
#define ISO_YSHIFT_RIGHT -0.05
|
|
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
|
|
|
|
/* You can add more billiard tables by adapting the functions */
|
|
/* xy_in_billiard and draw_billiard in sub_wave.c */
|
|
|
|
/* Physical patameters of wave equation */
|
|
|
|
// #define DT 0.00001
|
|
#define DT 0.000004
|
|
// #define DT 0.000002
|
|
// #define DT 0.00000002
|
|
// #define DT 0.000000005
|
|
#define VISCOSITY 10.0
|
|
#define T_OUT 2.0 /* outside temperature */
|
|
#define T_IN 0.0 /* inside temperature */
|
|
// #define T_OUT 0.0 /* outside temperature */
|
|
// #define T_IN 2.0 /* inside temperature */
|
|
#define SPEED 0.0 /* speed of drift to the right */
|
|
|
|
/* Boundary conditions, see list in global_pdes.c */
|
|
|
|
#define B_COND 1
|
|
|
|
/* Parameters for length and speed of simulation */
|
|
|
|
#define NSTEPS 1200 /* number of frames of movie */
|
|
#define NVID 30 /* number of iterations between images displayed on screen */
|
|
// #define NVID 100 /* number of iterations between images displayed on screen */
|
|
#define NSEG 100 /* number of segments of boundary */
|
|
#define BOUNDARY_WIDTH 1 /* width of billiard boundary */
|
|
#define DRAW_BILLIARD 0 /* set to 1 to draw billiard boundary */
|
|
|
|
#define PAUSE 100 /* number of frames after which to pause */
|
|
#define PSLEEP 1 /* sleep time during pause */
|
|
#define SLEEP1 2 /* initial sleeping time */
|
|
#define SLEEP2 1 /* final sleeping time */
|
|
|
|
/* For debugging purposes only */
|
|
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
|
|
#define VMAX 10.0 /* max value of wave amplitude */
|
|
|
|
/* Field representation */
|
|
|
|
#define FIELD_REP 0
|
|
|
|
#define F_INTENSITY 0 /* color represents intensity */
|
|
#define F_GRADIENT 1 /* color represents norm of gradient */
|
|
|
|
#define DRAW_FIELD_LINES 0 /* set to 1 to draw field lines */
|
|
#define FIELD_LINE_WIDTH 1 /* width of field lines */
|
|
#define N_FIELD_LINES 120 /* number of field lines */
|
|
#define FIELD_LINE_FACTOR 120 /* factor controlling precision when computing origin of field lines */
|
|
|
|
/* Color schemes, see list in global_pdes.c */
|
|
|
|
#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */
|
|
|
|
#define BLACK 1 /* black background */
|
|
|
|
#define COLOR_SCHEME 1 /* choice of color scheme */
|
|
|
|
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
|
|
// #define SLOPE 0.1 /* sensitivity of color on wave amplitude */
|
|
#define SLOPE 0.2 /* sensitivity of color on wave amplitude */
|
|
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
|
|
|
|
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
|
|
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
|
|
#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
|
|
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
|
|
// #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */
|
|
// #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */
|
|
#define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */
|
|
#define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */
|
|
// #define HUEMEAN 270.0 /* mean value of hue for color scheme C_HUE */
|
|
// #define HUEAMP -130.0 /* amplitude of variation of hue for color scheme C_HUE */
|
|
#define E_SCALE 100.0 /* scaling factor for energy representation */
|
|
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
|
|
#define LOG_SHIFT 0.0
|
|
|
|
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
|
|
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
|
|
#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */
|
|
#define ROTATE_COLOR_SCHEME 1 /* set to 1 to draw color scheme horizontally */
|
|
|
|
#include "global_pdes.c"
|
|
#include "sub_wave.c"
|
|
|
|
|
|
double courant2; /* Courant parameter squared */
|
|
double dx2; /* spatial step size squared */
|
|
double intstep; /* integration step */
|
|
double intstep1; /* integration step used in absorbing boundary conditions */
|
|
|
|
|
|
|
|
void init_gaussian(double x, double y, double mean, double amplitude, double scalex,
|
|
double *phi[NX], short int * xy_in[NX])
|
|
/* initialise field with gaussian at position (x,y) */
|
|
{
|
|
int i, j, in;
|
|
double xy[2], dist2, module, phase, scale2;
|
|
|
|
scale2 = scalex*scalex;
|
|
printf("Initialising field\n");
|
|
for (i=0; i<NX; i++)
|
|
for (j=0; j<NY; j++)
|
|
{
|
|
ij_to_xy(i, j, xy);
|
|
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
|
|
|
|
in = xy_in[i][j];
|
|
if (in == 1)
|
|
{
|
|
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
|
|
module = amplitude*exp(-dist2/scale2);
|
|
if (module < 1.0e-15) module = 1.0e-15;
|
|
|
|
phi[i][j] = mean + module/scalex;
|
|
} /* boundary temperatures */
|
|
else if (in >= 2) phi[i][j] = T_IN*pow(0.75, (double)(in-2));
|
|
// else if (in >= 2) phi[i][j] = T_IN*pow(1.0 - 0.5*(double)(in-2), (double)(in-2));
|
|
// else if (in >= 2) phi[i][j] = T_IN*(1.0 - (double)(in-2)/((double)MDEPTH))*(1.0 - (double)(in-2)/((double)MDEPTH));
|
|
else phi[i][j] = T_OUT;
|
|
}
|
|
}
|
|
|
|
void init_julia_set(double *phi[NX], short int * xy_in[NX])
|
|
/* change Julia set boundary condition */
|
|
{
|
|
int i, j, in;
|
|
double xy[2], dist2, module, phase, scale2;
|
|
|
|
// printf("Changing Julia set\n");
|
|
for (i=0; i<NX; i++)
|
|
for (j=0; j<NY; j++)
|
|
{
|
|
ij_to_xy(i, j, xy);
|
|
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
|
|
|
|
in = xy_in[i][j];
|
|
if (in >= 2) phi[i][j] = T_IN;
|
|
}
|
|
}
|
|
|
|
|
|
/*********************/
|
|
/* animation part */
|
|
/*********************/
|
|
|
|
|
|
void compute_gradient(double *phi[NX], double *nablax[NX], double *nablay[NX])
|
|
/* compute the gradient of the field */
|
|
{
|
|
int i, j, iplus, iminus, jplus, jminus;
|
|
double dx;
|
|
|
|
dx = (XMAX-XMIN)/((double)NX);
|
|
|
|
for (i=0; i<NX; i++)
|
|
for (j=0; j<NY; j++)
|
|
{
|
|
iplus = i+1; if (iplus == NX) iplus = NX-1;
|
|
iminus = i-1; if (iminus == -1) iminus = 0;
|
|
jplus = j+1; if (jplus == NX) jplus = NY-1;
|
|
jminus = j-1; if (jminus == -1) jminus = 0;
|
|
nablax[i][j] = (phi[iplus][j] - phi[iminus][j])/dx;
|
|
nablay[i][j] = (phi[i][jplus] - phi[i][jminus])/dx;
|
|
}
|
|
}
|
|
|
|
void draw_field_line(double x, double y, short int *xy_in[NX], double *nablax[NX],
|
|
double *nablay[NX], double delta, int nsteps)
|
|
/* draw a field line of the gradient, starting in (x,y) */
|
|
{
|
|
double x1, y1, x2, y2, pos[2], nabx, naby, norm2, norm;
|
|
int i = 0, ij[2], cont = 1;
|
|
|
|
glColor3f(1.0, 1.0, 1.0);
|
|
// glColor3f(0.0, 0.0, 0.0);
|
|
glLineWidth(FIELD_LINE_WIDTH);
|
|
x1 = x;
|
|
y1 = y;
|
|
|
|
// printf("Drawing field line \n");
|
|
|
|
glEnable(GL_LINE_SMOOTH);
|
|
glBegin(GL_LINE_STRIP);
|
|
xy_to_pos(x1, y1, pos);
|
|
glVertex2d(pos[0], pos[1]);
|
|
|
|
i = 0;
|
|
while ((cont)&&(i < nsteps))
|
|
{
|
|
xy_to_ij(x1, y1, ij);
|
|
|
|
if (ij[0] < 0) ij[0] = 0;
|
|
if (ij[0] > NX-1) ij[0] = NX-1;
|
|
if (ij[1] < 0) ij[1] = 0;
|
|
if (ij[1] > NY-1) ij[1] = NY-1;
|
|
|
|
nabx = nablax[ij[0]][ij[1]];
|
|
naby = nablay[ij[0]][ij[1]];
|
|
|
|
norm2 = nabx*nabx + naby*naby;
|
|
|
|
if (norm2 > 1.0e-14)
|
|
{
|
|
/* avoid too large step size */
|
|
if (norm2 < 1.0e-9) norm2 = 1.0e-9;
|
|
norm = sqrt(norm2);
|
|
x1 = x1 + delta*nabx/norm;
|
|
y1 = y1 + delta*naby/norm;
|
|
}
|
|
else cont = 0;
|
|
|
|
if (!xy_in[ij[0]][ij[1]]) cont = 0;
|
|
|
|
/* stop if the boundary is hit */
|
|
// if (xy_in[ij[0]][ij[1]] != 1) cont = 0;
|
|
|
|
// printf("x1 = %.3lg \t y1 = %.3lg \n", x1, y1);
|
|
|
|
xy_to_pos(x1, y1, pos);
|
|
glVertex2d(pos[0], pos[1]);
|
|
|
|
i++;
|
|
}
|
|
glEnd();
|
|
}
|
|
|
|
void draw_wave(double *phi[NX], short int *xy_in[NX], double scale, int time)
|
|
/* draw the field */
|
|
{
|
|
int i, j, iplus, iminus, jplus, jminus, ij[2], counter = 0;
|
|
static int first = 1;
|
|
double rgb[3], xy[2], x1, y1, x2, y2, dx, value, angle, dangle, intens, deltaintens, sum = 0.0;
|
|
double *nablax[NX], *nablay[NX];
|
|
static double linex[N_FIELD_LINES*FIELD_LINE_FACTOR], liney[N_FIELD_LINES*FIELD_LINE_FACTOR], distance[N_FIELD_LINES*FIELD_LINE_FACTOR], integral[N_FIELD_LINES*FIELD_LINE_FACTOR + 1];
|
|
|
|
for (i=0; i<NX; i++)
|
|
{
|
|
nablax[i] = (double *)malloc(NY*sizeof(double));
|
|
nablay[i] = (double *)malloc(NY*sizeof(double));
|
|
}
|
|
|
|
/* compute the gradient */
|
|
compute_gradient(phi, nablax, nablay);
|
|
|
|
/* compute the position of origins of field lines */
|
|
if ((first)&&(DRAW_FIELD_LINES))
|
|
{
|
|
first = 0;
|
|
|
|
printf("computing linex\n");
|
|
|
|
// x1 = LAMBDA + MU*1.01;
|
|
// y1 = 1.0;
|
|
x1 = 0.99*LAMBDA;
|
|
y1 = 0.0;
|
|
linex[0] = x1;
|
|
liney[0] = y1;
|
|
dangle = DPI/((double)(N_FIELD_LINES*FIELD_LINE_FACTOR));
|
|
|
|
for (i = 1; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
|
|
{
|
|
angle = (double)i*dangle;
|
|
// x2 = LAMBDA + MU*1.01*cos(angle);
|
|
// y2 = 0.5 + MU*1.01*sin(angle);
|
|
x2 = 0.99*LAMBDA*cos(angle);
|
|
y2 = 0.99*LAMBDA*sin(angle);
|
|
linex[i] = x2;
|
|
liney[i] = y2;
|
|
distance[i-1] = module2(x2-x1,y2-y1);
|
|
x1 = x2;
|
|
y1 = y2;
|
|
}
|
|
distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2- 0.99*LAMBDA,y2);
|
|
// distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-LAMBDA,y2-0.5);
|
|
}
|
|
|
|
dx = (XMAX-XMIN)/((double)NX);
|
|
glBegin(GL_QUADS);
|
|
|
|
for (i=0; i<NX; i++)
|
|
for (j=0; j<NY; j++)
|
|
{
|
|
if (FIELD_REP == F_INTENSITY) value = phi[i][j];
|
|
else if (FIELD_REP == F_GRADIENT)
|
|
{
|
|
value = module2(nablax[i][j], nablay[i][j]);
|
|
}
|
|
|
|
if (xy_in[i][j] == 1)
|
|
{
|
|
color_scheme(COLOR_SCHEME, value, scale, time, rgb);
|
|
glColor3f(rgb[0], rgb[1], rgb[2]);
|
|
}
|
|
else glColor3f(0.0, 0.0, 0.0);
|
|
|
|
glVertex2i(i, j);
|
|
glVertex2i(i+1, j);
|
|
glVertex2i(i+1, j+1);
|
|
glVertex2i(i, j+1);
|
|
}
|
|
glEnd ();
|
|
|
|
/* draw a field line */
|
|
if (DRAW_FIELD_LINES)
|
|
{
|
|
/* compute gradient norm along boundary and its integral */
|
|
for (i = 0; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
|
|
{
|
|
xy_to_ij(linex[i], liney[i], ij);
|
|
intens = module2(nablax[ij[0]][ij[1]], nablay[ij[0]][ij[1]])*distance[i];
|
|
if (i > 0) integral[i] = integral[i-1] + intens;
|
|
else integral[i] = intens;
|
|
}
|
|
deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
|
|
|
|
// printf("delta = %.5lg\n", deltaintens);
|
|
|
|
i = 0;
|
|
draw_field_line(linex[0], liney[0], xy_in, nablax, nablay, 0.00002, 100000);
|
|
for (j = 1; j < N_FIELD_LINES+1; j++)
|
|
{
|
|
while ((integral[i] <= j*deltaintens)&&(i < N_FIELD_LINES*FIELD_LINE_FACTOR)) i++;
|
|
draw_field_line(linex[i], liney[i], xy_in, nablax, nablay, 0.00002, 100000);
|
|
counter++;
|
|
}
|
|
printf("%i lines\n", counter);
|
|
}
|
|
|
|
|
|
for (i=0; i<NX; i++)
|
|
{
|
|
free(nablax[i]);
|
|
free(nablay[i]);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
void evolve_wave_half(double *phi_in[NX], double *phi_out[NX], short int *xy_in[NX])
|
|
/* time step of field evolution */
|
|
{
|
|
int i, j, iplus, iminus, jplus, jminus;
|
|
double delta1, delta2, x, y;
|
|
|
|
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
|
|
for (i=0; i<NX; i++){
|
|
for (j=0; j<NY; j++){
|
|
if (xy_in[i][j] == 1){
|
|
/* discretized Laplacian depending on boundary conditions */
|
|
if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
|
|
{
|
|
iplus = (i+1); if (iplus == NX) iplus = NX-1;
|
|
iminus = (i-1); if (iminus == -1) iminus = 0;
|
|
jplus = (j+1); if (jplus == NY) jplus = NY-1;
|
|
jminus = (j-1); if (jminus == -1) jminus = 0;
|
|
}
|
|
else if (B_COND == BC_PERIODIC)
|
|
{
|
|
iplus = (i+1) % NX;
|
|
iminus = (i-1) % NX;
|
|
if (iminus < 0) iminus += NX;
|
|
jplus = (j+1) % NY;
|
|
jminus = (j-1) % NY;
|
|
if (jminus < 0) jminus += NY;
|
|
}
|
|
|
|
delta1 = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
|
|
|
|
x = phi_in[i][j];
|
|
|
|
/* evolve phi */
|
|
if (B_COND != BC_ABSORBING)
|
|
{
|
|
phi_out[i][j] = x + intstep*(delta1 - SPEED*(phi_in[iplus][j] - phi_in[i][j]));
|
|
}
|
|
else /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
|
|
{
|
|
/* in the bulk */
|
|
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
|
|
{
|
|
phi_out[i][j] = x - intstep*delta2;
|
|
}
|
|
/* right border */
|
|
else if (i==NX-1)
|
|
{
|
|
phi_out[i][j] = x - intstep1*(x - phi_in[i-1][j]);
|
|
}
|
|
/* upper border */
|
|
else if (j==NY-1)
|
|
{
|
|
phi_out[i][j] = x - intstep1*(x - phi_in[i][j-1]);
|
|
}
|
|
/* left border */
|
|
else if (i==0)
|
|
{
|
|
phi_out[i][j] = x - intstep1*(x - phi_in[1][j]);
|
|
}
|
|
/* lower border */
|
|
else if (j==0)
|
|
{
|
|
phi_out[i][j] = x - intstep1*(x - phi_in[i][1]);
|
|
}
|
|
}
|
|
|
|
|
|
if (FLOOR)
|
|
{
|
|
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
|
|
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
|
|
}
|
|
|
|
void evolve_wave(double *phi[NX], double *phi_tmp[NX], short int *xy_in[NX])
|
|
/* time step of field evolution */
|
|
{
|
|
evolve_wave_half(phi, phi_tmp, xy_in);
|
|
evolve_wave_half(phi_tmp, phi, xy_in);
|
|
}
|
|
|
|
|
|
|
|
|
|
double compute_variance(double *phi[NX], short int * xy_in[NX])
|
|
/* compute the variance (total probability) of the field */
|
|
{
|
|
int i, j, n = 0;
|
|
double variance = 0.0;
|
|
|
|
for (i=1; i<NX; i++)
|
|
for (j=1; j<NY; j++)
|
|
{
|
|
if (xy_in[i][j])
|
|
{
|
|
n++;
|
|
variance += phi[i][j]*phi[i][j];
|
|
}
|
|
}
|
|
if (n==0) n=1;
|
|
return(variance/(double)n);
|
|
}
|
|
|
|
void renormalise_field(double *phi[NX], short int * xy_in[NX], double variance)
|
|
/* renormalise variance of field */
|
|
{
|
|
int i, j;
|
|
double stdv;
|
|
|
|
stdv = sqrt(variance);
|
|
|
|
for (i=1; i<NX; i++)
|
|
for (j=1; j<NY; j++)
|
|
{
|
|
if (xy_in[i][j])
|
|
{
|
|
phi[i][j] = phi[i][j]/stdv;
|
|
}
|
|
}
|
|
}
|
|
|
|
void print_level(int level)
|
|
{
|
|
double pos[2];
|
|
char message[50];
|
|
|
|
glColor3f(1.0, 1.0, 1.0);
|
|
sprintf(message, "Level %i", level);
|
|
xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
|
|
write_text(pos[0], pos[1], message);
|
|
}
|
|
|
|
|
|
void print_Julia_parameters()
|
|
{
|
|
double pos[2];
|
|
char message[50];
|
|
|
|
glColor3f(1.0, 1.0, 1.0);
|
|
if (julia_y >= 0.0) sprintf(message, "c = %.5f + %.5f i", julia_x, julia_y);
|
|
else sprintf(message, "c = %.5f %.5f i", julia_x, julia_y);
|
|
xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
|
|
write_text(pos[0], pos[1], message);
|
|
}
|
|
|
|
void set_Julia_parameters(int time, double *phi[NX], short int *xy_in[NX])
|
|
{
|
|
double jangle, cosj, sinj, radius = 0.15;
|
|
|
|
jangle = (double)time*DPI/(double)NSTEPS;
|
|
// jangle = (double)time*0.001;
|
|
// jangle = (double)time*0.0001;
|
|
|
|
cosj = cos(jangle);
|
|
sinj = sin(jangle);
|
|
julia_x = -0.9 + radius*cosj;
|
|
julia_y = radius*sinj;
|
|
init_julia_set(phi, xy_in);
|
|
|
|
printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
|
|
}
|
|
|
|
void set_Julia_parameters_cardioid(int time, double *phi[NX], short int *xy_in[NX])
|
|
{
|
|
double jangle, cosj, sinj, yshift;
|
|
|
|
jangle = pow(1.05 + (double)time*0.00003, 0.333);
|
|
yshift = 0.02*sin((double)time*PID*0.002);
|
|
// jangle = pow(1.0 + (double)time*0.00003, 0.333);
|
|
// jangle = pow(0.05 + (double)time*0.00003, 0.333);
|
|
// jangle = pow(0.1 + (double)time*0.00001, 0.333);
|
|
// yshift = 0.04*sin((double)time*PID*0.002);
|
|
|
|
cosj = cos(jangle);
|
|
sinj = sin(jangle);
|
|
julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
|
|
julia_y = 0.5*sinj*(1.0-cosj) + yshift;
|
|
// julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
|
|
// julia_y = 0.5*sinj*(1.0-cosj);
|
|
init_julia_set(phi, xy_in);
|
|
|
|
printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
|
|
}
|
|
|
|
void animation()
|
|
{
|
|
double time, scale, dx, var, jangle, cosj, sinj;
|
|
double *phi[NX], *phi_tmp[NX];
|
|
short int *xy_in[NX];
|
|
int i, j, s;
|
|
|
|
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
|
|
for (i=0; i<NX; i++)
|
|
{
|
|
phi[i] = (double *)malloc(NY*sizeof(double));
|
|
phi_tmp[i] = (double *)malloc(NY*sizeof(double));
|
|
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
|
|
}
|
|
|
|
npolyline = init_polyline(MDEPTH, polyline);
|
|
for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
|
|
|
|
dx = (XMAX-XMIN)/((double)NX);
|
|
intstep = DT/(dx*dx*VISCOSITY);
|
|
intstep1 = DT/(dx*VISCOSITY);
|
|
|
|
// julia_x = 0.1;
|
|
// julia_y = 0.6;
|
|
|
|
// set_Julia_parameters(0, phi, xy_in);
|
|
|
|
printf("Integration step %.3lg\n", intstep);
|
|
|
|
/* initialize wave wave function */
|
|
init_gaussian(-1.0, 0.0, 0.1, 0.0, 0.01, phi, xy_in);
|
|
// init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
|
|
|
|
if (SCALE)
|
|
{
|
|
var = compute_variance(phi, xy_in);
|
|
scale = sqrt(1.0 + var);
|
|
renormalise_field(phi, xy_in, var);
|
|
}
|
|
|
|
blank();
|
|
glColor3f(0.0, 0.0, 0.0);
|
|
|
|
|
|
glutSwapBuffers();
|
|
|
|
draw_wave(phi, xy_in, 1.0, 0);
|
|
if (DRAW_BILLIARD) draw_billiard(0, 1.0);
|
|
// print_Julia_parameters(i);
|
|
|
|
// print_level(MDEPTH);
|
|
|
|
glutSwapBuffers();
|
|
|
|
sleep(SLEEP1);
|
|
if (MOVIE) for (i=0; i<SLEEP1*25; i++) save_frame();
|
|
|
|
for (i=0; i<=NSTEPS; i++)
|
|
{
|
|
/* compute the variance of the field to adjust color scheme */
|
|
/* the color depends on the field divided by sqrt(1 + variance) */
|
|
if (SCALE)
|
|
{
|
|
var = compute_variance(phi, xy_in);
|
|
scale = sqrt(1.0 + var);
|
|
// printf("Norm: %5lg\t Scaling factor: %5lg\n", var, scale);
|
|
renormalise_field(phi, xy_in, var);
|
|
}
|
|
else scale = 1.0;
|
|
|
|
draw_wave(phi, xy_in, scale, i);
|
|
|
|
for (j=0; j<NVID; j++) evolve_wave(phi, phi_tmp, xy_in);
|
|
|
|
if (DRAW_BILLIARD) draw_billiard(0, 1.0);
|
|
|
|
// print_level(MDEPTH);
|
|
// print_Julia_parameters(i);
|
|
|
|
glutSwapBuffers();
|
|
|
|
/* modify Julia set */
|
|
// set_Julia_parameters(i, phi, xy_in);
|
|
|
|
if (MOVIE)
|
|
{
|
|
save_frame();
|
|
|
|
/* it seems that saving too many files too fast can cause trouble with the file system */
|
|
/* so this is to make a pause from time to time - parameter PAUSE may need adjusting */
|
|
if (i % PAUSE == PAUSE - 1)
|
|
{
|
|
printf("Making a short pause\n");
|
|
sleep(PSLEEP);
|
|
s = system("mv wave*.tif tif_heat/");
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
if (MOVIE)
|
|
{
|
|
for (i=0; i<20; i++) save_frame();
|
|
s = system("mv wave*.tif tif_heat/");
|
|
}
|
|
for (i=0; i<NX; i++)
|
|
{
|
|
free(phi[i]);
|
|
free(phi_tmp[i]);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
void display(void)
|
|
{
|
|
glPushMatrix();
|
|
|
|
blank();
|
|
glutSwapBuffers();
|
|
blank();
|
|
glutSwapBuffers();
|
|
|
|
animation();
|
|
sleep(SLEEP2);
|
|
|
|
glPopMatrix();
|
|
|
|
glutDestroyWindow(glutGetWindow());
|
|
|
|
}
|
|
|
|
|
|
int main(int argc, char** argv)
|
|
{
|
|
glutInit(&argc, argv);
|
|
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
|
|
glutInitWindowSize(WINWIDTH,WINHEIGHT);
|
|
glutCreateWindow("Heat equation in a planar domain");
|
|
|
|
init();
|
|
|
|
glutDisplayFunc(display);
|
|
|
|
glutMainLoop();
|
|
|
|
return 0;
|
|
}
|
|
|