Index: branches/numpy/anuga/fit_interpolate/fit.py
===================================================================
--- branches/numpy/anuga/fit_interpolate/fit.py	(revision 6533)
+++ branches/numpy/anuga/fit_interpolate/fit.py	(revision 6553)
@@ -280,5 +280,5 @@
             
             element_found, sigma0, sigma1, sigma2, k = \
-                           search_tree_of_vertices(self.root, self.mesh, x)
+                           search_tree_of_vertices(self.root, x)
             
             if element_found is True:
@@ -301,6 +301,4 @@
                         AtA[j,k] += sigmas[j]*sigmas[k]
             else:
-                # FIXME(Ole): This is the message referred to in ticket:314
-                
                 flag = is_inside_polygon(x,
                                          self.mesh_boundary_polygon,
@@ -310,4 +308,5 @@
                 assert flag is True, msg                
                 
+                # FIXME(Ole): This is the message referred to in ticket:314
                 minx = min(self.mesh_boundary_polygon[:,0])
                 maxx = max(self.mesh_boundary_polygon[:,0])                
@@ -317,7 +316,6 @@
                 msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
                     % (minx, maxx, miny, maxy)
-                
-
-                raise Exception(msg)
+                raise RuntimeError, msg
+
             self.AtA = AtA
 
@@ -375,4 +373,12 @@
                 self.build_fit_subset(points, z, verbose=verbose)
 
+                # FIXME(Ole): I thought this test would make sense here
+                # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
+                # Committed 11 March 2009
+                msg = 'Matrix AtA was not built'
+                assert self.AtA is not None, msg
+                
+                #print 'Matrix was built OK'
+
                 
             point_coordinates = None
@@ -382,8 +388,7 @@
         if point_coordinates is None:
             if verbose: print 'Warning: no data points in fit'
-            #assert self.AtA != None, 'no interpolation matrix'
-            #assert self.Atz != None
-            assert not self.AtA is None, 'no interpolation matrix'
-            assert not self.Atz is None
+            msg = 'No interpolation matrix.'
+            assert self.AtA is not None, msg
+            assert self.Atz is not None
             
             # FIXME (DSG) - do  a message
@@ -471,7 +476,4 @@
                 alpha=DEFAULT_ALPHA,
                 verbose=False,
-                acceptable_overshoot=1.01, 
-                # FIXME: Move to config - this value is assumed in caching test
-                # FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.
                 mesh_origin=None,
                 data_origin=None,
@@ -490,5 +492,4 @@
               'alpha': alpha,
               'verbose': verbose,
-              'acceptable_overshoot': acceptable_overshoot,
               'mesh_origin': mesh_origin,
               'data_origin': data_origin,
@@ -546,5 +547,4 @@
                  alpha=DEFAULT_ALPHA,
                  verbose=False,
-                 acceptable_overshoot=1.01,
                  mesh_origin=None,
                  data_origin=None,
@@ -572,17 +572,8 @@
           alpha: Smoothing parameter.
 
-          acceptable overshoot: NOT IMPLEMENTED
-          controls the allowed factor by which
-          fitted values
-          may exceed the value of input data. The lower limit is defined
-          as min(z) - acceptable_overshoot*delta z and upper limit
-          as max(z) + acceptable_overshoot*delta z
-          
-
           mesh_origin: A geo_reference object or 3-tuples consisting of
               UTM zone, easting and northing.
               If specified vertex coordinates are assumed to be
               relative to their respective origins.
-          
 
           point_attributes: Vector or array of data at the
Index: branches/numpy/anuga/fit_interpolate/interpolate.py
===================================================================
--- branches/numpy/anuga/fit_interpolate/interpolate.py	(revision 6533)
+++ branches/numpy/anuga/fit_interpolate/interpolate.py	(revision 6553)
@@ -484,7 +484,7 @@
             x = point_coordinates[i]
             element_found, sigma0, sigma1, sigma2, k = \
-                           search_tree_of_vertices(self.root, self.mesh, x)
-
-            # Update interpolation matrix A if necessary
+                           search_tree_of_vertices(self.root, x)
+            
+	    # Update interpolation matrix A if necessary
             if element_found is True:
                 # Assign values to matrix A
Index: branches/numpy/anuga/fit_interpolate/search_functions.py
===================================================================
--- branches/numpy/anuga/fit_interpolate/search_functions.py	(revision 6533)
+++ branches/numpy/anuga/fit_interpolate/search_functions.py	(revision 6553)
@@ -8,4 +8,5 @@
 import time
 
+from anuga.utilities.polygon import is_inside_triangle
 from anuga.utilities.numerical_tools import get_machine_precision
 from anuga.config import max_float
@@ -18,13 +19,14 @@
 search_more_cells_time = initial_search_value
 
-#FIXME test what happens if a 
-LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]),
-                        num.array([max_float,max_float]),
-                        num.array([max_float,max_float])),
-                       (num.array([1,1]),
-                        num.array([0,0]),
-                        num.array([-1.1,-1.1]))]]]
+# FIXME(Ole): Could we come up with a less confusing structure?
+LAST_TRIANGLE = [[-10,
+                   (num.array([[max_float, max_float],
+                               [max_float, max_float],
+                               [max_float, max_float]]),
+                    (num.array([1.,1.]),      
+                     num.array([0.,0.]),      
+                     num.array([-1.1,-1.1])))]]
 
-def search_tree_of_vertices(root, mesh, x):
+def search_tree_of_vertices(root, x):
     """
     Find the triangle (element) that the point x is in.
@@ -32,5 +34,4 @@
     Inputs:
         root: A quad tree of the vertices
-        mesh: The underlying mesh
         x:    The point being placed
     
@@ -47,60 +48,21 @@
     global search_more_cells_time
 
-    #Find triangle containing x:
-    element_found = False
-
-    # This will be returned if element_found = False
-    sigma2 = -10.0
-    sigma0 = -10.0
-    sigma1 = -10.0
-    k = -10.0
-
     # Search the last triangle first
-    try:
-        element_found, sigma0, sigma1, sigma2, k = \
-            _search_triangles_of_vertices(mesh, last_triangle, x)
-    except:
-        element_found = False
+    element_found, sigma0, sigma1, sigma2, k = \
+        _search_triangles_of_vertices(last_triangle, x)
                    
-    #print "last_triangle", last_triangle
     if element_found is True:
-        #print "last_triangle", last_triangle
         return element_found, sigma0, sigma1, sigma2, k
 
-    # This was only slightly faster than just checking the
-    # last triangle and it significantly slowed down
-    # non-gridded fitting
-    # If the last element was a dud, search its neighbours
-    #print "last_triangle[0][0]", last_triangle[0][0]
-    #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0])
-    #print "neighbours", neighbours
-    #neighbours = []
-  #   for k in neighbours:
-#         if k >= 0:
-#             tri = mesh.get_vertex_coordinates(k,
-#                                                    absolute=True)
-#             n0 = mesh.get_normal(k, 0)
-#             n1 = mesh.get_normal(k, 1)
-#             n2 = mesh.get_normal(k, 2) 
-#             triangle =[[k,(tri, (n0, n1, n2))]]
-#             element_found, sigma0, sigma1, sigma2, k = \
-#                            _search_triangles_of_vertices(mesh,
-#                                                          triangle, x)
-#             if element_found is True:
-#                 return element_found, sigma0, sigma1, sigma2, k
-            
-    #t0 = time.time()
+    
     # Get triangles in the cell that the point is in.
     # Triangle is a list, first element triangle_id,
     # second element the triangle
     triangles = root.search(x[0], x[1])
+    element_found, sigma0, sigma1, sigma2, k = \
+                   _search_triangles_of_vertices(triangles, x)
+
     is_more_elements = True
     
-    element_found, sigma0, sigma1, sigma2, k = \
-                   _search_triangles_of_vertices(mesh,
-                                                 triangles, x)
-    #search_one_cell_time += time.time()-t0
-    #print "search_one_cell_time",search_one_cell_time
-    #t0 = time.time()
     while not element_found and is_more_elements:
         triangles, branch = root.expand_search()
@@ -109,17 +71,16 @@
             # been searched.  This is the last try
             element_found, sigma0, sigma1, sigma2, k = \
-                           _search_triangles_of_vertices(mesh, triangles, x)
+                           _search_triangles_of_vertices(triangles, x)
             is_more_elements = False
         else:
             element_found, sigma0, sigma1, sigma2, k = \
-                       _search_triangles_of_vertices(mesh, triangles, x)
-        #search_more_cells_time += time.time()-t0
-    #print "search_more_cells_time", search_more_cells_time
+                       _search_triangles_of_vertices(triangles, x)
+                       
         
     return element_found, sigma0, sigma1, sigma2, k
 
 
-def _search_triangles_of_vertices(mesh, triangles, x):
-    """Search for triangle containing x amongs candidate_vertices in mesh
+def _search_triangles_of_vertices(triangles, x):
+    """Search for triangle containing x amongs candidate_vertices in triangles
 
     This is called by search_tree_of_vertices once the appropriate node
@@ -132,34 +93,40 @@
     global last_triangle
     
-    # these statments are needed if triangles is empty
-    #Find triangle containing x:
-    element_found = False
-
-    # This will be returned if element_found = False
+    # These statments are needed if triangles is empty
     sigma2 = -10.0
     sigma0 = -10.0
     sigma1 = -10.0
     k = -10
-    
-    #For all vertices in same cell as point x
+    # For all vertices in same cell as point x
+    element_found = False    
     for k, tri_verts_norms in triangles:
         tri = tri_verts_norms[0]
-        n0, n1, n2 = tri_verts_norms[1]
         # k is the triangle index
         # tri is a list of verts (x, y), representing a tringle
         # Find triangle that contains x (if any) and interpolate
-        element_found, sigma0, sigma1, sigma2 =\
-                       find_triangle_compute_interpolation(tri, n0, n1, n2, x)
-        if element_found is True:
+        
+        # Input check disabled to speed things up.
+        if is_inside_triangle(x, tri, 
+                              closed=True,
+                              check_inputs=False):
+            
+            n0, n1, n2 = tri_verts_norms[1]        
+            sigma0, sigma1, sigma2 =\
+                compute_interpolation_values(tri, n0, n1, n2, x)
+                
+            element_found = True    
+
             # Don't look for any other triangles in the triangle list
-            last_triangle = [[k,tri_verts_norms]]
-            break
+            last_triangle = [[k, tri_verts_norms]]
+            break            
+            
     return element_found, sigma0, sigma1, sigma2, k
 
 
             
-def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
-    """Compute linear interpolation of point x and triangle k in mesh.
-    It is assumed that x belongs to triangle k.max_float
+def compute_interpolation_values(triangle, n0, n1, n2, x):
+    """Compute linear interpolation of point x and triangle.
+    
+    n0, n1, n2 are normal to the tree edges.
     """
 
@@ -167,54 +134,9 @@
     xi0, xi1, xi2 = triangle
 
-    # this is where we can call some fast c code.
-      
-    # Integrity check - machine precision is too hard
-    # so we use hardwired single precision 
-    epsilon = 1.0e-6
-    
-    if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
-        # print "max(xi0[0], xi1[0], xi2[0])", max(xi0[0], xi1[0], xi2[0])
-        return False,0,0,0
-    if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:
-        return False,0,0,0
-    if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:
-        return False,0,0,0
-    if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:
-        return False,0,0,0
-    
-    # machine precision on some machines (e.g. nautilus)
-    epsilon = get_machine_precision() * 2
-    
-    # Compute interpolation - return as soon as possible
-    #  print "(xi0-xi1)", (xi0-xi1)
-    # print "n0", n0
-    # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)
-    
     sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
-    if sigma0 < -epsilon:
-        return False,0,0,0
     sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
-    if sigma1 < -epsilon:
-        return False,0,0,0
     sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
-    if sigma2 < -epsilon:
-        return False,0,0,0
-    
-    # epsilon = 1.0e-6
-    # we want to speed this up, so don't do assertions
-    #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
-    #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
-    #      %(delta, epsilon)
-    #assert delta < epsilon, msg
 
-
-    # Check that this triangle contains the data point
-    # Sigmas are allowed to get negative within
-    # machine precision on some machines (e.g. nautilus)
-    #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
-    #    element_found = True
-    #else:
-    #    element_found = False 
-    return True, sigma0, sigma1, sigma2
+    return sigma0, sigma1, sigma2
 
 def set_last_triangle():
Index: branches/numpy/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- branches/numpy/anuga/fit_interpolate/test_search_functions.py	(revision 6533)
+++ branches/numpy/anuga/fit_interpolate/test_search_functions.py	(revision 6553)
@@ -5,4 +5,5 @@
 from search_functions import search_tree_of_vertices, set_last_triangle
 from search_functions import _search_triangles_of_vertices
+from search_functions import compute_interpolation_values
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -10,5 +11,8 @@
 from anuga.utilities.polygon import is_inside_polygon
 from anuga.utilities.quad import build_quadtree, Cell
-
+from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle    
+
+import numpy as num
 
 class Test_search_functions(unittest.TestCase):
@@ -49,10 +53,12 @@
 
         x = [0.2, 0.7]
-        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+        found, s0, s1, s2, k = search_tree_of_vertices(root, x)
         assert k == 1 # Triangle one
         assert found is True
 
     def test_bigger(self):
-        """test_larger mesh
+        """test_bigger
+        
+        test larger mesh
         """
 
@@ -69,6 +75,7 @@
                   [0.1,0.9], [0.4,0.6], [0.9,0.1],
                   [10, 3]]:
-                
-            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+            
+            found, s0, s1, s2, k = search_tree_of_vertices(root, 
+                                                           ensure_numeric(x))
 
             if k >= 0:
@@ -102,5 +109,5 @@
                       [10, 3]]:
                 
-                found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+                found, s0, s1, s2, k = search_tree_of_vertices(root, x)
 
                 if k >= 0:
@@ -111,5 +118,6 @@
                     assert found is False                
 
-                
+            
+                if m == k == 0: return    
 
     def test_underlying_function(self):
@@ -124,11 +132,10 @@
 
         # One point
-        x = [0.5, 0.5]
+        x = ensure_numeric([0.5, 0.5])
         candidate_vertices = root.search(x[0], x[1])
 
         #print x, candidate_vertices
         found, sigma0, sigma1, sigma2, k = \
-               _search_triangles_of_vertices(mesh,
-                                             candidate_vertices,
+               _search_triangles_of_vertices(candidate_vertices,
                                              x)
 
@@ -151,7 +158,6 @@
             #print x, candidate_vertices
             found, sigma0, sigma1, sigma2, k = \
-                   _search_triangles_of_vertices(mesh,
-                                                 candidate_vertices,
-                                                 x)
+                   _search_triangles_of_vertices(candidate_vertices,
+                                                 ensure_numeric(x))
             if k >= 0:
                 V = mesh.get_vertex_coordinates(k) # nodes for triangle k
@@ -199,17 +205,54 @@
         x = [2.5, 1.5]
         element_found, sigma0, sigma1, sigma2, k = \
-                       search_tree_of_vertices(root, mesh, x)
+                       search_tree_of_vertices(root, x)
         # One point
         x = [3.00005, 2.999994]
         element_found, sigma0, sigma1, sigma2, k = \
-                       search_tree_of_vertices(root, mesh, x)
+                       search_tree_of_vertices(root, x)
         assert element_found is True
         assert k == 1
         
-################################################################################
-
+        
+    def test_compute_interpolation_values(self):
+        """test_compute_interpolation_values
+        
+        Test that interpolation values are correc
+        
+        This test used to check element_found as output from 
+        find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
+        and that failed before 18th March 2009.
+        
+        Now this function no longer returns this flag, so the test
+        is merely checknig the sigmas.
+        
+        """
+        
+        triangle = num.array([[306951.77151059, 6194462.14596986],
+                              [306952.58403545, 6194459.65001246],
+                              [306953.55109034, 6194462.0041216]])
+
+
+        n0 = [0.92499377, -0.37998227]
+        n1 = [0.07945684,  0.99683831]
+        n2 = [-0.95088404, -0.30954732]
+        
+        x = [306953.344, 6194461.5]
+        
+        # Test that point is indeed inside triangle
+        assert is_inside_polygon(x, triangle, 
+                                 closed=True, verbose=False)
+        assert is_inside_triangle(x, triangle, 
+                                  closed=True, verbose=False)                                 
+        
+        sigma0, sigma1, sigma2 = \
+            compute_interpolation_values(triangle, n0, n1, n2, x)
+            
+        msg = 'Point which is clearly inside triangle was not found'
+        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
+
+#-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_search_functions,'test')
-    runner = unittest.TextTestRunner() #verbosity=1)
+    suite = unittest.makeSuite(Test_search_functions, 'test')
+    runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)
     
