Index: inundation/ga/storm_surge/pyvolution/general_mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 321)
+++ inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 321)
@@ -0,0 +1,193 @@
+    
+class General_Mesh:
+    """Collection of triangular elements (purely geometric)
+
+    A triangular element is defined in terms of three vertex ids,
+    ordered counter clock-wise, 
+    each corresponding to a given coordinate set.
+    Vertices from different elements can point to the same
+    coordinate set.
+
+    Coordinate sets are implemented as an N x 2 Numeric array containing
+    x and y coordinates. 
+    
+
+    To instantiate:
+       Mesh(coordinates, triangles)
+
+    where
+
+      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
+      floats representing all x, y coordinates in the mesh.
+
+      triangles is either a list of 3-tuples or an Nx3 Numeric array of
+      integers representing indices of all vertices in the mesh.
+      Each vertex is identified by its index i in [0, M-1].
+
+
+    Example:
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        e = [2.0, 2.0]
+        
+        points = [a, b, c, e]
+        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        mesh = Mesh(points, triangles)
+    
+        #creates two triangles: bac and bce
+
+        This is a cut down version of mesh from pyvolution\mesh.py
+    """
+
+    def __init__(self, coordinates, triangles):
+        """
+        Build triangles from x,y coordinates (sequence of 2-tuples or
+        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
+        or Nx3 Numeric array of non-negative integers).
+        """
+
+        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
+
+        self.triangles = array(triangles).astype(Int)
+        self.coordinates = array(coordinates)
+
+        #Input checks
+        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.triangles.shape) == 2, msg
+
+        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.coordinates.shape) == 2, msg        
+
+        msg = 'Vertex indices reference non-existing coordinate sets'
+        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
+
+
+        #Register number of elements (N)
+        self.number_of_elements = N = self.triangles.shape[0]
+    
+        
+        self.normals = zeros((N, 6), Float)
+
+        #Get x,y coordinates for all triangles and store
+        self.vertex_coordinates = V = self.compute_vertex_coordinates()
+        
+        #Initialise each triangle
+        for i in range(N):
+            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
+            
+            x0 = V[i, 0]; y0 = V[i, 1]
+            x1 = V[i, 2]; y1 = V[i, 3]
+            x2 = V[i, 4]; y2 = V[i, 5]            
+
+            #Area
+            area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+
+            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
+            msg += ' is degenerate:  area == %f' %area
+            assert area > 0.0, msg
+        
+
+            #Normals
+            #The normal vectors
+            # - point outward from each edge
+            # - are orthogonal to the edge
+            # - have unit length
+            # - Are enumerated according to the opposite corner:
+            #   (First normal is associated with the edge opposite
+            #    the first vertex, etc)
+            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
+
+            n0 = array([x2 - x1, y2 - y1])
+	    l0 = sqrt(sum(n0**2))
+
+            n1 = array([x0 - x2, y0 - y2])
+            l1 = sqrt(sum(n1**2))
+
+            n2 = array([x1 - x0, y1 - y0])
+            l2 = sqrt(sum(n2**2))
+
+            #Normalise
+            n0 /= l0
+            n1 /= l1
+            n2 /= l2
+
+            #Compute and store
+            self.normals[i, :] = [n0[1], -n0[0],
+                                  n1[1], -n1[0],
+                                  n2[1], -n2[0]]
+            
+    def __len__(self):
+        return self.number_of_elements
+
+    def __repr__(self):
+        return 'Mesh: %d triangles, %d elements'\
+               %(self.coordinates.shape[0], len(self))
+
+    def get_normals(self):
+        """Return all normal vectors.
+        Return normal vectors for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.normals
+
+
+    def get_normal(self, i, j):
+        """Return normal vector j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.normals[i, 2*j:2*j+2]     
+    
+
+            
+    def get_vertex_coordinates(self):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.vertex_coordinates
+
+
+    def get_vertex_coordinate(self, i, j):
+        """Return coordinates for vertex j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.vertex_coordinates[i, 2*j:2*j+2]     
+
+
+    def compute_vertex_coordinates(self):
+        """Return vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+        """
+
+        #FIXME: Perhaps they should be ordered as in obj files??
+        #See quantity.get_vertex_values
+        
+        from Numeric import zeros, Float
+
+        N = self.number_of_elements
+        vertex_coordinates = zeros((N, 6), Float) 
+
+        for i in range(N):
+            for j in range(3):
+                k = self.triangles[i,j]  #Index of vertex 0
+                v_k = self.coordinates[k]
+                vertex_coordinates[i, 2*j+0] = v_k[0]
+                vertex_coordinates[i, 2*j+1] = v_k[1]
+
+        return vertex_coordinates
+  
+    def get_triangles(self, unique=False):
+        """Get connectivity
+        If unique is True give them only once as stored internally.
+        If unique is False return
+        """
+
+        if unique is True:
+            return self.triangles
+        else:
+	    from Numeric import reshape, array, Int
+	     	
+            m = len(self)  #Number of volumes
+            M = 3*m        #Total number of unique vertices
+	    return reshape(array(range(M)).astype(Int), (m,3))
Index: inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 320)
+++ inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 321)
@@ -13,227 +13,20 @@
 """
 
-#FIXME: What is the reason for duplication class Mesh
-#(or large parts of it at least) here?
-
-class Mesh:
-    """Collection of triangular elements (purely geometric)
-
-    A triangular element is defined in terms of three vertex ids,
-    ordered counter clock-wise, 
-    each corresponding to a given coordinate set.
-    Vertices from different elements can point to the same
-    coordinate set.
-
-    Coordinate sets are implemented as an N x 2 Numeric array containing
-    x and y coordinates. 
-    
-
-    To instantiate:
-       Mesh(coordinates, triangles)
-
-    where
-
-      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
-      floats representing all x, y coordinates in the mesh.
-
-      triangles is either a list of 3-tuples or an Nx3 Numeric array of
-      integers representing indices of all vertices in the mesh.
-      Each vertex is identified by its index i in [0, M-1].
+EPSILON = 1.0e-13
 
 
-    Example:
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        e = [2.0, 2.0]
-        
-        points = [a, b, c, e]
-        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
-        mesh = Mesh(points, triangles)
-    
-        #creates two triangles: bac and bce
-
-        This is a cut down version of mesh from pyvolution\mesh.py
-    """
-
-    def __init__(self, coordinates, triangles):
-        """
-        Build triangles from x,y coordinates (sequence of 2-tuples or
-        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
-        or Nx3 Numeric array of non-negative integers).
-        """
-
-        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
-
-        self.triangles = array(triangles).astype(Int)
-        self.coordinates = array(coordinates)
-
-        #Input checks
-        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
-        assert len(self.triangles.shape) == 2, msg
-
-        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
-        assert len(self.coordinates.shape) == 2, msg        
-
-        msg = 'Vertex indices reference non-existing coordinate sets'
-        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
-
-
-        #Register number of elements (N)
-        self.number_of_elements = N = self.triangles.shape[0]
-    
-        
-        self.normals = zeros((N, 6), Float)
-
-        #Get x,y coordinates for all triangles and store
-        self.vertex_coordinates = V = self.compute_vertex_coordinates()
-        
-        #Initialise each triangle
-        for i in range(N):
-            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
-            
-            x0 = V[i, 0]; y0 = V[i, 1]
-            x1 = V[i, 2]; y1 = V[i, 3]
-            x2 = V[i, 4]; y2 = V[i, 5]            
-
-            #Area
-            area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-
-            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
-            msg += ' is degenerate:  area == %f' %area
-            assert area > 0.0, msg
-        
-
-            #Normals
-            #The normal vectors
-            # - point outward from each edge
-            # - are orthogonal to the edge
-            # - have unit length
-            # - Are enumerated according to the opposite corner:
-            #   (First normal is associated with the edge opposite
-            #    the first vertex, etc)
-            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
-
-            n0 = array([x2 - x1, y2 - y1])
-	    l0 = sqrt(sum(n0**2))
-
-            n1 = array([x0 - x2, y0 - y2])
-            l1 = sqrt(sum(n1**2))
-
-            n2 = array([x1 - x0, y1 - y0])
-            l2 = sqrt(sum(n2**2))
-
-            #Normalise
-            n0 /= l0
-            n1 /= l1
-            n2 /= l2
-
-            #Compute and store
-            self.normals[i, :] = [n0[1], -n0[0],
-                                  n1[1], -n1[0],
-                                  n2[1], -n2[0]]             
-
-         
-        #Build vertex list
-        self.build_vertexlist()        
-
-
-        
-    def __len__(self):
-        return self.number_of_elements
-
-    def __repr__(self):
-        return 'Mesh: %d triangles, %d elements'\
-               %(self.coordinates.shape[0], len(self))
-
-
-
-    #DSG I don't know if I need this
-    def build_vertexlist(self):
-        pass
-    
-    # maybe add this latter..    
-    def check_integrity(self):
-        pass
-    
-
-    def get_normals(self):
-        """Return all normal vectors.
-        Return normal vectors for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
-        """
-        return self.normals
-
-
-    def get_normal(self, i, j):
-        """Return normal vector j of the i'th triangle.
-        Return value is the numeric array slice [x, y]
-        """
-        return self.normals[i, 2*j:2*j+2]     
-    
-
-            
-    def get_vertex_coordinates(self):
-        """Return all vertex coordinates.
-        Return all vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
-        """
-        return self.vertex_coordinates
-
-
-    def get_vertex_coordinate(self, i, j):
-        """Return coordinates for vertex j of the i'th triangle.
-        Return value is the numeric array slice [x, y]
-        """
-        return self.vertex_coordinates[i, 2*j:2*j+2]     
-
-
-    def compute_vertex_coordinates(self):
-        """Return vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-        """
-
-        #FIXME: Perhaps they should be ordered as in obj files??
-        #See quantity.get_vertex_values
-        
-        from Numeric import zeros, Float
-
-        N = self.number_of_elements
-        vertex_coordinates = zeros((N, 6), Float) 
-
-        for i in range(N):
-            for j in range(3):
-                k = self.triangles[i,j]  #Index of vertex 0
-                v_k = self.coordinates[k]
-                vertex_coordinates[i, 2*j+0] = v_k[0]
-                vertex_coordinates[i, 2*j+1] = v_k[1]
-
-        return vertex_coordinates
-  
-    def get_triangles(self, unique=False):
-        """Get connectivity
-        If unique is True give them only once as stored internally.
-        If unique is False return
-        """
-
-        if unique is True:
-            return self.triangles
-        else:
-	    from Numeric import reshape, array, Int
-	     	
-            m = len(self)  #Number of volumes
-            M = 3*m        #Total number of unique vertices
-	    return reshape(array(range(M)).astype(Int), (m,3))
-
-
-
+from general_mesh import General_Mesh
 from Numeric import zeros, array, Float, Int, dot,transpose
 from LinearAlgebra import solve_linear_equations
 
-EPSILON = 1.0e-13
 
-class Interpolation(Mesh):
+class Interpolation:
 
-    def __init__(self, vertex_coordinates, triangles, point_coordinates):
+    def __init__(self,
+                 vertex_coordinates,
+                 triangles,
+                 point_coordinates = None):
+
+        
         """ Build interpolation matrix mapping from
         function values at vertices to function values at data points
@@ -253,14 +46,22 @@
 
         #Convert input to Numeric arrays
-        point_coordinates = array(point_coordinates).astype(Float)
         vertex_coordinates = array(vertex_coordinates).astype(Float)
         triangles = array(triangles).astype(Int)                
         
         #Build underlying mesh
-        Mesh.__init__(self, vertex_coordinates, triangles)
+        self.mesh = General_Mesh(vertex_coordinates, triangles)
+
+        self.A_m = zeros((1,1), Float) #Not initialised
+
+        if point_coordinates:
+            self.build_A_m(point_coordinates)
         
+    def build_A_m(self, point_coordinates):
 
+        #Convert input to Numeric arrays
+        point_coordinates = array(point_coordinates).astype(Float)
+        
         #Build n x m interpolation matrix        
-        m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex)
+        m = self.mesh.coordinates.shape[0] #Number of basis functions (1/vertex)
         n = point_coordinates.shape[0]               #Number of data points         
 
@@ -272,17 +73,17 @@
 
             x = point_coordinates[i]
-            for k in range(len(self)):
+            for k in range(len(self.mesh)):
                 #For each triangle (brute force)
                 #FIXME: Real algorithm should only visit relevant triangles
 
                 #Get the three vertex_points
-                xi0 = self.get_vertex_coordinate(k, 0)
-                xi1 = self.get_vertex_coordinate(k, 1)
-                xi2 = self.get_vertex_coordinate(k, 2)                 
+                xi0 = self.mesh.get_vertex_coordinate(k, 0)
+                xi1 = self.mesh.get_vertex_coordinate(k, 1)
+                xi2 = self.mesh.get_vertex_coordinate(k, 2)                 
 
                 #Get the three normals
-                n0 = self.get_normal(k, 0)
-                n1 = self.get_normal(k, 1)
-                n2 = self.get_normal(k, 2)                
+                n0 = self.mesh.get_normal(k, 0)
+                n1 = self.mesh.get_normal(k, 1)
+                n2 = self.mesh.get_normal(k, 2)                
 
                 #Compute interpolation
@@ -299,16 +100,19 @@
 
                     #Assign values to A_m
-                    j = self.triangles[k,0] #Global vertex id
+                    j = self.mesh.triangles[k,0] #Global vertex id
                     self.A_m[i, j] = sigma0
 
-                    j = self.triangles[k,1] #Global vertex id
+                    j = self.mesh.triangles[k,1] #Global vertex id
                     self.A_m[i, j] = sigma1
 
-                    j = self.triangles[k,2] #Global vertex id
+                    j = self.mesh.triangles[k,2] #Global vertex id
                     self.A_m[i, j] = sigma2
+    
 
-    def smooth_to_mesh(self, z):
-        """ Linearly interpolate point data values to the vertices.
-        
+    def smooth_att_to_mesh(self, z):
+        """ Linearly interpolate many points with an attribute to the vertices.
+        Pre Condition:
+          A_m has been initialised
+          
         Inputs:
           z: Vector of data at the point_coordinates.  One value per point.
@@ -325,6 +129,32 @@
         return f
 
-    def smooth_to_points(self, f):
-        """ Linearly interpolate point data values to the vertices.
+    def smooth_atts_to_mesh(self, z_m):
+        """ Linearly interpolate many points with attributes to the vertices.
+        Pre Condition:
+          A_m has been initialised
+          
+        Inputs:
+          z_m: A Matrix of data at the point_coordinates.
+               One attribute per colum.
+          
+        """
+        #Convert input to Numeric arrays
+        z_m = array(z_m).astype(Float)
+
+        
+        #Build n x m interpolation matrix        
+        m = self.mesh.coordinates.shape[0] #Number of vertices
+        n = z_m.shape[1]               #Number of data points         
+
+        f_m = zeros((n,m), Float)
+        
+        for i in range(z_m.shape[1]):
+            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
+        
+        
+    def smooth_att_to_points(self, f):
+        """ Linearly interpolate vertices with an attribute to the points.
+        Pre Condition:
+          A_m has been initialised
         
         Inputs:
@@ -333,2 +163,26 @@
         """
         return dot(self.A_m,f)
+
+    def smooth_atts_to_points(self, f_m):
+        """ Linearly interpolate vertices with attributes to the points.
+        Pre Condition:
+          A_m has been initialised
+        
+        Inputs:
+          f_m: Matrix of data at the  vertices. 
+               One attribute per colum. 
+          
+        """
+        #Convert input to Numeric arrays
+        z_m = array(z_m).astype(Float)
+
+        
+        #Build n x m interpolation matrix        
+        m = self.mesh.coordinates.shape[0] #Number of vertices
+        n = z_m.shape[1]               #Number of data points         
+
+        f_m = zeros((n,m), Float)
+        
+        for i in range(z_m.shape[1]):
+            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
+        
Index: inundation/ga/storm_surge/pyvolution/test_least_squares.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 320)
+++ inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 321)
@@ -5,5 +5,5 @@
 
 from least_squares import *
-from config import epsilon
+from least_squares import EPSILON
 from Numeric import allclose, array
 
@@ -28,6 +28,6 @@
         data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point        
         
-        mesh = Interpolation(points, vertices, data)
-        assert allclose(mesh.A_m, [[1./3, 1./3, 1./3]]) 
+        interp = Interpolation(points, vertices, data)
+        assert allclose(interp.A_m, [[1./3, 1./3, 1./3]]) 
         
 
@@ -45,6 +45,6 @@
         data = points #Use data at vertices        
         
-        mesh = Interpolation(points, vertices, data)
-        assert allclose(mesh.A_m, [[1., 0., 0.],
+        interp = Interpolation(points, vertices, data)
+        assert allclose(interp.A_m, [[1., 0., 0.],
                                       [0., 1., 0.],
                                       [0., 0., 1.]]) 
@@ -65,7 +65,7 @@
         data = [ [0., 1.], [1., 0.], [1., 1.] ]
         
-        mesh = Interpolation(points, vertices, data)
-        
-        assert allclose(mesh.A_m, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
+        interp = Interpolation(points, vertices, data)
+        
+        assert allclose(interp.A_m, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
                                       [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
                                       [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
@@ -85,7 +85,7 @@
         data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
         
-        mesh = Interpolation(points, vertices, data)
-        
-        assert allclose(mesh.A_m, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0      
+        interp = Interpolation(points, vertices, data)
+        
+        assert allclose(interp.A_m, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0      
                                       [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
                                       [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
@@ -105,7 +105,7 @@
         data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
         
-        mesh = Interpolation(points, vertices, data)
-
-        assert allclose(sum(mesh.A_m, axis=1), 1.0)
+        interp = Interpolation(points, vertices, data)
+
+        assert allclose(sum(interp.A_m, axis=1), 1.0)
         
 
@@ -125,5 +125,5 @@
         data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
 
-        mesh = Interpolation(points, triangles, data) 
+        interp = Interpolation(points, triangles, data) 
         
         answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],  #Affects point d     
@@ -133,10 +133,10 @@
                   [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
                   [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f       
-        #print mesh.A_m
+        #print interp.A_m
         #print answer
         
-        assert allclose(mesh.A_m, answer)
-
-    def test_smooth_to_mesh(self):
+        assert allclose(interp.A_m, answer)
+
+    def test_smooth_att_to_mesh(self):
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -153,11 +153,11 @@
         data_coords = [ d1, d2, d3] #Use centroid as one data point        
         
-        mesh = Interpolation(points, triangles, data_coords)
+        interp = Interpolation(points, triangles, data_coords)
         z = [z1, z2, z3]
-        f =  mesh.smooth_to_mesh(z)
+        f =  interp.smooth_att_to_mesh(z)
         answer = [0, 5., 5.]
         assert allclose(f, answer)
         
-    def test_smooth_to_meshII(self):
+    def test_smooth_att_to_meshII(self):
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -171,10 +171,10 @@
         data_coords = [ d1, d2, d3]       
         z = self.linear_function(data_coords)
-        mesh = Interpolation(points, triangles, data_coords)
-        f =  mesh.smooth_to_mesh(z)
+        interp = Interpolation(points, triangles, data_coords)
+        f =  interp.smooth_att_to_mesh(z)
         answer = self.linear_function(points)
         assert allclose(f, answer)
     
-    def test_smooth_to_meshIII(self):
+    def test_smooth_att_to_meshIII(self):
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -202,6 +202,6 @@
         
         z = self.linear_function(point_coords)
-        mesh = Interpolation(vertices, triangles, point_coords)
-        f =  mesh.smooth_to_mesh(z)
+        interp = Interpolation(vertices, triangles, point_coords)
+        f =  interp.smooth_att_to_mesh(z)
         answer = self.linear_function(vertices)
         assert allclose(f, answer)
@@ -211,5 +211,5 @@
         return point[:,0]+point[:,1]
     
-    def test_smooth_to_points(self):
+    def test_smooth_att_to_points(self):
         v0 = [0.0, 0.0]
         v1 = [0.0, 5.0]
@@ -225,7 +225,7 @@
         point_coords = [ d0, d1, d2]
         
-        mesh = Interpolation(vertices, triangles, point_coords)
+        interp = Interpolation(vertices, triangles, point_coords)
         f = self.linear_function(vertices)
-        z =  mesh.smooth_to_points(f)
+        z =  interp.smooth_att_to_points(f)
         answer = self.linear_function(point_coords)
         #print "z",z
@@ -233,5 +233,5 @@
         assert allclose(z, answer)
        
-    def test_smooth_to_pointsII(self):
+    def test_smooth_att_to_pointsII(self):
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -258,11 +258,34 @@
                         [3.0,1.0]]
         
-        mesh = Interpolation(vertices, triangles, point_coords)
+        interp = Interpolation(vertices, triangles, point_coords)
         f = self.linear_function(vertices)
-        z =  mesh.smooth_to_points(f)
+        z =  interp.smooth_att_to_points(f)
         answer = self.linear_function(point_coords)
         #print "z",z
         #print "answer",answer
-        assert allclose(z, answer) 
+        assert allclose(z, answer)
+
+        
+    def test_smooth_atts_to_mesh(self):
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+         
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 3.0]
+        d3 = [3.0,1.0]
+        z1 = [2,4]
+        z2 = [4,8]
+        z3 = [4,8]
+        data_coords = [ d1, d2, d3] #Use centroid as one data point        
+        
+        interp = Interpolation(points, triangles, data_coords)
+        z = [z1, z2, z3]
+        f =  interp.smooth_att_to_mesh(z)
+        answer = [[0,0], [5., 10.], [5., 10.]]
+        assert allclose(f, answer)
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
