source: trunk/anuga_core/source/anuga/structures/Inlet.py @ 7972

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

Contains the Inlet class

File size: 4.4 KB
Line 
1from anuga.geometry.polygon import inside_polygon, is_inside_polygon
2from anuga.config import velocity_protection
3
4import numpy as num
5
6class Inlet:
7    """Contains information associated with each inlet
8    """
9
10    def __init__(self, domain, polygon, enquiry_point, inlet_vector):
11
12        self.domain = domain
13        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
14        self.polygon = polygon
15        self.enquiry_point = enquiry_point
16        self.inlet_vector = inlet_vector
17
18        # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
19        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_point)
20
21        self.compute_triangle_indices()
22        self.compute_area()
23
24
25    def compute_triangle_indices(self):
26
27        # Get boundary (in absolute coordinates)
28        bounding_polygon = self.domain_bounding_polygon
29        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
30
31        self.inlet_polygon = self.polygon
32
33        # Check that polygon lies within the mesh.
34        for point in self.inlet_polygon:
35                msg = 'Point %s in polygon for forcing term' %  str(point)
36                msg += ' did not fall within the domain boundary.'
37                assert is_inside_polygon(point, bounding_polygon), msg
38
39        self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon)
40
41        if len(self.triangle_indices) == 0:
42            region = 'Inlet polygon=%s' % (self.inlet_polygon)
43            msg = 'No triangles have been identified in '
44            msg += 'specified region: %s' % region
45            raise Exception, msg
46
47
48    def compute_area(self):
49       
50        # Compute inlet area as the sum of areas of triangles identified
51        # by polygon. Must be called after compute_inlet_triangle_indices().
52        if len(self.triangle_indices) == 0:
53            region = 'Inlet polygon=%s' % (self.inlet_polygon)
54            msg = 'No triangles have been identified in '
55            msg += 'specified region: %s' % region
56            raise Exception, msg
57       
58        self.area = 0.0
59        for j in self.triangle_indices:
60            self.area += self.domain.areas[j]
61
62        msg = 'Inlet exchange area has area = %f' % self.area
63        assert self.area > 0.0
64
65
66    def get_areas(self):
67       
68        # Must be called after compute_inlet_triangle_indices().
69        return self.domain.areas.take(self.triangle_indices)
70   
71       
72    def get_stages(self):
73       
74        return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices)
75       
76    def get_average_stage(self):
77       
78        return num.sum(self.get_stages())/self.triangle_indices.size
79       
80       
81    def get_elevations(self):   
82       
83        return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices)
84       
85    def get_average_elevation(self):
86       
87        return num.sum(self.get_elevations())/self.triangle_indices.size
88   
89   
90    def get_xmoms(self):
91   
92        return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices)
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_velocities(self):
134 
135            depths = self.get_stages() - self.get_elevations()
136            u = self.get_xmoms()/(depths + velocity_protection/depths)
137            v = self.get_ymoms()/(depths + velocity_protection/depths)
138           
139            return num.sum(u)/self.triangle_indices.size, num.sum(v)/self.triangle_indices.size
Note: See TracBrowser for help on using the repository browser.