Changeset 6185
 Timestamp:
 Jan 17, 2009, 2:56:59 PM (14 years ago)
 Location:
 anuga_core/source/anuga/fit_interpolate
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/fit_interpolate/interpolate.py
r6184 r6185 214 214 self._point_coordinates = None # FIXME (Ole): Probably obsolete 215 215 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 are237 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 as242 per underlying matrixmatrix multiplication243 point_coordinates: Interpolate polyline data to these positions.244 List of coordinate pairs [x, y] of245 data points or an nx2 Numeric array or a Geospatial_data object246 247 Output:248 Interpolated values at inputted points (z).249 """250 251 #print f252 #print vertex_coordinates253 #print gauge_neighbour_id254 #print point_coordinates255 256 257 # FIXME: There is an option of passing a tolerance into this258 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 check265 msg = 'point coordinates are not given (interpolate.py)'266 assert point_coordinates is not None, msg267 msg = 'function value must be specified at every interpolation node'268 assert f.shape[0]==vertex_coordinates.shape[0], msg269 msg = 'Must define function value at one or more nodes'270 assert f.shape[0]>0, msg271 272 n = f.shape[0]273 if n == 1:274 z = f*z275 msg = 'Polyline contained only one point. I need more. ', str(f)276 raise Exception, msg277 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 z285 216 286 217 … … 576 507 577 508 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. 517 def 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 matrixmatrix 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 578 577 579 578 def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id): … … 1051 1050 verbose=False) # No clutter 1052 1051 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=\ 1058 1056 self.interpolation_points) 1059 1057 1060 1058 #assert len(result), len(interpolation_points) 1061 1059 self.precomputed_values[name][i, :] = result 
anuga_core/source/anuga/fit_interpolate/test_interpolate.py
r6183 r6185 1853 1853 John Jakeman works. It has been exercised somewhat by tests of sts boundary, but never before separately. 1854 1854 """ 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] #31860 e = [1.0, 2.0]1861 f = [1.0, 2.0] #51862 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 afc1865 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 1882 1855 1883 1856 f = num.array([58.06150614, 58.06150614, 58.06150614]) … … 1896 1869 z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614] 1897 1870 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) 1899 1872 assert num.allclose(z, z_ref) 1900 1873 … … 1902 1875 f = num.array([58.06150614, 158.06150614, 258.06150614]) 1903 1876 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) 1905 1878 assert num.allclose(z, z_ref) 1906 1879 … … 1919 1892 gauge_neighbour_id = [1, 2, 1] 1920 1893 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) 1922 1895 z_ref = [1.1, 4.5, 5., 7.4, 11., 0.] 1923 1896 #print z
Note: See TracChangeset
for help on using the changeset viewer.