Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5728)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5729)
@@ -369,5 +369,5 @@
 
     def get_quantity(self, name, location='vertices', indices = None):
-        """Get quantity object.
+        """Get pointer to quantity object.
 
         name: Name of quantity
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 5728)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 5729)
@@ -1113,4 +1113,5 @@
         
 
+        
 
 class Triangle_intersection:
@@ -1147,2 +1148,25 @@
     
         return s
+
+def segment_midpoints(segments):
+    """Calculate midpoints of all segments
+    
+    Inputs:
+       segments: List of instances of class Segment
+       
+    Ouputs:
+       midpoints: List of points
+    """
+    
+    midpoints = []
+    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
+    for segment in segments:
+        assert isinstance(segment, Triangle_intersection), msg
+        
+        midpoint = sum(array(segment.segment))/2
+        midpoints.append(midpoint)
+
+    return midpoints
+    
+    
+    
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5728)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5729)
@@ -753,10 +753,10 @@
                                         verbose=False,
                                         use_cache=False):
-        # FIXME: Use this function for the time being. Later move code in here
+        """ Set values based on geo referenced geospatial data object.
+        """
 
         points = geospatial_data.get_data_points(absolute=False)
         values = geospatial_data.get_attributes()
         data_georef = geospatial_data.get_geo_reference()
-
 
 
@@ -1007,9 +1007,35 @@
 
     def get_interpolated_values(self, interpolation_points):
-
-        # Interpolation object based on internal (discontinuous triangles)
+        """ Get values at interpolation points
+        
+        If interpolation points have been given previously, the 
+        associated matrices will be reused to save time.
+        
+        The argument interpolation points must be given as either a list of absolute UTM coordinates or
+        a geospatial data object.
+        """
+
+        # FIXME (Ole): Could do with an input check (should be generalised and used widely)
+        # That will check that interpolation points is either a list of points, Nx2 array, or geospatial
+        
+        # Ensure points are converted to coordinates relative to mesh origin
+        # FIXME (Ole): This could all be refactored using the 'change_points_geo_ref' method 
+        # of Class geo_reference. The purpose is to make interpolation points relative
+        # to the mesh origin.
+        #
+        # Speed is also a consideration here.
+        
+        if isinstance(interpolation_points, Geospatial_data):        
+            # Ensure interpolation points are in absolute UTM coordinates
+            interpolation_points = interpolation_points.get_data_points(absolute=True)
+            
+        # Reconcile interpolation points with georeference of domain
+        interpolation_points = self.domain.geo_reference.get_relative(interpolation_points) 
+        interpolation_points = ensure_numeric(interpolation_points)
+        
+        # Get internal (discontinuous) triangles for use with the interpolation object.
         x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
                                                                 smooth=False)
-        # FIXME: This concat should roll into get_vertex_values
+        # FIXME (Ole): This concat should roll into get_vertex_values
         vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
                                          axis=1)
@@ -1020,6 +1046,8 @@
             I = self.interpolation_object
 
-            if allclose(interpolation_points, I._point_coordinates):
-                can_reuse = True
+            if allclose(interpolation_points.shape, 
+                        I._point_coordinates.shape):
+                if allclose(interpolation_points, I._point_coordinates):
+                    can_reuse = True
                 
 
@@ -1047,11 +1075,15 @@
 
         return X, Compatible list, Numeric array (see below)
-        interpolation_points: List of x, y coordinates where value is
-        sought (using interpolation). If points are given, values of
-        location and indices are ignored
-        
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroids
-                  and unique vertices. Default is 'vertices'
+        
+        Inputs:
+           interpolation_points: List of x, y coordinates where value is
+                                 sought (using interpolation). If points 
+                                 are given, values of location and indices 
+                                 are ignored. Assume either absolute UTM
+                                 coordinates or geospatial data object.
+        
+           location: Where values are to be stored.
+                     Permissible options are: vertices, edges, centroids
+                     and unique vertices. Default is 'vertices'
 
 
@@ -1075,4 +1107,5 @@
         internal ordering.
         """
+        
         from Numeric import take
 
@@ -1086,5 +1119,8 @@
         
         # FIXME (Ole): Consider deprecating 'edges' - but not if it is used
-        # elsewhere in ANUGA.
+        # elsewhere in ANUGA. 
+        # Edges have already been deprecated in set_values, see changeset:5521,
+        # but *might* be useful in get_values. Any thoughts anyone?
+        
         if location not in ['vertices', 'centroids', 'edges',
                             'unique vertices']:
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 5728)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 5729)
@@ -2354,4 +2354,58 @@
 
 
+    def test_get_interpolated_values_with_georef(self):
+    
+        zone = 56
+        xllcorner = 308500
+        yllcorner = 6189000
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices,
+                        geo_reference=Geo_reference(zone,xllcorner,yllcorner))
+
+        quantity = Quantity(domain)
+        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
+
+        #First pick one point (and turn it into absolute coordinates)
+        x, y = 2.0/3, 8.0/3
+        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
+        assert allclose(v, 6)
+        
+
+        # Then another to test that algorithm won't blindly
+        # reuse interpolation matrix 
+        x, y = 4.0/3, 4.0/3
+        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
+        assert allclose(v, 4)        
+        
+        # Try two points
+        pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 
+               [4.0/3 + xllcorner, 4.0/3 + yllcorner]]         
+        v = quantity.get_values(interpolation_points=pts)
+        assert allclose(v, [6, 4])               
+        
+        # Test it using the geospatial data format with absolute input points and default georef
+        pts = Geospatial_data(data_points=pts)
+        v = quantity.get_values(interpolation_points=pts)
+        assert allclose(v, [6, 4])                                
+        
+        
+        # Test it using the geospatial data format with relative input points
+        pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 
+                              geo_reference=Geo_reference(zone,xllcorner,yllcorner))
+        v = quantity.get_values(interpolation_points=pts)
+        assert allclose(v, [6, 4])                        
+        
+        
+        
 
     def test_getting_some_vertex_values(self):
Index: /anuga_core/source/anuga/geospatial_data/geospatial_data.py
===================================================================
--- /anuga_core/source/anuga/geospatial_data/geospatial_data.py	(revision 5728)
+++ /anuga_core/source/anuga/geospatial_data/geospatial_data.py	(revision 5729)
@@ -1543,5 +1543,8 @@
                     alpha = alpha)
 
-        points_geo=domain.geo_reference.change_points_geo_ref(points)
+        
+        # Convert points to geospatial data for use with get_values below
+        points_geo = Geospatial_data(points, domain.geo_reference)
+        
         #returns the predicted elevation of the points that were "split" out 
         #of the original data set for one particular alpha
Index: /anuga_core/source/anuga/geospatial_data/test_geospatial_data.py
===================================================================
--- /anuga_core/source/anuga/geospatial_data/test_geospatial_data.py	(revision 5728)
+++ /anuga_core/source/anuga/geospatial_data/test_geospatial_data.py	(revision 5729)
@@ -1892,5 +1892,5 @@
             os.remove(filename)
             
-    #        print value, alpha
+            # print value, alpha
             assert (alpha==0.01)
 
@@ -1976,6 +1976,6 @@
 
     #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')
-#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
-#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
+    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
+    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
     suite = unittest.makeSuite(Test_Geospatial_data, 'test')
     runner = unittest.TextTestRunner() #verbosity=2)
Index: /anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 5728)
+++ /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 5729)
@@ -100,4 +100,6 @@
 from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
      remove_lone_verts, sww2timeseries, get_centroid_values
+     
+from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
 from anuga.load_mesh.loadASCII import export_mesh_file
 from anuga.utilities.polygon import intersection
@@ -6370,19 +6372,12 @@
     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
     
-        
-    # Then store for each triangle the length of the intersecting segment(s),
-    # right hand normal(s) and midpoints as interpolation_points.
-    interpolation_points = []
-    for segment in segments:
-        midpoint = sum(array(segment.segment))/2
-        interpolation_points.append(midpoint)
-
-
+    # Get midpoints
+    interpolation_points = segment_midpoints(segments)       
+
+    # Interpolate
     if verbose:
         print 'Interpolating - ',        
         print 'total number of interpolation points = %d'\
               %len(interpolation_points)
-              
-              
     
     I = Interpolation_function(time,
@@ -6445,5 +6440,5 @@
     for t in time:
         total_flow=0
-        for i, p in enumerate(interpolation_points):
+        for i in range(len(interpolation_points)):
             elevation, stage, uh, vh = interpolation_function(t, point_id=i)
             normal = segments[i].normal
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5728)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5729)
@@ -83,5 +83,5 @@
 from Numeric import compress, arange
 
-
+from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
 from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
@@ -374,4 +374,68 @@
         return self.get_quantity('elevation').\
                get_maximum_location(indices=wet_elements)    
+               
+               
+               
+               
+    def get_flow_through_cross_section(self, polyline,
+                                       verbose=False):               
+        """Get the total flow through an arbitrary poly line.        
+        
+        This is a run-time equivalent of the function with same name in data_manager.py
+        
+        Input:
+            polyline: Representation of desired cross section - it may contain 
+                      multiple sections allowing for complex shapes. Assume 
+                      absolute UTM coordinates.
+                      Format [[x0, y0], [x1, y1], ...]        
+                  
+        Output:        
+            Q: Total flow [m^3/s] across given segments.
+        
+         
+        """       
+        
+        # Adjust polyline to mesh spatial origin
+        polyline = self.geo_reference.get_relative(polyline)
+        
+        # Find all intersections and associated triangles.
+        segments = self.get_intersecting_segments(polyline, verbose=verbose)
+
+        msg = 'No segments found'
+        assert len(segments) > 0, msg
+        
+        # Get midpoints
+        midpoints = segment_midpoints(segments)       
+        
+        # FIXME (Ole): HACK - need to make midpoints Geospatial instances
+        midpoints = self.geo_reference.get_absolute(midpoints)        
+        
+        # Compute flow        
+        if verbose: print 'Computing flow through specified cross section'
+        
+        # Get interpolated values
+        xmomentum = self.get_quantity('xmomentum')
+        ymomentum = self.get_quantity('ymomentum')        
+        
+        uh = xmomentum.get_values(interpolation_points=midpoints)
+        vh = ymomentum.get_values(interpolation_points=midpoints)        
+        
+        # Compute and sum flows across each segment
+        total_flow=0
+        for i in range(len(uh)):
+            
+            # Inner product of momentum vector with segment normal [m^2/s] 
+            normal = segments[i].normal
+            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 
+                
+            # Flow across this segment [m^3/s]
+            segment_flow = normal_momentum*segments[i].length
+
+            # Accumulate
+            total_flow += segment_flow
+
+            
+        return total_flow
+        
 
     def check_integrity(self):
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5728)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5729)
@@ -10073,7 +10073,7 @@
         u = 2 m/s
         h = 1 m
-        w = 5 m (width of channel)
-
-        q = u*h*w = 10 m^3/s
+        w = 3 m (width of channel)
+
+        q = u*h*w = 6 m^3/s
 
         #---------- First run without geo referencing        
@@ -10112,5 +10112,5 @@
 
         Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet: 
+        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet: 
 
 
@@ -10199,7 +10199,7 @@
         u = 2 m/s
         h = 2 m
-        w = 5 m (width of channel)
-
-        q = u*h*w = 20 m^3/s
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
 
 
@@ -10243,5 +10243,5 @@
 
         Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet: 
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
 
 
@@ -10300,7 +10300,7 @@
         u = 2 m/s
         h = 1 m
-        w = 5 m (width of channel)
-
-        q = u*h*w = 10 m^3/s
+        w = 3 m (width of channel)
+
+        q = u*h*w = 6 m^3/s
         Es = h + 0.5*v*v/g  # Specific energy head [m]
         Et = w + 0.5*v*v/g  # Total energy head [m]        
@@ -10346,5 +10346,5 @@
 
         Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet: 
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
 
         
Index: /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 5728)
+++ /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 5729)
@@ -1406,4 +1406,104 @@
             
         
+        
+    def test_get_flow_through_cross_section_with_geo(self):
+        """test_get_flow_through_cross_section(self):
+
+        Test that the total flow through a cross section can be
+        correctly obtained at run-time from the ANUGA domain.
+        
+        This test creates a flat bed with a known flow through it and tests
+        that the function correctly returns the expected flow.
+
+        The specifics are
+        e = -1 m
+        u = 2 m/s
+        h = 2 m
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
+
+
+        This run tries it with georeferencing and with elevation = -1
+        
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        # Setup
+        from mesh_factory import rectangular
+
+        # Create basic mesh (20m x 3m)
+        width = 3
+        length = 20
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width,
+                                                 length, width)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference=Geo_reference(56,308500,6189000))
+
+        domain.default_order = 2
+        domain.set_quantities_to_be_stored(None)
+
+
+        e = -1.0
+        w = 1.0
+        h = w-e
+        u = 2.0
+        uh = u*h
+
+        Br = Reflective_boundary(domain)     # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
+
+
+        # Initial conditions
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
+        domain.set_quantity('xmomentum', uh)
+        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+        
+        
+        # Interpolation points down the middle
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        interpolation_points = domain.geo_reference.get_absolute(I)        
+        
+        # Shortcuts to quantites
+        stage = domain.get_quantity('stage')        
+        xmomentum = domain.get_quantity('xmomentum')       
+        ymomentum = domain.get_quantity('ymomentum')               
+
+        for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
+            # Check that quantities are they should be in the interior
+
+            w_t = stage.get_values(interpolation_points)            
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)            
+            
+            assert allclose(w_t, w)
+            assert allclose(uh_t, uh)            
+            assert allclose(vh_t, 0.0)                        
+            
+            
+            # Check flows through the middle
+            for i in range(5):
+                x = length/2. + i*0.23674563 # Arbitrary
+                cross_section = [[x, 0], [x, width]]
+
+                cross_section = domain.geo_reference.get_absolute(cross_section)            
+                Q = domain.get_flow_through_cross_section(cross_section,
+                                                          verbose=False)
+
+                assert allclose(Q, uh*width)
+
+
+            
+        
+        
 
     def test_another_runup_example(self):
@@ -1431,8 +1531,8 @@
         domain = Domain(points, vertices, boundary) # Create domain
         domain.set_quantities_to_be_stored(None)
-        domain.set_maximum_allowed_speed(100) #
+        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
         
         # FIXME (Ole): Need tests where this is commented out
-        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
+        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
         domain.H0 = 0 # Backwards compatibility (6/2/7)
         domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
@@ -5805,7 +5905,7 @@
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Shallow_Water,'test')
-
-    #suite = unittest.makeSuite(Test_Shallow_Water,'test_bedslope_problem_first_order_moresteps')    
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test')
+
+    suite = unittest.makeSuite(Test_Shallow_Water,'test_get_flow_through_cross_section_with_geo')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
Index: /anuga_core/source/anuga/utilities/numerical_tools.py
===================================================================
--- /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 5728)
+++ /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 5729)
@@ -313,4 +313,6 @@
 
     return bins
+    
+
 
 def get_machine_precision():
