Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 6184)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 6185)
@@ -214,73 +214,4 @@
         self._point_coordinates = None # FIXME (Ole): Probably obsolete
         self.interpolation_matrices = {} # Store precomputed matrices
-
-
-    ##
-    # @brief Interpolate linearly from polyline nodes to midpoints of triangles.
-    # @param f The data on the polyline nodes.
-    # @param vertex_coordinates ??
-    # @param gauge_neighbour_id ??
-    # @param point_coordinates ??
-    # @param verbose True if this function is to be verbose.
-    def interpolate_polyline(self,
-                             f,
-                             vertex_coordinates,
-                             gauge_neighbour_id,
-                             point_coordinates=None,
-                             verbose=False):
-        """Interpolate linearly between values f on nodes (vertex coordinates)
-        of a polyline to midpoints of triangles of boundary.
-
-        f is the data on the polyline nodes.
-
-        The mesh values representing a smooth surface are
-        assumed to be specified in f.
-
-        Inputs:
-          f: Vector or array of data at the polyline nodes.
-              If f is an array, interpolation will be done for each column as
-              per underlying matrix-matrix multiplication
-          point_coordinates: Interpolate polyline data to these positions.
-              List of coordinate pairs [x, y] of
-	      data points or an nx2 Numeric array or a Geospatial_data object
-              
-	Output:
-	  Interpolated values at inputted points (z).
-        """
-        
-        #print f
-        #print vertex_coordinates
-        #print gauge_neighbour_id
-        #print point_coordinates
-
-
-        # FIXME: There is an option of passing a tolerance into this
-        
-        if isinstance(point_coordinates, Geospatial_data):
-            point_coordinates = point_coordinates.get_data_points(absolute=True)
- 
-        z = num.ones(len(point_coordinates), num.Float)
-
-        # input sanity check
-        msg = 'point coordinates are not given (interpolate.py)'
-        assert point_coordinates is not None, msg
-        msg = 'function value must be specified at every interpolation node'
-        assert f.shape[0]==vertex_coordinates.shape[0], msg
-        msg = 'Must define function value at one or more nodes'
-        assert f.shape[0]>0, msg
-
-        n = f.shape[0]
-        if n == 1:
-            z = f*z
-            msg = 'Polyline contained only one point. I need more. ', str(f)
-            raise Exception, msg
-            
-        # FIXME (John): add unit test for only One vertex point.
-        #               Exception should be thrown.
-        
-        elif n > 1:
-            _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
-            
-        return z
 
     
@@ -576,4 +507,72 @@
         
         
+
+##
+# @brief Interpolate linearly from polyline nodes to midpoints of triangles.
+# @param f The data on the polyline nodes.
+# @param vertex_coordinates ??
+# @param gauge_neighbour_id ??
+# @param point_coordinates ??
+# @param verbose True if this function is to be verbose.
+def interpolate_polyline(f,
+                         vertex_coordinates,
+                         gauge_neighbour_id,
+                         point_coordinates=None,
+                         verbose=False):
+    """Interpolate linearly between values f on nodes (vertex coordinates)
+    of a polyline to midpoints of triangles of boundary.
+
+    f is the data on the polyline nodes.
+
+    The mesh values representing a smooth surface are
+    assumed to be specified in f.
+
+    Inputs:
+      f: Vector or array of data at the polyline nodes.
+          If f is an array, interpolation will be done for each column as
+          per underlying matrix-matrix multiplication
+      point_coordinates: Interpolate polyline data to these positions.
+          List of coordinate pairs [x, y] of
+          data points or an nx2 Numeric array or a Geospatial_data object
+          
+    Output:
+      Interpolated values at inputted points (z).
+    """
+    
+    #print f
+    #print vertex_coordinates
+    #print gauge_neighbour_id
+    #print point_coordinates
+
+
+    # FIXME: There is an option of passing a tolerance into this
+    
+    if isinstance(point_coordinates, Geospatial_data):
+        point_coordinates = point_coordinates.get_data_points(absolute=True)
+
+    z = num.ones(len(point_coordinates), num.Float)
+
+    # input sanity check
+    msg = 'point coordinates are not given (interpolate.py)'
+    assert point_coordinates is not None, msg
+    msg = 'function value must be specified at every interpolation node'
+    assert f.shape[0]==vertex_coordinates.shape[0], msg
+    msg = 'Must define function value at one or more nodes'
+    assert f.shape[0]>0, msg
+
+    n = f.shape[0]
+    if n == 1:
+        z = f*z
+        msg = 'Polyline contained only one point. I need more. ', str(f)
+        raise Exception, msg
+        
+    # FIXME (John): add unit test for only One vertex point.
+    #               Exception should be thrown.
+    
+    elif n > 1:
+        _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
+        
+    return z
+
         
 def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id):
@@ -1051,11 +1050,10 @@
                                                       verbose=False) # No clutter
                     elif triangles is None and vertex_coordinates is not None:
-                        result = \
-                            interpol.interpolate_polyline(Q,
-                                                          vertex_coordinates,
-                                                          gauge_neighbour_id,
-                                                          point_coordinates=\
+                        result = interpolate_polyline(Q,
+                                                      vertex_coordinates,
+                                                      gauge_neighbour_id,
+                                                      point_coordinates=\
                                                           self.interpolation_points)
-
+                        
                     #assert len(result), len(interpolation_points)
                     self.precomputed_values[name][i, :] = result
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 6184)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 6185)
@@ -1853,31 +1853,4 @@
         John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
         """
-
-        a = [-1.0, 0.0]
-        b = [3.0, 4.0]
-        c = [4.0, 1.0]
-        d = [-3.0, 2.0] #3
-        e = [-1.0, -2.0]
-        f = [1.0, -2.0] #5
-
-        vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
-
-
-        point_coords = [[-2.0, 2.0],
-                        [-1.0, 1.0],
-                        [0.0, 2.0],
-                        [1.0, 1.0],
-                        [2.0, 1.0],
-                        [0.0, 0.0],
-                        [1.0, 0.0],
-                        [0.0, -1.0],
-                        [-0.2, -0.5],
-                        [-0.9, -1.5],
-                        [0.5, -1.9],
-                        [3.0, 1.0]]
-
-        interp = Interpolate(vertices, triangles)
-                
         
         f = num.array([58.06150614, 58.06150614, 58.06150614])
@@ -1896,5 +1869,5 @@
         z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
         
-        z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
         assert num.allclose(z, z_ref)
         
@@ -1902,5 +1875,5 @@
         f = num.array([58.06150614, 158.06150614, 258.06150614])
         z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]        
-        z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
         assert num.allclose(z, z_ref)        
                        
@@ -1919,5 +1892,5 @@
         gauge_neighbour_id = [1, 2, -1]
                                                
-        z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
+        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
         z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
         #print z
