Index: /inundation/pyvolution/least_squares.py
===================================================================
--- /inundation/pyvolution/least_squares.py	(revision 2584)
+++ /inundation/pyvolution/least_squares.py	(revision 2585)
@@ -216,21 +216,39 @@
     
     if verbose:
+    
+        point_coordinates = ensure_numeric(point_coordinates)
+        vertex_coordinates = ensure_numeric(vertex_coordinates)
+		
+        X = point_coordinates[:,0]
+        Y = point_coordinates[:,1]	
+        Z = point_attributes	
+	    
         print '+------------------------------------------------'
-        print 'least squares finished'
-        print 'points: %d points' %(len(point_attributes))
-        print '    x in [%f, %f]'%(min(point_coordinates[:,0]),
-                                   max(point_coordinates[:,0]))
-        print '    y in [%f, %f]'%(min(point_coordinates[:,1]),
-                                   max(point_coordinates[:,1]))
-        print '    z in [%f, %f]'%(min(point_attributes),
-                                   max(point_attributes))
-        print
-        print 'mesh: %d vertices' %(len(vertex_attributes))
-        print '    xi in [%f, %f]'%(min(vertex_coordinates[:,0]),
-                                   max(vertex_coordinates[:,0]))
-        print '    eta in [%f, %f]'%(min(vertex_coordinates[:,1]),
-                                   max(vertex_coordinates[:,1]))
-        print '    zeta in [%f, %f]'%(min(vertex_attributes),
-                                   max(vertex_attributes))
+        print 'Least squares statistics'
+        print '+------------------------------------------------'	
+        print 'points: %d points' %(len(Z))
+        print '    x in [%f, %f]'%(min(X), max(X))
+        print '    y in [%f, %f]'%(min(Y), max(Y))
+        print '    z in [%f, %f]'%(min(Z), max(Z))
+	print
+
+        indices = interp.point_indices
+	if indices is not None:
+            X = take(X, indices)
+            Y = take(Y, indices)	
+            Z = take(Z, indices)	
+            print 'Cropped points: %d points' %(len(Z))
+            print '    x in [%f, %f]'%(min(X), max(X))
+            print '    y in [%f, %f]'%(min(Y), max(Y))
+            print '    z in [%f, %f]'%(min(Z), max(Z))
+            print
+	
+	Xi = vertex_coordinates[:,0]
+        Eta = vertex_coordinates[:,1]	
+        Zeta = vertex_attributes		
+        print 'Mesh: %d vertices' %(len(Zeta))
+        print '    xi in [%f, %f]'%(min(Xi), max(Xi))
+        print '    eta in [%f, %f]'%(min(Eta), max(Eta))
+        print '    zeta in [%f, %f]'%(min(Zeta), max(Zeta))
         print '+------------------------------------------------'
 
Index: /inundation/pyvolution/test_least_squares.py
===================================================================
--- /inundation/pyvolution/test_least_squares.py	(revision 2584)
+++ /inundation/pyvolution/test_least_squares.py	(revision 2585)
@@ -725,5 +725,5 @@
         z = [2, 4, 4]
 
-        ref = fit_to_mesh(points, triangles, data_coords, z)
+        ref = fit_to_mesh(points, triangles, data_coords, z, verbose=False)
 
         #print attributes
@@ -943,5 +943,5 @@
         #Fit surface to mesh
         interp = Interpolation(points, triangles, data_points1, alpha=0.0,
-                               precrop = True)
+                               precrop = True, verbose=False)
         z = linear_function(data_points1) #Example z-values
         f = interp.fit(z)                 #Fitted values at vertices
@@ -1817,4 +1817,48 @@
         os.remove(point_file)
 
+	
+	
+    def test_fit_using_fit_to_mesh(self):
+        """Fit a surface to one set of points. Then interpolate that surface
+        using another set of points.
+        """
+        from mesh import Mesh
+
+
+        #Setup mesh used to represent fitted function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Datapoints to fit from
+        data_points1 = [[ 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],
+                        [ 15, -17],   #Outside mesh
+                        [ 1.0, 2.0],
+                        [ 1.0, 3.0],
+                        [ 2.0, 1.0],
+                        [ 3.0, 0.0],
+                        [ 3.0, 1.0]]
+
+        #Fit surface to mesh
+        z = linear_function(data_points1) #Example z-values
+	v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0, 
+	                precrop=True, verbose=False)
+	assert allclose(linear_function(points), v)
+
+	
 #-------------------------------------------------------------
 if __name__ == "__main__":
