Changeset 1832
- Timestamp:
- Sep 14, 2005, 6:22:00 PM (19 years ago)
- Location:
- production/karratha_2005
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
production/karratha_2005/create_mesh.py
r1830 r1832 6 6 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 7 7 8 from pmesh.mesh import * 8 9 9 from pyvolution.coordinate_transforms.geo_reference import Geo_reference 10 10 from pyvolution.coordinate_transforms.redfearn import redfearn 11 11 12 12 13 import project14 13 15 #------------------------------------------------------------- 16 if __name__ == "__main__": 14 def create_mesh_boundary(polygon, tags, resolution, filename = None): 15 """Create mesh from bounding polygon, tags for all segments and resolution. 16 17 Polygon is a list of points in latitude and longitudes - decimal degrees 17 18 18 19 20 Tags is a list of symbolic tags 21 22 Resolution is the maximal area 19 23 20 #Basic geometry21 24 22 south = project.south 23 north = project.north 24 west = project.west 25 east = project.east 25 FIXME - FUTURE: 26 Use class Point 27 Take multiple polygons 26 28 27 refzone = project.refzone 29 Accept deg, min, sec, e.g. 30 [ [(-20,30,55), (116, 20, 00)], ...] 31 32 33 """ 34 35 from pmesh.mesh import * 36 from pyvolution.coordinate_transforms.redfearn import redfearn 37 38 39 #Convert to UTM 40 import project 41 refzone = project.refzone #FIXME 28 42 29 43 30 #SW 31 zone, easting, northing = redfearn(south, west) 32 assert zone == refzone 33 point_sw = [easting, northing] 44 points = [] 45 for point in polygon: 46 zone, easting, northing = redfearn(point[0], point[1]) #FIXME: Use point.latitude etc 47 assert zone == refzone 48 points.append([easting, northing]) 34 49 35 #SE 36 zone, easting, northing = redfearn(south, east) 37 assert zone == refzone 38 point_se = [easting, northing] 50 polygon = points 39 51 40 #NW 41 zone, easting, northing = redfearn(north, west) 42 assert zone == refzone 43 point_nw = [easting, northing] 52 53 mesh_origin = project.mesh_origin 44 54 45 #NE 46 zone, easting, northing = redfearn(north, east) 47 assert zone == refzone 48 point_ne = [easting, northing] 55 geo = Geo_reference(xllcorner = mesh_origin[1], #From dem 56 yllcorner = mesh_origin[2], 57 zone = refzone) 58 59 60 # 61 dict = {} 62 dict['points'] = polygon 63 64 #Create segments 65 #E.g. [[0,1], [1,2], [2,3], [3,0]] 66 segments = [] 67 N = len(polygon) 68 for i in range(N): 69 lo = i 70 hi = (lo + 1) % N 71 segments.append( [lo, hi] ) 72 73 74 dict['segments'] = segments 75 76 77 #Create tags 78 #E.g. ['wall', 'wall', 'ocean', 'wall'] 79 segment_tags = [0]*N 80 for key in tags: 81 indices = tags[key] 82 for i in indices: 83 segment_tags[i] = key 49 84 85 dict['segment_tags'] = segment_tags 50 86 51 #geo = Geo_reference(xllcorner = 421544.35127423, #from dem52 # yllcorner = 7677669.5257159,53 # zone = refzone)54 87 55 geo = Geo_reference(xllcorner = 420468.31429902, #From dem56 yllcorner = 7677669.5257159,57 zone = refzone)58 88 59 #geo = Geo_reference(xllcorner = point_sw[0]+1000, #FIXME: BAD WHY60 # yllcorner = point_sw[1]+1000,61 # zone = refzone)62 63 64 89 print "***********************" 65 90 print "geo ref", geo … … 67 92 68 93 m = Mesh(geo_reference=geo) 69 70 #Boundary71 dict = {}72 dict['points'] = [point_sw,73 point_se,74 point_ne,75 point_nw]76 77 78 dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]79 dict['segment_tags'] = ['wall', 'wall', 'ocean', 'wall']80 94 81 95 m.addVertsSegs(dict) 82 96 83 97 84 base_resolution = 50000 98 #FIXME: Find centroid automatically 99 inner = m.addRegionEN(points[0][0]+1, points[0][1]+1) 100 inner.setMaxArea(resolution) 101 102 m.generateMesh('pzq28.0za1000000a') 103 104 m.export_mesh_file(filename) 105 85 106 86 inner = m.addRegionEN(point_sw[0]+1, point_sw[1]+1)87 inner.setMaxArea(base_resolution)88 107 89 m.generateMesh('pzq28.0za1000000a')90 108 91 m.export_mesh_file(project.meshname + '.msh') 109 92 110 -
production/karratha_2005/project.py
r1830 r1832 32 32 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees 33 33 34 #Origin of existing dem (FIXME: Temporary measure) 35 #mesh_origin = (50, 421544.35127423, 7677669.5257159) #250m 36 mesh_origin = (50, 420468.31429902, 7677669.5257159) #100m 37 38 39 #south = degminsec2decimal_degrees(-20,45,0) 40 #north = degminsec2decimal_degrees(-20,15,0) 41 #west = degminsec2decimal_degrees(116,30,0) 42 #east = degminsec2decimal_degrees(117,0,0) 43 34 44 south = degminsec2decimal_degrees(-20,45,0) 35 45 north = degminsec2decimal_degrees(-20,15,0) 36 west = degminsec2decimal_degrees(116,30,0) 37 east = degminsec2decimal_degrees(117,0,0) 46 west = degminsec2decimal_degrees(116,20,0) 47 east = degminsec2decimal_degrees(117,10,0) 48 49 p0 = [south, degminsec2decimal_degrees(116,32,0)] 50 p1 = [south, west] 51 p2 = [degminsec2decimal_degrees(-20,23,0), west] 52 p3 = [north, degminsec2decimal_degrees(116,45,0)] 53 p4 = [north, degminsec2decimal_degrees(117,0,0)] 54 p5 = [p2[0], degminsec2decimal_degrees(117,8,0)] 55 p6 = [degminsec2decimal_degrees(-20,30,0), east] 56 p7 = [degminsec2decimal_degrees(-20,38,0), east] 57 p8 = [south, east] 58 59 polygon = [p0, p1, p2, p3, p4, p5, p6, p7, p8] 38 60 39 61 refzone = 50 -
production/karratha_2005/run_karratha.py
r1830 r1832 7 7 8 8 9 from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary, Dirichlet_boundary, Time_boundary 9 from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary 10 10 from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\ 11 11 dem2pts, ferret2sww … … 41 41 east = project.east 42 42 43 #Origin of existing dem (FIXME: Temporary measure) 44 #mesh_origin = (50, 421544.35127423, 7677669.5257159) #250m 45 mesh_origin = (50, 420468.31429902, 7677669.5257159) #100m 43 44 #Mesh 45 from create_mesh import create_mesh_boundary 46 47 m = cache(create_mesh_boundary, 48 project.polygon, 49 {'tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]}, 50 'resolution': 52000, 51 'filename': project.meshname + '.msh'}, 52 verbose = True) 53 54 55 56 57 58 #tags = {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]}, 59 #resolution = 200000, 60 #filename = project.meshname + '.msh') 61 62 63 46 64 47 65 ferret2sww(source_dir+project.boundary_basename, … … 50 68 minlat=south-1, maxlat=north+1, 51 69 minlon=west-1, maxlon=east+1, 52 origin = mesh_origin,70 origin = project.mesh_origin, 53 71 zscale = 1, 54 72 fail_on_NaN = False, … … 68 86 domain.store = True 69 87 70 print "Number of triangles = ", len(domain)88 print 'Number of triangles = ', len(domain) 71 89 print 'The extent is ', domain.get_extent() 72 90 73 91 74 92 #Setup IC 75 domain.set_quantity('stage', 0) 93 tide = 0 94 domain.set_quantity('stage', tide) 76 95 domain.set_quantity('elevation', 77 96 filename = demname + '.pts', … … 83 102 print domain.get_boundary_tags() 84 103 85 #Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True) 104 Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True) 105 domain.starttime = 3000 #Obtained from MOST 106 86 107 87 108 Br = Reflective_boundary(domain) 88 Bd = Dirichlet_boundary([3,0,0]) 109 Bt = Transmissive_boundary(domain) 110 Bd = Dirichlet_boundary([tide,0,0]) 89 111 Bw = Time_boundary(domain=domain, 90 f=lambda t: [(60<t<660)*12, 0, 0]) 91 domain.set_boundary( {'wall': Br, 'ocean': Bw} ) 112 f=lambda t: [(60<t<660)*4, 0, 0]) 113 114 #domain.set_boundary( {'wall': Bd, 'ocean': Bf} ) 115 domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bf} ) 92 116 93 117 94 118 #Run 95 for t in domain.evolve(yieldstep = 15, finaltime = 4800): 119 120 for t in domain.evolve(yieldstep = 60, finaltime = 28610): 96 121 domain.write_time() 97 domain.write_boundary_statistics( quantities = 'stage')122 domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
Note: See TracChangeset
for help on using the changeset viewer.