Changeset 6185


Ignore:
Timestamp:
Jan 17, 2009, 2:56:59 PM (10 years ago)
Author:
ole
Message:

Moved interpolate_polyline out of Interpolation class - it had nothing to do with it anyway.

Location:
anuga_core/source/anuga/fit_interpolate
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r6184 r6185  
    214214        self._point_coordinates = None # FIXME (Ole): Probably obsolete
    215215        self.interpolation_matrices = {} # Store precomputed matrices
    216 
    217 
    218     ##
    219     # @brief Interpolate linearly from polyline nodes to midpoints of triangles.
    220     # @param f The data on the polyline nodes.
    221     # @param vertex_coordinates ??
    222     # @param gauge_neighbour_id ??
    223     # @param point_coordinates ??
    224     # @param verbose True if this function is to be verbose.
    225     def interpolate_polyline(self,
    226                              f,
    227                              vertex_coordinates,
    228                              gauge_neighbour_id,
    229                              point_coordinates=None,
    230                              verbose=False):
    231         """Interpolate linearly between values f on nodes (vertex coordinates)
    232         of a polyline to midpoints of triangles of boundary.
    233 
    234         f is the data on the polyline nodes.
    235 
    236         The mesh values representing a smooth surface are
    237         assumed to be specified in f.
    238 
    239         Inputs:
    240           f: Vector or array of data at the polyline nodes.
    241               If f is an array, interpolation will be done for each column as
    242               per underlying matrix-matrix multiplication
    243           point_coordinates: Interpolate polyline data to these positions.
    244               List of coordinate pairs [x, y] of
    245               data points or an nx2 Numeric array or a Geospatial_data object
    246              
    247         Output:
    248           Interpolated values at inputted points (z).
    249         """
    250        
    251         #print f
    252         #print vertex_coordinates
    253         #print gauge_neighbour_id
    254         #print point_coordinates
    255 
    256 
    257         # FIXME: There is an option of passing a tolerance into this
    258        
    259         if isinstance(point_coordinates, Geospatial_data):
    260             point_coordinates = point_coordinates.get_data_points(absolute=True)
    261  
    262         z = num.ones(len(point_coordinates), num.Float)
    263 
    264         # input sanity check
    265         msg = 'point coordinates are not given (interpolate.py)'
    266         assert point_coordinates is not None, msg
    267         msg = 'function value must be specified at every interpolation node'
    268         assert f.shape[0]==vertex_coordinates.shape[0], msg
    269         msg = 'Must define function value at one or more nodes'
    270         assert f.shape[0]>0, msg
    271 
    272         n = f.shape[0]
    273         if n == 1:
    274             z = f*z
    275             msg = 'Polyline contained only one point. I need more. ', str(f)
    276             raise Exception, msg
    277            
    278         # FIXME (John): add unit test for only One vertex point.
    279         #               Exception should be thrown.
    280        
    281         elif n > 1:
    282             _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
    283            
    284         return z
    285216
    286217   
     
    576507       
    577508       
     509
     510##
     511# @brief Interpolate linearly from polyline nodes to midpoints of triangles.
     512# @param f The data on the polyline nodes.
     513# @param vertex_coordinates ??
     514# @param gauge_neighbour_id ??
     515# @param point_coordinates ??
     516# @param verbose True if this function is to be verbose.
     517def interpolate_polyline(f,
     518                         vertex_coordinates,
     519                         gauge_neighbour_id,
     520                         point_coordinates=None,
     521                         verbose=False):
     522    """Interpolate linearly between values f on nodes (vertex coordinates)
     523    of a polyline to midpoints of triangles of boundary.
     524
     525    f is the data on the polyline nodes.
     526
     527    The mesh values representing a smooth surface are
     528    assumed to be specified in f.
     529
     530    Inputs:
     531      f: Vector or array of data at the polyline nodes.
     532          If f is an array, interpolation will be done for each column as
     533          per underlying matrix-matrix multiplication
     534      point_coordinates: Interpolate polyline data to these positions.
     535          List of coordinate pairs [x, y] of
     536          data points or an nx2 Numeric array or a Geospatial_data object
     537         
     538    Output:
     539      Interpolated values at inputted points (z).
     540    """
     541   
     542    #print f
     543    #print vertex_coordinates
     544    #print gauge_neighbour_id
     545    #print point_coordinates
     546
     547
     548    # FIXME: There is an option of passing a tolerance into this
     549   
     550    if isinstance(point_coordinates, Geospatial_data):
     551        point_coordinates = point_coordinates.get_data_points(absolute=True)
     552
     553    z = num.ones(len(point_coordinates), num.Float)
     554
     555    # input sanity check
     556    msg = 'point coordinates are not given (interpolate.py)'
     557    assert point_coordinates is not None, msg
     558    msg = 'function value must be specified at every interpolation node'
     559    assert f.shape[0]==vertex_coordinates.shape[0], msg
     560    msg = 'Must define function value at one or more nodes'
     561    assert f.shape[0]>0, msg
     562
     563    n = f.shape[0]
     564    if n == 1:
     565        z = f*z
     566        msg = 'Polyline contained only one point. I need more. ', str(f)
     567        raise Exception, msg
     568       
     569    # FIXME (John): add unit test for only One vertex point.
     570    #               Exception should be thrown.
     571   
     572    elif n > 1:
     573        _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id)
     574       
     575    return z
     576
    578577       
    579578def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id):
     
    10511050                                                      verbose=False) # No clutter
    10521051                    elif triangles is None and vertex_coordinates is not None:
    1053                         result = \
    1054                             interpol.interpolate_polyline(Q,
    1055                                                           vertex_coordinates,
    1056                                                           gauge_neighbour_id,
    1057                                                           point_coordinates=\
     1052                        result = interpolate_polyline(Q,
     1053                                                      vertex_coordinates,
     1054                                                      gauge_neighbour_id,
     1055                                                      point_coordinates=\
    10581056                                                          self.interpolation_points)
    1059 
     1057                       
    10601058                    #assert len(result), len(interpolation_points)
    10611059                    self.precomputed_values[name][i, :] = result
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r6183 r6185  
    18531853        John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately.
    18541854        """
    1855 
    1856         a = [-1.0, 0.0]
    1857         b = [3.0, 4.0]
    1858         c = [4.0, 1.0]
    1859         d = [-3.0, 2.0] #3
    1860         e = [-1.0, -2.0]
    1861         f = [1.0, -2.0] #5
    1862 
    1863         vertices = [a, b, c, d,e,f]
    1864         triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
    1865 
    1866 
    1867         point_coords = [[-2.0, 2.0],
    1868                         [-1.0, 1.0],
    1869                         [0.0, 2.0],
    1870                         [1.0, 1.0],
    1871                         [2.0, 1.0],
    1872                         [0.0, 0.0],
    1873                         [1.0, 0.0],
    1874                         [0.0, -1.0],
    1875                         [-0.2, -0.5],
    1876                         [-0.9, -1.5],
    1877                         [0.5, -1.9],
    1878                         [3.0, 1.0]]
    1879 
    1880         interp = Interpolate(vertices, triangles)
    1881                
    18821855       
    18831856        f = num.array([58.06150614, 58.06150614, 58.06150614])
     
    18961869        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
    18971870       
    1898         z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1871        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    18991872        assert num.allclose(z, z_ref)
    19001873       
     
    19021875        f = num.array([58.06150614, 158.06150614, 258.06150614])
    19031876        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]       
    1904         z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1877        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    19051878        assert num.allclose(z, z_ref)       
    19061879                       
     
    19191892        gauge_neighbour_id = [1, 2, -1]
    19201893                                               
    1921         z = interp.interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
     1894        z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates)
    19221895        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
    19231896        #print z
Note: See TracChangeset for help on using the changeset viewer.