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

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

Uses polygon_overlap method to find triangles that overlap the inlet.

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