Index: uga_work/development/flow_1d/channel_flow/channel_domain.py
===================================================================
--- /anuga_work/development/flow_1d/channel_flow/channel_domain.py	(revision 7831)
+++ 	(revision )
@@ -1,588 +1,0 @@
-"""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 = ['stage', 'xmomentum']
-        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
-        other_quantities = ['friction']
-        Generic_Domain.__init__(self,
-                                coordinates = coordinates,
-                                boundary = boundary,
-                                conserved_quantities = conserved_quantities,
-                                evolved_quantities = evolved_quantities,
-                                other_quantities = other_quantities,
-                                tagged_elements = tagged_elements)
-        
-        from config import minimum_allowed_height, g, h0
-        self.minimum_allowed_height = minimum_allowed_height
-        self.g = g
-        self.h0 = h0
-
-        #forcing terms not included in 1d domain ?WHy?
-        self.forcing_terms.append(gravity)
-        #self.forcing_terms.append(manning_friction)
-        #print "\nI have Removed forcing terms line 64 1dsw"
-
-        
-        #Stored output
-        self.store = True
-        self.format = 'sww'
-        self.smooth = True
-
-        
-        #Reduction operation for get_vertex_values
-        from util import mean
-        self.reduction = mean
-        #self.reduction = min  #Looks better near steep slopes
-
-        self.set_quantities_to_be_stored(['stage','xmomentum'])
-
-	self.__doc__ = 'shallow_water_domain'
-
-        self.check_integrity()
-
-
-    def check_integrity(self):
-
-        #Check that we are solving the shallow water wave equation
-
-        msg = 'First conserved quantity must be "stage"'
-        assert self.conserved_quantities[0] == 'stage', msg
-        msg = 'Second conserved quantity must be "xmomentum"'
-        assert self.conserved_quantities[1] == 'xmomentum', msg
-
-        msg = 'First evolved quantity must be "stage"'
-        assert self.evolved_quantities[0] == 'stage', msg
-        msg = 'Second evolved quantity must be "xmomentum"'
-        assert self.evolved_quantities[1] == 'xmomentum', msg
-        msg = 'Third evolved quantity must be "elevation"'
-        assert self.evolved_quantities[2] == 'elevation', msg
-        msg = 'Fourth evolved quantity must be "height"'
-        assert self.evolved_quantities[3] == 'height', msg
-        msg = 'Fifth evolved quantity must be "velocity"'
-        assert self.evolved_quantities[4] == 'velocity', msg
-
-        Generic_Domain.check_integrity(self)
-
-    def compute_fluxes(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        compute_fluxes_vel(self)
-
-    def distribute_to_vertices_and_edges(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        distribute_to_vertices_and_edges_limit_w_u(self)
-        
-
-#=============== End of Channel Domain ===============================
-#-----------------------------------
-# Compute fluxes interface
-#-----------------------------------
-def compute_fluxes(domain):
-    """
-    Python version of compute fluxes (local_compute_fluxes)
-    is available in test_shallow_water_vel_domain.py
-    """
- 
-
-    from Numeric import zeros, Float
-    import sys
-
-	
-    timestep = float(sys.maxint)
-
-    stage = domain.quantities['stage']
-    xmom = domain.quantities['xmomentum']
-    bed = domain.quantities['elevation']
-
-
-    from comp_flux_vel_ext import compute_fluxes_ext
-
-    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)
-
-
-#-----------------------------------
-# Compute flux definition with vel
-#-----------------------------------
-def compute_fluxes_vel(domain):
-    from Numeric import zeros, Float
-    import sys
-
-	
-    timestep = float(sys.maxint)
-
-    stage    = domain.quantities['stage']
-    xmom     = domain.quantities['xmomentum']
-    bed      = domain.quantities['elevation']
-    height   = domain.quantities['height']
-    velocity = domain.quantities['velocity']
-
-
-    from comp_flux_vel_ext import compute_fluxes_vel_ext
-
-    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)
-
-
-#--------------------------------------------------------------------------
-def distribute_to_vertices_and_edges_limit_w_u(domain):
-    """Distribution from centroids to vertices specific to the
-    shallow water wave
-    equation.
-
-    It will ensure that h (w-z) is always non-negative even in the
-    presence of steep bed-slopes by taking a weighted average between shallow
-    and deep cases.
-
-    In addition, all conserved quantities get distributed as per either a
-    constant (order==1) or a piecewise linear function (order==2).
-
-    FIXME: more explanation about removal of artificial variability etc
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
-
-    """
-
-    #from config import optimised_gradient_limiter
-
-    #Remove very thin layers of water
-    #protect_against_infinitesimal_and_negative_heights(domain)  
-
-    import sys
-    from Numeric import zeros, Float
-    from config import epsilon, h0
-	
-    N = domain.number_of_elements
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Bed = domain.quantities['elevation']
-    Height = domain.quantities['height']
-    Velocity = domain.quantities['velocity']
-
-    #Arrays   
-    w_C   = Stage.centroid_values    
-    uh_C  = Xmom.centroid_values    
-    z_C   = Bed.centroid_values
-    h_C   = Height.centroid_values
-    u_C   = Velocity.centroid_values
-	
-    #print id(h_C)
-    for i in range(N):
-	h_C[i] = w_C[i] - z_C[i]						
-        if h_C[i] <= 1.0e-12:
-	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
-	    h_C[i] = 0.0
-	    w_C[i] = z_C[i]
-	    #uh_C[i] = 0.0
-            
- #           u_C[i] = 0.0
- #       else:
- #           u_C[i] = uh_C[i]/h_C[i]
-		
-    h0 = 1.0e-12    
-    for i in range(len(h_C)):
-	if h_C[i] < 1.0e-12:
-	    u_C[i] = 0.0  #Could have been negative
-	    h_C[i] = 0.0
-	else:
-            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
-            #u_C[i]  = uh_C[i]/h_C[i]
-	
-    for name in [ 'velocity', 'stage' ]:
-        Q = domain.quantities[name]
-        if domain.order == 1:
-            Q.extrapolate_first_order()
-        elif domain.order == 2:
-            Q.extrapolate_second_order()
-        else:
-            raise 'Unknown order'
-
-    w_V  = domain.quantities['stage'].vertex_values			
-    z_V  = domain.quantities['elevation'].vertex_values	
-    h_V  = domain.quantities['height'].vertex_values
-    u_V  = domain.quantities['velocity'].vertex_values		
-    uh_V = domain.quantities['xmomentum'].vertex_values	
-
-    h_V[:]    = w_V - z_V
-    for i in range(len(h_C)):
-        for j in range(2):
-            if h_V[i,j] < 0.0 :
-                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
-                dh = h_V[i,j]
-                h_V[i,j] = 0.0
-                w_V[i,j] = z_V[i,j]
-                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
-                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
-                
-    uh_V[:] = u_V * h_V
-
-    
-    return
-
-#---------------------------------------------------------------------------
-def distribute_to_vertices_and_edges_limit_w_uh(domain):
-    """Distribution from centroids to vertices specific to the
-    shallow water wave equation.
-
-    In addition, all conserved quantities get distributed as per either a
-    constant (order==1) or a piecewise linear function (order==2).
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
-
-    """
-
-    import sys
-    from Numeric import zeros, Float
-    from config import epsilon, h0
-	
-    N = domain.number_of_elements
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Bed = domain.quantities['elevation']
-    Height = domain.quantities['height']
-    Velocity = domain.quantities['velocity']
-
-    #Arrays   
-    w_C   = Stage.centroid_values    
-    uh_C  = Xmom.centroid_values    
-    z_C   = Bed.centroid_values
-    h_C   = Height.centroid_values
-    u_C   = Velocity.centroid_values
-	
-
-    for i in range(N):
-	h_C[i] = w_C[i] - z_C[i]						
-        if h_C[i] <= 1.0e-6:
-	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
-	    h_C[i] = 0.0
-	    w_C[i] = z_C[i]
-	    uh_C[i] = 0.0
-            
-    for name in [ 'stage', 'xmomentum']:
-        Q = domain.quantities[name]
-        if domain.order == 1:
-            Q.extrapolate_first_order()
-        elif domain.order == 2:
-            Q.extrapolate_second_order()
-        else:
-            raise 'Unknown order'
-
-    w_V  = domain.quantities['stage'].vertex_values			
-    z_V  = domain.quantities['elevation'].vertex_values	
-    h_V  = domain.quantities['height'].vertex_values
-    u_V  = domain.quantities['velocity'].vertex_values		
-    uh_V = domain.quantities['xmomentum'].vertex_values	
-
-    h_V[:]    = w_V - z_V
-
-    for i in range(len(h_C)):
-        for j in range(2):
-            if h_V[i,j] < 0.0 :
-                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
-                dh = h_V[i,j]
-                h_V[i,j] = 0.0
-                w_V[i,j] = z_V[i,j]
-                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
-                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
-                u_V[i,j] = 0.0
-            if h_V[i,j] < h0:
-                u_V[i,j]
-            else:
-                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
-
-
-#--------------------------------------------------------
-#Boundaries - specific to the shallow_water_vel_domain
-#--------------------------------------------------------
-class Reflective_boundary(Boundary):
-    """Reflective boundary returns same conserved quantities as
-    those present in its neighbour volume but reflected.
-
-    This class is specific to the shallow water equation as it
-    works with the momentum quantities assumed to be the second
-    and third conserved quantities.
-    """
-
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for reflective boundary'
-            raise msg
-
-        #Handy shorthands
-        self.normals  = domain.normals
-        self.stage    = domain.quantities['stage'].vertex_values
-        self.xmom     = domain.quantities['xmomentum'].vertex_values
-        self.bed      = domain.quantities['elevation'].vertex_values
-        self.height   = domain.quantities['height'].vertex_values
-        self.velocity = domain.quantities['velocity'].vertex_values
-
-        from Numeric import zeros, Float
-        #self.conserved_quantities = zeros(3, Float)
-        self.evolved_quantities = zeros(5, Float)
-
-    def __repr__(self):
-        return 'Reflective_boundary'
-
-
-    def evaluate(self, vol_id, edge_id):
-        """Reflective boundaries reverses the outward momentum
-        of the volume they serve.
-        """
-
-        q = self.evolved_quantities
-        q[0] = self.stage[vol_id, edge_id]
-        q[1] = -self.xmom[vol_id, edge_id]
-        q[2] = self.bed[vol_id, edge_id]
-        q[3] = self.height[vol_id, edge_id]
-        q[4] = -self.velocity[vol_id, edge_id]
-
-        #print "In Reflective q ",q
-
-
-        return q
-
-class Dirichlet_boundary(Boundary):
-    """Dirichlet boundary returns constant values for the
-    conserved quantities
-    """
-
-
-    def __init__(self, evolved_quantities=None):
-        Boundary.__init__(self)
-
-        if evolved_quantities is None:
-            msg = 'Must specify one value for each evolved quantity'
-            raise msg
-
-        from Numeric import array, Float
-        self.evolved_quantities=array(evolved_quantities).astype(Float)
-
-    def __repr__(self):
-        return 'Dirichlet boundary (%s)' %self.evolved_quantities
-
-    def evaluate(self, vol_id=None, edge_id=None):
-        return self.evolved_quantities
-
-
-#----------------------------
-#Standard forcing terms:
-#---------------------------
-def gravity(domain):
-    """Apply gravitational pull in the presence of bed slope
-    """
-
-    from util import gradient
-    from Numeric import zeros, Float, array, sum
-
-
-
-    Stage     = domain.quantities['stage']
-    Xmom      = domain.quantities['xmomentum']
-    Elevation = domain.quantities['elevation']
-    Height    = domain.quantities['height']
-
-    xmom_ud  = Xmom.explicit_update
-    #stage_ud = Stage.explicit_update
-
-
-    #h = Stage.vertex_values - Elevation.vertex_values
-    h = Height.vertex_values
-    b = Elevation.vertex_values
-    w = Stage.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-
-    for k in range(domain.number_of_elements):
-        avg_h = 0.5*(h[k,0] + h[k,1])
-
-        #Compute bed slope
-        x0, x1 = x[k,:]
-        b0, b1 = b[k,:]
-        bx = gradient(x0, x1, b0, b1)
-        
-        #Update momentum (explicit update is reset to source values)
-        xmom_ud[k] += -g*bx*avg_h
-        #stage_ud[k] = 0.0
- 
- 
-def manning_friction(domain):
-    """Apply (Manning) friction to water momentum
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-    #vh = domain.quantities['ymomentum'].centroid_values
-    eta = domain.quantities['friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g
-
-    for k in range(N):
-        if eta[k] >= eps:
-            if h[k] >= eps:
-            	#S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
-                S = -g * eta[k]**2 * uh[k]
-            	S /= h[k]**(7.0/3)
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
-            	#ymom_update[k] += S*vh[k]
-
-def linear_friction(domain):
-    """Apply linear friction to water momentum
-
-    Assumes quantity: 'linear_friction' to be present
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-    tau = domain.quantities['linear_friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-
-    for k in range(N):
-        if tau[k] >= eps:
-            if h[k] >= eps:
-            	S = -tau[k]/h[k]
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
-
-
-
-def check_forcefield(f):
-    """Check that f is either
-    1: a callable object f(t,x,y), where x and y are vectors
-       and that it returns an array or a list of same length
-       as x and y
-    2: a scalar
-    """
-
-    from Numeric import ones, Float, array
-
-
-    if callable(f):
-        #N = 3
-        N = 2
-        #x = ones(3, Float)
-        #y = ones(3, Float)
-        x = ones(2, Float)
-        #y = ones(2, Float)
-        
-        try:
-            #q = f(1.0, x=x, y=y)
-            q = f(1.0, x=x)
-        except Exception, e:
-            msg = 'Function %s could not be executed:\n%s' %(f, e)
-	    #FIXME: Reconsider this semantics
-            raise msg
-
-        try:
-            q = array(q).astype(Float)
-        except:
-            msg = 'Return value from vector function %s could ' %f
-            msg += 'not be converted into a Numeric array of floats.\n'
-            msg += 'Specified function should return either list or array.'
-            raise msg
-
-        #Is this really what we want?
-        msg = 'Return vector from function %s ' %f
-        msg += 'must have same lenght as input vectors'
-        assert len(q) == N, msg
-
-    else:
-        try:
-            f = float(f)
-        except:
-            msg = 'Force field %s must be either a scalar' %f
-            msg += ' or a vector function'
-            raise msg
-    return f
-
Index: uga_work/development/flow_1d/channel_flow/channel_domain2.py
===================================================================
--- /anuga_work/development/flow_1d/channel_flow/channel_domain2.py	(revision 7831)
+++ 	(revision )
@@ -1,852 +1,0 @@
-"""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: uga_work/development/flow_1d/channel_flow/channel_domain_ext_backup.c
===================================================================
--- /anuga_work/development/flow_1d/channel_flow/channel_domain_ext_backup.c	(revision 7831)
+++ 	(revision )
@@ -1,510 +1,0 @@
-#include "Python.h"
-#include "Numeric/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-const double pi = 3.14159265358979;
-
-
-// Shared code snippets
-#include "util_ext.h"
-
-
-/* double max(double a, double b) { */
-/* 	double z; */
-/* 	z=(a>b)?a:b; */
-/* 	return z;} */
-
-/* double min(double a, double b) { */
-/* 	double z; */
-/* 	z=(a<b)?a:b; */
-/* 	return z;} */
-
-
-// Function to obtain speed from momentum and depth.
-// This is used by flux functions
-// Input parameters uh and h may be modified by this function.
-double _compute_speed(double *uh,
-		      double *h,
-		      double epsilon,
-		      double h0) {
-
-  double u;
-
-  if (*h < epsilon) {
-    *h = 0.0;  //Could have been negative
-     u = 0.0;
-  } else {
-    u = *uh/(*h + h0/ *h);
-  }
-
-
-  // Adjust momentum to be consistent with speed
-  *uh = u * *h;
-
-  return u;
-}
-
-
-
-//Innermost flux function (using w=z+h)
-int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
-			   double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed) {
-
-    int i;
-    double flux_left[2], flux_right[2];
-    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
-     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
- double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
- double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
-    double s_maxl, s_minl,s_maxr,s_minr, denom;
-    double zphalf,zmhalf,hleftstar,hrightstar;
-    double fluxtemp1,fluxtemp0,speedtemp;
-
-    zmhalf = max(q_leftm[2],q_leftp[2]);
-    zphalf = max(q_rightm[2],q_rightp[2]);   
-
-    a_leftm  = q_leftm[0];
-    d_leftm = q_leftm[1];
-    z_leftm  = q_leftm[2];
-    h_leftm  = max(0,q_leftm[3]+q_leftm[2]-zmhalf);
-    u_leftm  = q_leftm[4];
-    b_leftm  = q_leftm[5];
-    w_leftm  = h_leftm+z_leftm;
-
-    a_leftp  = q_leftp[0];
-    d_leftp = q_leftp[1];
-    z_leftp  = q_leftp[2];
-    h_leftp  = max(0,q_leftp[3]+q_leftp[2]-zmhalf);
-    u_leftp  = q_leftp[4];
-    b_leftp  = q_leftp[5];
-    w_leftp  = h_leftp+z_leftp;
-
-    a_rightm  = q_rightm[0];
-    d_rightm = q_rightm[1];
-    z_rightm  = q_rightm[2];
-    h_rightm  = max(0,q_rightm[3]+q_rightm[2]-zphalf);
-    u_rightm  = q_rightm[4];
-    b_rightm  = q_rightm[5];
-    w_rightm  = h_rightm+z_rightm;
-
-    a_rightp  = q_rightp[0];
-    d_rightp = q_rightp[1];
-    z_rightp  = q_rightp[2];
-    h_rightp  = max(0,q_rightp[3]+q_rightp[2]-zphalf);
-    u_rightp  = q_rightp[4];
-    b_rightp  = q_rightp[5];
-    w_rightp  = h_rightp+z_rightp;
-
-    hleftstar = q_leftp[3];
-    hrightstar = q_rightm[3];   
-
-
-
-
-    soundspeed_leftp = sqrt(g*h_leftp);
-    soundspeed_leftm = sqrt(g*h_leftm);
-    soundspeed_rightp = sqrt(g*h_rightp);
-    soundspeed_rightm = sqrt(g*h_rightm);
-
-    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
-    if (s_maxl < 0.0) s_maxl = 0.0;
-
-    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
-    if (s_minl > 0.0) s_minl = 0.0;
-
-    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
-    if (s_maxr < 0.0) s_maxr = 0.0;
-
-    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
-    if (s_minr > 0.0) s_minr = 0.0;
-
-    // Flux formulas for left hand side
-    flux_left[0] = d_leftm;
-    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
-
-    flux_right[0] = d_leftp;
-    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
-
-    
-    // Flux computation for left hand side
-    denom = s_maxl-s_minl;
-    if (denom < epsilon) {
-	for (i=0; i<2; i++) edgeflux[i] = 0.0;
-    } else {
-	edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
-//	edgeflux[0] += s_maxl*s_minl*(a_leftp-a_leftm);
-	edgeflux[0] /= denom;
-	edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
-	edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
-	edgeflux[1] /= denom;
-    
-
-    }
-    
-	fluxtemp0 = edgeflux[0];
-        fluxtemp1 = edgeflux[1];
-    
-
-    // Flux formulas for right hand side
-    flux_left[0] = d_rightm;
-    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm;
-
-    flux_right[0] = d_rightp;
-    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp;
-    
-   
-   
-    // Flux computation for right hand side
-    denom = s_maxr-s_minr;
-    if (denom < epsilon) {
-	for (i=0; i<2; i++) edgeflux[i] = 0.0;
-    } else {
-	edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
-//	    	edgeflux[0] += s_maxr*s_minr*(a_rightp-a_rightm);
-	edgeflux[0] /= denom;
-	edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
-	edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
-	edgeflux[1] /= denom;
-    
-
-    }
-    
-   
-    edgeflux[0]=edgeflux[0]-fluxtemp0;
-    edgeflux[1]=edgeflux[1]-fluxtemp1;
-    
-    
-     
-    edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
-
-    
-
-      
-
-     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
-         	// Maximal wavespeed
-	if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
-	    *max_speed = 0.0;
-	}else{
-	speedtemp = max(fabs(s_maxl),fabs(s_minl));
-        speedtemp = max(speedtemp,fabs(s_maxr));
-        speedtemp = max(speedtemp,fabs(s_minr));
-	*max_speed = speedtemp;
-	}
-
-//printf("%f\n",h_right);
-    return 0;
-}
-
-
-
-
-
-// Computational function for flux computation
-double _compute_fluxes_channel_ext(double cfl,
-			       double timestep,
-			       double epsilon,
-			       double g,
-			       double h0,
-			       long* neighbours,
-			       long* neighbour_vertices,
-			       double* normals,
-			       double* areas,
-			       double* area_edge_values,
-			       double* discharge_edge_values,
-			       double* bed_edge_values,
-			       double* height_edge_values,
-			       double* velocity_edge_values,
-			       double* width_edge_values,
-			       double* area_boundary_values,
-			       double* discharge_boundary_values,
-			       double* bed_boundary_values,
-			       double* height_boundary_values,
-			       double* velocity_boundary_values,
-			       double* width_boundary_values,
-			       double* area_explicit_update,
-			       double* discharge_explicit_update,
-			       int number_of_elements,
-			       double* max_speed_array) {
-
-    double flux[2], qlm[6],qlp[6], qrm[6],qrp[6], edgeflux[2];
-    double max_speed, normal;
-    int k, i, ki, n, m, nm=0;
-    double zstar;
-    for (k=0; k<number_of_elements; k++) {
-	flux[0] = 0.0;
-	flux[1] = 0.0;
-
-       
-	    ki = k*2;
-           
-
-         n = neighbours[ki];
-	    if (n<0) {
-		m = -n-1;
-
-	    qlm[0] = area_boundary_values[m];
-	    qlm[1] = discharge_boundary_values[m];
-	    qlm[2] = bed_boundary_values[m];
-	    qlm[3] = height_boundary_values[m];
-	    qlm[4] = velocity_boundary_values[m];
-            qlm[5] = width_boundary_values[m];
-
-	    }else{
-	m = neighbour_vertices[ki];
-		nm = n*2+m;
-	    
-
-	    qlm[0] = area_edge_values[nm];
-	    qlm[1] = discharge_edge_values[nm];
-	    qlm[2] = bed_edge_values[nm];
-	    qlm[3] = height_edge_values[nm];
-	    qlm[4] = velocity_edge_values[nm];
-            qlm[5] = width_edge_values[nm];
-    }
-	    qlp[0] = area_edge_values[ki];
-	    qlp[1] = discharge_edge_values[ki];
-	    qlp[2] = bed_edge_values[ki];
-	    qlp[3] = height_edge_values[ki];
-	    qlp[4] = velocity_edge_values[ki];
-            qlp[5] = width_edge_values[ki];
-
-         ki = k*2+1;
-           
-
-         n = neighbours[ki];
-	    if (n<0) {
-		m = -n-1;
-            qrp[0] = area_boundary_values[m];
-	    qrp[1] = discharge_boundary_values[m];
-	    qrp[2] = bed_boundary_values[m];
-	    qrp[3] = height_boundary_values[m];
-	    qrp[4] = velocity_boundary_values[m];
-            qrp[5] = width_boundary_values[m];
-          
-
-
-	    }else{
-	m = neighbour_vertices[ki];
-		nm = n*2+m;
-	    
-           
-            qrp[0] = area_edge_values[nm];
-	    qrp[1] = discharge_edge_values[nm];
-	    qrp[2] = bed_edge_values[nm];
-	    qrp[3] = height_edge_values[nm];
-	    qrp[4] = velocity_edge_values[nm];
-            qrp[5] = width_edge_values[nm];
-            }
-            qrm[0] = area_edge_values[ki];
-	    qrm[1] = discharge_edge_values[ki];
-	    qrm[2] = bed_edge_values[ki];
-	    qrm[3] = height_edge_values[ki];
-	    qrm[4] = velocity_edge_values[ki];
-            qrm[5] = width_edge_values[ki];
-
-             _flux_function_channel(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
-	    flux[0] -= edgeflux[0];
-	    flux[1] -= edgeflux[1];
-           
-	    // Update timestep based on edge i and possibly neighbour n
-	    if (max_speed > epsilon) {
-		// Original CFL calculation
-
-		timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
-		if (n>=0) {
-		    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
-		}
-	    }
-	 // End edge i (and neighbour n)
-	flux[0] /= areas[k];
-	area_explicit_update[k] = flux[0];
-	flux[1] /= areas[k];
-	discharge_explicit_update[k] = flux[1];
-	//Keep track of maximal speeds
-	max_speed_array[k]=max_speed;
-    }
-    return timestep;
-
-}
-
-
-//-------------------------------------------------------------
-// Old code
-//------------------------------------------------------------
-//Innermost flux function (using w=z+h)
-
-
-
-
-
-// Computational function for flux computation
-
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-
-
-//------------------------------------------------
-// New velocity based compute fluxes
-//------------------------------------------------
-
-PyObject *compute_fluxes_channel_ext(PyObject *self, PyObject *args) {
-
-    PyObject
-	*domain,
-	*area,
-	*discharge,
-	*bed,
-	*height,
-	*velocity,
-        *width;
-
-    PyArrayObject
-	*neighbours,
-	*neighbour_vertices,
-	*normals,
-	*areas,
-	*area_vertex_values,
-	*discharge_vertex_values,
-	*bed_vertex_values,
-	*height_vertex_values,
-	*velocity_vertex_values,
-	*width_vertex_values,
-	*area_boundary_values,
-	*discharge_boundary_values,
-	*bed_boundary_values,
-	*height_boundary_values,
-	*velocity_boundary_values,
-	*width_boundary_values,
-	*area_explicit_update,
-	*discharge_explicit_update,
-	*max_speed_array;
-
-  double timestep, epsilon, g, h0, cfl;
-  int number_of_elements;
-
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dOOOOOOO",
-			&timestep,
-			&domain,
-			&area,
-			&discharge,
-			&bed,
-			&height,
-			&velocity,
-                        &width)) {
-      PyErr_SetString(PyExc_RuntimeError, "comp_flux_channel_ext.c: compute_fluxes_channel_ext could not parse input");
-      return NULL;
-  }
-
-
-    epsilon           = get_python_double(domain,"epsilon");
-    g                 = get_python_double(domain,"g");
-    h0                = get_python_double(domain,"h0");
-    cfl               = get_python_double(domain,"CFL");
-
-
-    neighbours        = get_consecutive_array(domain, "neighbours");
-    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
-    normals           = get_consecutive_array(domain, "normals");
-    areas             = get_consecutive_array(domain, "areas");
-    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
-
-    area_vertex_values      = get_consecutive_array(area,    "vertex_values");
-    discharge_vertex_values       = get_consecutive_array(discharge,     "vertex_values");
-    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
-    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
-    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
-    width_vertex_values      = get_consecutive_array(width, "vertex_values");
-
-    area_boundary_values     = get_consecutive_array(area,     "boundary_values");
-    discharge_boundary_values      = get_consecutive_array(discharge,      "boundary_values");
-    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
-    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
-    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
-    width_boundary_values       = get_consecutive_array(width,     "boundary_values");
-
-
-    area_explicit_update = get_consecutive_array(area, "explicit_update");
-    discharge_explicit_update  = get_consecutive_array(discharge,  "explicit_update");
-
-    number_of_elements = area_vertex_values -> dimensions[0];
-
-    // Call underlying flux computation routine and update
-    // the explicit update arrays
-    timestep = _compute_fluxes_channel_ext(cfl,
-				       timestep,
-				       epsilon,
-				       g,
-				       h0,
-				       (long*) neighbours -> data,
-				       (long*) neighbour_vertices -> data,
-				       (double*) normals -> data,
-				       (double*) areas -> data,
-				       (double*) area_vertex_values -> data,
-				       (double*) discharge_vertex_values -> data,
-				       (double*) bed_vertex_values -> data,
-				       (double*) height_vertex_values -> data,
-				       (double*) velocity_vertex_values -> data,
-				       (double*) width_vertex_values -> data,
-				       (double*) area_boundary_values -> data,
-				       (double*) discharge_boundary_values -> data,
-				       (double*) bed_boundary_values -> data,
-				       (double*) height_boundary_values -> data,
-				       (double*) velocity_boundary_values -> data,
-				       (double*) width_boundary_values -> data,
-				       (double*) area_explicit_update -> data,
-				       (double*) discharge_explicit_update -> data,
-				       number_of_elements,
-				       (double*) max_speed_array -> data);
-
-
-    Py_DECREF(neighbours);
-    Py_DECREF(neighbour_vertices);
-    Py_DECREF(normals);
-    Py_DECREF(areas);
-    Py_DECREF(area_vertex_values);
-    Py_DECREF(discharge_vertex_values);
-    Py_DECREF(bed_vertex_values);
-    Py_DECREF(height_vertex_values);
-    Py_DECREF(velocity_vertex_values);
-    Py_DECREF(width_vertex_values);
-    Py_DECREF(area_boundary_values);
-    Py_DECREF(discharge_boundary_values);
-    Py_DECREF(bed_boundary_values);
-    Py_DECREF(height_boundary_values);
-    Py_DECREF(velocity_boundary_values);
-    Py_DECREF(width_boundary_values);
-    Py_DECREF(area_explicit_update);
-    Py_DECREF(discharge_explicit_update);
-    Py_DECREF(max_speed_array);
-
-
-    // Return updated flux timestep
-    return Py_BuildValue("d", timestep);
-}
-
-
-//-------------------------------
-// Method table for python module
-//-------------------------------
-
-static struct PyMethodDef MethodTable[] = {
-  {"compute_fluxes_channel_ext", compute_fluxes_channel_ext, METH_VARARGS, "Print out"},
-  {NULL}
-};
-
-/* // Module initialisation */
-/* void initcomp_flux_vel_ext(void){ */
-/*   Py_InitModule("comp_flux_vel_ext", MethodTable); */
-/*   import_array(); */
-/* } */
-
-void initchannel_domain_ext(void){
-  Py_InitModule("channel_domain_ext", MethodTable);
-  import_array();
-}
Index: uga_work/development/flow_1d/channel_flow/compile_all.py
===================================================================
--- /anuga_work/development/flow_1d/channel_flow/compile_all.py	(revision 7831)
+++ 	(revision )
@@ -1,15 +1,0 @@
-import os
-
-buildroot = os.getcwd()
-
-#Attempt to compile all extensions
-
-anuga_dir = '..' + os.sep + '..' + os.sep + '..'
-utilities_dir = anuga_dir + os.sep + 'trunk' + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities'
-
-execfile( utilities_dir + os.sep + 'compile.py')
-
-
-os.chdir(buildroot)    
-
-    
