1 /*
2 * $Id: demo2.i,v 1.1 2005-09-18 22:05:55 dhmunro Exp $
3 * Mesh plotting demo
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
demo2(which,time_limit)11 func demo2(which, time_limit)
12 /* DOCUMENT demo2
13 Exhibit quadrilateral mesh plots in 3 movies of a drumhead.
14 The drumhead is initially stationary, but has a bump near one
15 edge. Yorick is solving a 2D wave equation to compute the
16 evolution of this bump.
17
18 The first movie is a filled mesh plot with color "proportional"
19 to the height of the surface of the drum. A few well chosen
20 contour levels (here 3) add a lot to a filled mesh plot.
21
22 The second movie is a "3D" perspective plot of the height of the
23 drumhead. In this movie, the mesh lines are drawn, which is
24 slightly confusing since the cells are not all the same shape.
25
26 The second movie is a "3D" shaded plot of the height of the
27 drumhead. Yorick computes surface shading based on the angle
28 of each cell from a light source.
29
30 As you watch this, you might reflect on the two dimensionality
31 of your retina. What Yorick lacks by way of 3D graphics is
32 really just fancy hidden surface algorithms; the simple
33 painter's algorithm used here and in plwf.i is easy to
34 implement.
35
36 There are two optional arguments to demo2: the first is the
37 number fo the movie (1, 2, or 3) you want to watch; the second
38 is a time limit on the duration of each movie in seconds (default
39 is 60 seconds each).
40 */
41 {
42 require, "movie.i";
43
44 /* generate a 30-by-30 cell mesh on the [-1,1] square */
45 x= span(-1, 1, 31)(,-:1:31);
46 y= transpose(x);
47 /* map the square mesh into a mesh on the unit circle
48 this mesh has more nearly equal area cells than a polar
49 coordinate circle */
50 scale= max(abs(y),abs(x))/(abs(y,x)+1.e-30);
51 /* note that abs(y,x)=sqrt(x^2+y^2) */
52 x*= scale;
53 y*= scale;
54
55 f= f0= exp(-8.*abs(y+.67,x+.25)^2)*(1.-abs(y,x)^2);
56 fdot= 0.0*f(2:-1,2:-1);
57
58 af= abs(f(2:-1,2:-1));
59 lf= laplacian(f, y,x);
60 xdz= x(dif,zcen);
61 xzd= x(zcen,dif);
62 ydz= y(dif,zcen);
63 yzd= y(zcen,dif);
64 dt= 0.375*sqrt(min(abs(xdz*yzd - xzd*ydz)));
65
66 window, 0, wait=1, style="nobox.gs";
67 palette, "heat.gp";
68 limits, -1, 1, -1, 1;
69
70 /* roll the filled mesh movie */
71 if (which && which!=1) goto persp;
72 fc= f(zcen,zcen);
73 cmin= cmax= max(abs(fc));
74 cmin= -cmin;
75 level= cmax/4.;
76 movie, display_plf, time_limit;
77 write,format="%f frames of filled mesh drumhead completed in %f sec\n",
78 movie_timing(4), movie_timing(3);
79 write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
80 movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4),
81 movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
82
83 /* roll the perspective movie */
84 persp: if (which && which!=2) goto shade;
85 f= f0;
86 limits, -1, 1, -1, 1;
87 movie, display_plm, time_limit;
88 write,format="%f frames of wireframe surface drumhead completed in %f sec\n",
89 movie_timing(4), movie_timing(3);
90 write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
91 movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4),
92 movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
93
94 /* roll the shaded movie */
95 shade: if (which && which!=3) return;
96 f= f0;
97 limits, -1, 1, -1, 1;
98 movie, display_pl3, time_limit;
99 write,format="%f frames of filled surface drumhead completed in %f sec\n",
100 movie_timing(4), movie_timing(3);
101 write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
102 movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4),
103 movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
104
105 fma;
106 limits;
107 }
108
display_plf(i)109 func display_plf(i)
110 {
111 /* display first */
112 fc= f(zcen,zcen);
113 cmin= cmax= max(abs(fc));
114 cmin= -cmin;
115 plf, fc, -y, -x, cmin=cmin, cmax=cmax;
116 /* the 0 contour level is too noisy without some smoothing... */
117 fs= f(zcen,zcen)(pcen,pcen);
118 plc, fs, -y, -x, levs=0., marks=0, color="green", type="solid";
119 plc, f, -y, -x, levs=level, marks=0, color="black", type="dash";
120 plc, f, -y, -x, levs=-level, marks=0, color="green", type="dash";
121
122 /* then take a step forward in time */
123 lf= laplacian(f, y,x);
124 af= abs(f(2:-1,2:-1));
125 fdot+= lf*dt;
126 f(2:-1,2:-1)+= fdot*dt;
127
128 return i<200;
129 }
130
display_plm(i)131 func display_plm(i)
132 {
133 /* display first */
134 pl3d,0, f, y, x;
135
136 /* then take a step forward in time */
137 lf= laplacian(f, y,x);
138 af= abs(f(2:-1,2:-1));
139 fdot+= lf*dt;
140 f(2:-1,2:-1)+= fdot*dt;
141
142 return i<200;
143 }
144
display_pl3(i)145 func display_pl3(i)
146 {
147 /* display first */
148 pl3d,1, f, y, x;
149
150 /* then take a step forward in time */
151 lf= laplacian(f, y,x);
152 af= abs(f(2:-1,2:-1));
153 fdot+= lf*dt;
154 f(2:-1,2:-1)+= fdot*dt;
155
156 return i<200;
157 }
158
laplacian(f,y,x)159 func laplacian(f, y,x)
160 {
161 /* There are many ways to form the Laplacian as a finite difference.
162 This one is nice in Yorick. */
163 /* Start with the two median vectors across each zone. */
164 fdz= f(dif,zcen);
165 fzd= f(zcen,dif);
166 xdz= x(dif,zcen);
167 xzd= x(zcen,dif);
168 ydz= y(dif,zcen);
169 yzd= y(zcen,dif);
170
171 /* Estimate the gradient at the center of each cell. */
172 area= xdz*yzd - xzd*ydz;
173 gradfx= (fdz*yzd - fzd*ydz)/area;
174 gradfy= (xdz*fzd - xzd*fdz)/area;
175
176 /* Now consider the mesh formed by the center points of the original. */
177 x= x(zcen,zcen);
178 y= y(zcen,zcen);
179 xdz= x(dif,);
180 xzd= x(,dif);
181 ydz= y(dif,);
182 yzd= y(,dif);
183 area= xdz(,zcen)*yzd(zcen,) - xzd(zcen,)*ydz(,zcen);
184
185 return ((xdz*gradfy(zcen,)-ydz*gradfx(zcen,))(,dif) +
186 (yzd*gradfx(,zcen)-xzd*gradfy(,zcen))(dif,)) / area;
187 }
188
pl3d(shading,z,y,x)189 func pl3d(shading, z, y, x)
190 {
191 /* rotate so that (zp,yp) are screen (y,x) */
192 /* These orientations are cunningly chosen so that the painter's
193 algorithm correctly draws hidden surfaces first -- see help, plf
194 for a description of the order cells are drawn by plf. */
195 theta= 30. * pi/180.; /* angle of viewer above drumhead */
196 phi= 120. * pi/180;
197
198 ct= cos(phi);
199 st= sin(phi);
200 yp= y*ct - x*st;
201 xp= x*ct + y*st;
202
203 ct= cos(theta);
204 st= sin(theta);
205 zp= z*ct - xp*st;
206 xp= xp*ct + z*st;
207
208 if (!shading) {
209 color= [];
210 edges= 1;
211 } else {
212 /* compute the two median vectors for each cell */
213 m0x= xp(dif,zcen);
214 m0y= yp(dif,zcen);
215 m0z= zp(dif,zcen);
216 m1x= xp(zcen,dif);
217 m1y= yp(zcen,dif);
218 m1z= zp(zcen,dif);
219 /* define the normal vector to be their cross product */
220 nx= m0y*m1z - m0z*m1y;
221 ny= m0z*m1x - m0x*m1z;
222 nz= m0x*m1y - m0y*m1x;
223 n= abs(nx, ny, nz);
224 nx/= n;
225 ny/= n;
226 nz/= n;
227 color= bytscl(nx, cmin=0.0, cmax=1.0);
228 edges= 0;
229 }
230
231 plf, color, zp, yp, edges=edges;
232 }
233