Index: /anuga_work/development/anuga_1d/backup.py
===================================================================
--- /anuga_work/development/anuga_1d/backup.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/backup.py	(revision 6453)
@@ -0,0 +1,537 @@
+"""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
+
+This particular modification of the Domain class implements the ability to
+vary the width of the 1D channel that the water flows in. As a result the
+conserved variables are different than previous implementations and so are the
+equations.
+
+U_t + E_x = S
+
+where
+------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
+U = [A, Q]
+E = [Q, Q^2/A + 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,
+Padarn Wilson, Geoscience Australia, 2008
+"""
+
+
+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 = ['area', 'discharge']
+        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
+        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
+
+     
+        self.forcing_terms.append(gravity)
+        self.forcing_terms.append(boundary_stress)
+        #self.forcing_terms.append(manning_friction)
+       
+
+        
+        #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(['area','discharge'])
+
+	self.__doc__ = 'channel_domain_Ab'
+
+        self.check_integrity()
+
+
+    def check_integrity(self):
+
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "area"'
+        assert self.conserved_quantities[0] == 'area', msg
+        msg = 'Second conserved quantity must be "discharge"'
+        assert self.conserved_quantities[1] == 'discharge', msg
+
+        msg = 'First evolved quantity must be "area"'
+        assert self.evolved_quantities[0] == 'area', msg
+        msg = 'Second evolved quantity must be "discharge"'
+        assert self.evolved_quantities[1] == 'discharge', 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
+        msg = 'Fifth evolved quantity must be "width"'
+        assert self.evolved_quantities[5] == 'width', msg
+
+        Generic_Domain.check_integrity(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_channel(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_a_d(self)
+        
+
+#=============== End of Channel Domain ===============================
+
+#-----------------------------------
+# Compute flux definition with channel
+#-----------------------------------
+def compute_fluxes_channel(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    area    = domain.quantities['area']
+    discharge     = domain.quantities['discharge']
+    bed      = domain.quantities['elevation']
+    height   = domain.quantities['height']
+    velocity = domain.quantities['velocity']
+    width    = domain.quantities['width']
+
+    
+ 
+    from channel_domain_ext import compute_fluxes_channel_ext
+    domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
+
+#-----------------------------------------------------------------------
+# Distribute to verticies with stage reconstructed and then extrapolated
+#-----------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_a_d(domain):
+    
+    #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
+    Area = domain.quantities['area']
+    Discharge = domain.quantities['discharge']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+    Width = domain.quantities['width']
+    Stage = domain.quantities['stage']
+
+    #Arrays   
+    a_C   = Area.centroid_values    
+    d_C  = Discharge.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+    b_C   = Width.centroid_values
+    w_C   = Stage.centroid_values
+
+    h0 = 1.0e-12	
+    #print id(h_C)
+    for i in range(N):
+							
+        if a_C[i] <= h0:
+	    a_C[i] = 0.0
+	    h_C[i] = 0.0
+            d_C[i] = 0.0
+            u_C[i] = 0.0
+            w_C[i] = z_C[i]
+	   
+            
+
+        else:
+
+            if b_C[i]<=h0:
+                a_C[i] = 0.0
+                h_C[i] = 0.0
+                d_C[i] = 0.0
+                u_C[i] = 0.0
+                w_C[i] = z_C[i]
+
+            else:
+                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
+                w_C[i] = h_C[i]+z_C[i]
+                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
+		
+   
+	
+    for name in ['discharge','elevation','width','stage','height','area','velocity']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+    a_V  = domain.quantities['area'].vertex_values
+    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		
+    d_V = domain.quantities['discharge'].vertex_values
+    b_V = domain.quantities['width'].vertex_values
+
+    
+##     for i in range(len(h_C)):
+##         for j in range(2):
+##             if b_V[i,j] < h0 :
+##                 a_V[i,j]=0
+##                 h_V[i,j]=0
+##                 d_V[i,j]=0
+##                 u_V[i,j]=0
+##             else:                    
+##                 h_V[i,j] = w_V[i,j]-z_V[i,j]
+##                 if h_V[i,j]<h0:
+##                     h_V[i,j]=0
+##                     w_V[i,j]=z_V[i,j]
+##                 a_V[i,j] = b_V[i,j]*h_V[i,j]
+##                 if d_V[i,j]<h0 or a_V[i,j]<h0:
+##                     d_V[i,j]=0
+##                     u_V[i,j]=0
+##                 else:
+##                     u_V[i,j]=d_V[i,j]/(a_V[i,j]+h0/a_V[i,j])
+            
+                
+                
+    
+
+    
+    return
+
+
+#--------------------------------------------------------
+#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.area    = domain.quantities['area'].vertex_values
+        self.discharge     = domain.quantities['discharge'].vertex_values
+        self.bed      = domain.quantities['elevation'].vertex_values
+        self.height   = domain.quantities['height'].vertex_values
+        self.velocity = domain.quantities['velocity'].vertex_values
+        self.width    = domain.quantities['width'].vertex_values
+        self.stage    = domain.quantities['stage'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.evolved_quantities = zeros(7, 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.area[vol_id, edge_id]
+        q[1] = -self.discharge[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]
+        q[5] = self.width[vol_id,edge_id]
+        q[6] = self.stage[vol_id,edge_id]
+
+        #print "In Reflective q ",q
+
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+if k>5 and k<15:
+            print discharge_ud[k],-g*zx*avg_h*avg_b
+        discharge_ud[k] +=-g*zx*avg_h*avg_b    """
+
+
+    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
+
+
+
+    Area     = domain.quantities['area']
+    Discharge  = domain.quantities['discharge']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+    Width     = domain.quantities['width']
+
+    discharge_ud  = Discharge.explicit_update
+
+
+       
+    h = Height.vertex_values
+    b = Width.vertex_values
+    a = Area.vertex_values
+    z = Elevation.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])
+        avg_b = 0.5*(b[k,0] + b[k,1])
+        
+        #Compute bed slope
+        x0, x1 = x[k,:]
+        z0, z1 = z[k,:]
+        zx = gradient(x0, x1, z0, z1)
+        #Update momentum (explicit update is reset to source values)
+        discharge_ud[k]+= -g*zx*avg_h*avg_b
+        
+      
+
+def boundary_stress(domain):
+
+
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+
+
+
+    Area     = domain.quantities['area']
+    Discharge  = domain.quantities['discharge']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+    Width     = domain.quantities['width']
+
+    discharge_ud  = Discharge.explicit_update
+
+
+       
+    h = Height.vertex_values
+    b = Width.vertex_values
+    a = Area.vertex_values
+    z = Elevation.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)
+        discharge_ud[k] += 0.5*g*bx*avg_h*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/channel_domain.py
===================================================================
--- /anuga_work/development/anuga_1d/channel_domain.py	(revision 6452)
+++ /anuga_work/development/anuga_1d/channel_domain.py	(revision 6453)
@@ -6,11 +6,15 @@
 consisting of methods specific to the Shallow Water Wave Equation
 
+This particular modification of the Domain class implements the ability to
+vary the width of the 1D channel that the water flows in. As a result the
+conserved variables are different than previous implementations and so are the
+equations.
 
 U_t + E_x = S
 
 where
-
-U = [w, uh]
-E = [uh, u^2h + gh^2/2]
+------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
+U = [A, Q]
+E = [Q, Q^2/A + gh^2/2]
 S represents source terms forcing the system
 (e.g. gravity, friction, wind stress, ...)
@@ -32,5 +36,5 @@
 
 The conserved quantities are w, uh
-
+--------------------------------------------------------------------------
 For details see e.g.
 Christopher Zoppou and Stephen Roberts,
@@ -39,6 +43,6 @@
 
 
-John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-Geoscience Australia, 2006
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
+Padarn Wilson, Geoscience Australia, 2008
 """
 
Index: /anuga_work/development/anuga_1d/channel_domain2.py
===================================================================
--- /anuga_work/development/anuga_1d/channel_domain2.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/channel_domain2.py	(revision 6453)
@@ -0,0 +1,852 @@
+"""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 = ['area', 'discharge']
+        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'stage']
+        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(['area','discharge'])
+
+	self.__doc__ = 'shallow_water_domain'
+
+       # self.check_integrity()
+
+
+    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_channel(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_a_d(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)
+    domain.flux_timestep = .1
+
+#-----------------------------------
+# 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)
+
+#----------------------------------
+#  Compute fluxes channel
+#----------------------------------
+def compute_fluxes_channel(domain):
+    from Numeric import zeros, Float
+    import sys
+    
+    timestep = float(sys.maxint)
+
+    area      = domain.quantities['area']
+    discharge = domain.quantities['discharge']
+    bed       = domain.quantities['elevation']
+    height     = domain.quantities['height']
+    stage     = domain.quantities['stage']
+
+    #from channel_domain_ext import compute_fluxes_channel_ext
+
+   # domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,stage)
+   # domain.quantities['area'].explicit_update=ones(410,Float)
+    domain.flux_timestep = .1
+    #print area.vertex_values[0],area.vertex_values[1],area.vertex_values[2]
+    for i in range(len(domain.coordinates)-1):
+        fluxesl= channel_flux_func(domain,i-1)
+        fluxesr= channel_flux_func(domain,i)
+        #print domain.areas[i]
+        #print -1*(fluxesl[1]-fluxesr[1])
+        #print height.centroid_values[200],height.vertex_values[200]
+        #print fluxesl[0]
+        print -1*(fluxesl[0]-fluxesr[0])*domain.areas[i]
+        area.explicit_update[i]=-1*(fluxesl[0]-fluxesr[0])*domain.areas[i]
+        discharge.explicit_update[i]=-(fluxesl[1]-fluxesr[1])*domain.areas[i]
+        
+def channel_flux_func(domain, i):
+
+    area      = domain.quantities['area']
+    discharge = domain.quantities['discharge']
+    bed       = domain.quantities['elevation']
+    height     = domain.quantities['height']
+    stage     = domain.quantities['stage']
+
+    flux_left0=0
+    flux_left1=0
+    flux_right0=0
+    flux_right1=0
+    flux0=0
+    flux1=0
+  # Crude numerical flux calculation
+    from math import sqrt
+
+    if i==0 or i==(len(domain.coordinates)-1):
+        flux_left0=0
+        flux_left1=0
+        flux_right0=0
+        flux_right1=0
+    else:
+        g=9.8
+        a_left = area.centroid_values[i-1]
+        d_left = discharge.centroid_values[i-1]
+        z_left = bed.centroid_values[i-1]
+        h_left = height.centroid_values[i-1]
+        w_left = stage.centroid_values[i-1]
+
+        a_right = area.centroid_values[i]
+        d_right = discharge.centroid_values[i]
+        z_right = bed.centroid_values[i]
+        h_right = height.centroid_values[i]
+        w_right = stage.centroid_values[i]
+
+        z=(z_left+z_right)/2
+        #hbarr=0.5*(height.vertex_values[i][1]+height.vertex_values[i][0])
+        #hbarl=0.5*(height.vertex_values[i-1][1]+height.vertex_values[i-1][0])
+        hbarr=0
+        hbarl=0
+        #print hbarl,h_left
+        if a_left>1.0e-12:
+            u_left=d_left/a_left
+        else:
+            u_left=0
+        if a_right>1.0e-12:
+            u_right=d_right/a_right
+        else:
+            u_right=0
+
+        ## soundspeed_left = sqrt(g*h_left);
+##         soundspeed_right = sqrt(g*h_right);
+
+##         s_max = max(u_left+soundspeed_left,u_right+soundspeed_right)
+##         if s_max<0.0:
+##             s_max=0
+##         s_min = min(u_left-soundspeed_left,u_right-soundspeed_right)
+##         if s_min>0.0:
+##             s_min=0
+            
+        flux_left0 = d_left
+        flux_right1= d_right
+        if a_left<1.0e-12:
+            flux_left11=0
+        else:
+            flux_left11=d_left*d_left/a_left
+        flux_left1=flux_left11+g*(calculateI(h_left,0)-calculateI(hbarl,0))
+        if a_right<1.0e-12:
+            flux_right11=0
+        else:
+            flux_right11=d_right*d_right/a_right
+        flux_right1=flux_right11+g*(calculateI(h_right,0)-calculateI(hbarr,0))
+        #print g*(calculateI(h_right,0)-calculateI(hbarr,0)), g*(calculateI(h_left,0)-calculateI(hbarl,0)) 
+
+    flux0 = 0.5*(flux_left0+flux_right0)
+    flux1 = 0.5*(flux_left1+flux_right1)
+    return (flux0,flux1)
+            
+def calculateI(H,z0):
+    return H*H*0.5-z0*H+z0*z0*0.5
+
+    
+#--------------------------------------------------------------------------
+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])
+
+#---------------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_a_d(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
+    Area = domain.quantities['area']
+    Discharge = domain.quantities['discharge']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Stage = domain.quantities['stage']
+
+    #Arrays   
+    a_C   = Area.centroid_values    
+    d_C   = Discharge.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    w_C   = Stage.centroid_values
+	
+#work out stage
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]
+#make sure depth isn't zero        
+        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]
+	    d_C[i] = 0.0
+#distribute stage,discharge
+
+    for name in [ 'area', 'discharge']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    a_V  = domain.quantities['area'].vertex_values			
+    d_V  = domain.quantities['discharge'].vertex_values	
+    h_V  = domain.quantities['height'].vertex_values
+    z_V  = domain.quantities['elevation'].vertex_values		
+    w_V = domain.quantities['stage'].vertex_values	
+#height at verticies
+    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.
+        """
+#Commenting out some quantities not currently keeping track of
+##         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.stage[stage_id, stage_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
+
+#--------------------------------------------------------
+#Boundaries for channel - specific to the channel 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.area    = domain.quantities['area'].vertex_values
+        self.discharge     = domain.quantities['discharge'].vertex_values
+        self.bed      = domain.quantities['elevation'].vertex_values
+        self.height   = domain.quantities['height'].vertex_values
+        self.stage = domain.quantities['stage'].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.area[vol_id, edge_id]
+        q[1] = -self.discharge[vol_id, edge_id]
+        q[2] = self.bed[vol_id, edge_id]
+        q[3] = self.height[vol_id, edge_id]
+        q[4] = self.stage[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
+
+
+
+    Area     = domain.quantities['area']
+    Discharge      = domain.quantities['discharge']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+    Stage     = domain.quantities['stage']
+
+    discharge_ud  = Discharge.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)
+        discharge_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/channel_domain_Ab.py
===================================================================
--- /anuga_work/development/anuga_1d/channel_domain_Ab.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/channel_domain_Ab.py	(revision 6453)
@@ -0,0 +1,563 @@
+"""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
+
+This particular modification of the Domain class implements the ability to
+vary the width of the 1D channel that the water flows in. As a result the
+conserved variables are different than previous implementations and so are the
+equations.
+
+U_t + E_x = S
+
+where
+------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
+U = [A, Q]
+E = [Q, Q^2/A + 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,
+Padarn Wilson, Geoscience Australia, 2008
+"""
+
+
+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 = ['area', 'discharge']
+        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
+        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
+        self.setstageflag = False
+        self.discontinousb = False
+
+     
+        #self.forcing_terms.append(gravity)
+        #self.forcing_terms.append(boundary_stress)
+        #self.forcing_terms.append(manning_friction)
+       
+
+        
+        #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(['area','discharge'])
+
+	self.__doc__ = 'channel_domain_Ab'
+
+        self.check_integrity()
+
+
+    def check_integrity(self):
+
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "area"'
+        assert self.conserved_quantities[0] == 'area', msg
+        msg = 'Second conserved quantity must be "discharge"'
+        assert self.conserved_quantities[1] == 'discharge', msg
+
+        msg = 'First evolved quantity must be "area"'
+        assert self.evolved_quantities[0] == 'area', msg
+        msg = 'Second evolved quantity must be "discharge"'
+        assert self.evolved_quantities[1] == 'discharge', 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
+        msg = 'Fifth evolved quantity must be "width"'
+        assert self.evolved_quantities[5] == 'width', msg
+
+        Generic_Domain.check_integrity(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_channel(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_a_d(self)
+        
+
+#=============== End of Channel Domain ===============================
+
+#-----------------------------------
+# Compute flux definition with channel
+#-----------------------------------
+def compute_fluxes_channel(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    area    = domain.quantities['area']
+    discharge     = domain.quantities['discharge']
+    bed      = domain.quantities['elevation']
+    height   = domain.quantities['height']
+    velocity = domain.quantities['velocity']
+    width    = domain.quantities['width']
+
+
+    from channel_domain_ext import compute_fluxes_channel_ext
+    domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
+
+#-----------------------------------------------------------------------
+# Distribute to verticies with stage reconstructed and then extrapolated
+#-----------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_a_d(domain):
+    
+    #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
+   ##  linearb(domain)
+
+
+
+    
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Area = domain.quantities['area']
+    Discharge = domain.quantities['discharge']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+    Width = domain.quantities['width']
+    Stage = domain.quantities['stage']
+
+    #Arrays   
+    a_C   = Area.centroid_values    
+    d_C  = Discharge.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+    b_C   = Width.centroid_values
+    w_C   = Stage.centroid_values
+
+    if domain.setstageflag:
+        for i in range(len(a_C)):
+            a_C[i]=(w_C[i]-z_C[i])*b_C[i]
+     
+        domain.setstageflag = False
+
+    if domain.discontinousb:
+        domain.quantities['width'].extrapolate_second_order()
+        
+
+    h0 = 1.0e-12	
+    #print id(h_C)
+    for i in range(N):
+							
+        if a_C[i] <= h0:
+	    a_C[i] = 0.0
+	    h_C[i] = 0.0
+            d_C[i] = 0.0
+            u_C[i] = 0.0
+            w_C[i] = z_C[i]
+	   
+            
+
+        else:
+
+            if b_C[i]<=h0:
+                a_C[i] = 0.0
+                h_C[i] = 0.0
+                d_C[i] = 0.0
+                u_C[i] = 0.0
+                w_C[i] = z_C[i]
+
+            else:
+                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
+                w_C[i] = h_C[i]+z_C[i]
+                u_C[i] = d_C[i]/(a_C[i]+h0/a_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'
+    a_V  = domain.quantities['area'].vertex_values
+    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		
+    d_V = domain.quantities['discharge'].vertex_values
+    b_V = domain.quantities['width'].vertex_values
+
+   
+    
+    for i in range(len(h_C)):
+        for j in range(2):
+##             if b_V[i,j] < h0 :
+##                 a_V[i,j]=0
+##                 h_V[i,j]=0
+##                 d_V[i,j]=0
+##                 u_V[i,j]=0
+##             else:
+                
+            h_V[i,j] = w_V[i,j]-z_V[i,j]
+            if h_V[i,j]<h0:
+                h_V[i,j]=0
+                w_V[i,j]=z_V[i,j]
+            a_V[i,j] = b_V[i,j]*h_V[i,j]
+            d_V[i,j]=u_V[i,j]*a_V[i,j]
+            
+                
+                
+    
+
+    
+    return
+
+
+#--------------------------------------------------------
+#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.area    = domain.quantities['area'].vertex_values
+        self.discharge     = domain.quantities['discharge'].vertex_values
+        self.bed      = domain.quantities['elevation'].vertex_values
+        self.height   = domain.quantities['height'].vertex_values
+        self.velocity = domain.quantities['velocity'].vertex_values
+        self.width    = domain.quantities['width'].vertex_values
+        self.stage    = domain.quantities['stage'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.evolved_quantities = zeros(7, 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.area[vol_id, edge_id]
+        q[1] = -self.discharge[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]
+        q[5] = self.width[vol_id,edge_id]
+        q[6] = self.stage[vol_id,edge_id]
+
+        #print "In Reflective q ",q
+
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+if k>5 and k<15:
+            print discharge_ud[k],-g*zx*avg_h*avg_b
+        discharge_ud[k] +=-g*zx*avg_h*avg_b    """
+
+
+    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
+
+
+
+    Area     = domain.quantities['area']
+    Discharge  = domain.quantities['discharge']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+    Width     = domain.quantities['width']
+
+    discharge_ud  = Discharge.explicit_update
+
+
+       
+    h = Height.vertex_values
+    b = Width.vertex_values
+    a = Area.vertex_values
+    z = Elevation.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])
+        avg_b = 0.5*(b[k,0] + b[k,1])
+        
+        #Compute bed slope
+        x0, x1 = x[k,:]
+        z0, z1 = z[k,:]
+        zx = gradient(x0, x1, z0, z1)
+       
+        #Update momentum (explicit update is reset to source values)
+        discharge_ud[k]+= -g*zx*avg_h*avg_b
+      
+
+def boundary_stress(domain):
+
+
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+
+
+
+    Area     = domain.quantities['area']
+    Discharge  = domain.quantities['discharge']
+    Elevation = domain.quantities['elevation']
+    Height    = domain.quantities['height']
+    Width     = domain.quantities['width']
+
+    discharge_ud  = Discharge.explicit_update
+
+
+       
+    h = Height.vertex_values
+    b = Width.vertex_values
+    a = Area.vertex_values
+    z = Elevation.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)
+        discharge_ud[k] += 0.5*g*bx*avg_h*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
+
+def linearb(domain):
+
+    bC = domain.quantities['width'].vertex_values
+    
+    for i in range(len(bC)-1):
+        temp= 0.5*(bC[i,1]+bC[i+1,0])
+        bC[i,1]=temp
+        bC[i+1,0]=temp
+
+
+
Index: /anuga_work/development/anuga_1d/channel_flow_1_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/channel_flow_1_padarn.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/channel_flow_1_padarn.py	(revision 6453)
@@ -0,0 +1,109 @@
+import os
+from math import sqrt, pow, pi 
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Channel Flow 1 Padarn Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    return 1.4691*width(x)
+
+def width(x):
+    x1=(x/1000)*(x/1000)
+    x2=x1*(x/1000)
+    x3=x2*(x/1000)
+    return 10-64*(x1-2*x2+x3)
+
+def bed(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<525 and x[i]>475:
+            y[i]=1
+        else:
+            y[i]=0
+    return y
+
+def initial_discharge(x):
+    return 20
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [200]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+domain.set_quantity('elevation',bed)
+domain.set_quantity('discharge',initial_discharge)
+
+# Set boundry type, order, timestepping method and limiter
+domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+domain.order = 2
+domain.set_timestepping_method('rk2')
+domain.set_CFL(1.0)
+domain.set_limiter("vanleer")
+#domain.h0=0.0001
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,z,x,stage)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+    
Index: /anuga_work/development/anuga_1d/channel_flow_2_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/channel_flow_2_padarn.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/channel_flow_2_padarn.py	(revision 6453)
@@ -0,0 +1,110 @@
+import os
+from math import sqrt, pow, pi 
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Channel Flow 1 Padarn Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    return 1.4691*width(x)
+
+def width(x):
+    x1=(x/1000)*(x/1000)
+    x2=x1*(x/1000)
+    x3=x2*(x/1000)
+    return 10-64*(x1-2*x2+x3)
+
+def bed(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<525 and x[i]>475:
+            y[i]=0
+        else:
+            y[i]=0
+    return y
+ 
+
+def initial_discharge(x):
+    return 20
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [200]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+domain.set_quantity('elevation',bed)
+domain.set_quantity('discharge',initial_discharge)
+
+# Set boundry type, order, timestepping method and limiter
+domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+domain.order = 2
+domain.set_timestepping_method('rk2')
+domain.set_CFL(1.0)
+domain.set_limiter("vanleer")
+#domain.h0=0.0001
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,z,x,stage)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+    
Index: /anuga_work/development/anuga_1d/compile.py
===================================================================
--- /anuga_work/development/anuga_1d/compile.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/compile.py	(revision 6453)
@@ -0,0 +1,440 @@
+"""compile.py - compile Python C-extension
+
+   Commandline usage: 
+     python compile.py <filename>
+
+   Usage from within Python:
+     import compile
+     compile.compile(<filename>,..)
+
+   Ole Nielsen, Duncan Gray Oct 2001      
+"""     
+
+# FIXME (Ole): Although this script says it works with a range of compilers,
+# it has only really been used with gcc.
+
+import os, string, sys, types
+
+separation_line = '---------------------------------------'      
+ 
+def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
+  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
+  
+     Compile FN(s) using compiler CC (e.g. mpicc), 
+     Loader LD and shared flag SFLAG.
+     If CC is absent use default compiler dependent on platform
+     if LD is absent CC is used.
+     if SFLAG is absent platform default is used
+     FNs can be either one filename or a list of filenames
+     In the latter case, the first will be used to name so file.
+  """
+  
+  # Input check
+  #
+  assert not FNs is None, 'No filename provided'
+
+  if not type(FNs) == types.ListType:
+    FNs = [FNs]
+
+
+  libext = 'so' #Default extension (Unix)
+  libs = ''
+  version = sys.version[:3]
+  
+  # Determine platform and compiler
+  #
+  if sys.platform == 'sunos5':  #Solaris
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+      
+  elif sys.platform == 'osf1V5':  #Compaq AlphaServer
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'linux2':  #Linux
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'darwin':  #Mac OS X:
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'bundle -flat_namespace -undefined suppress'
+
+  elif sys.platform == 'cygwin':  #Cygwin (compilation same as linux)
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+    libext = 'dll'
+    libs = '/lib/python%s/config/libpython%s.dll.a' %(version,version)
+      
+  elif sys.platform == 'win32':  #Windows
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc.exe' #Some systems require this (a security measure?) 
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+     
+    # As of python2.5, .pyd is the extension for python extension
+    # modules.
+    if sys.version_info[0:2] >= (2, 5):
+      libext = 'pyd'
+    else:
+      libext = 'dll'
+
+    libs, is_found = set_python_dll_path()
+      
+  else:
+    if verbose: print "Unrecognised platform %s - revert to default"\
+                %sys.platform
+    
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = 'ld'
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+
+
+  # Verify that compiler can be executed
+  print 'Compiler: %s, version ' %compiler,
+  sys.stdout.flush()
+  s = '%s -dumpversion' %(compiler)
+  err = os.system(s)
+  print
+  
+  if err != 0:
+      msg = 'Unable to execute compiler: %s. ' %compiler
+      msg += 'Make sure it is available on the system path.\n'
+      msg += 'One way to check this is to run %s on ' %compiler
+      msg += 'the commandline.'
+      raise Exception, msg    
+
+  
+  # Find location of include files
+  #
+  if sys.platform == 'win32':  #Windows
+    python_include = os.path.join(sys.exec_prefix, 'include')    
+  else:  
+    python_include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
+                                  'python' + version)
+
+  # Check existence of Python.h
+  #
+  headerfile = python_include + os.sep + 'Python.h'
+  try:
+    open(headerfile, 'r')
+  except:
+    raise """Did not find Python header file %s.
+    Make sure files for Python C-extensions are installed. 
+    In debian linux, for example, you need to install a
+    package called something like python2.3-dev""" %headerfile
+
+
+
+  # Add Python path + utilities to includelist (see ticket:31)
+  # Assume there is only one 'utilities' dir under path dirs
+  
+  utilities_include_dir = None
+  for pathdir in sys.path:
+
+    utilities_include_dir = pathdir + os.sep + 'utilities'
+    #print pathdir
+    #print utilities_include_dir
+    try:
+      os.stat(utilities_include_dir)
+    except OSError:
+      pass
+    else:
+      #print 'Found %s to be used as include dir' %utilities_include_dir
+      break
+
+  # This is hacky since it
+  # assumes the location of the compile_all that determines buildroot
+  try:
+    utilities_include_dir = buildroot + os.sep + "source" + os.sep + "anuga" \
+                            + os.sep + 'utilities'
+  except:
+    # This will make compile work locally
+    utilities_include_dir = '.'
+    utilities_include_dir = '../utilities'
+
+
+    
+  try:
+    os.stat(utilities_include_dir)
+  except OSError: 
+    utilities_include_dir = buildroot + os.sep + 'utilities'
+  
+  
+  
+  # Check filename(s)
+  #
+  object_files = ''
+  for FN in FNs:        
+    root, ext = os.path.splitext(FN)
+    if ext == '':
+      FN = FN + '.c'
+    elif ext.lower() != '.c':
+      raise Exception, "Unrecognised extension: " + FN
+    
+    try:
+      open(FN, 'r')
+    except:
+      raise Exception, "Could not open: " + FN
+
+    if not object_files: root1 = root  #Remember first filename        
+    object_files += root + '.o '  
+  
+  
+    # Compile
+    #
+    if utilities_include_dir is None:    
+      s = '%s -c %s -I"%s" -o "%s.o" -Wall -O3'\
+          %(compiler, FN, python_include, root)
+    else:
+      if FN == "triangle.c" or FN == "mesh_engine_c_layer.c":
+        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -O3 -DTRILIBRARY=1 -DNO_TIMER=1'\
+            %(compiler, FN, python_include, utilities_include_dir, root)
+      else:
+        s = '%s -c %s -I"%s" -I"%s" -o "%s.o" -Wall -O3'\
+            %(compiler, FN, python_include, utilities_include_dir, root)
+
+    if os.name == 'posix' and os.uname()[4] == 'x86_64':
+      #Extra flags for 64 bit architectures
+      #Second clause will always fail on Win32 because uname is UNIX specific
+      #but won't get past first clause
+
+      #FIXME: Which one?
+      #s += ' -fPIC'
+      s += ' -fPIC -m64' 
+      
+      
+    if verbose:
+      print s
+
+    # Doesn't work on Windows anyway  
+    #else:
+    #  s = s + ' 2> /dev/null' #Suppress errors
+    
+    try:
+      err = os.system(s)
+      if err != 0:
+          raise 'Attempting to compile %s failed - please try manually' %FN 
+    except:
+      raise 'Could not compile %s - please try manually' %FN  
+
+  
+  # Make shared library (*.so or *.dll)
+  if libs is "":
+    s = '%s -%s %s -o %s.%s -lm' %(loader, sharedflag, object_files, root1, libext)
+  else:
+    s = '%s -%s %s -o %s.%s "%s" -lm' %(loader, sharedflag, object_files, root1, libext, libs)
+    
+  if verbose:
+    print s
+
+  # Doesn't work on Windows anyway      
+  #else:
+  #  s = s + ' 2> /dev/null' #Suppress warnings
+  
+  try:  
+    err=os.system(s)
+    if err != 0:        
+        raise 'Attempting to link %s failed - please try manually' %root1     
+  except:
+    raise 'Could not link %s - please try manually' %root1
+    
+
+def can_use_C_extension(filename):
+    """Determine whether specified C-extension
+    can and should be used.
+    """
+
+    from anuga.config import use_extensions
+
+    from os.path import splitext
+
+    root, ext = splitext(filename)
+    
+    C=False
+    if use_extensions:
+        try:
+            s = 'import %s' %root
+            #print s
+            exec(s)
+        except:
+            try:
+                open(filename)
+            except:
+                msg = 'C extension %s cannot be opened' %filename
+                print msg                
+            else:    
+                print '------- Trying to compile c-extension %s' %filename
+            
+                try:
+                    compile(filename)
+                except:
+                    print 'WARNING: Could not compile C-extension %s'\
+                          %filename
+                else:
+                    try:
+                        exec('import %s' %root)
+                    except:
+                        msg = 'C extension %s seems to compile OK, '
+                        msg += 'but it can still not be imported.'
+                        raise msg
+                    else:
+                        C=True
+        else:
+            C=True
+            
+    if not C:
+        pass
+        print 'NOTICE: C-extension %s not used' %filename
+
+    return C
+
+
+def set_python_dll_path():
+  """ On windows, find which of the two usual hiding places the python
+  dll is located.
+
+  If the file can't be found, return None.
+  """
+  import sys
+  from os import access, F_OK
+  
+  version = sys.version[:3]
+  v = version.replace('.','')
+  dllfilename = 'python%s.dll' %(v)
+  libs = os.path.join(sys.exec_prefix,dllfilename)
+
+  is_found = True   
+  if access(libs,F_OK) == 0 :
+    # Sometimes python dll is hidden in %WINDIR%/system32
+    libs = os.path.join(os.environ.get('WINDIR', 'C:\WINNT'), 'system32',
+                        dllfilename)
+    if access(libs,F_OK) == 0 :
+      # could not find the dll
+      libs = os.path.join(sys.exec_prefix,dllfilename)
+      is_found = False
+  return libs, is_found
+
+def check_python_dll():
+  libs, is_found = set_python_dll_path()
+  if not is_found:
+    print "%s not found.\nPlease install.\nIt is available on the web." \
+              %(libs)
+    import sys; sys.exit()
+    
+      
+if __name__ == '__main__':
+  from os.path import splitext
+  
+  if sys.platform == 'win32':
+    check_python_dll()
+  if len(sys.argv) > 1:
+      files = sys.argv[1:]
+      for filename in files:
+          root, ext = splitext(filename)
+
+          if ext <> '.c':
+              print 'WARNING (compile.py): Skipping %s. I only compile C-files.' %filename
+      
+  else:  
+      #path = os.path.split(sys.argv[0])[0] or os.getcwd()
+      path = '.'
+      files = os.listdir(path)
+      
+      
+
+  for filename in files:
+      root, ext = splitext(filename)
+
+      if ext == '.c':
+          for x in ['.dll', '.so']:
+              try:
+                  os.remove(root + x)
+              except:
+                  pass
+
+          print separation_line
+          print 'Trying to compile c-extension %s in %s'\
+                %(filename, os.getcwd())
+          try:
+            if filename == 'triang.c': 
+              compile(['triang.c','triangle.c'])
+            elif  filename == 'mesh_engine_c_layer.c': 
+              compile(['mesh_engine_c_layer.c','triangle.c'])
+            else:
+              compile(filename)
+          except Exception, e:
+              msg = 'Could not compile C extension %s\n' %filename
+              msg += str(e)
+              raise Exception, msg
+          else:
+              print 'C extension %s OK' %filename
+          print    
+        
+
Index: /anuga_work/development/anuga_1d/dam_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/dam_padarn.py	(revision 6452)
+++ /anuga_work/development/anuga_1d/dam_padarn.py	(revision 6453)
@@ -10,5 +10,6 @@
 
 def analytical_sol(C,t):
-    
+    if t==0:
+        t=0.0001
     #t  = 0.0     # time (s)
     # gravity (m/s^2)
@@ -90,8 +91,10 @@
 yieldstep = finaltime
 L = 2000.0     # Length of channel (m)
-number_of_cells = [200]
+number_of_cells = [810]
 k = 0
-widths = [1]
+widths = [1,2,5]
 heights= []
+velocities = []
+
 for i in range(len(widths)):
     k=widths[i]
@@ -118,5 +121,5 @@
         t0 = time.time()
 
-        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):            
             domain.write_time()
 
@@ -126,6 +129,6 @@
         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))
+        h_error = 1.0/(N)*sum(abs(h-HeightC))
+        u_error = 1.0/(N)*sum(abs(uh-DischargeC))
         #print "h_error %.10f" %(h_error[k])
         #print "uh_error %.10f"% (uh_error[k])
@@ -134,9 +137,11 @@
         X = domain.vertices
         heights.append(domain.quantities['height'].vertex_values)
-        VelocityQ = domain.quantities['velocity'].vertex_values
+        velocities.append( domain.quantities['velocity'].vertex_values)
         #stage = domain.quantities['stage'].vertex_values
         h, uh, u = analytical_sol(X.flat,domain.time)
         x = X.flat
-    
+            
+        print "Error in height", h_error
+        print "Error in xmom", u_error
         #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
 from pylab import *
@@ -144,21 +149,23 @@
 import matplotlib.axes3d as p3
 print 'test 2'
-hold(False)
+#hold(False)
 print 'test 3'
 plot1 = subplot(211)
 print 'test 4'
-plot(x,h,x,heights[0].flat)
+plot(x,heights[0].flat,x,h)
 print 'test 5'
 plot1.set_ylim([-1,11])
 xlabel('Position')
 ylabel('Stage')
-legend(('Analytical Solution', 'Numerical Solution'),
-           'upper right', shadow=True)
+legend(('Numerical Solution','Analytic Soltuion'), 'upper right', shadow=True)
 plot2 = subplot(212)
-plot(x,u,x,VelocityQ.flat)
+
+
+plot(x,velocities[0].flat,x,u)
 plot2.set_ylim([-35,35])
-    
 xlabel('Position')
 ylabel('Velocity')
+
+print heights[0].flat-heights[1].flat
    
 file = "dry_bed_"
@@ -168,5 +175,3 @@
 show()
     
-    
-#print "Error in height", h_error
-#print "Error in xmom", uh_error
+
Index: /anuga_work/development/anuga_1d/domain.py
===================================================================
--- /anuga_work/development/anuga_1d/domain.py	(revision 6452)
+++ /anuga_work/development/anuga_1d/domain.py	(revision 6453)
@@ -916,8 +916,8 @@
         # Update ghosts
         self.update_ghosts()
-
+        
         # Initial update of vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
+        self.distribute_to_vertices_and_edges()  
+        
         # Update extrema if necessary (for reporting)
         self.update_extrema()
Index: /anuga_work/development/anuga_1d/dry_dam_vel.py
===================================================================
--- /anuga_work/development/anuga_1d/dry_dam_vel.py	(revision 6452)
+++ /anuga_work/development/anuga_1d/dry_dam_vel.py	(revision 6453)
@@ -78,15 +78,9 @@
 for j in range(N+1):
     points[j] = j*cell_len
-
-boundary = {}
-boundary[0,0] = 'left'
-boundary[N-1,1] = 'right'
-
-domain = Domain(points, boundary = boundary)
-
+        
+domain = Domain(points)
+    
 domain.set_quantity('stage', stage)
-
-Br =  Reflective_boundary(domain)
-domain.set_boundary({'left':  Br, 'right': Br})
+domain.set_boundary({'exterior': Reflective_boundary(domain)})
 domain.order = 2
 domain.set_timestepping_method('euler')
Index: /anuga_work/development/anuga_1d/old_compile_all.py
===================================================================
--- /anuga_work/development/anuga_1d/old_compile_all.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/old_compile_all.py	(revision 6453)
@@ -0,0 +1,11 @@
+import os
+
+buildroot = os.getcwd()
+
+#Attempt to compile all extensions
+
+execfile('/home/steve/ga/anuga_core/source/anuga/utilities/compile.py')
+
+os.chdir(buildroot)    
+#execfile('test_all.py')
+    
Index: /anuga_work/development/anuga_1d/parabolic_cannal_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/parabolic_cannal_padarn.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/parabolic_cannal_padarn.py	(revision 6453)
@@ -0,0 +1,122 @@
+import os
+import random
+from math import sqrt, pow, pi ,sin, cos
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(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
+
+def width(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=1
+    return y
+
+def bed(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
+
+import time
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 2500.0   
+# Define the number of cells
+number_of_cells = [200]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N*4 # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+points = zeros(N+1,Float)
+for i in range(N+1):
+        points[i] = -2*L +i*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Define random array for width
+randomarray=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray[j]=random.normalvariate(5,1)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+domain.set_quantity('elevation',bed)
+
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+b = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+b
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,b,x,stage)
+ 
+plot1.set_ylim([-1,20])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
Index: /anuga_work/development/anuga_1d/quantity.py
===================================================================
--- /anuga_work/development/anuga_1d/quantity.py	(revision 6452)
+++ /anuga_work/development/anuga_1d/quantity.py	(revision 6453)
@@ -135,4 +135,5 @@
             #Use function specific method
             self.set_function_values(X, location)
+          
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
             if location == 'centroids':
@@ -163,14 +164,16 @@
         
         if location == 'centroids':
-            
+         
             P = self.domain.centroids
             self.set_values(f(P), location)
         else:
             #Vertices
+            
             P = self.domain.get_vertices()
             
             for i in range(2):
-                self.vertex_values[:,i] = f(P[:,i])        
-
+               
+                self.vertex_values[:,i] = f(P[:,i])
+               
     def set_array_values(self, values, location='vertices'):
         """Set values for quantity
Index: /anuga_work/development/anuga_1d/radial_dam_break.py
===================================================================
--- /anuga_work/development/anuga_1d/radial_dam_break.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/radial_dam_break.py	(revision 6453)
@@ -0,0 +1,165 @@
+import os
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Radial Dam Break Test"
+
+# Define functions for initial quantities
+
+def initial_area(x):
+    y=zeros(len(x),Float)
+    for i in range (len(x)):
+         if x[i]<=40:
+            y[i]=10*(x[i]*pi*2)
+         else:
+            y[i]=1*(x[i]*pi*2)
+    return y
+        
+
+def width(x):
+    return 2*pi*x
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 2.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 100.0   
+# Define the number of cells
+number_of_cells = [100]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+
+
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+## ################# REPEAT ABOVE WITH MORE CELLS
+## # Set final time and yield time for simulation
+## finaltime2 = 2.0
+## yieldstep2 = finaltime2
+
+## # Length of channel (m)
+## L2 = 100.0   
+## # Define the number of cells
+## number_of_cells2 = [100]
+
+## # Define cells for finite volume and their size
+## N2 = int(number_of_cells2[0])     
+## print "Evaluating domain with %d cells" %N2
+## cell_len2 = L2/N2 # Origin = 0.0
+## points2 = zeros(N2+1,Float)
+
+## # Define the centroid points
+## for j in range(N2+1):
+##     points2[j] = j*cell_len2
+
+## # Create domain with centroid points as defined above        
+## domain2 = Domain(points2)
+
+## # Set initial values of quantities - default to zero
+## domain2.set_quantity('area', initial_area)
+## domain2.set_quantity('width',width)
+
+
+## # Set boundry type, order, timestepping method and limiter
+## domain2.set_boundary({'exterior':Reflective_boundary(domain2)})
+## domain2.order = 2
+## domain2.set_timestepping_method('rk2')
+## domain2.set_CFL(1.0)
+## domain2.set_limiter("vanleer")
+## #domain.h0=0.0001
+
+## # Start timer
+## t02 = time.time()
+
+## for t in domain2.evolve(yieldstep = yieldstep2, finaltime = finaltime2):
+##     domain2.write_time()
+
+## N2 = float(N2)
+## HeightC2 = domain2.quantities['height'].centroid_values
+## DischargeC2 = domain2.quantities['discharge'].centroid_values
+## C2 = domain2.centroids
+## print 'That took %.2f seconds' %(time.time()-t02)
+## X2 = domain2.vertices
+## HeightQ2 = domain2.quantities['height'].vertex_values
+## VelocityQ2 = domain2.quantities['velocity'].vertex_values
+## x2 = X2.flat
+## z2 = domain2.quantities['elevation'].vertex_values.flat
+## stage2=HeightQ2.flat+z2
+
+## ################# REPEAT ABOVE
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+import pickle
+f=open('highresdam.txt','r')
+## Data=[[x,stage],[x,VelocityQ.flat]]
+## pickle.dump(Data,f)
+## f.close()
+[[x2,stage2],[x3,Velocity2]]=pickle.load(f)
+
+
+hold(False)
+plot1 = subplot(211)
+
+plot(x2,stage2,x,stage,'.')
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x3,Velocity2,x,VelocityQ.flat,'.')
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+
+
+
Index: /anuga_work/development/anuga_1d/shallow_water_domain_steve.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain_steve.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/shallow_water_domain_steve.py	(revision 6453)
@@ -0,0 +1,1402 @@
+"""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"
+
+        #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(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 basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
+                                       skip_initial_step):
+
+            #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_python(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(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_ext_steve import compute_fluxes_ext_short
+
+    domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed,height,velocity)
+
+
+
+# ###################################
+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.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.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.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]
+
+        #normal = self.normals[vol_id,edge_id]
+        
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, quantities=None):
+        Boundary.__init__(self)
+
+        if quantities is None:
+            msg = 'Must specify one value for each evolved quantity, w,uh,z,h,u'
+            raise msg
+
+        from Numeric import array, Float
+        self.quantities=array(quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.quantities
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.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/stead_state_varying_width.py
===================================================================
--- /anuga_work/development/anuga_1d/stead_state_varying_width.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/stead_state_varying_width.py	(revision 6453)
@@ -0,0 +1,158 @@
+import os
+from math import sqrt, pow, pi 
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Radial Dam Break Test"
+
+# Define functions for initial quantities
+
+def initial_area(x):
+    y=zeros(len(x),Float)
+    for i in range (len(x)):
+         if x[i]<=40:
+            y[i]=10*(x[i]*pi*2)
+         else:
+            y[i]=1*(x[i]*pi*2)
+    return y
+        
+
+def width(x):
+    return 2*pi*x
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 2.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 100.0   
+# Define the number of cells
+number_of_cells = [200]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+
+
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+################# REPEAT ABOVE WITH MORE CELLS
+# Set final time and yield time for simulation
+finaltime2 = 2.0
+yieldstep2 = finaltime2
+
+# Length of channel (m)
+L2 = 100.0   
+# Define the number of cells
+number_of_cells2 = [2000]
+
+# Define cells for finite volume and their size
+N2 = int(number_of_cells2[0])     
+print "Evaluating domain with %d cells" %N2
+cell_len2 = L2/N2 # Origin = 0.0
+points2 = zeros(N2+1,Float)
+
+# Define the centroid points
+for j in range(N2+1):
+    points2[j] = j*cell_len2
+
+# Create domain with centroid points as defined above        
+domain2 = Domain(points2)
+
+# Set initial values of quantities - default to zero
+domain2.set_quantity('area', initial_area)
+domain2.set_quantity('width',width)
+
+
+# Set boundry type, order, timestepping method and limiter
+domain2.set_boundary({'exterior':Reflective_boundary(domain2)})
+domain2.order = 2
+domain2.set_timestepping_method('rk2')
+domain2.set_CFL(1.0)
+domain2.set_limiter("vanleer")
+#domain.h0=0.0001
+
+# Start timer
+t02 = time.time()
+
+for t in domain2.evolve(yieldstep = yieldstep2, finaltime = finaltime2):
+    domain2.write_time()
+
+N2 = float(N2)
+HeightC2 = domain2.quantities['height'].centroid_values
+DischargeC2 = domain2.quantities['discharge'].centroid_values
+C2 = domain2.centroids
+print 'That took %.2f seconds' %(time.time()-t02)
+X2 = domain2.vertices
+HeightQ2 = domain2.quantities['height'].vertex_values
+VelocityQ2 = domain2.quantities['velocity'].vertex_values
+x2 = X2.flat
+z2 = domain2.quantities['elevation'].vertex_values.flat
+stage2=HeightQ2.flat+z2
+
+################# REPEAT ABOVE
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+plot1 = subplot(211)
+
+plot(x,stage,'ro',x2,stage2)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+
+
+
Index: /anuga_work/development/anuga_1d/straight_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/straight_padarn.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/straight_padarn.py	(revision 6453)
@@ -0,0 +1,165 @@
+import os
+from math import sqrt, pi
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+h1 = 5.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
+
+#def newLinePlot(title='Simple Plot'):
+#   import Gnuplot
+#    gg = Gnuplot.Gnuplot(persist=0)
+#    gg.terminal(postscript)
+#    gg.title(title)
+#    gg('set data style linespoints') 
+#    gg.xlabel('x')
+#    gg.ylabel('y')
+#    return gg
+
+#def linePlot(gg,x1,y1,x2,y2):
+#    import Gnuplot
+#    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
+#    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
+#    g.plot(plot1,plot2)
+
+h2=10.0
+
+print "TEST 1D-SOLUTION III -- DRY BED"
+
+def stage(x):
+    return 5
+
+def width(x):
+    return 1
+
+def bot(x):
+    return x/500
+ 
+import time
+
+finaltime = 20.0
+yieldstep = finaltime
+L = 2000.0     # Length of channel (m)
+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)
+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('area', stage)
+    domain.set_quantity('width',width)
+    domain.set_quantity('elevation',bot)
+    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)
+    HeightC = domain.quantities['height'].centroid_values
+    DischargeC = domain.quantities['discharge'].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))
+    #print "h_error %.10f" %(h_error[k])
+    #print "uh_error %.10f"% (uh_error[k])
+    k = k+1
+    print 'That took %.2f seconds' %(time.time()-t0)
+    X = domain.vertices
+    HeightQ = domain.quantities['height'].vertex_values
+    VelocityQ = domain.quantities['velocity'].vertex_values
+    #stage = domain.quantities['stage'].vertex_values
+    #h, uh, u = analytical_sol(X.flat,domain.time)
+    x = X.flat
+    z =bot(x)
+    w=HeightQ.flat+z
+    
+    #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    from pylab import *
+    import pylab as p
+    import matplotlib.axes3d as p3
+    print 'test 2'
+    hold(False)
+    print 'test 3'
+    plot1 = subplot(211)
+    print 'test 4'
+    plot(x,z,x,w)
+    print 'test 5'
+    plot1.set_ylim([-1,11])
+    xlabel('Position')
+    ylabel('Stage')
+    #legend(( 'Numerical Solution'),
+    #       'upper right', shadow=True)
+    plot2 = subplot(212)
+    plot(x,VelocityQ.flat)
+    plot2.set_ylim([-2,2])
+    
+    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/test_padarn.py
===================================================================
--- /anuga_work/development/anuga_1d/test_padarn.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/test_padarn.py	(revision 6453)
@@ -0,0 +1,130 @@
+import os
+import random
+from math import sqrt, pow, pi ,sin, cos
+from channel_domain import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_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 bed(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 initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=initial_stage([x[i]])-bed([x[i]])
+
+    return y
+
+def width(x):
+    return 1
+ 
+import time
+
+
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 2500.0
+# Define the number of cells
+number_of_cells = [50]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = 4*L/N # Origin = 0.0
+points = zeros(N+1,Float)
+for i in range(N+1):
+        points[i] = -2*L +i*cell_len
+
+
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('elevation',bed)
+domain.set_quantity('width',width)
+
+# Set boundry type, order, timestepping method and limiter
+domain.set_boundary({'exterior':Reflective_boundary(domain)})
+domain.order = 2
+domain.set_timestepping_method('euler')
+domain.set_CFL(1.0)
+domain.set_limiter("vanleer")
+#domain.h0=0.0001
+
+# Start timer
+t0 = time.time()
+
+
+
+
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['area'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,z,x,stage)
+ 
+plot1.set_ylim([-1,35])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
Index: /anuga_work/development/anuga_1d/test_shallow_water_domain_suggestion1.py
===================================================================
--- /anuga_work/development/anuga_1d/test_shallow_water_domain_suggestion1.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/test_shallow_water_domain_suggestion1.py	(revision 6453)
@@ -0,0 +1,533 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from shallow_water_domain_suggestion1 import *
+from shallow_water_domain_suggestion1 import flux_function as domain_flux_function
+
+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_domain_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
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+        #print edgeflux
+
+
+
+        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
+
+        #print domainedgeflux
+
+        assert allclose(domainedgeflux, edgeflux)
+        assert allclose(domainmaxspeed, maxspeed)
+
+    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)})
+
+        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 = 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)
Index: /anuga_work/development/anuga_1d/var_depth_only.py
===================================================================
--- /anuga_work/development/anuga_1d/var_depth_only.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/var_depth_only.py	(revision 6453)
@@ -0,0 +1,171 @@
+import os
+import random
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=(10-randomarray[i])
+
+   
+    return y
+
+def bed(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=randomarray[i]
+        randomarray[i]=random.normalvariate(3,1)
+    return y
+
+def stage(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=3+x[i]/2000
+    return y
+  
+
+
+def width(x):
+    return 1
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 1
+yieldstep = 0.1
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [50]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Define random array for width
+
+randomarray=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray[j]=random.normalvariate(3,1)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('stage',stage)
+#domain.set_quantity('area', initial_area)
+domain.set_quantity('elevation',bed)
+domain.set_quantity('width',width)
+domain.setstageflag = True
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+i=1
+
+print domain.quantities['elevation'].vertex_values
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+    N = float(N)
+    HeightC = domain.quantities['height'].centroid_values
+    DischargeC = domain.quantities['discharge'].centroid_values
+    C = domain.centroids
+    print 'That took %.2f seconds' %(time.time()-t0)
+    X = domain.vertices
+    HeightQ = domain.quantities['height'].vertex_values
+    VelocityQ = domain.quantities['velocity'].vertex_values
+    x = X.flat
+    z = domain.quantities['elevation'].vertex_values.flat
+    stage=HeightQ.flat+z
+
+
+
+
+    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
+    import pylab as p
+#hold(False)
+   
+    plot1=subplot(211)
+
+    plot(x,stage,x,z)
+ 
+    plot1.set_ylim([-1,11])
+    xlabel('Position')
+    ylabel('Stage')
+    legend(('Solution', 'Bed Height'),
+                'upper right', shadow=True)
+
+    plot2=subplot(212)
+    plot(x,VelocityQ.flat)
+    plot2.set_ylim([-5,5])
+    xlabel('Position') 
+    ylabel('Velocity')
+
+    savefig(str(i)+'.png')
+    i+=1
+    plot1.clear()
+
+print domain.quantities['elevation'].vertex_values
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+#hold(False)
+   
+plot1=subplot(211)
+
+plot(x,stage,x,z)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Solution', 'Bed Height'),
+           'upper right', shadow=True)
+
+plot2=subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-5,5])
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+
Index: /anuga_work/development/anuga_1d/var_depth_only2.py
===================================================================
--- /anuga_work/development/anuga_1d/var_depth_only2.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/var_depth_only2.py	(revision 6453)
@@ -0,0 +1,116 @@
+import os
+import random
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        #if x[i]>=500:
+        #    y[i]=0
+        #else:
+       # if x[i]<900:
+        y[i]=(5-bed([x[i]]))*bed([x[i]])
+       # else:
+         #   y[i]=10-bed([x[i]])
+    return y
+
+def bed(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]>100:
+            y[i]=2
+        else:
+            y[i]=1
+    return y
+
+def width(x):
+    return 1
+
+
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime =100.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [10]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('elevation',bed)
+domain.set_quantity('width',bed)
+
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+stage=HeightQ.flat+z
+b = domain.quantities['width'].vertex_values.flat
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,z,x,stage,x,b)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+## legend(('Analytical Solution', 'Numerical Solution'),
+##            'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
Index: /anuga_work/development/anuga_1d/var_width_depth.py
===================================================================
--- /anuga_work/development/anuga_1d/var_width_depth.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/var_width_depth.py	(revision 6453)
@@ -0,0 +1,143 @@
+import os
+import random
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=(10-randomarray[i])
+
+   
+    return y
+
+def width(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=randomarray[i]
+        randomarray[i]=random.normalvariate(3,1)
+    return y
+
+def bed(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=randomarray2[i]
+        randomarray2[i]=random.normalvariate(3,1)
+    return y
+
+
+
+def stage(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<200 and x[i]>150:
+            y[i]=8.0
+        else:
+            y[i]=8.0
+    return y
+  
+
+
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [100]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Define random array for width
+randomarray=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray[j]=random.normalvariate(3,1)
+randomarray2=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray2[j]=random.normalvariate(3,1)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('stage',stage,'centroids')
+domain.set_quantity('elevation', bed)
+domain.set_quantity('width',width)
+domain.setstageflag = True
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+print domain.quantities['elevation'].vertex_values
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+print domain.quantities['elevation'].vertex_values
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['elevation'].vertex_values.flat
+b = domain.quantities['width'].vertex_values.flat
+stage=HeightQ.flat+z
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+#hold(False)
+   
+plot1=subplot(211)
+
+plot(x,stage,x,z,x,b)
+ 
+plot1.set_ylim([-5,15])
+xlabel('Position')
+ylabel('Stage')
+legend(('Solution', 'Bed','Width'),
+            'upper right', shadow=True)
+
+plot2=subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+
Index: /anuga_work/development/anuga_1d/var_width_only.py
===================================================================
--- /anuga_work/development/anuga_1d/var_width_only.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/var_width_only.py	(revision 6453)
@@ -0,0 +1,127 @@
+import os
+import random
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=(10-randomarray[i])
+
+   
+    return y
+
+def width(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=randomarray[i]
+        randomarray[i]=random.normalvariate(3,1)
+    return y
+
+def stage(x):
+
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=8.0
+    return y
+  
+
+
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 10.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [50]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Define random array for width
+
+randomarray=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray[j]=random.normalvariate(3,1)
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('stage',stage)
+#domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+domain.setstageflag = True
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+print domain.quantities['elevation'].vertex_values
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+print domain.quantities['elevation'].vertex_values
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['width'].vertex_values.flat
+stage=HeightQ.flat+0
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+#hold(False)
+   
+plot1=subplot(211)
+
+plot(x,stage,x,z)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Solution', 'Width'),
+            'upper right', shadow=True)
+
+plot2=subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-5,5])
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
+
Index: /anuga_work/development/anuga_1d/var_width_only_random.py
===================================================================
--- /anuga_work/development/anuga_1d/var_width_only_random.py	(revision 6453)
+++ /anuga_work/development/anuga_1d/var_width_only_random.py	(revision 6453)
@@ -0,0 +1,114 @@
+import os
+import random
+from math import sqrt, pow, pi 
+from channel_domain_Ab import *
+from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+from config import g, epsilon
+
+
+print "Variable Width Only Test"
+
+# Define functions for initial quantities
+def initial_area(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=(10*randomarray[i])
+    return y
+
+def width(x):
+    y = zeros(len(x),Float)
+    for i in range(len(x)):
+        y[i]=randomarray[i]
+    return y
+
+
+ 
+import time
+
+# Set final time and yield time for simulation
+finaltime = 20.0
+yieldstep = finaltime
+
+# Length of channel (m)
+L = 1000.0   
+# Define the number of cells
+number_of_cells = [200]
+
+# Define cells for finite volume and their size
+N = int(number_of_cells[0])     
+print "Evaluating domain with %d cells" %N
+cell_len = L/N # Origin = 0.0
+points = zeros(N+1,Float)
+
+# Define the centroid points
+for j in range(N+1):
+    points[j] = j*cell_len
+
+# Create domain with centroid points as defined above        
+domain = Domain(points)
+
+# Define random array for width
+randomarray=zeros(len(points),Float)
+for j in range(N+1):
+    randomarray[j]=random.normalvariate(5,1)
+
+
+
+# Set initial values of quantities - default to zero
+domain.set_quantity('area', initial_area)
+domain.set_quantity('width',width)
+
+
+
+# Set boundry type, order, timestepping method and limiter
+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
+
+# Start timer
+t0 = time.time()
+
+
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+
+N = float(N)
+HeightC = domain.quantities['height'].centroid_values
+DischargeC = domain.quantities['discharge'].centroid_values
+C = domain.centroids
+print 'That took %.2f seconds' %(time.time()-t0)
+X = domain.vertices
+HeightQ = domain.quantities['height'].vertex_values
+VelocityQ = domain.quantities['velocity'].vertex_values
+x = X.flat
+z = domain.quantities['width'].vertex_values.flat
+stage=HeightQ.flat+0
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+
+plot(x,z,x,stage)
+ 
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+#legend(('Analytical Solution', 'Numerical Solution'),
+ #          'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
