Index: /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
===================================================================
--- /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5725)
+++ /anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c	(revision 5725)
@@ -0,0 +1,284 @@
+#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;} 
+
+
+
+
+
+//Innermost flux function (using w=z+h)
+int _flux_function_wellbalanced(double *q_left, double *q_right,
+        double z_left, double z_right,
+		double normals, double g, double epsilon,
+		double *edgeflux, double *max_speed) {
+		
+		int i;
+		double ql[2], qr[2], flux_left[2], flux_right[2];
+		double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
+		double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
+		double s_max, s_min, denom;
+		
+		if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
+		ql[0] = q_left[0];
+		ql[1] = q_left[1];
+		ql[1] = ql[1]*normals;
+	
+		qr[0] = q_right[0];
+		qr[1] = q_right[1];
+		qr[1] = qr[1]*normals;
+	  
+		//z = (z_left+z_right)/2.0;
+		z_star = max(z_left, z_right);		              											//equation(7)  
+				   
+		// Compute speeds in x-direction
+		w_left = ql[0];
+		h_left = w_left-z_left;		                                                                //This is used for computing the edgeflux[1].
+		h_left_star = max(0, w_left-z_star);														//(8)
+		uh_left = ql[1];
+		if (h_left_star < epsilon) {
+			u_left = 0.0; h_left_star = 0.0;
+		}
+		else {
+			u_left = uh_left/h_left_star;
+		}
+	
+		w_right = qr[0];
+		h_right = w_right-z_right;                                                                                                                                                                                                     
+		h_right_star = max(0, w_right-z_star);														//(9)
+		uh_right = qr[1];
+		if (h_right_star < epsilon) {
+			u_right = 0.0; h_right_star = 0.0; 
+		}
+		else {
+			u_right = uh_right/h_right_star;
+		}
+  
+		//printf("u_left = %g u_right = %g \n",u_left, u_right);
+		soundspeed_left = sqrt(g*h_left_star);
+		soundspeed_right = sqrt(g*h_right_star);
+		
+		s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+		if (s_max < 0.0) s_max = 0.0;
+	
+		s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+		if (s_min > 0.0) s_min = 0.0;
+		
+		
+		// Flux formulas
+		flux_left[0] = u_left*h_left_star;
+		flux_left[1] = u_left*u_left*h_left_star + 0.5*g*h_left_star*h_left_star;
+
+		flux_right[0] = u_right*h_right_star;
+		flux_right[1] = u_right*u_right*h_right_star + 0.5*g*h_right_star*h_right_star;
+
+		// Flux computation
+		denom = s_max-s_min;
+		if (denom < epsilon) {
+			for (i=0; i<2; i++) edgeflux[i] = 0.0;
+			*max_speed = 0.0;
+		} else {
+			edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
+			edgeflux[0] += s_max*s_min*(h_right_star-h_left_star);                   //s_max*s_min*(qr[0]-ql[0]);
+			edgeflux[0] /= denom;
+			edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
+			edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
+			edgeflux[1] /= denom;
+			edgeflux[1] *= normals;
+    
+		//edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star;
+		//edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star;
+		edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star;           //what about the h_right and h_right_star??????????????????????
+		// Maximal wavespeed
+        *max_speed = max(fabs(s_max), fabs(s_min));
+		}  
+		//printf("max_speed = %g \n",*max_speed);
+    return 0;		
+	}
+		
+		
+	
+	
+// Computational function for flux computation
+double _compute_fluxes_ext_wellbalanced(double timestep,
+		double epsilon,
+		double g,
+		long* neighbours,
+		long* neighbour_vertices,
+		double* normals,
+		double* areas,
+		double* stage_edge_values,
+		double* xmom_edge_values,
+		double* bed_edge_values,
+		double* stage_boundary_values,
+		double* xmom_boundary_values,
+		double* stage_explicit_update,
+		double* xmom_explicit_update,
+		int number_of_elements,
+		double* max_speed_array) {
+		
+		double flux[2], ql[2], qr[2], edgeflux[2];
+		double zl, zr, max_speed, normal;
+		int k, i, ki, n, m, nm=0;
+		
+		for (k=0; k<number_of_elements; k++) {
+			flux[0] = 0.0;
+			flux[1] = 0.0;
+			
+			for (i=0; i<2; i++) {
+				ki = k*2+i;
+				
+				ql[0] = stage_edge_values[ki];
+				ql[1] = xmom_edge_values[ki];
+				zl = bed_edge_values[ki];
+				
+				n = neighbours[ki];
+				if (n<0) {
+					m = -n-1;
+					qr[0] = stage_boundary_values[m];
+					qr[1] = xmom_boundary_values[m];
+					zr = zl;
+				} else {
+					m = neighbour_vertices[ki];
+					nm = n*2+m;
+					qr[0] = stage_edge_values[nm];
+					qr[1] = xmom_edge_values[nm];
+					zr = bed_edge_values[nm];				
+				}
+				
+				normal = normals[ki];
+				_flux_function_wellbalanced(ql, qr, zl, zr, normal, g, epsilon, 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*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
+					if (n>=0) {
+						timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
+					}
+				}
+            } // End edge i (and neighbour n)
+			flux[0] /= areas[k];
+			stage_explicit_update[k] = flux[0];
+			flux[1] /= areas[k];
+			xmom_explicit_update[k] = flux[1];
+			
+			//Keep track of maximal speeds
+			max_speed_array[k]=max_speed;
+		}
+		printf("timestep = %f \n",timestep);
+		return timestep;	
+	}
+	
+	
+	
+
+
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+PyObject *compute_fluxes_ext_wellbalanced(PyObject *self, PyObject *args) {
+  
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_vertices,
+	*normals, 
+	*areas,
+	*stage_edge_values,
+	*xmom_edge_values,
+	*bed_edge_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+    
+  double timestep, epsilon, g;
+  int number_of_elements;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
+			&timestep,
+			&epsilon,
+			&g,
+			&neighbours,
+			&neighbour_vertices,
+			&normals,
+			&areas,
+			&stage_edge_values,
+			&xmom_edge_values,
+			&bed_edge_values,
+			&stage_boundary_values,
+			&xmom_boundary_values,
+			&stage_explicit_update,
+			&xmom_explicit_update,
+			&number_of_elements,
+			&max_speed_array)) {
+    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
+    return NULL;
+  }
+
+ 
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_ext_wellbalanced(timestep,
+				 epsilon,
+				 g,
+				 (long*) neighbours -> data,
+				 (long*) neighbour_vertices -> data,
+				 (double*) normals -> data,
+				 (double*) areas -> data,
+				 (double*) stage_edge_values -> data,
+				 (double*) xmom_edge_values -> data,
+				 (double*) bed_edge_values -> data,
+				 (double*) stage_boundary_values -> data,
+				 (double*) xmom_boundary_values -> data,
+				 (double*) stage_explicit_update -> data,
+				 (double*) xmom_explicit_update -> data,
+				 number_of_elements,
+				 (double*) max_speed_array -> data);
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+
+static struct PyMethodDef MethodTable[] = {
+  {"compute_fluxes_ext_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initcomp_flux_ext_wellbalanced(void){
+  Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
+  import_array();
+}
+	
+	
Index: /anuga_work/development/anuga_1d/dam_h_elevation.py
===================================================================
--- /anuga_work/development/anuga_1d/dam_h_elevation.py	(revision 5725)
+++ /anuga_work/development/anuga_1d/dam_h_elevation.py	(revision 5725)
@@ -0,0 +1,171 @@
+import os
+from math import sqrt, sin, cos, pi, exp
+from shallow_water_domain_new import *
+from Numeric import zeros, Float
+#from analytic_dam_sudi import AnalyticDam
+
+h0=5.0
+h1=10.0
+
+#analytical_sol=AnalyticDam(h0,h1)
+
+"""
+def newLinePlot(title='Simple Plot'):
+    import Gnuplot
+    gg=Gnuplot.Gnuplot(persist=0)
+    gg.title(title)
+    gg('set data style linespoints')
+    gg.xlabel('x')
+    gg.ylabel('y')
+    return gg
+
+def linePlot(gg, x1, y1, x2, y2):
+    import Gnuplot
+    plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints")
+    plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3")
+    gg.plot(plot1, plot2)
+"""
+
+
+print "TEST 1D-SOLUTION I"
+
+L=2000.0
+N=800
+
+cell_len=L/N
+
+points=zeros(N+1, Float)
+for i in range(N+1):
+    points[i]=i*cell_len
+
+domain=Domain(points)
+
+domain.default_order = 2
+domain.default_time_order = 2
+domain.cfl = 1.0
+domain.limiter = "vanleer"
+
+
+
+def stage(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        if x[i]<=1000.0:
+            y[i]=10.0
+        else:
+            y[i]=5.0
+    return y
+
+def stage_perturb(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        if 1100.0 <= x[i] <=1200.0:
+            y[i]=10.0
+        else:
+            y[i]=10.0
+    return y
+
+def stage_sincos(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        y[i] = 5+exp(cos(pi*x[i]*0.001))+sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
+    return y
+
+
+def elevation_box(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        if 500.0 <= x[i] <= 1500.0:
+            y[i] = 2.5
+        else:
+            y[i] = 0.0
+    return y
+
+def elevation_parabol(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        if 1400.0 <= x[i] <= 1600.0:
+            y[i] = 1.25*(1+cos(0.01*pi*(x[i]-500)))
+        else:
+            y[i] = 0.0
+    return y
+
+
+def elevation_sincos(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        y[i] = sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
+    return y
+
+
+
+def xmom_sincos(x):
+    y=zeros(len(x), Float)
+    for i in range(len(x)):
+        y[i] = sin(cos(pi*x[i]*0.001))
+    return y
+
+domain.set_quantity('stage',stage_perturb)
+domain.set_quantity('elevation', elevation_parabol)
+#domain.set_quantity('xmomentum', xmom_sincos)
+domain.order=domain.default_order
+print "domain order", domain.order
+
+domain.set_boundary({'exterior':Reflective_boundary(domain)})
+
+X=domain.vertices
+C=domain.centroids
+#plot1x=newLinePlot("Height")
+#plot2x=newLinePlot("Momentum")
+
+
+import time
+t0=time.time()
+yieldstep=0.5 #30.0
+finaltime=10.0 #20.0
+print "integral", domain.quantities['stage'].get_integral()
+for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+    domain.write_time()
+    print "integral", domain.quantities['stage'].get_integral()
+    if t>0.0:
+        HeightQ=domain.quantities['stage'].vertex_values
+        MomentumQ=domain.quantities['xmomentum'].vertex_values
+        ElevationQ=domain.quantities['elevation'].vertex_values
+        #h, uh=analytical_sol(X.flat, domain.time)
+        #linePlot(plot1x, X, HeightQ, X, h)
+        #linePlot(plot2x, X, MomentumQ, X, uh)
+        #print "press return"
+        #pass
+    
+        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
+        #print 'Test1'
+        hold(False)
+        #print 'test 2'
+        plot1 = subplot(211)
+        #print 'test 3'
+    
+        plot(X,ElevationQ,X,HeightQ)
+        #print 'Test4'
+        plot1.set_ylim([-1,12])
+        xlabel('Position')
+        ylabel('Stage')
+        legend( ('Bed Elevation', 'Numerical Solution'),
+           'upper right', shadow=False)
+        #legend(('Analytical Solution', 'Numerical Solution'),
+        #   'lower right', shadow=False)
+        plot2 = subplot(212)
+        plot(X,MomentumQ)
+        #plot2.set_ylim([-5,35])
+        legend( ('Numerical Solution', 'for momentum'),
+           'upper right', shadow=False)
+    
+        xlabel('Position')
+        ylabel('Xmomentum')
+        
+        file = "dam_h_"
+        #file += str(number_of_cells[i])
+        file += ".eps"
+        #savefig(file)
+        show()
+
+print 'That took %.2f seconds'%(time.time()-t0)
Index: /anuga_work/development/anuga_1d/shallow_water_domain_new.py
===================================================================
--- /anuga_work/development/anuga_1d/shallow_water_domain_new.py	(revision 5725)
+++ /anuga_work/development/anuga_1d/shallow_water_domain_new.py	(revision 5725)
@@ -0,0 +1,932 @@
+"""Class Domain -
+1D interval domains for finite-volume computations of
+the shallow water wave equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the Shallow Water Wave Equation
+
+
+U_t + E_x = S
+
+where
+
+U = [w, uh]
+E = [uh, u^2h + gh^2/2]
+S represents source terms forcing the system
+(e.g. gravity, friction, wind stress, ...)
+
+and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
+
+The quantities are
+
+symbol    variable name    explanation
+x         x                horizontal distance from origin [m]
+z         elevation        elevation of bed on which flow is modelled [m]
+h         height           water height above z [m]
+w         stage            absolute water level, w = z+h [m]
+u                          speed in the x direction [m/s]
+uh        xmomentum        momentum in the x direction [m^2/s]
+
+eta                        mannings friction coefficient [to appear]
+nu                         wind stress coefficient [to appear]
+
+The conserved quantities are w, uh
+
+For details see e.g.
+Christopher Zoppou and Stephen Roberts,
+Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
+Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
+
+
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2006
+"""
+
+
+from domain import *
+Generic_Domain = Domain #Rename
+
+#Shallow water domain
+class Domain(Generic_Domain):
+
+    def __init__(self, coordinates, boundary = None, tagged_elements = None):
+
+        conserved_quantities = ['stage', 'xmomentum']
+        other_quantities = ['elevation', 'height', 'velocity']
+        Generic_Domain.__init__(self, coordinates, boundary,
+                                conserved_quantities, other_quantities,
+                                tagged_elements)
+        
+        from config import minimum_allowed_height, g
+        self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
+
+        #forcing terms not included in 1d domain ? Why?
+        self.forcing_terms.append(gravity)
+        
+        #Realtime visualisation
+        self.visualiser = None
+        self.visualise  = False
+        self.visualise_color_stage = False
+        self.visualise_stage_range = 1.0
+        self.visualise_timer = True
+        self.visualise_range_z = None
+        
+        #Stored output
+        self.store = True
+        self.format = 'sww'
+        self.smooth = True
+
+        #Evolve parametrs
+        self.cfl = 1.0
+        
+        #Reduction operation for get_vertex_values
+        from util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.quantities_to_be_stored = ['stage', 'xmomentum', 'elevation']
+
+
+
+
+    def set_quantities_to_be_stored(self, q):
+        """Specify which quantities will be stored in the sww file.
+
+        q must be either:
+          - the name of a quantity
+          - a list of quantity names
+          - None
+
+        In the two first cases, the named quantities will be stored at each
+        yieldstep
+        (This is in addition to the quantities elevation and friction)  
+        If q is None, storage will be switched off altogether.
+        """
+
+
+        if q is None:
+            self.quantities_to_be_stored = []            
+            self.store = False
+            return
+
+        if isinstance(q, basestring):
+            q = [q] # Turn argument into a list
+
+        #Check correcness    
+        for quantity_name in q:
+            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
+            assert quantity_name in self.conserved_quantities, msg 
+        
+        self.quantities_to_be_stored = q
+        
+
+    def initialise_visualiser(self,scale_z=1.0,rect=None):
+        #Realtime visualisation
+        if self.visualiser is None:
+            from realtime_visualisation_new import Visualiser
+            self.visualiser = Visualiser(self,scale_z,rect)
+            self.visualiser.setup['elevation']=True
+            self.visualiser.updating['stage']=True
+        self.visualise = True
+        if self.visualise_color_stage == True:
+            self.visualiser.coloring['stage'] = True
+            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
+        print 'initialise visualiser'
+        print self.visualiser.setup
+        print self.visualiser.updating
+
+    def check_integrity(self):
+        Generic_Domain.check_integrity(self)
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "stage"'
+        assert self.conserved_quantities[0] == 'stage', msg
+        msg = 'Second conserved quantity must be "xmomentum"'
+        assert self.conserved_quantities[1] == 'xmomentum', msg
+
+    def extrapolate_second_order_sw(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_sw(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_C_long(self) #compute_fluxes_C_wellbalanced(self)   #compute_fluxes_C_long(self) 		#compute_fluxes(self)    #compute_fluxes_C_short(self) 
+
+    def compute_timestep(self):
+        #Call correct module function
+        compute_timestep(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges(self)
+        
+    def evolve(self, yieldstep = None, finaltime = None,
+               skip_initial_step = False):
+        """Specialisation of basic evolve method from parent class
+        """
+
+        self.distribute_to_vertices_and_edges()
+        
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
+                                       skip_initial_step):
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+    def initialise_storage(self):
+        """Create and initialise self.writer object for storing data.
+        Also, save x and bed elevation
+        """
+
+        import data_manager
+
+        #Initialise writer
+        self.writer = data_manager.get_dataobject(self, mode = 'w')
+
+        #Store vertices and connectivity
+        self.writer.store_connectivity()
+
+
+    def store_timestep(self, name):
+        """Store named quantity and time.
+
+        Precondition:
+           self.write has been initialised
+        """
+        self.writer.store_timestep(name)
+
+
+#=============== End of Shallow Water Domain ===============================
+def compute_timestep(domain):
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                qr[0] = stage[n, m]
+                qr[1] =  xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+
+            #if domain.split == False:
+            #    edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            #elif domain.split == True:
+            #    edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
+            #Update optimal_timestep
+            #try:
+            #    timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            #except ZeroDivisionError:
+            #    pass
+
+    domain.timestep = timestep
+#
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+#The compute_fluxes is involving the flux_function above.
+def compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+#    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+
+        Stage.explicit_update[k] = flux[0]
+        Xmom.explicit_update[k] =  flux[1]
+        
+
+    domain.timestep = timestep
+    #print domain.quantities['stage'].centroid_values
+
+
+# Compute flux definition
+def compute_fluxes_C_long(domain):
+	from Numeric import zeros, Float
+	import sys
+
+	
+	timestep = float(sys.maxint)
+	epsilon = domain.epsilon
+	g = domain.g
+	neighbours = domain.neighbours
+	neighbour_vertices = domain.neighbour_vertices
+	normals = domain.normals
+	areas = domain.areas
+	stage_edge_values = domain.quantities['stage'].vertex_values
+	xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+	bed_edge_values = domain.quantities['elevation'].vertex_values
+	stage_boundary_values = domain.quantities['stage'].boundary_values
+	xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+	stage_explicit_update = domain.quantities['stage'].explicit_update
+	xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+	number_of_elements = len(stage_edge_values)
+	max_speed_array = domain.max_speed_array
+	        
+        from comp_flux_ext import compute_fluxes_ext
+	
+	domain.timestep = compute_fluxes_ext(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+
+# Compute flux definition
+def compute_fluxes_C_short(domain):
+	from Numeric import zeros, Float
+	import sys
+
+	
+	timestep = float(sys.maxint)
+
+        stage = domain.quantities['stage']
+        xmom = domain.quantities['xmomentum']
+        bed = domain.quantities['elevation']
+
+
+        from comp_flux_ext import compute_fluxes_ext_short
+
+ 	domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+
+
+def compute_fluxes_C_wellbalanced(domain):
+	from Numeric import zeros, Float
+	import sys
+
+	
+	timestep = float(sys.maxint)
+	epsilon = domain.epsilon
+	g = domain.g
+	neighbours = domain.neighbours
+	neighbour_vertices = domain.neighbour_vertices
+	normals = domain.normals
+	areas = domain.areas
+	stage_edge_values = domain.quantities['stage'].vertex_values
+	xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+	bed_edge_values = domain.quantities['elevation'].vertex_values
+	stage_boundary_values = domain.quantities['stage'].boundary_values
+	xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+	stage_explicit_update = domain.quantities['stage'].explicit_update
+	xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+	number_of_elements = len(stage_edge_values)
+	max_speed_array = domain.max_speed_array
+	        
+        from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced                  #from comp_flux_ext import compute_fluxes_ext
+	
+	domain.timestep = compute_fluxes_ext_wellbalanced(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+####################################
+
+
+
+
+
+
+# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+
+def distribute_to_vertices_and_edges(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    #from config import optimised_gradient_limiter
+
+    #Remove very thin layers of water
+    protect_against_infinitesimal_and_negative_heights(domain)
+    
+
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            #print "add extrapolate second order to shallow water"
+            #if name != 'height':
+            Q.extrapolate_second_order()
+            #Q.limit()
+        else:
+            raise 'Unknown order'
+
+    
+def protect_against_infinitesimal_and_negative_heights(domain):
+    """Protect against infinitesimal heights and associated high velocities
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    xmomc = domain.quantities['xmomentum'].centroid_values
+#    ymomc = domain.quantities['ymomentum'].centroid_values
+    hc = wc - zc  #Water depths at centroids
+
+    zv = domain.quantities['elevation'].vertex_values
+    wv = domain.quantities['stage'].vertex_values
+    hv = wv-zv
+    xmomv = domain.quantities['xmomentum'].vertex_values
+    #remove the above two lines and corresponding code below
+
+    #Update
+    for k in range(domain.number_of_elements):
+
+        if hc[k] < domain.minimum_allowed_height:
+            #Control stage
+            if hc[k] < domain.epsilon:
+                wc[k] = zc[k] # Contain 'lost mass' error
+                wv[k,0] = zv[k,0]
+                wv[k,1] = zv[k,1]
+                
+            xmomc[k] = 0.0
+            
+        
+    
+def h_limiter(domain):
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    hvbar = zeros(hv.shape, Float) #h-limited values
+
+    #Find min and max of this and neighbour's centroid values
+    hmax = zeros(hc.shape, Float)
+    hmin = zeros(hc.shape, Float)
+
+    for k in range(N):
+        hmax[k] = hmin[k] = hc[k]
+        #for i in range(3):
+        for i in range(2):   
+            n = domain.neighbours[k,i]
+            if n >= 0:
+                hn = hc[n] #Neighbour's centroid value
+
+                hmin[k] = min(hmin[k], hn)
+                hmax[k] = max(hmax[k], hn)
+
+
+    #Diffences between centroids and maxima/minima
+    dhmax = hmax - hc
+    dhmin = hmin - hc
+
+    #Deltas between vertex and centroid values
+    dh = zeros(hv.shape, Float)
+    #for i in range(3):
+    for i in range(2):
+        dh[:,i] = hv[:,i] - hc
+
+    #Phi limiter
+    for k in range(N):
+
+        #Find the gradient limiter (phi) across vertices
+        phi = 1.0
+        #for i in range(3):
+        for i in range(2):
+            r = 1.0
+            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
+            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
+
+            phi = min( min(r*beta_h, 1), phi )
+
+        #Then update using phi limiter
+        #for i in range(3):
+        for i in range(2):
+            hvbar[k,i] = hc[k] + phi*dh[k,i]
+
+    return hvbar
+
+def balance_deep_and_shallow(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters limiting using w, and stage computed by
+    gradient-limiters limiting using h (h-limiter).
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+
+    The h-limiter is always applied irrespective of the order.
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Limit h
+    hvbar = h_limiter(domain)
+
+    for k in range(domain.number_of_elements):
+        
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]))#,
+#                 abs(zv[k,2]-zc[k]))
+
+
+        hmin = min( hv[k,:] )
+
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed
+            alpha = 1.0
+
+        alpha  = 0.0
+	
+	if alpha < 1:
+
+            #for i in range(3):
+            for i in range(2):
+                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+#            ymomc = domain.quantities['ymomentum'].centroid_values
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+#           ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
+#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+
+
+###############################################
+#Boundaries - specific to the shallow water wave equation
+class Reflective_boundary(Boundary):
+    """Reflective boundary returns same conserved quantities as
+    those present in its neighbour volume but reflected.
+
+    This class is specific to the shallow water equation as it
+    works with the momentum quantities assumed to be the second
+    and third conserved quantities.
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for reflective boundary'
+            raise msg
+
+        #Handy shorthands
+        self.normals = domain.normals
+        self.stage   = domain.quantities['stage'].vertex_values
+        self.xmom    = domain.quantities['xmomentum'].vertex_values
+
+        from Numeric import zeros, Float
+        #self.conserved_quantities = zeros(3, Float)
+        self.conserved_quantities = zeros(2, Float)
+
+    def __repr__(self):
+        return 'Reflective_boundary'
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Reflective boundaries reverses the outward momentum
+        of the volume they serve.
+        """
+
+        q = self.conserved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = self.xmom[vol_id, edge_id]
+        normal = self.normals[vol_id,edge_id]
+
+        r = q
+        r[1] = normal*r[1]
+        r[1] = -r[1]
+        r[1] = normal*r[1]
+        q = r
+        #For start interval there is no outward momentum so do not need to
+        #reverse direction in this case
+
+        return q
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, conserved_quantities=None):
+        Boundary.__init__(self)
+
+        if conserved_quantities is None:
+            msg = 'Must specify one value for each conserved quantity'
+            raise msg
+
+        from Numeric import array, Float
+        self.conserved_quantities=array(conserved_quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.conserved_quantities
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.conserved_quantities
+
+
+#########################
+#Standard forcing terms:
+def gravity(domain):
+    """Apply gravitational pull in the presence of bed slope
+    """
+
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+
+    xmom = domain.quantities['xmomentum'].explicit_update
+    stage = domain.quantities['stage'].explicit_update
+#    ymom = domain.quantities['ymomentum'].explicit_update
+
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
+    #h = Stage.edge_values - Elevation.edge_values
+    h = Stage.vertex_values - Elevation.vertex_values
+    b = Elevation.vertex_values
+    w = Stage.vertex_values
+
+    x = domain.get_vertex_coordinates()
+    g = domain.g
+
+    for k in range(domain.number_of_elements):
+#        avg_h = sum( h[k,:] )/3
+        avg_h = sum( h[k,:] )/2
+
+        #Compute bed slope
+        #x0, y0, x1, y1, x2, y2 = x[k,:]
+        x0, x1 = x[k,:]
+        #z0, z1, z2 = v[k,:]
+        b0, b1 = b[k,:]
+
+        w0, w1 = w[k,:]
+        wx = gradient(x0, x1, w0, w1)
+
+        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+        bx = gradient(x0, x1, b0, b1)
+        
+        #Update momentum (explicit update is reset to source values)
+        xmom[k] += -g*bx*avg_h
+        #xmom[k] = -g*bx*avg_h
+        #stage[k] = 0.0
+ 
+ 
+
+def 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
