source: trunk/anuga_core/source/anuga/structures/inlet_enquiry.py @ 8150

Last change on this file since 8150 was 8073, checked in by steve, 14 years ago

Added in a unit test for inlet operator
File size: 3.1 KB
Line 
1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
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    def get_enquiry_position(self):
48
49        return self.domain.get_centroid_coordinates(absolute=True)[self.enquiry_index]
50
51    def get_enquiry_stage(self):
52
53        return self.domain.quantities['stage'].centroid_values[self.enquiry_index]
54
55
56    def get_enquiry_xmom(self):
57
58        return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index]
59
60    def get_enquiry_ymom(self):
61
62        return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index]
63
64
65    def get_enquiry_elevation(self):
66
67        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
68
69    def get_enquiry_depth(self):
70
71        return self.get_enquiry_stage() - self.get_enquiry_elevation()
72
73
74    def get_enquiry_velocity(self):
75
76            depth = self.get_enquiry_depth()
77            u = self.get_enquiry_xmom()/(depth + velocity_protection/depth)
78            v = self.get_enquiry_ymom()/(depth + velocity_protection/depth)
79
80            return u, v
81
82
83    def get_enquiry_xvelocity(self):
84
85            depth = self.get_enquiry_depth()
86            return self.get_enquiry_xmom()/(depth + velocity_protection/depth)
87
88    def get_enquiry_yvelocity(self):
89
90            depth = self.get_enquiry_depth()
91            return self.get_enquiry_ymom()/(depth + velocity_protection/depth)
92
93
94    def get_enquiry_speed(self):
95
96            u, v = self.get_enquiry_velocity()
97
98            return math.sqrt(u**2 + v**2)
99
100
101    def get_enquiry_velocity_head(self):
102
103        return 0.5*self.get_enquiry_speed()**2/g
104
105
106    def get_enquiry_total_energy(self):
107
108        return self.get_enquiry_velocity_head() + self.get_enquiry_stage()
109
110
111    def get_enquiry_specific_energy(self):
112
113        return self.get_enquiry_velocity_head() + self.get_enquiry_depth()
114
115
Note: See TracBrowser for help on using the repository browser.