Index: /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py	(revision 7717)
+++ /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py	(revision 7718)
@@ -23,7 +23,4 @@
 import profile , pstats
 from math import sqrt
-
-from anuga.fit_interpolate.search_functions import search_times, \
-     reset_search_times
 
 from anuga.fit_interpolate.interpolate import Interpolate
@@ -262,6 +259,6 @@
         # backing up.
         
-        search_one_cell_time, search_more_cells_time = search_times()
-        reset_search_times()
+        #search_one_cell_time, search_more_cells_time = search_times()
+        #reset_search_times()
         #print "bench - search_one_cell_time",search_one_cell_time
         #print "bench - search_more_cells_time", search_more_cells_time
Index: /anuga_core/source/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7717)
+++ /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7718)
@@ -3,5 +3,5 @@
 
 import unittest
-from mesh_quadtree import MeshQuadtree, compute_interpolation_values
+from mesh_quadtree import MeshQuadtree
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -228,41 +228,4 @@
         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
 
 #-------------------------------------------------------------
Index: /anuga_core/source/anuga/geometry/mesh_quadtree.py
===================================================================
--- /anuga_core/source/anuga/geometry/mesh_quadtree.py	(revision 7717)
+++ /anuga_core/source/anuga/geometry/mesh_quadtree.py	(revision 7718)
@@ -27,6 +27,4 @@
 
 
-# FIXME(Ole): Could we come up with a less confusing structure?
-# FIXME(James): remove this global var
 LAST_TRIANGLE = [[-10, -10,
                    (num.array([[max_float, max_float],
@@ -57,5 +55,5 @@
         N = len(mesh)
         self.mesh = mesh
-        self.last_triangle = LAST_TRIANGLE	       
+        self.set_last_triangle()  
 
         # Get x,y coordinates for all vertices for all triangles
@@ -91,16 +89,15 @@
 
         """
+        
+        x = ensure_numeric(x, num.float)
+                 
         if self.last_triangle[0][1] != -10:
             # check the last triangle found first
             element_found, sigma0, sigma1, sigma2, k = \
                        self._search_triangles_of_vertices(self.last_triangle, x)
-
             if element_found:
                 return True, sigma0, sigma1, sigma2, k
 
         branch = self.last_triangle[0][1]
-        
-        if branch == -10:
-            branch = self   
 
         # test neighbouring tris
@@ -112,22 +109,9 @@
             return True, sigma0, sigma1, sigma2, k       
 
-        # search to bottom of tree from last found leaf    
-        tri_data = branch.search(x)
-        triangles = self._trilist_from_data(tri_data)            
-        element_found, sigma0, sigma1, sigma2, k = \
-                    self._search_triangles_of_vertices(triangles, x)
-        if element_found:
-            return True, sigma0, sigma1, sigma2, k
-
         # search rest of tree
         element_found = False
-        while branch:
-            if not branch.parent:
-                # search from top of tree if we are at root
-                siblings = [self]
-            else:
-                siblings = branch.get_siblings()
-                
-            for sibling in siblings:
+        next_search = [branch]
+        while branch:               
+            for sibling in next_search:
                 tri_data = sibling.search(x)
                 triangles = self._trilist_from_data(tri_data)            
@@ -136,5 +120,6 @@
                 if element_found:
                     return True, sigma0, sigma1, sigma2, k
-                                        
+            
+            next_search = branch.get_siblings()                            
             branch = branch.parent
             if branch:
@@ -144,5 +129,5 @@
                             self._search_triangles_of_vertices(triangles, x)
                 if element_found:
-                    return True, sigma0, sigma1, sigma2, k        
+                    return True, sigma0, sigma1, sigma2, k      
 
         return element_found, sigma0, sigma1, sigma2, k
@@ -158,15 +143,6 @@
         This function is responsible for most of the compute time in
         fit and interpolate.
-        """
-        x = ensure_numeric(x, num.float)     
-        
-        # 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
-        element_found = False    
+        """  
+
         for k, node, tri_verts_norms in triangles:
             tri = tri_verts_norms[0]
@@ -178,27 +154,18 @@
             # Input check disabled to speed things up.    
             if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)):
-                
                 n0, n1, n2 = tri_verts_norms[1]        
-                sigma0, sigma1, sigma2 =\
-                    compute_interpolation_values(tri, n0, n1, n2, x)
-                    
-                element_found = True    
+                xi0, xi1, xi2 = tri
+
+                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)              
 
                 # Don't look for any other triangles in the triangle list
                 self.last_triangle = [[k, node, tri_verts_norms]]
-                break            
-                
-        return element_found, sigma0, sigma1, sigma2, k
+                return True, sigma0, sigma1, sigma2, k                      
+        return False, -1, -1, -1, -10
 
 
     def _trilist_from_data(self, indices):
-        """return a list of lists. For the inner lists,
-        The first element is the triangle index,
-        the second element is a list.for this list
-           the first element is a list of three (x, y) vertices,
-           the following elements are the three triangle normals.
-
-        """
-
         ret_list = []
         for i, node in indices:
@@ -212,19 +179,7 @@
     def set_last_triangle(self):
         self.last_triangle = LAST_TRIANGLE
+        self.last_triangle[0][1] = self # point at root by default          
 
     
-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.
-    """
 
-    # Get the three vertex_points of candidate triangle k
-    xi0, xi1, xi2 = triangle
-
-    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
                 
Index: /anuga_core/source/anuga/geometry/quad.py
===================================================================
--- /anuga_core/source/anuga/geometry/quad.py	(revision 7717)
+++ /anuga_core/source/anuga/geometry/quad.py	(revision 7718)
@@ -19,5 +19,5 @@
     """
   
-    def __init__(self, extents, parent, depth = 0,
+    def __init__(self, extents, parent,
          name = 'cell'):
   
@@ -27,6 +27,4 @@
         self.extents = extents
         self.parent = parent
-        self.searched = [-10, -15]
-        self.depth = depth
         
         # The points in this cell     
@@ -36,6 +34,6 @@
     
     def __repr__(self):
-        str = '(%d)%s: leaves: %d' \
-               % (self.depth, self.name , len(self.leaves))    
+        str = '%s: leaves: %d' \
+               % (self.name , len(self.leaves))    
         if self.children:
             str += ', children: %d' % (len(self.children))
@@ -84,12 +82,24 @@
             # try splitting this cell and see if we get a trivial in
             subregion1, subregion2 = self.extents.split()
-            if subregion1.is_trivial_in(new_region):
-                self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)]    
-                self.children[0]._insert(new_leaf)
-                return
-            elif subregion2.is_trivial_in(new_region):
-                self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)]    
-                self.children[1]._insert(new_leaf)
-                return                
+            
+            # option 1 - try splitting 4 ways
+            subregion11, subregion12 = subregion1.split()                                
+            subregion21, subregion22 = subregion2.split()
+            regions = [subregion11, subregion12, subregion21, subregion22]
+            for region in regions:
+                if region.is_trivial_in(new_region):
+                    self.children = [Cell(x, parent=self) for x in regions]
+                    self._insert(new_leaf)
+                    return               
+
+            # option 2 - try splitting 2 ways
+            #if subregion1.is_trivial_in(new_region):
+                #self.children = [Cell(subregion1, self), Cell(subregion2, self)]    
+                #self.children[0]._insert(new_leaf)
+                #return
+            #elif subregion2.is_trivial_in(new_region):
+                #self.children = [Cell(subregion1, self), Cell(subregion2, self)]    
+                #self.children[1]._insert(new_leaf)
+                #return                
     
         # recursion ended without finding a fit, so attach it as a leaf
@@ -128,14 +138,13 @@
         if depth == 0:
             log.critical() 
-        print '%s%s' % ('  '*depth, self.name), self.extents,' [', self.leaves, ']'
+        print '%s%s'  % ('  '*depth, self.name), self.extents,' [', self.leaves, ']'
         if self.children:
             log.critical()
             for child in self.children:
                 child.show(depth+1)
- 
 
     def search(self, x):
         """return a list of possible intersections with geometry"""
-        
+      
         intersecting_regions = self.test_leaves(x)
         
Index: /anuga_core/source/anuga/geometry/test_meshquad.py
===================================================================
--- /anuga_core/source/anuga/geometry/test_meshquad.py	(revision 7717)
+++ /anuga_core/source/anuga/geometry/test_meshquad.py	(revision 7718)
@@ -145,27 +145,38 @@
         self.assertEqual(idx, 2)
         
-    def NOtest_get_parent():
-        """ Get a child's parent.
+    def NOtest_num_visits(self):
+        """ Test optimisation code.
         """
-        cell = Cell(AABB(0, 6, 0, 6), 'cell')
+        a = [3, 7]
+        b = [5, 7]
+        c = [5, 5]
+        d = [7, 7]
+        e = [15, 15]
+        f = [15, 30]
+        g = [30, 10]
+        h = [30, 30]
 
-        p0 = [2,1]
-        p1 = [4,1]
-        p2 = [4,4]
-        p3 = [2,4]
-        p4 = [5,4]
-
-        points = [p0,p1,p2, p3, p4]
-        vertices = [[0,1,2],[0,2,3],[1,4,2]]
+        points = [a, b, c, d, e, f, g, h]
+        
+        #bac, bce, ecf, dbe, daf, dae
+        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
 
         mesh = Mesh(points, vertices)
 
-        Q = MeshQuadtree(mesh)
-        results = Q.search([4.5, 3])
-        assert len(results) == 1
-        self.assertEqual(results[0], 2)
-        results = Q.search([5,4.])
-        self.assertEqual(len(results),1)
-        self.assertEqual(results[0], 2)     
+        Q = MeshQuadtree(mesh)    
+
+
+        results = Q.search_fast([5.5, 5.5])
+        print 'visits: ', Q.count_visits()
+        
+        Q.clear_visits()
+        results = Q.search_fast([30, 10])
+        print 'visits: ', Q.count_visits()
+        
+        print 'second time:'
+
+        Q.clear_visits()        
+        results = Q.search_fast([5.5, 5.5])
+        print 'visits: ', Q.count_visits()
 ################################################################################
 
Index: /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7717)
+++ /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7718)
@@ -1812,7 +1812,7 @@
                 gauge_values[i].append(w[i])
 
-        #Reference (nautilus 26/6/2008)
-
-        G0 = [-0.20000000000000001, -0.20000000000000001, -0.20000000000000001, -0.19729755264559332, -0.18143006781428656, -0.17264603126739803, -0.16490767093739186, -0.15580226478343825, -0.15341783464373857, -0.14986976789238646, -0.1483294354194323, -0.19385112169384708, -0.19907249891161938, -0.19911053885072424, -0.19913592111915288, -0.1991524040791994, -0.19916445999802629, -0.19917790655716, -0.19919280664859376, -0.19920527965237472, -0.19921381524768964, -0.1992178338453264, -0.19921637002152079, -0.19920783637167225, -0.19921088191840772, -0.19919037392917549, -0.19918812396862254, -0.19918939056535054, -0.19919336800395021, -0.19919832209289323, -0.19920329732763048, -0.19920824490004829, -0.19921301385072837, -0.19921748460975705, -0.1992215655528298, -0.19922516461276893, -0.19922820250665138, -0.19923065047412888, -0.19923254166903823, -0.1992339741831877, -0.19923510970454786, -0.19923612584283332, -0.19923716297385921, -0.199238301846201, -0.19923956605537632, -0.19924092368339749, -0.19924231511786381, -0.19924367609907154, -0.19924495328100339, -0.19924611576332413, -0.19924716014245455]
+        # Captured data from code manually inspected for correctness 11/5/2010
+
+        G0 = [-0.20000000000000001, -0.20000000000000001, -0.19999999639817381, -0.19900000000000001, -0.19787180620524364, -0.18775357859677891, -0.17717419377267535, -0.17038758570131016, -0.1632648895200452, -0.16084615293651969, -0.15796677906889614, -0.15610187077034218, -0.19890877268012952, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19767101056475031, -0.19664020763833073, -0.19625551309686509, -0.19622800460530576, -0.19670420796840654, -0.19757556964453959, -0.19853453054821305, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001]
         G1 = [-0.29999999999999993, -0.29999999999999993, -0.29312012778165086, -0.26736324399674027, -0.23794021958953263, -0.21878933349240881, -0.2063216250007715, -0.19630760483269205, -0.18778429435995911, -0.18241377240055098, -0.17860448861463396, -0.17352695896444634, -0.16449192186501546, -0.18953645096562932, -0.20557835913635514, -0.21071182844586539, -0.21191011610076083, -0.21068538990761979, -0.20737339293005266, -0.20407505794124087, -0.20190406545920248, -0.19987295968833538, -0.19832486432402194, -0.19774281454937118, -0.19763203737080062, -0.19769585527535841, -0.19795493497088903, -0.19841633252351859, -0.19886909990019369, -0.19940203310127189, -0.19990812587783294, -0.20016182994265733, -0.20026318580583596, -0.20022855664242351, -0.20010463117095048, -0.199930824934102, -0.19973147245180795, -0.19957882377401115, -0.19948188719957555, -0.19944190810764481, -0.19945196574982568, -0.19949690579414533, -0.19955788251561279, -0.19961580710612042, -0.19965570581001513, -0.19966924072268322, -0.19966164565962002, -0.1996429044375452, -0.19962419292036049, -0.19961203640143024, -0.19960927744740292]
         G2 = [-0.40000000000000002, -0.39653752812451043, -0.32890250043784469, -0.29967800724628818, -0.27572700190320054, -0.25893877811152488, -0.2422572785756886, -0.22758771043596362, -0.21433402198450011, -0.20266876405969875, -0.19393588052745447, -0.18766846928826031, -0.18182722392707523, -0.16865402349317951, -0.18110962350013771, -0.19547215504627991, -0.20229670797148233, -0.20494160396870145, -0.20614855426298931, -0.20616358034144364, -0.2050749944506553, -0.20329491765401383, -0.20174840479413605, -0.20040155802388543, -0.19945041001985703, -0.19898209259951644, -0.19868944730468457, -0.19849574784901269, -0.19855116942355155, -0.19885365879108613, -0.19925381355621719, -0.19971339601933866, -0.20000203652933318, -0.20011337898513715, -0.2001188815225943, -0.20007116187977508, -0.20000646342088527, -0.19993034594397555, -0.19985512250322437, -0.19978530424723981, -0.19973800656091345, -0.19971854873923284, -0.19972327598147538, -0.19974189986453444, -0.19976101408429672, -0.19977705650936498, -0.19978787136068635, -0.19979206946039513, -0.19978865745584601, -0.19978181317853855, -0.19977875136946421]
