5 |
|
|
6 |
//#define float double |
//#define float double |
7 |
|
|
8 |
float eps = .02; |
#if 1 |
9 |
float Du = 2E-5; |
/* 3D map */ |
10 |
float Dv = 1E-5; |
/* Note: x, y, z \in [0, 2\pi] */ |
11 |
|
#define eps .02 |
12 |
|
|
13 |
|
//#define Du(i) 2E-5 / (eps*eps) |
14 |
|
#define Du(i) Du_arr[i] |
15 |
|
#define Du_init(x,y,z) exp(log(2E-5) + .7*triang2(z)) / (eps*eps) |
16 |
|
|
17 |
|
#define Dv(i) 1E-5 / (eps*eps) |
18 |
|
|
19 |
//#define F(i) 0.04 |
//#define F(i) 0.04 |
|
#define k(i) 0.06 |
|
20 |
#define F(i) F_arr[i] |
#define F(i) F_arr[i] |
21 |
//#define k(i) k_arr[i] |
#define F_init(x,y,z) .04 - .04 * triang2(y) |
22 |
|
|
23 |
|
//#define k(i) 0.06 |
24 |
|
#define k(i) k_arr[i] |
25 |
|
#define k_init(x,y,z) .05 - .02 * triang2(x) |
26 |
|
|
27 |
const int alln = 128; |
const int alln = 128; |
28 |
const int n0 = alln+2; |
const int n0 = alln+2; |
29 |
const int n1 = alln+2; |
const int n1 = alln+2; |
30 |
const int n2 = alln+2; |
const int n2 = alln+2; |
31 |
|
|
32 |
|
#else |
33 |
|
/* Xmorphia map parameters */ |
34 |
|
#define LAPLACIAN_2D |
35 |
|
|
36 |
|
#define eps .01 |
37 |
|
|
38 |
|
#define Du(i) 2E-5 / (eps*eps) |
39 |
|
#define Dv(i) 1E-5 / (eps*eps) |
40 |
|
|
41 |
|
#define F(i) F_arr[i] |
42 |
|
#define F_init(x,y,z) .04 - .04 * triang2(y) |
43 |
|
|
44 |
|
#define k(i) k_arr[i] |
45 |
|
#define k_init(x,y,z) .05 - .02 * triang2(x) |
46 |
|
|
47 |
|
const int n0 = 512+2; |
48 |
|
const int n1 = 512+2; |
49 |
|
const int n2 = 1+2; |
50 |
|
#endif |
51 |
|
|
52 |
const int e0 = 1; |
const int e0 = 1; |
53 |
const int e1 = n0; |
const int e1 = n0; |
54 |
const int e2 = n1 * n0; |
const int e2 = n1 * n0; |
59 |
float *V2; |
float *V2; |
60 |
float *F_arr; |
float *F_arr; |
61 |
float *k_arr; |
float *k_arr; |
62 |
|
float *Du_arr; |
63 |
|
float *Dv_arr; |
64 |
|
|
65 |
int iter; |
int iter; |
66 |
|
|
67 |
bool repaintRightaway = false; |
bool repaintRightaway = false; |
68 |
|
|
69 |
|
// priod 2 pi |
70 |
|
// /\ |
71 |
|
// / \ |
72 |
|
// \ / |
73 |
|
// \/ |
74 |
|
float triang(float x) { |
75 |
|
return asin(sin(x))/M_PI_2; |
76 |
|
} |
77 |
|
|
78 |
|
// \ |
79 |
|
// \ |
80 |
|
// \ |
81 |
|
// \ |
82 |
|
// |
83 |
|
float triang2(float x) { |
84 |
|
return 1-x/M_PI; |
85 |
|
} |
86 |
|
|
87 |
void init() { |
void init() { |
88 |
// Initialize arrays |
// Initialize arrays |
89 |
|
|
94 |
|
|
95 |
F_arr = new float[n0 * n1 * n2]; |
F_arr = new float[n0 * n1 * n2]; |
96 |
k_arr = new float[n0 * n1 * n2]; |
k_arr = new float[n0 * n1 * n2]; |
97 |
|
Du_arr = new float[n0 * n1 * n2]; |
98 |
|
Dv_arr = new float[n0 * n1 * n2]; |
99 |
|
|
100 |
for (int i = 0; i < n0 * n1 * n2; i++) { |
for (int i = 0; i < n0 * n1 * n2; i++) { |
101 |
|
|
106 |
U[i] = 0.5+0.5*sin(x + y + 3*z); |
U[i] = 0.5+0.5*sin(x + y + 3*z); |
107 |
V[i] = 0.5+0.5*sin(2 * y - 2*z); |
V[i] = 0.5+0.5*sin(2 * y - 2*z); |
108 |
|
|
109 |
F_arr[i] = +.05 + .05 * sin(x); |
#ifdef k_init |
110 |
k_arr[i] = +.04 + .04 * sin(y); |
k_arr[i] = k_init(x,y,z); |
111 |
|
#endif |
112 |
|
#ifdef F_init |
113 |
|
F_arr[i] = F_init(x,y,z); |
114 |
|
#endif |
115 |
|
#ifdef Du_init |
116 |
|
Du_arr[i] = Du_init(x,y,z); |
117 |
|
#endif |
118 |
|
#ifdef Dv_init |
119 |
|
Dv_arr[i] = Dv_init(x,y,z); |
120 |
|
#endif |
121 |
} |
} |
122 |
|
|
123 |
iter = 0; |
iter = 0; |
160 |
copy2d(V, 0, e0 * (n0-2), e1, n1, e2, n2); |
copy2d(V, 0, e0 * (n0-2), e1, n1, e2, n2); |
161 |
copy2d(V, e0 * (n0-1), e0, e1, n1, e2, n2); |
copy2d(V, e0 * (n0-1), e0, e1, n1, e2, n2); |
162 |
|
|
163 |
|
#ifndef LAPLACIAN_2D |
164 |
#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]) |
165 |
|
#else |
166 |
|
#define laplacian(a) (a[i-e0]+a[i+e0]+a[i-e1]+a[i+e1]-4*a[i]) |
167 |
|
#endif |
168 |
|
|
169 |
// Evaluate |
// Evaluate |
|
float Du = ::Du / (eps * eps); |
|
|
float Dv = ::Dv / (eps * eps); |
|
|
|
|
170 |
int n = n0 * n1 * n2; |
int n = n0 * n1 * n2; |
171 |
for (int i = e2; i < n - e2; i++) { |
for (int i = e2; i < n - e2; i++) { |
172 |
U2[i] = U[i] + (Du * laplacian(U) - U[i] * V[i]*V[i] + F(i) * (1 - U[i])); |
U2[i] = U[i] + (Du(i) * laplacian(U) - U[i] * V[i]*V[i] + F(i) * (1 - U[i])); |
173 |
V2[i] = V[i] + (Dv * laplacian(V) + U[i] * V[i]*V[i] - (F(i)+k(i)) * V[i]); |
V2[i] = V[i] + (Dv(i) * laplacian(V) + U[i] * V[i]*V[i] - (F(i)+k(i)) * V[i]); |
174 |
} |
} |
175 |
|
|
176 |
// Swap |
// Swap |
210 |
} |
} |
211 |
|
|
212 |
|
|
213 |
float scale = .8f; |
float scale = 1.f; |
214 |
float rotate[4] = {1,0,0,0}; |
float rotate[4] = {1,0,0,0}; |
215 |
|
|
216 |
|
|
263 |
|
|
264 |
glTranslatef(-1, -1, -1); |
glTranslatef(-1, -1, -1); |
265 |
glScalef(2./(n0 - 2), 2./(n1-2), 2./(n2-2)); |
glScalef(2./(n0 - 2), 2./(n1-2), 2./(n2-2)); |
266 |
|
glTranslatef(-1, -1, -1); |
267 |
|
|
268 |
glPointSize(5); |
glPointSize(5); |
269 |
int i = 0; |
int i = 0; |
272 |
|
|
273 |
static int through; |
static int through; |
274 |
through = through%(n2-2)+1; |
through = through%(n2-2)+1; |
275 |
for (int i2 = 0; i2 < n2; i2++) |
for (int i2 = 1; i2 < n2-1; i2++) |
276 |
{ |
{ |
277 |
int i = i2 * e2; |
glBegin(GL_POINTS); |
278 |
glBegin(GL_POINTS); |
for (int i1 = 1; i1 < n1-1; i1++) { |
279 |
for (int i1 = 0; i1 < n1; i1++) |
int i = i1 * e1 + i2 * e2 + e0; |
280 |
for (int i0 = 0; i0 < n0; i0++) { |
for (int i0 = 1; i0 < n0-1; i0++) { |
281 |
if(U[i] < 0.4 || i2 == through) { |
if(U[i] < 0.4 || i2 == through) { |
282 |
bool stripe = |
bool stripe = |
283 |
(i1 % 16 == 0) || |
(i1 % 16 == 0) || |
288 |
} |
} |
289 |
i++; |
i++; |
290 |
} |
} |
291 |
glEnd(); |
} |
292 |
|
glEnd(); |
293 |
} |
} |
294 |
|
|
295 |
glPopMatrix(); |
glPopMatrix(); |