Index: /anuga_core/source/anuga/fit_interpolate/test_search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 4665)
+++ /anuga_core/source/anuga/fit_interpolate/test_search_functions.py	(revision 4666)
@@ -182,11 +182,9 @@
         mesh = Mesh(points, vertices)
 
-        # Don't do this, want to control the mx and mins
+        # Don't do this, want to control the max and mins
         #root = build_quadtree(mesh, max_points_per_cell=4)
     
-        #Initialise
-        Cell.initialise(mesh)
-
-        root = Cell(-3, 9, -3, 9,
+
+        root = Cell(-3, 9, -3, 9, mesh,
                     max_points_per_cell = 4)
         #Insert indices of all vertices
Index: /anuga_core/source/anuga/shallow_water/benchmark_sww2dem.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/benchmark_sww2dem.py	(revision 4665)
+++ /anuga_core/source/anuga/shallow_water/benchmark_sww2dem.py	(revision 4666)
@@ -22,4 +22,6 @@
 import profile , pstats
 import tempfile
+import gc
+
 from Scientific.IO.NetCDF import NetCDFFile
 from Numeric import array
@@ -27,5 +29,4 @@
 from anuga.fit_interpolate.interpolate import Interpolate
 from anuga.fit_interpolate.fit import Fit
-from anuga.pmesh.mesh import Mesh
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 from anuga.shallow_water import Domain
@@ -56,92 +57,81 @@
 
 
-class Benchmarksww2dem:
 
+def trial(vert_rows,
+          vert_columns,
+          output_quantities=['elevation'],
+          save=False,
+          verbose=False,
+          run_profile=False):
+    import time
+  
+    sww_file_name = build_sww(vert_rows, vert_columns, save=save)
+        
+    #Initial time and memory
+    t0 = time.time()
+    #m0 = None on windows
+    m0 = mem_usage()
+    #print "m0", m0
+
+    fileout = export_grid(sww_file_name, quantities=output_quantities)
+        
+    time_taken_sec = (time.time()-t0)
+    m1 = mem_usage()
+    #print "m1", m1
+    if m0 is None or m1 is None:
+        memory_used = None
+    else:
+        memory_used = (m1 - m0)
+    return time_taken_sec, memory_used, m0, m1
+
+def build_sww(vert_rows, vert_columns, save):
     """
-
-    Note(DSG-DSG): If you are interested in benchmarking fitting, before
-    and after blocking O:\1\dgray\before_blocking_subsandpit is before blocking
-
+    Build an sww file we can read. 
+    verts_rows and vert_columns have to be greater that 1
     """
-
-    def __init__(self):
-        pass
-
-    def trial(self,
-              vert_rows,
-              vert_columns,
-              output_quantities=['elevation'],
-              save=False,
-              verbose=False,
-              run_profile=False):
-  
-        import time
+    from anuga.pmesh.mesh import Mesh
+    m = Mesh()
+    m.build_grid(vert_rows, vert_columns)
+    if save is True:
+            m.export_mesh_file("aaaa.tsh")
+    mesh_dict =  m.Mesh2IOTriangulationDict()
+    
+    sww_fileName = tempfile.mktemp(".sww" )
+    # sww_fileName = "aa.sww" 
+    elevation = array(range(len(mesh_dict["vertices"])))
+    stage = elevation
+    ymomentum = elevation
+    xmomentum = elevation
         
-        sww_file_name = self._build_sww(vert_rows, vert_columns, save=save)
+    # NetCDF file definition
+    fid = NetCDFFile(sww_fileName, 'w')
+    sww = Write_sww()
+    sww.header(fid, 0,
+               len(mesh_dict['triangles']),
+               len(mesh_dict["vertices"]))
+    sww.triangulation(fid,
+                      mesh_dict["vertices"], mesh_dict['triangles'],
+                      elevation)
         
-        #Initial time and memory
-        t0 = time.time()
-        #m0 = None on windows
-        m0 = mem_usage()
-        #print "m0", m0
-
-        fileout = export_grid(sww_file_name, quantities=output_quantities)
-        
-        time_taken_sec = (time.time()-t0)
-        m1 = mem_usage()
-        #print "m1", m1
-        if m0 is None or m1 is None:
-            memory_used = None
-        else:
-            memory_used = (m1 - m0)
-        #print 'That took %.2f seconds' %time_taken_sec
-        return time_taken_sec, memory_used, m0, m1
-
-    def _build_sww(self, vert_rows, vert_columns, save):
-        """
-        Build an sww file we can read. 
-        verts_rows and vert_columns have to be greater that 1
-        """
-        m = Mesh()
-        m.build_grid(vert_rows, vert_columns)
-        if save is True:
-            m.export_mesh_file("aaaa.tsh")
-        mesh_dict =  m.Mesh2IOTriangulationDict()
-        
-        sww_fileName = tempfile.mktemp(".sww" )
-        # sww_fileName = "aa.sww" 
-        elevation = array(range(len(mesh_dict["vertices"])))
-        stage = elevation
-        ymomentum = elevation
-        xmomentum = elevation
-        
-        # NetCDF file definition
-        fid = NetCDFFile(sww_fileName, 'w')
-        sww = Write_sww()
-        sww.header(fid, 0,
-                   len(mesh_dict['triangles']),
-                   len(mesh_dict["vertices"]))
-        sww.triangulation(fid,
-                   mesh_dict["vertices"], mesh_dict['triangles'],
-                          elevation)
-        
-        for time_step in range(10):
-            sww.quantities(fid, 
-                           time=time_step,
-                           stage=stage,
-                           xmomentum=xmomentum,
-                           ymomentum=ymomentum)     
-        #Close
-        fid.close()
-        return sww_fileName
+    for time_step in range(10):
+        sww.quantities(fid, 
+                       time=time_step,
+                       stage=stage,
+                       xmomentum=xmomentum,
+                       ymomentum=ymomentum)     
+    #Close
+    fid.close()
+    del mesh_dict
+    del Mesh
+    return sww_fileName
     
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    b = Benchmarksww2dem() 
     delimiter = ','
 
     ofile = 'lbm_resultsII.csv'
     run_profile = False #True
-    size_list = [[4,5],[50,40]]
+    size_list = [[4,5],[5,40]]
+    #size_list = [[5,4]]
 
     quantitiy_list = [['elevation'],
@@ -165,5 +155,7 @@
             #sys.gettotalrefcount()
             #sys.getobjects() 
-            time, mem, m0, m1 = b.trial(size[0]
+            gc.collect()
+            #print " gc.get_objects() ", gc.get_objects() 
+            time, mem, m0, m1 = trial(size[0]
                                         ,size[1]
                                         ,output_quantities=quantitiy
@@ -172,4 +164,7 @@
             print "time",time
             print "mem", mem
+            gc.collect()
+            print "mem", mem
+            print " gc.get_objects() ", gc.get_objects() 
             #sys.gettotalrefcount() 
             fd.write(str(size[0]*size[1]) + delimiter +
Index: /anuga_core/source/anuga/utilities/quad.py
===================================================================
--- /anuga_core/source/anuga/utilities/quad.py	(revision 4665)
+++ /anuga_core/source/anuga/utilities/quad.py	(revision 4666)
@@ -27,5 +27,5 @@
     """
   
-    def __init__(self, southern, northern, western, eastern, 
+    def __init__(self, southern, northern, western, eastern, mesh,
 		 name = 'cell',
     	         max_points_per_cell = 4):
@@ -39,4 +39,5 @@
         self.western = round(western,5)    
         self.eastern = round(eastern,5)
+        self.mesh = mesh
 
         # The points in this cell     
@@ -63,11 +64,12 @@
         cn = self.northern
         cw = self.western    
-        ce = self.eastern
+        ce = self.eastern   
+        mesh = self.mesh
 
         # create 4 child cells
-        self.AddChild(Cell((cn+cs)/2,cn,cw,(cw+ce)/2,self.name+'_nw'))
-        self.AddChild(Cell((cn+cs)/2,cn,(cw+ce)/2,ce,self.name+'_ne'))
-        self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,self.name+'_se'))
-        self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,self.name+'_sw'))
+        self.AddChild(Cell((cn+cs)/2,cn,cw,(cw+ce)/2,mesh,self.name+'_nw'))
+        self.AddChild(Cell((cn+cs)/2,cn,(cw+ce)/2,ce,mesh,self.name+'_ne'))
+        self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,mesh,self.name+'_se'))
+        self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,mesh,self.name+'_sw'))
         
  
@@ -140,5 +142,5 @@
 	if len(args) == 2:
 	    point_id = int(args[1])
-            x, y = self.__class__.mesh.get_nodes()[point_id]
+            x, y = self.mesh.get_nodes()[point_id]
 
             #print point_id, x, y
@@ -181,6 +183,6 @@
                     self.store(point)
                 else:
-                    x = self.__class__.mesh.coordinates[point][0]
-                    y = self.__class__.mesh.coordinates[point][1]
+                    x = self.mesh.coordinates[point][0]
+                    y = self.mesh.coordinates[point][1]
                     print "(" + str(x) + "," + str(y) + ")"
                     raise 'point not in region: %s' %str(point)
@@ -223,12 +225,12 @@
             # print "verts", verts
             for vert in verts:
-                triangle_list = self.__class__.mesh.get_triangles_and_vertices_per_node(vert)
+                triangle_list = self.mesh.get_triangles_and_vertices_per_node(vert)
                 for k, _ in triangle_list:
                     if not triangles.has_key(k):
                         # print 'k',k
-                        tri = self.__class__.mesh.get_vertex_coordinates(k)
-                        n0 = self.__class__.mesh.get_normal(k, 0)
-                        n1 = self.__class__.mesh.get_normal(k, 1)
-                        n2 = self.__class__.mesh.get_normal(k, 2) 
+                        tri = self.mesh.get_vertex_coordinates(k)
+                        n0 = self.mesh.get_normal(k, 0)
+                        n1 = self.mesh.get_normal(k, 1)
+                        n2 = self.mesh.get_normal(k, 2) 
                         triangles[k]=(tri, (n0, n1, n2))
             self.triangles = triangles.items()
@@ -418,20 +420,21 @@
 	
 
-    #Class initialisation method	
-    def initialise(cls, mesh):
-        cls.mesh = mesh
-
-    initialise = classmethod(initialise)
+    #Class initialisation method
+    # this is bad.  It adds a huge memory structure to the class.
+    # When the instance is deleted the mesh hangs round (leaks).
+    #def initialise(cls, mesh):
+    #    cls.mesh = mesh
+
+    #initialise = classmethod(initialise)
 
 def build_quadtree(mesh, max_points_per_cell = 4):
     """Build quad tree for mesh.
 
-    All vertices in mesh are stored in quadtree and a reference to the root is returned.
+    All vertices in mesh are stored in quadtree and a reference
+    to the root is returned.
     """
 
     from Numeric import minimum, maximum
 
-    #Initialise
-    Cell.initialise(mesh)
 
     #Make root cell
@@ -445,8 +448,9 @@
 
     
-    #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
+    # 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
@@ -462,5 +466,5 @@
     
     #FIXME: Use mesh.filename if it exists
-    root = Cell(ymin, ymax, xmin, xmax,
+    root = Cell(ymin, ymax, xmin, xmax,mesh,
                 #name = .... 
                 max_points_per_cell = max_points_per_cell)
Index: /anuga_core/source/anuga/utilities/test_quad.py
===================================================================
--- /anuga_core/source/anuga/utilities/test_quad.py	(revision 4665)
+++ /anuga_core/source/anuga/utilities/test_quad.py	(revision 4666)
@@ -12,5 +12,4 @@
 
     def setUp(self):
-        self.cell = Cell(100, 140, 0, 40, 'cell')
 
         a = [3, 107]
@@ -30,5 +29,5 @@
 	mesh = Mesh(points, vertices)
         self.mesh = mesh
-        Cell.initialise(mesh)
+        self.cell = Cell(100, 140, 0, 40, mesh, 'cell')
 
     def tearDown(self):
