Changeset 6186
 Timestamp:
 Jan 17, 2009, 3:24:28 PM (10 years ago)
 Location:
 anuga_core/source/anuga/fit_interpolate
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/fit_interpolate/interpolate.py
r6185 r6186 516 516 # @param verbose True if this function is to be verbose. 517 517 def interpolate_polyline(f, 518 vertex_coordinates,518 polyline_nodes, 519 519 gauge_neighbour_id, 520 520 point_coordinates=None, 521 521 verbose=False): 522 """Interpolate linearly between values f on nodes (vertex coordinates)522 """Interpolate linearly between values f on polyline nodes 523 523 of a polyline to midpoints of triangles of boundary. 524 524 … … 539 539 Interpolated values at inputted points (z). 540 540 """ 541 542 #print f543 #print vertex_coordinates544 #print gauge_neighbour_id545 #print point_coordinates546 547 548 541 # FIXME: There is an option of passing a tolerance into this 549 542 … … 551 544 point_coordinates = point_coordinates.get_data_points(absolute=True) 552 545 553 z = num.ones(len(point_coordinates), num.Float) 554 555 # input sanity check 546 z = num.zeros(len(point_coordinates), num.Float) 547 548 f = ensure_numeric(f, num.Float) 549 polyline_nodes = ensure_numeric(polyline_nodes, num.Float) 550 point_coordinates = ensure_numeric(point_coordinates, num.Float) 551 gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.Int) 552 553 n = polyline_nodes.shape[0] # Number of nodes in polyline 554 # Input sanity check 556 555 msg = 'point coordinates are not given (interpolate.py)' 557 556 assert point_coordinates is not None, msg 558 557 msg = 'function value must be specified at every interpolation node' 559 assert f.shape[0]== vertex_coordinates.shape[0], msg558 assert f.shape[0]==polyline_nodes.shape[0], msg 560 559 msg = 'Must define function value at one or more nodes' 561 560 assert f.shape[0]>0, msg 562 561 563 n = f.shape[0] 562 564 563 if n == 1: 565 z = f*z 566 msg = 'Polyline contained only one point. I need more. ', str(f) 564 msg = 'Polyline contained only one point. I need more. ' + str(f) 567 565 raise Exception, msg 568 569 # FIXME (John): add unit test for only One vertex point.570 # Exception should be thrown.571 572 566 elif n > 1: 573 _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id) 567 _interpolate_polyline_aux(point_coordinates, len(point_coordinates), 568 polyline_nodes, len(polyline_nodes), 569 f, 570 gauge_neighbour_id, 571 z) 574 572 575 573 return z 576 574 577 575 578 def _interpolate_polyline_aux(n, z, f, point_coordinates, vertex_coordinates, gauge_neighbour_id): 579 """Auxiliary function 576 def _interpolate_polyline_aux(point_coordinates, number_of_points, 577 polyline_nodes, number_of_nodes, 578 f, gauge_neighbour_id, z): 579 """Auxiliary function used by interpolate_polyline 580 580 """ 581 for i in range(len(point_coordinates)): 582 583 x2 = point_coordinates[i][0] 584 y2 = point_coordinates[i][1] 585 581 582 for i in range(number_of_points): 583 584 x2, y2 = point_coordinates[i,:] 586 585 found = False 587 for j in range(n): 586 587 for j in range(number_of_nodes): 588 588 589 589 neighbour_id = gauge_neighbour_id[j] 590 590 if neighbour_id >= 0: 591 x0 = vertex_coordinates[j][0] 592 y0 = vertex_coordinates[j][1] 593 x1 = vertex_coordinates[neighbour_id][0] 594 y1 = vertex_coordinates[neighbour_id][1] 591 x0, y0 = polyline_nodes[j,:] 592 x1, y1 = polyline_nodes[neighbour_id,:] 595 593 596 if point_on_line([x2, y2],594 if point_on_line([x2, y2], 597 595 [[x0, y0], [x1, y1]], 598 596 rtol=1.0e6): 
anuga_core/source/anuga/fit_interpolate/test_interpolate.py
r6185 r6186 1896 1896 #print z 1897 1897 assert num.allclose(z, z_ref) 1898 1899 # Test exception thrown for one point 1900 f = num.array([5]) 1901 vertex_coordinates = num.array([[4., 4.]]) 1902 try: 1903 z = interpolate_polyline(f, vertex_coordinates, gauge_neighbour_id, point_coordinates) 1904 except Exception: 1905 pass 1906 else: 1907 raise Exception, 'One point should have raised exception' 1898 1908 1899 1909 #
Note: See TracChangeset
for help on using the changeset viewer.