Changeset 2684

Timestamp:

Apr 10, 2006, 3:15:05 PM (17 years ago)

Author:

duncan

Message:

Moving methods around, comments, a new test (to use when interpolate_function can handle discontinuous meshes)

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (6 diffs)
• test_interpolate.py (modified) (2 diffs)

inundation/fit_interpolate/interpolate.py

@@ -13 +13 @@
 
 TO DO
-* remove points outside the mesh ?(in interpolate_block)?
 * geo-ref (in interpolate_block)
-* add functional interpolate interface - in mesh and points, out interp data 
 """
 
@@ -68 +66 @@
 
         """
-
-        # why no alpha in parameter list?
-        
         # Initialise variabels
         self._A_can_be_reused = False
@@ -91 +86 @@
         # points, geo-ref and triangles, rather than just points and geo-ref
         # this info will usually be coming from a domain instance.
+        
         if mesh_origin is None:
             geo = None
@@ -104 +100 @@
         
         
-    def _build_interpolation_matrix_A(self,
-                                     point_coordinates,
-                                     verbose = False):
-        """Build n x m interpolation matrix, where
-        n is the number of data points and
-        m is the number of basis functions phi_k (one per vertex)
-
-        This algorithm uses a quad tree data structure for fast binning
-        of data points
-        origin is a 3-tuple consisting of UTM zone, easting and northing.
-        If specified coordinates are assumed to be relative to this origin.
-
-        This one will override any data_origin that may be specified in
-        instance interpolation
-
-        Preconditions
-        Point_coordindates and mesh vertices have the same origin.
-        """
-
-
-        
-        #Convert point_coordinates to Numeric arrays, in case it was a list.
-        point_coordinates = ensure_numeric(point_coordinates, Float)
-        
-        if verbose: print 'Getting indices inside mesh boundary'
-        self.inside_poly_indices, self.outside_poly_indices  = \
-                     in_and_outside_polygon(point_coordinates,
-                                            self.mesh.get_boundary_polygon(),
-                                            closed = True, verbose = verbose)
-        
-        #Build n x m interpolation matrix
-        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
-        n = point_coordinates.shape[0]     #Nbr of data points
-
-        if verbose: print 'Number of datapoints: %d' %n
-        if verbose: print 'Number of basis functions: %d' %m
-
-        A = Sparse(n,m)
-
-        #Compute matrix elements for points inside the mesh
-        for i in self.inside_poly_indices:
-            #For each data_coordinate point
-            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
-            x = point_coordinates[i]
-            element_found, sigma0, sigma1, sigma2, k = \
-                           search_tree_of_vertices(self.root, self.mesh, x)
-            #Update interpolation matrix A if necessary
-            if element_found is True:
-                #Assign values to matrix A
-
-                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
-                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
-                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
-
-                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
-                js     = [j0,j1,j2]
-
-                for j in js:
-                    A[i,j] = sigmas[j]
-            else:
-                msg = 'Could not find triangle for point', x 
-                raise Exception(msg)
-        return A
-
 
      # FIXME: What is a good start_blocking_len value?
@@ -278 +210 @@
 
 
+    def _build_interpolation_matrix_A(self,
+                                     point_coordinates,
+                                     verbose = False):
+        """Build n x m interpolation matrix, where
+        n is the number of data points and
+        m is the number of basis functions phi_k (one per vertex)
+
+        This algorithm uses a quad tree data structure for fast binning
+        of data points
+        origin is a 3-tuple consisting of UTM zone, easting and northing.
+        If specified coordinates are assumed to be relative to this origin.
+
+        This one will override any data_origin that may be specified in
+        instance interpolation
+
+        Preconditions
+        Point_coordindates and mesh vertices have the same origin.
+        """
+
+
+        
+        #Convert point_coordinates to Numeric arrays, in case it was a list.
+        point_coordinates = ensure_numeric(point_coordinates, Float)
+        
+        if verbose: print 'Getting indices inside mesh boundary'
+        #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() 
+        self.inside_poly_indices, self.outside_poly_indices  = \
+                     in_and_outside_polygon(point_coordinates,
+                                            self.mesh.get_boundary_polygon(),
+                                            closed = True, verbose = verbose)
+        #print "self.inside_poly_indices",self.inside_poly_indices
+        #print "self.outside_poly_indices",self.outside_poly_indices 
+        #Build n x m interpolation matrix
+        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        n = point_coordinates.shape[0]     #Nbr of data points
+
+        if verbose: print 'Number of datapoints: %d' %n
+        if verbose: print 'Number of basis functions: %d' %m
+
+        A = Sparse(n,m)
+
+        #Compute matrix elements for points inside the mesh
+        for i in self.inside_poly_indices:
+            #For each data_coordinate point
+            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
+            x = point_coordinates[i]
+            element_found, sigma0, sigma1, sigma2, k = \
+                           search_tree_of_vertices(self.root, self.mesh, x)
+            #Update interpolation matrix A if necessary
+            if element_found is True:
+                #Assign values to matrix A
+
+                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
+                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
+                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
+
+                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
+                js     = [j0,j1,j2]
+
+                for j in js:
+                    A[i,j] = sigmas[j]
+            else:
+                msg = 'Could not find triangle for point', x 
+                raise Exception(msg)
+        return A
+
 class Interpolation_function:
     """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
@@ -499 +497 @@
         #Compute interpolated values
         q = zeros(len(self.quantity_names), Float)
-
+        #print "self.precomputed_values", self.precomputed_values
         for i, name in enumerate(self.quantity_names):
             Q = self.precomputed_values[name]

inundation/fit_interpolate/test_interpolate.py

@@ -798 +798 @@
 
 
+    def fix_test_interpolation_precompute_points(self):
+        # looking at a discrete mesh
+        #
+    
+        #Three timesteps
+        time = [0.0, 60.0]    
+
+        #Setup mesh used to represent fitted function
+        points = [[ 15., -20.],
+                  [ 15.,  10.],
+                  [  0., -20.],
+                  [  0.,  10.],
+                  [  0., -20.],
+                  [ 15.,  10.]]
+        
+        triangles = [[0, 1, 2],
+                     [3, 4, 5]]
+
+        #New datapoints where interpolated values are sought
+        interpolation_points = [[ 1.,  0.], [0.,1.]]
+
+        #One quantity
+        Q = zeros( (2,6), Float )
+
+        #Linear in time and space
+        for i, t in enumerate(time):
+            Q[i, :] = t*linear_function(points)
+        #print "Q", Q
+
+
+        
+        interp = Interpolate(points, triangles)
+        f = array([linear_function(points),2*linear_function(points) ])
+        f = transpose(f)
+        #print "f",f
+        z = interp.interpolate(f, interpolation_points)
+        answer = [linear_function(interpolation_points),
+                  2*linear_function(interpolation_points) ]
+        answer = transpose(answer)
+        print "z",z
+        print "answer",answer
+        assert allclose(z, answer)
+
+
+        #Check interpolation of one quantity using interpolaton points)
+        I = Interpolation_function(time, Q,
+                                   vertex_coordinates = points,
+                                   triangles = triangles, 
+                                   interpolation_points = interpolation_points,
+                                   verbose = False)
+        print "I.precomputed_values", I.precomputed_values
+        
+        self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
+                        ' failed')
+        
     def test_interpolation_function_outside_point(self):
         # Test spatio-temporal interpolation
@@ -911 +966 @@
 
     suite = unittest.makeSuite(Test_Interpolate,'test')
-    #suite = unittest.makeSuite(Test_Interpolate,'test_points_outside_the_polygon')
+    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)