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

Last change on this file since 7984 was 7984, checked in by steve, 13 years ago

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