Index: /anuga_work/development/anuga_1d/dry_dam_vel.py
===================================================================
--- /anuga_work/development/anuga_1d/dry_dam_vel.py	(revision 5845)
+++ /anuga_work/development/anuga_1d/dry_dam_vel.py	(revision 5845)
@@ -0,0 +1,147 @@
+import os
+from math import sqrt, pi
+from shallow_water_vel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+h1 = 10.0
+h0 = 0.0
+
+def analytical_sol(C,t):
+    
+    #t  = 0.0     # time (s)
+    # gravity (m/s^2)
+    #h1 = 10.0    # depth upstream (m)
+    #h0 = 0.0     # depth downstream (m)
+    L = 2000.0   # length of stream/domain (m)
+    n = len(C)    # number of cells
+
+    u = zeros(n,Float)
+    h = zeros(n,Float)
+    x = C-3*L/4.0
+    
+
+    for i in range(n):
+        # Calculate Analytical Solution at time t > 0
+        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
+        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
+        u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
+        h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1))
+
+        if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
+            u[i] = 0.0
+            h[i] = h0
+        elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
+            u[i] = u3_
+            h[i] = h3_
+
+        elif ( x[i] <= -t*sqrt(g*h1) ):
+            u[i] = 0.0 
+            h[i] = h1 
+        elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
+            u[i] = u3 
+            h[i] = h3 
+        else:
+            u[i] = 0.0 
+            h[i] = h0 
+
+    return h , u*h, u
+
+print "TEST 1D-SOLUTION III -- DRY BED"
+
+def stage(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<=L/4.0:
+            y[i] = h0
+        elif x[i]<=3*L/4.0:
+            y[i] = h1
+        else:
+            y[i] = h0
+    return y
+
+
+import time
+
+finaltime = 10.0
+yieldstep = finaltime
+L = 2000.0     # Length of channel (m)
+number_of_cells = [800]#,200,500,1000,2000,5000,10000,20000]
+h_error = zeros(len(number_of_cells),Float)
+uh_error = zeros(len(number_of_cells),Float)
+u_error = zeros(len(number_of_cells),Float)
+k = 0
+for i in range(len(number_of_cells)):
+    N = int(number_of_cells[i])
+    print "Evaluating domain with %d cells" %N
+    cell_len = L/N # Origin = 0.0
+    points = zeros(N+1,Float)
+    for j in range(N+1):
+        points[j] = j*cell_len
+        
+    domain = Domain(points)
+    
+    domain.set_quantity('stage', stage)
+    domain.set_boundary({'exterior': Reflective_boundary(domain)})
+    domain.order = 2
+    domain.set_timestepping_method('rk2')
+    domain.set_CFL(1.0)
+    domain.set_limiter("vanleer")
+    #domain.h0=0.0001
+
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+        domain.write_time()
+
+    N = float(N)
+    StageC = domain.quantities['stage'].centroid_values
+    XmomC = domain.quantities['xmomentum'].centroid_values
+    u_C = domain.quantities['velocity'].centroid_values
+    C = domain.centroids
+    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))
+    u_error[k] = 1.0/(N)*sum(abs(u-u_C))
+    print "h_error %.10f" %(h_error[k])
+    print "uh_error %.10f"% (uh_error[k])
+    print "u_error %.10f"% (u_error[k])
+    k = k+1
+    print 'That took %.2f seconds' %(time.time()-t0)
+    X = domain.vertices
+    StageQ = domain.quantities['stage'].vertex_values
+    XmomQ = domain.quantities['xmomentum'].vertex_values
+    velQ = domain.quantities['velocity'].vertex_values
+    
+    h, uh, u = analytical_sol(X.flat,domain.time)
+    x = X.flat
+    
+    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    print 'test 2'
+    hold(False)
+    print 'test 3'
+    plot1 = subplot(211)
+    print 'test 4'
+    plot(x,h,x,StageQ.flat)
+    print 'test 5'
+    plot1.set_ylim([-1,11])
+    xlabel('Position')
+    ylabel('Stage')
+    legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+    plot2 = subplot(212)
+    plot(x,u,x,velQ.flat)
+    plot2.set_ylim([-35,35])
+    
+    xlabel('Position')
+    ylabel('Velocity')
+   
+    file = "dry_bed_"
+    file += str(number_of_cells[i])
+    file += ".eps"
+    #savefig(file)
+    show()
+    
+print "Error in height", h_error
+print "Error in xmom", uh_error
Index: /anuga_work/development/anuga_1d/parabolic_cannal_vel.py
===================================================================
--- /anuga_work/development/anuga_1d/parabolic_cannal_vel.py	(revision 5845)
+++ /anuga_work/development/anuga_1d/parabolic_cannal_vel.py	(revision 5845)
@@ -0,0 +1,144 @@
+import os
+from math import sqrt, pi, sin, cos
+from shallow_water_vel_domain 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
+
+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 = Domain(points)
+domain.order = 2
+domain.set_timestepping_method('euler')
+domain.set_CFL(1.0)
+domain.beta = 1.0
+domain.set_limiter("vanleer")
+
+    
+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 = 10.0 #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:
+        
+        x = X.flat
+        print 'domain.time=',domain.time
+        
+        w_V = domain.quantities['stage'].vertex_values.flat
+        uh_V = domain.quantities['xmomentum'].vertex_values.flat
+        u_V = domain.quantities['velocity'].vertex_values.flat
+        h_V = domain.quantities['height'].vertex_values.flat
+        
+        u,h,z,z_b,T = analytic_cannal(x,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)
+
+        #print 'size X',X.shape
+        #print 'size w',w.shape
+        signh = (h_V>0)
+
+        plot1 = subplot(311)
+        #plot(x,w,  x,w_V,  x,z_b)
+        plot(x,signh)
+        plot1.set_xlim([-6000,6000])
+        xlabel('Position')
+        ylabel('Stage')
+        legend(('Analytic Solution', 'Numerical Solution', 'Bed'),
+               'upper center', shadow=True)
+
+        
+        plot2 = subplot(312)
+        plot(x,u*h,x,uh_V)
+        #plot2.set_ylim([-1,25])
+        xlabel('Position')
+        ylabel('Xmomentum')
+
+        print u_V
+        
+        plot3 = subplot(313)
+        plot(x,u, x,u_V)
+        #plot2.set_ylim([-1,25])
+        xlabel('Position')
+        ylabel('Velocity')
+        show()
+raw_input("Press the return key!")
Index: /anuga_work/development/anuga_1d/shallow_water_vel_domain.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_vel_domain.py	(revision 5845)
+++ /anuga_work/development/anuga_1d/shallow_water_vel_domain.py	(revision 5845)
@@ -0,0 +1,615 @@
+"""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']
+        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
+        other_quantities = ['friction']
+        Generic_Domain.__init__(self,
+                                coordinates = coordinates,
+                                boundary = boundary,
+                                conserved_quantities = conserved_quantities,
+                                evolved_quantities = evolved_quantities,
+                                other_quantities = other_quantities,
+                                tagged_elements = 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"
+
+        
+        #Stored output
+        self.store = True
+        self.format = 'sww'
+        self.smooth = True
+
+        
+        #Reduction operation for get_vertex_values
+        from util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.set_quantities_to_be_stored(['stage','xmomentum'])
+
+	self.__doc__ = 'shallow_water_domain'
+
+        self.check_integrity()
+
+
+    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 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
+
+        msg = 'First evolved quantity must be "stage"'
+        assert self.evolved_quantities[0] == 'stage', msg
+        msg = 'Second evolved quantity must be "xmomentum"'
+        assert self.evolved_quantities[1] == 'xmomentum', msg
+        msg = 'Third evolved quantity must be "elevation"'
+        assert self.evolved_quantities[2] == 'elevation', msg
+        msg = 'Fourth evolved quantity must be "height"'
+        assert self.evolved_quantities[3] == 'height', msg
+        msg = 'Fifth evolved quantity must be "velocity"'
+        assert self.evolved_quantities[4] == 'velocity', msg
+
+        Generic_Domain.check_integrity(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_vel(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges_limit_w_u(self)
+        
+
+#=============== End of Shallow Water Domain ===============================
+#-----------------------------------
+# Compute fluxes interface
+#-----------------------------------
+def compute_fluxes(domain):
+    """
+    Python version of compute fluxes (local_compute_fluxes)
+    is available in test_shallow_water_vel_domain.py
+    """
+ 
+
+    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_vel_ext import compute_fluxes_ext
+
+    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)
+
+
+#-----------------------------------
+# Compute flux definition with vel
+#-----------------------------------
+def compute_fluxes_vel(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    stage    = domain.quantities['stage']
+    xmom     = domain.quantities['xmomentum']
+    bed      = domain.quantities['elevation']
+    height   = domain.quantities['height']
+    velocity = domain.quantities['velocity']
+
+
+    from comp_flux_vel_ext import compute_fluxes_vel_ext
+
+    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)
+
+
+#--------------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_w_u(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] <= 1.0e-12:
+	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
+	    h_C[i] = 0.0
+	    w_C[i] = z_C[i]
+	    #uh_C[i] = 0.0
+            
+ #           u_C[i] = 0.0
+ #       else:
+ #           u_C[i] = uh_C[i]/h_C[i]
+		
+    h0 = 1.0e-12    
+    for i in range(len(h_C)):
+	if h_C[i] < 1.0e-12:
+	    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])
+            #u_C[i]  = uh_C[i]/h_C[i]
+	
+    for name in [ 'velocity', 'stage' ]:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    w_V  = domain.quantities['stage'].vertex_values			
+    z_V  = domain.quantities['elevation'].vertex_values	
+    h_V  = domain.quantities['height'].vertex_values
+    u_V  = domain.quantities['velocity'].vertex_values		
+    uh_V = domain.quantities['xmomentum'].vertex_values	
+
+    h_V[:]    = w_V - z_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
+                w_V[i,j] = z_V[i,j]
+                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
+                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
+                
+    uh_V[:] = u_V * h_V
+
+    
+    return
+
+#---------------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_w_uh(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave equation.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    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
+	
+
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]						
+        if h_C[i] <= 1.0e-6:
+	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
+	    h_C[i] = 0.0
+	    w_C[i] = z_C[i]
+	    uh_C[i] = 0.0
+            
+    for name in [ 'stage', 'xmomentum']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    w_V  = domain.quantities['stage'].vertex_values			
+    z_V  = domain.quantities['elevation'].vertex_values	
+    h_V  = domain.quantities['height'].vertex_values
+    u_V  = domain.quantities['velocity'].vertex_values		
+    uh_V = domain.quantities['xmomentum'].vertex_values	
+
+    h_V[:]    = w_V - z_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
+                w_V[i,j] = z_V[i,j]
+                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
+                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
+                u_V[i,j] = 0.0
+            if h_V[i,j] < h0:
+                u_V[i,j]
+            else:
+                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
+
+
+#--------------------------------------------------------
+#Boundaries - specific to the shallow_water_vel_domain
+#--------------------------------------------------------
+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.normals  = domain.normals
+        self.stage    = domain.quantities['stage'].vertex_values
+        self.xmom     = domain.quantities['xmomentum'].vertex_values
+        self.bed      = domain.quantities['elevation'].vertex_values
+        self.height   = domain.quantities['height'].vertex_values
+        self.velocity = domain.quantities['velocity'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.evolved_quantities = zeros(5, 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.evolved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = -self.xmom[vol_id, edge_id]
+        q[2] = self.bed[vol_id, edge_id]
+        q[3] = self.height[vol_id, edge_id]
+        q[4] = -self.velocity[vol_id, edge_id]
+
+        #print "In Reflective q ",q
+
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, evolved_quantities=None):
+        Boundary.__init__(self)
+
+        if evolved_quantities is None:
+            msg = 'Must specify one value for each evolved quantity'
+            raise msg
+
+        from Numeric import array, Float
+        self.evolved_quantities=array(evolved_quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.evolved_quantities
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.evolved_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
+
+
+
+    Stage     = domain.quantities['stage']
+    Xmom      = domain.quantities['xmomentum']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+
+    xmom_ud  = Xmom.explicit_update
+    #stage_ud = Stage.explicit_update
+
+
+    #h = Stage.vertex_values - Elevation.vertex_values
+    h = Height.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 = 0.5*(h[k,0] + h[k,1])
+
+        #Compute bed slope
+        x0, x1 = x[k,:]
+        b0, b1 = b[k,:]
+        bx = gradient(x0, x1, b0, b1)
+        
+        #Update momentum (explicit update is reset to source values)
+        xmom_ud[k] += -g*bx*avg_h
+        #stage_ud[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
+    tau = domain.quantities['linear_friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+
+    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]
+
+
+
+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
+
Index: /anuga_work/development/anuga_1d/test_shallow_water_vel_domain.py
===================================================================
--- /anuga_work/development/anuga_1d/test_shallow_water_vel_domain.py	(revision 5845)
+++ /anuga_work/development/anuga_1d/test_shallow_water_vel_domain.py	(revision 5845)
@@ -0,0 +1,517 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from shallow_water_vel_domain import *
+
+
+from Numeric import allclose, array, ones, Float, maximum, zeros
+
+
+class Test_Shallow_Water(unittest.TestCase):
+    def setUp(self):
+        self.points = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        domain = Domain(self.points)
+        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
+
+    def test_reflective_boundary(self):
+        """
+        Test setting reflective boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test reflective'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [-6.0, -6.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [0.0, 0.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [-3.0, -3.0])
+
+        
+
+    def test_dirichlet_boundary(self):
+        """
+        Test setting dirichlet boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Dirichlet_boundary([3.0, 8.0, 1.0, 2.0, 4.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test dirichlet'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [3.0, 3.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [8.0, 8.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [1.0, 1.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [4.0, 4.0])
+
+
+    def test_compute_fluxes(self):
+        """
+        Compare shallow_water_domain flux calculation against a previous
+        Python implementation (defined in this file)
+        """
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_boundary({'exterior' : Dirichlet_boundary([0.0, 0.0, 0.0, 0.0, 0.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        stage_ud, xmom_ud = local_compute_fluxes(domain)
+
+        domain.compute_fluxes()
+        
+        print domain.quantities['stage'].vertex_values
+        print domain.quantities['xmomentum'].vertex_values
+        print domain.quantities['elevation'].vertex_values
+        print domain.quantities['height'].vertex_values
+        print domain.quantities['velocity'].vertex_values
+
+        print domain.quantities['stage'].boundary_values
+        print domain.quantities['xmomentum'].boundary_values
+        print domain.quantities['elevation'].boundary_values
+        print domain.quantities['height'].boundary_values
+        print domain.quantities['velocity'].boundary_values
+        
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
+
+        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
+        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        #This assumes h0 = 1.0e-3!!
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        #print maxspeed
+        #print edgeflux
+        
+        assert allclose(array([2.0, 8.9],Float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+
+        #print maxspeed
+        #print edgeflux        
+        
+        assert allclose(array([-2.0, -8.9],Float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+
+        def slope_one(x):
+            return x
+
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
+
+
+##     def test_evolve_first_order(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+        
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 1
+##         domain.set_timestepping_method('euler')
+        
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         ##	print domain.quantities['stage'].vertex_values 
+##         ##	print domain.quantities['elevation'].vertex_values 
+##         ##        print domain.quantities['xmomentum'].vertex_values
+##         ##  
+##         ##
+##         ##        print domain.quantities['stage'].centroid_values  
+##         ##	print domain.quantities['elevation'].centroid_values    
+##         ##        print domain.quantities['xmomentum'].centroid_values    
+
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+##     def test_evolve_euler_second_order_space(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk2')
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+##     def test_evolve_second_order_space_time(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk3')
+
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+
+#==============================================================================
+
+def local_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
+
+    tmp0 = zeros(N,Float)
+    tmp1 = zeros(N,Float)
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    #    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
+    #   ymom =  Ymom.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    #print 'stage edge values', stage
+    #print 'xmom edge values', xmom
+    #print 'bed values', bed
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    #print 'stage_bdry',stage_bdry
+    #print 'xmom_bdry', xmom_bdry
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, Float) #Work array for summing up fluxes
+    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]
+                #print i, ' ' , m
+                #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 'ql',ql
+            #print 'qr',qr
+            
+
+            edgeflux, max_speed = local_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]
+        tmp0[k] = flux[0]
+        tmp1[k] = flux[1]
+
+
+    return tmp0, tmp1
+
+
+def local_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
+    
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    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**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    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
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')    
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
+
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
+    #print "restricted test"
+    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
