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

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