Index: /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py	(revision 4650)
+++ /anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py	(revision 4651)
@@ -28,5 +28,7 @@
 from anuga.pmesh.mesh import Mesh
 from anuga.geospatial_data.geospatial_data import Geospatial_data
-#from anuga.shallow_water import Domain
+from anuga.shallow_water import Domain
+from anuga.fit_interpolate.fit import Fit, fit_to_mesh
+from anuga.fit_interpolate.interpolate import benchmark_interpolate
 
 def mem_usage():
@@ -96,5 +98,8 @@
                        "PPC" + str(max_points_per_cell) + \
                        ".txt"
-                       
+                    
+        
+        domain = Domain(mesh_dict['vertices'], mesh_dict['triangles'],
+                        use_cache=False, verbose=False)
         #Initial time and memory
         t0 = time.time()
@@ -102,25 +107,21 @@
         m0 = mem_usage()
         
-        #domain = Domain(mesh_dict['vertices'], mesh_dict['triangles'],
-         #               use_cache=False, verbose=False)
         if is_fit is True:
-            from anuga.fit_interpolate.fit import Fit, fit_to_mesh
 
             print "Fit in Fit"
+            geospatial = Geospatial_data(points_dict['points'],
+                                     points_dict['point_attributes'])
             if use_file_type == None:
-                points = points_dict['points']
-                point_attributes = points_dict['point_attributes']
+                points = geospatial
+                filename = None
             else:
                 #check that the type
                 fileName = tempfile.mktemp("." + use_file_type)
-                G1 = Geospatial_data(points_dict['points'],
-                                     points_dict['point_attributes'])
-                G1.export_points_file(fileName, absolute=True)
-                points = fileName
-                point_attributes = None
-                
+                geospatial.export_points_file(fileName, absolute=True)
+                points = None
+                filename = fileName
             if run_profile is True:
                     
-                s = """fit_to_mesh(mesh_dict['vertices'], mesh_dict['triangles'],points,point_attributes)"""
+                s = """domain.set_quantity('elevation',points,filename=filename,use_cache=False)"""
                 pobject = profile.Profile()
                 presult = pobject.runctx(s,
@@ -140,8 +141,6 @@
                 os.remove(prof_file)
             else:
-                fit_to_mesh(mesh_dict['vertices'],
-                            mesh_dict['triangles'],
-                            points, # this can also be a points file name
-                            point_attributes)
+                domain.set_quantity('elevation',points,filename=filename,
+                                    use_cache=False)
             if not use_file_type == None:
                 os.remove(fileName)
@@ -151,11 +150,6 @@
             print "Interpolate!"
             
-            interp = Interpolate(mesh_dict['vertices'],
-                                 mesh_dict['triangles'], 
-                                 max_vertices_per_cell = max_points_per_cell)
-            
             if run_profile:
-                s="""calc=interp.interpolate(mesh_dict['vertex_attributes']
-                ,points_dict['points'],start_blocking_len=blocking_len)"""
+                s="""benchmark_interpolate(mesh_dict['vertices'],mesh_dict['vertex_attributes'],mesh_dict['triangles'],points_dict['points'],max_points_per_cell=max_points_per_cell)"""
                 pobject = profile.Profile()
                 presult = pobject.runctx(s,
@@ -176,8 +170,9 @@
                     
             else:
-                calc = interp.interpolate(mesh_dict['vertex_attributes']
-                                          ,points_dict['points']
-                                          ,start_blocking_len = 500000)
-            
+                 benchmark_interpolate(mesh_dict['vertices'],
+                                       mesh_dict['vertex_attributes'],
+                                       mesh_dict['triangles'],
+                                       points_dict['points'],
+                                       max_points_per_cell=max_points_per_cell)
         time_taken_sec = (time.time()-t0)
         m1 = mem_usage()
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 4650)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 4651)
@@ -291,4 +291,24 @@
         return A
 
+def benchmark_interpolate(vertices,
+                          vertex_attributes,
+                          triangles, points,
+                          max_points_per_cell=None,
+                          start_blocking_len=500000):
+    """
+    No test for this yet.
+    Note, this has no time the input data has no time dimension.  Which is
+    different from most of the data we interpolate, eg sww info.
+     
+	Output:
+	  Interpolated values at inputted points.
+    """
+    interp = Interpolate(vertices,
+                         triangles, 
+                         max_vertices_per_cell=max_points_per_cell)
+            
+    calc = interp.interpolate(vertex_attributes
+                              ,points
+                              ,start_blocking_len=start_blocking_len)
 def interpolate_sww2csv(sww_file,
                         points,
Index: /anuga_core/source/anuga/fit_interpolate/search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 4650)
+++ /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 4651)
@@ -126,4 +126,9 @@
 
     # 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
+    
     xmax = max(xi0[0], xi1[0], xi2[0])
     xmin = min(xi0[0], xi1[0], xi2[0])
@@ -131,7 +136,4 @@
     ymin = min(xi0[1], xi1[1], xi2[1])
 
-    # Integrity check - machine precision is too hard
-    # so we use hardwired single precision 
-    epsilon = 1.0e-6
     
     if  x[0] > xmax + epsilon:
@@ -143,5 +145,4 @@
     if  x[1] < ymin - epsilon:
         return False,0,0,0,0
-    
     
     # Get the three normals 
Index: /anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py	(revision 4650)
+++ /anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py	(revision 4651)
@@ -14,12 +14,11 @@
 ofile = 'lbm_resultsII.csv'
 delimiter = ','
-run_profile = True
-use_least_squares_list = [False]
+run_profile = True #False #True
 is_fit_list = [True, False]
 num_of_points_list = [3, 200, 600, 2000, 6000, 10000, 20000] 
 maxArea_list = [ 0.008, 0.0016, 0.0008]
 max_points_per_cell_list = [2,4,8,16,30,64]
-use_file_type_list = ['pts']
-num_of_points_list = [10]
+use_file_type_list = ['pts'] #'pts'
+#num_of_points_list = [10]
 maxArea_list = [ 0.008]
 max_points_per_cell_list = [30]
