Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1014)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1015)
@@ -526,5 +526,5 @@
         else:
             raise 'Unknown order'
-
+	
     #Take bed elevation into account when water heights are small
     balance_deep_and_shallow(domain)
@@ -567,4 +567,6 @@
     """
     
+    #FIXME: Look here for error
+    
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
@@ -574,4 +576,6 @@
     hc = wc - zc  #Water depths at centroids    
 
+    #print zc
+    #print '1', wc	
     #Update
     for k in range(domain.number_of_elements):
@@ -584,4 +588,5 @@
             xmomc[k] = ymomc[k] = 0.0
         
+    #print '2', wc		
 
 
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1014)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1015)
@@ -1624,10 +1624,10 @@
 	    #print xmom
 	    #assert allclose (xmom, initial_xmom) 
-        
-	
+
+	    
     def test_conservation_3(self):
         """Test that stage is conserved globally
 	
-	This one uses a convoluted bed, reflective bdries and a suitable 
+	This one uses a larger grid, convoluted bed, reflective bdries and a suitable 
 	initial condition
 	"""
@@ -1638,10 +1638,10 @@
 
         #Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = rectangular(2, 1)
 
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
         domain.smooth = False
-        domain.default_order=1  #FIXME: 2 will fail
+        domain.default_order=2
 	domain.beta_h = 0
 
@@ -1655,6 +1655,5 @@
                     z[i] = -0.5
                 if 0.5 <= x[i] < 0.7:  
-                    z[i] = 0.4       
-                    #z[i] = 0.7      		    
+                    z[i] = 0.39     		    
                 if 0.7 <= x[i]:  
                     z[i] = x[i]/3		    
@@ -1678,4 +1677,80 @@
 	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				
+	domain.distribute_to_vertices_and_edges()
+	
+	
+	#print domain.quantities['stage'].centroid_values
+	assert allclose(domain.quantities['stage'].centroid_values,
+	                ref_centroid_values)		
+	
+	
+	#Check that initial limiter doesn't violate cons quan
+	assert allclose (domain.quantities['stage'].get_integral(),
+	                 initial_volume) 
+	
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5):
+            volume =  domain.quantities['stage'].get_integral()         
+            #print t, 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 
+	initial condition
+	"""
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary,\
+             Dirichlet_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=2
+	domain.beta_h = 0.0
+
+	#IC
+        def x_slope(x, y):
+	    z = 0*x
+	    for i in range(len(x)):
+                if x[i] < 0.3:  
+                    z[i] = x[i]/3
+                if 0.3 <= x[i] < 0.5:  
+                    z[i] = -0.5
+                if 0.5 <= x[i] < 0.7:  
+                    z[i] = 0.34 #OK
+                    #z[i] = 0.35	#Fails	    
+                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	
+	
+        # Boundary conditions (reflective everywhere)
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, '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()      
+	
 	
 	#Let surface settle according to limiter 
@@ -1684,4 +1759,14 @@
 	#needs to respect the law of the limiters?
 	domain.distribute_to_vertices_and_edges()
+	
+	
+	#FIXME: At this point study the centroid of individual triangles.
+	
+	#Check that initial limiter doesn't violate cons quan
+	assert allclose (domain.quantities['stage'].get_integral(),
+	                 initial_volume) 
+	
+			 
+			 	
 	initial_volume = domain.quantities['stage'].get_integral()         
 		
@@ -2950,5 +3035,4 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
     runner = unittest.TextTestRunner()
     runner.run(suite)
