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

Last change on this file since 7975 was 7975, checked in by steve, 12 years ago

Got a working culvert!

File size: 5.4 KB
Line 
1
2from anuga.geometry.polygon import inside_polygon, is_inside_polygon
3from anuga.config import velocity_protection, g
4import math
5
6import numpy as num
7
8class Inlet:
9    """Contains information associated with each inlet
10    """
11
12    def __init__(self, domain, polygon, enquiry_point, inlet_vector):
13
14        self.domain = domain
15        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
16        self.polygon = polygon
17        self.enquiry_point = enquiry_point
18        self.inlet_vector = inlet_vector
19
20        # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
21        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_point)
22
23        self.compute_triangle_indices()
24        self.compute_area()
25
26
27    def compute_triangle_indices(self):
28
29        # Get boundary (in absolute coordinates)
30        bounding_polygon = self.domain_bounding_polygon
31        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
32
33        self.inlet_polygon = self.polygon
34
35        # Check that polygon lies within the mesh.
36        for point in self.inlet_polygon:
37                msg = 'Point %s in polygon for forcing term' %  str(point)
38                msg += ' did not fall within the domain boundary.'
39                assert is_inside_polygon(point, bounding_polygon), msg
40
41        self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon)
42
43        if len(self.triangle_indices) == 0:
44            region = 'Inlet polygon=%s' % (self.inlet_polygon)
45            msg = 'No triangles have been identified in '
46            msg += 'specified region: %s' % region
47            raise Exception, msg
48
49
50    def compute_area(self):
51       
52        # Compute inlet area as the sum of areas of triangles identified
53        # by polygon. Must be called after compute_inlet_triangle_indices().
54        if len(self.triangle_indices) == 0:
55            region = 'Inlet polygon=%s' % (self.inlet_polygon)
56            msg = 'No triangles have been identified in '
57            msg += 'specified region: %s' % 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_areas(self):
69       
70        # Must be called after compute_inlet_triangle_indices().
71        return self.domain.areas.take(self.triangle_indices)
72   
73       
74    def get_stages(self):
75       
76        return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices)
77       
78       
79    def get_average_stage(self):
80       
81        return num.sum(self.get_stages())/self.triangle_indices.size
82       
83       
84    def get_elevations(self):   
85       
86        return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices)
87       
88    def get_average_elevation(self):
89       
90        return num.sum(self.get_elevations())/self.triangle_indices.size
91   
92   
93    def get_xmoms(self):
94   
95        return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices)
96       
97       
98    def get_average_xmom(self):
99       
100        return num.sum(self.get_xmoms())/self.triangle_indices.size
101       
102   
103    def get_ymoms(self):
104       
105        return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices)
106 
107 
108    def get_average_ymom(self):
109       
110        return num.sum(self.get_ymoms())/self.triangle_indices.size
111 
112   
113    def get_heights(self):
114   
115        return self.get_stages() - self.get_elevations()
116   
117   
118    def get_total_water_volume(self):
119       
120       return num.sum(self.get_heights()*self.get_areas())
121   
122   
123    def get_average_height(self):
124   
125        return self.get_total_water_volume()/self.area
126       
127       
128    def get_velocities(self):
129       
130            depths = self.get_stages() - self.get_elevations()
131            u = self.get_xmoms()/(depths + velocity_protection/depths)
132            v = self.get_ymoms()/(depths + velocity_protection/depths)
133           
134            return u, v
135           
136           
137    def get_average_speed(self):
138 
139            u, v = self.get_velocities()
140           
141            average_u = num.sum(u)/self.triangle_indices.size
142            average_v = num.sum(v)/self.triangle_indices.size
143           
144            return math.sqrt(average_u**2 + average_v**2)
145
146
147
148    def get_average_velocity_head(self):
149
150        return 0.5*self.get_average_speed()**2/g
151
152
153    def get_average_total_energy(self):
154       
155        return self.get_average_velocity_head() + self.get_average_stage()
156       
157   
158    def get_average_specific_energy(self):
159       
160        return self.get_average_velocity_head() + self.get_average_height()
161
162
163
164    def set_heights(self,height):
165
166        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height)
167
168
169    def set_stages(self,stage):
170
171        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
172
173    def set_xmoms(self,xmom):
174
175        self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
176
177
178    def set_ymoms(self,ymom):
179
180        self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
181
182
183    def set_elevations(self,elevation):
184
185        self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation)
Note: See TracBrowser for help on using the repository browser.