Ignore:
Timestamp:
Jan 11, 2006, 5:32:23 PM (19 years ago)
Author:
duncan
Message:

comments

Location:
inundation/fit_interpolate
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/fit_interpolate/interpolate.py

    r2192 r2201  
    1515* remove points outside the mesh ?(in interpolate_block)?
    1616* geo-ref (in interpolate_block)
     17* add functional interpolate interface - in mesh and points, out interp data
    1718"""
    1819
     
    3738                 vertex_coordinates,
    3839                 triangles,
    39                  #point_coordinates = None,
    4040                 mesh_origin = None,
    4141                 verbose=False,
    42                  max_points_per_cell=30):
     42                 max_vertices_per_cell=30):
    4343
    4444
     
    4949
    5050          vertex_coordinates: List of coordinate pairs [xi, eta] of
    51           points constituting mesh (or a an m x 2 Numeric array)
    52           Points may appear multiple times
    53           (e.g. if vertices have discontinuities)
     51              points constituting a mesh (or an m x 2 Numeric array)
     52              Points may appear multiple times
     53              (e.g. if vertices have discontinuities)
    5454
    5555          triangles: List of 3-tuples (or a Numeric array) of
    56           integers representing indices of all vertices in the mesh.
    57 
    58           point_coordinates: List of coordinate pairs [x, y] of
    59           data points (or an nx2 Numeric array)
    60           If point_coordinates is absent, only smoothing matrix will
    61           be built
    62 
    63           alpha: Smoothing parameter
    64 
    65           data_origin and mesh_origin are 3-tuples consisting of
    66           UTM zone, easting and northing. If specified
    67           point coordinates and vertex coordinates are assumed to be
    68           relative to their respective origins.
    69 
    70         """
    71        
    72         from pyvolution.util import ensure_numeric
    73 
     56              integers representing indices of all vertices in the mesh.
     57
     58          mesh_origin: 3-tuples consisting of
     59              UTM zone, easting and northing.
     60              If specified vertex coordinates are assumed to be
     61              relative to their respective origins.
     62
     63          max_vertices_per_cell: Number of vertices in a quad tree cell
     64          at which the cell is split into 4.
     65
     66        """
     67       
    7468        # Initialise variabels
    75         self.A_can_be_reused = False
    76         self.point_coordinates = None
     69        self._A_can_be_reused = False
     70        self._point_coordinates = None
    7771       
    7872        #Convert input to Numeric arrays
     
    10296
    10397        self.root = build_quadtree(self.mesh,
    104                               max_points_per_cell = max_points_per_cell)
     98                              max_points_per_cell = max_vertices_per_cell)
    10599       
    106100       
     
    148142            #For each data_coordinate point
    149143            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
    150 
    151144            x = point_coordinates[i]
    152            
    153145            element_found, sigma0, sigma1, sigma2, k = \
    154146                           search_tree_of_vertices(self.root, self.mesh, x)
     
    168160            else:
    169161                print 'Could not find triangle for point', x
    170                
    171162        return A
    172163
     
    281272              per underlying matrix-matrix multiplication
    282273          point_coordinates: Interpolate mesh data to these positions.
     274              List of coordinate pairs [x, y] of
     275              data points (or an nx2 Numeric array)
     276              If point_coordinates is absent, the points inputted last time
     277              this method was called are used, if possible.
    283278          start_blocking_len: If the # of points is more or greater than this,
    284279              start blocking
     
    290285        # Can I interpolate, based on previous point_coordinates?
    291286        if point_coordinates is None:
    292             if self.A_can_be_reused is True and \
    293             len(self.point_coordinates) < start_blocking_len:
     287            if self._A_can_be_reused is True and \
     288            len(self._point_coordinates) < start_blocking_len:
    294289                z = self._get_point_data_z(f,
    295290                                           verbose=verbose)
    296             elif self.point_coordinates is not None:
     291            elif self._point_coordinates is not None:
    297292                #     if verbose, give warning
    298293                if verbose:
    299294                    print 'WARNING: Recalculating A matrix, due to blocking.'
    300                 point_coordinates = self.point_coordinates
     295                point_coordinates = self._point_coordinates
    301296            else:
    302297                #There are no good point_coordinates. import sys; sys.exit()
     
    305300           
    306301        if point_coordinates is not None:
    307             self.point_coordinates = point_coordinates
     302            self._point_coordinates = point_coordinates
    308303            if len(point_coordinates) < start_blocking_len or \
    309304                   start_blocking_len == 0:
    310                 self.A_can_be_reused = True
     305                self._A_can_be_reused = True
    311306                z = self.interpolate_block(f, point_coordinates,
    312307                                           verbose=verbose)
    313308            else:
    314309                #Handle blocking
    315                 self.A_can_be_reused = False
     310                self._A_can_be_reused = False
    316311                start=0
    317312                z = self.interpolate_block(f, point_coordinates[0:0])
     
    337332        """
    338333        if point_coordinates is not None:
    339             self.A =self._build_interpolation_matrix_A(point_coordinates,
     334            self._A =self._build_interpolation_matrix_A(point_coordinates,
    340335                                                       verbose=verbose)
    341336        return self._get_point_data_z(f)
    342337
    343338    def _get_point_data_z(self, f, verbose=False):
    344         return self.A * f
     339        return self._A * f
    345340       
    346341#-------------------------------------------------------------
  • inundation/fit_interpolate/search_functions.py

    r2190 r2201  
    77"""
    88from Numeric import dot
    9 
    10 
    11 
    12 import time
    13 
    14 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
    15      ArrayType, allclose, take
    16 
    17 from pyvolution.mesh import Mesh
    18 from pyvolution.sparse import Sparse, Sparse_CSR
    19 from pyvolution.cg_solve import conjugate_gradient, VectorShapeError
    20 from coordinate_transforms.geo_reference import Geo_reference
    21 from pyvolution.quad import build_quadtree
    22 from utilities.numerical_tools import ensure_numeric
    23 from utilities.polygon import inside_polygon
    249
    2510
  • inundation/fit_interpolate/test_interpolate.py

    r2190 r2201  
    6767        data = [ [4,4] ]
    6868        interp = Interpolate(points, triangles,
    69                                max_points_per_cell = 4)
     69                               max_vertices_per_cell = 4)
    7070        #print "PDSG - interp.get_A()", interp.get_A()
    7171        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
     
    507507        #print "answer",answer
    508508        assert allclose(z, answer)
    509         assert allclose(interp.A_can_be_reused, True)
     509        assert allclose(interp._A_can_be_reused, True)
    510510
    511511        z = interp.interpolate(f)
     
    515515        z = interp.interpolate(f, start_blocking_len=10)
    516516        assert allclose(z, answer)
    517         assert allclose(interp.A_can_be_reused, False)
     517        assert allclose(interp._A_can_be_reused, False)
    518518
    519519        #A is recalculated
    520520        z = interp.interpolate(f)
    521521        assert allclose(z, answer)
    522         assert allclose(interp.A_can_be_reused, True)
    523 
     522        assert allclose(interp._A_can_be_reused, True)
     523       
    524524        interp = Interpolate(vertices, triangles)
    525525        #Must raise an exception, no points specified
Note: See TracChangeset for help on using the changeset viewer.