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

Last change on this file since 8007 was 8003, checked in by habili, 14 years ago

fix for is_inside_polygon

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