source: anuga_work/development/profile_okushiri/create_okushiri.py @ 7860

Last change on this file since 7860 was 4657, checked in by duncan, 17 years ago

Used to benchmark fitting, Okushiri example.

File size: 4.0 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#-------------------------------------------------------------
18if __name__ == "__main__":
19
20
21    #--------------------------------------------------------------------------
22    # Create the triangular mesh based on overall clipping polygon with a
23    # tagged
24    # boundary and interior regions defined in project.py along with
25    # resolutions (maximal area of per triangle) for each polygon
26    #--------------------------------------------------------------------------
27
28
29    base_resolution = 100 # Use this to coarsen or refine entire mesh.
30
31    # Basic geometry and bounding polygon
32    xleft   = 0
33    xright  = 5.448
34    ybottom = 0
35    ytop    = 3.402
36
37    point_sw = [xleft, ybottom]
38    point_se = [xright, ybottom]
39    point_nw = [xleft, ytop]   
40    point_ne = [xright, ytop]
41
42    bounding_polygon = [point_se,
43                        point_ne,
44                        point_nw,
45                        point_sw]
46   
47
48    # Localised refined area for gulleys
49    xl = 4.8
50    xr = 5.3
51    yb = 1.6
52    yt = 2.3
53    p0 = [xl, yb]
54    p1 = [xr, yb]
55    p2 = [xr, yt]
56    p3 = [xl, yt]
57   
58    gulleys = [p0, p1, p2, p3]
59   
60
61    # Island area and drawdown region (original)
62    #island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
63    #island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
64    #island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
65    #island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
66    #island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
67    #island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
68    #island_6 = [xl-.01, yb]  #OK
69    #island_7 = [xl-.01, yt]  #OK     
70
71
72    # Island area and drawdown region
73    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
74    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
75    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
76    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
77    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
78    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
79    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
80    island_7 = [xl-.01, yt]
81 
82    island = [island_0, island_1, island_2,
83              island_3, island_4, island_5,
84              island_6, island_7]
85
86
87    # Region spanning half right hand side of domain just inside boundary (org)
88    #rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
89    #rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
90    #rhs_se = [xright-0.02, ybottom+0.02]
91    #rhs_ne = [xright-0.02, ytop-0.02]
92
93    # Region spanning half right hand side of domain just inside boundary
94    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
95    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
96    rhs_se = [xright-0.1, ybottom+0.2]
97    rhs_ne = [xright-0.1, ytop-0.2]       
98
99    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
100
101   
102    # Interior regions and creation of mesh
103    interior_regions = [[rhs_region, 0.0005],
104                        [island, 0.0002*base_resolution],
105                        [gulleys, 0.00002*base_resolution]]   
106
107    meshname = project.mesh_filename + '.msh'
108    m = create_mesh_from_regions(bounding_polygon,
109                                 boundary_tags={'wall': [0, 1, 3],
110                                                'wave': [2]},     
111                                 maximum_triangle_area=0.1*base_resolution,
112                                 interior_regions=interior_regions,
113                                 filename=project.mesh_filename,
114                                 verbose=True)
115
Note: See TracBrowser for help on using the repository browser.