#include #include #include #include #include #include using std::cout; using std::cerr; //#define float double #if 1 float eps = .015; float teps = .5; float Du = 2E-5; float Dv = 1E-5; float F = .05; float k = .065; int maxiter = 10000; #define SCALE 128 const int alln = 128; const int n0 = alln+2; const int n1 = alln+2; const int n2 = alln+2; #else /* Xmorphia map parameters */ #define LAPLACIAN_2D #define eps .01 #define Du_f(x,y,z) (2E-5 / (eps*eps)) #define Dv_f(x,y,z) (1E-5 / (eps*eps)) #define F_f(x,y,z) (.08 * y) #define k_f(x,y,z) (.03 + .04 * x) #define SCALE 512 const int n0 = 512+2; const int n1 = 512+2; const int n2 = 1+2; #endif const int e0 = 1; const int e1 = n0; const int e2 = n1 * n0; float *U; float *V; float *U2; float *V2; int iter; bool repaintRightaway = false; void init() { // Initialize arrays U = new float[n0 * n1 * n2]; V = new float[n0 * n1 * n2]; U2 = new float[n0 * n1 * n2]; V2 = new float[n0 * n1 * n2]; for (int i = 0; i < n0 * n1 * n2; i++) { float z = (float)(i / e2 - 1) / (n2 - 2) * 2 * M_PI; float y = (float)(i % e2 / e1 - 1) / (n1 - 2) * 2 * M_PI; float x = (float)(i % e2 % e1 - 1) / (n0 - 2) * 2 * M_PI; float xp = x + 1.5*sin(2*z+3*y); float yp = y + 1.8*cos(x+z); float zp = z + 2.5*sin(3*x-2*y); U[i] = 0.5+(0.25*sin(xp + yp + zp)+0.25*cos(2*x+z)*sin(yp)); V[i] = 0.5+(0.25*sin(2*yp - 2*zp)+0.25*cos(y+2*z)*sin(zp)); } iter = 0; } void copy2d(float *data, int dst, int src, int e0, int n0, int e1, int n1) { for (int i1 = 0; i1 < n1; i1++) { for (int i0 = 0; i0 < n0; i0++) { int i = e0 * i0 + e1 * i1; data[i + dst] = data[i + src]; } } } void copy2d2(float *data1, float *data2, int dst, int src, int e0, int n0, int e1, int n1) { copy2d(data1, dst, src, e0, n0, e1, n1); copy2d(data2, dst, src, e0, n0, e1, n1); } void evolve() { float *__restrict__ U = ::U; float *__restrict__ V = ::V; float *__restrict__ U2 = ::U2; float *__restrict__ V2 = ::V2; #if 1 /* Periodic boundary conditions */ copy2d2(U, V, 0, e2 * (n2-2), e0, n0, e1, n1); copy2d2(U, V, e2 * (n2-1), e2, e0, n0, e1, n1); copy2d2(U, V, 0, e1 * (n1-2), e0, n0, e2, n2); copy2d2(U, V, e1 * (n1-1), e1, e0, n0, e2, n2); copy2d2(U, V, 0, e0 * (n0-2), e1, n1, e2, n2); copy2d2(U, V, e0 * (n0-1), e0, e1, n1, e2, n2); #endif #ifndef LAPLACIAN_2D #define laplacian6(a) (a[i-e0]+a[i+e0]+a[i-e1]+a[i+e1]+a[i-e2]+a[i+e2]-6*a[i]) #define laplacian8(a) (a[i-e0-e1-e2]+\ a[i-e0-e1+e2]+\ a[i-e0+e1-e2]+\ a[i-e0+e1+e2]+\ a[i+e0-e1-e2]+\ a[i+e0-e1+e2]+\ a[i+e0+e1-e2]+\ a[i+e0+e1+e2]-8*a[i])*(1./3) #define laplacian(a) (.5*laplacian6(a)+.5*laplacian8(a)) #else #define laplacian(a) (a[i-e0]+a[i+e0]+a[i-e1]+a[i+e1]-4*a[i]) #endif float Du = ::Du / (eps * eps); float Dv = ::Dv / (eps * eps); // Evaluate int i = e0 + e1 + e2; for (int i2 = 1; i2 < n2 - 1; i2++) { for (int i1 = 1; i1 < n1 - 1; i1++) { int i_end = i + n0 - 2; for (; i < i_end; i++) { U2[i] = U[i] + teps * (Du * laplacian(U) - U[i] * V[i]*V[i] + F * (1 - U[i])); V2[i] = V[i] + teps * (Dv * laplacian(V) + U[i] * V[i]*V[i] - (F+k) * V[i]); } i += 2 * e0; } i += 2 * e1; } // Swap ::U = U2; ::U2 = U; ::V = V2; ::V2 = V; iter++; } void writeout() { char file[1000]; sprintf(file, "map3d_F%G_k%G_eps%G_teps%G_iter%d.dat", F, k, eps, teps, iter); FILE *f = fopen(file, "w"); if (!f) { perror(file); exit(1); } fprintf(f, "%d %d %d 2\n", n2-2, n1-2, n0-2); for (int i2 = 1; i2 < n2 - 1; i2++) for (int i1 = 1; i1 < n1 - 1; i1++) for (int i0 = 1; i0 < n0 - 1; i0++) { int i = e0 * i0 + e1 * i1 + e2 * i2; fwrite(&U[i], sizeof(U[i]), 1, f); fwrite(&V[i], sizeof(V[i]), 1, f); //fprintf(f, "%G %G\n", U[i], V[i]); } fclose(f); unlink("map3d.dat"); link(file, "map3d.dat"); } void stats() { int n = n0 * n1 * n2; fprintf(stderr, "Iteration %d: ", iter); fflush(stderr); int num=0; float umin=1E30, umax=-1E30, usum=0; float vmin=1E30, vmax=-1E30, vsum=0; for (int i2 = 1; i2 < n2 - 1; i2++) for (int i1 = 1; i1 < n1 - 1; i1++) for (int i0 = 1; i0 < n0 - 1; i0++) { int i = e0 * i0 + e1 * i1 + e2 * i2; if (!finite(U[i]) || !finite(V[i])) { fprintf(stderr, "Not finite %d\n", i); abort(); } num++; if (U[i] < umin) umin = U[i]; if (U[i] > umax) umax = U[i]; usum += U[i]; if (V[i] < vmin) vmin = V[i]; if (V[i] > vmax) vmax = V[i]; vsum += V[i]; } fprintf(stderr, "\tu: %G %G %G\tv: %G %G %G\n", umin, usum/num, umax, vmin, vsum/num, vmax); fflush(stderr); } float scale = 1.f; float rotate[4] = {1,0,0,0}; void quatmul(float a[], float b[], float c[]) { float x = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]; float y = a[0]*b[1] + a[1]*b[0] + a[2]*b[3] - a[3]*b[2]; float z = a[0]*b[2] - a[1]*b[3] + a[2]*b[0] + a[3]*b[1]; float w = a[0]*b[3] + a[1]*b[2] - a[2]*b[1] + a[3]*b[0]; c[0] = x; c[1] = y; c[2] = z; c[3] = w; } float threshold = 51./128; void display() { glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); glPushMatrix(); glScalef(scale, scale, scale); // Compute rotation matrix from the quaternion float aa = rotate[0]*rotate[0]; float ab = rotate[0]*rotate[1]; float ac = rotate[0]*rotate[2]; float ad = rotate[0]*rotate[3]; float bb = rotate[1]*rotate[1]; float bc = rotate[1]*rotate[2]; float bd = rotate[1]*rotate[3]; float cc = rotate[2]*rotate[2]; float cd = rotate[2]*rotate[3]; float dd = rotate[3]*rotate[3]; float m = 1. / (aa + bb + cc + dd); GLfloat mat[16] = { m*(aa+bb-cc-dd), m*2*(ad+bc), m*2*(-ac+bd), 0, m*2*(-ad+bc), m*(aa-bb+cc-dd), m*2*(ab+cd), 0, m*2*(ac+bd), m*2*(-ab+cd), m*(aa-bb-cc+dd), 0, 0, 0, 0, 1 }; glMultMatrixf(mat); glScalef(2./SCALE, 2./SCALE, 2./SCALE); glTranslatef(-n0/2, -n1/2, -n2/2); glPointSize(5); int i = 0; glEnable(GL_DEPTH_TEST); static int through; through = through%(n2-2)+1; for (int i2 = 1; i2 < n2-1; i2++) { glBegin(GL_POINTS); for (int i1 = 1; i1 < n1-1; i1++) { int i = i1 * e1 + i2 * e2 + e0; for (int i0 = 1; i0 < n0-1; i0++) { if(U[i] < threshold || i2 == through) { bool stripe = (i1 % 16 == 0) || (i0 % 16 == 0) || (i2 % 16 == 0) ; glColor3f(U[i], V[i], stripe); glVertex3f(i0, i1, i2); } i++; } } glEnd(); } glPopMatrix(); glutSwapBuffers(); glFlush(); glFinish(); } void idle() { static int i; if(!repaintRightaway) { evolve(); i++; if (iter % 100 == 0) stats(); if (iter % 1000 == 0) writeout(); if (iter %20 == 0) repaintRightaway = true; if (iter >= maxiter) exit(0); } else { repaintRightaway = false; glutPostRedisplay(); } } void keypress(unsigned char key, int x, int y) { cout << key << " pressed\n"; switch (key) { case '+': cout << (threshold += 1./128) << "\n"; break; case '-': cout << (threshold -= 1./128) << "\n"; break; } } bool buttons[3] = {0,0,0}; float buttonx[3] = {0,0,0}; float buttony[3] = {0,0,0}; void mouse(int button, int state, int x, int y) { repaintRightaway = true; if (button > 2) return; buttonx[button] = x; buttony[button] = y; buttons[button] = (state == GLUT_DOWN); } void motion(int x, int y) { if (buttons[0]) { float a = .01*(buttonx[0] - x); float b = .01*(buttony[0] - y); buttonx[0] = x; buttony[0] = y; float q0[4] = { cos(.5*a), 0, sin(.5*a), 0 }; float q1[4] = { cos(.5*b), sin(.5*b), 0, 0 }; quatmul(q0, rotate, rotate); quatmul(q1, rotate, rotate); } if (buttons[2]) { float cx = .5 * glutGet(GLUT_WINDOW_WIDTH); float cy = .5 * glutGet(GLUT_WINDOW_HEIGHT); float a = atan2(x-cx, y-cy) - atan2(buttonx[2]-cx, buttony[2]-cy); buttonx[2] = x; buttony[2] = y; float q0[4] = { cos(.5*a), 0, 0, sin(.5*a) }; quatmul(q0, rotate, rotate); } if (buttons[1]) { scale *= exp(.01 * (y - buttony[1])); buttonx[1] = x; buttony[1] = y; } glutPostRedisplay(); } int main(int argc, char *argv[]) { int i; while ((i = getopt(argc, argv, "F:k:e:t:i:")) != EOF) switch (i) { case '?': case 'h': break; case 'F': F = atof(optarg); break; case 'k': k = atof(optarg); break; case 'e': eps = atof(optarg); break; case 't': teps = atof(optarg); break; case 'i': maxiter = atoi(optarg); break; } glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH); glutCreateWindow("rd1"); glClearColor(0, 0, 0, 0); init(); glutDisplayFunc(display); glutIdleFunc(idle); glutKeyboardFunc(keypress); glutMouseFunc(mouse); glutMotionFunc(motion); glutMainLoop(); }