/[gzz]/gzz/lava/gfx/rd1.cxx
ViewVC logotype

Diff of /gzz/lava/gfx/rd1.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by jvk, Tue Aug 27 13:40:17 2002 UTC revision 1.5 by jvk, Fri Aug 30 10:32:35 2002 UTC
# Line 2  Line 2 
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;
# Line 35  const int n2 = alln+2; Line 39  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;
# Line 66  int iter; Line 68  int iter;
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    
# Line 92  void init() { Line 76  void init() {
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;
# Line 135  void copy2d(float *data, Line 105  void copy2d(float *data,
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])
# Line 167  void evolve() { Line 137  void evolve() {
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;
# Line 180  void evolve() { Line 178  void evolve() {
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);
# Line 188  void stats() { Line 212  void stats() {
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);
# Line 225  void quatmul(float a[], float b[], float Line 253  void quatmul(float a[], float b[], float
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);
# Line 261  void display() { Line 290  void display() {
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;
# Line 278  void display() { Line 306  void display() {
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) ||
# Line 305  void idle() { Line 333  void idle() {
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;
# Line 313  void idle() { Line 342  void idle() {
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};
# Line 369  int main(int argc, char *argv[]) { Line 406  int main(int argc, char *argv[]) {
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();

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26