Changeset 5403
 Timestamp:
 Jun 13, 2008, 3:48:16 PM (14 years ago)
 Location:
 anuga_core/source/anuga
 Files:

 5 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 matrixmatrix 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.0e6): 123 153 found=True 124 154 x0=vertex_coordinates[j][0] 
anuga_core/source/anuga/shallow_water/data_manager.py
r5391 r5403 5413 5413 #Title 5414 5414 fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \ 5415 5415 'Northing' +d+ 'depth m' + "\n") 5416 5416 i = 0 5417 5417 for depth, point_utm, lat, long in map(None, depths, 
anuga_core/source/anuga/utilities/polygon.py
r5372 r5403 17 17 from anuga.geospatial_data.geospatial_data import ensure_absolute 18 18 19 def point_on_line_py(point, line): 20 from Numeric import fabs 21 point = ensure_numeric(point) 22 line = ensure_numeric(line) 23 24 25 x=point[0];y=point[1] 26 x0=line[0,0];y0=line[0,1] 27 x1=line[1,0];y1=line[1,1] 28 #from pylab import plot,show 29 #plot(line[:,0],line[:,1]) 30 #plot([x],[y],'o') 31 #show() 32 a0 = x  x0; 33 a1 = y  y0; 34 35 a_normal0 = a1; 36 a_normal1 = a0; 37 38 b0 = x1  x0; 39 b1 = y1  y0; 40 41 #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ): 42 eps=200 43 #print 'normal',a_normal0*b0 + a_normal1*b1 44 if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ): 45 #for some reason (perhaps catastrophic cancellation) urs_boundary.py 46 #example only works for eps=2 47 #point is somewhere on the infinite extension of the line 48 # FIXME (Ole): Perhaps add a tolerance here instead of 0.0 49 50 len_a = sqrt(a0*a0 + a1*a1); 51 len_b = sqrt(b0*b0 + b1*b1); 52 53 if (a0*b0 + a1*b1 >= 0 and len_a <= len_b): 54 return 1 55 else: 56 return 0 57 else: 58 return 0 59 60 61 def point_on_line(point, line): 19 20 def point_on_line(point, line, rtol=0.0, atol=0.0): 62 21 """Determine whether a point is on a line segment 63 22 … … 68 27 69 28 Output: 70 71 """ 29 30 Note: Line can be degenerate and function still works to discern coinciding points from noncoinciding. 31 """ 32 33 # FIXME(Ole): Perhaps make defaults as in allclose: rtol=1.0e5, atol=1.0e8 72 34 73 35 point = ensure_numeric(point) 74 36 line = ensure_numeric(line) 75 37 76 77 38 res = _point_on_line(point[0], point[1], 78 39 line[0,0], line[0,1], 79 line[1,0], line[1,1]) 40 line[1,0], line[1,1], 41 rtol, atol) 80 42 81 43 return bool(res) 
anuga_core/source/anuga/utilities/polygon_ext.c
r5321 r5403 21 21 int __point_on_line(double x, double y, 22 22 double x0, double y0, 23 double x1, double y1) { 23 double x1, double y1, 24 double rtol, 25 double atol) { 24 26 /*Determine whether a point is on a line segment 25 27 … … 31 33 32 34 double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b; 35 double nominator, denominator; 36 int is_parallel; 33 37 34 38 a0 = x  x0; … … 41 45 b1 = y1  y0; 42 46 43 if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) { 44 //Point is somewhere on the infinite extension of the line 45 // FIXME (Ole): Perhaps add a tolerance here instead of 0.0 47 nominator = fabs(a_normal0*b0 + a_normal1*b1); 48 denominator = b0*b0 + b1*b1; 49 50 // Determine if line is parallel to point vector up to a tolerance 51 is_parallel = 0; 52 if (denominator == 0.0) { 53 // Use absolute tolerance 54 if (nominator <= atol) { 55 is_parallel = 1; 56 } 57 } else { 58 // Denominator is positive  use relative tolerance 59 if (nominator/denominator <= rtol) { 60 is_parallel = 1; 61 } 62 } 63 64 if (is_parallel) { 65 // Point is somewhere on the infinite extension of the line 66 // subject to specified absolute tolerance 46 67 47 68 len_a = sqrt(a0*a0 + a1*a1); … … 173 194 int verbose) { 174 195 175 double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j ;196 double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0; 176 197 int i, j, k, outside_index, inside_index, inside; 177 198 … … 223 244 224 245 //Check for case where point is contained in line segment 225 if (__point_on_line(x, y, px_i, py_i, px_j, py_j )) {246 if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) { 226 247 if (closed == 1) { 227 248 inside = 1; … … 260 281 // 261 282 262 double x, y, x0, y0, x1, y1 ;283 double x, y, x0, y0, x1, y1, rtol, atol; 263 284 int res; 264 285 PyObject *result; 265 286 266 287 // Convert Python arguments to C 267 if (!PyArg_ParseTuple(args, "dddddd ", &x, &y, &x0, &y0, &x1, &y1)) {288 if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) { 268 289 PyErr_SetString(PyExc_RuntimeError, 269 290 "point_on_line could not parse input"); … … 273 294 274 295 // Call underlying routine 275 res = __point_on_line(x, y, x0, y0, x1, y1 );296 res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol); 276 297 277 298 // Return values a and b 
anuga_core/source/anuga/utilities/test_polygon.py
r5288 r5403 165 165 assert not point_on_line( [40,19], [[40,20], [40,40]] ) 166 166 167 # Degenerate line 168 assert point_on_line( [40,19], [[40,19], [40,19]] ) 169 assert not point_on_line( [40,19], [[40,40], [40,40]] ) 167 170 168 171
Note: See TracChangeset
for help on using the changeset viewer.