Index: /anuga_core/source/anuga/fit_interpolate/fit.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/fit.py	(revision 7709)
+++ /anuga_core/source/anuga/fit_interpolate/fit.py	(revision 7710)
@@ -57,6 +57,5 @@
                  mesh_origin=None,
                  alpha = None,
-                 verbose=False,
-                 max_vertices_per_cell=None):
+                 verbose=False):
 
 
@@ -79,7 +78,4 @@
               If specified vertex coordinates are assumed to be
               relative to their respective origins.
-
-          max_vertices_per_cell: Number of vertices in a quad tree cell
-          at which the cell is split into 4.
 
           Note: Don't supply a vertex coords as a geospatial object and
@@ -103,6 +99,5 @@
                  mesh,
                  mesh_origin,
-                 verbose,
-                 max_vertices_per_cell)
+                 verbose)
         
         m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
Index: /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 7709)
+++ /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 7710)
@@ -54,6 +54,5 @@
                  mesh=None,
                  mesh_origin=None,
-                 verbose=False,
-                 max_vertices_per_cell=None):
+                 verbose=False):
 
 
@@ -82,13 +81,8 @@
               relative to their respective origins.
 
-          max_vertices_per_cell: Number of vertices in a quad tree cell
-          at which the cell is split into 4.
-
           Note: Don't supply a vertex coords as a geospatial object and
               a mesh origin, since geospatial has its own mesh origin.
         """
         global build_quadtree_time
-        if max_vertices_per_cell == None:
-            max_vertices_per_cell = MAX_VERTICES_PER_CELL
         if mesh is None:
             if vertex_coordinates is not None and  triangles is not None:
@@ -113,6 +107,5 @@
             #This stores indices of vertices
             t0 = time.time()
-            self.root = build_quadtree(self.mesh,
-                                       max_points_per_cell = max_vertices_per_cell)
+            self.root = build_quadtree(self.mesh)
         
             build_quadtree_time =  time.time()-t0
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7709)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7710)
@@ -57,5 +57,4 @@
 # @param interpolation_points Points to interpolate to.
 # @param mesh_origin A geo_ref object or 3-tuples of UTMzone, easting, northing.
-# @param max_vertices_per_cell Max number of vertices before splitting cell.
 # @param start_blocking_len Block if # of points greater than this.
 # @param use_cache If True, cache.
@@ -66,5 +65,4 @@
                 interpolation_points,
                 mesh_origin=None,
-                max_vertices_per_cell=None,
                 start_blocking_len=500000,
                 use_cache=False,
@@ -103,7 +101,4 @@
                  relative to their respective origins.
 
-    max_vertices_per_cell: Number of vertices in a quad tree cell
-                           at which the cell is split into 4.
-
                            Note: Don't supply a vertex coords as a geospatial
                            object and a mesh origin, since geospatial has its
@@ -136,5 +131,4 @@
             ensure_numeric(triangles))
     kwargs = {'mesh_origin': mesh_origin,
-              'max_vertices_per_cell': max_vertices_per_cell,
               'verbose': verbose}
 
@@ -175,6 +169,5 @@
                  triangles,
                  mesh_origin=None,
-                 verbose=False,
-                 max_vertices_per_cell=None):
+                 verbose=False):
 
         """ Build interpolation matrix mapping from
@@ -209,6 +202,5 @@
                                 triangles=triangles,
                                 mesh_origin=mesh_origin,
-                                verbose=verbose,
-                                max_vertices_per_cell=max_vertices_per_cell)
+                                verbose=verbose)
 
         # Initialise variables
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7709)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 7710)
@@ -307,6 +307,5 @@
 
         data = [ [4,4] ]
-        interp = Interpolate(points, triangles,
-                               max_vertices_per_cell = 4)
+        interp = Interpolate(points, triangles)
         #print "PDSG - interp.get_A()", interp.get_A()
         answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
Index: /anuga_core/source/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7709)
+++ /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 7710)
@@ -50,5 +50,5 @@
         mesh.check_integrity()
 
-        root = build_quadtree(mesh, max_points_per_cell = 1)
+        root = build_quadtree(mesh)
         set_last_triangle()
 
@@ -58,6 +58,6 @@
         found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, [0, 0])
         assert found is True
-		
-		
+        
+        
     def test_small(self):
         """test_small: Two triangles
@@ -70,5 +70,5 @@
         mesh.check_integrity()
 
-        root = build_quadtree(mesh, max_points_per_cell = 1)
+        root = build_quadtree(mesh)
         set_last_triangle()
 
@@ -76,6 +76,6 @@
         found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
         assert k == 1 # Triangle one
-        assert found is True		
-		
+        assert found is True        
+        
     def test_bigger(self):
         """test_bigger
@@ -90,5 +90,5 @@
         mesh.check_integrity()
 
-        root = build_quadtree(mesh, max_points_per_cell = 4)
+        root = build_quadtree(mesh)
         set_last_triangle()
 
@@ -98,6 +98,6 @@
             
             found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,
-                                                           ensure_numeric(x))								   
-														   
+                                                           ensure_numeric(x))                                   
+                                                           
             if k >= 0:
                 V = mesh.get_vertex_coordinates(k) # nodes for triangle k
@@ -121,25 +121,25 @@
 
         
-        for m in range(8):
-            root = build_quadtree(mesh, max_points_per_cell = m)
-            set_last_triangle()
-            #print m, root.show()
-
-            for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.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)
-
-                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                
-
+
+        root = build_quadtree(mesh)
+        set_last_triangle()
+        #print m, root.show()
+
+        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
+                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
+                  [10, 3]]:
             
-                if m == k == 0: return    
+            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 k == 0: return    
 
     def test_underlying_function(self):
@@ -150,5 +150,5 @@
         mesh = Mesh(points, vertices, boundary)
 
-        root = build_quadtree(mesh, max_points_per_cell = 4)
+        root = build_quadtree(mesh)
         set_last_triangle()
 
@@ -157,5 +157,5 @@
 
         triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
-	
+    
         found, sigma0, sigma1, sigma2, k = \
                _search_triangles_of_vertices(triangles, x)
Index: /anuga_core/source/anuga/utilities/aabb.py
===================================================================
--- /anuga_core/source/anuga/utilities/aabb.py	(revision 7710)
+++ /anuga_core/source/anuga/utilities/aabb.py	(revision 7710)
@@ -0,0 +1,72 @@
+# Allow children to be slightly bigger than their parents to prevent straddling of a boundary
+SPLIT_BORDER_RATIO    = 0.55
+
+class AABB:
+    """Axially-aligned bounding box class.
+    """
+    
+    def __init__(self, xmin, xmax, ymin, ymax):
+        """ Define axially-algned bounding box.
+            xmin is minimum x
+            xmax is maximum x (absolute coord, ie, not size)
+            ymin is minimum y
+            ymax is maximum y (absolute coord, ie, not size)
+        """
+        self.xmin = xmin    
+        self.xmax = xmax
+        self.ymin = ymin    
+        self.ymax = ymax
+
+
+    def __repr__(self):
+        return 'AABB(xmin:%f, xmax:%f, ymin:%f, ymax:%f)' \
+               % (round(self.xmin,1), round(self.xmax,1), round(self.ymin,1), round(self.ymax, 1)) 
+
+
+    def grow(self, amount):
+        """ Expand region by given amount.
+            amount is a multiplier, ie 1.1 will expand border by 10%.
+        """
+        self.ymax += (self.ymax-self.ymin)*amount
+        self.xmax += (self.xmax-self.xmin)*amount
+        self.ymin -= (self.ymax-self.ymin)*amount
+        self.xmin -= (self.xmax-self.xmin)*amount    
+
+        
+    def size(self):
+        """return size as (w,h)"""
+        return self.xmax - self.xmin, self.ymax - self.ymin
+
+        
+    def split(self, border=SPLIT_BORDER_RATIO):
+        """Split along shorter axis.
+           return 2 subdivided AABBs.
+        """
+        
+        width, height = self.size()
+        assert width >= 0 and height >= 0
+        
+        if (width > height):
+            # split vertically
+            return AABB(self.xmin, self.xmin+width*border, self.ymin, self.ymax), \
+                   AABB(self.xmax-width*border, self.xmax, self.ymin, self.ymax)
+        else:
+            # split horizontally       
+            return AABB(self.xmin, self.xmax, self.ymin, self.ymin+height*border), \
+                   AABB(self.xmin, self.xmax, self.ymax-height*border, self.ymax)    
+
+    
+    def is_trivial_in(self, test):
+        """ Is trivial in.
+            test an x,y point to test against the bounding box
+            return True if the test point falls within the bounding box
+        """
+        if (test.xmin < self.xmin) or (test.xmax > self.xmax):
+            return False        
+        if (test.ymin < self.ymin) or (test.ymax > self.ymax):
+            return False        
+        return True
+ 
+    def contains(self, x, y):
+        return (self.xmin <= x <= self.xmax) and (self.ymin <= y <= self.ymax)
+        
Index: /anuga_core/source/anuga/utilities/quad.py
===================================================================
--- /anuga_core/source/anuga/utilities/quad.py	(revision 7709)
+++ /anuga_core/source/anuga/utilities/quad.py	(revision 7710)
@@ -9,52 +9,6 @@
 import string, types, sys
 import anuga.utilities.log as log
+from anuga.utilities.aabb import AABB
 
-# Allow children to be slightly bigger than their parents to prevent straddling of a boundary
-SPLIT_BORDER_RATIO    = 0.55
-
-class AABB:
-    """Axially-aligned bounding box class.
-    """
-    
-    def __init__(self, xmin, xmax, ymin, ymax):
-        self.xmin = xmin    
-        self.xmax = xmax
-        self.ymin = ymin    
-        self.ymax = ymax
-
-    def __repr__(self):
-        return '(xmin:%f, xmax:%f, ymin:%f, ymax:%f)' \
-               % (round(self.xmin,1), round(self.xmax,1), round(self.ymin,1), round(self.ymax, 1)) 
-        
-    def size(self):
-        """return size as (w,h)"""
-        return self.xmax - self.xmin, self.ymax - self.ymin
-        
-    def split(self, border=SPLIT_BORDER_RATIO):
-        """Split along shorter axis.
-           return 2 subdivided AABBs.
-        """
-        
-        width, height = self.size()
-        assert width >= 0 and height >= 0
-        
-        if (width > height):
-            # split vertically
-            return AABB(self.xmin, self.xmin+width*border, self.ymin, self.ymax), \
-                   AABB(self.xmax-width*border, self.xmax, self.ymin, self.ymax)
-        else:
-            # split horizontally       
-            return AABB(self.xmin, self.xmax, self.ymin, self.ymin+height*border), \
-                   AABB(self.xmin, self.xmax, self.ymax-height*border, self.ymax)    
-    
-    def is_trivial_in(self, test):
-        if (test.xmin < self.xmin) or (test.xmax > self.xmax):
-            return False        
-        if (test.ymin < self.ymin) or (test.ymax > self.ymax):
-            return False        
-        return True
- 
-    def contains(self, x, y):
-        return (self.xmin <= x <= self.xmax) and (self.ymin <= y <= self.ymax)
             
 class Cell(TreeNode):
@@ -63,12 +17,4 @@
     One cell in the plane delimited by southern, northern,
     western, eastern boundaries.
-
-    Public Methods:
-        insert(point)
-        search(x, y)
-        split()
-        store()
-        retrieve()
-        count()
     """
   
@@ -88,5 +34,5 @@
     def __repr__(self):
         str = '%s: leaves: %d' \
-               % (self.name , len(self.leaves))	
+               % (self.name , len(self.leaves))    
         if self.children:
             str += ', children: %d' % (len(self.children))
@@ -211,35 +157,13 @@
 #from anuga.pmesh.mesh import Mesh
     
-def build_quadtree(mesh, max_points_per_cell = 4):
+def build_quadtree(mesh):
     """Build quad tree for mesh.
 
-    All vertices in mesh are stored in quadtree and a reference
+    All vertex indices in mesh are stored in a quadtree and a reference
     to the root is returned.
     """
-
-
-    #Make root cell
-    #print mesh.coordinates
-
-    xmin, xmax, ymin, ymax = mesh.get_extent(absolute=True)
-    
-    # Ensure boundary points are fully contained in region
-    # It is a property of the cell structure that
-    # points on xmax or ymax of any given cell
-    # belong to the neighbouring cell.
-    # Hence, the root cell needs to be expanded slightly
-    ymax += (ymax-ymin)/10
-    xmax += (xmax-xmin)/10
-
-    # To avoid round off error
-    ymin -= (ymax-ymin)/10
-    xmin -= (xmax-xmin)/10   
-
-    #print "xmin", xmin 
-    #print "xmax", xmax
-    #print "ymin", ymin 
-    #print "ymax", ymax
-    
-    root = Cell(AABB(xmin, xmax, ymin, ymax))
+    extents = AABB(*mesh.get_extent(absolute=True))   
+    extents.grow(1.001) # To avoid round off error
+    root = Cell(extents)
     
     N = len(mesh)
@@ -247,5 +171,5 @@
     # Get x,y coordinates for all vertices for all triangles
     V = mesh.get_vertex_coordinates(absolute=True)
-	
+    
     # Check each triangle
     for i in range(N):
