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