1 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect |
---|
2 | from anuga.config import velocity_protection, g |
---|
3 | import math |
---|
4 | |
---|
5 | import numpy as num |
---|
6 | |
---|
7 | class Inlet: |
---|
8 | """Contains information associated with each inlet |
---|
9 | """ |
---|
10 | |
---|
11 | def __init__(self, domain, line, verbose=False): |
---|
12 | |
---|
13 | self.domain = domain |
---|
14 | self.domain_bounding_polygon = self.domain.get_boundary_polygon() |
---|
15 | self.line = line |
---|
16 | self.verbose = verbose |
---|
17 | |
---|
18 | self.compute_triangle_indices() |
---|
19 | self.compute_area() |
---|
20 | |
---|
21 | |
---|
22 | def compute_triangle_indices(self): |
---|
23 | |
---|
24 | # Get boundary (in absolute coordinates) |
---|
25 | bounding_polygon = self.domain_bounding_polygon |
---|
26 | domain_centroids = self.domain.get_centroid_coordinates(absolute=True) |
---|
27 | vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True) |
---|
28 | |
---|
29 | # Check that line lies within the mesh. |
---|
30 | for point in self.line: |
---|
31 | msg = 'Point %s ' % str(point) |
---|
32 | msg += ' did not fall within the domain boundary.' |
---|
33 | assert is_inside_polygon(point, bounding_polygon), msg |
---|
34 | |
---|
35 | |
---|
36 | |
---|
37 | self.triangle_indices = line_intersect(vertex_coordinates, self.line) |
---|
38 | |
---|
39 | if len(self.triangle_indices) == 0: |
---|
40 | msg = 'Inlet line=%s ' % (self.line) |
---|
41 | msg += 'No triangles intersecting line ' |
---|
42 | raise Exception, msg |
---|
43 | |
---|
44 | |
---|
45 | |
---|
46 | def compute_area(self): |
---|
47 | |
---|
48 | # Compute inlet area as the sum of areas of triangles identified |
---|
49 | # by line. Must be called after compute_inlet_triangle_indices(). |
---|
50 | if len(self.triangle_indices) == 0: |
---|
51 | region = 'Inlet line=%s' % (self.inlet_line) |
---|
52 | msg = 'No triangles have been identified in region ' |
---|
53 | raise Exception, msg |
---|
54 | |
---|
55 | self.area = 0.0 |
---|
56 | for j in self.triangle_indices: |
---|
57 | self.area += self.domain.areas[j] |
---|
58 | |
---|
59 | msg = 'Inlet exchange area has area = %f' % self.area |
---|
60 | assert self.area > 0.0 |
---|
61 | |
---|
62 | |
---|
63 | def get_area(self): |
---|
64 | |
---|
65 | return self.area |
---|
66 | |
---|
67 | |
---|
68 | def get_areas(self): |
---|
69 | |
---|
70 | # Must be called after compute_inlet_triangle_indices(). |
---|
71 | return self.domain.areas.take(self.triangle_indices) |
---|
72 | |
---|
73 | |
---|
74 | def get_stages(self): |
---|
75 | |
---|
76 | return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices) |
---|
77 | |
---|
78 | |
---|
79 | def get_average_stage(self): |
---|
80 | |
---|
81 | return num.sum(self.get_stages()*self.get_areas())/self.area |
---|
82 | |
---|
83 | def get_elevations(self): |
---|
84 | |
---|
85 | return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices) |
---|
86 | |
---|
87 | def get_average_elevation(self): |
---|
88 | |
---|
89 | return num.sum(self.get_elevations()*self.get_areas())/self.area |
---|
90 | |
---|
91 | |
---|
92 | def get_xmoms(self): |
---|
93 | |
---|
94 | return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices) |
---|
95 | |
---|
96 | |
---|
97 | def get_average_xmom(self): |
---|
98 | |
---|
99 | return num.sum(self.get_xmoms()*self.get_areas())/self.area |
---|
100 | |
---|
101 | |
---|
102 | def get_ymoms(self): |
---|
103 | |
---|
104 | return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices) |
---|
105 | |
---|
106 | |
---|
107 | def get_average_ymom(self): |
---|
108 | |
---|
109 | return num.sum(self.get_ymoms()*self.get_areas())/self.area |
---|
110 | |
---|
111 | |
---|
112 | def get_heights(self): |
---|
113 | |
---|
114 | return self.get_stages() - self.get_elevations() |
---|
115 | |
---|
116 | |
---|
117 | def get_total_water_volume(self): |
---|
118 | |
---|
119 | return num.sum(self.get_heights()*self.get_areas()) |
---|
120 | |
---|
121 | |
---|
122 | def get_average_height(self): |
---|
123 | |
---|
124 | return self.get_total_water_volume()/self.area |
---|
125 | |
---|
126 | |
---|
127 | def get_velocities(self): |
---|
128 | |
---|
129 | heights = self.get_heights() |
---|
130 | u = self.get_xmoms()/(heights + velocity_protection/heights) |
---|
131 | v = self.get_ymoms()/(heights + velocity_protection/heights) |
---|
132 | |
---|
133 | return u, v |
---|
134 | |
---|
135 | |
---|
136 | def get_xvelocities(self): |
---|
137 | |
---|
138 | heights = self.get_heights() |
---|
139 | return self.get_xmoms()/(heights + velocity_protection/heights) |
---|
140 | |
---|
141 | def get_yvelocities(self): |
---|
142 | |
---|
143 | heights = self.get_heights() |
---|
144 | return self.get_ymoms()/(heights + velocity_protection/heights) |
---|
145 | |
---|
146 | |
---|
147 | def get_average_speed(self): |
---|
148 | |
---|
149 | u, v = self.get_velocities() |
---|
150 | |
---|
151 | average_u = num.sum(u*self.get_areas())/self.area |
---|
152 | average_v = num.sum(v*self.get_areas())/self.area |
---|
153 | |
---|
154 | return math.sqrt(average_u**2 + average_v**2) |
---|
155 | |
---|
156 | |
---|
157 | def get_average_velocity_head(self): |
---|
158 | |
---|
159 | return 0.5*self.get_average_speed()**2/g |
---|
160 | |
---|
161 | |
---|
162 | def get_average_total_energy(self): |
---|
163 | |
---|
164 | return self.get_average_velocity_head() + self.get_average_stage() |
---|
165 | |
---|
166 | |
---|
167 | def get_average_specific_energy(self): |
---|
168 | |
---|
169 | return self.get_average_velocity_head() + self.get_average_height() |
---|
170 | |
---|
171 | |
---|
172 | |
---|
173 | def set_heights(self,height): |
---|
174 | |
---|
175 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height) |
---|
176 | |
---|
177 | |
---|
178 | def set_stages(self,stage): |
---|
179 | |
---|
180 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage) |
---|
181 | |
---|
182 | |
---|
183 | def set_xmoms(self,xmom): |
---|
184 | |
---|
185 | self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom) |
---|
186 | |
---|
187 | |
---|
188 | def set_ymoms(self,ymom): |
---|
189 | |
---|
190 | self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom) |
---|
191 | |
---|
192 | |
---|
193 | def set_elevations(self,elevation): |
---|
194 | |
---|
195 | self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation) |
---|