Index: /anuga_core/source/anuga/shallow_water_balanced/__init__.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/__init__.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/__init__.py	(revision 7734)
@@ -0,0 +1,21 @@
+"""Make directory available as a Python package
+"""
+
+# Add path of package to PYTHONPATH to allow C-extensions to be loaded
+import sys
+sys.path += __path__
+
+# Make selected classes available directly
+from swb_domain import Domain,\
+     Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary, Time_boundary, File_boundary,\
+     Transmissive_momentum_set_stage_boundary,\
+     Dirichlet_discharge_boundary,\
+     Field_boundary,\
+     Transmissive_stage_zero_momentum_boundary,\
+     Transmissive_n_momentum_zero_t_momentum_set_stage_boundary,\
+     AWI_boundary
+
+
+
+ 
Index: /anuga_core/source/anuga/shallow_water_balanced/run_boundary_condition_time.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/run_boundary_condition_time.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/run_boundary_condition_time.py	(revision 7734)
@@ -0,0 +1,72 @@
+
+
+import numpy as num
+from swb_domain import *
+
+
+"""test_boundary_condition_time
+
+This tests that boundary conditions are evaluated
+using the right time from domain.
+"""
+
+# Setup computational domain
+from anuga.abstract_2d_finite_volumes.mesh_factory \
+     import rectangular_cross
+
+#--------------------------------------------------------------
+# Setup computational domain
+#--------------------------------------------------------------
+N = 20
+points, vertices, boundary = rectangular_cross(N, N)
+domain = Domain(points, vertices, boundary)
+
+#--------------------------------------------------------------
+# Setup initial conditions
+#--------------------------------------------------------------
+domain.set_quantity('elevation', 0.0)
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('stage', 0.0)
+
+#--------------------------------------------------------------
+# Setup boundary conditions
+#--------------------------------------------------------------
+# Time dependent boundary
+Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
+
+Br = Reflective_boundary(domain)              # Reflective wall
+
+domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+
+interactive_visualisation = True
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("elevation", zScale =1.0, dynamic=True)    
+    vis.render_quantity_height("stage", zScale =1.0, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+
+for t in domain.evolve(yieldstep = 0.02, finaltime = 20.0):
+    #print domain.get_time()
+    #q = Bt.evaluate()
+    
+    # FIXME (Ole): This test would not have passed in
+    # changeset:5846.
+    #msg = 'Time boundary not evaluated correctly'
+    #assert num.allclose(t, q[0]), msg
+
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+if interactive_visualisation:
+    vis.evolveFinished()    
+    
Index: /anuga_core/source/anuga/shallow_water_balanced/run_step.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/run_step.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/run_step.py	(revision 7734)
@@ -0,0 +1,123 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave which showed diffusive effects of
+thefirst order and standard second order method. Problem resolved if "rk2" timestepping
+and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from swb_domain import *
+
+from math import cos
+import numpy as num
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'wave_'+time
+output_file = 'wave'
+
+#copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+interactive_visualisation = True
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1000.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = False
+#domain.use_centroid_velocities = False
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',-100.0)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)
+
+def stage(x,y):
+    z = num.zeros_like(x)
+    num.putmask(z, x < L/2 , 1.0)
+    return z
+
+domain.set_quantity('stage', stage  )            
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+## Sine wave
+                  f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+if interactive_visualisation:
+    vis.evolveFinished()
+
Index: /anuga_core/source/anuga/shallow_water_balanced/run_wave.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/run_wave.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/run_wave.py	(revision 7734)
@@ -0,0 +1,117 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave which showed diffusive effects of
+thefirst order and standard second order method. Problem resolved if "rk2" timestepping
+and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from swb_domain import *
+
+from math import cos
+import numpy as num
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'wave_'+time
+output_file = 'wave'
+
+#copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+interactive_visualisation = True
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1000.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = False
+#domain.use_centroid_velocities = False
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',-100.0)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)            
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+wave_length = 1000.0
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+## Sine wave
+                  f=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+if interactive_visualisation:
+    vis.evolveFinished()
+
Index: /anuga_core/source/anuga/shallow_water_balanced/run_wave_shelf.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/run_wave_shelf.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/run_wave_shelf.py	(revision 7734)
@@ -0,0 +1,139 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave which showed diffusive effects of
+thefirst order and standard second order method. Problem resolved if "rk2" timestepping
+and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from swb_domain import *
+
+from math import cos
+import numpy as num
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'wave_shelf_'+time
+output_file = 'wave_shelf'
+
+#copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+interactive_visualisation = True
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1000.
+dy = dx
+L = 150000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = False
+#domain.use_centroid_velocities = False
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+depth = -100.0
+def topography(x,y):
+    z = (x*0.0)+depth
+
+    N = len(x)
+    for i in range(N):
+        if  40000.0 < x[i] < 60000.0:
+            z[i] = depth*( 1.0 - (x[i] - 40000.0)/20000.0)
+
+        if x[i] > 58000.0:
+            z[i] = 0.1*depth
+
+        if x[i] > 100000.0:
+            z[i] = 0.1*depth*(1.0 - (x[i] - 75000.0)/50000)
+    return z
+
+domain.set_quantity('elevation',topography)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)            
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+wave_length = 2000.0
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+## Sine wave
+                  f=lambda t: [(amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bw, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+min = 60.0
+hr = 60*min
+
+for t in domain.evolve(yieldstep = min, finaltime = 5*hr):
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+    if domain.get_time()> 2*hr:
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+if interactive_visualisation:
+    vis.evolveFinished()
+
Index: /anuga_core/source/anuga/shallow_water_balanced/swb_boundary_conditions.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/swb_boundary_conditions.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/swb_boundary_conditions.py	(revision 7734)
@@ -0,0 +1,33 @@
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Boundary
+
+class Transmissive_boundary(Boundary):
+    """Transmissive boundary returns same conserved quantities as
+    those present in its neighbour volume.
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for transmissive boundary'
+            raise Exception, msg
+
+        self.domain = domain
+
+    def __repr__(self):
+        return 'Transmissive_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Transmissive boundaries return the edge values
+        of the volume they serve.
+        """
+
+        if self.domain.get_centroid_transmissive_bc() :
+            q = self.domain.get_evolved_quantities(vol_id)
+        else:
+            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+        return q
Index: /anuga_core/source/anuga/shallow_water_balanced/swb_domain_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/swb_domain_ext.c	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/swb_domain_ext.c	(revision 7734)
@@ -0,0 +1,668 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.3):
+//  gcc -c swb_domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
+//  gcc -shared swb_domain_ext.o  -o swb_domain_ext.so
+//
+// or use python compile.py
+//
+// See the module swb_domain.py for more documentation on 
+// how to use this module
+//
+//
+// Stephen Roberts, ANU 2009
+// Ole Nielsen, GA 2004
+
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+//#include "numpy_shim.h"
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+const double pi = 3.14159265358979;
+
+// Computational function for rotation
+int _rotate(double *q, double n1, double n2) {
+  /*Rotate the last  2 coordinates of q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector (n1, n2).
+
+    Result is returned in array 2x1 r
+    To rotate in opposite direction, call rotate with (q, n1, -n2)
+
+    Contents of q are changed by this function */
+
+
+  double q1, q2;
+
+  // Shorthands
+  q1 = q[1];  // x coordinate
+  q2 = q[2];  // y coordinate
+
+  // Rotate
+  q[1] =  n1*q1 + n2*q2;
+  q[2] = -n2*q1 + n1*q2;
+
+  return 0;
+}
+
+// Innermost flux function
+int _flux_function(double *q_inside, double *q_outside,
+		   double n1, double n2, 
+		   double g,
+		   double *edgeflux, double *local_max_speed) 
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', 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
+  */
+
+  int i;
+
+  double z_star, h_inside_star, h_outside_star;
+  double w_inside,  h_inside,  u_inside,  v_inside,  z_inside;
+  double w_outside, h_outside, u_outside, v_outside, z_outside;
+  double s_min, s_max, soundspeed_inside, soundspeed_outside;
+  double epsilon, denom, inverse_denominator;
+
+  // Workspace
+  double cons_inside[3], cons_outside[3], flux_outside[3], flux_inside[3];
+
+  epsilon = 1.0e-15;
+
+  // Setup aliases
+  w_inside = q_inside[0]; 
+  h_inside = q_inside[3];
+  z_inside = q_inside[4];
+  u_inside = q_inside[5];
+  v_inside = q_inside[6];
+
+  w_outside = q_outside[0]; 
+  h_outside = q_outside[3];
+  z_outside = q_outside[4];
+  u_outside = q_outside[5];
+  v_outside = q_outside[6];
+
+
+  // Rotate the velocity terms and update. Setup cons_ to hold conservative rotated
+  // variables 
+  cons_inside[0] = w_inside;
+  cons_inside[1] = u_inside;
+  cons_inside[2] = v_inside;
+
+  _rotate(cons_inside, n1, n2);
+
+  u_inside = cons_inside[1];
+  v_inside = cons_inside[2];
+  cons_inside[1] = u_inside*h_inside;
+  cons_inside[2] = v_inside*h_inside;
+
+  cons_outside[0] = w_outside;
+  cons_outside[1] = u_outside;
+  cons_outside[2] = v_outside;
+
+  _rotate(cons_outside, n1, n2);
+
+  u_outside = cons_outside[1];
+  v_outside = cons_outside[2];
+  cons_outside[1] = u_outside*h_outside;
+  cons_outside[2] = v_outside*h_outside;
+
+
+
+  // Deal with discontinuous z
+  z_star = max(z_inside, z_outside);
+  h_inside_star  = max(0.0, w_inside-z_star);
+  h_outside_star = max(0.0, w_outside-z_star);
+	           
+
+  // Maximal and minimal wave speeds
+  soundspeed_inside  = sqrt(g*h_inside_star);
+  soundspeed_outside = sqrt(g*h_outside_star);  
+  
+  s_max = max(u_inside + soundspeed_inside, u_outside + soundspeed_outside);
+  s_max = max(0.0, s_max);
+
+  s_min = min(u_inside - soundspeed_inside, u_outside - soundspeed_outside);
+  s_min = min(0.0, s_min);
+  
+  // Flux formulas
+  flux_inside[0] = u_inside*h_inside_star;
+  flux_inside[1] = u_inside*u_inside*h_inside_star + 0.5*g*h_inside_star*h_inside_star;
+  flux_inside[2] = u_inside*v_inside*h_inside_star;
+
+  flux_outside[0] = u_outside*h_outside_star;
+  flux_outside[1] = u_outside*u_outside*h_outside_star + 0.5*g*h_outside_star*h_outside_star;
+  flux_outside[2] = u_outside*v_outside*h_outside_star;
+
+  // Flux computation
+  denom = s_max - s_min;
+
+  memset(edgeflux, 0, 3*sizeof(double));
+  
+  if (denom < epsilon) 
+    {     
+      *local_max_speed = 0.0;
+      //printf("here %g\n",h_inside);
+      // Add in edge pressure term due to discontinuous bed
+      if ( h_inside > 0.0 ) edgeflux[1] = 0.5*g*h_inside*h_inside ;
+    } 
+  else 
+    {
+      inverse_denominator = 1.0/denom;
+      for (i = 0; i < 3; i++) 
+	{
+	  edgeflux[i] = s_max*flux_inside[i] - s_min*flux_outside[i];
+	  edgeflux[i] += s_max*s_min*(cons_outside[i] - cons_inside[i]);
+	  edgeflux[i] *= inverse_denominator;
+	}
+
+      // Balance the pressure term
+      edgeflux[1] += 0.5*g*h_inside*h_inside - 0.5*g*h_inside_star*h_inside_star;
+      
+    }
+  
+// Maximal wavespeed
+  *local_max_speed = max(fabs(s_max), fabs(s_min));
+  //printf("local speed % g  h_inside %g \n",*local_max_speed, h_inside);
+  // Rotate back
+  _rotate(edgeflux, n1, -n2);
+  
+  return 0;
+}
+
+
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+
+//========================================================================
+// Compute fluxes
+//========================================================================
+
+PyObject *compute_fluxes(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function_central"
+
+    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.
+
+    Python call:
+    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, height, bed, xvel, yvel)
+
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      returns timestep which is the largest step satisfying all volumes.
+
+
+  */
+
+  PyObject 
+    *domain,
+    *stage, 
+    *xmom, 
+    *ymom,
+    *height,
+    *bed,
+    *xvel,
+    *yvel;
+
+    PyArrayObject 
+      *num_neighbours          ,
+      *num_neighbour_edges     ,
+      *num_normals             ,
+      *num_edgelengths         ,      
+      *num_radii               ,   
+      *num_areas               ,
+      *num_tri_full_flag       ,
+      *num_max_speed           ,
+      *num_stage_edge_values   ,
+      *num_xmom_edge_values    ,
+      *num_ymom_edge_values    ,
+      *num_height_edge_values  ,
+      *num_bed_edge_values     ,
+      *num_xvel_edge_values    ,
+      *num_yvel_edge_values    ,
+      *num_stage_boundary_values ,
+      *num_xmom_boundary_values  ,
+      *num_ymom_boundary_values  ,
+      *num_height_boundary_values,
+      *num_bed_boundary_values   ,
+      *num_xvel_boundary_values  ,
+      *num_yvel_boundary_values  ,
+      *num_stage_explicit_update ,
+      *num_xmom_explicit_update  ,
+      *num_ymom_explicit_update ;
+
+    long
+      *neighbours          ,
+      *neighbour_edges     ,
+      *tri_full_flag       ; 
+
+    double
+      *normals             ,
+      *edgelengths         ,      
+      *radii               ,   
+      *areas               ,
+      *max_speed           ,
+      *stage_edge_values   ,
+      *xmom_edge_values    ,
+      *ymom_edge_values    ,
+      *height_edge_values  ,
+      *bed_edge_values     ,
+      *xvel_edge_values    ,
+      *yvel_edge_values    ,
+      *stage_boundary_values ,
+      *xmom_boundary_values  ,
+      *ymom_boundary_values  ,
+      *height_boundary_values,
+      *bed_boundary_values   ,
+      *xvel_boundary_values  ,
+      *yvel_boundary_values  ,
+      *stage_explicit_update ,
+      *xmom_explicit_update  ,
+      *ymom_explicit_update ;
+
+    
+
+    double flux_timestep, g;
+    double local_max_speed, length, inv_area;
+
+    int k, i, m, n;
+    int ki, nm=0, ki2; // Index shorthands
+   
+    double ql[7], qr[7]; // Work arrays for storing local evolving quantities
+    double edgeflux[3];  // Work array for summing up fluxes
+
+    
+    //======================================================================
+    // Convert Python arguments to C
+    //======================================================================
+    if (!PyArg_ParseTuple(args, "dOOOOOOOO", 
+			  &flux_timestep, 
+			  &domain, 
+			  &stage, 
+			  &xmom, 
+			  &ymom, 
+			  &height,
+			  &bed,
+			  &xvel, 
+			  &yvel)) {
+      PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+      return NULL;
+    }
+    
+    //======================================================================
+    // Extract data from python objects
+    //======================================================================
+    g = get_python_double(domain,"g");
+    
+    num_neighbours             = get_consecutive_array(domain, "neighbours");
+    num_neighbour_edges        = get_consecutive_array(domain, "neighbour_edges"); 
+    num_normals                = get_consecutive_array(domain, "normals");
+    num_edgelengths            = get_consecutive_array(domain, "edgelengths");    
+    num_radii                  = get_consecutive_array(domain, "radii");    
+    num_areas                  = get_consecutive_array(domain, "areas");    
+    num_tri_full_flag          = get_consecutive_array(domain, "tri_full_flag");
+    num_max_speed              = get_consecutive_array(domain, "max_speed");
+    
+    num_stage_edge_values      = get_consecutive_array(stage, "edge_values");    
+    num_xmom_edge_values       = get_consecutive_array(xmom, "edge_values");    
+    num_ymom_edge_values       = get_consecutive_array(ymom, "edge_values"); 
+    num_height_edge_values     = get_consecutive_array(height, "edge_values"); 
+    num_bed_edge_values        = get_consecutive_array(bed, "edge_values");       
+    num_xvel_edge_values       = get_consecutive_array(xvel, "edge_values");    
+    num_yvel_edge_values       = get_consecutive_array(yvel, "edge_values");    
+    
+    num_stage_boundary_values  = get_consecutive_array(stage, "boundary_values");    
+    num_xmom_boundary_values   = get_consecutive_array(xmom, "boundary_values");    
+    num_ymom_boundary_values   = get_consecutive_array(ymom, "boundary_values"); 
+    num_height_boundary_values = get_consecutive_array(height, "boundary_values"); 
+    num_bed_boundary_values    = get_consecutive_array(bed, "boundary_values"); 
+    num_xvel_boundary_values   = get_consecutive_array(xvel, "boundary_values"); 
+    num_yvel_boundary_values   = get_consecutive_array(yvel, "boundary_values"); 
+    
+    num_stage_explicit_update  = get_consecutive_array(stage, "explicit_update");    
+    num_xmom_explicit_update   = get_consecutive_array(xmom, "explicit_update");    
+    num_ymom_explicit_update   = get_consecutive_array(ymom, "explicit_update"); 
+    
+    //---------------------------------------------------------------------------
+    
+    neighbours         = (long *) num_neighbours-> data;
+    neighbour_edges    = (long *) num_neighbour_edges-> data;
+    normals            = (double *) num_normals -> data;
+    edgelengths        = (double *) num_edgelengths -> data;
+    radii              = (double *) num_radii -> data;
+    areas              = (double *) num_areas -> data; 
+    tri_full_flag      = (long *) num_tri_full_flag -> data; 
+    max_speed          = (double *) num_max_speed -> data;
+    
+    stage_edge_values       = (double *) num_stage_edge_values -> data;   
+    xmom_edge_values        = (double *) num_xmom_edge_values -> data;
+    ymom_edge_values        = (double *) num_ymom_edge_values -> data;
+    height_edge_values      = (double *) num_height_edge_values -> data;
+    bed_edge_values         = (double *) num_bed_edge_values -> data;
+    xvel_edge_values        = (double *) num_xvel_edge_values -> data;
+    yvel_edge_values        = (double *) num_yvel_edge_values -> data;
+    
+    stage_boundary_values      = (double *) num_stage_boundary_values -> data; 
+    xmom_boundary_values       = (double *) num_xmom_boundary_values -> data;
+    ymom_boundary_values       = (double *) num_ymom_boundary_values -> data;
+    height_boundary_values     = (double *) num_height_boundary_values -> data;
+    bed_boundary_values        = (double *) num_bed_boundary_values -> data;
+    xvel_boundary_values       = (double *) num_xvel_boundary_values -> data;
+    yvel_boundary_values       = (double *) num_yvel_boundary_values -> data;
+    
+    stage_explicit_update      = (double *) num_stage_explicit_update -> data;
+    xmom_explicit_update       = (double *) num_xmom_explicit_update -> data;
+    ymom_explicit_update       = (double *) num_ymom_explicit_update -> data;
+    
+    
+    int number_of_elements = num_stage_edge_values -> dimensions[0];
+
+    //======================================================================
+    // Flux computation routine and update the explicit update arrays 
+    //======================================================================
+
+    // Set explicit_update to zero for all conserved_quantities.
+    // This assumes compute_fluxes called before forcing terms
+
+    memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
+    memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
+    memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));    
+  
+    // For all triangles
+    for (k = 0; k < number_of_elements; k++) 
+      {  
+	max_speed[k] = 0.0;
+
+	// Loop through neighbours and compute edge flux for each  epsilon
+	for (i = 0; i < 3; i++) 
+	  {
+	    ki = k*3 + i; // Linear index to edge i of triangle k
+      
+	    // Get left hand side values from triangle k, edge i 
+	    ql[0] = stage_edge_values[ki];
+	    ql[1] = xmom_edge_values[ki];
+	    ql[2] = ymom_edge_values[ki];
+	    ql[3] = height_edge_values[ki];
+	    ql[4] = bed_edge_values[ki];
+	    ql[5] = xvel_edge_values[ki];
+	    ql[6] = yvel_edge_values[ki];
+
+	    // Get right hand side values either from neighbouring triangle
+	    // or from boundary array (Quantities at neighbour on nearest face).
+	    n = neighbours[ki];
+	    if (n < 0) 
+	      {
+		// Neighbour is a boundary condition
+		m = -n - 1; // Convert negative flag to boundary index
+		
+		qr[0] = stage_boundary_values[m];
+		qr[1] = xmom_boundary_values[m];
+		qr[2] = ymom_boundary_values[m];
+		qr[3] = height_boundary_values[m];
+		qr[4] = bed_boundary_values[m];
+		qr[5] = xvel_boundary_values[m];
+		qr[6] = yvel_boundary_values[m];
+		
+	      } 
+	    else 
+	      {
+		// Neighbour is a real triangle
+		m = neighbour_edges[ki];
+		nm = n*3 + m; // Linear index (triangle n, edge m)
+		
+		qr[0] = stage_edge_values[nm];
+		qr[1] = xmom_edge_values[nm];
+		qr[2] = ymom_edge_values[nm];
+		qr[3] = height_edge_values[nm];
+		qr[4] = bed_edge_values[nm];
+		qr[5] = xvel_edge_values[nm];
+		qr[6] = yvel_edge_values[nm];
+		
+	      }
+	    
+	    // Now we have values for this edge - both from left and right side.
+
+	    // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
+	    ki2 = 2*ki; //k*6 + i*2
+
+	    // Edge flux computation (triangle k, edge i)
+	    _flux_function(ql,
+			   qr,
+			   normals[ki2], 
+			   normals[ki2+1],
+			   g,
+			   edgeflux, 
+			   &local_max_speed);
+	    
+	    
+	    // Multiply edgeflux by edgelength
+	    length = edgelengths[ki];
+	    edgeflux[0] *= length;            
+	    edgeflux[1] *= length;            
+	    edgeflux[2] *= length;                        
+      
+      
+	    // Update triangle k with flux from edge i
+	    stage_explicit_update[k] -= edgeflux[0];
+	    xmom_explicit_update[k]  -= edgeflux[1];
+	    ymom_explicit_update[k]  -= edgeflux[2];
+      
+
+	    // Update flux_timestep based on edge i and possibly neighbour n
+	    if (tri_full_flag[k] == 1) 
+	      {
+		if (local_max_speed > 1.0e-15) 
+		  {
+		    // Calculate safe timestep based on local_max_speed and size of triangle k
+		    // A reduction of timestep based on domain.CFL is applied in update_timestep 
+		    flux_timestep = min(flux_timestep, radii[k]/local_max_speed);
+		    
+		  }
+	      }
+	    
+	  } // End edge i 
+	
+	
+	// Normalise triangle k by area and store for when all conserved
+	// quantities get updated
+	inv_area = 1.0/areas[k];
+	stage_explicit_update[k] *= inv_area;
+	xmom_explicit_update[k]  *= inv_area;
+	ymom_explicit_update[k]  *= inv_area;
+	
+	
+	// Keep track of maximal speeds
+	max_speed[k] = max(max_speed[k],local_max_speed);
+	
+      } // End triangle k
+    
+    
+    //======================================================================
+    // Cleanup
+    //======================================================================
+
+    Py_DECREF(num_neighbours            );
+    Py_DECREF(num_neighbour_edges       );
+    Py_DECREF(num_normals               );
+    Py_DECREF(num_edgelengths           );      
+    Py_DECREF(num_radii                 );   
+    Py_DECREF(num_areas                 );
+    Py_DECREF(num_tri_full_flag         );
+    Py_DECREF(num_max_speed             );
+    Py_DECREF(num_stage_edge_values     );
+    Py_DECREF(num_xmom_edge_values      );
+    Py_DECREF(num_ymom_edge_values      );
+    Py_DECREF(num_height_edge_values    );
+    Py_DECREF(num_bed_edge_values       );
+    Py_DECREF(num_xvel_edge_values      );
+    Py_DECREF(num_yvel_edge_values      );
+    Py_DECREF(num_stage_boundary_values );
+    Py_DECREF(num_xmom_boundary_values  );
+    Py_DECREF(num_ymom_boundary_values  );
+    Py_DECREF(num_height_boundary_values);
+    Py_DECREF(num_bed_boundary_values   );
+    Py_DECREF(num_xvel_boundary_values  );
+    Py_DECREF(num_yvel_boundary_values  );
+    Py_DECREF(num_stage_explicit_update );
+    Py_DECREF(num_xmom_explicit_update  );
+    Py_DECREF(num_ymom_explicit_update  );
+
+    //======================================================================
+    // Return updated flux flux_timestep
+    //======================================================================
+    return Py_BuildValue("d", flux_timestep);
+}
+
+//========================================================================
+// Flux Function
+//========================================================================
+
+PyObject *flux_function(PyObject *self, PyObject *args) {
+  //
+  // Gateway to innermost flux function.
+  // This is only used by the unit tests as the c implementation is
+  // normally called by compute_fluxes in this module.
+
+
+  PyArrayObject *normal, *ql, *qr,  *edgeflux;
+  double g, local_max_speed;
+
+  if (!PyArg_ParseTuple(args, "OOOOd",
+            &normal, &ql, &qr, &edgeflux, &g)) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "swb_domain_ext.c: flux_function could not parse input arguments");
+    return NULL;
+  }
+
+  
+  _flux_function((double*) ql -> data, 
+		 (double*) qr -> data, 
+		 ((double*) normal -> data)[0],
+		 ((double*) normal -> data)[1],          
+		 g,
+		 (double*) edgeflux -> data, 
+		 &local_max_speed);
+  
+  return Py_BuildValue("d", local_max_speed);  
+}
+
+
+//========================================================================
+// Gravity
+//========================================================================
+
+PyObject *gravity(PyObject *self, PyObject *args) {
+  //
+  //  gravity(g, h, v, x, xmom, ymom)
+  //
+  
+  
+  PyArrayObject *h, *v, *x, *xmom, *ymom;
+  int k, N, k3, k6;
+  double g, avg_h, zx, zy;
+  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
+  //double epsilon;
+  
+  if (!PyArg_ParseTuple(args, "dOOOOO",
+			&g, &h, &v, &x,
+			&xmom, &ymom)) {
+    //&epsilon)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
+    return NULL;
+  }
+  
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(h);
+  CHECK_C_CONTIG(v);
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+  
+  N = h -> dimensions[0];
+  for (k=0; k<N; k++) {
+    k3 = 3*k;  // base index
+    
+    // Get bathymetry
+    z0 = ((double*) v -> data)[k3 + 0];
+    z1 = ((double*) v -> data)[k3 + 1];
+    z2 = ((double*) v -> data)[k3 + 2];
+    
+    // Optimise for flat bed
+    // Note (Ole): This didn't produce measurable speed up.
+    // Revisit later
+    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
+    //  continue;
+    //} 
+    
+    // Get average depth from centroid values
+    avg_h = ((double *) h -> data)[k];
+    
+    // Compute bed slope
+    k6 = 6*k;  // base index
+    
+    x0 = ((double*) x -> data)[k6 + 0];
+    y0 = ((double*) x -> data)[k6 + 1];
+    x1 = ((double*) x -> data)[k6 + 2];
+    y1 = ((double*) x -> data)[k6 + 3];
+    x2 = ((double*) x -> data)[k6 + 4];
+    y2 = ((double*) x -> data)[k6 + 5];
+    
+    
+    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+    
+    // Update momentum
+    ((double*) xmom -> data)[k] += -g*zx*avg_h;
+    ((double*) ymom -> data)[k] += -g*zy*avg_h;
+  }
+  
+  return Py_BuildValue("");
+}
+
+
+
+//========================================================================
+// Method table for python module
+//========================================================================
+
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+  {"compute_fluxes_c", compute_fluxes,     METH_VARARGS, "Print out"},
+  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
+  {"flux_function_c",  flux_function,      METH_VARARGS, "Print out"},
+  {NULL, NULL, 0, NULL}
+};
+
+// Module initialisation
+void initswb_domain_ext(void){
+  Py_InitModule("swb_domain_ext", MethodTable);
+
+  import_array(); // Necessary for handling of NumPY structures
+}
Index: /anuga_core/source/anuga/shallow_water_balanced/test_swb_domain_ext.py
===================================================================
--- /anuga_core/source/anuga/shallow_water_balanced/test_swb_domain_ext.py	(revision 7734)
+++ /anuga_core/source/anuga/shallow_water_balanced/test_swb_domain_ext.py	(revision 7734)
@@ -0,0 +1,495 @@
+#!/usr/bin/env python
+
+import unittest, os
+import os.path
+from math import pi, sqrt
+
+from swb_domain import *
+
+from anuga.config import g
+
+import numpy as num
+
+# Get gateway to C implementation of flux function for direct testing
+from swb_domain_ext import flux_function_c
+from swb_domain_ext import gravity_c
+from swb_domain_ext import compute_fluxes_c
+
+
+class Test_swb_domain_ext(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+    def test_compute_fluxes0(self):
+        # Do a full triangle and check that fluxes cancel out for
+        # the constant stage case
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
+        domain.check_integrity()
+        domain.distribute_to_vertices_and_edges()
+        
+        assert num.allclose(domain.neighbours,
+                            [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
+        assert num.allclose(domain.neighbour_edges,
+                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
+        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
+
+        #print ql
+        #print qr
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(2, 2)
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
+
+        assert num.allclose(edgeflux0 + edgeflux, 0.)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
+        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
+        #print ql
+        #print qr
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(3, 0)
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
+
+        assert num.allclose(edgeflux1 + edgeflux, 0.)
+
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
+        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
+        #print ql
+        #print qr
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(0, 1)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
+        assert num.allclose(edgeflux2 + edgeflux, 0.)
+
+        # Scale by edgelengths, add up anc check that total flux is zero
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
+
+        # Now check that compute_flux yields zeros as well
+        w  = domain.quantities['stage']
+        uh = domain.quantities['xmomentum']
+        vh = domain.quantities['ymomentum']
+        h  = domain.quantities['height']
+        z  = domain.quantities['elevation']
+        u  = domain.quantities['xvelocity']
+        v  = domain.quantities['yvelocity']
+        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
+
+        assert num.allclose(timestep,  0.1064794275)
+
+
+        for name in ['stage', 'xmomentum', 'ymomentum']:
+            assert num.allclose(domain.quantities[name].explicit_update[1], 0.0)
+
+    def test_compute_fluxes1(self):
+        #Use values from previous version
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        val0 = 2. + 2.0/3
+        val1 = 4. + 4.0/3
+        val2 = 8. + 2.0/3
+        val3 = 2. + 8.0/3
+
+        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
+                                      [val2, val2, val2], [val3, val3, val3]])
+        domain.check_integrity()
+        domain.set_default_order(1)
+        domain.distribute_to_vertices_and_edges()
+
+
+        w  = domain.quantities['stage']
+        uh = domain.quantities['xmomentum']
+        vh = domain.quantities['ymomentum']
+        h  = domain.quantities['height']
+        z  = domain.quantities['elevation']
+        u  = domain.quantities['xvelocity']
+        v  = domain.quantities['yvelocity']
+
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
+        assert num.allclose(normal, [1, 0])
+
+  
+        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
+        #print ql, val1
+        assert num.allclose(ql, [val1, 0, 0, val1, 0, 0, 0])
+
+        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
+        assert num.allclose(qr, [val2, 0, 0, val2, 0, 0, 0])
+
+        local_max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
+
+        #print edgeflux0
+        #print local_max_speed
+        
+        # Flux across edge in the east direction (as per normal vector)
+        assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
+        assert num.allclose(local_max_speed, 9.21592824046)
+
+        #Flux across edge in the west direction (opposite sign for xmomentum)
+        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
+        assert num.allclose(normal_opposite, [-1, 0])
+
+        max_speed = flux_function_c(normal_opposite, ql, qr, edgeflux, g)
+        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
+
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
+        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
+
+        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
+        assert num.allclose(max_speed, 7.22956891292)
+
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
+        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
+
+        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
+        assert num.allclose(max_speed, 7.22956891292)
+
+        #Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
+
+
+        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
+        assert num.allclose(timestep, 0.0511510624314)
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+
+        assert num.allclose(domain.quantities['stage'].explicit_update,
+                            [0., -0.68218178, -111.77316251, -35.68522449])
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
+                            [-69.68888889, -166.6, 69.68888889, 0])
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update,
+                            [-69.68888889, -35.93333333, 0., 69.68888889])
+
+    def test_compute_fluxes2(self):
+        #Random values, incl momentum
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        domain.set_default_order(1)
+
+        w  = domain.quantities['stage']
+        uh = domain.quantities['xmomentum']
+        vh = domain.quantities['ymomentum']
+        h  = domain.quantities['height']
+        z  = domain.quantities['elevation']
+        u  = domain.quantities['xvelocity']
+        v  = domain.quantities['yvelocity']
+
+        zl = 10    # Assume flat bed
+        val0 = zl + 2. + 2.0/3
+        val1 = zl + 4. + 4.0/3
+        val2 = zl + 8. + 2.0/3
+        val3 = zl + 2. + 8.0/3
+
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+
+        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
+
+        domain.set_quantity('stage', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
+
+        domain.check_integrity()
+        domain.distribute_to_vertices_and_edges()
+
+        
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
+        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
+        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
+
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
+        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
+
+        # Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
+        assert num.allclose(timestep, 0.0529903533177)
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+    # FIXME (Ole): Need test like this for fluxes in very shallow water.
+    def test_compute_fluxes3(self):
+        #Random values, incl momentum
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        domain.set_default_order(1)
+
+
+        w  = domain.quantities['stage']
+        uh = domain.quantities['xmomentum']
+        vh = domain.quantities['ymomentum']
+        h  = domain.quantities['height']
+        z  = domain.quantities['elevation']
+        u  = domain.quantities['xvelocity']
+        v  = domain.quantities['yvelocity']
+        
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl = -3.75    # Assume constant bed (must be less than stage)
+        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+
+
+        domain.set_quantity('stage', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
+
+        domain.check_integrity()
+        domain.distribute_to_vertices_and_edges()        
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
+        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
+
+        #print max_speed, edgeflux0
+        assert num.allclose(edgeflux0, [-10.51926294, 577.81462335, -3.64772873], rtol=1.0e-2)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
+        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
+
+        #print max_speed, edgeflux1
+        assert num.allclose(edgeflux1, [ 5.93945967, 19.14204647,  377.48851296], rtol=1.0e-2)
+        
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
+        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
+
+        #print max_speed, edgeflux2
+        assert num.allclose(edgeflux2, [17.21047971, -192.90578267, -203.05771337], rtol=1.0e-2)
+
+        
+        # Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        timestep = compute_fluxes_c(100.0,domain,w,uh,vh,h,z,u,v)
+
+        assert num.allclose(timestep, 0.0438142267244, rtol=1.0e-2)
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+
+
+    def test_compute_fluxes4(self):
+        #Random values, incl momentum
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        domain.set_default_order(1)
+
+
+        w  = domain.quantities['stage']
+        uh = domain.quantities['xmomentum']
+        vh = domain.quantities['ymomentum']
+        h  = domain.quantities['height']
+        z  = domain.quantities['elevation']
+        u  = domain.quantities['xvelocity']
+        v  = domain.quantities['yvelocity']
+        
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl = -3.75    # Assume constant bed (must be less than stage)
+        domain.set_quantity('elevation', [0.0, -2.0, 0.0 ,0.0], location='centroids') 
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+
+
+        domain.set_quantity('stage', -1.0)
+
+        domain.set_quantity('xmomentum',0.0)
+
+        domain.set_quantity('ymomentum',0.0)
+
+        domain.check_integrity()
+        domain.distribute_to_vertices_and_edges()        
+
+
+        print w.edge_values
+        print u.edge_values
+        print v.edge_values
+        print h.edge_values
+        print z.edge_values
+
+
+        domain.compute_fluxes()
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            print domain.quantities[name].explicit_update[1]
+            #assert num.allclose(total_flux[i],
+            #                    domain.quantities[name].explicit_update[1])
+
+
+#################################################################################
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_swb_domain_ext, 'test')
+    runner = unittest.TextTestRunner(verbosity=1)
+    runner.run(suite)
