# source:trunk/anuga_core/source/anuga/structures/inlet.py@8049

Last change on this file since 8049 was 8049, checked in by habili, 13 years ago

Changed from using polygons to polylines. The inlets are now represented by triangles that intersect or contain a polyline - The enquiry gap is now computed as: gap = (self.apron + self.enquiry_gap)*self.culvert_vector

File size: 7.6 KB
Line
1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, polyline_overlap
2from anuga.config import velocity_protection, g
3import math
4
5import numpy as num
6
7class Inlet:
8    """Contains information associated with each inlet
9    """
10
11    def __init__(self, domain, polyline, enquiry_pt,  outward_culvert_vector=None, verbose=False):
12
13        self.domain = domain
14        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
15        self.polyline = polyline
16        self.enquiry_pt = enquiry_pt
17        self.outward_culvert_vector = outward_culvert_vector
18        self.verbose = verbose
19
20        self.compute_indices()
21        self.compute_area()
22
23
24    def compute_indices(self):
25
26        # Get boundary (in absolute coordinates)
27        bounding_polygon = self.domain_bounding_polygon
28        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
29        vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
30
31        # Check that polyline lies within the mesh.
32        for point in self.polyline:
33                msg = 'Point %s ' %  str(point)
34                msg += ' did not fall within the domain boundary.'
35                assert is_inside_polygon(point, bounding_polygon), msg
36
37        point = self.enquiry_pt
38        msg = 'Enquiry Point %s ' %  str(point)
39        msg += ' did not fall within the domain boundary.'
40        assert is_inside_polygon(point, bounding_polygon), msg
41
42        self.triangle_indices = polyline_overlap(vertex_coordinates, self.polyline)
43
44        if len(self.triangle_indices) == 0:
45            msg = 'Inlet polyline=%s ' % (self.polyline)
46            msg += 'No triangles intersecting polyline '
47            raise Exception, msg
48
49        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt)
50
51        if self.enquiry_index in self.triangle_indices:
52            msg = 'Enquiry point %s' % (self.enquiry_pt)
53            msg += 'is in an inlet triangle'
54            raise Exception, msg
55
56    def compute_area(self):
57
58        # Compute inlet area as the sum of areas of triangles identified
59        # by polyline. Must be called after compute_inlet_triangle_indices().
60        if len(self.triangle_indices) == 0:
61            region = 'Inlet polyline=%s' % (self.inlet_polyline)
62            msg = 'No triangles have been identified in region '
63            raise Exception, msg
64
65        self.area = 0.0
66        for j in self.triangle_indices:
67            self.area += self.domain.areas[j]
68
69        msg = 'Inlet exchange area has area = %f' % self.area
70        assert self.area > 0.0
71
72
73    def get_area(self):
74
75        return self.area
76
77
78    def get_areas(self):
79
80        # Must be called after compute_inlet_triangle_indices().
81        return self.domain.areas.take(self.triangle_indices)
82
83
84    def get_stages(self):
85
86        return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices)
87
88
89    def get_average_stage(self):
90
91        return num.sum(self.get_stages()*self.get_areas())/self.area
92
93    def get_elevations(self):
94
95        return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices)
96
97    def get_average_elevation(self):
98
99        return num.sum(self.get_elevations()*self.get_areas())/self.area
100
101
102    def get_xmoms(self):
103
104        return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices)
105
106
107    def get_average_xmom(self):
108
109        return num.sum(self.get_xmoms()*self.get_areas())/self.area
110
111
112    def get_ymoms(self):
113
114        return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices)
115
116
117    def get_average_ymom(self):
118
119        return num.sum(self.get_ymoms()*self.get_areas())/self.area
120
121
122    def get_heights(self):
123
124        return self.get_stages() - self.get_elevations()
125
126
127    def get_total_water_volume(self):
128
129       return num.sum(self.get_heights()*self.get_areas())
130
131
132    def get_average_height(self):
133
134        return self.get_total_water_volume()/self.area
135
136
137    def get_velocities(self):
138
139            heights = self.get_heights()
140            u = self.get_xmoms()/(heights + velocity_protection/heights)
141            v = self.get_ymoms()/(heights + velocity_protection/heights)
142
143            return u, v
144
145
146    def get_xvelocities(self):
147
148            heights = self.get_heights()
149            return self.get_xmoms()/(heights + velocity_protection/heights)
150
151    def get_yvelocities(self):
152
153            heights = self.get_heights()
154            return self.get_ymoms()/(heights + velocity_protection/heights)
155
156
157    def get_average_speed(self):
158
159            u, v = self.get_velocities()
160
161            average_u = num.sum(u*self.get_areas())/self.area
162            average_v = num.sum(v*self.get_areas())/self.area
163
164            return math.sqrt(average_u**2 + average_v**2)
165
166
168
169        return 0.5*self.get_average_speed()**2/g
170
171
172    def get_average_total_energy(self):
173
174        return self.get_average_velocity_head() + self.get_average_stage()
175
176
177    def get_average_specific_energy(self):
178
179        return self.get_average_velocity_head() + self.get_average_height()
180
181
182    def get_enquiry_stage(self):
183
184        return self.domain.quantities['stage'].centroid_values[self.enquiry_index]
185
186
187    def get_enquiry_xmom(self):
188
189        return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index]
190
191    def get_enquiry_ymom(self):
192
193        return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index]
194
195
196    def get_enquiry_elevation(self):
197
198        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
199
200    def get_enquiry_height(self):
201
202        return self.get_enquiry_stage() - self.get_enquiry_elevation()
203
204
205    def get_enquiry_velocity(self):
206
207            height = self.get_enquiry_height()
208            u = self.get_enquiry_xmom()/(height + velocity_protection/height)
209            v = self.get_enquiry_ymom()/(height + velocity_protection/height)
210
211            return u, v
212
213
214    def get_enquiry_xvelocity(self):
215
216            height = self.get_enquiry_height()
217            return self.get_enquiry_xmom()/(height + velocity_protection/height)
218
219    def get_enquiry_yvelocity(self):
220
221            height = self.get_enquiry_height()
222            return self.get_enquiry_ymom()/(height + velocity_protection/height)
223
224
225    def get_enquiry_speed(self):
226
227            u, v = self.get_enquiry_velocity()
228
229            return math.sqrt(u**2 + v**2)
230
231
233
234        return 0.5*self.get_enquiry_speed()**2/g
235
236
237    def get_enquiry_total_energy(self):
238
239        return self.get_enquiry_velocity_head() + self.get_enquiry_stage()
240
241
242    def get_enquiry_specific_energy(self):
243
244        return self.get_enquiry_velocity_head() + self.get_enquiry_height()
245
246
247    def set_heights(self,height):
248
249        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height)
250
251
252    def set_stages(self,stage):
253
254        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
255
256
257    def set_xmoms(self,xmom):
258
259        self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
260
261
262    def set_ymoms(self,ymom):
263
264        self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
265
266
267    def set_elevations(self,elevation):
268
269        self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation)
Note: See TracBrowser for help on using the repository browser.