1 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon |
---|
2 | from anuga.config import velocity_protection |
---|
3 | |
---|
4 | import numpy as num |
---|
5 | |
---|
6 | class 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 |
---|