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

Last change on this file since 7986 was 7986, checked in by habili, 14 years ago

deal with manning

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