source: anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py @ 5959

Last change on this file since 5959 was 4783, checked in by duncan, 17 years ago

modified to aid in debugging

File size: 4.7 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
14
15
16
17#--------------------------------------------------------------------------
18# Create the triangular mesh based on overall clipping polygon with a
19# tagged
20# boundary and interior regions defined in project.py along with
21# resolutions (maximal area of per triangle) for each polygon
22#--------------------------------------------------------------------------
23
24
25base_resolution = 1 # Use this to coarsen or refine entire mesh.
26
27# Basic geometry and bounding polygon
28xleft   = 0
29xright  = 5.448
30ybottom = 0
31ytop    = 3.402
32
33point_sw = [xleft, ybottom]
34point_se = [xright, ybottom]
35point_nw = [xleft, ytop]   
36point_ne = [xright, ytop]
37
38bounding_polygon = [point_se,
39                    point_ne,
40                    point_nw,
41                    point_sw]
42
43
44# Localised refined area for gulleys
45xl = 4.8
46xr = 5.3
47yb = 1.6
48yt = 2.3
49p0 = [xl, yb]
50p1 = [xr, yb]
51p2 = [xr, yt]
52p3 = [xl, yt]
53
54gulleys = [p0, p1, p2, p3]
55
56
57
58# Island area and drawdown region
59island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
60island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
61island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
62island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
63island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
64island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
65island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
66island_7 = [xl-.01, yt]
67
68island = [island_0, island_1, island_2,
69          island_3, island_4, island_5,
70          island_6, island_7]
71
72
73
74# Region spanning half right hand side of domain just inside boundary
75rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
76rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
77rhs_se = [xright-0.1, ybottom+0.2]
78rhs_ne = [xright-0.1, ytop-0.2]       
79
80rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
81
82
83# Interior regions and creation of mesh
84interior_regions = [[rhs_region, 0.0005],
85                    [island, 0.0003*base_resolution],
86                    [gulleys, 0.00003*base_resolution]]   
87
88
89
90
91def prepare_timeboundary(filename):
92    """Convert benchmark 2 time series to NetCDF tms file.
93    This is a 'throw-away' code taylor made for files like
94    'Benchmark_2_input.txt' from the LWRU2 benchmark
95    """
96
97    from Scientific.IO.NetCDF import NetCDFFile
98    from Numeric import array
99
100
101    print 'Creating', filename
102
103    # Read the ascii (.txt) version of this file
104    fid = open(filename[:-4] + '.txt')
105
106    # Skip first line
107    line = fid.readline()
108
109    # Read remaining lines
110    lines = fid.readlines()
111    fid.close()
112
113
114    N = len(lines)
115    T = zeros(N, Float)  #Time
116    Q = zeros(N, Float)  #Values
117
118    for i, line in enumerate(lines):
119        fields = line.split()
120
121        T[i] = float(fields[0])
122        Q[i] = float(fields[1])
123
124
125    # Create tms NetCDF file
126
127    fid = NetCDFFile(filename, 'w')
128    fid.institution = 'Geoscience Australia'
129    fid.description = 'Input wave for Benchmark 2'
130    fid.starttime = 0.0
131    fid.createDimension('number_of_timesteps', len(T))
132    fid.createVariable('time', Float, ('number_of_timesteps',))
133    fid.variables['time'][:] = T
134
135    fid.createVariable('stage', Float, ('number_of_timesteps',))
136    fid.variables['stage'][:] = Q[:]
137
138    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
139    fid.variables['xmomentum'][:] = 0.0
140
141    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
142    fid.variables['ymomentum'][:] = 0.0
143
144    fid.close()
145
146
147def create_mesh(elevation_in_mesh=False):
148    # Prepare time boundary
149    prepare_timeboundary(project.boundary_filename)
150
151
152
153    meshname = project.mesh_filename + '.msh'
154    m = create_mesh_from_regions(bounding_polygon,
155                                 boundary_tags={'wall': [0, 1, 3],
156                                                'wave': [2]},     
157                                 maximum_triangle_area=0.1*base_resolution,
158                                 interior_regions=interior_regions,
159                                 filename=project.mesh_filename,
160                                 use_cache=False,
161                                 verbose=True)
162   
163    if elevation_in_mesh is True:
164        fit_to_mesh_file(project.mesh_filename, project.bathymetry_filename,
165                         project.mesh_filename)
166
167
168#-------------------------------------------------------------
169if __name__ == "__main__":
170    create_mesh()
171
172
Note: See TracBrowser for help on using the repository browser.