Index: /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py
===================================================================
--- /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py	(revision 1294)
+++ /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py	(revision 1295)
@@ -4,5 +4,5 @@
    Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
    Geoscience Australia, ANU
-   
+
 Specific methods pertaining to the 2D shallow water equation
 are imported from shallow_water
@@ -25,5 +25,5 @@
 from pmesh2domain import pmesh_to_domain_instance
 from util import file_function, Polygon_function, read_polygon
-from Numeric import zeros, Float 
+from Numeric import zeros, Float
 
 #-------
@@ -36,7 +36,7 @@
 
 #------------------------------------------
-# Reduction operation for get_vertex_values              
+# Reduction operation for get_vertex_values
 from util import mean
-domain.reduction = mean 
+domain.reduction = mean
 
 #-----------------
@@ -176,5 +176,5 @@
 # Initial water surface elevation
 domain.set_quantity('stage', 0.0)
-    
+
 #----------------------------------------------------------
 # Decide which quantities are to be stored at each timestep
@@ -191,6 +191,8 @@
 domain.default_order = 1
 domain.visualise = True
+domain.visualise_color_stage = True
+domain.visualise_stage_range = 1.0
 domain.smooth = True
-         
+
 #---------
 #Evolution
@@ -199,6 +201,6 @@
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
-    
+
 print 'Done'
 
-    
+
Index: /inundation/ga/storm_surge/pyvolution/netherlands.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/netherlands.py	(revision 1294)
+++ /inundation/ga/storm_surge/pyvolution/netherlands.py	(revision 1295)
@@ -63,4 +63,8 @@
                 z[i] = -x[i]/20+0.4
 
+            if (x[i] - 0.72)**2 + (y[i] - 0.4)**2 < 0.05**2:# or\
+                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
+                z[i] = -x[i]/20+0.4
+
             #Wall
             if x[i] > 0.995:
@@ -79,5 +83,5 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 30
+N = 40
 
 #N = 15
@@ -100,4 +104,8 @@
 print "Number of triangles = ", len(domain)
 
+#Set bed-slope and friction
+inflow_stage = 0.4
+manning = 0.02
+Z = Weir(inflow_stage)
 
 if N > 150:
@@ -112,12 +120,9 @@
 else:
     domain.visualise = True
+    domain.visualise_color_stage = True
+    domain.visualise_timer = True
     domain.checkpoint = False
     domain.store = False
 
-
-#Set bed-slope and friction
-inflow_stage = 0.1
-manning = 0.02
-Z = Weir(inflow_stage)
 
 print 'Field values'
@@ -134,5 +139,5 @@
 
 #Constant inflow
-Bd = Dirichlet_boundary([2*inflow_stage, 0.0, 0.0])
+Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
 
 
@@ -151,5 +156,5 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.2, finaltime = 10.0):
+for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
     domain.write_time()
 
Index: /inundation/ga/storm_surge/pyvolution/realtime_visualisation.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/realtime_visualisation.py	(revision 1294)
+++ /inundation/ga/storm_surge/pyvolution/realtime_visualisation.py	(revision 1295)
@@ -11,5 +11,5 @@
         #n = vector(0,0,0)
 
-    return n    
+    return n
 
 
@@ -20,16 +20,16 @@
 
         #Outward normal
-        
+
         self.normal=normal(v1-v0, v2-v0)
 
         self.panel = faces(frame = frame)
-        self.backpanel = faces(frame = frame)        
+        self.backpanel = faces(frame = frame)
         for v in (v0,v1,v2):
             self.panel.append(pos=v, normal = self.normal, color=color)
 
         for v in (v2,v1,v0):
-            self.backpanel.append(pos=v, normal = -self.normal, color=color)            
-
-            
+            self.backpanel.append(pos=v, normal = -self.normal, color=color)
+
+
     def move(self,v):
         """Move panel in direction given by vector v
@@ -37,5 +37,5 @@
 
         for vertex in self.panel.pos:
-            vertex += v            
+            vertex += v
 
     def update_height(self,d,fix_baseline=False):
@@ -47,6 +47,6 @@
         v0 = self.panel.pos[0]
         v1 = self.panel.pos[1]
-        v2 = self.panel.pos[2]        
-        
+        v2 = self.panel.pos[2]
+
         n = normal(v1-v0, self.normal)
 
@@ -54,5 +54,5 @@
             self.panel.pos[2] -= d*n
         else:
-            self.panel.pos[:2] += d*n        
+            self.panel.pos[:2] += d*n
 
     def update_color(self, c):
@@ -61,5 +61,5 @@
 
         self.panel.color = c
-        
+
 
     def set_vertexheights(self, heights, floor_heights = None):
@@ -70,5 +70,5 @@
             floor_heights = zeros(heights.shape, Float)
 
-        all_vertices_below_threshold = True    
+        all_vertices_below_threshold = True
         for k in range(3):
             w = heights[k]
@@ -77,5 +77,5 @@
             if w-z >= hmin:
                 all_vertices_below_threshold = False
-            
+
             vertex = self.panel.pos[k]
             vertex[2] = w
@@ -85,8 +85,8 @@
 
 
-        #Update color to visualise dry areas 
-        if all_vertices_below_threshold: 
+        #Update color to visualise dry areas
+        if all_vertices_below_threshold:
             self.panel.color = self.bottom_color
-            self.backpanel.color = self.bottom_color                
+            self.backpanel.color = self.bottom_color
         else:
             self.panel.color = self.top_color
@@ -94,5 +94,5 @@
 
 
-        #update normal    
+        #update normal
         v0 = self.panel.pos[0]
         v1 = self.panel.pos[1]
@@ -104,5 +104,5 @@
         self.backpanel.normal=-n
 
-        
+
 
 
@@ -120,5 +120,5 @@
     except:
         bed = False
-        
+
     N = Q.vertex_values.shape[0]
 
@@ -128,10 +128,10 @@
         z0 = Q.vertex_values[i, 0]
         z1 = Q.vertex_values[i, 1]
-        z2 = Q.vertex_values[i, 2]        
+        z2 = Q.vertex_values[i, 2]
 
         x0, y0 = domain.get_vertex_coordinate(i,0)
         x1, y1 = domain.get_vertex_coordinate(i,1)
-        x2, y2 = domain.get_vertex_coordinate(i,2)        
-        
+        x2, y2 = domain.get_vertex_coordinate(i,2)
+
         V0 = vector(x0, y0, z0)
         V1 = vector(x1, y1, z1)
@@ -149,9 +149,9 @@
 
         col = (c0,c1,c2)
-            
+
         visual_top = Triangle(V0,V1,V2,color=col,frame=fr)
         visual_top.top_color = col
 
-    
+
         #Bottom surface
         #v0, v1, v2 = volume.vertices
@@ -168,5 +168,5 @@
         V0 = vector(x0, y0, z0)
         V1 = vector(x1, y1, z1)
-        V2 = vector(x2, y2, z2)        
+        V2 = vector(x2, y2, z2)
 
         c0 = 0.3
@@ -176,17 +176,17 @@
             c2 = 0.4*c2   #To show triangles in slightly different shades
 
-        col = (c0,c1,c2)   
+        col = (c0,c1,c2)
         visual_bottom = Triangle(V0, V1, V2, color=col,frame=fr)
         visual_top.bottom_color=col
-        
+
         domain.visuals.append( (visual_top, visual_bottom) )
-        
+
     update(domain)
-    #print 'Scale', scene.scale    
-
-
-def update(domain):    
+    #print 'Scale', scene.scale
+
+
+def update(domain):
     """Update vertex heights.
-    The argument index refers to which conserved quantity to visualise.    
+    The argument index refers to which conserved quantity to visualise.
     If domain.smooth is set True, vertex values will be averaged
     yielding a smoother surface.
@@ -198,5 +198,5 @@
     Q = domain.quantities['stage']
     N = Q.vertex_values.shape[0]
-    
+
     #print scene.forward
     #FIXME: Use smoother from pyvolution instead
@@ -218,13 +218,13 @@
 
 
-    for k in range(N):    
+    for k in range(N):
         if domain.smooth:
             #The averages
             x = zeros(3, Float)
-            for i in range(3):            
+            for i in range(3):
                 vertex = domain.triangles[k, i]
                 A = array(vertex_heights[vertex])
                 x[i] = sum(A)/len(A)
-        else:        
+        else:
             #The true values
             x = [Q.vertex_values[k, 0],
@@ -232,8 +232,8 @@
                  Q.vertex_values[k, 2]]
 
-            
+
         #Do it
         floor_heights = array([pos[2] for pos in domain.visuals[k][1].panel.pos])
-        
+
         domain.visuals[k][0].set_vertexheights(x, floor_heights)
 
@@ -246,8 +246,8 @@
 #Original
 scene.center = (0.5,0.5,0)
-scene.forward = vector(-0.1, 0.5, -0.5)
+scene.forward = vector(0.0, 0.5, -0.5)
 
 #Temporary (for bedslope)
-#scene.forward = vector(0.0006, 0.7, -0.03)    
+#scene.forward = vector(0.0006, 0.7, -0.03)
 
 
@@ -258,5 +258,5 @@
 #Temporary for hackett - end
 
-    
+
 scene.ambient = 0.4
 scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
@@ -266,3 +266,2 @@
 
 
-
Index: /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py	(revision 1294)
+++ /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py	(revision 1295)
@@ -59,5 +59,5 @@
         self.border_model = curve(frame = self.frame,
                                   pos=[(0,0),(0,1),(1,1),(1,0),(0,0)])
-        self.timer=label(pos=(0.75,0.5,0.5),text='Time=%10.5e'%self.domain.time)
+        self.timer=label(pos=(0.75,0.75,0.5),text='Time=%10.5e'%self.domain.time)
 
         #scene.autoscale=0
@@ -135,14 +135,14 @@
 
 
-        try:
-            self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor)
+        #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
+        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
 
 
@@ -189,10 +189,10 @@
         colour   = self.colour
 
-        try:
-            update_arrays_color_weave(vertices,quantity,col,pos,normals,colour,N,
+        #try:
+        update_arrays_color_python(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)
+        #except:
+        #    update_arrays_color_python(vertices,quantity,col,pos,normals,colour,N,
+        #                         min_x,min_y,range_xy,range_z,scale_z)
 
 
@@ -359,12 +359,22 @@
         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
+        q0 = quantity[i,0]/0.4
+        q1 = quantity[i,1]/0.4
+        q2 = quantity[i,2]/0.4
+
+        q0r = min(q0,0.0)
+        q1r = min(q1,0.0)
+        q2r = min(q2,0.0)
+        q0b = max(q0,0.0)
+        q1b = max(q1,0.0)
+        q2b = max(q2,0.0)
+
+
+        colour[6*i  ,:] = s*array([0.3 - q0r, 0.3 - q0, 0.3 + q0b])
+        colour[6*i+1,:] = s*array([0.3 - q1r, 0.3 - q1, 0.3 + q1b])
+        colour[6*i+2,:] = s*array([0.3 - q2r, 0.3 - q2, 0.3 + q2b])
+        colour[6*i+3,:] = s*array([0.3 - q0r, 0.3 - q0, 0.3 + q0b])
+        colour[6*i+4,:] = s*array([0.3 - q2r, 0.3 - q2, 0.3 + q2b])
+        colour[6*i+5,:] = s*array([0.3 - q1r, 0.3 - q1, 0.3 + q1b])
 
         s =  15.0/8.0 - s
@@ -417,5 +427,5 @@
             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++) {
@@ -427,14 +437,4 @@
                 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);
@@ -445,4 +445,13 @@
                 }
 
+            s = 0.2+fabs((v(0,2)+v(1,2)+v(2,2))/3.0);
+            for (int j=0; j<3; 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);
+                }
             s =  15.0/8.0 - s;
         }
@@ -459,5 +468,5 @@
 #Original
 scene.center = (0.5,0.5,0.0)
-scene.forward = vector(-0.1, 0.5, -0.5)
+scene.forward = vector(-0.0, 0.5, -0.8)
 
 #Temporary (for bedslope)
@@ -492,4 +501,9 @@
 
     surface = domain.surface
-    surface.update_quantity('stage')
-    surface.update_timer()
+    if domain.visualise_color_stage:
+        surface.update_quantity_color('stage')
+    else:
+        surface.update_quantity('stage')
+
+    if domain.visualise_timer:
+        surface.update_timer()
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1294)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1295)
@@ -80,4 +80,7 @@
         #Realtime visualisation
         self.visualise = False
+        self.visualise_color_stage = False
+        self.visualise_stage_range = 1.0
+        self.visualise_timer = False
 
         #Stored output
