Show rd1.cxx syntax highlighted
/*
rd1.cxx
*
* Copyright (c) 2003, Janne Kujala and Tuomas J. Lukka
*
* This file is part of Gzz.
*
* Gzz is free software; you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* Gzz is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
* Public License for more details.
*
* You should have received a copy of the GNU Lesser General
* Public License along with Gzz; if not, write to the Free
* Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*
*
*/
/*
* Written by Janne Kujala and Tuomas J. Lukka
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
#include <iostream>
#include <unistd.h>
using std::cout;
using std::cerr;
//#define float double
#if 1
/* 3D map */
#define eps .02
#define teps .5
#define Du_f(x,y,z) (2E-5 / (eps*eps))
#define Dv_f(x,y,z) 1E-5 / (eps*eps)
//#define Dv_f(x,y,z) (( .75E-5 + .5E-5 * z ) / (eps*eps))
#define F_f(x,y,z) 0.05
//#define F_f(x,y,z) (.08 * y)
//#define F_f(x,y,z) (.04 + .01 * y)
#define k_f(x,y,z) 0.065
//#define k_f(x,y,z) (.03 + .04 * x)
//#define k_f(x,y,z) (.06 + .0075 * x)
#define SCALE 64
const int alln = 64;
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;
float *F_arr;
float *k_arr;
float *Du_arr;
float *Dv_arr;
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 + 2*sin(2*z);
float yp = y + 3*cos(x);
float zp = z + sin(x) - sin(2*y);
U[i] = 0.5+0.5*sin(xp + yp + 3*zp);
V[i] = 0.5+0.5*sin(2 * yp - 2*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
// Evaluate
const float dx = 1.0 / (n0 - 2);
const float dy = 1.0 / (n1 - 2);
const float dz = 1.0 / (n2 - 2);
const float dk = dx * (k_f(1, 0, 0) - k_f(0, 0, 0));
const float dF = dx * (F_f(1, 0, 0) - F_f(0, 0, 0));
const float dDu = dx * (Du_f(1, 0, 0) - Du_f(0, 0, 0));
const float dDv = dx * (Dv_f(1, 0, 0) - Dv_f(0, 0, 0));
int i = e0 + e1 + e2;
for (int i2 = 1; i2 < n2 - 1; i2++) {
for (int i1 = 1; i1 < n1 - 1; i1++) {
float z = dz * (i2 - 1);
float y = dy * (i1 - 1);
float k = k_f(0, y, z);
float F = F_f(0, y, z);
float Du = Du_f(0, y, z);
float Dv = Dv_f(0, y, z);
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]);
F += dF;
k += dk;
Du += dDu;
Dv += dDv;
}
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_f(0,0,0), k_f(0,0,0), 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;
} 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[]) {
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();
}
See more files for this project here