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

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

Added in an inlet_operator class to implement inflow "boundary condition" by setting up a inlet defined by a line close to the boundary. Look at testing_culvert_inlet.py as an example.

File size: 5.4 KB
RevLine
[8049]1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, polyline_overlap
[7975]2from anuga.config import velocity_protection, g
[7973]3import math
[7972]4
5import numpy as num
6
7class Inlet:
8    """Contains information associated with each inlet
9    """
10
[8050]11    def __init__(self, domain, polyline, verbose=False):
[7972]12
13        self.domain = domain
14        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
[8048]15        self.polyline = polyline
[7984]16        self.verbose = verbose
[7972]17
[8050]18        self.compute_triangle_indices()
[7972]19        self.compute_area()
20
21
[8050]22    def compute_triangle_indices(self):
[7972]23
24        # Get boundary (in absolute coordinates)
25        bounding_polygon = self.domain_bounding_polygon
26        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
[8042]27        vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
[7972]28
[8048]29        # Check that polyline lies within the mesh.
30        for point in self.polyline:
[7984]31                msg = 'Point %s ' %  str(point)
[7972]32                msg += ' did not fall within the domain boundary.'
33                assert is_inside_polygon(point, bounding_polygon), msg
[8003]34
[7972]35
[8050]36
[8048]37        self.triangle_indices = polyline_overlap(vertex_coordinates, self.polyline)
[7972]38
39        if len(self.triangle_indices) == 0:
[8049]40            msg = 'Inlet polyline=%s ' % (self.polyline)
41            msg += 'No triangles intersecting polyline '
[7972]42            raise Exception, msg
43
[7980]44
[8050]45
[7972]46    def compute_area(self):
47
48        # Compute inlet area as the sum of areas of triangles identified
[8048]49        # by polyline. Must be called after compute_inlet_triangle_indices().
[7972]50        if len(self.triangle_indices) == 0:
[8048]51            region = 'Inlet polyline=%s' % (self.inlet_polyline)
[7984]52            msg = 'No triangles have been identified in region '
[7972]53            raise Exception, msg
54
55        self.area = 0.0
56        for j in self.triangle_indices:
57            self.area += self.domain.areas[j]
58
59        msg = 'Inlet exchange area has area = %f' % self.area
60        assert self.area > 0.0
61
62
[7980]63    def get_area(self):
64
65        return self.area
66
67
[7972]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
[7973]78
[7972]79    def get_average_stage(self):
[7980]80
81        return num.sum(self.get_stages()*self.get_areas())/self.area
[7972]82
83    def get_elevations(self):
84
85        return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices)
86
87    def get_average_elevation(self):
[7980]88
89        return num.sum(self.get_elevations()*self.get_areas())/self.area
[7972]90
91
92    def get_xmoms(self):
93
94        return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices)
95
[7973]96
[7972]97    def get_average_xmom(self):
[7980]98
99        return num.sum(self.get_xmoms()*self.get_areas())/self.area
[7972]100
101
102    def get_ymoms(self):
103
104        return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices)
105
106
107    def get_average_ymom(self):
108
[7980]109        return num.sum(self.get_ymoms()*self.get_areas())/self.area
[7972]110
[7984]111
[7972]112    def get_heights(self):
113
114        return self.get_stages() - self.get_elevations()
115
116
117    def get_total_water_volume(self):
118
[7975]119       return num.sum(self.get_heights()*self.get_areas())
[7980]120
121
[7972]122    def get_average_height(self):
123
124        return self.get_total_water_volume()/self.area
125
126
127    def get_velocities(self):
128
[7980]129            heights = self.get_heights()
130            u = self.get_xmoms()/(heights + velocity_protection/heights)
131            v = self.get_ymoms()/(heights + velocity_protection/heights)
[7972]132
133            return u, v
[7980]134
135
136    def get_xvelocities(self):
137
138            heights = self.get_heights()
139            return self.get_xmoms()/(heights + velocity_protection/heights)
140
141    def get_yvelocities(self):
142
143            heights = self.get_heights()
144            return self.get_ymoms()/(heights + velocity_protection/heights)
[7972]145
146
[7975]147    def get_average_speed(self):
[7972]148
[7973]149            u, v = self.get_velocities()
[7972]150
[7980]151            average_u = num.sum(u*self.get_areas())/self.area
152            average_v = num.sum(v*self.get_areas())/self.area
[7973]153
[7975]154            return math.sqrt(average_u**2 + average_v**2)
[7973]155
156
158
159        return 0.5*self.get_average_speed()**2/g
160
161
162    def get_average_total_energy(self):
[7973]163
[7973]165
[7975]166
167    def get_average_specific_energy(self):
[7973]168
[7973]170
[7975]171
[7984]172
[7975]173    def set_heights(self,height):
174
175        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height)
176
177
178    def set_stages(self,stage):
179
180        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
181
[7977]182
[7975]183    def set_xmoms(self,xmom):
184
[7984]185        self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
[7975]186
187
188    def set_ymoms(self,ymom):
189
[7984]190        self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
[7975]191
192
193    def set_elevations(self,elevation):
194
[7976]195        self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation)
Note: See TracBrowser for help on using the repository browser.