Changeset 6184
 Timestamp:
 Jan 17, 2009, 2:48:28 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/fit_interpolate/interpolate.py
r6173 r6184 38 38 from anuga.abstract_2d_finite_volumes.util import file_function 39 39 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 40 from utilities.polygon import point_on_line 41 40 42 41 43 import Numeric as num … … 246 248 Interpolated values at inputted points (z). 247 249 """ 250 251 #print f 252 #print vertex_coordinates 253 #print gauge_neighbour_id 254 #print point_coordinates 255 248 256 249 257 # FIXME: There is an option of passing a tolerance into this … … 252 260 point_coordinates = point_coordinates.get_data_points(absolute=True) 253 261 254 from utilities.polygon import point_on_line255 256 262 z = num.ones(len(point_coordinates), num.Float) 257 263 … … 274 280 275 281 elif n > 1: 276 for i in range(len(point_coordinates)): 277 found = False 278 for j in range(n): 279 if gauge_neighbour_id[j] >= 0: 280 if point_on_line(point_coordinates[i], 281 [vertex_coordinates[j], 282 vertex_coordinates[gauge_neighbour_id[j]]], 283 rtol=1.0e6): 284 found = True 285 x0 = vertex_coordinates[j][0] 286 y0 = vertex_coordinates[j][1] 287 x1 = vertex_coordinates[gauge_neighbour_id[j]][0] 288 y1 = vertex_coordinates[gauge_neighbour_id[j]][1] 289 x2 = point_coordinates[i][0] 290 y2 = point_coordinates[i][1] 291 292 segment_len = sqrt((x1x0)**2 + (y1y0)**2) 293 dist = sqrt((x2x0)**2 + (y2y0)**2) 294 z[i] = (f[gauge_neighbour_id[j]]  f[j]) \ 295 / segment_len*dist + f[j] 296 break 297 298 if not found: 299 z[i] = 0.0 282 _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id) 283 300 284 return z 301 285 286 302 287 303 288 ## … … 361 346 point_coordinates = self._point_coordinates 362 347 else: 363 # There are no good point_coordinates. import sys; sys.exit()348 # There are no good point_coordinates. import sys; sys.exit() 364 349 msg = 'ERROR (interpolate.py): No point_coordinates inputted' 365 350 raise Exception(msg) … … 373 358 verbose=verbose) 374 359 else: 375 # Handle blocking360 # Handle blocking 376 361 self._A_can_be_reused = False 377 362 start = 0 … … 588 573 589 574 575 576 577 578 579 def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id): 580 """Auxiliary function 581 """ 582 for i in range(len(point_coordinates)): 583 584 x2 = point_coordinates[i][0] 585 y2 = point_coordinates[i][1] 586 587 found = False 588 for j in range(n): 589 590 neighbour_id = gauge_neighbour_id[j] 591 if neighbour_id >= 0: 592 x0 = vertex_coordinates[j][0] 593 y0 = vertex_coordinates[j][1] 594 x1 = vertex_coordinates[neighbour_id][0] 595 y1 = vertex_coordinates[neighbour_id][1] 596 597 if point_on_line([x2,y2], 598 [[x0, y0], [x1, y1]], 599 rtol=1.0e6): 600 601 found = True 602 segment_len = sqrt((x1x0)**2 + (y1y0)**2) 603 dist = sqrt((x2x0)**2 + (y2y0)**2) 604 z[i] = (f[neighbour_id]  f[j])*dist/segment_len + f[j] 605 break 606 607 if not found: 608 z[i] = 0.0 609 610 590 611 ## 591 612 # @brief ??
Note: See TracChangeset
for help on using the changeset viewer.