Index: /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 7712)
+++ /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 7713)
@@ -36,8 +36,4 @@
 import numpy as num
 
-
-# tests fail if 2 is used
-MAX_VERTICES_PER_CELL = 13 # A value of 8 or lower can cause problems,
-                           # if a vert has 9 triangles.
                            
 build_quadtree_time = 0
Index: /anuga_core/source/anuga/fit_interpolate/search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 7712)
+++ /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 7713)
@@ -33,5 +33,5 @@
 # FIXME(Ole): Could we come up with a less confusing structure?
 # FIXME(James): remove this global var
-LAST_TRIANGLE = [[-10,
+LAST_TRIANGLE = [[-10, -10,
                    (num.array([[max_float, max_float],
                                [max_float, max_float],
@@ -64,35 +64,63 @@
     global search_more_cells_time
 
-    # Search the last triangle first
+    if last_triangle[0][1] != -10:
+        # check the last triangle found first
+        element_found, sigma0, sigma1, sigma2, k = \
+                   _search_triangles_of_vertices(last_triangle, x)
+
+        if element_found:
+            return True, sigma0, sigma1, sigma2, k
+
+    # search to bottom of tree from last found leaf
+    branch = last_triangle[0][1]
+    
+    if branch == -10:
+        branch = root   
+    
+    tri_data = branch.search(x[0], x[1])
+    triangles = _trilist_from_data(mesh, tri_data)            
     element_found, sigma0, sigma1, sigma2, k = \
-        _search_triangles_of_vertices(last_triangle, x)
-                   
-    if element_found is True:
-        return element_found, sigma0, sigma1, sigma2, k
-
-    
-    # Get triangles in the cell that the point is in.
-    tri_indices = root.search(x[0], x[1])
-    triangles = _trilist_from_indices(mesh, tri_indices)
-    
-    element_found, sigma0, sigma1, sigma2, k = \
-                   _search_triangles_of_vertices(triangles, x)
-
-    is_more_elements = True
-    
-    # while not element_found and is_more_elements:
-        # triangles, branch = root.expand_search()
-        # if branch == []:
-            # # Searching all the verts from the root cell that haven't
-            # # been searched.  This is the last try
-            # element_found, sigma0, sigma1, sigma2, k = \
-                           # _search_triangles_of_vertices(triangles, x)
-            # is_more_elements = False
-        # else:
-            # element_found, sigma0, sigma1, sigma2, k = \
-                       # _search_triangles_of_vertices(triangles, x)
-                       
-        
+                _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 = [root]
+        else:
+            siblings = branch.get_siblings()
+            
+        for sibling in siblings:
+            tri_data = sibling.search(x[0], x[1])
+            triangles = _trilist_from_data(mesh, tri_data)            
+            element_found, sigma0, sigma1, sigma2, k = \
+                        _search_triangles_of_vertices(triangles, x)
+            if element_found:
+                return True, sigma0, sigma1, sigma2, k
+                                    
+        branch = branch.parent
+        if branch:
+            tri_data = branch.test_leaves(x[0], x[1])
+            triangles = _trilist_from_data(mesh, tri_data)            
+            element_found, sigma0, sigma1, sigma2, k = \
+                        _search_triangles_of_vertices(triangles, x)
+            if element_found:
+                return True, sigma0, sigma1, sigma2, k        
+
     return element_found, sigma0, sigma1, sigma2, k
+
+
+    ## Get triangles in the cell that the point is in.
+    #tri_data = root.search(x[0], x[1])
+    #triangles = _trilist_from_data(mesh, tri_data)
+    
+    #element_found, sigma0, sigma1, sigma2, k = \
+                   #_search_triangles_of_vertices(triangles, x)
+
+    #if element_found:
+        #return True, sigma0, sigma1, sigma2, k
 
 
@@ -116,7 +144,8 @@
     sigma1 = -10.0
     k = -10
+    
     # For all vertices in same cell as point x
     element_found = False    
-    for k, tri_verts_norms in triangles:
+    for k, node, tri_verts_norms in triangles:
         tri = tri_verts_norms[0]
         tri = ensure_numeric(tri)        
@@ -135,5 +164,5 @@
 
             # Don't look for any other triangles in the triangle list
-            last_triangle = [[k, tri_verts_norms]]
+            last_triangle = [[k, node, tri_verts_norms]]
             break            
             
@@ -141,5 +170,5 @@
 
 
-def _trilist_from_indices(mesh, indices):
+def _trilist_from_data(mesh, indices):
     """return a list of lists. For the inner lists,
     The first element is the triangle index,
@@ -151,10 +180,10 @@
 
     ret_list = []
-    for i in indices:
+    for i, node in indices:
         vertices = mesh.get_vertex_coordinates(triangle_id=i, absolute=True)
         n0 = mesh.get_normal(i, 0)
         n1 = mesh.get_normal(i, 1)
         n2 = mesh.get_normal(i, 2) 
-        ret_list.append([i, [vertices, (n0, n1, n2)]])
+        ret_list.append([i, node, [vertices, (n0, n1, n2)]])
     return ret_list
                 
@@ -196,4 +225,6 @@
 class MeshQuadtree(Cell):
     """ A quadtree constructed from the given mesh.
+        This class is actually the root node of a tree,
+        and derives from a Cell.
     """
     def __init__(self, mesh):
@@ -205,5 +236,5 @@
         extents = AABB(*mesh.get_extent(absolute=True))   
         extents.grow(1.001) # To avoid round off error
-        Cell.__init__(self, extents)
+        Cell.__init__(self, extents, None)  # root has no parent
         
         N = len(mesh)
Index: /anuga_core/source/anuga/fit_interpolate/test_meshquad.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_meshquad.py	(revision 7712)
+++ /anuga_core/source/anuga/fit_interpolate/test_meshquad.py	(revision 7713)
@@ -43,8 +43,8 @@
 
         # test a point that falls within a triangle
-        result = Q.search(10, 10, get_vertices=True)
+        result = Q.search(10, 10)
         assert type(result) in [types.ListType,types.TupleType],\
                             'should be a list'
-        pos, index = result[0]
+        index, parent = result[0]
         self.assertEqual(index, 1)
 
@@ -138,8 +138,10 @@
         results = Q.search(4.5, 3)
         assert len(results) == 1
-        self.assertEqual(results[0], 2)
+        idx, _ = results[0]
+        self.assertEqual(idx, 2)
         results = Q.search(5,4.)
         self.assertEqual(len(results),1)
-        self.assertEqual(results[0], 2)
+        idx, _ = results[0]
+        self.assertEqual(idx, 2)
         
     def NOtest_get_parent():
Index: /anuga_core/source/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7712)
+++ /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7713)
@@ -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 _trilist_from_indices
+from search_functions import _trilist_from_data
 from search_functions import compute_interpolation_values, MeshQuadtree
 
@@ -156,5 +156,5 @@
         x = ensure_numeric([0.5, 0.5])
 
-        triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
+        triangles = _trilist_from_data(mesh, root.search(x[0], x[1]))
     
         found, sigma0, sigma1, sigma2, k = \
@@ -175,5 +175,5 @@
                   [10, 3]]:
                 
-            triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
+            triangles = _trilist_from_data(mesh, root.search(x[0], x[1]))
 
             #print x, candidate_vertices
Index: /anuga_core/source/anuga/geometry/quad.py
===================================================================
--- /anuga_core/source/anuga/geometry/quad.py	(revision 7712)
+++ /anuga_core/source/anuga/geometry/quad.py	(revision 7713)
@@ -19,5 +19,5 @@
     """
   
-    def __init__(self, extents,
+    def __init__(self, extents, parent, depth = 0,
          name = 'cell'):
   
@@ -26,4 +26,7 @@
     
         self.extents = extents
+        self.parent = parent
+        self.searched = [-10, -15]
+        self.depth = depth
         
         # The points in this cell     
@@ -33,6 +36,6 @@
     
     def __repr__(self):
-        str = '%s: leaves: %d' \
-               % (self.name , len(self.leaves))    
+        str = '(%d)%s: leaves: %d' \
+               % (self.depth, self.name , len(self.leaves))    
         if self.children:
             str += ', children: %d' % (len(self.children))
@@ -82,9 +85,9 @@
             subregion1, subregion2 = self.extents.split()
             if subregion1.is_trivial_in(new_region):
-                self.children = [Cell(subregion1), Cell(subregion2)]    
+                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), Cell(subregion2)]    
+                self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)]    
                 self.children[1]._insert(new_leaf)
                 return                
@@ -132,7 +135,23 @@
  
 
-    def search(self, x, y, get_vertices = False):
+    def search(self, x, y):
         """return a list of possible intersections with geometry"""
+
+ #       if self.searched == [x,y]:
+ #           print 'ERROR: already searched at ', [x,y]
+ #       self.searched = [x,y]
         
+        intersecting_regions = self.test_leaves(x, y)
+        
+        # recurse down into nodes that the point passes through
+        if self.children:
+            for child in self.children:    
+                if child.extents.contains(x, y):
+                    intersecting_regions.extend(child.search(x, y))
+             
+        return intersecting_regions
+ 
+ 
+    def test_leaves(self, x, y):
         intersecting_regions = []
         
@@ -141,15 +160,17 @@
             aabb, data = leaf
             if aabb.contains(x, y):
-                if get_vertices:
-                    intersecting_regions.append(leaf)
-                else:
-                    intersecting_regions.append(data)
-        
-        # recurse down into nodes that the point passes through
-        if self.children:
-            for child in self.children:    
-                if child.extents.contains(x, y):
-                    intersecting_regions.extend(child.search(x, y, get_vertices))
-             
-        return intersecting_regions
-        
+                intersecting_regions.append((data, self))
+                
+        return intersecting_regions                
+ 
+ 
+    def get_siblings(self):
+        """ return parent and siblings of this node.
+        """
+        if not self.parent:
+            return []
+         
+        siblings = list(self.parent.children)
+        siblings.remove(self)
+        return siblings
+                
Index: /anuga_core/source/anuga/geometry/test_geometry.py
===================================================================
--- /anuga_core/source/anuga/geometry/test_geometry.py	(revision 7712)
+++ /anuga_core/source/anuga/geometry/test_geometry.py	(revision 7713)
@@ -59,5 +59,5 @@
         
     def test_add_data(self):
-        cell = Cell(AABB(0,10, 0,5))
+        cell = Cell(AABB(0,10, 0,5), None)
         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
@@ -70,17 +70,27 @@
         
     def test_search(self):
-        test_region = (AABB(8,9, 1, 2), 222)
-        cell = Cell(AABB(0,10, 0,5))
-        cell.insert([(AABB(1,3, 1, 3), 111), test_region,  \
+        test_tag = 222
+        cell = Cell(AABB(0,10, 0,5), None)
+        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), test_tag),  \
                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
 
-        result =  cell.search(x = 8.5, y = 1.5, get_vertices=True)
+        result = cell.search(x = 8.5, y = 1.5)
         assert type(result) in [types.ListType,types.TupleType],\
                             'should be a list'
-        self.assertEqual(result, [test_region], 'only 1 point should intersect')
+        assert(len(result) == 1)
+        data, node = result[0]
+        self.assertEqual(data, test_tag, 'only 1 point should intersect')
 
+    def test_get_siblings(self):
+        """ Make sure children know their parent. """
+        cell = Cell(AABB(0,10, 0,5), None)
+        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222)])
+        assert len(cell.children) == 2
+        assert cell.parent == None
+        assert cell.children[0].parent == cell
+        assert cell.children[1].parent == cell
 
     def test_clear_1(self):
-        cell = Cell(AABB(0,10, 0,5))    
+        cell = Cell(AABB(0,10, 0,5), None)    
         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
@@ -94,5 +104,5 @@
 
 if __name__ == "__main__":
-    mysuite = unittest.makeSuite(Test_Quad,'test')
+    mysuite = unittest.makeSuite(Test_Geometry, 'test')
     runner = unittest.TextTestRunner()
     runner.run(mysuite)
