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

Last change on this file since 8272 was 8250, checked in by steve, 13 years ago

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