source: inundation/ga/storm_surge/pyvolution/least_squares.py @ 303

"""Least squares smooting and interpolation.
   Very first cut - just to check a few things.
   Architecture is not set in concrete at all.

   O
"""


from mesh import Mesh


class Interpolation(Mesh):

    def __init__(self, vertex_coordinates, triangles, data):
        """ Build interpolation matrix mapping from
        function values at vertices to function values at data points

        Inputs:
        
          vertex_coordinates: List of coordinate pairs [xi, eta] of points
          constituting mesh (or a an m x 2 Numeric array)
        
          triangles: List of 3-tuples (or a Numeric array) of
          integers representing indices of all vertices in the mesh.

          data: List of coordinate pairs [x, y] of data points
          (or an nx2 Numeric array)
          
        """

        from Numeric import zeros, array, Float, Int, dot
        from config import epsilon

        #Convert input to Numeric arrays
        data = array(data).astype(Float)
        vertex_coordinates = array(vertex_coordinates).astype(Float)
        triangles = array(triangles).astype(Int)                
        
        #Build underlying mesh
        Mesh.__init__(self, vertex_coordinates, triangles)
        

        #Build n x m interpolation matrix        
        m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex)
        n = data.shape[0]               #Number of data points         

        self.matrix = zeros((n,m), Float)

        #Compute matrix elements
        for i in range(n):
            #For each data point

            x = data[i]
            for k in range(len(self)):
                #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)                 

                #Get the three normals
                n0 = self.get_normal(k, 0)
                n1 = self.get_normal(k, 1)
                n2 = self.get_normal(k, 2)                

                #Compute interpolation
                sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
                sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
                sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)

                #FIXME: Maybe move out to test or something
                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon

                
                #Check that this trinagle contains data point
                if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:

                    #Assign values to matrix
                    v_id = self.vertices[k,0]
                    self.matrix[i, v_id] = sigma0

                    v_id = self.vertices[k,1]
                    self.matrix[i, v_id] = sigma1

                    v_id = self.vertices[k,2]
                    self.matrix[i, v_id] = sigma2