Changeset 1290

Timestamp:

May 6, 2005, 6:24:03 PM (19 years ago)

Author:

steve

Message:

visualisation stuff

Location:

inundation/ga/storm_surge

Files:

• analytical solutions/Hobart.py (modified) (1 diff)
• pyvolution/domain.py (modified) (3 diffs)
• pyvolution/mesh.py (modified) (4 diffs)
• pyvolution/netherlands-chris.py (added)
• pyvolution/netherlands.py (modified) (1 diff)
• pyvolution/quantity.py (modified) (60 diffs)
• pyvolution/realtime_visualisation_new.py (modified) (7 diffs)
• pyvolution/view_tsh.py (modified) (1 diff)
• steve/Eden_tide_Sept03.dat (deleted)
• steve/Merimbula_17540_refined.tsh (deleted)
• steve/Merimbula_Weather_data_Sept03_m_per_s.dat (deleted)
• steve/Merimbula_gauge_results.txt (deleted)
• steve/XYmerimbula_bathymetry.xya (deleted)
• steve/meri0.xya (deleted)
• steve/merimbula_10785.tsh (deleted)
• steve/merimbula_interpolated.tsh (deleted)
• steve/run_least_squares_merimbula.py (deleted)
• steve/run_merimbula_lake.py (modified) (3 diffs)
• steve/weed_zone.000 (deleted)
• steve/weed_zone.001 (deleted)
• steve/weed_zone.002 (deleted)
• steve/weed_zone.003 (deleted)
• steve/weed_zone.004 (deleted)
• steve/weed_zone.005 (deleted)
• steve/weed_zone.006 (deleted)
• steve/weed_zone.007 (deleted)
• steve/weed_zone.008 (deleted)
• steve/weed_zone.009 (deleted)
• steve/weed_zone.010 (deleted)
• steve/weed_zone.011 (deleted)
• steve/weed_zone.012 (deleted)
• steve/weed_zone.013 (deleted)
• steve/weed_zone.014 (deleted)
• steve/weed_zone.015 (deleted)
• steve/weed_zone.016 (deleted)
• steve/weed_zone.017 (deleted)
• steve/weed_zone.018 (deleted)
• steve/weed_zone.019 (deleted)
• steve/weed_zone.020 (deleted)
• steve/weed_zone.021 (deleted)
• steve/weed_zone.022 (deleted)
• steve/weed_zone.023 (deleted)
• steve/weed_zone.024 (deleted)
• steve/weed_zone.025 (deleted)
• steve/weed_zone.026 (deleted)
• steve/weed_zone.027 (deleted)
• steve/weed_zone.028 (deleted)
• steve/weed_zone.029 (deleted)
• steve/weed_zone.030 (deleted)
• steve/weed_zone.031 (deleted)
• steve/weed_zone.032 (deleted)
• steve/weed_zone.033 (deleted)
• steve/weed_zone.034 (deleted)
• steve/weed_zone.035 (deleted)
• steve/weed_zone.036 (deleted)
• steve/weed_zone.037 (deleted)
• steve/weed_zone.038 (deleted)
• steve/weed_zone.039 (deleted)
• steve/weed_zone.040 (deleted)
• steve/weed_zone.041 (deleted)
• steve/weed_zone.042 (deleted)
• steve/weed_zone.043 (deleted)
• steve/weed_zone.044 (deleted)
• steve/weed_zone.045 (deleted)
• steve/weed_zone.046 (deleted)
• steve/weed_zone.047 (deleted)
• steve/weed_zone.all (deleted)
• zeus/anuga-workspace.zwi (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)
• zeus/steve.zpi (modified) (1 diff)

inundation/ga/storm_surge/analytical solutions/Hobart.py

@@ -42 +42 @@
 domain.default_order = 1
 domain.smooth = True
+domain.visualise = True
 
 

inundation/ga/storm_surge/pyvolution/domain.py

@@ -141 +141 @@
         X: Compatible list, Numeric array, const or function (see below)
         location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
-
-        In case of location == 'centroid' the dimension values must
+                  Permissible options are: vertices, edges, centroids
+
+        In case of location == 'centroids' the dimension values must
         be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3.
@@ -152 +152 @@
         internal ordering.
 
-        #FIXME (Ole): I suggest the following interface
-        set_quantity(name, X, location, region)
-        where
+        #FIXME (Ole): I suggest the following interface
+        set_quantity(name, X, location, region)
+        where
             name: Name of quantity
             X:
@@ -194 +194 @@
         name: Name of quantity
 
-        In case of location == 'centroid' the dimension values must
+        In case of location == 'centroids' the dimension values must
         be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3.

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -67 +67 @@
         """
 
-        
+
 
         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
@@ -113 +113 @@
                 #a triangle to the midpoint of the side of the triangle
                 #closest to the centroid
-                d0 = sqrt(sum( (centroid-m0)**2 ))
+                d0 = sqrt(sum( (centroid-m0)**2 ))
                 d1 = sqrt(sum( (centroid-m1)**2 ))
                 d2 = sqrt(sum( (centroid-m2)**2 ))
@@ -129 +129 @@
 
                 self.radii[i]=self.areas[i]/(2*(a+b+c))
-                
+
 
             #Initialise Neighbours (-1 means that it is a boundary neighbour)
@@ -446 +446 @@
             #Note: Could be arbitrary, but nice to have
             #a unique way of selecting
-            dist_A = sqrt(sum( (A-pmin)**2 ))
-            dist_B = sqrt(sum( (B-pmin)**2 ))
+            dist_A = sqrt(sum( (A-pmin)**2 ))
+            dist_B = sqrt(sum( (B-pmin)**2 ))
 
             #Find minimal point

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -151 +151 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.01, finaltime = 10.0):
+for t in domain.evolve(yieldstep = 0.2, finaltime = 10.0):
     domain.write_time()
 

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -6 +6 @@
 
    domain: Associated domain structure. Required.
-   
+
    vertex_values: N x 3 array of values at each vertex for each element.
                   Default None
@@ -28 +28 @@
 
         if vertex_values is None:
-            N = domain.number_of_elements            
+            N = domain.number_of_elements
             self.vertex_values = zeros((N, 3), Float)
-        else:    
+        else:
             self.vertex_values = array(vertex_values).astype(Float)
 
@@ -36 +36 @@
             assert V == 3,\
                    'Three vertex values per element must be specified'
-            
+
 
             msg = 'Number of vertex values (%d) must be consistent with'\
@@ -42 +42 @@
             msg += 'number of elements in specified domain (%d).'\
                    %domain.number_of_elements
-            
+
             assert N == domain.number_of_elements, msg
 
-        self.domain = domain    
+        self.domain = domain
 
         #Allocate space for other quantities
         self.centroid_values = zeros(N, Float)
-        self.edge_values = zeros((N, 3), Float)        
+        self.edge_values = zeros((N, 3), Float)
 
         #Intialise centroid and edge_values
@@ -56 +56 @@
     def __len__(self):
         return self.centroid_values.shape[0]
-    
+
     def interpolate(self):
         """Compute interpolated values at edges and centroid
@@ -62 +62 @@
         """
 
-        N = self.vertex_values.shape[0]        
+        N = self.vertex_values.shape[0]
         for i in range(N):
             v0 = self.vertex_values[i, 0]
             v1 = self.vertex_values[i, 1]
             v2 = self.vertex_values[i, 2]
-            
+
             self.centroid_values[i] = (v0 + v1 + v2)/3
 
@@ -77 +77 @@
         #(either from this module or C-extension)
         interpolate_from_vertices_to_edges(self)
-            
-            
+
+
     def set_values(self, X, location='vertices', indexes = None):
         """Set values for quantity
@@ -84 +84 @@
         X: Compatible list, Numeric array (see below), constant or function
         location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
+                  Permissible options are: vertices, edges, centroids
                   Default is "vertices"
 
-        In case of location == 'centroid' the dimension values must
+        In case of location == 'centroids' the dimension values must
         be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3
@@ -96 +96 @@
         If values are described a function, it will be evaluated at
         specified points
-        
+
         If indexex is not 'unique vertices' Indexes is the set of element ids
         that the operation applies to.
         If indexex is 'unique vertices' Indexes is the set of vertex ids
         that the operation applies to.
-        
-        
+
+
         If selected location is vertices, values for centroid and edges
         will be assigned interpolated values.
@@ -115 +115 @@
         if X is None:
             msg = 'Given values are None'
-            raise msg            
+            raise msg
 
         import types, Numeric
@@ -122 +122 @@
                                  'Indices must be a list or None'
 
-        
+
         if callable(X):
             #Use function specific method
-            self.set_function_values(X, location, indexes = indexes)        
+            self.set_function_values(X, location, indexes = indexes)
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
             if location == 'centroids':
@@ -134 +134 @@
                     for i in indexes:
                         self.centroid_values[i,:] = X
-                        
+
             elif location == 'edges':
                 if (indexes ==  None):
@@ -142 +142 @@
                     for i in indexes:
                         self.edge_values[i,:] = X
-                        
+
             elif location == 'unique vertices':
                 if (indexes ==  None):
@@ -151 +151 @@
                     for unique_vert_id in indexes:
                         triangles = self.domain.vertexlist[unique_vert_id]
-                        
+
                         #In case there are unused points
-                        if triangles is None: continue 
+                        if triangles is None: continue
 
                         #Go through all triangle, vertex pairs
@@ -159 +159 @@
                         for triangle_id, vertex_id in triangles:
                             self.vertex_values[triangle_id, vertex_id] = X
-                
+
                         #Intialise centroid and edge_values
                         self.interpolate()
@@ -168 +168 @@
                     #Brute force
                     for i_vertex in indexes:
-                        self.vertex_values[i_vertex,:] = X     
+                        self.vertex_values[i_vertex,:] = X
 
         else:
@@ -177 +177 @@
             #Intialise centroid and edge_values
             self.interpolate()
-            
+
         if location == 'centroids':
             #Extrapolate 1st order - to capture notion of area being specified
-            self.extrapolate_first_order()            
-          
+            self.extrapolate_first_order()
+
     def get_values(self, location='vertices', indexes = None):
         """get values for quantity
@@ -197 +197 @@
         (N if indexes = None).  Each value will be a list of the three
         vertex values for this quantity.
-        
+
         Indexes is the set of element ids that the operation applies to.
-        
+
         """
         from Numeric import take
@@ -205 +205 @@
         if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
             msg = 'Invalid location: %s' %location
-            raise msg         
-        
+            raise msg
+
         import types, Numeric
         assert type(indexes) in [types.ListType, types.NoneType,
                                  Numeric.ArrayType],\
                                  'Indices must be a list or None'
-            
+
         if location == 'centroids':
             if (indexes ==  None):
                 indexes = range(len(self))
-            return take(self.centroid_values,indexes)                       
+            return take(self.centroid_values,indexes)
         elif location == 'edges':
             if (indexes ==  None):
                 indexes = range(len(self))
-            return take(self.edge_values,indexes)                      
+            return take(self.edge_values,indexes)
         elif location == 'unique vertices':
             if (indexes ==  None):
@@ -227 +227 @@
             for unique_vert_id in indexes:
                 triangles = self.domain.vertexlist[unique_vert_id]
-                        
+
                 #In case there are unused points
-                if triangles is None: 
-                    msg = 'Unique vertex not associated with triangles' 
-                    raise msg 
+                if triangles is None:
+                    msg = 'Unique vertex not associated with triangles'
+                    raise msg
 
                 # Go through all triangle, vertex pairs
@@ -280 +280 @@
             raise 'Not implemented: %s' %location
 
-            
+
     def set_array_values(self, values, location='vertices', indexes = None):
         """Set values for quantity
@@ -288 +288 @@
         Permissible options are: vertices, edges, centroid, unique vertices
         Default is "vertices"
-        
+
         indexes - if this action is carried out on a subset of
         elements or unique vertices
         The element/unique vertex indexes are specified here.
-        
+
         In case of location == 'centroid' the dimension values must
         be a list of a Numerical array of length N, N being the number
@@ -315 +315 @@
             indexes = array(indexes).astype(Int)
             msg = 'Number of values must match number of indexes'
-            assert values.shape[0] == indexes.shape[0], msg 
-            
+            assert values.shape[0] == indexes.shape[0], msg
+
         N = self.centroid_values.shape[0]
-        
+
         if location == 'centroids':
             assert len(values.shape) == 1, 'Values array must be 1d'
@@ -325 +325 @@
                 msg = 'Number of values must match number of elements'
                 assert values.shape[0] == N, msg
-            
+
                 self.centroid_values = values
             else:
                 msg = 'Number of values must match number of indexes'
                 assert values.shape[0] == indexes.shape[0], msg
-                
+
                 #Brute force
                 for i in range(len(indexes)):
                     self.centroid_values[indexes[i]] = values[i]
-                    
+
         elif location == 'edges':
             assert len(values.shape) == 2, 'Values array must be 2d'
@@ -340 +340 @@
             msg = 'Number of values must match number of elements'
             assert values.shape[0] == N, msg
-            
-            msg = 'Array must be N x 3'            
+
+            msg = 'Array must be N x 3'
             assert values.shape[1] == 3, msg
-            
+
             self.edge_values = values
-            
+
         elif location == 'unique vertices':
             assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
@@ -360 +360 @@
                  #   self.set_vertex_values(values)
                 #else:
-                #    for element_index, value in map(None, indexes, values): 
+                #    for element_index, value in map(None, indexes, values):
                 #        self.vertex_values[element_index, :] = value
-                        
+
             elif len(values.shape) == 2:
                 #Vertex values are given as a triplet for each triangle
-                
-                msg = 'Array must be N x 3'            
+
+                msg = 'Array must be N x 3'
                 assert values.shape[1] == 3, msg
-                
+
                 if indexes == None:
                     self.vertex_values = values
                 else:
-                    for element_index, value in map(None, indexes, values): 
+                    for element_index, value in map(None, indexes, values):
                         self.vertex_values[element_index] = value
-            else:    
+            else:
                 msg = 'Values array must be 1d or 2d'
                 raise msg
-            
+
 
     # FIXME have a get_vertex_values as well, so the 'stage' quantity can be
@@ -387 +387 @@
         one value for each row in vertexlist.
 
-        indexes is the list of vertex_id's that will be set. 
+        indexes is the list of vertex_id's that will be set.
 
         Note: Functions not allowed
@@ -409 +409 @@
         for i_index,unique_vert_id in enumerate(vertex_list):
             triangles = self.domain.vertexlist[unique_vert_id]
-                
+
             if triangles is None: continue #In case there are unused points
 
@@ -416 +416 @@
             for triangle_id, vertex_id in triangles:
                 self.vertex_values[triangle_id, vertex_id] = A[i_index]
-                
+
         #Intialise centroid and edge_values
         self.interpolate()
- 
+
     def smooth_vertex_values(self, value_array='field_values',
                              precision = None):
@@ -432 +432 @@
         from Numeric import concatenate, zeros, Float, Int, array, reshape
 
-        
+
         A,V = self.get_vertex_values(xy=False,
                                      value_array=value_array,
@@ -450 +450 @@
         elif value_array == 'conserved_quantities':
             Volume.interpolate_conserved_quantities()
-        
-    
+
+
     #Method for outputting model results
     #FIXME: Split up into geometric and numeric stuff.
@@ -457 +457 @@
     #FIXME: STill remember to move XY to mesh
     def get_vertex_values(self,
-                          xy=True, 
+                          xy=True,
                           smooth = None,
                           precision = None,
@@ -465 +465 @@
         The vertex values are returned as one sequence in the 1D float array A.
         If requested the coordinates will be returned in 1D arrays X and Y.
-        
+
         The connectivity is represented as an integer array, V, of dimension
         M x 3, where M is the number of volumes. Each row has three indices
@@ -474 +474 @@
         reduction operator. In this case vertex coordinates will be
         de-duplicated.
-        
+
         If no smoothings is required, vertex coordinates and values will
         be aggregated as a concatenation of values at
         vertices 0, vertices 1 and vertices 2
 
-        
+
         Calling convention
         if xy is True:
            X,Y,A,V = get_vertex_values
         else:
-           A,V = get_vertex_values         
-           
+           A,V = get_vertex_values
+
         """
 
@@ -499 +499 @@
         if reduction is None:
             reduction = self.domain.reduction
-            
+
         #Create connectivity
-                        
+
         if smooth == True:
 
@@ -511 +511 @@
             for k in range(N):
                 L = self.domain.vertexlist[k]
-                
+
                 #Go through all triangle, vertex pairs
                 #contributing to vertex k and register vertex value
 
                 if L is None: continue #In case there are unused points
-                
+
                 contributions = []
                 for volume_id, vertex_id in L:
@@ -522 +522 @@
                     contributions.append(v)
 
-                A[k] = reduction(contributions)    
+                A[k] = reduction(contributions)
 
 
@@ -528 +528 @@
                 X = self.domain.coordinates[:,0].astype(precision)
                 Y = self.domain.coordinates[:,1].astype(precision)
-                
+
                 return X, Y, A, V
             else:
@@ -540 +540 @@
             M = 3*m        #Total number of unique vertices
             V = reshape(array(range(M)).astype(Int), (m,3))
-            
+
             A = self.vertex_values.flat
 
-            #Do vertex coordinates    
-            if xy is True:                
+            #Do vertex coordinates
+            if xy is True:
                 C = self.domain.get_vertex_coordinates()
 
                 X = C[:,0:6:2].copy()
-                Y = C[:,1:6:2].copy()                
-
-                return X.flat, Y.flat, A, V           
+                Y = C[:,1:6:2].copy()
+
+                return X.flat, Y.flat, A, V
             else:
                 return A, V
 
-            
+
     def extrapolate_first_order(self):
         """Extrapolate conserved quantities from centroid to
@@ -560 +560 @@
         first order scheme.
         """
-        
+
         qc = self.centroid_values
         qv = self.vertex_values
@@ -566 +566 @@
         for i in range(3):
             qv[:,i] = qc
-            
+
 
     def get_integral(self):
@@ -574 +574 @@
         for k in range(self.domain.number_of_elements):
             area = self.domain.areas[k]
-            qc = self.centroid_values[k]     
+            qc = self.centroid_values[k]
             integral += qc*area
-            
-        return integral    
-        
+
+        return integral
+
 
 class Conserved_quantity(Quantity):
@@ -590 +590 @@
     def __init__(self, domain, vertex_values=None):
         Quantity.__init__(self, domain, vertex_values)
-        
+
         from Numeric import zeros, Float
 
@@ -603 +603 @@
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
-                
+
 
     def update(self, timestep):
@@ -615 +615 @@
         #(either from this module or C-extension)
         return compute_gradients(self)
-            
+
 
     def limit(self):
@@ -621 +621 @@
         #(either from this module or C-extension)
         limit(self)
-        
-        
+
+
     def extrapolate_second_order(self):
         #Call correct module function
@@ -629 +629 @@
 
 
-def update(quantity, timestep):    
+def update(quantity, timestep):
     """Update centroid values based on values stored in
     explicit_update and semi_implicit_update as well as given timestep
@@ -642 +642 @@
         domain.quantities['ymomentum'].explicit_update = ...
 
-        
+
 
     Explicit terms must have the form
-    
+
         G(q, t)
-    
+
     and explicit scheme is
-    
+
        q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
 
 
     Semi implicit forcing terms are assumed to have the form
-    
+
        G(q, t) = H(q, t) q
-    
+
     and the semi implicit scheme will then be
-    
+
       q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
 
-    
+
     """
-    
+
     from Numeric import sum, equal, ones, Float
-        
+
     N = quantity.centroid_values.shape[0]
 
@@ -673 +673 @@
 
     for k in range(N):
-        x = quantity.centroid_values[k] 
+        x = quantity.centroid_values[k]
         if x == 0.0:
             #FIXME: Is this right
-            quantity.semi_implicit_update[k] = 0.0            
-        else:
-            quantity.semi_implicit_update[k] /= x             
-            
+            quantity.semi_implicit_update[k] = 0.0
+        else:
+            quantity.semi_implicit_update[k] /= x
+
     #Explicit updates
     quantity.centroid_values += timestep*quantity.explicit_update
-            
+
     #Semi implicit updates
     denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
@@ -702 +702 @@
         q1 = quantity.vertex_values[k, 1]
         q2 = quantity.vertex_values[k, 2]
-            
+
         quantity.edge_values[k, 0] = 0.5*(q1+q2)
-        quantity.edge_values[k, 1] = 0.5*(q0+q2)  
+        quantity.edge_values[k, 1] = 0.5*(q0+q2)
         quantity.edge_values[k, 2] = 0.5*(q0+q1)
 
 
 
-def extrapolate_second_order(quantity):        
+def extrapolate_second_order(quantity):
     """Extrapolate conserved quantities from centroid to
     vertices for each volume using
     second order scheme.
     """
-        
+
     a, b = quantity.compute_gradients()
 
@@ -720 +720 @@
     qc = quantity.centroid_values
     qv = quantity.vertex_values
-    
+
     #Check each triangle
     for k in range(quantity.domain.number_of_elements):
-        #Centroid coordinates            
+        #Centroid coordinates
         x, y = quantity.domain.centroid_coordinates[k]
-        
+
         #vertex coordinates
         x0, y0, x1, y1, x2, y2 = X[k,:]
-        
+
         #Extrapolate
         qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
         qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
-
-
-def compute_gradients(quantity):                    
+        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
+
+
+def compute_gradients(quantity):
     """Compute gradients of triangle surfaces defined by centroids of
     neighbouring volumes.
@@ -744 +744 @@
     from Numeric import zeros, Float
     from util import gradient
-    
+
     centroid_coordinates = quantity.domain.centroid_coordinates
-    surrogate_neighbours = quantity.domain.surrogate_neighbours    
-    centroid_values = quantity.centroid_values    
-    number_of_boundaries = quantity.domain.number_of_boundaries     
-    
+    surrogate_neighbours = quantity.domain.surrogate_neighbours
+    centroid_values = quantity.centroid_values
+    number_of_boundaries = quantity.domain.number_of_boundaries
+
     N = centroid_values.shape[0]
 
     a = zeros(N, Float)
     b = zeros(N, Float)
-    
+
     for k in range(N):
         if number_of_boundaries[k] < 2:
             #Two or three true neighbours
 
-            #Get indices of neighbours (or self when used as surrogate)     
+            #Get indices of neighbours (or self when used as surrogate)
             k0, k1, k2 = surrogate_neighbours[k,:]
 
-            #Get data       
+            #Get data
             q0 = centroid_values[k0]
             q1 = centroid_values[k1]
-            q2 = centroid_values[k2]                    
+            q2 = centroid_values[k2]
 
             x0, y0 = centroid_coordinates[k0] #V0 centroid
@@ -773 +773 @@
             #Gradient
             a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        
+
         elif number_of_boundaries[k] == 2:
             #One true neighbour
@@ -789 +789 @@
 
             x0, y0 = centroid_coordinates[k0] #V0 centroid
-            x1, y1 = centroid_coordinates[k1] #V1 centroid        
+            x1, y1 = centroid_coordinates[k1] #V1 centroid
 
             #Gradient
@@ -798 +798 @@
 
         else:
-            #No true neighbours -       
+            #No true neighbours -
             #Fall back to first order scheme
             pass
-        
-    
+
+
     return a, b
-        
-                    
-
-def limit(quantity):        
+
+
+
+def limit(quantity):
     """Limit slopes for each volume to eliminate artificial variance
     introduced by e.g. second order extrapolator
-    
+
     This is an unsophisticated limiter as it does not take into
     account dependencies among quantities.
-    
+
     precondition:
     vertex values are estimated from gradient
@@ -823 +823 @@
 
     N = quantity.domain.number_of_elements
-    
+
     beta_w = quantity.domain.beta_w
-        
+
     qc = quantity.centroid_values
     qv = quantity.vertex_values
-        
+
     #Find min and max of this and neighbour's centroid values
     qmax = zeros(qc.shape, Float)
-    qmin = zeros(qc.shape, Float)        
-    
+    qmin = zeros(qc.shape, Float)
+
     for k in range(N):
         qmax[k] = qmin[k] = qc[k]
@@ -839 +839 @@
             if n >= 0:
                 qn = qc[n] #Neighbour's centroid value
-                
+
                 qmin[k] = min(qmin[k], qn)
                 qmax[k] = max(qmax[k], qn)
-      
-    
+
+
     #Diffences between centroids and maxima/minima
-    dqmax = qmax - qc 
+    dqmax = qmax - qc
     dqmin = qmin - qc
-        
+
     #Deltas between vertex and centroid values
     dq = zeros(qv.shape, Float)
@@ -853 +853 @@
         dq[:,i] = qv[:,i] - qc
 
-    #Phi limiter    
+    #Phi limiter
     for k in range(N):
-        
+
         #Find the gradient limiter (phi) across vertices
         phi = 1.0
@@ -862 +862 @@
             if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
             if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
-            
-            phi = min( min(r*beta_w, 1), phi )    
+
+            phi = min( min(r*beta_w, 1), phi )
 
         #Then update using phi limiter
-        for i in range(3):            
+        for i in range(3):
             qv[k,i] = qc[k] + phi*dq[k,i]
 
@@ -874 +874 @@
 if compile.can_use_C_extension('quantity_ext.c'):
     #Replace python version with c implementations
-        
+
     from quantity_ext import limit, compute_gradients,\
-    extrapolate_second_order, interpolate_from_vertices_to_edges, update 
-
+    extrapolate_second_order, interpolate_from_vertices_to_edges, update

inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

@@ -1 +1 @@
 from visual import *
+
+elevation_color = (0.3,0.4,0.5)
+stage_color = (0.1,0.4,0.99)
 
 class Surface:
@@ -8 +11 @@
         self.frame  = frame()
 
-        self.bed_model   = faces(frame=self.frame)
-        self.stage_model = faces(frame=self.frame)
         self.domain = domain
         self.scale_z  = scale_z
-
         self.vertices = domain.vertex_coordinates
 
-        self.stage = domain.quantities['stage'].vertex_values
-
-        try:
-            self.bed   = domain.quantities['elevation'].vertex_values
-            self.xmomentum = domain.quantities['xmomentum'].vertex_values
-            self.ymomentum = domain.quantities['ymomentum'].vertex_values
-        except:
-            self.bed = None
-            self.xmomentum = None
-            self.ymomentum = None
-
-
-        self.max_x = max(max(self.vertices[:,0],self.vertices[:,2],self.vertices[:,4]))
-        self.min_x = min(min(self.vertices[:,0],self.vertices[:,2],self.vertices[:,4]))
-        self.max_y = max(max(self.vertices[:,1],self.vertices[:,3],self.vertices[:,5]))
-        self.min_y = min(min(self.vertices[:,1],self.vertices[:,3],self.vertices[:,5]))
+        # models for each quantity
+        self.z_models = {}
+        keys = self.domain.quantities.keys()
+        #print keys
+        for key in keys:
+            self.z_models[key] = faces(frame= self.frame)
+
+
+
+        self.max_x = max(max(self.vertices[:,0]),max(self.vertices[:,2]),max(self.vertices[:,4]))
+        self.min_x = min(min(self.vertices[:,0]),min(self.vertices[:,2]),min(self.vertices[:,4]))
+        self.max_y = max(max(self.vertices[:,1]),max(self.vertices[:,3]),max(self.vertices[:,5]))
+        self.min_y = min(min(self.vertices[:,1]),min(self.vertices[:,3]),min(self.vertices[:,5]))
         self.range_x = self.max_x - self.min_x
         self.range_y = self.max_y - self.min_y
-
         self.range_xy = max(self.range_x, self.range_y)
 
-        if self.bed != None and range_z == None:
-            self.max_z = max(max(self.bed))
-            self.min_z = min(min(self.bed))
-            self.range_z = max(self.max_z - self.min_z, 1.0e-10)
-
-            print 'min_bed=',self.min_z
-            print 'max_bed=',self.max_z
-            print 'range_z=',self.range_z
-        else:
-            self.max_z = range_z/2.0
-            self.min_z = -range_z/2.0
-            self.range_z = max(self.max_z - self.min_z, 1.0e-10)
+#        print 'min_x=',self.min_x
+#        print 'max_x=',self.max_x
+#        print 'range_x=',self.range_x
+#        print 'min_y=',self.min_y
+#        print 'max_y=',self.max_y
+#        print 'range_y',self.range_y
+
+        self.range_z = 1.0
+        self.max_z = self.range_z/2.0
+        self.min_z = -self.range_z/2.0
+        self.range_z = max(self.max_z - self.min_z, 1.0e-10)
 
 
@@ -68 +63 @@
         #scene.autoscale=0
         #self.update_all()
+        #scene.uniform=0
+        self.setup_range_z()
+
+    def setup_range_z(self,qname1='elevation',qname2='stage'):
+
+        #print qname1
+        #print qname2
+        self.range_z = 1.0e-10
+        try:
+            q1 = self.domain.quantities[qname1].vertex_values
+            max_z = max(max(q1))
+            min_z = min(min(q1))
+            print max_z, min_z
+            self.range_z = max(self.range_z, max_z - min_z)
+        except:
+            print 'could not find range of '+qname1
+            pass
+
+        try:
+            q2 = self.domain.quantities[qname2].vertex_values
+            max_z = max(max(q2))
+            min_z = min(min(q2))
+            print max_z, min_z
+            self.range_z = max(self.range_z, max_z - min_z)
+        except:
+            print 'Visualisation: could not find range of '+qname2
+            pass
 
     def update_all(self):
 
-        self.update_bed()
-        self.update_stage()
+        self.update_quantity('elevation',elevation_color)
+        self.update_quantity('stage',stage_color)
 
     def update_timer(self):
         self.timer.text='Time=%10.5e'%self.domain.time
 
-    def update_bed(self):
-
-        #print 'update bed arrays'
-        c0 = 0.3
-        c1 = 0.3
-        c2 = 0.3
-        if self.bed != None:
-            self.update_arrays(self.bed,  (c0,c1,c2) )
+
+    def update_quantity(self,qname,qcolor=None):
+
+        #print 'update '+qname+' arrays'
+
+        if qcolor is None:
+            if qname=='elevation':
+                qcolor = elevation_color
+
+            if qname=='stage':
+                qcolor = stage_color
+
+
+        try:
+            self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor)
 
             #print 'update bed image'
-            self.bed_model.pos    = self.pos
-            self.bed_model.color  = self.colour
-            self.bed_model.normal = self.normals
-
-
-    def update_stage(self):
-        c0 = 0.1
-        c1 = 0.4
-        c2 = 0.99
-        #print 'update stage arrays'
-        if self.stage != None:
-            self.update_arrays(self.stage, (c0,c1,c2) )
-            #print 'update stage image'
-            self.stage_model.pos    = self.pos
-            self.stage_model.color  = self.colour
-            self.stage_model.normal = self.normals
-
-    def update_xmomentum(self):
-        c0 = 0.1
-        c1 = 0.4
-        c2 = 0.99
-        #print 'update stage arrays'
-        self.update_arrays(self.xmomentum, (c0,c1,c2) )
-        #print 'update stage image'
-        self.stage_model.pos    = self.pos
-        self.stage_model.color  = self.colour
-        self.stage_model.normal = self.normals
-
-    def update_ymomentum(self):
-        c0 = 0.1
-        c1 = 0.4
-        c2 = 0.99
-        #print 'update stage arrays'
-        self.update_arrays(self.ymomentum, (c0,c1,c2) )
-        #print 'update stage image'
-        self.stage_model.pos    = self.pos
-        self.stage_model.color  = self.colour
-        self.stage_model.normal = self.normals
+            self.z_models[qname].pos    = self.pos
+            self.z_models[qname].color  = self.colour
+            self.z_models[qname].normal = self.normals
+        except:
+            print 'Visualisation: Could not update quantity '+qname
+            pass
+
+    def update_quantity_color(self,qname,qcolor=None):
+
+        #print 'update '+qname+' arrays'
+
+        if qcolor is None:
+            if qname=='elevation':
+                qcolor = elevation_color
+
+            if qname=='stage':
+                qcolor = stage_color
+
+
+        try:
+            self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor)
+
+            #print 'update bed image'
+            self.z_models[qname].pos    = self.pos
+            self.z_models[qname].color  = self.colour
+            self.z_models[qname].normal = self.normals
+        except:
+            print 'Visualisation: Could not update quantity '+qname
+            pass
+
 
     def update_arrays(self,quantity,qcolor):
@@ -151 +172 @@
 
 
+    def update_arrays_color(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
+        range_z = self.range_z
+
+        vertices = self.vertices
+        pos      = self.pos
+        normals  = self.normals
+        colour   = self.colour
+
+        try:
+            update_arrays_color_weave(vertices,quantity,col,pos,normals,colour,N,
+                                min_x,min_y,range_xy,range_z,scale_z)
+        except:
+            update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
+                                 min_x,min_y,range_xy,range_z,scale_z)
+
+
+#==================================================================================
 
 def  update_arrays_python(vertices,quantity,col,pos,normals,colour,N,
@@ -277 +324 @@
                  type_converters = converters.blitz, compiler='gcc');
 
+def  update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
+                          min_x,min_y,range_xy,range_z,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]/range_z*scale_z*0.5
+
+        v10 = v[1,:]-v[0,:]
+        v20 = v[2,:]-v[0,:]
+
+        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]
+
+        norm = sqrt( normal[0]**2 + normal[1]**2 + normal[2]**2)
+
+        normal[0] = normal[0]/norm
+        normal[1] = normal[1]/norm
+        normal[2] = normal[2]/norm
+
+        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,:]
+
+        s = abs(v(2,j))
+
+        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_color_weave(vertices,quantity,col,pos,normals,colour,
+                         N,min_x,min_y,range_xy,range_z,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)/range_z*scale_z*0.5;
+            }
+
+
+            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;
+
+            s = 0.2+fabs((v(0,2)+v(1,2)+v(2,2))/3.0);
+
+            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','range_z','scale_z','v','normal','v10','v20'],
+                 type_converters = converters.blitz, compiler='gcc');
 
 scene.width = 1000
@@ -282 +458 @@
 
 #Original
-scene.center = (0.5,0.5,0)
+scene.center = (0.5,0.5,0.0)
 scene.forward = vector(-0.1, 0.5, -0.5)
 
@@ -304 +480 @@
 
 
-def create_surface(domain,range_z=None):
-
-    surface = Surface(domain,range_z=range_z)
+def create_surface(domain,scale_z=0.5,range_z=None):
+
+    surface = Surface(domain,scale_z=0.5,range_z=range_z)
     domain.surface = surface
 
-    surface.update_bed()
-    surface.update_stage()
-    #surface.update_xmomentum()
-    #surface.update_ymomentum
+    surface.update_quantity('elevation')
+    surface.update_quantity('stage')
     surface.update_timer()
 
 def update(domain):
 
-    #print 'start visual update'
     surface = domain.surface
-    surface.update_stage()
-    #surface.update_xmomentum()
-    #surface.update_ymomentum()
+    surface.update_quantity('stage')
     surface.update_timer()
-    #print 'end visual update'

inundation/ga/storm_surge/pyvolution/view_tsh.py

@@ -40 +40 @@
 
     surface =  Surface(domain, scale_z = scale_z)
-    surface.update_bed()
-
-
+    surface.update_quantity('elevation')

inundation/ga/storm_surge/steve/run_merimbula_lake.py

@@ -29 +29 @@
 #-------
 # Domain
-filename = 'merimbula_interpolated.tsh'  
+filename = 'merimbula_interpolated.tsh'
+filename = 'merimbula_10785_1.tsh'
 print 'Creating domain from', filename
 domain = pmesh_to_domain_instance(filename, Domain)
@@ -162 +163 @@
 #   All other boundaries are reflective
 Br = Reflective_boundary(domain)
-domain.set_boundary({'external': Br, 'open': Bf})
+domain.set_boundary({'exterior': Br, 'open': Bf})
 
 #-----------
@@ -174 +175 @@
 #--------------------------------
 # Initial water surface elevation
-domain.set_quantity('stage', 1.0)
+domain.set_quantity('stage', 0.0)
     
 #----------------------------------------------------------

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

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

inundation/ga/storm_surge/zeus/steve.zpi

@@ -75 +75 @@
     <file>..\steve\get_triangle_data.py</file>
     <file>..\steve\get_triangle_data_new.py</file>
-    <file>..\steve\run_least_squares_merimbula.py</file>
     <file>..\pyvolution-parallel\vtk-pylab.py</file>
     <folder name="Header Files" />