Index: trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain.py	(revision 8153)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain.py	(revision 8153)
@@ -0,0 +1,474 @@
+"""Class Domain -
+1D interval domains for finite-volume computations of the shallow water wave equation.
+
+This module contains a specialisation of class Generic_domain from module generic_domain.py consisting of methods specific to mixed rectangular Pipe/Channel flow using the Shallow Water Wave Equation.
+
+This particular modification of the Domain class implements the ability to
+vary the width of the 1D pipe that the water flows in. As a result the
+conserved variables are different than previous implementations and so are the
+equations. In particular, the conserved quantites generalise the open channel flow conserved quantities.
+
+U_t + E_x = S
+
+where
+
+U = [A, Q] = [b*h, u*b*h]
+E = [Q, Q^2/A + g*b*h^2/2]
+S represents source terms forcing the system
+(e.g. gravity, boundary_stree, friction, wind stress, ...)
+gravity = -g*b*h*z_x
+boundary_stress = 1/2*g*b_x*h^2
+
+and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
+
+The quantities are
+
+symbol    variable name    explanation
+A         area             equivalent Wetted area = b*h
+Q         discharge        flux of water = u*b*h
+x         x                horizontal distance from origin [m]
+z         elevation        elevation of bottom of pipe [m]
+h         height           water height above z [m]
+w         stage            equivalent 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]
+b         width            width of pipe
+eta                        mannings friction coefficient [to appear]
+
+The conserved quantities are A, Q
+--------------------------------------------------------------------------
+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
+
+Add paper details here
+
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
+Padarn Wilson, Geoscience Australia, 2008
+"""
+
+
+from anuga_1d.base.generic_domain import *
+import numpy
+
+
+#Shallow water pipe 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 anuga_1d.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
+
+
+        # forcing terms gravity and boundary stress are included in the flux calculation
+        #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 anuga_1d.base.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'
+
+        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 = 'Sixth evolved quantity must be "width"'
+        assert self.evolved_quantities[5] == 'width', msg
+        msg = 'Seventh evolved quantity must be "stage"'
+        assert self.evolved_quantities[6] == 'stage', msg
+
+        Generic_domain.check_integrity(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_pipe(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 Pipe Domain ===============================
+
+#-----------------------------------
+# Compute flux definition with pipe
+#-----------------------------------
+def compute_fluxes_pipe(domain):
+    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 anuga_1d.pipe.pipe_domain_ext import compute_fluxes_pipe_ext
+    domain.flux_timestep = compute_fluxes_pipe_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 anuga_1d.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
+
+    if domain.setstageflag:
+        a_C[:,] = (w_C-z_C)*b_C
+        domain.setstageflag = False
+
+    if domain.discontinousb:
+        width.extrapolate_second_order()
+        
+
+    h0 = 1.0e-12
+
+    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
+    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
+
+    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
+    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
+    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
+
+    w_C[:] = h_C + z_C
+    
+		
+
+    bed.extrapolate_second_order()
+    
+    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  = area.vertex_values
+    w_V  = stage.vertex_values
+    z_V  = bed.vertex_values
+    h_V  = height.vertex_values
+    u_V  = velocity.vertex_values
+    d_V  = discharge.vertex_values
+    b_V  = width.vertex_values
+
+
+    h_V[:] = w_V-z_V
+
+    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
+    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
+    a_V[:] = b_V*h_V
+    d_V[:] = u_V*a_V
+
+    
+    return
+
+
+#--------------------------------------------------------
+#Boundaries - specific to the pipe_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
+
+        self.evolved_quantities = numpy.zeros(7, numpy.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]
+
+        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
+
+        assert len(evolved_quantities) == 7
+
+        self.evolved_quantities=numpy.array(evolved_quantities,numpy.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
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].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 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: trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c	(revision 8153)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c	(revision 8153)
@@ -0,0 +1,428 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+const double pi = 3.14159265358979;
+
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+// Hard coded flux computation choice
+#define _flux_function_pipe _flux_function_pipe_naive
+
+
+// Formulae for fluxes
+double flux_formula0(double d, double u, double h, double b, double g) {
+  return 1.0; //d;
+}
+
+double flux_formula1(double d, double u, double h, double b, double g) {
+  return 1.0; // u*d + 0.5*g*h*h*b;
+}
+
+// Compute fluxes
+void compute_flux(double flux_left[], double flux_right[], double *edgeflux) {
+  edgeflux[0] = 0.5*(flux_left[0]+flux_right[0]);
+  edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]);
+}
+
+// Compute forcing terms
+double forcingterm(double h_rightm, double h_leftp, double b_rightm, double b_leftp, double z_rightm, double z_leftp, double g) {
+  return 1.0; //-0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightm-z_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightm-b_leftp)*g;
+}
+
+// NAIVE VERSION
+int _flux_function_pipe_naive(double *q_leftm,double *q_leftp, double *q_rightm,
+			      double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed){
+
+  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;
+  double fluxtemp1,fluxtemp0,speedtemp;
+
+  a_leftm  = q_leftm[0];
+  d_leftm = q_leftm[1];
+  z_leftm  = q_leftm[2];
+  h_leftm  = q_leftm[3];
+  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  = q_leftp[3];
+  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  = q_rightm[3];
+  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  = q_rightp[3];
+  u_rightp  = q_rightp[4];
+  b_rightp  = q_rightp[5];
+  w_rightp  = h_rightp+z_rightp;
+
+  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] = flux_formula0(d_leftm, u_leftm, h_leftm, b_leftm, g);
+  flux_left[1] = flux_formula1(d_leftm, u_leftm, h_leftm, b_leftm, g);
+
+  flux_right[0] = flux_formula0(d_leftp, u_leftp, h_leftp, b_leftp, g);
+  flux_right[1] = flux_formula1(d_leftp, u_leftp, h_leftp, b_leftp, g);
+
+  // Flux computation for left hand side
+  compute_flux(flux_left, flux_right, edgeflux);
+  fluxtemp0 = edgeflux[0];
+  fluxtemp1 = edgeflux[1];
+
+  // Flux formulas for right hand side
+  flux_left[0] = flux_formula0(d_rightm, u_rightm, h_rightm, b_rightm, g);
+  flux_left[1] = flux_formula1(d_rightm, u_rightm, h_rightm, b_rightm, g);
+
+  flux_right[0] = flux_formula0(d_rightp, u_rightp, h_rightp, b_rightp, g);
+  flux_right[1] = flux_formula1(d_rightp, u_rightp, h_rightp, b_rightp, g);
+
+  // Flux computation for right hand side
+  compute_flux(flux_left, flux_right, edgeflux);
+  edgeflux[0] = edgeflux[0]-fluxtemp0;
+  edgeflux[1] = edgeflux[1]-fluxtemp1;
+
+  // Forcing terms  
+  edgeflux[1] -= forcingterm(h_rightm, h_leftp, b_rightm, b_leftp, z_rightm, z_leftp, g);
+
+  // 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;
+  }
+
+  edgeflux[0] = 100.0;
+  edgeflux[1] = -10.0;
+
+  return 0;
+}
+
+
+// Computational function for flux computation
+double _compute_fluxes_pipe_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;
+  int k, ki, n, m, nm=0;
+  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_pipe(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;
+
+}
+
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+
+
+//------------------------------------------------
+// New velocity based compute fluxes
+//------------------------------------------------
+
+PyObject *compute_fluxes_pipe_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_pipe_ext.c: compute_fluxes_pipe_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_pipe_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_pipe_ext", compute_fluxes_pipe_ext, METH_VARARGS, "Print out"},
+  {NULL}
+};
+
+/* // Module initialisation */
+void initpipe_domain_ext(void){
+  Py_InitModule("pipe_domain_ext", MethodTable);
+  import_array();
+}
Index: trunk/anuga_work/development/2010-projects/anuga_1d/pipe/test_pipe.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/pipe/test_pipe.py	(revision 8153)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/pipe/test_pipe.py	(revision 8153)
@@ -0,0 +1,129 @@
+import os
+import random
+from math import sqrt, pow, pi ,sin, cos
+import numpy
+
+from anuga_1d.pipe.pipe_domain import *
+
+from anuga_1d.config import g, epsilon
+
+
+
+print "Pipe 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 = numpy.zeros(len(x),numpy.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 = numpy.zeros(N,numpy.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 = numpy.zeros(len(x),numpy.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 = numpy.zeros(N+1,numpy.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
+Z          = domain.quantities['elevation'].vertex_values
+Stage      = HeightQ + Z
+
+
+
+from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+
+hold(False)
+   
+plot1 = subplot(211)
+plot(X.flat,Z.flat, X.flat,Stage.flat)
+ 
+plot1.set_ylim([-1,35])
+xlabel('Position')
+ylabel('Stage')
+legend(('Bed', 'Stage'), 'upper right', shadow=True)
+plot2 = subplot(212)
+plot(X.flat,VelocityQ.flat)
+plot2.set_ylim([-10,10])
+    
+xlabel('Position')
+ylabel('Velocity')
+   
+show()
