Index: /inundation/ga/storm_surge/pyvolution/general_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 1011)
@@ -267,4 +267,8 @@
 
 
-
-
+    def get_area(self):
+    	"""Return total area of mesh
+	"""
+
+	return sum(self.areas)		
+        
Index: /inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 1011)
@@ -414,7 +414,4 @@
         return polygon
         
-            
-
-        
 
     def check_integrity(self):
Index: /inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 1011)
@@ -568,4 +568,15 @@
             
 
+    def get_integral(self):
+        """Compute the integral of quantity across entire domain
+	"""
+	integral = 0
+	for k in range(self.domain.number_of_elements):
+            area = self.domain.areas[k]
+            qc = self.centroid_values[k]     
+            integral += qc*area
+	    
+	return integral    
+        
 
 class Conserved_quantity(Quantity):
Index: /inundation/ga/storm_surge/pyvolution/test_general_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_general_mesh.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/test_general_mesh.py	(revision 1011)
@@ -34,4 +34,15 @@
         assert domain.get_vertices([0,4]) == [domain.triangles[0],
                                               domain.triangles[4]]
+    def test_areas(self):
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        domain = Domain(points, vertices, boundary)
+
+	assert domain.get_area() == 1.0	        
+
 
     def test_get_unique_vertex_values(self):
Index: /inundation/ga/storm_surge/pyvolution/test_quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 1011)
@@ -16,8 +16,8 @@
         a = [0.0, 0.0]
         b = [0.0, 2.0]
-        c = [2.0,0.0]
+        c = [2.0, 0.0]
         d = [0.0, 4.0]
         e = [2.0, 2.0]
-        f = [4.0,0.0]
+        f = [4.0, 0.0]
 
         points = [a, b, c, d, e, f]
@@ -196,4 +196,28 @@
         assert allclose(quantity.centroid_values,
                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
+
+
+    def test_integral(self):
+        quantity = Quantity(self.mesh4)
+
+	#Try constants first	
+	const = 5
+        quantity.set_values(const, location = 'vertices')	
+	#print 'Q', quantity.get_integral()
+		
+	assert allclose(quantity.get_integral(), 
+			self.mesh4.get_area() * const)
+	
+	#Try with a linear function
+        def f(x, y):
+            return x+y
+
+        quantity.set_values(f, location = 'vertices')
+	
+	
+	ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
+        
+	assert allclose (quantity.get_integral(), ref_integral)	
+
 
 
@@ -900,5 +924,5 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test_limiter2')
+    suite = unittest.makeSuite(TestCase,'test')
     #print "restricted test"
     #suite = unittest.makeSuite(TestCase,'test_set_vertex_values_subset')
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1010)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1011)
@@ -1522,5 +1522,184 @@
 
         
-        
+    def test_conservation_1(self):
+        """Test that stage is conserved globally
+	
+	This one uses a flat 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=1
+
+	#IC
+        def x_slope(x, y):
+            return x/3
+
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('friction', 0)	
+        domain.set_quantity('stage', x_slope)	
+	
+        # 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()      
+	
+	#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) 
+	 
+            #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()   
+	    #print 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 
+	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=1
+
+	#IC
+        def x_slope(x, y):
+            return x/3
+
+        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()      
+	
+	#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) 
+	 
+            #FIXME: What would we expect from momentum
+            #xmom = domain.quantities['xmomentum'].get_integral()   
+	    #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 
+	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=1
+
+	#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.4
+                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 
+	#FIXME: What exactly is the significance of this?
+	#Is it a problem that an arbitrary initial condition 
+	#needs to respect the law of the limiters?
+	domain.distribute_to_vertices_and_edges()
+	initial_volume = domain.quantities['stage'].get_integral()         
+		
+        #print initial_volume
+	
+	#print initial_xmom
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
+            volume =  domain.quantities['stage'].get_integral()         
+            #print t, volume, initial_volume
+	    
+	    assert allclose (volume, initial_volume) 
+	 
+            #FIXME: What would we expect from momentum
+            #xmom = domain.quantities['xmomentum'].get_integral()   
+	    #print xmom
+	    #assert allclose (xmom, initial_xmom) 
+        
+	
+	        
     def test_second_order_flat_bed_onestep(self):
 
