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

Last change on this file since 8049 was 8049, checked in by habili, 12 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
167    def get_average_velocity_head(self):
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
232    def get_enquiry_velocity_head(self):
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.