# source:trunk/anuga_core/source/anuga/structures/inlet.py@7977

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

Refactoring of routines.

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