1 | """Create mesh and time boundary for the Okushiri island validation |
---|
2 | """ |
---|
3 | |
---|
4 | |
---|
5 | from Numeric import array, zeros, Float, allclose |
---|
6 | |
---|
7 | from anuga.pmesh.mesh import * |
---|
8 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
9 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
10 | from anuga.geospatial_data import Geospatial_data |
---|
11 | |
---|
12 | import project |
---|
13 | |
---|
14 | |
---|
15 | |
---|
16 | |
---|
17 | #------------------------------------------------------------- |
---|
18 | if __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 | |
---|