Changeset 5372
- Timestamp:
- May 28, 2008, 2:47:30 PM (16 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r5361 r5372 15 15 16 16 from anuga.utilities.numerical_tools import ensure_numeric 17 from Numeric import arange, choose, zeros, Float, array 17 from Numeric import arange, choose, zeros, Float, array, allclose, take, compress 18 18 19 19 from anuga.geospatial_data.geospatial_data import ensure_absolute … … 33 33 time_thinning=1, 34 34 verbose=False, 35 use_cache=False): 35 use_cache=False, 36 boundary_polygon=None): 36 37 """Read time history of spatial data from NetCDF file and return 37 38 a callable object. … … 114 115 'domain_starttime': domain_starttime, 115 116 'time_thinning': time_thinning, 116 'verbose': verbose} 117 'verbose': verbose, 118 'boundary_polygon':boundary_polygon} 117 119 118 120 … … 130 132 dependencies=[filename], 131 133 compression=False, 132 verbose=verbose) 134 verbose=verbose, 135 boundary_polygon=boundary_polygon) 133 136 134 137 else: … … 166 169 domain_starttime=None, 167 170 time_thinning=1, 168 verbose=False): 171 verbose=False, 172 boundary_polygon=None): 169 173 """Internal function 170 174 … … 193 197 domain_starttime, 194 198 time_thinning=time_thinning, 195 verbose=verbose) 199 verbose=verbose, 200 boundary_polygon=boundary_polygon) 196 201 else: 197 202 raise 'Must be a NetCDF File' … … 204 209 domain_starttime=None, 205 210 time_thinning=1, 206 verbose=False): 211 verbose=False, 212 boundary_polygon=None): 207 213 """Read time history of spatial data from NetCDF sww file and 208 214 return a callable object f(t,x,y) … … 282 288 raise msg 283 289 290 if filename[-3:] == 'sts' and boundary_polygon is None: 291 #What if mux file only contains one point 292 msg = 'Files of type sts require boundary polygon' 293 raise msg 294 284 295 # Get first timestep 285 296 try: … … 308 319 y = reshape(y, (len(y),1)) 309 320 vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array 321 322 if boundary_polygon is not None: 323 #removes sts points that do not lie on boundary 324 boundary_polygon=ensure_numeric(boundary_polygon) 325 boundary_polygon[:,0] -= xllcorner 326 boundary_polygon[:,1] -= yllcorner 327 temp=[] 328 boundary_id=[] 329 gauge_id=[] 330 for i in range(len(boundary_polygon)): 331 for j in range(len(x)): 332 if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4): 333 #FIX ME: 334 #currently gauges lat and long is stored as float and 335 #then cast to double. This cuases slight repositioning 336 #of vertex_coordinates. 337 temp.append(boundary_polygon[i]) 338 gauge_id.append(j) 339 boundary_id.append(i) 340 break 341 gauge_neighbour_id=[] 342 for i in range(len(boundary_id)-1): 343 if boundary_id[i]+1==boundary_id[i+1]: 344 gauge_neighbour_id.append(i+1) 345 else: 346 gauge_neighbour_id.append(-1) 347 if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 and boundary_id[0]==0: 348 gauge_neighbour_id.append(0) 349 else: 350 gauge_neighbour_id.append(-1) 351 gauge_neighbour_id=ensure_numeric(gauge_neighbour_id) 352 if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1: 353 msg='incorrect number of segments' 354 raise msg 355 vertex_coordinates=ensure_numeric(temp) 356 if len(vertex_coordinates)==0: 357 msg = 'None of the sts gauges fall on the boundary' 358 raise msg 359 else: 360 gauge_neighbour_id=None 310 361 311 362 if interpolation_points is not None: … … 342 393 for i, name in enumerate(quantity_names): 343 394 quantities[name] = fid.variables[name][:] 395 if boundary_polygon is not None: 396 #removes sts points that do not lie on boundary 397 quantities[name] = take(quantities[name],gauge_id,1) 344 398 fid.close() 345 346 399 347 400 from anuga.fit_interpolate.interpolate import Interpolation_function … … 362 415 interpolation_points, 363 416 time_thinning=time_thinning, 364 verbose=verbose ), starttime417 verbose=verbose,gauge_neighbour_id=gauge_neighbour_id), starttime 365 418 366 419 # NOTE (Ole): Caching Interpolation function is too slow as -
anuga_core/source/anuga/fit_interpolate/interpolate.py
r5361 r5372 93 93 f, 94 94 vertex_coordinates, 95 gauge_neighbour_id, 95 96 point_coordinates=None, 96 start_blocking_len=500000,97 97 verbose=False): 98 98 … … 112 112 assert f.shape[0]>0,msg 113 113 114 #print point_on_line_py(point_coordinates[3],[vertex_coordinates[0],vertex_coordinates[1]])115 #print vertex_coordinates[0],vertex_coordinates[1]116 117 #print point_coordinates118 #print vertex_coordinates119 114 n=f.shape[0] 120 115 if n==1: 121 116 z=f*z 122 117 elif n>1: 123 index=0#index of segment on which last point was found124 k=0#number of points on line125 118 for i in range(len(point_coordinates)): 126 119 found = False 127 #find the segment the cetnroid lies on 128 #Start with the segment the last point was found on 129 #For n points there will be n-1 segments 130 for j in range(n-1): 131 #print 'searcing segment', index+j 132 #reached last segment look at first segment 133 if (index+j)==n-2: 134 index=0 135 if point_on_line_py(point_coordinates[i],[vertex_coordinates[(j)+index],vertex_coordinates[(j+1)+index]]): 136 found=True 137 x0=vertex_coordinates[j][0];y0=vertex_coordinates[j][1] 138 x1=vertex_coordinates[j+1][0];y1=vertex_coordinates[j+1][1] 139 x2=point_coordinates[i][0];y2=point_coordinates[i][1] 140 141 segment_len=sqrt((x1-x0)**2+(y1-y0)**2) 142 dist=sqrt((x2-x0)**2+(y2-y0)**2) 143 z[i]=(f[j+1]-f[j])/segment_len*dist+f[j] 144 #print 'element found on segment',j+index 145 index=j 146 break 147 148 120 for j in range(n): 121 if gauge_neighbour_id[j]>=0: 122 if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]): 123 found=True 124 x0=vertex_coordinates[j][0] 125 y0=vertex_coordinates[j][1] 126 x1=vertex_coordinates[gauge_neighbour_id[j]][0] 127 y1=vertex_coordinates[gauge_neighbour_id[j]][1] 128 x2=point_coordinates[i][0] 129 y2=point_coordinates[i][1] 130 131 segment_len=sqrt((x1-x0)**2+(y1-y0)**2) 132 dist=sqrt((x2-x0)**2+(y2-y0)**2) 133 z[i]=(f[gauge_neighbour_id[j]]-f[j])/segment_len*dist+f[j] 134 #print 'element found on segment' 135 break 136 149 137 if not found: 150 138 z[i]=0.0 151 139 #print 'point not on urs boundary' 152 else:153 k+=1154 #print 'number of midpoints on urs boundary',k155 #print z156 140 return z 157 141 … … 550 534 interpolation_points=None, 551 535 time_thinning=1, 552 verbose=False): 536 verbose=False, 537 gauge_neighbour_id=None): 553 538 """Initialise object and build spatial interpolation if required 554 539 … … 749 734 verbose=False) # Don't clutter 750 735 elif triangles is None and vertex_coordinates is not None: 751 result=interpol.interpolate_polyline(Q,vertex_coordinates, point_coordinates=self.interpolation_points)736 result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points) 752 737 753 738 #assert len(result), len(interpolation_points) -
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5369 r5372 6441 6441 domain_fbound = pmesh_to_domain_instance(meshname, Domain) 6442 6442 domain_fbound.set_quantity('stage', tide) 6443 Bf = File_boundary(sts_file+'.sts', domain_fbound )6443 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon) 6444 6444 Br = Reflective_boundary(domain_fbound) 6445 6445 Bd=Dirichlet_boundary([2.0,220,-220]) … … 6466 6466 i+=1 6467 6467 6468 assert temp_fbound-temp_drchlt<epsilon 6468 assert allclose(temp_fbound,temp_drchlt) 6469 os.remove(sts_file+'.sts') 6470 os.remove(meshname) 6471 6472 def test_file_boundary_sts2(self): 6473 # mux2 file has points not included in boundary 6474 # mux2 gauges are not stored with the same order as they are 6475 # found in bounding_polygon 6476 from anuga.shallow_water import Domain 6477 from anuga.shallow_water import Reflective_boundary 6478 from anuga.shallow_water import Dirichlet_boundary 6479 from anuga.shallow_water import File_boundary 6480 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6481 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 6482 bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]] 6483 tide=0. 6484 time_step_count = 5 6485 time_step = 2 6486 lat_long_points=bounding_polygon[0:2] 6487 lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1]) 6488 lat_long_points.insert(0,[6.0,97.01]) 6489 n=len(lat_long_points) 6490 first_tstep=ones(n,Int) 6491 last_tstep=(time_step_count)*ones(n,Int) 6492 gauge_depth=20*ones(n,Float) 6493 ha=2*ones((n,time_step_count),Float) 6494 ua=10*ones((n,time_step_count),Float) 6495 va=-10*ones((n,time_step_count),Float) 6496 base_name, files = self.write_mux2(lat_long_points, 6497 time_step_count, time_step, 6498 first_tstep, last_tstep, 6499 depth=gauge_depth, 6500 ha=ha, 6501 ua=ua, 6502 va=va) 6503 6504 sts_file=base_name 6505 urs2sts(base_name,sts_file,mean_stage=tide,verbose=False) 6506 self.delete_mux(files) 6507 6508 #print 'start create mesh from regions' 6509 for i in range(len(bounding_polygon)): 6510 zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1]) 6511 extent_res=1000000 6512 meshname = 'urs_test_mesh' + '.tsh' 6513 interior_regions=None 6514 boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]} 6515 #have to change boundary tags from last example because now bounding 6516 #polygon starts in different place. 6517 create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags, 6518 maximum_triangle_area=extent_res,filename=meshname, 6519 interior_regions=interior_regions,verbose=False) 6520 6521 domain_fbound = pmesh_to_domain_instance(meshname, Domain) 6522 domain_fbound.set_quantity('stage', tide) 6523 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon) 6524 Br = Reflective_boundary(domain_fbound) 6525 Bd=Dirichlet_boundary([2.0,220,-220]) 6526 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 6527 finaltime=time_step*(time_step_count-1) 6528 yieldstep=time_step 6529 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 6530 i=0 6531 for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 6532 skip_initial_step = False): 6533 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 6534 i+=1 6535 6536 domain_drchlt = pmesh_to_domain_instance(meshname, Domain) 6537 domain_drchlt.set_quantity('stage', tide) 6538 Br = Reflective_boundary(domain_drchlt) 6539 Bd=Dirichlet_boundary([2.0,220,-220]) 6540 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 6541 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 6542 i=0 6543 for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 6544 skip_initial_step = False): 6545 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 6546 i+=1 6547 6548 assert allclose(temp_fbound,temp_drchlt) 6469 6549 os.remove(sts_file+'.sts') 6470 6550 os.remove(meshname) -
anuga_core/source/anuga/utilities/polygon.py
r5321 r5372 16 16 from anuga.utilities.numerical_tools import ensure_numeric 17 17 from anuga.geospatial_data.geospatial_data import ensure_absolute 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 18 60 19 61 def point_on_line(point, line):
Note: See TracChangeset
for help on using the changeset viewer.