Index: /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1026)
+++ /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1027)
@@ -37,4 +37,5 @@
 
          fileName = tempfile.mktemp(".tsh")
+         print fileName
          file = open(fileName,"w")
          file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
@@ -163,3 +164,2 @@
 
 
-
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1026)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1027)
@@ -15,15 +15,15 @@
     Low speeds at center and edges
     """
-    
+
     from math import sqrt, exp, cos, pi
 
     x = array(x)
-    y = array(y)    
+    y = array(y)
 
     N = len(x)
     s = 0*x  #New array
-    
+
     for k in range(N):
-        
+
         r = sqrt(x[k]**2 + y[k]**2)
 
@@ -49,12 +49,12 @@
 
     x = array(x)
-    y = array(y)    
+    y = array(y)
 
     N = len(x)
     a = 0*x  #New array
 
-    for k in range(N):    
-        r = sqrt(x[k]**2 + y[k]**2)    
-    
+    for k in range(N):
+        r = sqrt(x[k]**2 + y[k]**2)
+
         angle = atan(y[k]/x[k])
 
@@ -64,26 +64,26 @@
         #Take normal direction
         angle -= pi/2
-    
-        #Ensure positive radians    
+
+        #Ensure positive radians
         if angle < 0:
             angle += 2*pi
 
         a[k] = angle/pi*180
-        
+
     return a
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_Shallow_Water(unittest.TestCase):
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
 
     def test_rotate(self):
-        normal = array([0.0,-1.0])        
+        normal = array([0.0,-1.0])
 
         q = array([1.0,2.0,3.0])
-      
+
         r = rotate(q, normal, direction = 1)
         assert r[0] == 1
@@ -91,5 +91,5 @@
         assert r[2] == 2
 
-        w = rotate(r, normal, direction = -1)        
+        w = rotate(r, normal, direction = -1)
         assert allclose(w, q)
 
@@ -114,35 +114,35 @@
     def test_flux_constants(self):
         w = 2.0
-       
-        normal = array([1.,0])        
+
+        normal = array([1.,0])
         ql = array([w, 0, 0])
         qr = array([w, 0, 0])
         zl = zr = 0.
         h = w - (zl+zr)/2
-        
+
         flux, max_speed = flux_function(normal, ql, qr, zl, zr)
 
         assert allclose(flux, [0., 0.5*g*h**2, 0.])
         assert max_speed == sqrt(g*h)
-        
+
     #def test_flux_slope(self):
     #    #FIXME: TODO
     #    w = 2.0
-    #   
-    #    normal = array([1.,0])        
+    #
+    #    normal = array([1.,0])
     #    ql = array([w, 0, 0])
     #    qr = array([w, 0, 0])
     #    zl = zr = 0.
     #    h = w - (zl+zr)/2
-    #    
+    #
     #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
     #
     #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
     #    assert max_speed == sqrt(g*h)
-        
+
 
     def test_flux1(self):
         #Use data from previous version of pyvolution
-        normal = array([1.,0])        
+        normal = array([1.,0])
         ql = array([-0.2, 2, 3])
         qr = array([-0.2, 2, 3])
@@ -156,5 +156,5 @@
     def test_flux2(self):
         #Use data from previous version of pyvolution
-        normal = array([0., -1.])        
+        normal = array([0., -1.])
         ql = array([-0.075, 2, 3])
         qr = array([-0.075, 2, 3])
@@ -166,6 +166,6 @@
 
     def test_flux3(self):
-        #Use data from previous version of pyvolution        
-        normal = array([-sqrt(2)/2, sqrt(2)/2])        
+        #Use data from previous version of pyvolution
+        normal = array([-sqrt(2)/2, sqrt(2)/2])
         ql = array([-0.075, 2, 3])
         qr = array([-0.075, 2, 3])
@@ -177,6 +177,6 @@
 
     def test_flux4(self):
-        #Use data from previous version of pyvolution                
-        normal = array([-sqrt(2)/2, sqrt(2)/2])        
+        #Use data from previous version of pyvolution
+        normal = array([-sqrt(2)/2, sqrt(2)/2])
         ql = array([-0.34319278, 0.10254161, 0.07273855])
         qr = array([-0.30683287, 0.1071986, 0.05930515])
@@ -185,5 +185,5 @@
 
         assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
-        assert allclose(max_speed, 1.31414103233)                                
+        assert allclose(max_speed, 1.31414103233)
 
     def test_sw_domain_simple(self):
@@ -196,8 +196,8 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
-        
-        domain = Domain(points, vertices)        
+
+        domain = Domain(points, vertices)
         domain.check_integrity()
 
@@ -207,9 +207,9 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.    
+        assert domain.get_conserved_quantities(0, edge=1) == 0.
 
 
     def test_boundary_conditions(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -227,6 +227,6 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Third'}                                         
-                     
+                     (3, 2): 'Third'}
+
 
         domain = Domain(points, vertices, boundary)
@@ -241,5 +241,5 @@
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -248,5 +248,5 @@
         R = Reflective_boundary(domain)
         domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
-        
+
         domain.update_boundary()
 
@@ -254,5 +254,5 @@
         assert domain.quantities['stage'].boundary_values[0] == 2.5
         assert domain.quantities['stage'].boundary_values[0] ==\
-               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5) 
+               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
         assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
         assert domain.quantities['stage'].boundary_values[2] ==\
@@ -275,5 +275,5 @@
                domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
         assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
-        
+
         #Ymomentum
         assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
@@ -286,5 +286,5 @@
 
     def test_boundary_conditionsII(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -303,6 +303,6 @@
                      (3, 1): 'Second',
                      (3, 2): 'Third',
-                     (0, 1): 'Internal'}                                         
-                     
+                     (0, 1): 'Internal'}
+
 
         domain = Domain(points, vertices, boundary)
@@ -317,5 +317,5 @@
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -333,5 +333,5 @@
         #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]
@@ -344,6 +344,6 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
-        domain = Domain(points, vertices)        
+
+        domain = Domain(points, vertices)
         domain.set_quantity('stage', [[2,2,2], [2,2,2],
                                       [2,2,2], [2,2,2]])
@@ -351,42 +351,42 @@
 
         assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
-        assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])        
+        assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
 
         zl=zr=0. #Assume flat bed
-        
+
         #Flux across right edge of volume 1
-        normal = domain.get_normal(1,0)        
+        normal = domain.get_normal(1,0)
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
         flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
         #Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(2,2)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
         assert allclose(flux + flux0, 0.)
-    
+
         #Flux across upper edge of volume 1
-        normal = domain.get_normal(1,1)        
+        normal = domain.get_normal(1,1)
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
         flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
         #Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(3,0)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
         assert allclose(flux + flux1, 0.)
-        
+
         #Flux across lower left hypotenuse of volume 1
-        normal = domain.get_normal(1,2)        
+        normal = domain.get_normal(1,2)
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
         flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
         #Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(0,1)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
         assert allclose(flux + flux2, 0.)
 
@@ -396,5 +396,5 @@
         e1 = domain.edgelengths[1, 1]
         e2 = domain.edgelengths[1, 2]
-        
+
         assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
 
@@ -410,5 +410,5 @@
     def test_compute_fluxes1(self):
         #Use values from previous version
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -421,5 +421,5 @@
         #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
@@ -427,5 +427,5 @@
         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]])
@@ -435,5 +435,5 @@
 
         #Flux across right edge of volume 1
-        normal = domain.get_normal(1,0)        
+        normal = domain.get_normal(1,0)
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
@@ -444,5 +444,5 @@
 
         #Flux across upper edge of volume 1
-        normal = domain.get_normal(1,1)        
+        normal = domain.get_normal(1,1)
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
@@ -452,5 +452,5 @@
 
         #Flux across lower left hypotenuse of volume 1
-        normal = domain.get_normal(1,2)        
+        normal = domain.get_normal(1,2)
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
@@ -459,5 +459,5 @@
         assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
         assert allclose(max_speed, 7.22956891292)
-        
+
         #Scale, add up and check that compute_fluxes is correct for vol 1
         e0 = domain.edgelengths[1, 0]
@@ -487,16 +487,16 @@
                         [-69.68888889, -166.6, 69.68888889, 0])
         assert allclose(domain.quantities['ymomentum'].explicit_update,
-                        [-69.68888889, -35.93333333, 0., 69.68888889])        
-
-        
+                        [-69.68888889, -35.93333333, 0., 69.68888889])
+
+
         #assert allclose(domain.quantities[name].explicit_update
 
-        
-        
+
+
 
 
     def test_compute_fluxes2(self):
         #Random values, incl momentum
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -509,5 +509,5 @@
         #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
@@ -519,6 +519,6 @@
 
         domain.set_quantity('elevation', zl*ones( (4,3) ))
-        
-        
+
+
         domain.set_quantity('stage', [[val0, val0-1, val0-2],
                                       [val1, val1+1, val1],
@@ -532,7 +532,7 @@
         domain.set_quantity('ymomentum',
                             [[1, -1, 0], [0, -3, 2],
-                             [0, 1, 0], [-1, 2, 2]])                
-
-        
+                             [0, 1, 0], [-1, 2, 2]])
+
+
         domain.check_integrity()
 
@@ -540,5 +540,5 @@
 
         #Flux across right edge of volume 1
-        normal = domain.get_normal(1,0)        
+        normal = domain.get_normal(1,0)
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
@@ -546,5 +546,5 @@
 
         #Flux across upper edge of volume 1
-        normal = domain.get_normal(1,1)        
+        normal = domain.get_normal(1,1)
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
@@ -552,9 +552,9 @@
 
         #Flux across lower left hypotenuse of volume 1
-        normal = domain.get_normal(1,2)        
+        normal = domain.get_normal(1,2)
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
         flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
         #Scale, add up and check that compute_fluxes is correct for vol 1
         e0 = domain.edgelengths[1, 0]
@@ -564,15 +564,15 @@
         total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
 
-        
+
         domain.compute_fluxes()
         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
             assert allclose(total_flux[i],
-                            domain.quantities[name].explicit_update[1])        
+                            domain.quantities[name].explicit_update[1])
         #assert allclose(total_flux, domain.explicit_update[1,:])
-        
+
 
     def test_compute_fluxes3(self):
         #Random values, incl momentum
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -585,5 +585,5 @@
         #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
@@ -594,6 +594,6 @@
         zl=zr=-3.75 #Assume constant bed (must be less than stage)
         domain.set_quantity('elevation', zl*ones( (4,3) ))
-        
-        
+
+
         domain.set_quantity('stage', [[val0, val0-1, val0-2],
                                       [val1, val1+1, val1],
@@ -607,7 +607,7 @@
         domain.set_quantity('ymomentum',
                             [[1, -1, 0], [0, -3, 2],
-                             [0, 1, 0], [-1, 2, 2]])                
-
-        
+                             [0, 1, 0], [-1, 2, 2]])
+
+
         domain.check_integrity()
 
@@ -615,5 +615,5 @@
 
         #Flux across right edge of volume 1
-        normal = domain.get_normal(1,0)        
+        normal = domain.get_normal(1,0)
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
@@ -621,5 +621,5 @@
 
         #Flux across upper edge of volume 1
-        normal = domain.get_normal(1,1)        
+        normal = domain.get_normal(1,1)
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
@@ -627,9 +627,9 @@
 
         #Flux across lower left hypotenuse of volume 1
-        normal = domain.get_normal(1,2)        
+        normal = domain.get_normal(1,2)
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
         flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
         #Scale, add up and check that compute_fluxes is correct for vol 1
         e0 = domain.edgelengths[1, 0]
@@ -638,5 +638,5 @@
 
         total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
-        
+
         domain.compute_fluxes()
         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
@@ -647,5 +647,5 @@
 
     def test_catching_negative_heights(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -658,5 +658,5 @@
         #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
@@ -681,5 +681,5 @@
 
 
-    #####################################################    
+    #####################################################
     def test_initial_condition(self):
         """Test that initial condition is output at time == 0
@@ -688,5 +688,5 @@
         from config import g
         import copy
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -699,8 +699,8 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
 
-        #Set up for a gradient of (3,0) at mid triangle         
+        #Set up for a gradient of (3,0) at mid triangle
         def slope(x, y):
             return 3*x
@@ -725,14 +725,14 @@
             else:
                 assert not allclose(stage, initial_stage)
-                
-      
-        
-
-    #####################################################    
+
+
+
+
+    #####################################################
     def test_gravity(self):
         #Assuming no friction
 
         from config import g
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -745,8 +745,8 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
 
-        #Set up for a gradient of (3,0) at mid triangle         
+        #Set up for a gradient of (3,0) at mid triangle
         def slope(x, y):
             return 3*x
@@ -762,15 +762,15 @@
             assert allclose(domain.quantities[name].explicit_update, 0)
             assert allclose(domain.quantities[name].semi_implicit_update, 0)
-            
-        domain.compute_forcing_terms()                                    
+
+        domain.compute_forcing_terms()
 
         assert allclose(domain.quantities['stage'].explicit_update, 0)
         assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
         assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
-        
-        
+
+
     def test_manning_friction(self):
         from config import g
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -783,8 +783,8 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
 
-        #Set up for a gradient of (3,0) at mid triangle         
+        #Set up for a gradient of (3,0) at mid triangle
         def slope(x, y):
             return 3*x
@@ -797,11 +797,11 @@
         domain.set_quantity('elevation', slope)
         domain.set_quantity('stage', stage)
-        domain.set_quantity('friction', eta)        
+        domain.set_quantity('friction', eta)
 
         for name in domain.conserved_quantities:
             assert allclose(domain.quantities[name].explicit_update, 0)
             assert allclose(domain.quantities[name].semi_implicit_update, 0)
-            
-        domain.compute_forcing_terms()                                    
+
+        domain.compute_forcing_terms()
 
         assert allclose(domain.quantities['stage'].explicit_update, 0)
@@ -811,5 +811,5 @@
         assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
         assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
-        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)        
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
 
         #Create some momentum for friction to work with
@@ -817,16 +817,16 @@
         S = -g * eta**2 / h**(7.0/3)
 
-        domain.compute_forcing_terms()                                   
+        domain.compute_forcing_terms()
         assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
         assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
-        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)        
-
-        #A more complex example 
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
+
+        #A more complex example
         domain.quantities['stage'].semi_implicit_update[:] = 0.0
         domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
-        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0        
-       
+        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
+
         domain.set_quantity('xmomentum', 3)
-        domain.set_quantity('ymomentum', 4)        
+        domain.set_quantity('ymomentum', 4)
 
         S = -g * eta**2 * 5 / h**(7.0/3)
@@ -837,5 +837,5 @@
         assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
         assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
-        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)        
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
 
     def test_constant_wind_stress(self):
@@ -854,5 +854,5 @@
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        
+
         domain = Domain(points, vertices)
 
@@ -860,5 +860,5 @@
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)        
+        domain.set_quantity('friction', 0)
 
         Br = Reflective_boundary(domain)
@@ -870,13 +870,13 @@
         domain.forcing_terms = []
         domain.forcing_terms.append( Wind_stress(s, phi) )
-        
+
         domain.compute_forcing_terms()
 
 
         const = eta_w*rho_a/rho_w
-        
+
         #Convert to radians
         phi = phi*pi/180
-        
+
         #Compute velocity vector (u, v)
         u = s*cos(phi)
@@ -911,5 +911,5 @@
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)        
+        domain.set_quantity('friction', 0)
 
         Br = Reflective_boundary(domain)
@@ -924,10 +924,10 @@
         domain.forcing_terms = []
         domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
-        
+
         domain.compute_forcing_terms()
 
         #Compute reference solution
         const = eta_w*rho_a/rho_w
-        
+
         N = domain.number_of_elements
 
@@ -939,5 +939,5 @@
         s_vec = speed(t,x,y)
         phi_vec = angle(t,x,y)
-        
+
 
         for k in range(N):
@@ -945,5 +945,5 @@
             phi = phi_vec[k]*pi/180
             s = s_vec[k]
-        
+
             #Compute velocity vector (u, v)
             u = s*cos(phi)
@@ -966,5 +966,5 @@
         from util import file_function
         import time
-        
+
 
         a = [0.0, 0.0]
@@ -984,5 +984,5 @@
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)        
+        domain.set_quantity('friction', 0)
 
         Br = Reflective_boundary(domain)
@@ -995,5 +995,5 @@
         #Take x=1 and y=0
         filename = 'test_windstress_from_file.txt'
-        start = time.mktime(time.strptime('2000', '%Y'))        
+        start = time.mktime(time.strptime('2000', '%Y'))
         fid = open(filename, 'w')
         dt = 1  #One second interval
@@ -1006,18 +1006,18 @@
                                       angle(t,[1],[0])[0]))
             t += dt
-    
+
         fid.close()
 
         #Setup wind stress
         F = file_function(filename)
-        W = Wind_stress(F) 
+        W = Wind_stress(F)
         domain.forcing_terms = []
         domain.forcing_terms.append(W)
-        
+
         domain.compute_forcing_terms()
 
         #Compute reference solution
         const = eta_w*rho_a/rho_w
-        
+
         N = domain.number_of_elements
 
@@ -1026,9 +1026,9 @@
         s = speed(t,[1],[0])[0]
         phi = angle(t,[1],[0])[0]
-        
+
         #Convert to radians
         phi = phi*pi/180
 
-        
+
         #Compute velocity vector (u, v)
         u = s*cos(phi)
@@ -1049,5 +1049,5 @@
         are wrong - e.g. returns a scalar
         """
-        
+
         from config import rho_a, rho_w, eta_w
         from math import pi, cos, sin, sqrt
@@ -1069,5 +1069,5 @@
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)        
+        domain.set_quantity('friction', 0)
 
         Br = Reflective_boundary(domain)
@@ -1088,5 +1088,5 @@
             msg = 'Should have raised exception'
             raise msg
-        
+
 
         try:
@@ -1098,5 +1098,5 @@
             msg = 'Should have raised exception'
             raise msg
-        
+
         try:
             domain.forcing_terms.append(Wind_stress(s = speed,
@@ -1107,9 +1107,9 @@
             msg = 'Should have raised exception'
             raise msg
-        
+
 
     #####################################################
     def test_first_order_extrapolator_const_z(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1122,5 +1122,5 @@
         #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
@@ -1129,5 +1129,5 @@
         val3 = 2.+8.0/3
 
-        zl=zr=-3.75 #Assume constant bed (must be less than stage)        
+        zl=zr=-3.75 #Assume constant bed (must be less than stage)
         domain.set_quantity('elevation', zl*ones( (4,3) ))
         domain.set_quantity('stage', [[val0, val0-1, val0-2],
@@ -1143,7 +1143,7 @@
         #Check that centroid values were distributed to vertices
         C = domain.quantities['stage'].centroid_values
-        for i in range(3): 
+        for i in range(3):
             assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
-        
+
 
     def test_first_order_limiter_variable_z(self):
@@ -1151,5 +1151,5 @@
         from Numeric import alltrue, greater_equal
         from config import epsilon
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1162,5 +1162,5 @@
         #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
@@ -1177,13 +1177,13 @@
 
         E = domain.quantities['elevation'].vertex_values
-        L = domain.quantities['stage'].vertex_values        
-
-        
+        L = domain.quantities['stage'].vertex_values
+
+
         #Check that some stages are not above elevation (within eps)
         #- so that the limiter has something to work with
-        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))                
+        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
 
         domain.order = 1
-        domain.distribute_to_vertices_and_edges()                            
+        domain.distribute_to_vertices_and_edges()
 
         #Check that all stages are above elevation (within eps)
@@ -1191,9 +1191,9 @@
 
 
-    #####################################################    
+    #####################################################
     def test_distribute_basic(self):
-        #Using test data generated by pyvolution-2 
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1206,7 +1206,7 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
-	
+
         val0 = 2.
         val1 = 4.
@@ -1216,21 +1216,21 @@
         domain.set_quantity('stage', [val0, val1, val2, val3], 'centroids')
         L = domain.quantities['stage'].vertex_values
-	
+
 	#First order
-        domain.order = 1	
+        domain.order = 1
         domain.distribute_to_vertices_and_edges()
 	assert allclose(L[1], val1)
-	
+
 	#Second order
-        domain.order = 2	
+        domain.order = 2
         domain.distribute_to_vertices_and_edges()
 	assert allclose(L[1], [2.2, 4.9, 4.9])
-		
-
-	
-    def test_distribute_away_from_bed(self):		
-        #Using test data generated by pyvolution-2 
+
+
+
+    def test_distribute_away_from_bed(self):
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1243,31 +1243,31 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
-        L = domain.quantities['stage'].vertex_values	
-        
+        L = domain.quantities['stage'].vertex_values
+
 	def stage(x,y):
             return x**2
-        
+
 	domain.set_quantity('stage', stage, 'centroids')
-	
-        a, b = domain.quantities['stage'].compute_gradients()	    	
+
+        a, b = domain.quantities['stage'].compute_gradients()
 	assert allclose(a[1], 3.33333334)
 	assert allclose(b[1], 0.0)
-	
+
         domain.order = 1
         domain.distribute_to_vertices_and_edges()
         assert allclose(L[1], 1.77777778)
-	
+
         domain.order = 2
         domain.distribute_to_vertices_and_edges()
-        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])	
-
-
-	
-    def test_distribute_away_from_bed1(self):		
-        #Using test data generated by pyvolution-2 
+        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
+
+
+
+    def test_distribute_away_from_bed1(self):
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1280,30 +1280,30 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
-        L = domain.quantities['stage'].vertex_values	
-        
+        L = domain.quantities['stage'].vertex_values
+
 	def stage(x,y):
             return x**4+y**2
-        
+
 	domain.set_quantity('stage', stage, 'centroids')
-        #print domain.quantities['stage'].centroid_values		
-	
-        a, b = domain.quantities['stage'].compute_gradients()	    	
+        #print domain.quantities['stage'].centroid_values
+
+        a, b = domain.quantities['stage'].compute_gradients()
 	assert allclose(a[1], 25.18518519)
 	assert allclose(b[1], 3.33333333)
-	
+
         domain.order = 1
         domain.distribute_to_vertices_and_edges()
         assert allclose(L[1], 4.9382716)
-	
+
         domain.order = 2
         domain.distribute_to_vertices_and_edges()
-        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])	
-
-    
-    
+        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
+
+
+
     def test_distribute_near_bed(self):
-        #Using test data generated by pyvolution-2 
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
 
@@ -1318,9 +1318,9 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
 
 
-        #Set up for a gradient of (3,0) at mid triangle         
+        #Set up for a gradient of (3,0) at mid triangle
         def slope(x, y):
             return 10*x
@@ -1333,22 +1333,22 @@
         domain.set_quantity('stage', stage, 'centroids')
 
-        #print domain.quantities['elevation'].centroid_values	
+        #print domain.quantities['elevation'].centroid_values
         #print domain.quantities['stage'].centroid_values
-	
+
         E = domain.quantities['elevation'].vertex_values
         L = domain.quantities['stage'].vertex_values
 
         #print E
-        domain.order = 1	
+        domain.order = 1
         domain.distribute_to_vertices_and_edges()
 	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
 	assert allclose(L[1], [0.1, 20.1, 20.1])
-        
-        domain.order = 2	
+
+        domain.order = 2
         domain.distribute_to_vertices_and_edges()
 	assert allclose(L[1], [0.1, 20.1, 20.1])
 
     def test_distribute_near_bed1(self):
-        #Using test data generated by pyvolution-2 
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
 
@@ -1363,9 +1363,9 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
+
         domain = Domain(points, vertices)
 
 
-        #Set up for a gradient of (3,0) at mid triangle         
+        #Set up for a gradient of (3,0) at mid triangle
         def slope(x, y):
             return x**4+y**2
@@ -1378,7 +1378,7 @@
         domain.set_quantity('stage', stage)
 
-        #print domain.quantities['elevation'].centroid_values	
+        #print domain.quantities['elevation'].centroid_values
         #print domain.quantities['stage'].centroid_values
-	
+
         E = domain.quantities['elevation'].vertex_values
         L = domain.quantities['stage'].vertex_values
@@ -1388,14 +1388,14 @@
         domain.distribute_to_vertices_and_edges()
 	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
-	assert allclose(L[1], [4.1, 16.1, 20.1])        
-
-        domain.order = 2	
+	assert allclose(L[1], [4.1, 16.1, 20.1])
+
+        domain.order = 2
         domain.distribute_to_vertices_and_edges()
 	assert allclose(L[1], [4.1, 16.1, 20.1])
 
- 
+
 
     def test_second_order_distribute_real_data(self):
-        #Using test data generated by pyvolution-2 
+        #Using test data generated by pyvolution-2
         #Assuming no friction and flat bed (0.0)
 
@@ -1411,5 +1411,5 @@
         #bae, efb, cbf, feg
         vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
-        
+
         domain = Domain(points, vertices)
 
@@ -1418,22 +1418,22 @@
 
         domain.set_quantity('elevation', slope)
-        domain.set_quantity('stage', 
-			    [0.01298164, 0.00365611, 0.01440365, -0.0381856437096], 
+        domain.set_quantity('stage',
+			    [0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
 			    'centroids')
-        domain.set_quantity('xmomentum', 
+        domain.set_quantity('xmomentum',
 			    [0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
 			    'centroids')
-        domain.set_quantity('ymomentum', 			    
+        domain.set_quantity('ymomentum',
 			    [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
-			    'centroids')			    
-	
+			    'centroids')
+
         E = domain.quantities['elevation'].vertex_values
         L = domain.quantities['stage'].vertex_values
-        X = domain.quantities['xmomentum'].vertex_values	
-        Y = domain.quantities['ymomentum'].vertex_values		
+        X = domain.quantities['xmomentum'].vertex_values
+        Y = domain.quantities['ymomentum'].vertex_values
 
         #print E
         domain.order = 2
-        domain.beta_h = 0.0 #Use first order in h-limiter        
+        domain.beta_h = 0.0 #Use first order in h-limiter
         domain.distribute_to_vertices_and_edges()
 
@@ -1446,11 +1446,11 @@
 	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
 
-	
-	
+
+
     def test_balance_deep_and_shallow(self):
         """Test that balanced limiters preserve conserved quantites.
 	"""
         import copy
-	
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1465,5 +1465,5 @@
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
 
-        mesh = Domain(points, elements)        
+        mesh = Domain(points, elements)
         mesh.check_integrity()
 
@@ -1471,8 +1471,8 @@
         mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
         mesh.set_quantity('elevation', 0) #Flat bed
-	stage = mesh.quantities['stage']     	
-	
+	stage = mesh.quantities['stage']
+
 	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
-    
+
         #Limit
 	mesh.distribute_to_vertices_and_edges()
@@ -1483,26 +1483,26 @@
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
-        
+
 
         #Now try with a non-flat bed - closely hugging initial stage in places
 	#This will create alphas in the range [0, 0.478260, 1]
         mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-        mesh.set_quantity('elevation', [[0,0,0], 
-					[1.8,1.9,5.9], 
-					[4.6,0,0], 
+        mesh.set_quantity('elevation', [[0,0,0],
+					[1.8,1.9,5.9],
+					[4.6,0,0],
 					[0,2,4]])
-	stage = mesh.quantities['stage']     	
-	
+	stage = mesh.quantities['stage']
+
 	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
-	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy	
-    
+	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
+
         #Limit
 	mesh.distribute_to_vertices_and_edges()
 
-	
+
         #Assert that all vertex quantities have changed
         for k in range(mesh.number_of_elements):
-	    #print ref_vertex_values[k,:], stage.vertex_values[k,:]	
-            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])		
+	    #print ref_vertex_values[k,:], stage.vertex_values[k,:]
+            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
         #and assert that quantities are still conserved
         from Numeric import sum
@@ -1510,6 +1510,6 @@
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
-        
-    
+
+
         #Also check that Python and C version produce the same
         assert allclose (stage.vertex_values,
@@ -1518,12 +1518,12 @@
 			  [6.93333333, 4.53333333, 4.53333333],
 			  [5.33333333, 5.33333333, 5.33333333]])
-			  
-			 
-
-        
+
+
+
+
     def test_conservation_1(self):
         """Test that stage is conserved globally
-	
-	This one uses a flat bed, reflective bdries and a suitable 
+
+	This one uses a flat bed, reflective bdries and a suitable
 	initial condition
 	"""
@@ -1546,7 +1546,7 @@
 
         domain.set_quantity('elevation', 0)
-        domain.set_quantity('friction', 0)	
-        domain.set_quantity('stage', x_slope)	
-	
+        domain.set_quantity('friction', 0)
+        domain.set_quantity('stage', x_slope)
+
         # Boundary conditions (reflective everywhere)
         Br = Reflective_boundary(domain)
@@ -1555,28 +1555,28 @@
         domain.check_integrity()
 
-	#domain.visualise = True #If you want to take a sticky beak  
-	
-	initial_volume = domain.quantities['stage'].get_integral()         
-	initial_xmom = domain.quantities['xmomentum'].get_integral()      
-	
+	#domain.visualise = True #If you want to take a sticky beak
+
+	initial_volume = domain.quantities['stage'].get_integral()
+	initial_xmom = domain.quantities['xmomentum'].get_integral()
+
 	#print initial_xmom
 
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
-            volume =  domain.quantities['stage'].get_integral()         
-	    assert allclose (volume, initial_volume) 
-	 
+            volume =  domain.quantities['stage'].get_integral()
+	    assert allclose (volume, initial_volume)
+
             #I don't believe that the total momentum should be the same
 	    #It starts with zero and ends with zero though
-            #xmom = domain.quantities['xmomentum'].get_integral()   
+            #xmom = domain.quantities['xmomentum'].get_integral()
 	    #print xmom
-	    #assert allclose (xmom, initial_xmom) 
-        
-
-	
+	    #assert allclose (xmom, initial_xmom)
+
+
+
     def test_conservation_2(self):
         """Test that stage is conserved globally
-	
-	This one uses a slopy bed, reflective bdries and a suitable 
+
+	This one uses a slopy bed, reflective bdries and a suitable
 	initial condition
 	"""
@@ -1599,7 +1599,7 @@
 
         domain.set_quantity('elevation', x_slope)
-        domain.set_quantity('friction', 0)	
-        domain.set_quantity('stage', 0.4) #Steady	
-	
+        domain.set_quantity('friction', 0)
+        domain.set_quantity('stage', 0.4) #Steady
+
         # Boundary conditions (reflective everywhere)
         Br = Reflective_boundary(domain)
@@ -1608,26 +1608,26 @@
         domain.check_integrity()
 
-	#domain.visualise = True #If you want to take a sticky beak  
-	
-	initial_volume = domain.quantities['stage'].get_integral()         
-	initial_xmom = domain.quantities['xmomentum'].get_integral()      
-	
+	#domain.visualise = True #If you want to take a sticky beak
+
+	initial_volume = domain.quantities['stage'].get_integral()
+	initial_xmom = domain.quantities['xmomentum'].get_integral()
+
 	#print initial_xmom
 
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
-            volume =  domain.quantities['stage'].get_integral()         
-	    assert allclose (volume, initial_volume) 
-	 
+            volume =  domain.quantities['stage'].get_integral()
+	    assert allclose (volume, initial_volume)
+
             #FIXME: What would we expect from momentum
-            #xmom = domain.quantities['xmomentum'].get_integral()   
+            #xmom = domain.quantities['xmomentum'].get_integral()
 	    #print xmom
-	    #assert allclose (xmom, initial_xmom) 
-
-	    
+	    #assert allclose (xmom, initial_xmom)
+
+
     def test_conservation_3(self):
         """Test that stage is conserved globally
-	
-	This one uses a larger grid, convoluted bed, reflective bdries and a suitable 
+
+	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
 	initial condition
 	"""
@@ -1650,20 +1650,20 @@
 	    z = 0*x
 	    for i in range(len(x)):
-                if x[i] < 0.3:  
+                if x[i] < 0.3:
                     z[i] = x[i]/3
-                if 0.3 <= x[i] < 0.5:  
+                if 0.3 <= x[i] < 0.5:
                     z[i] = -0.5
-                if 0.5 <= x[i] < 0.7:  
-                    z[i] = 0.39     		    
-                if 0.7 <= x[i]:  
-                    z[i] = x[i]/3		    
-            return z		    		    
-		    		    
-		    
+                if 0.5 <= x[i] < 0.7:
+                    z[i] = 0.39
+                if 0.7 <= x[i]:
+                    z[i] = x[i]/3
+            return z
+
+
 
         domain.set_quantity('elevation', x_slope)
-        domain.set_quantity('friction', 0)	
-        domain.set_quantity('stage', 0.4) #Steady	
-	
+        domain.set_quantity('friction', 0)
+        domain.set_quantity('stage', 0.4) #Steady
+
         # Boundary conditions (reflective everywhere)
         Br = Reflective_boundary(domain)
@@ -1672,37 +1672,37 @@
         domain.check_integrity()
 
-	#domain.visualise = True #If you want to take a sticky beak  
-	
-	initial_volume = domain.quantities['stage'].get_integral()         
-	initial_xmom = domain.quantities['xmomentum'].get_integral()      
-	
+	#domain.visualise = True #If you want to take a sticky beak
+
+	initial_volume = domain.quantities['stage'].get_integral()
+	initial_xmom = domain.quantities['xmomentum'].get_integral()
+
 	import copy
 	ref_centroid_values =\
 		 copy.copy(domain.quantities['stage'].centroid_values)
-		 		
-	#print 'ORG', domain.quantities['stage'].centroid_values				
+
+	#print 'ORG', domain.quantities['stage'].centroid_values
 	domain.distribute_to_vertices_and_edges()
-	
-	
+
+
 	#print domain.quantities['stage'].centroid_values
 	assert allclose(domain.quantities['stage'].centroid_values,
-	                ref_centroid_values)		
-	
-	
+	                ref_centroid_values)
+
+
 	#Check that initial limiter doesn't violate cons quan
 	assert allclose (domain.quantities['stage'].get_integral(),
-	                 initial_volume) 
-	
+	                 initial_volume)
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
-            volume =  domain.quantities['stage'].get_integral()         
+            volume =  domain.quantities['stage'].get_integral()
             #print t, volume, initial_volume
-	    assert allclose (volume, initial_volume) 	    
-	            
-	
+	    assert allclose (volume, initial_volume)
+
+
     def test_conservation_4(self):
         """Test that stage is conserved globally
-	
-	This one uses a larger grid, convoluted bed, reflective bdries and a suitable 
+
+	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
 	initial condition
 	"""
@@ -1725,23 +1725,23 @@
 	    z = 0*x
 	    for i in range(len(x)):
-                if x[i] < 0.3:  
+                if x[i] < 0.3:
                     z[i] = x[i]/3
-                if 0.3 <= x[i] < 0.5:  
+                if 0.3 <= x[i] < 0.5:
                     z[i] = -0.5
                 if 0.5 <= x[i] < 0.7:
-                    #z[i] = 0.3 #OK with beta == 0.2		  
+                    #z[i] = 0.3 #OK with beta == 0.2
                     z[i] = 0.34 #OK with beta == 0.0
-                    #z[i] = 0.35#Fails after 80 timesteps with an error 
-		                #of the order 1.0e-5    
-                if 0.7 <= x[i]:  
-                    z[i] = x[i]/3		    
-            return z		    		    
-		    		    
-		    
+                    #z[i] = 0.35#Fails after 80 timesteps with an error
+		                #of the order 1.0e-5
+                if 0.7 <= x[i]:
+                    z[i] = x[i]/3
+            return z
+
+
 
         domain.set_quantity('elevation', x_slope)
-        domain.set_quantity('friction', 0)	
-        domain.set_quantity('stage', 0.4) #Steady	
-	
+        domain.set_quantity('friction', 0)
+        domain.set_quantity('stage', 0.4) #Steady
+
         # Boundary conditions (reflective everywhere)
         Br = Reflective_boundary(domain)
@@ -1750,47 +1750,47 @@
         domain.check_integrity()
 
-	#domain.visualise = True #If you want to take a sticky beak  
-	
-	initial_volume = domain.quantities['stage'].get_integral()         
-	initial_xmom = domain.quantities['xmomentum'].get_integral()      
-	
+	#domain.visualise = True #If you want to take a sticky beak
+
+	initial_volume = domain.quantities['stage'].get_integral()
+	initial_xmom = domain.quantities['xmomentum'].get_integral()
+
 	import copy
 	ref_centroid_values =\
 		 copy.copy(domain.quantities['stage'].centroid_values)
-		 		
+
 	#Test limiter by itself
 	domain.distribute_to_vertices_and_edges()
-	
+
 	#Check that initial limiter doesn't violate cons quan
 	assert allclose (domain.quantities['stage'].get_integral(),
-	                 initial_volume) 
-        #NOTE: This would fail if any initial stage was less than the 
+	                 initial_volume)
+        #NOTE: This would fail if any initial stage was less than the
 	#corresponding bed elevation - but that is reasonable.
-			 
-		
+
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
-            volume =  domain.quantities['stage'].get_integral()         
+            volume =  domain.quantities['stage'].get_integral()
 
             #print t, volume, initial_volume
-    
-	    
-	    #if not allclose (volume, initial_volume): 
+
+
+	    #if not allclose (volume, initial_volume):
 	    #    print 't==4.05'
             #    for k in range(domain.number_of_elements):
 	    # 	    pass
 	    #	    print domain.quantities['stage'].centroid_values[k] -\
-            #             ref_centroid_values[k] 
-            
-	    assert allclose (volume, initial_volume) 
-        
-
-	
-	
-	
+            #             ref_centroid_values[k]
+
+	    assert allclose (volume, initial_volume)
+
+
+
+
+
     def test_conservation_5(self):
-        """Test that momentum is conserved globally in 
+        """Test that momentum is conserved globally in
 	steady state scenario
-	
+
 	This one uses a slopy bed, dirichlet and reflective bdries
 	"""
@@ -1813,42 +1813,42 @@
 
         domain.set_quantity('elevation', x_slope)
-        domain.set_quantity('friction', 0)	
-        domain.set_quantity('stage', 0.4) #Steady	
-	
+        domain.set_quantity('friction', 0)
+        domain.set_quantity('stage', 0.4) #Steady
+
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)	
+        Br = Reflective_boundary(domain)
         Bleft = Dirichlet_boundary([0.5,0,0])
-        Bright = Dirichlet_boundary([0.1,0,0])	
-        domain.set_boundary({'left': Bleft, 'right': Bright, 
+        Bright = Dirichlet_boundary([0.1,0,0])
+        domain.set_boundary({'left': Bleft, 'right': Bright,
 			     'top': Br, 'bottom': Br})
 
         domain.check_integrity()
 
-	#domain.visualise = True #If you want to take a sticky beak  
-	
-	initial_volume = domain.quantities['stage'].get_integral()         
-	initial_xmom = domain.quantities['xmomentum'].get_integral()      
-	
+	#domain.visualise = True #If you want to take a sticky beak
+
+	initial_volume = domain.quantities['stage'].get_integral()
+	initial_xmom = domain.quantities['xmomentum'].get_integral()
+
 
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
-            stage =  domain.quantities['stage'].get_integral()         
-            xmom = domain.quantities['xmomentum'].get_integral()   	    
-            ymom = domain.quantities['ymomentum'].get_integral()   	    	    
-	    
+            stage =  domain.quantities['stage'].get_integral()
+            xmom = domain.quantities['xmomentum'].get_integral()
+            ymom = domain.quantities['ymomentum'].get_integral()
+
 	    if allclose(t, 6):  #Steady state reached
                 steady_xmom = domain.quantities['xmomentum'].get_integral()
                 steady_ymom = domain.quantities['ymomentum'].get_integral()
-                steady_stage = domain.quantities['stage'].get_integral()		
-		
+                steady_stage = domain.quantities['stage'].get_integral()
+
 	    if t > 6:
                 #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
-                assert allclose(xmom, steady_xmom)		
+                assert allclose(xmom, steady_xmom)
                 assert allclose(ymom, steady_ymom)
-                assert allclose(stage, steady_stage)		
-						
-	    
-		
-	        
+                assert allclose(stage, steady_stage)
+
+
+
+
     def test_second_order_flat_bed_onestep(self):
 
@@ -1879,5 +1879,5 @@
         #Data from earlier version of pyvolution
         assert allclose(domain.min_timestep, 0.0396825396825)
-        assert allclose(domain.max_timestep, 0.0396825396825) 
+        assert allclose(domain.max_timestep, 0.0396825396825)
 
         assert allclose(domain.quantities['stage'].centroid_values[:12],
@@ -1897,5 +1897,5 @@
                          0.00024152, 0.02656103, 0.00024152, 0.02656103,
                          0.00024152, 0.02656103, 0.00040506, 0.0272623])
-        
+
         assert allclose(domain.quantities['stage'].vertex_values[:12,2],
                         [0.00182037, 0.02656103, 0.00676264,
@@ -1919,5 +1919,5 @@
 
 
-        
+
     def test_second_order_flat_bed_moresteps(self):
 
@@ -1948,5 +1948,5 @@
         #Data from earlier version of pyvolution
         #assert allclose(domain.min_timestep, 0.0396825396825)
-        #assert allclose(domain.max_timestep, 0.0396825396825) 
+        #assert allclose(domain.max_timestep, 0.0396825396825)
         #print domain.quantities['stage'].centroid_values
 
@@ -1972,9 +1972,9 @@
         Br = Reflective_boundary(domain)
         Bd = Dirichlet_boundary([0.2,0.,0.])
-        
+
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
         domain.check_integrity()
 
-  
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
@@ -1982,5 +1982,5 @@
             #domain.write_time()
 
-        #FIXME: These numbers were from version before 25/10        
+        #FIXME: These numbers were from version before 25/10
         #assert allclose(domain.min_timestep, 0.0140413643926)
         #assert allclose(domain.max_timestep, 0.0140947355753)
@@ -1994,6 +1994,6 @@
 
             #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
-            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
-        
+            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
+
     def test_flatbed_second_order(self):
         from mesh_factory import rectangular
@@ -2020,5 +2020,5 @@
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
         domain.check_integrity()
-  
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
@@ -2027,5 +2027,5 @@
 
         assert allclose(domain.min_timestep, 0.0210448446782)
-        assert allclose(domain.max_timestep, 0.0210448446782) 
+        assert allclose(domain.max_timestep, 0.0210448446782)
 
         #print domain.quantities['stage'].vertex_values[:4,0]
@@ -2036,5 +2036,5 @@
         #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
         #                [0.00101913,0.05352143,0.00104852,0.05354394])
-        
+
         assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
                         [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
@@ -2042,10 +2042,10 @@
         #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
         #                [0.00090581,0.03685719,0.00088303,0.03687313])
-        
+
         assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
                         [-0.00139463,0.0006156,-0.00060364,0.00061827])
 
 
-        
+
     def test_flatbed_second_order_distribute(self):
         #Use real data from pyvolution 2
@@ -2078,5 +2078,5 @@
         for V in [False, True]:
             if V:
-                #Set centroids as if system had been evolved            
+                #Set centroids as if system had been evolved
                 L = zeros(2*N*N, Float)
                 L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
@@ -2091,5 +2091,5 @@
                           4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
                           0.00000000e+000, 5.57305948e-005]
-        
+
                 X = zeros(2*N*N, Float)
                 X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
@@ -2118,11 +2118,11 @@
                         0.00000000e+000, -2.57635178e-005]
 
-        
+
                 domain.set_quantity('stage', L, 'centroids')
                 domain.set_quantity('xmomentum', X, 'centroids')
                 domain.set_quantity('ymomentum', Y, 'centroids')
-        
+
                 domain.check_integrity()
-            else:    
+            else:
                 #Evolution
                 for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
@@ -2130,5 +2130,5 @@
                 assert allclose(domain.min_timestep, 0.0210448446782)
                 assert allclose(domain.max_timestep, 0.0210448446782)
-        
+
 
             #Centroids were correct but not vertices.
@@ -2149,13 +2149,13 @@
             ##print domain.quantities['xmomentum'].centroid_values[17], V
             ##print
-            if not V:            
+            if not V:
                 assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
-            else:  
-                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)            
+            else:
+                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
 
             import copy
             XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
             assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
-        
+
             domain.distribute_to_vertices_and_edges()
 
@@ -2187,10 +2187,10 @@
             #print domain.quantities['xmomentum'].vertex_values[:4,0]
             #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
-            #                [0.00090581,0.03685719,0.00088303,0.03687313])        
-
-
-
-        
-        
+            #                [0.00090581,0.03685719,0.00088303,0.03687313])
+
+
+
+
+
     def test_bedslope_problem_first_order(self):
 
@@ -2228,5 +2228,5 @@
         assert allclose(domain.max_timestep, 0.050010003001)
 
-        
+
     def test_bedslope_problem_first_order_moresteps(self):
 
@@ -2264,5 +2264,5 @@
         #Data from earlier version of pyvolution
         #print domain.quantities['stage'].centroid_values
-        
+
         assert allclose(domain.quantities['stage'].centroid_values,
                         [-0.02998628, -0.01520652, -0.03043492,
@@ -2290,6 +2290,6 @@
                         -0.14175896, -0.14560798, -0.13911658,
                         -0.14439383, -0.13924047, -0.14829043])
-                        
-        
+
+
     def test_bedslope_problem_second_order_one_step(self):
 
@@ -2349,13 +2349,13 @@
         #print domain.quantities['stage'].extrapolate_second_order()
         #domain.distribute_to_vertices_and_edges()
-        #print domain.quantities['stage'].vertex_values[:,0]        
-        
+        #print domain.quantities['stage'].vertex_values[:,0]
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
-            #domain.write_time()            
+            #domain.write_time()
             pass
 
 
-        #print domain.quantities['stage'].centroid_values 
+        #print domain.quantities['stage'].centroid_values
 	assert allclose(domain.quantities['stage'].centroid_values,
 		 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
@@ -2377,5 +2377,5 @@
 		  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
 
-	
+
 
     def test_bedslope_problem_second_order_two_steps(self):
@@ -2392,5 +2392,5 @@
         domain.smooth = False
         domain.default_order=2
-        domain.beta_h = 0.0 #Use first order in h-limiter        
+        domain.beta_h = 0.0 #Use first order in h-limiter
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2437,13 +2437,13 @@
         #print domain.quantities['stage'].extrapolate_second_order()
         #domain.distribute_to_vertices_and_edges()
-        #print domain.quantities['stage'].vertex_values[:,0]        
-        
+        #print domain.quantities['stage'].vertex_values[:,0]
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
             pass
 
-        
+
         #Data from earlier version of pyvolution ft=0.1
-        assert allclose(domain.min_timestep, 0.0376895634803) 
+        assert allclose(domain.min_timestep, 0.0376895634803)
         assert allclose(domain.max_timestep, 0.0415635655309)
 
@@ -2469,7 +2469,7 @@
                          -0.25031463, -0.24454764, -0.24885323, -0.24286438])
 
-	
-	
-	        
+
+
+
     def test_bedslope_problem_second_order_two_yieldsteps(self):
 
@@ -2485,5 +2485,5 @@
         domain.smooth = False
         domain.default_order=2
-        domain.beta_h = 0.0 #Use first order in h-limiter        
+        domain.beta_h = 0.0 #Use first order in h-limiter
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2530,14 +2530,14 @@
         #print domain.quantities['stage'].extrapolate_second_order()
         #domain.distribute_to_vertices_and_edges()
-        #print domain.quantities['stage'].vertex_values[:,0]        
-        
+        #print domain.quantities['stage'].vertex_values[:,0]
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
-            #domain.write_time()            
+            #domain.write_time()
             pass
 
 
-	
-	assert allclose(domain.quantities['stage'].centroid_values, 
+
+	assert allclose(domain.quantities['stage'].centroid_values,
 		 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
 		  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
@@ -2574,5 +2574,5 @@
         domain.smooth = False
         domain.default_order=2
-        domain.beta_h = 0.0 #Use first order in h-limiter        
+        domain.beta_h = 0.0 #Use first order in h-limiter
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2619,11 +2619,11 @@
         #print domain.quantities['stage'].extrapolate_second_order()
         #domain.distribute_to_vertices_and_edges()
-        #print domain.quantities['stage'].vertex_values[:,0]        
-        
+        #print domain.quantities['stage'].vertex_values[:,0]
+
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
             pass
 
-        
+
         assert allclose(domain.quantities['stage'].centroid_values,
      [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
@@ -2631,5 +2631,5 @@
       -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
       -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
-      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174, 
+      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
       -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
       -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
@@ -2644,16 +2644,16 @@
        0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
        0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
-       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113, 
+       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
        0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
-       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039, 
+       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
        0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
        0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
-       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247, 
+       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
        0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
        0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
        0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
 
-	
-        assert allclose(domain.quantities['ymomentum'].centroid_values,	
+
+        assert allclose(domain.quantities['ymomentum'].centroid_values,
      [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
       -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
@@ -2674,7 +2674,7 @@
       -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
       -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
-	
-
-	
+
+
+
 
     def test_temp_play(self):
@@ -2691,5 +2691,5 @@
         domain.smooth = False
         domain.default_order=2
-        domain.beta_h = 0.0 #Use first order in h-limiter        
+        domain.beta_h = 0.0 #Use first order in h-limiter
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2713,14 +2713,14 @@
 	assert allclose(domain.quantities['stage'].centroid_values[:4],
 			[0.00206836, 0.01296714, 0.00363415, 0.01438924])
-	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],			
+	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
 			[0.01360154, 0.00671133, 0.01264578, 0.00648503])
 	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
-			[-1.19201077e-003, -7.23647546e-004, 
+			[-1.19201077e-003, -7.23647546e-004,
 			-6.39083123e-005, 6.29815168e-005])
-	
-	
+
+
     def test_complex_bed(self):
         #No friction is tested here
-        
+
         from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
              Transmissive_boundary, Time_boundary,\
@@ -2729,5 +2729,5 @@
         from mesh_factory import rectangular
         from Numeric import array
-    
+
         N = 12
         points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
@@ -2740,5 +2740,5 @@
         domain.default_order=2
 
-        
+
         inflow_stage = 0.1
         Z = Weir(inflow_stage)
@@ -2757,5 +2757,5 @@
         #print domain.quantities['stage'].centroid_values
 
-        #FIXME: These numbers were from version before 25/10        
+        #FIXME: These numbers were from version before 25/10
         #assert allclose(domain.quantities['stage'].centroid_values,
 # [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
@@ -2810,8 +2810,8 @@
         """
         import time
-       
-	#Create sww file of simple propagation from left to right 
-	#through rectangular domain
-	
+
+        #Create sww file of simple propagation from left to right
+        #through rectangular domain
+
         from mesh_factory import rectangular
 
@@ -2821,22 +2821,22 @@
         #Create shallow water domain
         domain1 = Domain(points, vertices, boundary)
-	
+
         from util import mean
-        domain1.reduction = mean 
+        domain1.reduction = mean
         domain1.smooth = False  #Exact result
 
         domain1.default_order = 2
-	domain1.store = True	
+        domain1.store = True
         domain1.set_datadir('.')
         domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
-	
-	#FIXME: This is extremely important!
+
+        #FIXME: This is extremely important!
         #How can we test if they weren't stored?
         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
-		
+
         #Bed-slope and friction at vertices (and interpolated elsewhere)
         domain1.set_quantity('elevation', 0)
-        domain1.set_quantity('friction', 0)	
+        domain1.set_quantity('friction', 0)
 
         # Boundary conditions
@@ -2855,41 +2855,40 @@
 
         cv1 = domain1.quantities['stage'].centroid_values
-	
-	
-	#Create an triangle shaped domain (reusing coordinates from domain 1), 
-	#formed from the lower and right hand  boundaries and 
-	#the sw-ne diagonal 
-	#from domain 1. Call it domain2
-	
-	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
-	           [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
-		   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
-	
-	vertices = [ [1,2,0], 
-	             [3,4,1], [2,1,4], [4,5,2], 
-	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]	
+
+
+        #Create an triangle shaped domain (reusing coordinates from domain 1),
+        #formed from the lower and right hand  boundaries and
+        #the sw-ne diagonal
+        #from domain 1. Call it domain2
+
+        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
+                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
+
+        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
+                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
 
         boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
 	             (4,2):'right', (6,2):'right', (8,2):'right',
-		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
-		     
+		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
+
         domain2 = Domain(points, vertices, boundary)
-	
-	domain2.reduction = domain1.reduction
+
+        domain2.reduction = domain1.reduction
         domain2.smooth = False
         domain2.default_order = 2
-	
+
         #Bed-slope and friction at vertices (and interpolated elsewhere)
         domain2.set_quantity('elevation', 0)
         domain2.set_quantity('friction', 0)
-        domain2.set_quantity('stage', 0)        
+        domain2.set_quantity('stage', 0)
 
         # Boundary conditions
-        Br = Reflective_boundary(domain2)        
-	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
+        Br = Reflective_boundary(domain2)
+        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
 				      domain2)
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
-        domain2.check_integrity()        
-	
+        domain2.check_integrity()
+
 
 
@@ -2898,39 +2897,39 @@
             pass
 
-	
-	#Use output from domain1 as spatio-temporal boundary for domain2
-        #and verify that results at right hand side are close.	
+
+        #Use output from domain1 as spatio-temporal boundary for domain2
+        #and verify that results at right hand side are close.
 
         cv2 = domain2.quantities['stage'].centroid_values
 
         #print take(cv1, (12,14,16))  #Right
-	#print take(cv2, (4,6,8))
-        #print take(cv1, (0,6,12))  #Bottom	
-	#print take(cv2, (0,1,4))
-        #print take(cv1, (0,8,16))  #Diag	
-	#print take(cv2, (0,3,8))
-
-	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
-	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
+        #print take(cv2, (4,6,8))
+        #print take(cv1, (0,6,12))  #Bottom
+        #print take(cv2, (0,1,4))
+        #print take(cv1, (0,8,16))  #Diag
+        #print take(cv2, (0,3,8))
+
+        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
+        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
         assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
 
         #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)        
-	
-	
+        os.remove(domain1.filename + '.' + domain1.format)
+
+
 
     def test_spatio_temporal_boundary_2(self):
         """Test that boundary values can be read from file and interpolated
         in both time and space.
-	This is a more basic test, verifying that boundary object 
-	produces the expected results
+        This is a more basic test, verifying that boundary object
+        produces the expected results
 
 
         """
         import time
-       
-	#Create sww file of simple propagation from left to right 
-	#through rectangular domain
-	
+
+        #Create sww file of simple propagation from left to right
+        #through rectangular domain
+
         from mesh_factory import rectangular
 
@@ -2940,26 +2939,26 @@
         #Create shallow water domain
         domain1 = Domain(points, vertices, boundary)
-	
+
         from util import mean
-        domain1.reduction = mean 
+        domain1.reduction = mean
         domain1.smooth = True #To mimic MOST output
 
         domain1.default_order = 2
-	domain1.store = True	
+        domain1.store = True
         domain1.set_datadir('.')
         domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
-	
-	#FIXME: This is extremely important!
+
+        #FIXME: This is extremely important!
         #How can we test if they weren't stored?
         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
-		
+
         #Bed-slope and friction at vertices (and interpolated elsewhere)
         domain1.set_quantity('elevation', 0)
-        domain1.set_quantity('friction', 0)	
+        domain1.set_quantity('friction', 0)
 
         # Boundary conditions
         Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3,0,0])	
+        Bd = Dirichlet_boundary([0.3,0,0])
         domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
         #Initial condition
@@ -2973,129 +2972,130 @@
             #domain1.write_time()
 
-	
-	#Create an triangle shaped domain (coinciding with some
-        #coordinates from domain 1), 
-	#formed from the lower and right hand  boundaries and 
-	#the sw-ne diagonal 
-	#from domain 1. Call it domain2
-	
-	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+
+        #Create an triangle shaped domain (coinciding with some
+        #coordinates from domain 1),
+        #formed from the lower and right hand  boundaries and
+        #the sw-ne diagonal
+        #from domain 1. Call it domain2
+
+        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
 	           [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
 		   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
-	
-	vertices = [ [1,2,0], 
-	             [3,4,1], [2,1,4], [4,5,2], 
-	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]	
+
+        vertices = [ [1,2,0],
+	             [3,4,1], [2,1,4], [4,5,2],
+	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
 
         boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
 	             (4,2):'right', (6,2):'right', (8,2):'right',
-		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
-		     
+		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
+
         domain2 = Domain(points, vertices, boundary)
-	
-	domain2.reduction = domain1.reduction
+
+        domain2.reduction = domain1.reduction
         domain2.smooth = False
         domain2.default_order = 2
-	
+
         #Bed-slope and friction at vertices (and interpolated elsewhere)
         domain2.set_quantity('elevation', 0)
         domain2.set_quantity('friction', 0)
-        domain2.set_quantity('stage', 0)        
-
-	
-	#Read results for specific timesteps t=1 and t=2
-        from Scientific.IO.NetCDF import NetCDFFile	
-	fid = NetCDFFile(domain1.filename + '.' + domain1.format)
-	
-	x = fid.variables['x'][:]
-	y = fid.variables['y'][:]	
-	s1 = fid.variables['stage'][1,:]
-	s2 = fid.variables['stage'][2,:]	
-	fid.close()
-
-	from Numeric import take, reshape, concatenate	
-	shp = (len(x), 1)
-	points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)	
-	#The diagonal points of domain 1 are 0, 5, 10, 15
-        
-	#print points[0], points[5], points[10], points[15]
-	assert allclose( take(points, [0,5,10,15]), 
-	                 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])  
-
-	
+        domain2.set_quantity('stage', 0)
+
+
+        #Read results for specific timesteps t=1 and t=2
+        from Scientific.IO.NetCDF import NetCDFFile
+        fid = NetCDFFile(domain1.filename + '.' + domain1.format)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        s1 = fid.variables['stage'][1,:]
+        s2 = fid.variables['stage'][2,:]
+        fid.close()
+
+        from Numeric import take, reshape, concatenate
+        shp = (len(x), 1)
+        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
+        #The diagonal points of domain 1 are 0, 5, 10, 15
+
+        #print points[0], points[5], points[10], points[15]
+        assert allclose( take(points, [0,5,10,15]),
+	                 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
+
+
         # Boundary conditions
-        Br = Reflective_boundary(domain2)        
-	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
+        Br = Reflective_boundary(domain2)
+        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
 				      domain2)
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
-        domain2.check_integrity()        
-
-	#Test that interpolation points are the mid points of the all boundary
-	#segments	
-
-	boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
+        domain2.check_integrity()
+
+        #Test that interpolation points are the mid points of the all boundary
+        #segments
+
+        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
 			      [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
 			      [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
-			      
+
         boundary_midpoints.sort()
-	R = Bf.F.interpolation_points.tolist()	
-	R.sort()
-	assert allclose(boundary_midpoints, R)  
-	
-	#Check spatially interpolated output at time == 1
-	domain2.time = 1
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(0,0)
-	assert allclose(R0[0], (s1[0] + s1[5])/2)
-	
-	#Second diagonal midpoint
-	R0 = Bf.evaluate(3,0)
-	assert allclose(R0[0], (s1[5] + s1[10])/2)
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(8,0)
-	assert allclose(R0[0], (s1[10] + s1[15])/2)		
-		
-	#Check spatially interpolated output at time == 2
-	domain2.time = 2
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(0,0)
-	assert allclose(R0[0], (s2[0] + s2[5])/2)
-	
-	#Second diagonal midpoint
-	R0 = Bf.evaluate(3,0)
-	assert allclose(R0[0], (s2[5] + s2[10])/2)
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(8,0)
-	assert allclose(R0[0], (s2[10] + s2[15])/2)		
-	
-	
-	#Now check temporal interpolation
+        R = Bf.F.interpolation_points.tolist()
+        R.sort()
+        assert allclose(boundary_midpoints, R)
+
+        #Check spatially interpolated output at time == 1
+        domain2.time = 1
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], (s1[0] + s1[5])/2)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], (s1[5] + s1[10])/2)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], (s1[10] + s1[15])/2)
+
+        #Check spatially interpolated output at time == 2
+        domain2.time = 2
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], (s2[0] + s2[5])/2)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], (s2[5] + s2[10])/2)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], (s2[10] + s2[15])/2)
+
+
+        #Now check temporal interpolation
 
         domain2.time = 1 + 2.0/3
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(0,0)
-	assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
-	
-	#Second diagonal midpoint
-	R0 = Bf.evaluate(3,0)
-	assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
-	
-	#First diagonal midpoint
-	R0 = Bf.evaluate(8,0)
-	assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)		
-		
-	
-		
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
+
+
+
         #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)        	    
-#-------------------------------------------------------------
+        os.remove(domain1.filename + '.' + domain1.format)
+
+
+        #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Shallow_Water,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
-
