source: trunk/anuga_core/source/anuga/structures/culvert.py @ 7991

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

Finally resolved conflicts between nariman and steve

File size: 4.4 KB
Line 
1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
3import anuga.utilities.log as log
4import inlet
5import numpy as num
6import math
7
8class Culvert:
9    """Culvert flow - transfer water from one rectangular box to another.
10    Sets up the geometry of problem
11   
12    This is the base class for culverts. Inherit from this class (and overwrite
13    compute_discharge method for specific subclasses)
14   
15    Input: Two points, pipe_size (either diameter or width, height),
16    mannings_rougness,
17    """ 
18
19    def __init__(self,
20                 domain,
21                 end_points, 
22                 width,
23                 height,
24                 apron,
25                 manning,
26                 enquiry_gap,
27                 use_velocity_head,
28                 verbose):
29       
30        # Input check
31       
32        self.domain = domain
33        self.end_points = end_points
34        self.width  = width
35        self.height = height
36        self.apron  = apron
37        self.manning = manning
38       
39        self.enquiry_gap = enquiry_gap
40        self.use_velocity_head = use_velocity_head
41        self.verbose=verbose
42       
43        # Create the fundamental culvert polygons and create inlet objects
44        self.__create_culvert_polygons()
45
46        #FIXME (SR) Put this into a foe loop to deal with more inlets
47
48        self.inlets = []
49        polygon0 = self.inlet_polygons[0]
50        ep0 = self.inlet_equiry_pts[0]
51        outward_vector0 = self.culvert_vector
52        self.inlets.append(inlet.Inlet(self.domain, polygon0, ep0, outward_vector0))
53
54        polygon1 = self.inlet_polygons[1]
55        ep1 = self.inlet_equiry_pts[1]
56        outward_vector1  = - self.culvert_vector
57        self.inlets.append(inlet.Inlet(self.domain, polygon1, ep1, outward_vector1))
58
59
60    def __call__(self):
61        msg = 'Method __call__ must be overridden by Culvert subclass'
62        raise Exception, msg
63
64    def __create_culvert_polygons(self):
65
66        """Create polygons at the end of a culvert inlet and outlet.
67        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
68        for enquiring the total energy at both ends of the culvert and transferring flow.
69        """
70
71        # Calculate geometry
72        x0, y0 = self.end_points[0]
73        x1, y1 = self.end_points[1]
74
75        dx = x1 - x0
76        dy = y1 - y0
77
78        self.culvert_vector = num.array([dx, dy])
79        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
80        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
81
82        # Unit direction vector and normal
83        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
84        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
85
86        # Short hands
87        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
88        h = self.apron*self.culvert_vector    # Vector of length=height in the
89                             # direction of the culvert
90
91        gap = (1 + self.enquiry_gap)*h
92
93        self.inlet_polygons = []
94        self.inlet_equiry_pts = []
95
96        # Build exchange polygon and enquiry point
97        for i in [0, 1]:
98            i0 = (2*i-1)
99            p0 = self.end_points[i] + w
100            p1 = self.end_points[i] - w
101            p2 = p1 + i0*h
102            p3 = p0 + i0*h
103            ep = self.end_points[i] + i0*gap
104
105            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
106            self.inlet_equiry_pts.append(ep)
107
108        # Check that enquiry points are outside inlet polygons
109        for i in [0,1]:
110            polygon = self.inlet_polygons[i]
111            ep = self.inlet_equiry_pts[i]
112           
113            area = polygon_area(polygon)
114           
115            msg = 'Polygon %s ' %(polygon)
116            msg += ' has area = %f' % area
117            assert area > 0.0, msg
118
119            msg = 'Enquiry point falls inside an exchange polygon.'
120            assert not inside_polygon(ep, polygon), msg
121   
122    def get_inlets(self):
123       
124        return self.inlets
125       
126       
127    def get_culvert_length(self):
128       
129        return self.culvert_length
130       
131       
132    def get_culvert_width(self):
133       
134        return self.width
135       
136       
137    def get_culvert_height(self):
138   
139        return self.height
140
141
142    def get_culvert_apron(self):
143
144        return self.apron
145       
Note: See TracBrowser for help on using the repository browser.