source: trunk/anuga_core/source/anuga/structures/structure_operator.py @ 7993

Last change on this file since 7993 was 7993, checked in by habili, 12 years ago

Consolidation of the culvert routines into structure classes

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