source: anuga_validation/okushiri_2005/create_okushiri.py @ 5848

Last change on this file since 5848 was 3915, checked in by ole, 17 years ago

Okushiri demo

File size: 6.5 KB
Line 
1"""Create mesh and time boundary for the Okushiri island validation
2"""
3
4
5from Numeric import array, zeros, Float, allclose
6
7from anuga.pmesh.mesh import *
8from anuga.pmesh.mesh_interface import create_mesh_from_regions
9from anuga.coordinate_transforms.geo_reference import Geo_reference
10from anuga.geospatial_data import Geospatial_data
11
12import project
13
14def prepare_bathymetry(filename):
15    """Convert benchmark 2 bathymetry to NetCDF pts file.
16    This is a 'throw-away' code taylor made for files like
17    'Benchmark_2_bathymetry.txt' from the LWRU2 benchmark
18    """
19
20    print 'Creating', filename
21   
22    # Read the ascii (.txt) version of this file,
23    # make it comma separated and invert the bathymetry
24    # (Below mean sea level should be negative)
25    infile = open(filename[:-4] + '.txt')
26
27    points = []
28    attribute = []
29    for line in infile.readlines()[1:]: #Skip first line (the header)
30        fields = line.strip().split()
31
32        x = float(fields[0])
33        y = float(fields[1])
34        z = -float(fields[2]) # Bathymetry is inverted in original file
35       
36        points.append([x,y])
37        attribute.append(z)
38    infile.close()
39
40    # Convert to geospatial data and store as NetCDF
41    G = Geospatial_data(data_points=points,
42                        attributes=attribute)
43    G.export_points_file(filename)
44   
45
46
47def prepare_timeboundary(filename):
48    """Convert benchmark 2 time series to NetCDF tms file.
49    This is a 'throw-away' code taylor made for files like
50    'Benchmark_2_input.txt' from the LWRU2 benchmark
51    """
52
53    from Scientific.IO.NetCDF import NetCDFFile
54    from Numeric import array
55
56
57    print 'Creating', filename
58
59    # Read the ascii (.txt) version of this file
60    fid = open(filename[:-4] + '.txt')
61
62    # Skip first line
63    line = fid.readline()
64
65    # Read remaining lines
66    lines = fid.readlines()
67    fid.close()
68
69
70    N = len(lines)
71    T = zeros(N, Float)  #Time
72    Q = zeros(N, Float)  #Values
73
74    for i, line in enumerate(lines):
75        fields = line.split()
76
77        T[i] = float(fields[0])
78        Q[i] = float(fields[1])
79
80
81    # Create tms NetCDF file
82
83    fid = NetCDFFile(filename, 'w')
84    fid.institution = 'Geoscience Australia'
85    fid.description = 'Input wave for Benchmark 2'
86    fid.starttime = 0.0
87    fid.createDimension('number_of_timesteps', len(T))
88    fid.createVariable('time', Float, ('number_of_timesteps',))
89    fid.variables['time'][:] = T
90
91    fid.createVariable('stage', Float, ('number_of_timesteps',))
92    fid.variables['stage'][:] = Q[:]
93
94    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
95    fid.variables['xmomentum'][:] = 0.0
96
97    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
98    fid.variables['ymomentum'][:] = 0.0
99
100    fid.close()
101
102
103#-------------------------------------------------------------
104if __name__ == "__main__":
105
106
107    # Prepare bathymetry
108    prepare_bathymetry(project.bathymetry_filename)
109
110    # Prepare time boundary
111    prepare_timeboundary(project.boundary_filename)
112
113
114    #--------------------------------------------------------------------------
115    # Create the triangular mesh based on overall clipping polygon with a
116    # tagged
117    # boundary and interior regions defined in project.py along with
118    # resolutions (maximal area of per triangle) for each polygon
119    #--------------------------------------------------------------------------
120
121
122    base_resolution = 1 # Use this to coarsen or refine entire mesh.
123
124    # Basic geometry and bounding polygon
125    xleft   = 0
126    xright  = 5.448
127    ybottom = 0
128    ytop    = 3.402
129
130    point_sw = [xleft, ybottom]
131    point_se = [xright, ybottom]
132    point_nw = [xleft, ytop]   
133    point_ne = [xright, ytop]
134
135    bounding_polygon = [point_se,
136                        point_ne,
137                        point_nw,
138                        point_sw]
139   
140
141    # Localised refined area for gulleys
142    xl = 4.8
143    xr = 5.3
144    yb = 1.6
145    yt = 2.3
146    p0 = [xl, yb]
147    p1 = [xr, yb]
148    p2 = [xr, yt]
149    p3 = [xl, yt]
150   
151    gulleys = [p0, p1, p2, p3]
152   
153
154    # Island area and drawdown region (original)
155    #island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
156    #island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
157    #island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
158    #island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
159    #island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
160    #island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
161    #island_6 = [xl-.01, yb]  #OK
162    #island_7 = [xl-.01, yt]  #OK     
163
164
165    # Island area and drawdown region
166    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
167    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
168    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
169    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
170    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
171    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
172    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
173    island_7 = [xl-.01, yt]
174 
175    island = [island_0, island_1, island_2,
176              island_3, island_4, island_5,
177              island_6, island_7]
178
179
180    # Region spanning half right hand side of domain just inside boundary (org)
181    #rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
182    #rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
183    #rhs_se = [xright-0.02, ybottom+0.02]
184    #rhs_ne = [xright-0.02, ytop-0.02]
185
186    # Region spanning half right hand side of domain just inside boundary
187    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
188    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
189    rhs_se = [xright-0.1, ybottom+0.2]
190    rhs_ne = [xright-0.1, ytop-0.2]       
191
192    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
193
194   
195    # Interior regions and creation of mesh
196    interior_regions = [[rhs_region, 0.0005],
197                        [island, 0.0002*base_resolution],
198                        [gulleys, 0.00002*base_resolution]]   
199
200    meshname = project.mesh_filename + '.msh'
201    m = create_mesh_from_regions(bounding_polygon,
202                                 boundary_tags={'wall': [0, 1, 3],
203                                                'wave': [2]},     
204                                 maximum_triangle_area=0.1*base_resolution,
205                                 interior_regions=interior_regions,
206                                 filename=project.mesh_filename,
207                                 verbose=True)
208
Note: See TracBrowser for help on using the repository browser.