Changeset 4797
- Timestamp:
- Nov 7, 2007, 2:19:42 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_validation/conical_island/run_circular.py
r4777 r4797 10 10 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary 11 11 from anuga.shallow_water import Dirichlet_boundary, Time_boundary 12 from anuga.shallow_water.data_manager import get_maximum_inundation_data 12 13 from anuga.abstract_2d_finite_volumes.util import file_function 13 14 from anuga.pmesh.mesh_interface import create_mesh_from_regions 15 from anuga.utilities.polygon import read_polygon, plot_polygons 14 16 #from cmath import cos,sin,cosh,pi 15 17 from math import cos,pi,sin,tan,sqrt … … 22 24 #------------------------- 23 25 26 def get_xy(x=0,y=0,r=1,angle=1): 27 x1 = x + cos(angle)*r 28 y1 = y + sin(-angle)*r 29 # print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle) 30 return x1,y1, 24 31 25 32 def prepare_timeboundary(filename): … … 50 57 for i, line in enumerate(lines): 51 58 fields = line.split() 52 print 'fields',i,line59 # print 'fields',i,line 53 60 54 61 T[i] = float(fields[0]) … … 82 89 83 90 84 85 poly_all= [[0,0],[0,25],[30,25],[30,0]] 86 poly=[[9,10],[9,17.6],[16.6,17.6],[16.6,10]] 87 internal_poly=[[poly,1]] 91 angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5, 92 95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 93 247.5, 270.0, 292.5, 315.0, 337.5) 94 angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 95 247.5, 270.0, 292.5, 315.0, 337.5) 96 97 x = 12.96 98 y = 13.80 99 r1 = 1.5 100 r2 = 2.5 101 d=0.75 102 res=0.0005 103 104 poly = [] 105 internal_poly=[] 106 for i,angle in enumerate(angles): 107 angle = ((angle-180)/180.0)*pi 108 # p1 = get_xy(x,y,r1,angle-0.01745*d) 109 # p2 = get_xy(x,y,r2,angle-0.01745*d) 110 # p3 = get_xy(x,y,r2,angle+0.01745*d) 111 # p4 = get_xy(x,y,r1,angle+0.01745*d) 112 # poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]] 113 114 p = get_xy(x,y,4,angle) 115 poly.insert(0,[p[0],p[1]]) 116 internal_poly.insert(0,(poly,1)) 117 118 poly1 = [] 119 for i,angle in enumerate(angles): 120 angle = ((angle-180)/180.0)*pi 121 p = get_xy(x,y,2.,angle) 122 poly1.insert(0,[p[0],p[1]]) 123 124 internal_poly.insert(1,(poly1,.1)) 125 126 poly2 = [] 127 for i,angle in enumerate(angles): 128 angle = ((angle-180)/180.0)*pi 129 p = get_xy(x,y,1.,angle) 130 poly2.insert(0,[p[0],p[1]]) 131 132 internal_poly.insert(2,(poly2,1)) 133 134 print internal_poly 135 136 plot_polygons(internal_poly,'poly.png',) 137 138 poly_all= [[0,7],[0,25],[30,25],[30,7]] 139 #poly1=[[9,10],[9,17.6],[16.6,17.6],[16.6,10]] 140 #internal_poly.insert(0,(poly1,0.025)) 88 141 create_mesh_from_regions(poly_all, 89 boundary_tags={'wall': [0, 1,2],'wave': [3]},90 maximum_triangle_area=1 0,142 boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]}, 143 maximum_triangle_area=1, 91 144 interior_regions=internal_poly, 92 145 filename='temp.msh', … … 101 154 #------------------------- 102 155 domain.set_quantity('friction', 0.0) 103 water_height = 0.32 156 water_height = 0.322 104 157 domain.set_quantity('stage', water_height) 105 158 … … 121 174 z=0.625 122 175 123 print 'x,y,f',i,j,a,b,z176 # print 'x,y,f',i,j,a,b,z 124 177 file.write("%s, %s, %s\n" %(i, j, z)) 125 178 file.close() … … 133 186 # Set simulation parameters 134 187 #------------------------- 135 domain.set_name('test ') # Name of output sww file188 domain.set_name('test1') # Name of output sww file 136 189 domain.set_default_order(1) # Apply second order scheme 137 190 domain.set_all_limiters(0.9) # Max second order scheme (old lim) … … 160 213 161 214 # Create and assign boundary objects 162 boundary_filename="ts2 a_edit_input_32"215 boundary_filename="ts2cnew1_input" 163 216 prepare_timeboundary(boundary_filename+'.txt') 164 217 … … 174 227 Br = Reflective_boundary(domain) 175 228 Bd = Dirichlet_boundary([water_height,0,0]) 176 domain.set_boundary({'wave': Bts, 'wall': B d})229 domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd}) 177 230 178 231 #------------------------- … … 182 235 t0 = time.time() 183 236 184 for t in domain.evolve(yieldstep = 5, finaltime = 75):237 for t in domain.evolve(yieldstep = 1, finaltime = 32): 185 238 domain.write_time() 186 239 # domain.write_time(track_speeds=False) 187 240 241 for t in domain.evolve(yieldstep = 0.1, finaltime = 36 242 ,skip_initial_step = True): 243 domain.write_time() 244 245 for t in domain.evolve(yieldstep = 1, finaltime = 75 246 ,skip_initial_step = True): 247 domain.write_time() 248 188 249 print 'That took %.2f seconds' %(time.time()-t0) 250 251 252 253 #angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5, 254 # 95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 255 # 247.5, 270.0, 292.5, 315.0, 337.5) 256 #x = 12.96 257 #y = 13.80 258 #r1 = 1.1 259 #r2 = 3.6 260 #d=2 261 #poly = {} 262 for i,angle in enumerate(angles): 263 # angle = ((angle-180)/180.0)*pi 264 ## angle = angle1*3.14157 265 # p1 = get_xy(x,y,r1,angle-0.01745*d) 266 # p2 = get_xy(x,y,r2,angle-0.01745*d) 267 # p3 = get_xy(x,y,r2,angle+0.01745*d) 268 # p4 = get_xy(x,y,r1,angle+0.01745*d) 269 # poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]] 270 ## print i,poly[i] 271 # 272 273 run_up, x_y = get_maximum_inundation_data(filename='test.sww',polygon=poly[i], verbose=False) 274 275 print 'maximum_inundation_data',angle, run_up, x_y 276 277 278 279 280 281
Note: See TracChangeset
for help on using the changeset viewer.