Changeset 1266

Timestamp:

May 2, 2005, 12:02:12 AM (20 years ago)

Author:

steve

Message:

Played with weave to convert update of new visualisation code.

Location:

inundation/ga/storm_surge

Files:

• pyvolution/netherlands.py (modified) (12 diffs)
• pyvolution/realtime_visualisation_new.py (modified) (9 diffs)
• pyvolution/view_tsh.py (added)
• steve/get_triangle_data_new.py (added)
• zeus/anuga-workspace.zwi (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -7 +7 @@
 
 ######################
-# Module imports 
+# Module imports
 #
 from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
@@ -14 +14 @@
 from mesh_factory import rectangular
 from Numeric import array
-    
+
 
 class Weir:
@@ -20 +20 @@
     x,y are assumed to be in the unit square
     """
-    
+
     def __init__(self, stage):
         self.inflow_stage = stage
 
-    def __call__(self, x, y):    
-        from Numeric import zeros, Float 
-    
+    def __call__(self, x, y):
+        from Numeric import zeros, Float
+
         N = len(x)
-        assert N == len(y)    
+        assert N == len(y)
 
         z = zeros(N, Float)
@@ -38 +38 @@
                 #z[i] = -x[i]/5
                 z[i] = -x[i]/20
-                
-            
+
+
             #Weir
             if x[i] > 0.3 and x[i] < 0.4:
@@ -53 +53 @@
                 #z[i] = -x[i]/5
                 z[i] = -x[i]/20
-            
+
             #Poles
             #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
@@ -65 +65 @@
             #Wall
             if x[i] > 0.995:
-                z[i] = -x[i]/20+0.3                
+                z[i] = -x[i]/20+0.3
 
         return z/2
-        
-    
-    
+
+
+
 ######################
 # Domain
@@ -79 +79 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 30
+N = 50
 
 #N = 15
@@ -89 +89 @@
 #Create shallow water domain
 domain = Domain(points, vertices, boundary)
-    
+
 domain.check_integrity()
 domain.default_order = 2
@@ -101 +101 @@
 
 
-if N > 40:
+if N > 150:
     domain.visualise = False
     domain.checkpoint = False
@@ -113 +113 @@
     domain.visualise = True
     domain.checkpoint = False
-    domain.store = False    
+    domain.store = False
 
-    
+
 #Set bed-slope and friction
-inflow_stage = 0.08
+inflow_stage = 0.1
 manning = 0.02
 Z = Weir(inflow_stage)
@@ -139 +139 @@
 #Set boundary conditions
 domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
-                  
+
 
 ######################
@@ -151 +151 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.1, finaltime = 7.0):
+for t in domain.evolve(yieldstep = 0.005, finaltime = 10.0):
     domain.write_time()
-    
-print 'That took %.2f seconds' %(time.time()-t0)    
-    
+
+
+print 'That took %.2f seconds' %(time.time()-t0)

inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

@@ -2 +2 @@
 
 class Surface:
-    def __init__(self,domain,scale=1.0/1.0):
+    def __init__(self,domain,scale_z=1.0):
         """Create visualisation of domain
         """
@@ -10 +10 @@
         self.stage_model = faces(frame=self.frame)
         self.domain = domain
-        self.scale  = scale
+        self.scale_z  = scale_z
 
         self.vertices = domain.vertex_coordinates
@@ -29 +29 @@
         self.range_bed = self.max_bed - self.min_bed
 
-        print self.max_x, self.min_x, self.max_y, self.min_y, self.min_bed
-        print self.max_bed, self.range_x, self.range_y, self.range_xy, self.range_bed
-
-        print 'Calculating triangles'
+        #print self.max_x, self.min_x, self.max_y, self.min_y,
+        #self.min_bed print self.max_bed, self.range_x, self.range_y,
+        #self.range_xy, self.range_bed
+        #print 'Calculating triangles'
 
         self.pos     = zeros( (6*len(domain),3), Float)
@@ -38 +38 @@
         self.normals = zeros( (6*len(domain),3), Float)
 
-        print 'keys',self.domain.quantities.keys()
-        print 'shape of stage',shape(self.stage)
+        #print 'keys',self.domain.quantities.keys()
+        #print 'shape of stage',shape(self.stage)
 
         self.update_all()
@@ -50 +50 @@
     def update_bed(self):
 
-        print 'update bed arrays'
+        #print 'update bed arrays'
         c0 = 0.3
         c1 = 0.3
         c2 = 0.3
         self.update_arrays(self.bed,  (c0,c1,c2) )
-        print 'update bed image'
+
+        #print 'update bed image'
         self.bed_model.pos    = self.pos
         self.bed_model.color  = self.colour
@@ -65 +66 @@
         c1 = 0.4
         c2 = 0.99
-        print 'update stage arrays'
+        #print 'update stage arrays'
         self.update_arrays(self.stage, (c0,c1,c2) )
-        print 'update stage image'
+        #print 'update stage image'
         self.stage_model.pos    = self.pos
         self.stage_model.color  = self.colour
         self.stage_model.normal = self.normals
 
-    def calc_normal(self,v):
-
-        from math import sqrt
-
-        cross = zeros(3,Float)
-        normal = zeros(3,Float)
+
+
+    def update_arrays(self,quantity,qcolor):
+
+        col = asarray(qcolor)
+
+        N = len(self.domain)
+
+        scale_z = self.scale_z
+        min_x = self.min_x
+        min_y = self.min_y
+        range_xy = self.range_xy
+
+        vertices = self.vertices
+        pos      = self.pos
+        normals  = self.normals
+        colour   = self.colour
+
+        try:
+            update_arrays_weave(vertices,quantity,col,pos,normals,colour,N,min_x,min_y,range_xy,scale_z)
+        except:
+            update_arrays_python(vertices,quantity,col,pos,normals,colour,N,min_x,min_y,range_xy,scale_z)
+
+
+
+def  update_arrays_python(vertices,quantity,col,pos,normals,colour,N,min_x,min_y,range_xy,scale_z):
+
+    from math import sqrt
+    normal = zeros( 3, Float)
+    v = zeros( (3,3), Float)
+    s = 1.0
+
+    for i in range( N ):
+
+        for j in range(3):
+            v[j,0] = (vertices[i,2*j  ]-min_x)/range_xy
+            v[j,1] = (vertices[i,2*j+1]-min_y)/range_xy
+            v[j,2] = quantity[i,j]*scale_z
 
         v10 = v[1,:]-v[0,:]
@@ -92 +125 @@
         normal[2] = normal[2]/norm
 
-        return normal
-
-
-
-    def update_arrays(self, quantity, qcolor):
-
-        col = asarray(qcolor)
-        v = zeros( (3,3), Float)
-        normal = zeros( 3, Float)
-        N = len(self.domain)
-        scale = self.scale
-
-
-
-        for i in range( N ):
-
-            for j in range(3):
-                v[j,0] = 2.0*(self.vertices[i,2*j  ]-self.min_x)/self.range_xy-1.0
-                v[j,1] = 2.0*(self.vertices[i,2*j+1]-self.min_y)/self.range_xy-1.0
-                v[j,2] = quantity[i,j]*scale
-
-
-
-            normal = self.calc_normal(v)
-
-            self.pos[6*i  ,:] = v[0,:]
-            self.pos[6*i+1,:] = v[1,:]
-            self.pos[6*i+2,:] = v[2,:]
-            self.pos[6*i+3,:] = v[0,:]
-            self.pos[6*i+4,:] = v[2,:]
-            self.pos[6*i+5,:] = v[1,:]
-
-
-
-            self.colour[6*i  ,:] = col
-            self.colour[6*i+1,:] = col
-            self.colour[6*i+2,:] = col
-            self.colour[6*i+3,:] = col
-            self.colour[6*i+4,:] = col
-            self.colour[6*i+5,:] = col
-
-            self.normals[6*i  ,:] = normal
-            self.normals[6*i+1,:] = normal
-            self.normals[6*i+2,:] = normal
-            self.normals[6*i+3,:] = -normal
-            self.normals[6*i+4,:] = -normal
-            self.normals[6*i+5,:] = -normal
-
-
-
-
-#scene.width = 1000
-#scene.height = 800
+        pos[6*i  ,:] = v[0,:]
+        pos[6*i+1,:] = v[1,:]
+        pos[6*i+2,:] = v[2,:]
+        pos[6*i+3,:] = v[0,:]
+        pos[6*i+4,:] = v[2,:]
+        pos[6*i+5,:] = v[1,:]
+
+
+
+        colour[6*i  ,:] = s*col
+        colour[6*i+1,:] = s*col
+        colour[6*i+2,:] = s*col
+        colour[6*i+3,:] = s*col
+        colour[6*i+4,:] = s*col
+        colour[6*i+5,:] = s*col
+
+        s =  15.0/8.0 - s
+
+        normals[6*i  ,:] = normal
+        normals[6*i+1,:] = normal
+        normals[6*i+2,:] = normal
+        normals[6*i+3,:] = -normal
+        normals[6*i+4,:] = -normal
+        normals[6*i+5,:] = -normal
+
+
+
+def  update_arrays_weave(vertices,quantity,col,pos,normals,colour,N,min_x,min_y,range_xy,scale_z):
+
+    import weave
+    from weave import converters
+
+    from math import sqrt
+    normal = zeros( 3, Float)
+    v = zeros( (3,3), Float)
+    v10 = zeros( 3, Float)
+    v20 = zeros( 3, Float)
+
+    code1 = """
+        double s = 1.0;
+
+        for (int i=0; i<N ; i++ ) {
+            for (int j=0; j<3 ; j++) {
+                v(j,0) = (vertices(i,2*j  )-min_x)/range_xy;
+                v(j,1) = (vertices(i,2*j+1)-min_y)/range_xy;
+                v(j,2) = quantity(i,j)*scale_z;
+            }
+
+
+            for (int j=0; j<3; j++) {
+                v10(j) = v(1,j)-v(0,j);
+                v20(j) = v(2,j)-v(0,j);
+            }
+
+            normal(0) = v10(1)*v20(2) - v20(1)*v10(2);
+            normal(1) = v10(2)*v20(0) - v20(2)*v10(0);
+            normal(2) = v10(0)*v20(1) - v20(0)*v10(1);
+
+            double norm =  sqrt(normal(0)*normal(0) + normal(1)*normal(1) + normal(2)*normal(2));
+
+            normal(0) = normal(0)/norm;
+            normal(1) = normal(1)/norm;
+            normal(2) = normal(2)/norm;
+
+            for (int j=0; j<3; j++) {
+                pos(6*i  ,j) = v(0,j);
+                pos(6*i+1,j) = v(1,j);
+                pos(6*i+2,j) = v(2,j);
+                pos(6*i+3,j) = v(0,j);
+                pos(6*i+4,j) = v(2,j);
+                pos(6*i+5,j) = v(1,j);
+
+                colour(6*i  ,j) = s*col(j);
+                colour(6*i+1,j) = s*col(j);
+                colour(6*i+2,j) = s*col(j);
+                colour(6*i+3,j) = s*col(j);
+                colour(6*i+4,j) = s*col(j);
+                colour(6*i+5,j) = s*col(j);
+
+                normals(6*i  ,j) = normal(j);
+                normals(6*i+1,j) = normal(j);
+                normals(6*i+2,j) = normal(j);
+                normals(6*i+3,j) = -normal(j);
+                normals(6*i+4,j) = -normal(j);
+                normals(6*i+5,j) = -normal(j);
+                }
+
+            s =  15.0/8.0 - s;
+        }
+        """
+
+    weave.inline(code1, ['vertices','quantity','col','pos','normals','colour','N',
+                         'min_x','min_y','range_xy','scale_z','v','normal','v10','v20'],
+                 type_converters = converters.blitz, compiler='gcc');
+
+
+scene.width = 1000
+scene.height = 800
 
 #Original
-#scene.center = (0.5,0.5,0)
-#scene.forward = vector(-0.1, 0.5, -0.5)
+scene.center = (0.5,0.5,0)
+scene.forward = vector(-0.1, 0.5, -0.5)
 
 #Temporary (for bedslope)
@@ -161 +239 @@
 
 
-#scene.ambient = 0.4
+scene.ambient = 0.4
 #scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
 #               (-0.2, 0.2, 0.1)]
 
+scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4)]
+
+
 
 def create_surface(domain):
@@ -173 +254 @@
 def update(domain):
 
+    #print 'start visual update'
     surface = domain.surface
     surface.update_stage()
-
-
-
+    #print 'end visual update'
+

inundation/ga/storm_surge/zeus/anuga-workspace.zwi

@@ -2 +2 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>steve</active>
+    <active>pyvolution</active>
     <project name="pyvolution">pyvolution.zpi</project>
     <project name="steve">steve.zpi</project>