2 |
#include <stdlib.h> |
#include <stdlib.h> |
3 |
#include <math.h> |
#include <math.h> |
4 |
#include <GL/glut.h> |
#include <GL/glut.h> |
5 |
|
#include <iostream> |
6 |
|
|
7 |
|
using std::cout; |
8 |
|
using std::cerr; |
9 |
|
|
10 |
//#define float double |
//#define float double |
11 |
|
|
12 |
#if 1 |
#if 1 |
13 |
/* 3D map */ |
/* 3D map */ |
14 |
/* Note: x, y, z \in [0, 2\pi] */ |
#define eps .015 |
15 |
#define eps .02 |
#define teps .5 |
16 |
|
|
17 |
//#define Du(i) 2E-5 / (eps*eps) |
#define Du_f(x,y,z) (2E-5 / (eps*eps)) |
|
#define Du(i) Du_arr[i] |
|
|
#define Du_init(x,y,z) exp(log(2E-5) + .7*triang2(z)) / (eps*eps) |
|
18 |
|
|
19 |
#define Dv(i) 1E-5 / (eps*eps) |
#define Dv_f(x,y,z) 1E-5 / (eps*eps) |
20 |
|
//#define Dv_f(x,y,z) (( .75E-5 + .5E-5 * z ) / (eps*eps)) |
21 |
|
|
22 |
//#define F(i) 0.04 |
#define F_f(x,y,z) 0.015 |
23 |
#define F(i) F_arr[i] |
//#define F_f(x,y,z) (.08 * y) |
24 |
#define F_init(x,y,z) .04 - .04 * triang2(y) |
//#define F_f(x,y,z) (.04 + .01 * y) |
25 |
|
|
26 |
//#define k(i) 0.06 |
#define k_f(x,y,z) 0.04 |
27 |
#define k(i) k_arr[i] |
//#define k_f(x,y,z) (.03 + .04 * x) |
28 |
#define k_init(x,y,z) .05 - .02 * triang2(x) |
//#define k_f(x,y,z) (.06 + .0075 * x) |
29 |
|
|
30 |
const int alln = 128; |
#define SCALE 81 |
31 |
|
const int alln = 81; |
32 |
const int n0 = alln+2; |
const int n0 = alln+2; |
33 |
const int n1 = alln+2; |
const int n1 = alln+2; |
34 |
const int n2 = alln+2; |
const int n2 = alln+2; |
39 |
|
|
40 |
#define eps .01 |
#define eps .01 |
41 |
|
|
42 |
#define Du(i) 2E-5 / (eps*eps) |
#define Du_f(x,y,z) (2E-5 / (eps*eps)) |
43 |
#define Dv(i) 1E-5 / (eps*eps) |
#define Dv_f(x,y,z) (1E-5 / (eps*eps)) |
|
|
|
|
#define F(i) F_arr[i] |
|
|
#define F_init(x,y,z) .04 - .04 * triang2(y) |
|
44 |
|
|
45 |
#define k(i) k_arr[i] |
#define F_f(x,y,z) (.08 * y) |
46 |
#define k_init(x,y,z) .05 - .02 * triang2(x) |
#define k_f(x,y,z) (.03 + .04 * x) |
47 |
|
|
48 |
|
#define SCALE 512 |
49 |
const int n0 = 512+2; |
const int n0 = 512+2; |
50 |
const int n1 = 512+2; |
const int n1 = 512+2; |
51 |
const int n2 = 1+2; |
const int n2 = 1+2; |
68 |
|
|
69 |
bool repaintRightaway = false; |
bool repaintRightaway = false; |
70 |
|
|
|
// priod 2 pi |
|
|
// /\ |
|
|
// / \ |
|
|
// \ / |
|
|
// \/ |
|
|
float triang(float x) { |
|
|
return asin(sin(x))/M_PI_2; |
|
|
} |
|
|
|
|
|
// \ |
|
|
// \ |
|
|
// \ |
|
|
// \ |
|
|
// |
|
|
float triang2(float x) { |
|
|
return 1-x/M_PI; |
|
|
} |
|
|
|
|
71 |
void init() { |
void init() { |
72 |
// Initialize arrays |
// Initialize arrays |
73 |
|
|
76 |
U2 = new float[n0 * n1 * n2]; |
U2 = new float[n0 * n1 * n2]; |
77 |
V2 = new float[n0 * n1 * n2]; |
V2 = new float[n0 * n1 * n2]; |
78 |
|
|
|
F_arr = new float[n0 * n1 * n2]; |
|
|
k_arr = new float[n0 * n1 * n2]; |
|
|
Du_arr = new float[n0 * n1 * n2]; |
|
|
Dv_arr = new float[n0 * n1 * n2]; |
|
|
|
|
79 |
for (int i = 0; i < n0 * n1 * n2; i++) { |
for (int i = 0; i < n0 * n1 * n2; i++) { |
80 |
|
|
81 |
float z = (float)(i / e2) / n2 * 2 * M_PI; |
float z = (float)(i / e2 - 1) / (n2 - 2) * 2 * M_PI; |
82 |
float y = (float)(i % e2 / e1) / n1 * 2 * M_PI; |
float y = (float)(i % e2 / e1 - 1) / (n1 - 2) * 2 * M_PI; |
83 |
float x = (float)(i % e2 % e1) / n0 * 2 * M_PI; |
float x = (float)(i % e2 % e1 - 1) / (n0 - 2) * 2 * M_PI; |
84 |
|
|
85 |
U[i] = 0.5+0.5*sin(x + y + 3*z); |
float xp = x + 2*sin(2*z); |
86 |
V[i] = 0.5+0.5*sin(2 * y - 2*z); |
float yp = y + 3*cos(x); |
87 |
|
float zp = z + sin(x) - sin(2*y); |
88 |
#ifdef k_init |
|
89 |
k_arr[i] = k_init(x,y,z); |
U[i] = 0.5+0.5*sin(xp + yp + 3*zp); |
90 |
#endif |
V[i] = 0.5+0.5*sin(2 * yp - 2*zp*zp); |
|
#ifdef F_init |
|
|
F_arr[i] = F_init(x,y,z); |
|
|
#endif |
|
|
#ifdef Du_init |
|
|
Du_arr[i] = Du_init(x,y,z); |
|
|
#endif |
|
|
#ifdef Dv_init |
|
|
Dv_arr[i] = Dv_init(x,y,z); |
|
|
#endif |
|
91 |
} |
} |
92 |
|
|
93 |
iter = 0; |
iter = 0; |
105 |
} |
} |
106 |
} |
} |
107 |
|
|
108 |
|
void copy2d2(float *data1, float *data2, |
109 |
|
int dst, int src, |
110 |
|
int e0, int n0, |
111 |
|
int e1, int n1) { |
112 |
|
copy2d(data1, dst, src, e0, n0, e1, n1); |
113 |
|
copy2d(data2, dst, src, e0, n0, e1, n1); |
114 |
|
} |
115 |
|
|
116 |
void evolve() { |
void evolve() { |
117 |
float *__restrict__ U = ::U; |
float *__restrict__ U = ::U; |
118 |
float *__restrict__ V = ::V; |
float *__restrict__ V = ::V; |
119 |
float *__restrict__ U2 = ::U2; |
float *__restrict__ U2 = ::U2; |
120 |
float *__restrict__ V2 = ::V2; |
float *__restrict__ V2 = ::V2; |
|
|
|
|
// Copy borders for proper tiling |
|
|
copy2d(U, 0, e2 * (n2-2), e0, n0, e1, n1); |
|
|
copy2d(U, e2 * (n2-1), e2, e0, n0, e1, n1); |
|
121 |
|
|
122 |
copy2d(V, 0, e2 * (n2-2), e0, n0, e1, n1); |
#if 1 /* Periodic boundary conditions */ |
123 |
copy2d(V, e2 * (n2-1), e2, e0, n0, e1, n1); |
copy2d2(U, V, 0, e2 * (n2-2), e0, n0, e1, n1); |
124 |
|
copy2d2(U, V, e2 * (n2-1), e2, e0, n0, e1, n1); |
125 |
|
|
126 |
copy2d(U, 0, e1 * (n1-2), e0, n0, e2, n2); |
copy2d2(U, V, 0, e1 * (n1-2), e0, n0, e2, n2); |
127 |
copy2d(U, e1 * (n1-1), e1, e0, n0, e2, n2); |
copy2d2(U, V, e1 * (n1-1), e1, e0, n0, e2, n2); |
128 |
|
|
129 |
copy2d(V, 0, e1 * (n1-2), e0, n0, e2, n2); |
copy2d2(U, V, 0, e0 * (n0-2), e1, n1, e2, n2); |
130 |
copy2d(V, e1 * (n1-1), e1, e0, n0, e2, n2); |
copy2d2(U, V, e0 * (n0-1), e0, e1, n1, e2, n2); |
131 |
|
#endif |
|
copy2d(U, 0, e0 * (n0-2), e1, n1, e2, n2); |
|
|
copy2d(U, e0 * (n0-1), e0, e1, n1, e2, n2); |
|
|
|
|
|
copy2d(V, 0, e0 * (n0-2), e1, n1, e2, n2); |
|
|
copy2d(V, e0 * (n0-1), e0, e1, n1, e2, n2); |
|
132 |
|
|
133 |
#ifndef LAPLACIAN_2D |
#ifndef LAPLACIAN_2D |
134 |
#define laplacian(a) (a[i-e0]+a[i+e0]+a[i-e1]+a[i+e1]+a[i-e2]+a[i+e2]-6*a[i]) |
#define laplacian(a) (a[i-e0]+a[i+e0]+a[i-e1]+a[i+e1]+a[i-e2]+a[i+e2]-6*a[i]) |
137 |
#endif |
#endif |
138 |
|
|
139 |
// Evaluate |
// Evaluate |
140 |
int n = n0 * n1 * n2; |
const float dx = 1.0 / (n0 - 2); |
141 |
for (int i = e2; i < n - e2; i++) { |
const float dy = 1.0 / (n1 - 2); |
142 |
U2[i] = U[i] + (Du(i) * laplacian(U) - U[i] * V[i]*V[i] + F(i) * (1 - U[i])); |
const float dz = 1.0 / (n2 - 2); |
143 |
V2[i] = V[i] + (Dv(i) * laplacian(V) + U[i] * V[i]*V[i] - (F(i)+k(i)) * V[i]); |
|
144 |
|
const float dk = dx * (k_f(1, 0, 0) - k_f(0, 0, 0)); |
145 |
|
const float dF = dx * (F_f(1, 0, 0) - F_f(0, 0, 0)); |
146 |
|
const float dDu = dx * (Du_f(1, 0, 0) - Du_f(0, 0, 0)); |
147 |
|
const float dDv = dx * (Dv_f(1, 0, 0) - Dv_f(0, 0, 0)); |
148 |
|
|
149 |
|
int i = e0 + e1 + e2; |
150 |
|
for (int i2 = 1; i2 < n2 - 1; i2++) { |
151 |
|
for (int i1 = 1; i1 < n1 - 1; i1++) { |
152 |
|
float z = dz * (i2 - 1); |
153 |
|
float y = dy * (i1 - 1); |
154 |
|
|
155 |
|
float k = k_f(0, y, z); |
156 |
|
float F = F_f(0, y, z); |
157 |
|
float Du = Du_f(0, y, z); |
158 |
|
float Dv = Dv_f(0, y, z); |
159 |
|
|
160 |
|
int i_end = i + n0 - 2; |
161 |
|
for (; i < i_end; i++) { |
162 |
|
U2[i] = U[i] + teps * (Du * laplacian(U) - U[i] * V[i]*V[i] + F * (1 - U[i])); |
163 |
|
V2[i] = V[i] + teps * (Dv * laplacian(V) + U[i] * V[i]*V[i] - (F+k) * V[i]); |
164 |
|
F += dF; |
165 |
|
k += dk; |
166 |
|
Du += dDu; |
167 |
|
Dv += dDv; |
168 |
|
} |
169 |
|
i += 2 * e0; |
170 |
|
} |
171 |
|
i += 2 * e1; |
172 |
} |
} |
173 |
|
|
174 |
// Swap |
// Swap |
175 |
::U = U2; ::U2 = U; |
::U = U2; ::U2 = U; |
176 |
::V = V2; ::V2 = V; |
::V = V2; ::V2 = V; |
178 |
iter++; |
iter++; |
179 |
} |
} |
180 |
|
|
181 |
|
void writeout() { |
182 |
|
const char *filenew = "map3d.dat.new"; |
183 |
|
char *file = "map3d.dat"; |
184 |
|
|
185 |
|
|
186 |
|
FILE *f = fopen(filenew, "w"); |
187 |
|
if (!f) { |
188 |
|
perror(file); |
189 |
|
exit(1); |
190 |
|
} |
191 |
|
|
192 |
|
fprintf(f, "%d %d %d 2\n", n2-2, n1-2, n0-2); |
193 |
|
|
194 |
|
for (int i2 = 1; i2 < n2 - 1; i2++) |
195 |
|
for (int i1 = 1; i1 < n1 - 1; i1++) |
196 |
|
for (int i0 = 1; i0 < n0 - 1; i0++) { |
197 |
|
int i = e0 * i0 + e1 * i1 + e2 * i2; |
198 |
|
fwrite(&U[i], sizeof(U[i]), 1, f); |
199 |
|
fwrite(&V[i], sizeof(V[i]), 1, f); |
200 |
|
//fprintf(f, "%G %G\n", U[i], V[i]); |
201 |
|
} |
202 |
|
fclose(f); |
203 |
|
rename(filenew, file); |
204 |
|
|
205 |
|
} |
206 |
|
|
207 |
void stats() { |
void stats() { |
208 |
int n = n0 * n1 * n2; |
int n = n0 * n1 * n2; |
209 |
fprintf(stderr, "Iteration %d: ", iter); |
fprintf(stderr, "Iteration %d: ", iter); |
212 |
int num=0; |
int num=0; |
213 |
float umin=1E30, umax=-1E30, usum=0; |
float umin=1E30, umax=-1E30, usum=0; |
214 |
float vmin=1E30, vmax=-1E30, vsum=0; |
float vmin=1E30, vmax=-1E30, vsum=0; |
215 |
for (int i = e2; i < n - e2; i++) { |
for (int i2 = 1; i2 < n2 - 1; i2++) |
216 |
if (!finite(U2[i]) || !finite(V2[i])) { |
for (int i1 = 1; i1 < n1 - 1; i1++) |
217 |
fprintf(stderr, "Not finite %d\n", i); |
for (int i0 = 1; i0 < n0 - 1; i0++) { |
218 |
abort(); |
int i = e0 * i0 + e1 * i1 + e2 * i2; |
219 |
} |
if (!finite(U[i]) || !finite(V[i])) { |
220 |
|
fprintf(stderr, "Not finite %d\n", i); |
221 |
num++; |
abort(); |
222 |
if (U2[i] < umin) umin = U2[i]; |
} |
223 |
if (U2[i] > umax) umax = U2[i]; |
|
224 |
usum += U2[i]; |
num++; |
225 |
|
if (U[i] < umin) umin = U[i]; |
226 |
if (V2[i] < vmin) vmin = V2[i]; |
if (U[i] > umax) umax = U[i]; |
227 |
if (V2[i] > vmax) vmax = V2[i]; |
usum += U[i]; |
228 |
vsum += V2[i]; |
|
229 |
} |
if (V[i] < vmin) vmin = V[i]; |
230 |
|
if (V[i] > vmax) vmax = V[i]; |
231 |
|
vsum += V[i]; |
232 |
|
} |
233 |
|
|
234 |
fprintf(stderr, "\tu: %G %G %G\tv: %G %G %G\n", |
fprintf(stderr, "\tu: %G %G %G\tv: %G %G %G\n", |
235 |
umin, usum/num, umax, |
umin, usum/num, umax, |
236 |
vmin, vsum/num, vmax); |
vmin, vsum/num, vmax); |
253 |
c[3] = w; |
c[3] = w; |
254 |
} |
} |
255 |
|
|
256 |
|
float threshold = 51./128; |
257 |
|
|
258 |
void display() { |
void display() { |
259 |
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); |
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); |
290 |
glMultMatrixf(mat); |
glMultMatrixf(mat); |
291 |
|
|
292 |
|
|
293 |
glTranslatef(-1, -1, -1); |
glScalef(2./SCALE, 2./SCALE, 2./SCALE); |
294 |
glScalef(2./(n0 - 2), 2./(n1-2), 2./(n2-2)); |
glTranslatef(-n0/2, -n1/2, -n2/2); |
|
glTranslatef(-1, -1, -1); |
|
295 |
|
|
296 |
glPointSize(5); |
glPointSize(5); |
297 |
int i = 0; |
int i = 0; |
306 |
for (int i1 = 1; i1 < n1-1; i1++) { |
for (int i1 = 1; i1 < n1-1; i1++) { |
307 |
int i = i1 * e1 + i2 * e2 + e0; |
int i = i1 * e1 + i2 * e2 + e0; |
308 |
for (int i0 = 1; i0 < n0-1; i0++) { |
for (int i0 = 1; i0 < n0-1; i0++) { |
309 |
if(U[i] < 0.4 || i2 == through) { |
if(U[i] < threshold || i2 == through) { |
310 |
bool stripe = |
bool stripe = |
311 |
(i1 % 16 == 0) || |
(i1 % 16 == 0) || |
312 |
(i0 % 16 == 0) || |
(i0 % 16 == 0) || |
333 |
if(!repaintRightaway) { |
if(!repaintRightaway) { |
334 |
evolve(); i++; |
evolve(); i++; |
335 |
if (iter % 100 == 0) stats(); |
if (iter % 100 == 0) stats(); |
336 |
|
if (iter % 1000 == 0) writeout(); |
337 |
if (iter %20 == 0) repaintRightaway = true; |
if (iter %20 == 0) repaintRightaway = true; |
338 |
} else { |
} else { |
339 |
repaintRightaway = false; |
repaintRightaway = false; |
342 |
} |
} |
343 |
|
|
344 |
|
|
345 |
|
void keypress(unsigned char key, int x, int y) { |
346 |
|
cout << key << " pressed\n"; |
347 |
|
switch (key) { |
348 |
|
case '+': cout << (threshold += 1./128) << "\n"; break; |
349 |
|
case '-': cout << (threshold -= 1./128) << "\n"; break; |
350 |
|
} |
351 |
|
} |
352 |
|
|
353 |
|
|
354 |
bool buttons[3] = {0,0,0}; |
bool buttons[3] = {0,0,0}; |
355 |
float buttonx[3] = {0,0,0}; |
float buttonx[3] = {0,0,0}; |
406 |
|
|
407 |
glutDisplayFunc(display); |
glutDisplayFunc(display); |
408 |
glutIdleFunc(idle); |
glutIdleFunc(idle); |
409 |
|
glutKeyboardFunc(keypress); |
410 |
glutMouseFunc(mouse); |
glutMouseFunc(mouse); |
411 |
glutMotionFunc(motion); |
glutMotionFunc(motion); |
412 |
glutMainLoop(); |
glutMainLoop(); |