# source:trunk/anuga_core/source/anuga/structures/inlet_enquiry.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: 3.0 KB
Line
1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, polyline_overlap
2from anuga.config import velocity_protection, g
3import math
4
5import numpy as num
6
7import inlet
8
9class Inlet_enquiry(inlet.Inlet):
10    """Contains information associated with each inlet plus an enquiry point
11    """
12
13    def __init__(self, domain, polyline, enquiry_pt,  outward_culvert_vector=None, verbose=False):
14
15
16
17        inlet.Inlet.__init__(self, domain, polyline, verbose)
18
19
20        self.enquiry_pt = enquiry_pt
21        self.outward_culvert_vector = outward_culvert_vector
22
23
24        self.compute_enquiry_index()
25
26
27    def compute_enquiry_index(self):
28
29        # Get boundary (in absolute coordinates)
30        bounding_polygon = self.domain_bounding_polygon
31        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
32        vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
33
34
35        point = self.enquiry_pt
36        msg = 'Enquiry Point %s ' %  str(point)
37        msg += ' did not fall within the domain boundary.'
38        assert is_inside_polygon(point, bounding_polygon), msg
39
40        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt)
41
42        if self.enquiry_index in self.triangle_indices:
43            msg = 'Enquiry point %s' % (self.enquiry_pt)
44            msg += 'is in an inlet triangle'
45            raise Exception, msg
46
47
48
49    def get_enquiry_stage(self):
50
51        return self.domain.quantities['stage'].centroid_values[self.enquiry_index]
52
53
54    def get_enquiry_xmom(self):
55
56        return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index]
57
58    def get_enquiry_ymom(self):
59
60        return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index]
61
62
63    def get_enquiry_elevation(self):
64
65        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
66
67    def get_enquiry_height(self):
68
69        return self.get_enquiry_stage() - self.get_enquiry_elevation()
70
71
72    def get_enquiry_velocity(self):
73
74            height = self.get_enquiry_height()
75            u = self.get_enquiry_xmom()/(height + velocity_protection/height)
76            v = self.get_enquiry_ymom()/(height + velocity_protection/height)
77
78            return u, v
79
80
81    def get_enquiry_xvelocity(self):
82
83            height = self.get_enquiry_height()
84            return self.get_enquiry_xmom()/(height + velocity_protection/height)
85
86    def get_enquiry_yvelocity(self):
87
88            height = self.get_enquiry_height()
89            return self.get_enquiry_ymom()/(height + velocity_protection/height)
90
91
92    def get_enquiry_speed(self):
93
94            u, v = self.get_enquiry_velocity()
95
96            return math.sqrt(u**2 + v**2)
97
98
100
101        return 0.5*self.get_enquiry_speed()**2/g
102
103
104    def get_enquiry_total_energy(self):
105
106        return self.get_enquiry_velocity_head() + self.get_enquiry_stage()
107
108
109    def get_enquiry_specific_energy(self):
110
111        return self.get_enquiry_velocity_head() + self.get_enquiry_height()
112
113
Note: See TracBrowser for help on using the repository browser.