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

Last change on this file since 7978 was 7978, checked in by habili, 12 years ago

Refactoring of the culvert code.

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