Index: /inundation/pyvolution/data_manager.py
===================================================================
--- /inundation/pyvolution/data_manager.py	(revision 1752)
+++ /inundation/pyvolution/data_manager.py	(revision 1753)
@@ -963,5 +963,5 @@
 
     import os 
-    from Scientific.IO.NetCDF import NetCDFFile
+    #from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import Float, arrayrange, concatenate
 
@@ -1010,6 +1010,22 @@
 
     if verbose: print 'Store to NetCDF file %s' %ptsname
+    write_ptsfile(ptsname, points, attribute, attribute_name)
+
+
+def write_ptsfile(ptsname, points, attribute, attribute_name = None):
+    """Write points and associated attribute to pts (NetCDF) format
+    """
+
+    from Numeric import Float
+    
+    if attribute_name is None:
+        attribute_name = 'attribute'        
+    
+
+    from Scientific.IO.NetCDF import NetCDFFile
+    
     # NetCDF file definition
     outfile = NetCDFFile(ptsname, 'w')
+
 
     #Create new file
@@ -1040,7 +1056,7 @@
     nc_attribute[:] = attribute
 
-    infile.close()
     outfile.close()
     
+
 def dem2pts(basename_in, basename_out=None, verbose=False,
             easting_min=None, easting_max=None,
@@ -1063,4 +1079,6 @@
     elevation: N Float array
     """
+
+    #FIXME: Can this be written feasibly using write_pts?
 
     import os
Index: /inundation/pyvolution/domain.py
===================================================================
--- /inundation/pyvolution/domain.py	(revision 1752)
+++ /inundation/pyvolution/domain.py	(revision 1753)
@@ -16,5 +16,6 @@
     def __init__(self, coordinates, vertices, boundary = None,
                  conserved_quantities = None, other_quantities = None,
-                 tagged_elements = None, geo_reference = None, use_inscribed_circle=False):
+                 tagged_elements = None, geo_reference = None,
+                 use_inscribed_circle=False):
 
         Mesh.__init__(self, coordinates, vertices, boundary,
@@ -141,8 +142,8 @@
 
 
-    def set_quantity(self, name,
-                     X,
-                     location = 'vertices',
-                     indices = None):
+    def set_quantity(self, name, *args, **kwargs):
+                     #X,
+                     #location = 'vertices',
+                     #indices = None):
         """Set values for named quantity
 
@@ -196,8 +197,8 @@
         #Set value
         #FIXME: use **kwargs
-        self.quantities[name].set_values(X,
-                                         location = location,
-                                         indices = indices)
-
+        #self.quantities[name].set_values(X,
+        #                                 location = location,
+        #                                 indices = indices)
+        self.quantities[name].set_values(*args, **kwargs)
 
 
Index: /inundation/pyvolution/quantity.py
===================================================================
--- /inundation/pyvolution/quantity.py	(revision 1752)
+++ /inundation/pyvolution/quantity.py	(revision 1753)
@@ -91,5 +91,6 @@
                    location = 'vertices',                   
                    indices = None,
-                   verbose = None):
+                   verbose = None,
+                   use_cache = False):
         
         """Set values for quantity based on different sources.
@@ -155,4 +156,11 @@
                   will be left undefined.
 
+        verbose: True means that output to stdout is generated
+
+        use_cache: True means that caching of intermediate results is
+                   attempted for least squares fit.
+                   
+                   
+        
 
         Exactly one of the arguments
@@ -211,8 +219,10 @@
             assert values is not None, msg
             self.set_values_from_points(points, values, alpha,
-                                        location, indices, verbose)
+                                        location, indices, verbose,
+                                        use_cache)
         elif filename is not None:
             self.set_values_from_file(filename, attribute_name, alpha,
-                                      location, indices, verbose)
+                                      location, indices, verbose,
+                                      use_cache)
         else:
             raise 'This can\'t happen :-)'
@@ -437,10 +447,9 @@
 
     def set_values_from_points(self, points, values, alpha,
-                               location, indices, verbose):
-        """
-        """
-
-        
-        #FIXME: Needs unit test         
+                               location, indices, verbose, use_cache):
+        """Set quantity values from arbitray data points using least squares
+        """
+
+        from Numeric import Float 
         from util import ensure_numeric
         from least_squares import fit_to_mesh
@@ -457,11 +466,25 @@
         triangles = self.domain.triangles
 
-        #FIXME Pass and use caching here
-        vertex_attributes = fit_to_mesh(coordinates,
-                                        triangles,
-                                        points,
-                                        values,
-                                        alpha = alpha,
-                                        verbose = verbose)
+        if use_cache is True:
+            try:
+                from caching import cache
+            except:
+                msg = 'Caching was requested, but caching module'+\
+                      'could not be imported'
+                raise msg
+            
+            args = (coordinates, triangles, points, values)
+            kwargs = {'alpha': alpha, 'verbose': verbose}
+            vertex_attributes = cache(fit_to_mesh,
+                                      args, kwargs,
+                                      verbose = verbose)
+        else:    
+            vertex_attributes = fit_to_mesh(coordinates,
+                                            triangles,
+                                            points,
+                                            values,
+                                            alpha = alpha,
+                                            verbose = verbose)
+            
         
         self.set_values_from_array(vertex_attributes,
@@ -473,5 +496,5 @@
     
     def set_values_from_file(self, filename, attribute_name, alpha,
-                             location, indices, verbose):
+                             location, indices, verbose, use_cache):
         """Set quantity based on arbitrary points in .pts file
         using least_squares attribute_name selects name of attribute
@@ -481,5 +504,4 @@
 
 
-        #FIXME: Needs unit test 
         from types import StringType
         msg = 'Filename must be a text string'
@@ -513,5 +535,5 @@
         #Call least squares method    
         self.set_values_from_points(points, z, alpha,
-                                    location, indices, verbose)
+                                    location, indices, verbose, use_cache)
 
 
@@ -523,5 +545,5 @@
         location: Where values are to be stored.
                   Permissible options are: vertices, edges, centroid
-                  Default is "vertices"
+                  and unique vertices. Default is 'vertices'
 
         In case of location == 'centroids' the dimension values must
Index: /inundation/pyvolution/test_advection.py
===================================================================
--- /inundation/pyvolution/test_advection.py	(revision 1752)
+++ /inundation/pyvolution/test_advection.py	(revision 1753)
@@ -77,5 +77,5 @@
         domain.check_integrity()
 
-        domain.set_quantity('stage', [1.0], 'centroids')
+        domain.set_quantity('stage', [1.0], location='centroids')
 
 	domain.distribute_to_vertices_and_edges()
@@ -136,5 +136,5 @@
         #Populate boundary array with dirichlet conditions.
         domain.neighbours = array([[1,-1,-2], [0,-3,-4]])
-        domain.set_quantity('stage', [1.0, 0.0], 'centroids')
+        domain.set_quantity('stage', [1.0, 0.0], location='centroids')
         domain.distribute_to_vertices_and_edges()
 
Index: /inundation/pyvolution/test_data_manager.py
===================================================================
--- /inundation/pyvolution/test_data_manager.py	(revision 1752)
+++ /inundation/pyvolution/test_data_manager.py	(revision 1753)
@@ -492,4 +492,56 @@
         #Cleanup
         os.remove(sww.filename)
+
+
+    def test_write_pts(self):
+        #Get (enough) datapoints
+
+        from Numeric import array
+        points = array([[ 0.66666667, 0.66666667],
+                        [ 1.33333333, 1.33333333],
+                        [ 2.66666667, 0.66666667],
+                        [ 0.66666667, 2.66666667],
+                        [ 0.0, 1.0],
+                        [ 0.0, 3.0],
+                        [ 1.0, 0.0],
+                        [ 1.0, 1.0],
+                        [ 1.0, 2.0],
+                        [ 1.0, 3.0],
+                        [ 2.0, 1.0],
+                        [ 3.0, 0.0],
+                        [ 3.0, 1.0]])
+
+        z = points[:,0] + 2*points[:,1]
+
+        ptsfile = 'testptsfile.pts'
+        write_ptsfile(ptsfile, points, z,
+                      attribute_name = 'linear_combination')
+
+        #Check contents
+        #Get NetCDF
+        from Scientific.IO.NetCDF import NetCDFFile        
+        fid = NetCDFFile(ptsfile, 'r')
+
+        # Get the variables
+        #print fid.variables.keys()
+        points1 = fid.variables['points']
+        z1 = fid.variables['linear_combination']
+
+        #Check values
+
+        #print points[:]
+        #print ref_points
+        assert allclose(points, points1)
+
+        #print attributes[:]
+        #print ref_elevation
+        assert allclose(z, z1)
+
+        #Cleanup
+        fid.close()
+
+        import os
+        os.remove(ptsfile)
+        
 
 
Index: /inundation/pyvolution/test_domain.py
===================================================================
--- /inundation/pyvolution/test_domain.py	(revision 1752)
+++ /inundation/pyvolution/test_domain.py	(revision 1753)
@@ -290,7 +290,7 @@
 
 
-        domain.set_quantity('stage', [1,2,3,4], 'centroids')
-        domain.set_quantity('xmomentum', [1,2,3,4], 'centroids')
-        domain.set_quantity('ymomentum', [1,2,3,4], 'centroids')
+        domain.set_quantity('stage', [1,2,3,4], location='centroids')
+        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
+        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
 
 
Index: /inundation/pyvolution/test_quantity.py
===================================================================
--- /inundation/pyvolution/test_quantity.py	(revision 1752)
+++ /inundation/pyvolution/test_quantity.py	(revision 1753)
@@ -276,7 +276,8 @@
 
 
+        
 
     def test_set_vertex_values_using_least_squares(self):
-        from least_squares import Interpolation, fit_to_mesh
+
         
         quantity = Quantity(self.mesh4)
@@ -297,27 +298,29 @@
                        [ 3.0, 1.0]]
 
-
-        interp = Interpolation(quantity.domain.coordinates, quantity.domain.triangles,
-                               data_points, alpha=0.0)
-
         z = linear_function(data_points)
-        answer = linear_function(quantity.domain.coordinates)
-
-        f = interp.fit(z)
-
+        
+        #Obsoleted bit
+        #interp = Interpolation(quantity.domain.coordinates,
+        #                       quantity.domain.triangles,
+        #                       data_points, alpha=0.0)
+        #
+        #answer = linear_function(quantity.domain.coordinates)
+        #
+        #f = interp.fit(z)
+        #
         #print "f",f
         #print "answer",answer
-        assert allclose(f, answer)
-
-
-        quantity.set_values(f)
-
-
+        #assert allclose(f, answer)
+
+
+        #Use built-in least squares fit
+        quantity.set_values(points = data_points, values = z, alpha = 0)
         answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
         #print quantity.vertex_values, answer
         assert allclose(quantity.vertex_values.flat, answer)
 
-        #Now try using the general interface
-
+
+        #Now try by setting the same values directly
+        from least_squares import fit_to_mesh        
         vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
                                         quantity.domain.triangles,
@@ -329,6 +332,50 @@
         #print vertex_attributes
         quantity.set_values(vertex_attributes)
+        assert allclose(quantity.vertex_values.flat, answer)
+
+
+    def test_set_values_from_file(self):
+        quantity = Quantity(self.mesh4)
+
+        #Get (enough) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667],
+                       [ 0.0, 1.0],
+                       [ 0.0, 3.0],
+                       [ 1.0, 0.0],
+                       [ 1.0, 1.0],
+                       [ 1.0, 2.0],
+                       [ 1.0, 3.0],
+                       [ 2.0, 1.0],
+                       [ 3.0, 0.0],
+                       [ 3.0, 1.0]]
+
+        z = linear_function(data_points)
+        
+
+        #Create pts file
+        from data_manager import write_ptsfile
+        ptsfile = 'testptsfile.pts'
+        att = 'spam_and_eggs'
+        write_ptsfile(ptsfile, data_points, z, attribute_name = att)
+
+        #Check that values can be set from file
+        quantity.set_values(filename = ptsfile,
+                            attribute_name = att, alpha = 0)
+        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
+        #print quantity.vertex_values, answer
+        assert allclose(quantity.vertex_values.flat, answer)
+
+
+        #Check that values can be set from file using default attribute
+        quantity.set_values(filename = ptsfile, alpha = 0)
         assert allclose(quantity.vertex_values.flat, answer)        
 
+        #Cleanup
+        import os
+        os.remove(ptsfile)
+        
 
 
@@ -832,6 +879,6 @@
         indices = [1]
         quantity.set_values(value,
-                                  location = 'centroids',
-                                  indices = indices)
+                            location = 'centroids',
+                            indices = indices)
         #print "quantity.centroid_values",quantity.centroid_values
         assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
Index: /inundation/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/pyvolution/test_shallow_water.py	(revision 1752)
+++ /inundation/pyvolution/test_shallow_water.py	(revision 1753)
@@ -1238,5 +1238,6 @@
         val3 = 2.
 
-        domain.set_quantity('stage', [val0, val1, val2, val3], 'centroids')
+        domain.set_quantity('stage', [val0, val1, val2, val3],
+                            location='centroids')
         L = domain.quantities['stage'].vertex_values
 
@@ -1274,5 +1275,5 @@
             return x**2
 
-	domain.set_quantity('stage', stage, 'centroids')
+	domain.set_quantity('stage', stage, location='centroids')
 
         a, b = domain.quantities['stage'].compute_gradients()
@@ -1311,5 +1312,5 @@
             return x**4+y**2
 
-	domain.set_quantity('stage', stage, 'centroids')
+	domain.set_quantity('stage', stage, location='centroids')
         #print domain.quantities['stage'].centroid_values
 
@@ -1355,5 +1356,5 @@
 
         domain.set_quantity('elevation', slope)
-        domain.set_quantity('stage', stage, 'centroids')
+        domain.set_quantity('stage', stage, location='centroids')
 
         #print domain.quantities['elevation'].centroid_values
@@ -1443,12 +1444,15 @@
         domain.set_quantity('elevation', slope)
         domain.set_quantity('stage',
-			    [0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
-			    'centroids')
+			    [0.01298164, 0.00365611,
+                             0.01440365, -0.0381856437096],
+			    location='centroids')
         domain.set_quantity('xmomentum',
-			    [0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
-			    'centroids')
+			    [0.00670439, 0.01263789,
+                             0.00647805, 0.0178180740668],
+			    location='centroids')
         domain.set_quantity('ymomentum',
-			    [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
-			    'centroids')
+			    [-7.23510980e-004, -6.30413883e-005,
+                             6.30413883e-005, 0.000200907255866],
+			    location='centroids')
 
         E = domain.quantities['elevation'].vertex_values
@@ -1466,7 +1470,13 @@
 	#print Y[1,:]
 
-        assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
-	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
-	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
+        assert allclose(L[1,:], [-0.00825735775384,
+                                 -0.00801881482869,
+                                 0.0272445025825])
+	assert allclose(X[1,:], [0.0143507718962,
+                                 0.0142502147066,
+                                 0.00931268339717])
+	assert allclose(Y[1,:], [-0.000117062180693,
+                                 7.94434448109e-005,
+                                 -0.000151505429018])
 
 
@@ -2211,7 +2221,7 @@
 
 
-                domain.set_quantity('stage', L, 'centroids')
-                domain.set_quantity('xmomentum', X, 'centroids')
-                domain.set_quantity('ymomentum', Y, 'centroids')
+                domain.set_quantity('stage', L, location='centroids')
+                domain.set_quantity('xmomentum', X, location='centroids')
+                domain.set_quantity('ymomentum', Y, location='centroids')
 
                 domain.check_integrity()
Index: /inundation/validation/LWRU2/lwru2.py
===================================================================
--- /inundation/validation/LWRU2/lwru2.py	(revision 1752)
+++ /inundation/validation/LWRU2/lwru2.py	(revision 1753)
@@ -11,6 +11,6 @@
 
 
-read_mesh = True  #Use large unstructured mesh generated from create_mesh.py
-#read_mesh = False  #Use small structured mesh
+#read_mesh = True  #Use large unstructured mesh generated from create_mesh.py
+read_mesh = False  #Use small structured mesh
 
 import sys
@@ -128,29 +128,12 @@
 domain.filename, _ = os.path.splitext(base)
 
-#LS code to be included in set_quantity
-print 'Reading points'
-from pyvolution import util, least_squares
-#points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
-points, attributes = cache(util.read_xya,
-                           filenames.bathymetry_filename[:-4] + '.pts')
+print 'Initial values'
 
-elevation = attributes['elevation']
-
-#Fit attributes to mesh
-#vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
-#                                              domain.triangles,
-#                                              points,
-#                                              attributes['elevation'],
-#                                              alpha = 0.0001,
-#                                              verbose=True)
-
-vertex_attributes = cache(least_squares.fit_to_mesh,
-                          (domain.coordinates, domain.triangles,
-                           points, elevation) ,
-                          {'alpha': 0.0001, 'verbose': True})
-
-print 'Initial values'
-domain.set_quantity('elevation', vertex_attributes)
-#domain.set_quantity('elevation', points=points, attributes=attributes['elevation'], alpha=0.0001)
+domain.set_quantity('elevation',
+                    filename=filenames.bathymetry_filename[:-4] + '.pts',
+                    alpha = 0.0001,
+                    verbose = True,
+                    use_cache = True)
+                    
 domain.set_quantity('friction', 0.0)
 domain.set_quantity('stage', 0.0)
