1 | """ This Script shows how a bad mesh can be generated causing |
---|
2 | excessive values for stage and momentum. |
---|
3 | |
---|
4 | """ |
---|
5 | |
---|
6 | #------------------------------------------------------------------------------ |
---|
7 | # Import necessary modules |
---|
8 | #------------------------------------------------------------------------------ |
---|
9 | |
---|
10 | # Standard modules |
---|
11 | import os |
---|
12 | import time |
---|
13 | import sys |
---|
14 | |
---|
15 | # Related major packages |
---|
16 | from anuga.shallow_water import Domain |
---|
17 | from anuga.shallow_water import Dirichlet_boundary |
---|
18 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
19 | from anuga.pmesh.mesh import importMeshFromFile |
---|
20 | from anuga.utilities.polygon import read_polygon, plot_polygons |
---|
21 | |
---|
22 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
23 | |
---|
24 | |
---|
25 | |
---|
26 | #------------------------------------------------------------------------------ |
---|
27 | # Filenames |
---|
28 | #------------------------------------------------------------------------------ |
---|
29 | |
---|
30 | outname ='mesh' |
---|
31 | meshname ='mesh.msh' |
---|
32 | |
---|
33 | # Create a points file |
---|
34 | ptfile = 'points_from_a_file.csv' |
---|
35 | file = open(ptfile,"w") |
---|
36 | file.write("x,y\n\ |
---|
37 | 304200.0, 6185300.0\n\ |
---|
38 | 304250.0, 6185350.0\n\ |
---|
39 | 304300.0, 6185400.0\n") |
---|
40 | file.close() |
---|
41 | |
---|
42 | #------------------------------------------------------------------------------ |
---|
43 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
44 | # boundary and interior regions defined in project.py along with |
---|
45 | # resolutions (maximal area of per triangle) for each polygon |
---|
46 | #------------------------------------------------------------------------------ |
---|
47 | |
---|
48 | W=304180.0 |
---|
49 | S=6185270.0 |
---|
50 | E=307650.0 |
---|
51 | N=6189040.0 |
---|
52 | |
---|
53 | bounding_polygon = [[W, S], [E, S], [E, N], [W, N]] |
---|
54 | |
---|
55 | |
---|
56 | |
---|
57 | create_mesh_from_regions(bounding_polygon, |
---|
58 | boundary_tags={'south': [0], |
---|
59 | 'east': [1], |
---|
60 | 'north': [2], |
---|
61 | 'west': [3]}, |
---|
62 | maximum_triangle_area=10000, #background mesh size |
---|
63 | filename=meshname, |
---|
64 | use_cache=False, |
---|
65 | verbose=True) |
---|
66 | |
---|
67 | m = importMeshFromFile(meshname) |
---|
68 | |
---|
69 | m.add_vertices([(W +500,S +500),(W +1000,S +1000) ]) |
---|
70 | |
---|
71 | # load points from a points file |
---|
72 | m.add_vertices(Geospatial_data(ptfile)) |
---|
73 | |
---|
74 | p1 = (W +500,N - 500) |
---|
75 | p2 = (W +1000,N - 1000) |
---|
76 | p3 = (W +1000,N - 1500) |
---|
77 | m.add_points_and_segments([p1,p2,p3], [[0,1],[1,2]]) |
---|
78 | |
---|
79 | # add points and segments, where the segments are automatically added. |
---|
80 | # This will work on the next release of ANUGA. |
---|
81 | # It doesn't work on the current sourceforge release. |
---|
82 | p1 = (E -1000,N - 1000) |
---|
83 | p2 = (E -1000,N - 1500) |
---|
84 | #m.add_points_and_segments([p1,p2]) |
---|
85 | |
---|
86 | m.generate_mesh() |
---|
87 | m.export_mesh_file('end.msh') |
---|
88 | |
---|
89 | import sys; sys.exit() |
---|
90 | #------------------------------------------------------------------------------ |
---|
91 | # Setup computational domain |
---|
92 | #------------------------------------------------------------------------------ |
---|
93 | |
---|
94 | domain = Domain(meshname, use_cache=False, verbose=True) |
---|
95 | domain.set_name(outname) |
---|
96 | print domain.statistics() |
---|
97 | |
---|
98 | #------------------------------------------------------------------------------ |
---|
99 | # Setup initial conditions |
---|
100 | #------------------------------------------------------------------------------ |
---|
101 | # Setting up initial conditions for a Flood model includes: |
---|
102 | # Setting A Base topography over the model domain |
---|
103 | # Setting the variation of Roughness over the terrain, which may be |
---|
104 | # varied dependent on depth say |
---|
105 | # Setting initial water levels at BOundaries and in Lakes and Water Bodies. |
---|
106 | |
---|
107 | |
---|
108 | tide = 1.5 # 100yr tidal level as per WCC drainage design code |
---|
109 | base_friction = 0.015 |
---|
110 | |
---|
111 | domain.set_quantity('friction', base_friction) |
---|
112 | domain.set_quantity('stage', tide) |
---|
113 | domain.set_quantity('elevation', 0) |
---|
114 | |
---|
115 | |
---|
116 | #------------------------------------------------------------------------------ |
---|
117 | # Setup boundary conditions |
---|
118 | #------------------------------------------------------------------------------ |
---|
119 | |
---|
120 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
121 | |
---|
122 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
123 | |
---|
124 | # Boundary conditions |
---|
125 | domain.set_boundary({'west': Bd, |
---|
126 | 'south': Bd, |
---|
127 | 'north': Bd, |
---|
128 | 'east': Bd}) |
---|
129 | |
---|
130 | |
---|
131 | #------------------------------------------------------------------------------ |
---|
132 | # Evolve system through time |
---|
133 | #------------------------------------------------------------------------------ |
---|
134 | |
---|
135 | import time |
---|
136 | t0 = time.time() |
---|
137 | |
---|
138 | for t in domain.evolve(yieldstep = 1, finaltime = 2): |
---|
139 | print domain.timestepping_statistics() |
---|
140 | |
---|
141 | print 'Finished' |
---|
142 | |
---|