# source:trunk/anuga_core/source/anuga/structures/box_culvert.py@7978

Last change on this file since 7978 was 7978, checked in by habili, 13 years ago

Refactoring of the culvert code.

File size: 3.7 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 Box_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=None,
23                 height=None,
24                 verbose=False):
25
26        # Input check
27
28        self.domain = domain
29
30        self.end_points = end_points
31
32        self.width = width
33        self.height = height
34
35        self.verbose=verbose
36
37        # Create the fundamental culvert polygons and create inlet objects
38        self.__create_culvert_polygons()
39
40        #FIXME (SR) Put this into a foe loop to deal with more inlets
41        self.inlets = []
42        polygon0 = self.inlet_polygons
43        inlet0_vector = self.culvert_vector
44        self.inlets.append(inlet.Inlet(self.domain, polygon0))
45
46        polygon1 = self.inlet_polygons
47        inlet1_vector = - self.culvert_vector
48        self.inlets.append(inlet.Inlet(self.domain, polygon1))
49
50
51    def __create_culvert_polygons(self):
52
53        """Create polygons at the end of a culvert inlet and outlet.
54        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
55        for enquiring the total energy at both ends of the culvert and transferring flow.
56        """
57
58        # Calculate geometry
59        x0, y0 = self.end_points
60        x1, y1 = self.end_points
61
62        dx = x1 - x0
63        dy = y1 - y0
64
65        self.culvert_vector = num.array([dx, dy])
66        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
67        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
68
69        # Unit direction vector and normal
70        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
71        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
72
73        # Short hands
74        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
75        h = self.height*self.culvert_vector    # Vector of length=height in the
76                             # direction of the culvert
77
78        self.inlet_polygons = []
79
80        # Build exchange polygon and enquiry points 0 and 1
81        for i in [0, 1]:
82            i0 = (2*i-1)
83            p0 = self.end_points[i] + w
84            p1 = self.end_points[i] - w
85            p2 = p1 + i0*h
86            p3 = p0 + i0*h
87            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
88
89        # Check that enquiry points are outside inlet polygons
90        for i in [0,1]:
91            polygon = self.inlet_polygons[i]
92            # FIXME (SR) Probably should calculate the area of all the triangles
93            # associated with this polygon, as there is likely to be some
94            # inconsistency between triangles and ploygon
95            area = polygon_area(polygon)
96
97            msg = 'Polygon %s ' %(polygon)
98            msg += ' has area = %f' % area
99            assert area > 0.0, msg
100
101
102    def get_inlets(self):
103
104        return self.inlets
105
106
107    def get_culvert_length(self):
108
109        return self.culvert_length
110
111
112    def get_culvert_width(self):
113
114        return self.width
115
116
117    def get_culvert_height(self):
118
119        return self.height
120
Note: See TracBrowser for help on using the repository browser.