source: inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py @ 1295

Last change on this file since 1295 was 1295, checked in by steve, 20 years ago

Playing with coloured surfaces

  • Property svn:executable set to *
File size: 15.4 KB
Line 
1from visual import *
2
3elevation_color = (0.3,0.4,0.5)
4stage_color = (0.1,0.4,0.99)
5
6class Surface:
7    def __init__(self,domain,scale_z=1.0,range_z=None):
8        """Create visualisation of domain
9        """
10
11        self.frame  = frame()
12
13        self.domain = domain
14        self.scale_z  = scale_z
15        self.vertices = domain.vertex_coordinates
16
17        # models for each quantity
18        self.z_models = {}
19        keys = self.domain.quantities.keys()
20        #print keys
21        for key in keys:
22            self.z_models[key] = faces(frame= self.frame)
23
24
25
26        self.max_x = max(max(self.vertices[:,0]),max(self.vertices[:,2]),max(self.vertices[:,4]))
27        self.min_x = min(min(self.vertices[:,0]),min(self.vertices[:,2]),min(self.vertices[:,4]))
28        self.max_y = max(max(self.vertices[:,1]),max(self.vertices[:,3]),max(self.vertices[:,5]))
29        self.min_y = min(min(self.vertices[:,1]),min(self.vertices[:,3]),min(self.vertices[:,5]))
30        self.range_x = self.max_x - self.min_x
31        self.range_y = self.max_y - self.min_y
32        self.range_xy = max(self.range_x, self.range_y)
33
34#        print 'min_x=',self.min_x
35#        print 'max_x=',self.max_x
36#        print 'range_x=',self.range_x
37#        print 'min_y=',self.min_y
38#        print 'max_y=',self.max_y
39#        print 'range_y',self.range_y
40
41        self.range_z = 1.0
42        self.max_z = self.range_z/2.0
43        self.min_z = -self.range_z/2.0
44        self.range_z = max(self.max_z - self.min_z, 1.0e-10)
45
46
47        #print self.max_x, self.min_x, self.max_y, self.min_y,
48        #self.min_bed print self.max_bed, self.range_x, self.range_y,
49        #self.range_xy, self.range_bed
50        #print 'Calculating triangles'
51
52        self.pos     = zeros( (6*len(domain),3), Float)
53        self.colour  = zeros( (6*len(domain),3), Float)
54        self.normals = zeros( (6*len(domain),3), Float)
55
56        #print 'keys',self.domain.quantities.keys()
57        #print 'shape of stage',shape(self.stage)
58
59        self.border_model = curve(frame = self.frame,
60                                  pos=[(0,0),(0,1),(1,1),(1,0),(0,0)])
61        self.timer=label(pos=(0.75,0.75,0.5),text='Time=%10.5e'%self.domain.time)
62
63        #scene.autoscale=0
64        #self.update_all()
65        #scene.uniform=0
66        self.setup_range_z()
67
68    def setup_range_z(self,qname1='elevation',qname2='stage'):
69
70        #print qname1
71        #print qname2
72        self.range_z = 1.0e-10
73        try:
74            q1 = self.domain.quantities[qname1].vertex_values
75            max_z = max(max(q1))
76            min_z = min(min(q1))
77            print max_z, min_z
78            self.range_z = max(self.range_z, max_z - min_z)
79        except:
80            print 'could not find range of '+qname1
81            pass
82
83        try:
84            q2 = self.domain.quantities[qname2].vertex_values
85            max_z = max(max(q2))
86            min_z = min(min(q2))
87            print max_z, min_z
88            self.range_z = max(self.range_z, max_z - min_z)
89        except:
90            print 'Visualisation: could not find range of '+qname2
91            pass
92
93    def update_all(self):
94
95        self.update_quantity('elevation',elevation_color)
96        self.update_quantity('stage',stage_color)
97
98    def update_timer(self):
99        self.timer.text='Time=%10.5e'%self.domain.time
100
101
102    def update_quantity(self,qname,qcolor=None):
103
104        #print 'update '+qname+' arrays'
105
106        if qcolor is None:
107            if qname=='elevation':
108                qcolor = elevation_color
109
110            if qname=='stage':
111                qcolor = stage_color
112
113
114        try:
115            self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor)
116
117            #print 'update bed image'
118            self.z_models[qname].pos    = self.pos
119            self.z_models[qname].color  = self.colour
120            self.z_models[qname].normal = self.normals
121        except:
122            print 'Visualisation: Could not update quantity '+qname
123            pass
124
125    def update_quantity_color(self,qname,qcolor=None):
126
127        #print 'update '+qname+' arrays'
128
129        if qcolor is None:
130            if qname=='elevation':
131                qcolor = elevation_color
132
133            if qname=='stage':
134                qcolor = stage_color
135
136
137        #try:
138        self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor)
139
140            #print 'update bed image'
141        self.z_models[qname].pos    = self.pos
142        self.z_models[qname].color  = self.colour
143        self.z_models[qname].normal = self.normals
144        #except:
145        #    print 'Visualisation: Could not update quantity '+qname
146        #    pass
147
148
149    def update_arrays(self,quantity,qcolor):
150
151        col = asarray(qcolor)
152
153        N = len(self.domain)
154
155        scale_z = self.scale_z
156        min_x = self.min_x
157        min_y = self.min_y
158        range_xy = self.range_xy
159        range_z = self.range_z
160
161        vertices = self.vertices
162        pos      = self.pos
163        normals  = self.normals
164        colour   = self.colour
165
166        try:
167            update_arrays_weave(vertices,quantity,col,pos,normals,colour,N,
168                                min_x,min_y,range_xy,range_z,scale_z)
169        except:
170            update_arrays_python(vertices,quantity,col,pos,normals,colour,N,
171                                 min_x,min_y,range_xy,range_z,scale_z)
172
173
174    def update_arrays_color(self,quantity,qcolor):
175
176        col = asarray(qcolor)
177
178        N = len(self.domain)
179
180        scale_z = self.scale_z
181        min_x = self.min_x
182        min_y = self.min_y
183        range_xy = self.range_xy
184        range_z = self.range_z
185
186        vertices = self.vertices
187        pos      = self.pos
188        normals  = self.normals
189        colour   = self.colour
190
191        #try:
192        update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
193                                min_x,min_y,range_xy,range_z,scale_z)
194        #except:
195        #    update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
196        #                         min_x,min_y,range_xy,range_z,scale_z)
197
198
199#==================================================================================
200
201def  update_arrays_python(vertices,quantity,col,pos,normals,colour,N,
202                          min_x,min_y,range_xy,range_z,scale_z):
203
204    from math import sqrt
205    normal = zeros( 3, Float)
206    v = zeros( (3,3), Float)
207    s = 1.0
208
209    for i in range( N ):
210
211        for j in range(3):
212            v[j,0] = (vertices[i,2*]-min_x)/range_xy
213            v[j,1] = (vertices[i,2*j+1]-min_y)/range_xy
214            v[j,2] = quantity[i,j]/range_z*scale_z*0.5
215
216        v10 = v[1,:]-v[0,:]
217        v20 = v[2,:]-v[0,:]
218
219        normal[0] = v10[1]*v20[2] - v20[1]*v10[2]
220        normal[1] = v10[2]*v20[0] - v20[2]*v10[0]
221        normal[2] = v10[0]*v20[1] - v20[0]*v10[1]
222
223        norm = sqrt( normal[0]**2 + normal[1]**2 + normal[2]**2)
224
225        normal[0] = normal[0]/norm
226        normal[1] = normal[1]/norm
227        normal[2] = normal[2]/norm
228
229        pos[6*,:] = v[0,:]
230        pos[6*i+1,:] = v[1,:]
231        pos[6*i+2,:] = v[2,:]
232        pos[6*i+3,:] = v[0,:]
233        pos[6*i+4,:] = v[2,:]
234        pos[6*i+5,:] = v[1,:]
235
236
237
238        colour[6*,:] = s*col
239        colour[6*i+1,:] = s*col
240        colour[6*i+2,:] = s*col
241        colour[6*i+3,:] = s*col
242        colour[6*i+4,:] = s*col
243        colour[6*i+5,:] = s*col
244
245        s =  15.0/8.0 - s
246
247        normals[6*,:] = normal
248        normals[6*i+1,:] = normal
249        normals[6*i+2,:] = normal
250        normals[6*i+3,:] = -normal
251        normals[6*i+4,:] = -normal
252        normals[6*i+5,:] = -normal
253
254
255
256def  update_arrays_weave(vertices,quantity,col,pos,normals,colour,
257                         N,min_x,min_y,range_xy,range_z,scale_z):
258
259    import weave
260    from weave import converters
261
262    from math import sqrt
263    normal = zeros( 3, Float)
264    v = zeros( (3,3), Float)
265    v10 = zeros( 3, Float)
266    v20 = zeros( 3, Float)
267
268    code1 = """
269        double s = 1.0;
270
271        for (int i=0; i<N ; i++ ) {
272            for (int j=0; j<3 ; j++) {
273                v(j,0) = (vertices(i,2*j  )-min_x)/range_xy;
274                v(j,1) = (vertices(i,2*j+1)-min_y)/range_xy;
275                v(j,2) = quantity(i,j)/range_z*scale_z*0.5;
276            }
277
278
279            for (int j=0; j<3; j++) {
280                v10(j) = v(1,j)-v(0,j);
281                v20(j) = v(2,j)-v(0,j);
282            }
283
284            normal(0) = v10(1)*v20(2) - v20(1)*v10(2);
285            normal(1) = v10(2)*v20(0) - v20(2)*v10(0);
286            normal(2) = v10(0)*v20(1) - v20(0)*v10(1);
287
288            double norm =  sqrt(normal(0)*normal(0) + normal(1)*normal(1) + normal(2)*normal(2));
289
290            normal(0) = normal(0)/norm;
291            normal(1) = normal(1)/norm;
292            normal(2) = normal(2)/norm;
293
294            for (int j=0; j<3; j++) {
295                pos(6*i  ,j) = v(0,j);
296                pos(6*i+1,j) = v(1,j);
297                pos(6*i+2,j) = v(2,j);
298                pos(6*i+3,j) = v(0,j);
299                pos(6*i+4,j) = v(2,j);
300                pos(6*i+5,j) = v(1,j);
301
302                colour(6*i  ,j) = s*col(j);
303                colour(6*i+1,j) = s*col(j);
304                colour(6*i+2,j) = s*col(j);
305                colour(6*i+3,j) = s*col(j);
306                colour(6*i+4,j) = s*col(j);
307                colour(6*i+5,j) = s*col(j);
308
309                normals(6*i  ,j) = normal(j);
310                normals(6*i+1,j) = normal(j);
311                normals(6*i+2,j) = normal(j);
312                normals(6*i+3,j) = -normal(j);
313                normals(6*i+4,j) = -normal(j);
314                normals(6*i+5,j) = -normal(j);
315                }
316
317            s =  15.0/8.0 - s;
318        }
319
320        """
321
322    weave.inline(code1, ['vertices','quantity','col','pos','normals','colour','N',
323                         'min_x','min_y','range_xy','range_z','scale_z','v','normal','v10','v20'],
324                 type_converters = converters.blitz, compiler='gcc');
325
326def  update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
327                          min_x,min_y,range_xy,range_z,scale_z):
328
329    from math import sqrt
330    normal = zeros( 3, Float)
331    v = zeros( (3,3), Float)
332    s = 1.0
333
334    for i in range( N ):
335
336        for j in range(3):
337            v[j,0] = (vertices[i,2*]-min_x)/range_xy
338            v[j,1] = (vertices[i,2*j+1]-min_y)/range_xy
339            v[j,2] = quantity[i,j]/range_z*scale_z*0.5
340
341        v10 = v[1,:]-v[0,:]
342        v20 = v[2,:]-v[0,:]
343
344        normal[0] = v10[1]*v20[2] - v20[1]*v10[2]
345        normal[1] = v10[2]*v20[0] - v20[2]*v10[0]
346        normal[2] = v10[0]*v20[1] - v20[0]*v10[1]
347
348        norm = sqrt( normal[0]**2 + normal[1]**2 + normal[2]**2)
349
350        normal[0] = normal[0]/norm
351        normal[1] = normal[1]/norm
352        normal[2] = normal[2]/norm
353
354        pos[6*,:] = v[0,:]
355        pos[6*i+1,:] = v[1,:]
356        pos[6*i+2,:] = v[2,:]
357        pos[6*i+3,:] = v[0,:]
358        pos[6*i+4,:] = v[2,:]
359        pos[6*i+5,:] = v[1,:]
360
361        q0 = quantity[i,0]/0.4
362        q1 = quantity[i,1]/0.4
363        q2 = quantity[i,2]/0.4
364
365        q0r = min(q0,0.0)
366        q1r = min(q1,0.0)
367        q2r = min(q2,0.0)
368        q0b = max(q0,0.0)
369        q1b = max(q1,0.0)
370        q2b = max(q2,0.0)
371
372
373        colour[6*,:] = s*array([0.3 - q0r, 0.3 - q0, 0.3 + q0b])
374        colour[6*i+1,:] = s*array([0.3 - q1r, 0.3 - q1, 0.3 + q1b])
375        colour[6*i+2,:] = s*array([0.3 - q2r, 0.3 - q2, 0.3 + q2b])
376        colour[6*i+3,:] = s*array([0.3 - q0r, 0.3 - q0, 0.3 + q0b])
377        colour[6*i+4,:] = s*array([0.3 - q2r, 0.3 - q2, 0.3 + q2b])
378        colour[6*i+5,:] = s*array([0.3 - q1r, 0.3 - q1, 0.3 + q1b])
379
380        s =  15.0/8.0 - s
381
382        normals[6*,:] = normal
383        normals[6*i+1,:] = normal
384        normals[6*i+2,:] = normal
385        normals[6*i+3,:] = -normal
386        normals[6*i+4,:] = -normal
387        normals[6*i+5,:] = -normal
388
389
390
391def  update_arrays_color_weave(vertices,quantity,col,pos,normals,colour,
392                         N,min_x,min_y,range_xy,range_z,scale_z):
393
394    import weave
395    from weave import converters
396
397    from math import sqrt
398    normal = zeros( 3, Float)
399    v = zeros( (3,3), Float)
400    v10 = zeros( 3, Float)
401    v20 = zeros( 3, Float)
402
403    code1 = """
404        double s = 1.0;
405
406        for (int i=0; i<N ; i++ ) {
407            for (int j=0; j<3 ; j++) {
408                v(j,0) = (vertices(i,2*j  )-min_x)/range_xy;
409                v(j,1) = (vertices(i,2*j+1)-min_y)/range_xy;
410                v(j,2) = quantity(i,j)/range_z*scale_z*0.5;
411            }
412
413
414            for (int j=0; j<3; j++) {
415                v10(j) = v(1,j)-v(0,j);
416                v20(j) = v(2,j)-v(0,j);
417            }
418
419            normal(0) = v10(1)*v20(2) - v20(1)*v10(2);
420            normal(1) = v10(2)*v20(0) - v20(2)*v10(0);
421            normal(2) = v10(0)*v20(1) - v20(0)*v10(1);
422
423            double norm =  sqrt(normal(0)*normal(0) + normal(1)*normal(1) + normal(2)*normal(2));
424
425            normal(0) = normal(0)/norm;
426            normal(1) = normal(1)/norm;
427            normal(2) = normal(2)/norm;
428
429
430
431            for (int j=0; j<3; j++) {
432                pos(6*i  ,j) = v(0,j);
433                pos(6*i+1,j) = v(1,j);
434                pos(6*i+2,j) = v(2,j);
435                pos(6*i+3,j) = v(0,j);
436                pos(6*i+4,j) = v(2,j);
437                pos(6*i+5,j) = v(1,j);
438
439                normals(6*i  ,j) = normal(j);
440                normals(6*i+1,j) = normal(j);
441                normals(6*i+2,j) = normal(j);
442                normals(6*i+3,j) = -normal(j);
443                normals(6*i+4,j) = -normal(j);
444                normals(6*i+5,j) = -normal(j);
445                }
446
447            s = 0.2+fabs((v(0,2)+v(1,2)+v(2,2))/3.0);
448            for (int j=0; j<3; j++) {
449                colour(6*i  ,j) = s*col(j);
450                colour(6*i+1,j) = s*col(j);
451                colour(6*i+2,j) = s*col(j);
452                colour(6*i+3,j) = s*col(j);
453                colour(6*i+4,j) = s*col(j);
454                colour(6*i+5,j) = s*col(j);
455                }
456            s =  15.0/8.0 - s;
457        }
458
459        """
460
461    weave.inline(code1, ['vertices','quantity','col','pos','normals','colour','N',
462                         'min_x','min_y','range_xy','range_z','scale_z','v','normal','v10','v20'],
463                 type_converters = converters.blitz, compiler='gcc');
464
465scene.width = 1000
466scene.height = 800
467
468#Original
469scene.center = (0.5,0.5,0.0)
470scene.forward = vector(-0.0, 0.5, -0.8)
471
472#Temporary (for bedslope)
473#scene.forward = vector(0.0006, 0.7, -0.03)
474
475
476#Temporary for hackett - begin
477#scene.autoscale = 0
478#scene.scale = (0.002, 0.002, 0.01) #Scale z so that countours stand out more
479#scene.center = (300.0,500.0,-10)
480#Temporary for hackett - end
481
482
483scene.ambient = 0.4
484#scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
485#               (-0.2, 0.2, 0.1)]
486
487scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4)]
488
489
490
491def create_surface(domain,scale_z=0.5,range_z=None):
492
493    surface = Surface(domain,scale_z=0.5,range_z=range_z)
494    domain.surface = surface
495
496    surface.update_quantity('elevation')
497    surface.update_quantity('stage')
498    surface.update_timer()
499
500def update(domain):
501
502    surface = domain.surface
503    if domain.visualise_color_stage:
504        surface.update_quantity_color('stage')
505    else:
506        surface.update_quantity('stage')
507
508    if domain.visualise_timer:
509        surface.update_timer()
Note: See TracBrowser for help on using the repository browser.