Index: /anuga_work/development/anuga_1d/analytic_dam_sudi.py
===================================================================
--- /anuga_work/development/anuga_1d/analytic_dam_sudi.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/analytic_dam_sudi.py	(revision 5827)
@@ -26,5 +26,5 @@
             else:
                 zmax=z
-
+            #print 'func=',func
         if( abs(z) > 99.0):
             print 'no convergence'
Index: /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
===================================================================
--- /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5826)
+++ /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5827)
@@ -26,9 +26,9 @@
 //Innermost flux function (using w=z+h)
 int _flux_function_wellbalanced(double *q_left, double *q_right,
-        double z_left, double z_right,
 		double normals, double g, double epsilon,
 		double *edgeflux, double *max_speed) {
 		
 		int i;
+        double z_left, z_right;
 		double ql[2], qr[2], flux_left[2], flux_right[2];
 		double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
@@ -36,21 +36,23 @@
 		double s_max, s_min, denom;
 		
+		w_left = q_left[0];
+		uh_left = q_left[1];
+		uh_left = uh_left*normals;
+        z_left = q_left[2];
+	
+		w_right = q_right[0];
+		uh_right = q_right[1];
+		uh_right = uh_right*normals;
+        z_right = q_right[2];
+	  
 		if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
-		ql[0] = q_left[0];
-		ql[1] = q_left[1];
-		ql[1] = ql[1]*normals;
-	
-		qr[0] = q_right[0];
-		qr[1] = q_right[1];
-		qr[1] = qr[1]*normals;
-	  
+        
 		//z = (z_left+z_right)/2.0;
 		z_star = max(z_left, z_right);		              											//equation(7)  
 				   
 		// Compute speeds in x-direction
-		w_left = ql[0];
 		h_left = w_left-z_left;		                                                                //This is used for computing the edgeflux[1].
 		h_left_star = max(0, w_left-z_star);														//(8)
-		uh_left = ql[1];
+    
 		if (h_left_star < epsilon) {
 			u_left = 0.0; h_left_star = 0.0;
@@ -60,8 +62,7 @@
 		}
 	
-		w_right = qr[0];
 		h_right = w_right-z_right;                                                                                                                                                                                                     
 		h_right_star = max(0, w_right-z_star);														//(9)
-		uh_right = qr[1];
+
 		if (h_right_star < epsilon) {
 			u_right = 0.0; h_right_star = 0.0; 
@@ -134,5 +135,5 @@
 		double* max_speed_array) {
 		
-		double flux[2], ql[2], qr[2], edgeflux[2];
+		double flux[2], ql[3], qr[3], edgeflux[2];
 		double zl, zr, max_speed, normal;
 		int k, i, ki, n, m, nm=0;
@@ -147,6 +148,8 @@
 				ql[0] = stage_edge_values[ki];
 				ql[1] = xmom_edge_values[ki];
-				zl = bed_edge_values[ki];
-				
+				ql[2] = bed_edge_values[ki];
+				//ql[3] = velocity_edge_values[ki];
+                //ql[4] = height_edge_values[ki];
+                
 				n = neighbours[ki];
 				if (n<0) {
@@ -154,5 +157,6 @@
 					qr[0] = stage_boundary_values[m];
 					qr[1] = xmom_boundary_values[m];
-					zr = zl;
+					qr[2] = ql[2];
+                    
 				} else {
 					m = neighbour_vertices[ki];
@@ -160,9 +164,9 @@
 					qr[0] = stage_edge_values[nm];
 					qr[1] = xmom_edge_values[nm];
-					zr = bed_edge_values[nm];				
+					qr[2] = bed_edge_values[nm];				
 				}
 				
 				normal = normals[ki];
-				_flux_function_wellbalanced(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
+				_flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
 				flux[0] -= edgeflux[0];
 				flux[1] -= edgeflux[1];
@@ -186,5 +190,5 @@
 			max_speed_array[k]=max_speed;
 		}
-		printf("timestep = %f \n",timestep);
+		//printf("timestep = %f \n",timestep);
 		return timestep;	
 	}
Index: /anuga_work/development/anuga_1d/config.py
===================================================================
--- /anuga_work/development/anuga_1d/config.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/config.py	(revision 5827)
@@ -3,5 +3,5 @@
 
 epsilon = 1.0e-12
-h0 = 1.0e-6
+h0 = 1.0e-12
 
 default_boundary_tag = 'exterior'
@@ -39,5 +39,5 @@
 max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order
 
-manning = 0.3  #Manning's friction coefficient
+manning = 0.0  #Manning's friction coefficient
 g = 9.8       #Gravity
 #g(phi) = 9780313 * (1 + 0.0053024 sin(phi)**2 - 0.000 0059 sin(2*phi)**2) micro m/s**2, where phi is the latitude
@@ -109,4 +109,4 @@
 
 #Specific to shallow water W.E.
-minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0
-maximum_allowed_speed = 100.0
+minimum_allowed_height = 1.0e-6 #Water depth below which it is considered to be 0
+maximum_allowed_speed = 1000.0
Index: /anuga_work/development/anuga_1d/dam_h_elevation.py
===================================================================
--- /anuga_work/development/anuga_1d/dam_h_elevation.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/dam_h_elevation.py	(revision 5827)
@@ -31,5 +31,5 @@
 
 L=2000.0
-N=200
+N=50
 
 cell_len=L/N
@@ -123,6 +123,6 @@
 import time
 t0=time.time()
-yieldstep=45.0 #30.0
-finaltime=45.0 #20.0
+yieldstep=10.0 #30.0
+finaltime=10.0 #20.0
 print "integral", domain.quantities['stage'].get_integral()
 for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
@@ -152,5 +152,5 @@
 
         plot(X,ElevationQ,X,HeightQ)
-        #plot1.set_ylim([9.99,10.01])
+        plot1.set_ylim([-1.0,11.0])
         xlabel('Position')
         ylabel('Stage')
Index: /anuga_work/development/anuga_1d/dry_dam_sudi.py
===================================================================
--- /anuga_work/development/anuga_1d/dry_dam_sudi.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/dry_dam_sudi.py	(revision 5827)
@@ -1,5 +1,5 @@
 import os
 from math import sqrt, pi
-from shallow_water_domain import *
+from shallow_water_domain_suggestion1 import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon
@@ -8,5 +8,5 @@
     
     #t  = 0.0     # time (s)
-    g  = 9.81    # gravity (m/s^2)
+    # gravity (m/s^2)
     h1 = 10.0    # depth upstream (m)
     h0 = 0.0     # depth downstream (m)
@@ -43,5 +43,5 @@
             h[i] = h0 
 
-    return h , u*h
+    return h , u*h, u
 
 #def newLinePlot(title='Simple Plot'):
@@ -65,14 +65,4 @@
 print "TEST 1D-SOLUTION III -- DRY BED"
 
-L = 2000.0     # Length of channel (m)
-N = 800        # Number of computational cells
-cell_len = L/N # Origin = 0.0
-
-points = zeros(N+1,Float)
-for i in range(N+1):
-    points[i] = i*cell_len
-
-domain = Domain(points)
-
 def stage(x):
     h1 = 10.0
@@ -92,7 +82,7 @@
 
 finaltime = 20.0
-yieldstep = 1.0
+yieldstep = 20.0
 L = 2000.0     # Length of channel (m)
-number_of_cells = [200]#,200,500,1000,2000,5000,10000,20000]
+number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
 h_error = zeros(len(number_of_cells),Float)
 uh_error = zeros(len(number_of_cells),Float)
@@ -110,8 +100,8 @@
     domain.set_quantity('stage', stage)
     domain.set_boundary({'exterior': Reflective_boundary(domain)})
-    domain.default_order = 2
-    domain.default_time_order = 2
+    domain.order = 2
+    domain.set_timestepping_method('rk2')
     domain.cfl = 1.0
-    domain.limiter = "minmod"
+    domain.limiter = "vanleer"
 
     t0 = time.time()
@@ -124,5 +114,5 @@
     XmomC = domain.quantities['xmomentum'].centroid_values
     C = domain.centroids
-    h, uh = analytical_sol(C,domain.time)
+    h, uh, u = analytical_sol(C,domain.time)
     h_error[k] = 1.0/(N)*sum(abs(h-StageC))
     uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
@@ -134,5 +124,7 @@
     StageQ = domain.quantities['stage'].vertex_values
     XmomQ = domain.quantities['xmomentum'].vertex_values
-    h, uh = analytical_sol(X.flat,domain.time)
+    velQ = domain.quantities['velocity'].vertex_values
+    
+    h, uh, u = analytical_sol(X.flat,domain.time)
     x = X.flat
     
@@ -151,9 +143,9 @@
            'upper right', shadow=True)
     plot2 = subplot(212)
-    plot(x,uh,x,XmomQ.flat)
+    plot(x,u,x,velQ.flat)
     plot2.set_ylim([-35,35])
     
     xlabel('Position')
-    ylabel('Xmomentum')
+    ylabel('Velocity')
    
     file = "dry_bed_"
Index: /anuga_work/development/anuga_1d/parabolic_cannal.py
===================================================================
--- /anuga_work/development/anuga_1d/parabolic_cannal.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/parabolic_cannal.py	(revision 5827)
@@ -1,8 +1,8 @@
 import os
 from math import sqrt, pi, sin, cos
-from shallow_water_1d import *
+from shallow_water_domain import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon
-
+"""
 def newLinePlot(title='Simple Plot'):
     import Gnuplot
@@ -20,5 +20,5 @@
     plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1")
     g.plot(plot1,plot2,plot3)
-
+"""
 
 def analytic_cannal(C,t):
@@ -54,5 +54,5 @@
 
 L_x = 2500.0     # Length of channel (m)
-N = 8         # Number of compuational cells
+N = 100         # Number of compuational cells
 cell_len = 4*L_x/N   # Origin = 0.0
 
@@ -63,7 +63,6 @@
 domain = Domain(points)
 
-domain.order = 2 #make this unnecessary
-domain.default_order = 2
-domain.default_time_order = 2
+domain.order = 2
+domain.set_timestepping_method('rk2')
 domain.cfl = 1.0
 domain.beta = 1.0
@@ -72,5 +71,6 @@
 #domain.limiter = "vanalbada"
 #domain.limiter = "superbee"
-domain.limiter = "steve_minmod"
+#domain.limiter = "steve_minmod"
+domain.limiter = "pyvolution"
 
 def stage(x):
@@ -113,18 +113,18 @@
 domain.set_quantity('elevation',elevation)
 #domain.set_quantity('height',height)
-
 domain.set_boundary({'exterior': Reflective_boundary(domain)})
  
 import time
 t0 = time.time()
-finaltime = 1122.0*0.75
-yieldstep = finaltime
-yieldstep = 10.0
-plot1 = newLinePlot("Stage")
-plot2 = newLinePlot("Xmomentum")
+#finaltime = 1122.0*0.75
+yieldstep = finaltime = 10.0
+#yieldstep = 10.0
+#plot1 = newLinePlot("Stage")
+#plot2 = newLinePlot("Xmomentum")
 C = domain.centroids
 X = domain.vertices
 StageQ = domain.quantities['stage'].vertex_values
 XmomQ = domain.quantities['xmomentum'].vertex_values
+#Bed = domain.quantities['elevation'].vertex_values
 import time
 t0 = time.time()
@@ -133,6 +133,6 @@
     domain.write_time()
     u,h,z,z_b = analytic_cannal(X.flat,t)
-    linePlot(plot1,X,z,X,z_b,X,StageQ)
-    linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
+    #linePlot(plot1,X,z,X,z_b,X,StageQ)
+    #linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
     HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
     u,hc,z,z_b = analytic_cannal(C,t)
@@ -142,27 +142,26 @@
     error = 1.0/(N)*sum(abs(hc-HeightQ))
     print 'Error measured at centroids', error
-    #raw_input()         
-
-#from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-#import Gnuplot  
-X = domain.vertices
-u,h,z,z_b = analytic_cannal(X.flat,domain.time)
-StageQ = domain.quantities['stage'].vertex_values
-XmomQ = domain.quantities['xmomentum'].vertex_values
-hold(False)
-plot1 = subplot(211)
-plot(X,h,X,StageQ)
-#plot1.set_ylim([4,11])
-#title('Free Surface Elevation of a Dry Dam-Break')
-xlabel('Position')
-ylabel('Stage')
-legend(('Analytical Solution', 'Numerical Solution'),
-       'upper right', shadow=True)
-plot2 = subplot(212)
-plot(X,u*h,X,XmomQ)
-#plot2.set_ylim([-1,25])
-#title('Xmomentum Profile of a Dry Dam-Break')
-xlabel('Position')
-ylabel('Xmomentum')
-show()
+    # #raw_input()         
+    # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    # #import Gnuplot  
+    X = domain.vertices
+    u,h,z,z_b = analytic_cannal(X.flat,domain.time)
+    StageQ = domain.quantities['stage'].vertex_values
+    XmomQ = domain.quantities['xmomentum'].vertex_values
+    hold(False)
+    plot1 = subplot(211)
+    plot(X,h,X,StageQ,X,z_b)
+    #plot1.set_ylim([4,11])
+    #title('?????????')
+    xlabel('Position')
+    ylabel('Stage')
+    legend(('Analytical Solution', 'Numerical Solution', 'Bed'),
+           'upper right', shadow=True)
+    plot2 = subplot(212)
+    plot(X,u*h,X,XmomQ)
+    #plot2.set_ylim([-1,25])
+    #title('?????????????????????????')
+    xlabel('Position')
+    ylabel('Xmomentum')
+    show()
Index: /anuga_work/development/anuga_1d/parabolic_cannal_sudi.py
===================================================================
--- /anuga_work/development/anuga_1d/parabolic_cannal_sudi.py	(revision 5827)
+++ /anuga_work/development/anuga_1d/parabolic_cannal_sudi.py	(revision 5827)
@@ -0,0 +1,126 @@
+import os
+from math import sqrt, pi, sin, cos
+from shallow_water_domain_suggestion1 import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+def analytic_cannal(C,t):
+    N = len(C)
+    u = zeros(N,Float)    ## water velocity
+    h = zeros(N,Float)    ## water depth
+    x = C
+    g = 9.81
+    ## Define Basin Bathymetry
+    z_b = zeros(N,Float) ## elevation of basin
+    z = zeros(N,Float)   ## elevation of water surface
+    z_infty = 10.0       ## max equilibrium water depth at lowest point.
+    L_x = 2500.0         ## width of channel
+    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+    for i in range(N):
+        z_b[i] = z_infty*(x[i]**2/L_x**2)
+        u[i] = -A0*omega*sin(omega*t)
+        z[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
+    h = z-z_b
+    T = 2.0*pi/omega
+    return u,h,z,z_b, T
+
+
+L_x = 2500.0     # Length of channel (m)
+N = 400         # Number of compuational cells
+cell_len = 4*L_x/N   # Origin = 0.0
+
+points = zeros(N+1,Float)
+for i in range(N+1):
+        points[i] = -2*L_x +i*cell_len
+
+domain = Domain(points)
+domain.order = 2
+domain.set_timestepping_method('rk2')
+domain.cfl = 1.0
+domain.beta = 1.0
+domain.limiter = "vanleer"
+
+def stage(x):
+    z_infty = 10.0       ## max equilibrium water depth at lowest point.
+    L_x = 2500.0         ## width of channel
+    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+    t=0.0
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
+        #y[i] = 12.0
+    return y
+
+def elevation(x):
+    N = len(x)
+    z_infty = 10.0
+    z = zeros(N,Float)
+    L_x = 2500.0
+    A0 = 0.5*L_x
+    omega = sqrt(2*g*z_infty)/L_x
+    for i in range(N):
+        z[i] = z_infty*(x[i]**2/L_x**2)
+    return z
+
+def height(x):
+    z_infty = 10.0       ## max equilibrium water depth at lowest point.
+    L_x = 2500.0         ## width of channel
+    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+    t=0.0
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i] = max(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))-z_infty*(x[i]**2/L_x**2),0.0)
+    return y
+
+    
+domain.set_quantity('stage', stage)
+domain.set_quantity('elevation',elevation)
+domain.set_boundary({'exterior': Reflective_boundary(domain)})
+ 
+X = domain.vertices
+u,h,z,z_b,T = analytic_cannal(X.flat,domain.time)
+print 'T = ',T
+
+yieldstep = finaltime = T/2.0
+StageQ = domain.quantities['stage'].vertex_values
+XmomQ = domain.quantities['xmomentum'].vertex_values
+
+import time
+t0 = time.time()
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+    print "integral", domain.quantities['stage'].get_integral()
+    if t>0.0:
+        
+        print 'X=',X
+        print 'domain.time=',domain.time
+        
+        StageQ = domain.quantities['stage'].vertex_values
+        XmomQ = domain.quantities['xmomentum'].vertex_values
+        
+        u,h,z,z_b,T = analytic_cannal(X.flat,domain.time)
+        w = z
+        for k in range(len(h)):
+            if h[k] < 0.0:
+                h[k] = 0.0
+                w[k] = z_b[k]
+        
+        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+        hold(False)
+        plot1 = subplot(211)
+        plot(X,w,  X,StageQ,  X,z_b)
+        #plot1.set_ylim([4,11])
+        xlabel('Position')
+        ylabel('Stage')
+        legend(('Analytic Solution', 'Numerical Solution', 'Bed'),
+               'upper center', shadow=True)
+        plot2 = subplot(212)
+        plot(X,u*h,X,XmomQ)
+        #plot2.set_ylim([-1,25])
+        xlabel('Position')
+        ylabel('Xmomentum')
+        show()
+raw_input("Press the return key!")
Index: /anuga_work/development/anuga_1d/quantity.py
===================================================================
--- /anuga_work/development/anuga_1d/quantity.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/quantity.py	(revision 5827)
@@ -764,6 +764,6 @@
             elif n1 < 0:
                 phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
-            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
-                phi = 0.0
+            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
+            #    phi = 0.0
             else:
                 if limiter == "minmod":
Index: /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py	(revision 5827)
+++ /anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py	(revision 5827)
@@ -0,0 +1,1450 @@
+"""Class Domain -
+1D interval domains for finite-volume computations of
+the shallow water wave equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the Shallow Water Wave Equation
+
+
+U_t + E_x = S
+
+where
+
+U = [w, uh]
+E = [uh, u^2h + gh^2/2]
+S represents source terms forcing the system
+(e.g. gravity, friction, wind stress, ...)
+
+and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
+
+The quantities are
+
+symbol    variable name    explanation
+x         x                horizontal distance from origin [m]
+z         elevation        elevation of bed on which flow is modelled [m]
+h         height           water height above z [m]
+w         stage            absolute water level, w = z+h [m]
+u                          speed in the x direction [m/s]
+uh        xmomentum        momentum in the x direction [m^2/s]
+
+eta                        mannings friction coefficient [to appear]
+nu                         wind stress coefficient [to appear]
+
+The conserved quantities are w, uh
+
+For details see e.g.
+Christopher Zoppou and Stephen Roberts,
+Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
+Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
+
+
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2006
+"""
+
+
+from domain import *
+Generic_Domain = Domain #Rename
+
+#Shallow water domain
+class Domain(Generic_Domain):
+
+    def __init__(self, coordinates, boundary = None, tagged_elements = None):
+
+        conserved_quantities = ['stage', 'xmomentum']
+        other_quantities = ['elevation', 'friction', 'height', 'velocity']
+        Generic_Domain.__init__(self, coordinates, boundary,
+                                conserved_quantities, other_quantities,
+                                tagged_elements)
+        
+        from config import minimum_allowed_height, g, h0
+        self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
+        self.h0 = h0
+
+        #forcing terms not included in 1d domain ?WHy?
+        self.forcing_terms.append(gravity)
+        #self.forcing_terms.append(manning_friction)
+        #print "\nI have Removed forcing terms line 64 1dsw"
+
+        #Realtime visualisation
+        self.visualiser = None
+        self.visualise  = False
+        self.visualise_color_stage = False
+        self.visualise_stage_range = 1.0
+        self.visualise_timer = True
+        self.visualise_range_z = None
+        
+        #Stored output
+        self.store = True
+        self.format = 'sww'
+        self.smooth = True
+
+        #Evolve parametrs
+        self.cfl = 1.0
+        
+        #Reduction operation for get_vertex_values
+        from util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.quantities_to_be_stored = ['stage','xmomentum']
+
+	self.__doc__ = 'shallow_water_domain'
+
+
+    def set_quantities_to_be_stored(self, q):
+        """Specify which quantities will be stored in the sww file.
+
+        q must be either:
+          - the name of a quantity
+          - a list of quantity names
+          - None
+
+        In the two first cases, the named quantities will be stored at each
+        yieldstep
+        (This is in addition to the quantities elevation and friction)  
+        If q is None, storage will be switched off altogether.
+        """
+
+
+        if q is None:
+            self.quantities_to_be_stored = []            
+            self.store = False
+            return
+
+        if isinstance(q, basestring):
+            q = [q] # Turn argument into a list
+
+        #Check correcness    
+        for quantity_name in q:
+            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
+            assert quantity_name in self.conserved_quantities, msg 
+        
+        self.quantities_to_be_stored = q
+        
+
+    def initialise_visualiser(self,scale_z=1.0,rect=None):
+        #Realtime visualisation
+        if self.visualiser is None:
+            from realtime_visualisation_new import Visualiser
+            self.visualiser = Visualiser(self,scale_z,rect)
+            self.visualiser.setup['elevation']=True
+            self.visualiser.updating['stage']=True
+        self.visualise = True
+        if self.visualise_color_stage == True:
+            self.visualiser.coloring['stage'] = True
+            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
+        print 'initialise visualiser'
+        print self.visualiser.setup
+        print self.visualiser.updating
+
+    def check_integrity(self):
+        Generic_Domain.check_integrity(self)
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "stage"'
+        assert self.conserved_quantities[0] == 'stage', msg
+        msg = 'Second conserved quantity must be "xmomentum"'
+        assert self.conserved_quantities[1] == 'xmomentum', msg
+
+    def extrapolate_second_order_sw(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_sw(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_C_wellbalanced(self)
+
+    def compute_timestep(self):
+        #Call correct module function
+        compute_timestep(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges(self)
+        
+    def evolve(self, yieldstep = None, finaltime = None, duration = None,
+               skip_initial_step = False):
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call check integrity here rather than from user scripts
+        #self.check_integrity()
+
+        #msg = 'Parameter beta_h must be in the interval [0, 1)'
+        #assert 0 <= self.beta_h < 1.0, msg
+        #msg = 'Parameter beta_w must be in the interval [0, 1)'
+        #assert 0 <= self.beta_w < 1.0, msg
+
+
+        #Initial update of vertex and edge values before any storage
+        #and or visualisation
+        
+        #self.distribute_to_vertices_and_edges ?????????????????????????????????
+        
+        #Initialise real time viz if requested
+        #if self.visualise is True and self.time == 0.0:
+        #    if self.visualiser is None:
+        #        self.initialise_visualiser()
+        #
+        #    self.visualiser.update_timer()
+        #    self.visualiser.setup_all()
+
+        #Store model data, e.g. for visualisation
+        #if self.store is True and self.time == 0.0:
+        #    self.initialise_storage()
+        #    #print 'Storing results in ' + self.writer.filename
+        #else:
+        #    pass
+        #    #print 'Results will not be stored.'
+        #    #print 'To store results set domain.store = True'
+        #    #FIXME: Diagnostic output should be controlled by
+        #    # a 'verbose' flag living in domain (or in a parent class)
+
+        #Call basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
+                                       skip_initial_step):
+            #Real time viz
+        #    if self.visualise is True:
+        #        self.visualiser.update_all()
+        #        self.visualiser.update_timer()
+
+
+            #Store model data, e.g. for subsequent visualisation
+        #    if self.store is True:
+        #        self.store_timestep(self.quantities_to_be_stored)
+
+            #FIXME: Could maybe be taken from specified list
+            #of 'store every step' quantities
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+    def initialise_storage(self):
+        """Create and initialise self.writer object for storing data.
+        Also, save x and bed elevation
+        """
+
+        import data_manager
+
+        #Initialise writer
+        self.writer = data_manager.get_dataobject(self, mode = 'w')
+
+        #Store vertices and connectivity
+        self.writer.store_connectivity()
+
+
+    def store_timestep(self, name):
+        """Store named quantity and time.
+
+        Precondition:
+           self.write has been initialised
+        """
+        self.writer.store_timestep(name)
+
+
+#=============== End of Shallow Water Domain ===============================
+
+#Rotation of momentum vector
+def rotate(q, normal, direction = 1):
+    """Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+
+    If direction is negative the rotation is inverted.
+
+    Input vector is preserved
+
+    This function is specific to the shallow water wave equation
+    """
+
+    from Numeric import zeros, Float
+
+    assert len(q) == 3,\
+           'Vector of conserved quantities must have length 3'\
+           'for 2D shallow water equation'
+
+    try:
+        l = len(normal)
+    except:
+        raise 'Normal vector must be an Numeric array'
+
+    assert l == 2, 'Normal vector must have 2 components'
+
+
+    n1 = normal[0]
+    n2 = normal[1]
+
+    r = zeros(len(q), Float) #Rotated quantities
+    r[0] = q[0]              #First quantity, height, is not rotated
+
+    if direction == -1:
+        n2 = -n2
+
+
+    r[1] =  n1*q[1] + n2*q[2]
+    r[2] = -n2*q[1] + n1*q[2]
+
+    return r
+
+
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+	
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/(h_left +  h0/h_left)
+
+
+    uh_left = u_left*h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/(h_right + h0/h_right)
+
+    uh_right = u_right*h_right
+	
+
+    #We have got   h   and    u   at vertex, then  the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+# 
+def compute_fluxes(domain):
+    """
+	Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+
+    domain.distribute_to_vertices_and_edges()
+    domain.update_boundary()
+    
+    N = domain.number_of_elements
+    Stage = domain.quantities['stage']
+    Xmom  = domain.quantities['xmomentum']
+    Bed   = domain.quantities['elevation']
+
+    stage = Stage.vertex_values
+    xmom  = Xmom.vertex_values
+    bed   = Bed.vertex_values
+    
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+
+
+	
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+     
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+
+        Stage.explicit_update[k] = flux[0]
+        Xmom.explicit_update[k] =  flux[1]
+        #Ymom.explicit_update[k] = flux[2]
+        #print "flux cell",k,flux[0]
+
+    domain.flux_timestep = timestep
+    #print domain.quantities['stage'].centroid_values
+#
+def compute_timestep(domain):
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                qr[0] = stage[n, m]
+                qr[1] =  xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #Update optimal_timestep
+            try:
+                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+    domain.timestep = timestep
+
+# Compute flux definition
+def compute_fluxes_C_long(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+    #print 'timestep=',timestep
+    #print 'The type of timestep is',type(timestep)
+	
+    epsilon = domain.epsilon
+    #print 'epsilon=',epsilon
+    #print 'The type of epsilon is',type(epsilon)
+	
+    g = domain.g
+    #print 'g=',g
+    #print 'The type of g is',type(g)
+	
+    neighbours = domain.neighbours
+    #print 'neighbours=',neighbours
+    #print 'The type of neighbours is',type(neighbours)
+	
+    neighbour_vertices = domain.neighbour_vertices
+    #print 'neighbour_vertices=',neighbour_vertices
+    #print 'The type of neighbour_vertices is',type(neighbour_vertices)
+	
+    normals = domain.normals
+    #print 'normals=',normals
+    #print 'The type of normals is',type(normals)
+	
+    areas = domain.areas
+    #print 'areas=',areas
+    #print 'The type of areas is',type(areas)
+	
+    stage_edge_values = domain.quantities['stage'].vertex_values
+    #print 'stage_edge_values=',stage_edge_values
+    #print 'The type of stage_edge_values is',type(stage_edge_values)
+	
+    xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+    #print 'xmom_edge_values=',xmom_edge_values
+    #print 'The type of xmom_edge_values is',type(xmom_edge_values)
+	
+    bed_edge_values = domain.quantities['elevation'].vertex_values
+    #print 'bed_edge_values=',bed_edge_values
+    #print 'The type of bed_edge_values is',type(bed_edge_values)
+	
+    stage_boundary_values = domain.quantities['stage'].boundary_values
+    #print 'stage_boundary_values=',stage_boundary_values
+    #print 'The type of stage_boundary_values is',type(stage_boundary_values)
+	
+    xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+    #print 'xmom_boundary_values=',xmom_boundary_values
+    #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
+	
+    stage_explicit_update = domain.quantities['stage'].explicit_update
+    #print 'stage_explicit_update=',stage_explicit_update
+    #print 'The type of stage_explicit_update is',type(stage_explicit_update)
+	
+    xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+    #print 'xmom_explicit_update=',xmom_explicit_update
+    #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
+	
+    number_of_elements = len(stage_edge_values)
+    #print 'number_of_elements=',number_of_elements
+    #print 'The type of number_of_elements is',type(number_of_elements)
+	
+    max_speed_array = domain.max_speed_array
+    #print 'max_speed_array=',max_speed_array
+    #print 'The type of max_speed_array is',type(max_speed_array)
+        
+        
+    from comp_flux_ext import compute_fluxes_ext
+	
+    domain.flux_timestep = compute_fluxes_ext(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+
+# Compute flux definition
+def compute_fluxes_C_short(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    stage = domain.quantities['stage']
+    xmom = domain.quantities['xmomentum']
+    bed = domain.quantities['elevation']
+
+
+    from comp_flux_ext import compute_fluxes_ext_short
+
+    domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+
+
+
+# ###################################
+def compute_fluxes_C_wellbalanced(domain):
+    #from Numeric import zeros, Float
+    #import sys
+
+	
+    #timestep = float(sys.maxint)
+    #epsilon = domain.epsilon
+    #g = domain.g
+    #neighbours = domain.neighbours
+    #neighbour_vertices = domain.neighbour_vertices
+    #normals = domain.normals
+    #areas = domain.areas
+    #stage_edge_values = domain.quantities['stage'].vertex_values
+    #xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+    #bed_edge_values = domain.quantities['elevation'].vertex_values
+    #stage_boundary_values = domain.quantities['stage'].boundary_values
+    #xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+    #stage_explicit_update = domain.quantities['stage'].explicit_update
+    #xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+    #number_of_elements = len(stage_edge_values)
+    #max_speed_array = domain.max_speed_array
+
+    import sys
+    from Numeric import zeros, Float
+    
+    N = domain.number_of_elements
+    timestep = float(sys.maxint)
+    epsilon = domain.epsilon
+    g = domain.g
+    neighbours = domain.neighbours
+    neighbour_vertices = domain.neighbour_vertices
+    normals = domain.normals
+    areas = domain.areas
+    
+    Stage = domain.quantities['stage']
+    Xmom  = domain.quantities['xmomentum']
+    Bed   = domain.quantities['elevation']
+
+    stage_boundary_values = Stage.boundary_values
+    xmom_boundary_values =  Xmom.boundary_values
+    stage_explicit_update = Stage.explicit_update
+    xmom_explicit_update = Xmom.explicit_update
+    max_speed_array = domain.max_speed_array
+    
+    domain.distribute_to_vertices_and_edges()
+    domain.update_boundary()
+    stage_V = Stage.vertex_values
+    xmom_V  = Xmom.vertex_values
+    bed_V   = Bed.vertex_values
+    #h_V     = Height.vertex_values
+    #u_V     = Velocity.vertex_values
+
+    number_of_elements = len(stage_V)
+
+    #flux = zeros(2, Float) #Work array for summing up fluxes
+    #ql = zeros(2, Float)
+    #qr = zeros(2, Float)
+	        
+    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
+	
+    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_V,
+                                  xmom_V,
+                                  bed_V,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+# ###################################
+
+
+
+
+
+
+# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+
+def distribute_to_vertices_and_edges(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    #from config import optimised_gradient_limiter
+
+    #Remove very thin layers of water
+    #protect_against_infinitesimal_and_negative_heights(domain)  
+
+    import sys
+    from Numeric import zeros, Float
+    from config import epsilon, h0
+	
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+
+    #Arrays   
+    w_C   = Stage.centroid_values    
+    uh_C  = Xmom.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+	
+    #print id(h_C)
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]						
+        if h_C[i] <= 0.0:
+	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
+	    h_C[i] = 1.0e-15
+	    w_C[i] = z_C[i]
+	    uh_C[i] = 0.0
+		
+    
+##    for i in range(len(h_C)):
+##	if h_C[i] < epsilon:
+##	    u_C[i] = 0.0  #Could have been negative
+##	    h_C[i] = 0.0
+##	else:
+
+    u_C[:]  = uh_C/(h_C + h0/h_C)
+	
+    for name in [ 'velocity', 'stage' ]:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            #print "add extrapolate second order to shallow water"
+            #if name != 'height':
+            Q.extrapolate_second_order()
+            #Q.limit()
+        else:
+            raise 'Unknown order'
+
+    stage_V = domain.quantities['stage'].vertex_values			
+    bed_V   = domain.quantities['elevation'].vertex_values	
+    h_V    = domain.quantities['height'].vertex_values
+    u_V    = domain.quantities['velocity'].vertex_values		
+    xmom_V = domain.quantities['xmomentum'].vertex_values	
+
+    h_V[:]    = stage_V - bed_V
+    for i in range(len(h_C)):
+        for j in range(2):
+            if h_V[i,j] < 0.0 :
+                print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
+                dh = h_V[i,j]
+                h_V[i,j] = 0.0
+                stage_V[i,j] = bed_V[i,j]
+                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
+                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
+                
+    xmom_V[:] = u_V * h_V
+    
+    return
+#
+
+
+
+
+
+
+
+#
+def protect_against_infinitesimal_and_negative_heights(domain):
+    """Protect against infinitesimal heights and associated high velocities
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    xmomc = domain.quantities['xmomentum'].centroid_values
+#    ymomc = domain.quantities['ymomentum'].centroid_values
+    hc = wc - zc  #Water depths at centroids
+
+    zv = domain.quantities['elevation'].vertex_values
+    wv = domain.quantities['stage'].vertex_values
+    hv = wv-zv
+    xmomv = domain.quantities['xmomentum'].vertex_values
+    #remove the above two lines and corresponding code below
+
+    #Update
+    for k in range(domain.number_of_elements):
+
+        if hc[k] < domain.minimum_allowed_height:
+            #Control stage
+            if hc[k] < domain.epsilon:
+                wc[k] = zc[k] # Contain 'lost mass' error
+                wv[k,0] = zv[k,0]
+                wv[k,1] = zv[k,1]
+                
+            xmomc[k] = 0.0
+            
+        #N = domain.number_of_elements
+        #if (k == 0) | (k==N-1):
+        #    wc[k] = zc[k] # Contain 'lost mass' error
+        #    wv[k,0] = zv[k,0]
+        #    wv[k,1] = zv[k,1]
+    
+def h_limiter(domain):
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    hvbar = zeros(hv.shape, Float) #h-limited values
+
+    #Find min and max of this and neighbour's centroid values
+    hmax = zeros(hc.shape, Float)
+    hmin = zeros(hc.shape, Float)
+
+    for k in range(N):
+        hmax[k] = hmin[k] = hc[k]
+        #for i in range(3):
+        for i in range(2):   
+            n = domain.neighbours[k,i]
+            if n >= 0:
+                hn = hc[n] #Neighbour's centroid value
+
+                hmin[k] = min(hmin[k], hn)
+                hmax[k] = max(hmax[k], hn)
+
+
+    #Diffences between centroids and maxima/minima
+    dhmax = hmax - hc
+    dhmin = hmin - hc
+
+    #Deltas between vertex and centroid values
+    dh = zeros(hv.shape, Float)
+    #for i in range(3):
+    for i in range(2):
+        dh[:,i] = hv[:,i] - hc
+
+    #Phi limiter
+    for k in range(N):
+
+        #Find the gradient limiter (phi) across vertices
+        phi = 1.0
+        #for i in range(3):
+        for i in range(2):
+            r = 1.0
+            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
+            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
+
+            phi = min( min(r*beta_h, 1), phi )
+
+        #Then update using phi limiter
+        #for i in range(3):
+        for i in range(2):
+            hvbar[k,i] = hc[k] + phi*dh[k,i]
+
+    return hvbar
+
+def balance_deep_and_shallow(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters limiting using w, and stage computed by
+    gradient-limiters limiting using h (h-limiter).
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+
+    The h-limiter is always applied irrespective of the order.
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Limit h
+    hvbar = h_limiter(domain)
+
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]))#,
+#                 abs(zv[k,2]-zc[k]))
+
+
+        hmin = min( hv[k,:] )
+
+        #Create alpha in [0,1], where alpha==0 means using the h-limited
+        #stage and alpha==1 means using the w-limited stage as
+        #computed by the gradient limiter (both 1st or 2nd order)
+
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed
+            alpha = 1.0
+
+        alpha  = 0.0
+	#Let
+	#
+	#  wvi be the w-limited stage (wvi = zvi + hvi)
+	#  wvi- be the h-limited state (wvi- = zvi + hvi-)
+	#
+	#
+	#where i=0,1,2 denotes the vertex ids
+	#
+        #Weighted balance between w-limited and h-limited stage is
+        #
+	#  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+	#
+        #It follows that the updated wvi is
+        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
+        #
+	# Momentum is balanced between constant and limited
+
+
+        #for i in range(3):
+        #    wv[k,i] = zv[k,i] + hvbar[k,i]
+
+	#return
+
+	if alpha < 1:
+
+            #for i in range(3):
+            for i in range(2):
+                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+#            ymomc = domain.quantities['ymomentum'].centroid_values
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+#           ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
+#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+
+
+###############################################
+#Boundaries - specific to the shallow water wave equation
+class Reflective_boundary(Boundary):
+    """Reflective boundary returns same conserved quantities as
+    those present in its neighbour volume but reflected.
+
+    This class is specific to the shallow water equation as it
+    works with the momentum quantities assumed to be the second
+    and third conserved quantities.
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for reflective boundary'
+            raise msg
+
+        #Handy shorthands
+        #self.stage   = domain.quantities['stage'].edge_values
+        #self.xmom    = domain.quantities['xmomentum'].edge_values
+        #self.ymom    = domain.quantities['ymomentum'].edge_values
+        self.normals = domain.normals
+        self.stage   = domain.quantities['stage'].vertex_values
+        self.xmom    = domain.quantities['xmomentum'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.conserved_quantities = zeros(2, Float)
+
+    def __repr__(self):
+        return 'Reflective_boundary'
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Reflective boundaries reverses the outward momentum
+        of the volume they serve.
+        """
+
+        q = self.conserved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = self.xmom[vol_id, edge_id]
+        #q[2] = self.ymom[vol_id, edge_id]
+        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
+        normal = self.normals[vol_id,edge_id]
+
+        #r = rotate(q, normal, direction = 1)
+        #r[1] = -r[1]
+        #q = rotate(r, normal, direction = -1)
+        r = q
+        r[1] = normal*r[1]
+        r[1] = -r[1]
+        r[1] = normal*r[1]
+        q = r
+        #For start interval there is no outward momentum so do not need to
+        #reverse direction in this case
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, conserved_quantities=None):
+        Boundary.__init__(self)
+
+        if conserved_quantities is None:
+            msg = 'Must specify one value for each conserved quantity'
+            raise msg
+
+        from Numeric import array, Float
+        self.conserved_quantities=array(conserved_quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.conserved_quantities
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.conserved_quantities
+
+
+#########################
+#Standard forcing terms:
+#
+def gravity(domain):
+    """Apply gravitational pull in the presence of bed slope
+    """
+
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+
+    xmom = domain.quantities['xmomentum'].explicit_update
+    stage = domain.quantities['stage'].explicit_update
+#    ymom = domain.quantities['ymomentum'].explicit_update
+
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
+    #h = Stage.edge_values - Elevation.edge_values
+    h = Stage.vertex_values - Elevation.vertex_values
+    b = Elevation.vertex_values
+    w = Stage.vertex_values
+
+    x = domain.get_vertex_coordinates()
+    g = domain.g
+
+    for k in range(domain.number_of_elements):
+#        avg_h = sum( h[k,:] )/3
+        avg_h = sum( h[k,:] )/2
+
+        #Compute bed slope
+        #x0, y0, x1, y1, x2, y2 = x[k,:]
+        x0, x1 = x[k,:]
+        #z0, z1, z2 = v[k,:]
+        b0, b1 = b[k,:]
+
+        w0, w1 = w[k,:]
+        wx = gradient(x0, x1, w0, w1)
+
+        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+        bx = gradient(x0, x1, b0, b1)
+        
+        #Update momentum (explicit update is reset to source values)
+        xmom[k] += -g*bx*avg_h
+        #xmom[k] = -g*bx*avg_h
+        #stage[k] = 0.0
+ 
+ 
+def manning_friction(domain):
+    """Apply (Manning) friction to water momentum
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+    #vh = domain.quantities['ymomentum'].centroid_values
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g
+
+    for k in range(N):
+        if eta[k] >= eps:
+            if h[k] >= eps:
+            	#S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+                S = -g * eta[k]**2 * uh[k]
+            	S /= h[k]**(7.0/3)
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+            	#ymom_update[k] += S*vh[k]
+
+def linear_friction(domain):
+    """Apply linear friction to water momentum
+
+    Assumes quantity: 'linear_friction' to be present
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+#    vh = domain.quantities['ymomentum'].centroid_values
+    tau = domain.quantities['linear_friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g #Not necessary? Why was this added?
+
+    for k in range(N):
+        if tau[k] >= eps:
+            if h[k] >= eps:
+            	S = -tau[k]/h[k]
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+ #           	ymom_update[k] += S*vh[k]
+
+
+
+def check_forcefield(f):
+    """Check that f is either
+    1: a callable object f(t,x,y), where x and y are vectors
+       and that it returns an array or a list of same length
+       as x and y
+    2: a scalar
+    """
+
+    from Numeric import ones, Float, array
+
+
+    if callable(f):
+        #N = 3
+        N = 2
+        #x = ones(3, Float)
+        #y = ones(3, Float)
+        x = ones(2, Float)
+        #y = ones(2, Float)
+        
+        try:
+            #q = f(1.0, x=x, y=y)
+            q = f(1.0, x=x)
+        except Exception, e:
+            msg = 'Function %s could not be executed:\n%s' %(f, e)
+	    #FIXME: Reconsider this semantics
+            raise msg
+
+        try:
+            q = array(q).astype(Float)
+        except:
+            msg = 'Return value from vector function %s could ' %f
+            msg += 'not be converted into a Numeric array of floats.\n'
+            msg += 'Specified function should return either list or array.'
+            raise msg
+
+        #Is this really what we want?
+        msg = 'Return vector from function %s ' %f
+        msg += 'must have same lenght as input vectors'
+        assert len(q) == N, msg
+
+    else:
+        try:
+            f = float(f)
+        except:
+            msg = 'Force field %s must be either a scalar' %f
+            msg += ' or a vector function'
+            raise msg
+    return f
+
+class Wind_stress:
+    """Apply wind stress to water momentum in terms of
+    wind speed [m/s] and wind direction [degrees]
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Initialise windfield from wind speed s [m/s]
+        and wind direction phi [degrees]
+
+        Inputs v and phi can be either scalars or Python functions, e.g.
+
+        W = Wind_stress(10, 178)
+
+        #FIXME - 'normal' degrees are assumed for now, i.e. the
+        vector (1,0) has zero degrees.
+        We may need to convert from 'compass' degrees later on and also
+        map from True north to grid north.
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def speed(t,x,y):
+            ...
+            return s
+
+        def angle(t,x,y):
+            ...
+            return phi
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        W = Wind_stress(speed, angle)
+
+        The instantiated object W can be appended to the list of
+        forcing_terms as in
+
+        Alternatively, one vector valued function for (speed, angle)
+        can be applied, providing both quantities simultaneously.
+        As in
+        W = Wind_stress(F), where returns (speed, angle) for each t.
+
+        domain.forcing_terms.append(W)
+        """
+
+        from config import rho_a, rho_w, eta_w
+        from Numeric import array, Float
+
+        if len(args) == 2:
+            s = args[0]
+            phi = args[1]
+        elif len(args) == 1:
+            #Assume vector function returning (s, phi)(t,x,y)
+            vector_function = args[0]
+            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+            s = lambda t,x: vector_function(t,x=x)[0]
+            phi = lambda t,x: vector_function(t,x=x)[1]
+        else:
+           #Assume info is in 2 keyword arguments
+
+           if len(kwargs) == 2:
+               s = kwargs['s']
+               phi = kwargs['phi']
+           else:
+               raise 'Assumes two keyword arguments: s=..., phi=....'
+
+        print 'phi', phi
+        self.speed = check_forcefield(s)
+        self.phi = check_forcefield(phi)
+
+        self.const = eta_w*rho_a/rho_w
+
+
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import Float, ones, ArrayType
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        #ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        N = domain.number_of_elements
+        t = domain.time
+
+        if callable(self.speed):
+            xc = domain.get_centroid_coordinates()
+            #s_vec = self.speed(t, xc[:,0], xc[:,1])
+            s_vec = self.speed(t, xc)
+        else:
+            #Assume s is a scalar
+
+            try:
+                s_vec = self.speed * ones(N, Float)
+            except:
+                msg = 'Speed must be either callable or a scalar: %s' %self.s
+                raise msg
+
+
+        if callable(self.phi):
+            xc = domain.get_centroid_coordinates()
+            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
+            phi_vec = self.phi(t, xc)
+        else:
+            #Assume phi is a scalar
+
+            try:
+                phi_vec = self.phi * ones(N, Float)
+            except:
+                msg = 'Angle must be either callable or a scalar: %s' %self.phi
+                raise msg
+
+        #assign_windfield_values(xmom_update, ymom_update,
+        #                        s_vec, phi_vec, self.const)
+        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
+
+
+#def assign_windfield_values(xmom_update, ymom_update,
+#                            s_vec, phi_vec, const):
+def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
+    """Python version of assigning wind field to update vectors.
+    A c version also exists (for speed)
+    """
+    from math import pi, cos, sin, sqrt
+
+    N = len(s_vec)
+    for k in range(N):
+        s = s_vec[k]
+        phi = phi_vec[k]
+
+        #Convert to radians
+        phi = phi*pi/180
+
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress
+        #S = const * sqrt(u**2 + v**2)
+        S = const * u
+        xmom_update[k] += S*u
+        #ymom_update[k] += S*v
Index: /anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py	(revision 5827)
+++ /anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py	(revision 5827)
@@ -0,0 +1,1431 @@
+"""Class Domain -
+1D interval domains for finite-volume computations of
+the shallow water wave equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the Shallow Water Wave Equation
+
+
+U_t + E_x = S
+
+where
+
+U = [w, uh]
+E = [uh, u^2h + gh^2/2]
+S represents source terms forcing the system
+(e.g. gravity, friction, wind stress, ...)
+
+and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
+
+The quantities are
+
+symbol    variable name    explanation
+x         x                horizontal distance from origin [m]
+z         elevation        elevation of bed on which flow is modelled [m]
+h         height           water height above z [m]
+w         stage            absolute water level, w = z+h [m]
+u                          speed in the x direction [m/s]
+uh        xmomentum        momentum in the x direction [m^2/s]
+
+eta                        mannings friction coefficient [to appear]
+nu                         wind stress coefficient [to appear]
+
+The conserved quantities are w, uh
+
+For details see e.g.
+Christopher Zoppou and Stephen Roberts,
+Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
+Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
+
+
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2006
+"""
+
+
+from domain import *
+Generic_Domain = Domain #Rename
+
+#Shallow water domain
+class Domain(Generic_Domain):
+
+    def __init__(self, coordinates, boundary = None, tagged_elements = None):
+
+        conserved_quantities = ['stage', 'xmomentum']
+        other_quantities = ['elevation', 'friction', 'height', 'velocity']
+        Generic_Domain.__init__(self, coordinates, boundary,
+                                conserved_quantities, other_quantities,
+                                tagged_elements)
+        
+        from config import minimum_allowed_height, g, h0
+        self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
+        self.h0 = h0
+
+        #forcing terms not included in 1d domain ?WHy?
+        self.forcing_terms.append(gravity)
+        #self.forcing_terms.append(manning_friction)
+        #print "\nI have Removed forcing terms line 64 1dsw"
+
+        #Realtime visualisation
+        self.visualiser = None
+        self.visualise  = False
+        self.visualise_color_stage = False
+        self.visualise_stage_range = 1.0
+        self.visualise_timer = True
+        self.visualise_range_z = None
+        
+        #Stored output
+        self.store = True
+        self.format = 'sww'
+        self.smooth = True
+
+        #Evolve parametrs
+        self.cfl = 1.0
+        
+        #Reduction operation for get_vertex_values
+        from util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.quantities_to_be_stored = ['stage','xmomentum']
+
+	self.__doc__ = 'shallow_water_domain'
+
+
+    def set_quantities_to_be_stored(self, q):
+        """Specify which quantities will be stored in the sww file.
+
+        q must be either:
+          - the name of a quantity
+          - a list of quantity names
+          - None
+
+        In the two first cases, the named quantities will be stored at each
+        yieldstep
+        (This is in addition to the quantities elevation and friction)  
+        If q is None, storage will be switched off altogether.
+        """
+
+
+        if q is None:
+            self.quantities_to_be_stored = []            
+            self.store = False
+            return
+
+        if isinstance(q, basestring):
+            q = [q] # Turn argument into a list
+
+        #Check correcness    
+        for quantity_name in q:
+            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
+            assert quantity_name in self.conserved_quantities, msg 
+        
+        self.quantities_to_be_stored = q
+        
+
+    def initialise_visualiser(self,scale_z=1.0,rect=None):
+        #Realtime visualisation
+        if self.visualiser is None:
+            from realtime_visualisation_new import Visualiser
+            self.visualiser = Visualiser(self,scale_z,rect)
+            self.visualiser.setup['elevation']=True
+            self.visualiser.updating['stage']=True
+        self.visualise = True
+        if self.visualise_color_stage == True:
+            self.visualiser.coloring['stage'] = True
+            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
+        print 'initialise visualiser'
+        print self.visualiser.setup
+        print self.visualiser.updating
+
+    def check_integrity(self):
+        Generic_Domain.check_integrity(self)
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "stage"'
+        assert self.conserved_quantities[0] == 'stage', msg
+        msg = 'Second conserved quantity must be "xmomentum"'
+        assert self.conserved_quantities[1] == 'xmomentum', msg
+
+    def extrapolate_second_order_sw(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_sw(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_C_wellbalanced(self)
+
+    def compute_timestep(self):
+        #Call correct module function
+        compute_timestep(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges(self)
+        
+    def evolve(self, yieldstep = None, finaltime = None, duration = None,
+               skip_initial_step = False):
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call check integrity here rather than from user scripts
+        #self.check_integrity()
+
+        #msg = 'Parameter beta_h must be in the interval [0, 1)'
+        #assert 0 <= self.beta_h < 1.0, msg
+        #msg = 'Parameter beta_w must be in the interval [0, 1)'
+        #assert 0 <= self.beta_w < 1.0, msg
+
+
+        #Initial update of vertex and edge values before any storage
+        #and or visualisation
+        
+        #self.distribute_to_vertices_and_edges() ?????????????????????????????????
+        
+        #Initialise real time viz if requested
+        #if self.visualise is True and self.time == 0.0:
+        #    if self.visualiser is None:
+        #        self.initialise_visualiser()
+        #
+        #    self.visualiser.update_timer()
+        #    self.visualiser.setup_all()
+
+        #Store model data, e.g. for visualisation
+        #if self.store is True and self.time == 0.0:
+        #    self.initialise_storage()
+        #    #print 'Storing results in ' + self.writer.filename
+        #else:
+        #    pass
+        #    #print 'Results will not be stored.'
+        #    #print 'To store results set domain.store = True'
+        #    #FIXME: Diagnostic output should be controlled by
+        #    # a 'verbose' flag living in domain (or in a parent class)
+
+        #Call basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
+                                       skip_initial_step):
+            #Real time viz
+        #    if self.visualise is True:
+        #        self.visualiser.update_all()
+        #        self.visualiser.update_timer()
+
+
+            #Store model data, e.g. for subsequent visualisation
+        #    if self.store is True:
+        #        self.store_timestep(self.quantities_to_be_stored)
+
+            #FIXME: Could maybe be taken from specified list
+            #of 'store every step' quantities
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+    def initialise_storage(self):
+        """Create and initialise self.writer object for storing data.
+        Also, save x and bed elevation
+        """
+
+        import data_manager
+
+        #Initialise writer
+        self.writer = data_manager.get_dataobject(self, mode = 'w')
+
+        #Store vertices and connectivity
+        self.writer.store_connectivity()
+
+
+    def store_timestep(self, name):
+        """Store named quantity and time.
+
+        Precondition:
+           self.write has been initialised
+        """
+        self.writer.store_timestep(name)
+
+
+#=============== End of Shallow Water Domain ===============================
+
+#Rotation of momentum vector
+def rotate(q, normal, direction = 1):
+    """Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+
+    If direction is negative the rotation is inverted.
+
+    Input vector is preserved
+
+    This function is specific to the shallow water wave equation
+    """
+
+    from Numeric import zeros, Float
+
+    assert len(q) == 3,\
+           'Vector of conserved quantities must have length 3'\
+           'for 2D shallow water equation'
+
+    try:
+        l = len(normal)
+    except:
+        raise 'Normal vector must be an Numeric array'
+
+    assert l == 2, 'Normal vector must have 2 components'
+
+
+    n1 = normal[0]
+    n2 = normal[1]
+
+    r = zeros(len(q), Float) #Rotated quantities
+    r[0] = q[0]              #First quantity, height, is not rotated
+
+    if direction == -1:
+        n2 = -n2
+
+
+    r[1] =  n1*q[1] + n2*q[2]
+    r[2] = -n2*q[1] + n1*q[2]
+
+    return r
+
+
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+	
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/(h_left +  h0/h_left)
+
+
+    uh_left = u_left*h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/(h_right + h0/h_right)
+
+    uh_right = u_right*h_right
+	
+
+    #We have got   h   and    u   at vertex, then  the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+# 
+def compute_fluxes(domain):
+    """
+	Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    
+    """
+
+    import sys
+    from Numeric import zeros, Float
+    N = domain.number_of_elements
+    Stage = domain.quantities['stage']
+    Xmom  = domain.quantities['xmomentum']
+    Bed   = domain.quantities['elevation']
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    stage, xmom, bed, h_V, u_V = distribute_to_vertices_and_edges(domain)
+    domain.update_boundary()
+	
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+     
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+
+        Stage.explicit_update[k] = flux[0]
+        Xmom.explicit_update[k] =  flux[1]
+        #Ymom.explicit_update[k] = flux[2]
+        #print "flux cell",k,flux[0]
+
+    domain.flux_timestep = timestep
+    #print domain.quantities['stage'].centroid_values
+#
+def compute_timestep(domain):
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                qr[0] = stage[n, m]
+                qr[1] =  xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #Update optimal_timestep
+            try:
+                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+    domain.timestep = timestep
+
+# Compute flux definition
+def compute_fluxes_C_long(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+    #print 'timestep=',timestep
+    #print 'The type of timestep is',type(timestep)
+	
+    epsilon = domain.epsilon
+    #print 'epsilon=',epsilon
+    #print 'The type of epsilon is',type(epsilon)
+	
+    g = domain.g
+    #print 'g=',g
+    #print 'The type of g is',type(g)
+	
+    neighbours = domain.neighbours
+    #print 'neighbours=',neighbours
+    #print 'The type of neighbours is',type(neighbours)
+	
+    neighbour_vertices = domain.neighbour_vertices
+    #print 'neighbour_vertices=',neighbour_vertices
+    #print 'The type of neighbour_vertices is',type(neighbour_vertices)
+	
+    normals = domain.normals
+    #print 'normals=',normals
+    #print 'The type of normals is',type(normals)
+	
+    areas = domain.areas
+    #print 'areas=',areas
+    #print 'The type of areas is',type(areas)
+	
+    stage_edge_values = domain.quantities['stage'].vertex_values
+    #print 'stage_edge_values=',stage_edge_values
+    #print 'The type of stage_edge_values is',type(stage_edge_values)
+	
+    xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+    #print 'xmom_edge_values=',xmom_edge_values
+    #print 'The type of xmom_edge_values is',type(xmom_edge_values)
+	
+    bed_edge_values = domain.quantities['elevation'].vertex_values
+    #print 'bed_edge_values=',bed_edge_values
+    #print 'The type of bed_edge_values is',type(bed_edge_values)
+	
+    stage_boundary_values = domain.quantities['stage'].boundary_values
+    #print 'stage_boundary_values=',stage_boundary_values
+    #print 'The type of stage_boundary_values is',type(stage_boundary_values)
+	
+    xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+    #print 'xmom_boundary_values=',xmom_boundary_values
+    #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
+	
+    stage_explicit_update = domain.quantities['stage'].explicit_update
+    #print 'stage_explicit_update=',stage_explicit_update
+    #print 'The type of stage_explicit_update is',type(stage_explicit_update)
+	
+    xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+    #print 'xmom_explicit_update=',xmom_explicit_update
+    #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
+	
+    number_of_elements = len(stage_edge_values)
+    #print 'number_of_elements=',number_of_elements
+    #print 'The type of number_of_elements is',type(number_of_elements)
+	
+    max_speed_array = domain.max_speed_array
+    #print 'max_speed_array=',max_speed_array
+    #print 'The type of max_speed_array is',type(max_speed_array)
+        
+        
+    from comp_flux_ext import compute_fluxes_ext
+	
+    domain.flux_timestep = compute_fluxes_ext(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+
+# Compute flux definition
+def compute_fluxes_C_short(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    stage = domain.quantities['stage']
+    xmom = domain.quantities['xmomentum']
+    bed = domain.quantities['elevation']
+
+
+    from comp_flux_ext import compute_fluxes_ext_short
+
+    domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+
+
+
+# ###################################
+def compute_fluxes_C_wellbalanced(domain):
+    import sys
+    from Numeric import zeros, Float
+    from config import epsilon, h0
+    
+    N = domain.number_of_elements
+    timestep = float(sys.maxint)
+    #epsilon = domain.epsilon
+    g = domain.g
+    neighbours = domain.neighbours
+    neighbour_vertices = domain.neighbour_vertices
+    normals = domain.normals
+    areas = domain.areas
+    
+    Stage = domain.quantities['stage']
+    Xmom  = domain.quantities['xmomentum']
+    Bed   = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+
+    w_C   = Stage.centroid_values    
+    uh_C  = Xmom.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]						
+        if h_C[i] < epsilon :
+	    h_C[i] = 0.0
+	    w_C[i] = z_C[i]
+	    uh_C[i] = 0.0
+
+    for i in range(len(h_C)):
+	if h_C[i] < epsilon:
+	    u_C[i] = 0.0  #Could have been negative
+	    h_C[i] = 0.0
+	else:
+	    u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
+
+    for name in ['height', 'velocity']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            #print "add extrapolate second order to shallow water"
+            #if name != 'height':
+            Q.extrapolate_second_order()
+            #Q.limit()
+        else:
+            raise 'Unknown order'
+
+    w_V = domain.quantities['stage'].vertex_values			
+    h_V = domain.quantities['height'].vertex_values
+    u_V = domain.quantities['velocity'].vertex_values		
+    uh_V = domain.quantities['xmomentum'].vertex_values
+    z_V   = w_V - h_V
+
+    
+    stage_boundary_values = Stage.boundary_values
+    xmom_boundary_values =  Xmom.boundary_values
+    stage_explicit_update = Stage.explicit_update
+    xmom_explicit_update = Xmom.explicit_update
+    max_speed_array = domain.max_speed_array
+    #number_of_elements = len(w_V)
+    
+    #domain.distribute_to_vertices_and_edges()
+    #domain.update_boundary()
+    #stage_V = Stage.vertex_values
+    #xmom_V  = Xmom.vertex_values
+    #bed_V   = Bed.vertex_values
+    
+        
+    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
+	
+    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  w_V,
+                                  uh_V,
+                                  z_V,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  N,
+                                  max_speed_array)
+
+# ###################################
+
+
+
+
+
+
+# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+
+def distribute_to_vertices_and_edges(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    #from config import optimised_gradient_limiter
+
+    #Remove very thin layers of water
+    protect_against_infinitesimal_and_negative_heights(domain)  
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+
+    #Arrays   
+    w_C   = Stage.centroid_values    
+    uh_C  = Xmom.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+	
+    #print id(h_C)
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]						#This is h at centroids:  OK
+    
+    from config import epsilon, h0
+	
+    for i in range(len(h_C)):
+	if h_C[i] < epsilon:
+	    u_C[i] = 0.0  #Could have been negative
+	    h_C[i] = 0.0
+	else:
+	    u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
+	
+    for name in ['stage', 'velocity', 'height' ]:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            #print "add extrapolate second order to shallow water"
+            #if name != 'height':
+            Q.extrapolate_second_order()
+            #Q.limit()
+        else:
+            raise 'Unknown order'
+
+    stage  = domain.quantities['stage'].vertex_values			#This is w at vertices:   OK
+    h_V    = domain.quantities['height'].vertex_values
+    bed    = stage - h_V
+    u_V    = domain.quantities['velocity'].vertex_values		
+    xmom_V = domain.quantities['xmomentum'].vertex_values	
+	
+    #h_V[:]    = stage - bed
+    xmom_V[:] = u_V * h_V
+    return stage, xmom_V, bed, h_V, u_V
+#
+def protect_against_infinitesimal_and_negative_heights(domain):
+    """Protect against infinitesimal heights and associated high velocities
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    xmomc = domain.quantities['xmomentum'].centroid_values
+#    ymomc = domain.quantities['ymomentum'].centroid_values
+    hc = wc - zc  #Water depths at centroids
+
+    zv = domain.quantities['elevation'].vertex_values
+    wv = domain.quantities['stage'].vertex_values
+    hv = wv-zv
+    xmomv = domain.quantities['xmomentum'].vertex_values
+    #remove the above two lines and corresponding code below
+
+    #Update
+    for k in range(domain.number_of_elements):
+
+        if hc[k] < domain.minimum_allowed_height:
+            #Control stage
+            if hc[k] < domain.epsilon:
+                wc[k] = zc[k] # Contain 'lost mass' error
+                wv[k,0] = zv[k,0]
+                wv[k,1] = zv[k,1]
+                
+            xmomc[k] = 0.0
+            
+        #N = domain.number_of_elements
+        #if (k == 0) | (k==N-1):
+        #    wc[k] = zc[k] # Contain 'lost mass' error
+        #    wv[k,0] = zv[k,0]
+        #    wv[k,1] = zv[k,1]
+    
+def h_limiter(domain):
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    hvbar = zeros(hv.shape, Float) #h-limited values
+
+    #Find min and max of this and neighbour's centroid values
+    hmax = zeros(hc.shape, Float)
+    hmin = zeros(hc.shape, Float)
+
+    for k in range(N):
+        hmax[k] = hmin[k] = hc[k]
+        #for i in range(3):
+        for i in range(2):   
+            n = domain.neighbours[k,i]
+            if n >= 0:
+                hn = hc[n] #Neighbour's centroid value
+
+                hmin[k] = min(hmin[k], hn)
+                hmax[k] = max(hmax[k], hn)
+
+
+    #Diffences between centroids and maxima/minima
+    dhmax = hmax - hc
+    dhmin = hmin - hc
+
+    #Deltas between vertex and centroid values
+    dh = zeros(hv.shape, Float)
+    #for i in range(3):
+    for i in range(2):
+        dh[:,i] = hv[:,i] - hc
+
+    #Phi limiter
+    for k in range(N):
+
+        #Find the gradient limiter (phi) across vertices
+        phi = 1.0
+        #for i in range(3):
+        for i in range(2):
+            r = 1.0
+            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
+            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
+
+            phi = min( min(r*beta_h, 1), phi )
+
+        #Then update using phi limiter
+        #for i in range(3):
+        for i in range(2):
+            hvbar[k,i] = hc[k] + phi*dh[k,i]
+
+    return hvbar
+
+def balance_deep_and_shallow(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters limiting using w, and stage computed by
+    gradient-limiters limiting using h (h-limiter).
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+
+    The h-limiter is always applied irrespective of the order.
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Limit h
+    hvbar = h_limiter(domain)
+
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]))#,
+#                 abs(zv[k,2]-zc[k]))
+
+
+        hmin = min( hv[k,:] )
+
+        #Create alpha in [0,1], where alpha==0 means using the h-limited
+        #stage and alpha==1 means using the w-limited stage as
+        #computed by the gradient limiter (both 1st or 2nd order)
+
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed
+            alpha = 1.0
+
+        alpha  = 0.0
+	#Let
+	#
+	#  wvi be the w-limited stage (wvi = zvi + hvi)
+	#  wvi- be the h-limited state (wvi- = zvi + hvi-)
+	#
+	#
+	#where i=0,1,2 denotes the vertex ids
+	#
+        #Weighted balance between w-limited and h-limited stage is
+        #
+	#  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+	#
+        #It follows that the updated wvi is
+        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
+        #
+	# Momentum is balanced between constant and limited
+
+
+        #for i in range(3):
+        #    wv[k,i] = zv[k,i] + hvbar[k,i]
+
+	#return
+
+	if alpha < 1:
+
+            #for i in range(3):
+            for i in range(2):
+                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+#            ymomc = domain.quantities['ymomentum'].centroid_values
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+#           ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
+#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+
+
+###############################################
+#Boundaries - specific to the shallow water wave equation
+class Reflective_boundary(Boundary):
+    """Reflective boundary returns same conserved quantities as
+    those present in its neighbour volume but reflected.
+
+    This class is specific to the shallow water equation as it
+    works with the momentum quantities assumed to be the second
+    and third conserved quantities.
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for reflective boundary'
+            raise msg
+
+        #Handy shorthands
+        #self.stage   = domain.quantities['stage'].edge_values
+        #self.xmom    = domain.quantities['xmomentum'].edge_values
+        #self.ymom    = domain.quantities['ymomentum'].edge_values
+        self.normals = domain.normals
+        self.stage   = domain.quantities['stage'].vertex_values
+        self.xmom    = domain.quantities['xmomentum'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.conserved_quantities = zeros(2, Float)
+
+    def __repr__(self):
+        return 'Reflective_boundary'
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Reflective boundaries reverses the outward momentum
+        of the volume they serve.
+        """
+
+        q = self.conserved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = self.xmom[vol_id, edge_id]
+        #q[2] = self.ymom[vol_id, edge_id]
+        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
+        normal = self.normals[vol_id,edge_id]
+
+        #r = rotate(q, normal, direction = 1)
+        #r[1] = -r[1]
+        #q = rotate(r, normal, direction = -1)
+        r = q
+        r[1] = normal*r[1]
+        r[1] = -r[1]
+        r[1] = normal*r[1]
+        q = r
+        #For start interval there is no outward momentum so do not need to
+        #reverse direction in this case
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, conserved_quantities=None):
+        Boundary.__init__(self)
+
+        if conserved_quantities is None:
+            msg = 'Must specify one value for each conserved quantity'
+            raise msg
+
+        from Numeric import array, Float
+        self.conserved_quantities=array(conserved_quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.conserved_quantities
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.conserved_quantities
+
+
+#########################
+#Standard forcing terms:
+#
+def gravity(domain):
+    """Apply gravitational pull in the presence of bed slope
+    """
+
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+
+    xmom = domain.quantities['xmomentum'].explicit_update
+    stage = domain.quantities['stage'].explicit_update
+#    ymom = domain.quantities['ymomentum'].explicit_update
+
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
+    #h = Stage.edge_values - Elevation.edge_values
+    h = Stage.vertex_values - Elevation.vertex_values
+    b = Elevation.vertex_values
+    w = Stage.vertex_values
+
+    x = domain.get_vertex_coordinates()
+    g = domain.g
+
+    for k in range(domain.number_of_elements):
+#        avg_h = sum( h[k,:] )/3
+        avg_h = sum( h[k,:] )/2
+
+        #Compute bed slope
+        #x0, y0, x1, y1, x2, y2 = x[k,:]
+        x0, x1 = x[k,:]
+        #z0, z1, z2 = v[k,:]
+        b0, b1 = b[k,:]
+
+        w0, w1 = w[k,:]
+        wx = gradient(x0, x1, w0, w1)
+
+        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+        bx = gradient(x0, x1, b0, b1)
+        
+        #Update momentum (explicit update is reset to source values)
+        xmom[k] += -g*bx*avg_h
+        #xmom[k] = -g*bx*avg_h
+        #stage[k] = 0.0
+ 
+ 
+def manning_friction(domain):
+    """Apply (Manning) friction to water momentum
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+    #vh = domain.quantities['ymomentum'].centroid_values
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g
+
+    for k in range(N):
+        if eta[k] >= eps:
+            if h[k] >= eps:
+            	#S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+                S = -g * eta[k]**2 * uh[k]
+            	S /= h[k]**(7.0/3)
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+            	#ymom_update[k] += S*vh[k]
+
+def linear_friction(domain):
+    """Apply linear friction to water momentum
+
+    Assumes quantity: 'linear_friction' to be present
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+#    vh = domain.quantities['ymomentum'].centroid_values
+    tau = domain.quantities['linear_friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g #Not necessary? Why was this added?
+
+    for k in range(N):
+        if tau[k] >= eps:
+            if h[k] >= eps:
+            	S = -tau[k]/h[k]
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+ #           	ymom_update[k] += S*vh[k]
+
+
+
+def check_forcefield(f):
+    """Check that f is either
+    1: a callable object f(t,x,y), where x and y are vectors
+       and that it returns an array or a list of same length
+       as x and y
+    2: a scalar
+    """
+
+    from Numeric import ones, Float, array
+
+
+    if callable(f):
+        #N = 3
+        N = 2
+        #x = ones(3, Float)
+        #y = ones(3, Float)
+        x = ones(2, Float)
+        #y = ones(2, Float)
+        
+        try:
+            #q = f(1.0, x=x, y=y)
+            q = f(1.0, x=x)
+        except Exception, e:
+            msg = 'Function %s could not be executed:\n%s' %(f, e)
+	    #FIXME: Reconsider this semantics
+            raise msg
+
+        try:
+            q = array(q).astype(Float)
+        except:
+            msg = 'Return value from vector function %s could ' %f
+            msg += 'not be converted into a Numeric array of floats.\n'
+            msg += 'Specified function should return either list or array.'
+            raise msg
+
+        #Is this really what we want?
+        msg = 'Return vector from function %s ' %f
+        msg += 'must have same lenght as input vectors'
+        assert len(q) == N, msg
+
+    else:
+        try:
+            f = float(f)
+        except:
+            msg = 'Force field %s must be either a scalar' %f
+            msg += ' or a vector function'
+            raise msg
+    return f
+
+class Wind_stress:
+    """Apply wind stress to water momentum in terms of
+    wind speed [m/s] and wind direction [degrees]
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Initialise windfield from wind speed s [m/s]
+        and wind direction phi [degrees]
+
+        Inputs v and phi can be either scalars or Python functions, e.g.
+
+        W = Wind_stress(10, 178)
+
+        #FIXME - 'normal' degrees are assumed for now, i.e. the
+        vector (1,0) has zero degrees.
+        We may need to convert from 'compass' degrees later on and also
+        map from True north to grid north.
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def speed(t,x,y):
+            ...
+            return s
+
+        def angle(t,x,y):
+            ...
+            return phi
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        W = Wind_stress(speed, angle)
+
+        The instantiated object W can be appended to the list of
+        forcing_terms as in
+
+        Alternatively, one vector valued function for (speed, angle)
+        can be applied, providing both quantities simultaneously.
+        As in
+        W = Wind_stress(F), where returns (speed, angle) for each t.
+
+        domain.forcing_terms.append(W)
+        """
+
+        from config import rho_a, rho_w, eta_w
+        from Numeric import array, Float
+
+        if len(args) == 2:
+            s = args[0]
+            phi = args[1]
+        elif len(args) == 1:
+            #Assume vector function returning (s, phi)(t,x,y)
+            vector_function = args[0]
+            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+            s = lambda t,x: vector_function(t,x=x)[0]
+            phi = lambda t,x: vector_function(t,x=x)[1]
+        else:
+           #Assume info is in 2 keyword arguments
+
+           if len(kwargs) == 2:
+               s = kwargs['s']
+               phi = kwargs['phi']
+           else:
+               raise 'Assumes two keyword arguments: s=..., phi=....'
+
+        print 'phi', phi
+        self.speed = check_forcefield(s)
+        self.phi = check_forcefield(phi)
+
+        self.const = eta_w*rho_a/rho_w
+
+
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import Float, ones, ArrayType
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        #ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        N = domain.number_of_elements
+        t = domain.time
+
+        if callable(self.speed):
+            xc = domain.get_centroid_coordinates()
+            #s_vec = self.speed(t, xc[:,0], xc[:,1])
+            s_vec = self.speed(t, xc)
+        else:
+            #Assume s is a scalar
+
+            try:
+                s_vec = self.speed * ones(N, Float)
+            except:
+                msg = 'Speed must be either callable or a scalar: %s' %self.s
+                raise msg
+
+
+        if callable(self.phi):
+            xc = domain.get_centroid_coordinates()
+            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
+            phi_vec = self.phi(t, xc)
+        else:
+            #Assume phi is a scalar
+
+            try:
+                phi_vec = self.phi * ones(N, Float)
+            except:
+                msg = 'Angle must be either callable or a scalar: %s' %self.phi
+                raise msg
+
+        #assign_windfield_values(xmom_update, ymom_update,
+        #                        s_vec, phi_vec, self.const)
+        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
+
+
+#def assign_windfield_values(xmom_update, ymom_update,
+#                            s_vec, phi_vec, const):
+def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
+    """Python version of assigning wind field to update vectors.
+    A c version also exists (for speed)
+    """
+    from math import pi, cos, sin, sqrt
+
+    N = len(s_vec)
+    for k in range(N):
+        s = s_vec[k]
+        phi = phi_vec[k]
+
+        #Convert to radians
+        phi = phi*pi/180
+
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress
+        #S = const * sqrt(u**2 + v**2)
+        S = const * u
+        xmom_update[k] += S*u
+        #ymom_update[k] += S*v
Index: /anuga_work/development/anuga_1d/steady_flow.py
===================================================================
--- /anuga_work/development/anuga_1d/steady_flow.py	(revision 5826)
+++ /anuga_work/development/anuga_1d/steady_flow.py	(revision 5827)
@@ -1,6 +1,6 @@
 import os
-import Gnuplot
+#import Gnuplot
 from math import sqrt, pi
-from shallow_water_1d import *
+from shallow_water_domain import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon
@@ -21,5 +21,5 @@
 h_1 = 0.5
 u_1 = 0.6
-g = 9.81
+
 
 #h_c = (q**2/g)**(1/3)
@@ -103,5 +103,5 @@
 
 L = 50.0
-N = 50 #800
+N = 400
 cell_len = L/N
 points = zeros(N+1,Float)
@@ -134,46 +134,59 @@
 wc, uc, hc = analytical_sol(C)
 
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
-hold(False)
-rc('text', usetex=True)
+#from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
+#hold(False)
+#rc('text', usetex=True)
 h_error = zeros(1,Float)
 uh_error = zeros(1,Float)
 k = 0
-yieldstep = 25.0
-finaltime = 25.0
+yieldstep = 1.0
+finaltime = 1.0
 domain.limiter = "vanleer"
-#domain.limiter = "steve_minmod"
-domain.default_order = 2
-domain.default_time_order = 2#check order is not working
+domain.order = 2
+domain.set_timestepping_method('rk2')
+domain.cfl = 1.0
+print 'THE DOMAIN LIMITER is', domain.limiter
+import time
+t0=time.time()
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
-print 'Test 1'    
-plot1 = subplot(211)
-print 'Test 2'
-plot(X,w,X,StageQ,X,ElevationQ)
-print 'Test 3'
-#xlabel('Position')
-#ylabel('Stage (m)')
-#legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
-#       'upper right', shadow=True)
-plot1.set_ylim([-0.1,1.0])
-print 'Test 4'
-plot2 = subplot(212)
-print 'Test 5'
-plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
-print 'Test 6'
-#xlabel('x (m)')
-#ylabel(r'X-momentum ($m^2/s$)')
-h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
-uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
-print "h_error %.10f" %(h_error[k])
-print "uh_error %.10f"% (uh_error[k])
-print "h_max %.10f"%max(wc-StageC)
-print "uh max  %.10f"%max(uc*hc-XmomC)
-filename = "steady_flow_"
-filename += domain.limiter
-filename += str(400)
-filename += ".ps"
-#savefig(filename)
-show()
-
+    print "integral", domain.quantities['stage'].get_integral()
+    if t>0.0:
+        print 'Test 1'
+        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
+        #ion()
+        hold(False)
+        #rc('text', usetex=True)
+        clf()
+        
+        plot1 = subplot(211)
+        print 'Test 2'
+        plot(X,w,X,StageQ,X,ElevationQ)
+        print 'Test 3'
+        xlabel('Position')
+        ylabel('Stage (m)')
+        legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
+               'upper right', shadow=True)
+        plot1.set_ylim([-0.1,1.0])
+        print 'Test 4'
+        plot2 = subplot(212)
+        print 'Test 5'
+        plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
+        print 'Test 6'
+        xlabel('x (m)')
+        ylabel(r'X-momentum (m^2/s)')
+        plot2.set_ylim([0.297,0.301])
+        h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
+        uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
+        print "h_error %.10f" %(h_error[k])
+        print "uh_error %.10f"% (uh_error[k])
+        print "h_max %.10f"%max(wc-StageC)
+        print "uh max  %.10f"%max(uc*hc-XmomC)
+        #filename = "steady_flow_"
+        #filename += domain.limiter
+        #filename += str(400)
+        #filename += ".ps"
+        #savefig(filename)
+        show()
+print 'That took %.2f seconds'%(time.time()-t0)
+raw_input("Press return key!")
