Index: /anuga_core/source/anuga/fit_interpolate/search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 6543)
+++ /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 6544)
@@ -19,11 +19,11 @@
 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.Int),      #array default#
-                        num.array([0,0], num.Int),      #array default#
-                        num.array([-1.1,-1.1]))]]]
+LAST_TRIANGLE = [[-10,
+                   (num.array([[max_float,max_float],
+                               [max_float,max_float],
+                               [max_float,max_float]]),
+                    (num.array([1,1], num.Int),      #array default#
+                     num.array([0,0], num.Int),      #array default#
+                     num.array([-1.1,-1.1])))]]
 
 def search_tree_of_vertices(root, mesh, x):
@@ -48,8 +48,6 @@
     global search_more_cells_time
 
-    #Find triangle containing x:
-    element_found = False
-
     # This will be returned if element_found = False
+    element_found = False    
     sigma2 = -10.0
     sigma0 = -10.0
@@ -61,35 +59,13 @@
         element_found, sigma0, sigma1, sigma2, k = \
             _search_triangles_of_vertices(mesh, last_triangle, x)
+            
     except:
+        print 'This should never happen:', last_triangle
         element_found = False
                    
-    #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,
@@ -97,5 +73,4 @@
     triangles = root.search(x[0], x[1])
     is_more_elements = True
-    
     element_found, sigma0, sigma1, sigma2, k = \
                    _search_triangles_of_vertices(mesh,
@@ -133,15 +108,12 @@
     global last_triangle
     
-    # these statments are needed if triangles is empty
-    #Find triangle containing x:
+    # These statments are needed if triangles is empty
     element_found = False
-
-    # This will be returned if element_found = False
     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
     for k, tri_verts_norms in triangles:
         tri = tri_verts_norms[0]
@@ -150,23 +122,28 @@
         # tri is a list of verts (x, y), representing a tringle
         # Find triangle that contains x (if any) and interpolate
-        if is_inside_triangle(x, tri, closed=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]        
-            element_found, sigma0, sigma1, sigma2 =\
-                find_triangle_compute_interpolation(tri, n0, n1, n2, x)
-        else:
-            element_found = False
+            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            
 
             
-        if element_found is True:
-            # Don't look for any other triangles in the triangle list
-            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.
     """
 
@@ -174,60 +151,11 @@
     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
+    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
+    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
+    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
+
+    return sigma0, sigma1, sigma2
+
     
-    #if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
-    #    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:
-        #print 'sigma0', sigma0
-        #sigma0 = 0.0
-        #return False,0,0,0
-    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
-    #if sigma1 < -epsilon:
-        #print 'sigma1', sigma1    
-        #sigma1 = 0.0        
-        #return False,0,0,0
-    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
-    #if sigma2 < -epsilon:
-        #print 'sigma2', sigma2    
-        #sigma2 = 0.0        
-        #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
-
 def set_last_triangle():
     global last_triangle
Index: /anuga_core/source/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 6543)
+++ /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 6544)
@@ -5,5 +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 find_triangle_compute_interpolation
+from search_functions import compute_interpolation_values
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -11,4 +11,5 @@
 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    
 
@@ -52,10 +53,14 @@
 
         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, 
+                                                       mesh, 
+                                                       ensure_numeric(x))
         assert k == 1 # Triangle one
         assert found is True
 
     def test_bigger(self):
-        """test_larger mesh
+        """test_bigger
+        
+        test larger mesh
         """
 
@@ -72,6 +77,9 @@
                   [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, 
+                                                           mesh, 
+                                                           ensure_numeric(x))
 
             if k >= 0:
@@ -95,5 +103,4 @@
         mesh.check_integrity()
 
-        
         for m in range(8):
             root = build_quadtree(mesh, max_points_per_cell = m)
@@ -105,14 +112,20 @@
                       [10, 3]]:
                 
-                found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+                found, s0, s1, s2, k = search_tree_of_vertices(root, 
+                                                               mesh, 
+                                                               x)
 
                 if k >= 0:
                     V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+
+                    assert is_inside_triangle(x, V, closed=True)
                     assert is_inside_polygon(x, V)
+                    
                     assert found is True
                 else:
                     assert found is False                
 
-                
+            
+                if m == k == 0: return    
 
     def test_underlying_function(self):
@@ -127,8 +140,8 @@
 
         # 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
+        # print x, candidate_vertices
         found, sigma0, sigma1, sigma2, k = \
                _search_triangles_of_vertices(mesh,
@@ -156,5 +169,5 @@
                    _search_triangles_of_vertices(mesh,
                                                  candidate_vertices,
-                                                 x)
+                                                 ensure_numeric(x))
             if k >= 0:
                 V = mesh.get_vertex_coordinates(k) # nodes for triangle k
@@ -211,8 +224,16 @@
         
         
-    def test_triangle_compute_interpolation(self):
-        """test_triangle_compute_interpolation
-        
-        Test that triangle can be found if point is inside it
+    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.
+        
         """
         
@@ -232,16 +253,11 @@
                                  closed=True, verbose=False)
         assert is_inside_triangle(x, triangle, 
-                                  closed=True, verbose=False)                                 
-        
-        element_found, sigma0, sigma1, sigma2 = \
-            find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
+                                  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 element_found is True, msg
-        
-        #print sigma0, sigma1, sigma2
-        #print sigma0 + sigma1 + sigma2
-                
-        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6
+        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
 
 #-------------------------------------------------------------
Index: /anuga_core/source/anuga/utilities/polygon.py
===================================================================
--- /anuga_core/source/anuga/utilities/polygon.py	(revision 6543)
+++ /anuga_core/source/anuga/utilities/polygon.py	(revision 6544)
@@ -232,12 +232,20 @@
     """
 
+    triangle = ensure_numeric(triangle)        
+    point = ensure_numeric(point, num.Float)    
+    
     if check_inputs is True:
-        triangle = ensure_numeric(triangle)    
-
-        point = ensure_numeric(point, num.Float) # Convert point to array of points
         msg = 'is_inside_triangle must be invoked with one point only'
         assert num.allclose(point.shape, [2]), msg
     
-        
+    
+    # Use C-implementation
+    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
+    
+    
+
+    # FIXME (Ole): The rest of this function has been made 
+    # obsolete by the C extension.
+    
     # Quickly reject points that are clearly outside
     if point[0] < min(triangle[:,0]): return False 
@@ -245,12 +253,4 @@
     if point[1] < min(triangle[:,1]): return False
     if point[1] > max(triangle[:,1]): return False        
-
-    # Use C-implementation
-    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
-    
-    
-
-    # FIXME (Ole): The rest of this function has been made 
-    # obsolete by the C extension.
         
     # Start search    
Index: /anuga_core/source/anuga/utilities/polygon_ext.c
===================================================================
--- /anuga_core/source/anuga/utilities/polygon_ext.c	(revision 6543)
+++ /anuga_core/source/anuga/utilities/polygon_ext.c	(revision 6544)
@@ -254,6 +254,29 @@
   double a00, a10, a01, a11, b0, b1;
   double denom, alpha, beta;
+  
+  double x, y; // Point coordinates
   int i, j, res;
 
+  x = point[0];
+  y = point[1];
+  
+  // Quickly reject points that are clearly outside
+  if ((x < triangle[0]) && 
+      (x < triangle[2]) && 
+      (x < triangle[4])) return 0;       
+      
+  if ((x > triangle[0]) && 
+      (x > triangle[2]) && 
+      (x > triangle[4])) return 0;             
+  
+  if ((y < triangle[1]) && 
+      (y < triangle[3]) && 
+      (y < triangle[5])) return 0;       
+      
+  if ((y > triangle[1]) && 
+      (y > triangle[3]) && 
+      (y > triangle[5])) return 0;             
+  
+  
   // v0 = C-A 
   v0x = triangle[4]-triangle[0]; 
@@ -274,6 +297,6 @@
   if (fabs(denom) > 0.0) {
     // v = point-A  
-    vx = point[0] - triangle[0]; 
-    vy = point[1] - triangle[1];     
+    vx = x - triangle[0]; 
+    vy = y - triangle[1];     
     
     b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)        
@@ -291,5 +314,5 @@
     for (i=0; i<3; i++) {
       j = (i+1) % 3; // Circular index into triangle array
-      res = __point_on_line(point[0], point[1],
+      res = __point_on_line(x, y,
                             triangle[2*i], triangle[2*i+1], 
                             triangle[2*j], triangle[2*j+1], 			    
Index: /anuga_core/source/anuga/utilities/test_polygon.py
===================================================================
--- /anuga_core/source/anuga/utilities/test_polygon.py	(revision 6543)
+++ /anuga_core/source/anuga/utilities/test_polygon.py	(revision 6544)
@@ -1896,5 +1896,11 @@
         assert res is False
         
-        
+
+        res = is_inside_triangle([10, 3], 
+                                 [[0.1, 0.],
+                                  [0.1, 0.08333333],
+                                  [0., 0.]])
+        assert res is False
+                
                 
 #-------------------------------------------------------------
