Changeset 1027

Timestamp:

Mar 7, 2005, 5:13:39 PM (20 years ago)

Author:

steve

Message:

Changed test files to define class with a useful name (shows up when using zeus editor)

Location:

inundation/ga/storm_surge/pyvolution

Files:

• test_pmesh2domain.py (modified) (2 diffs)
• test_shallow_water.py (modified) (151 diffs)

inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py

@@ -37 +37 @@
 
          fileName = tempfile.mktemp(".tsh")
+         print fileName
          file = open(fileName,"w")
          file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
@@ -163 +164 @@
 
 
-

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -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 +49 @@
 
     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 +64 @@
         #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 +91 @@
         assert r[2] == 2
 
-        w = rotate(r, normal, direction = -1)        
+        w = rotate(r, normal, direction = -1)
         assert allclose(w, q)
 
@@ -114 +114 @@
     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 +156 @@
     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 +166 @@
 
     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 +177 @@
 
     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 +185 @@
 
         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 +196 @@
 
         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 +207 @@
 
 
-        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 +227 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Third'}                                         
-                     
+                     (3, 2): 'Third'}
+
 
         domain = Domain(points, vertices, boundary)
@@ -241 +241 @@
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -248 +248 @@
         R = Reflective_boundary(domain)
         domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
-        
+
         domain.update_boundary()
 
@@ -254 +254 @@
         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 +275 @@
                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 +286 @@
 
     def test_boundary_conditionsII(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -303 +303 @@
                      (3, 1): 'Second',
                      (3, 2): 'Third',
-                     (0, 1): 'Internal'}                                         
-                     
+                     (0, 1): 'Internal'}
+
 
         domain = Domain(points, vertices, boundary)
@@ -317 +317 @@
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -333 +333 @@
         #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 +344 @@
         #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 +351 @@
 
         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 +396 @@
         e1 = domain.edgelengths[1, 1]
         e2 = domain.edgelengths[1, 2]
-        
+
         assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
 
@@ -410 +410 @@
     def test_compute_fluxes1(self):
         #Use values from previous version
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -421 +421 @@
         #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 +427 @@
         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 +435 @@
 
         #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 +444 @@
 
         #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 +452 @@
 
         #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 +459 @@
         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 +487 @@
                         [-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 +509 @@
         #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 +519 @@
 
         domain.set_quantity('elevation', zl*ones( (4,3) ))
-        
-        
+
+
         domain.set_quantity('stage', [[val0, val0-1, val0-2],
                                       [val1, val1+1, val1],
@@ -532 +532 @@
         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 +540 @@
 
         #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 +546 @@
 
         #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 +552 @@
 
         #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 +564 @@
         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 +585 @@
         #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 +594 @@
         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 +607 @@
         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 +615 @@
 
         #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 +621 @@
 
         #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 +627 @@
 
         #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 +638 @@
 
         total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
-        
+
         domain.compute_fluxes()
         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
@@ -647 +647 @@
 
     def test_catching_negative_heights(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -658 +658 @@
         #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 +681 @@
 
 
-    #####################################################    
+    #####################################################
     def test_initial_condition(self):
         """Test that initial condition is output at time == 0
@@ -688 +688 @@
         from config import g
         import copy
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -699 +699 @@
         #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 +725 @@
             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 +745 @@
         #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 +762 @@
             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 +783 @@
         #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 +797 @@
         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 +811 @@
         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 +817 @@
         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 +837 @@
         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 +854 @@
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        
+
         domain = Domain(points, vertices)
 
@@ -860 +860 @@
         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 +870 @@
         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 +911 @@
         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 +924 @@
         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 +939 @@
         s_vec = speed(t,x,y)
         phi_vec = angle(t,x,y)
-        
+
 
         for k in range(N):
@@ -945 +945 @@
             phi = phi_vec[k]*pi/180
             s = s_vec[k]
-        
+
             #Compute velocity vector (u, v)
             u = s*cos(phi)
@@ -966 +966 @@
         from util import file_function
         import time
-        
+
 
         a = [0.0, 0.0]
@@ -984 +984 @@
         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 +995 @@
         #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 +1006 @@
                                       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 +1026 @@
         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 +1049 @@
         are wrong - e.g. returns a scalar
         """
-        
+
         from config import rho_a, rho_w, eta_w
         from math import pi, cos, sin, sqrt
@@ -1069 +1069 @@
         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 +1088 @@
             msg = 'Should have raised exception'
             raise msg
-        
+
 
         try:
@@ -1098 +1098 @@
             msg = 'Should have raised exception'
             raise msg
-        
+
         try:
             domain.forcing_terms.append(Wind_stress(s = speed,
@@ -1107 +1107 @@
             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 +1122 @@
         #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 +1129 @@
         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 +1143 @@
         #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 +1151 @@
         from Numeric import alltrue, greater_equal
         from config import epsilon
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -1162 +1162 @@
         #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 +1177 @@
 
         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 +1191 @@
 
 
-    #####################################################    
+    #####################################################
     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 +1206 @@
         #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 +1216 @@
         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 +1243 @@
         #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 +1280 @@
         #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 +1318 @@
         #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 +1333 @@
         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 +1363 @@
         #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 +1378 @@
         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 +1388 @@
         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 +1411 @@
         #bae, efb, cbf, feg
         vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
-        
+
         domain = Domain(points, vertices)
 
@@ -1418 +1418 @@
 
         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 +1446 @@
         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 +1465 @@
         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 +1471 @@
         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 +1483 @@
             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 +1510 @@
             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 +1518 @@
                           [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 +1546 @@
 
         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 +1555 @@
         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 +1599 @@
 
         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 +1608 @@
         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 +1650 @@
             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 +1672 @@
         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 +1725 @@
             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 +1750 @@
         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 +1813 @@
 
         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 +1879 @@
         #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 +1897 @@
                          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 +1919 @@
 
 
-        
+
     def test_second_order_flat_bed_moresteps(self):
 
@@ -1948 +1948 @@
         #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 +1972 @@
         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 +1982 @@
             #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 +1994 @@
 
             #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 +2020 @@
         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 +2027 @@
 
         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 +2036 @@
         #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 +2042 @@
         #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 +2078 @@
         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 +2091 @@
                           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 +2118 @@
                         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 +2130 @@
                 assert allclose(domain.min_timestep, 0.0210448446782)
                 assert allclose(domain.max_timestep, 0.0210448446782)
-        
+
 
             #Centroids were correct but not vertices.
@@ -2149 +2149 @@
             ##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 +2187 @@
             #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 +2228 @@
         assert allclose(domain.max_timestep, 0.050010003001)
 
-        
+
     def test_bedslope_problem_first_order_moresteps(self):
 
@@ -2264 +2264 @@
         #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 +2290 @@
                         -0.14175896, -0.14560798, -0.13911658,
                         -0.14439383, -0.13924047, -0.14829043])
-                        
-        
+
+
     def test_bedslope_problem_second_order_one_step(self):
 
@@ -2349 +2349 @@
         #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 +2377 @@
                   -0.25865535, -0.24776246, -0.25865535, -0.24521113])
 
-        
+
 
     def test_bedslope_problem_second_order_two_steps(self):
@@ -2392 +2392 @@
         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 +2437 @@
         #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 +2469 @@
                          -0.25031463, -0.24454764, -0.24885323, -0.24286438])
 
-        
-        
-                
+
+
+
     def test_bedslope_problem_second_order_two_yieldsteps(self):
 
@@ -2485 +2485 @@
         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 +2530 @@
         #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 +2574 @@
         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 +2619 @@
         #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 +2631 @@
       -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 +2644 @@
        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 +2674 @@
       -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 +2691 @@
         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 +2713 @@
         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 +2729 @@
         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 +2740 @@
         domain.default_order=2
 
-        
+
         inflow_stage = 0.1
         Z = Weir(inflow_stage)
@@ -2757 +2757 @@
         #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 +2810 @@
         """
         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 +2821 @@
         #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 +2855 @@
 
         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 +2897 @@
             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 +2939 @@
         #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 +2972 @@
             #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)
-