Changeset 5403 for anuga_core/source/anuga/fit_interpolate/interpolate.py
- Timestamp:
- Jun 13, 2008, 3:48:16 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/interpolate.py
r5372 r5403 96 96 point_coordinates=None, 97 97 verbose=False): 98 98 """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles 99 of boundary. 100 101 f is the data on the polyline nodes. 102 103 The mesh values representing a smooth surface are 104 assumed to be specified in f. 105 106 Inputs: 107 f: Vector or array of data at the polyline nodes. 108 If f is an array, interpolation will be done for each column as 109 per underlying matrix-matrix multiplication 110 point_coordinates: Interpolate polyline data to these positions. 111 List of coordinate pairs [x, y] of 112 data points or an nx2 Numeric array or a Geospatial_data object 113 114 Output: 115 Interpolated values at inputted points (z). 116 """ 117 118 119 # FIXME: There is an option of passing a tolerance into this 120 99 121 if isinstance(point_coordinates, Geospatial_data): 100 122 point_coordinates = point_coordinates.get_data_points( \ 101 123 absolute = True) 102 124 103 from utilities.polygon import point_on_line ,point_on_line_py125 from utilities.polygon import point_on_line 104 126 from Numeric import ones 105 127 z=ones(len(point_coordinates),Float) … … 108 130 assert point_coordinates is not None, msg 109 131 msg='function value must be specified at every interpolation node' 110 assert f.shape[0]==vertex_coordinates.shape[0], msg132 assert f.shape[0]==vertex_coordinates.shape[0], msg 111 133 msg='Must define function value at one or more nodes' 112 assert f.shape[0]>0, msg134 assert f.shape[0]>0, msg 113 135 114 136 n=f.shape[0] 115 137 if n==1: 116 138 z=f*z 139 msg = 'Polyline contained only one point. I need more. ', str(f) 140 raise Exception, msg 141 142 # FIXME (John): add unit test for only One vertex point. Exception should be thrown. 143 144 117 145 elif n>1: 118 146 for i in range(len(point_coordinates)): … … 120 148 for j in range(n): 121 149 if gauge_neighbour_id[j]>=0: 122 if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]): 150 if point_on_line(point_coordinates[i], 151 [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]], 152 rtol=1.0e-6): 123 153 found=True 124 154 x0=vertex_coordinates[j][0]
Note: See TracChangeset
for help on using the changeset viewer.