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

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

Confirming if the correct Culvert Velocity is being used, that is the resultant of Delta_Total_Energy

File size: 6.6 KB
RevLine 
[7993]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:
[7995]8    """Structure Operator - transfer water from one rectangular box to another.
[7993]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,
[7996]27                 description,
[7993]28                 verbose):
29       
30        self.domain = domain
31        self.domain.set_fractional_step_operator(self)
32        self.end_points = [end_point0, end_point1]
33       
34        if height is None:
35            height = width
36
37        if apron is None:
38            apron = width
39
40        self.width  = width
41        self.height = height
42        self.apron  = apron
43        self.manning = manning
44        self.enquiry_gap = enquiry_gap
[7996]45        self.description = description
[7993]46        self.verbose = verbose
[7995]47
48        self.discharge = 0.0
49        self.velocity = 0.0
[7996]50        self.delta_total_energy = 0.0
51        self.driving_energy = 0.0
[7993]52       
53        self.__create_exchange_polygons()
54
55        self.inlets = []
56        polygon0 = self.inlet_polygons[0]
57        enquiry_point0 = self.inlet_equiry_points[0]
58        outward_vector0 = self.culvert_vector
59        self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_point0, outward_vector0))
60
61        polygon1 = self.inlet_polygons[1]
62        exchange_polygon1 = self.inlet_equiry_points[1]
63        outward_vector1  = - self.culvert_vector
64        self.inlets.append(inlet.Inlet(self.domain, polygon1, exchange_polygon1, outward_vector1))
65
66    def __call__(self):
67
68        pass
69
70    def __create_exchange_polygons(self):
71
72        """Create polygons at the end of a culvert inlet and outlet.
73        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
74        for enquiring the total energy at both ends of the culvert and transferring flow.
75        """
76
77        # Calculate geometry
78        x0, y0 = self.end_points[0]
79        x1, y1 = self.end_points[1]
80
81        dx = x1 - x0
82        dy = y1 - y0
83
84        self.culvert_vector = num.array([dx, dy])
85        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
86        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
87
88        # Unit direction vector and normal
89        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
90        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
91
92        # Short hands
93        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
94        h = self.apron*self.culvert_vector    # Vector of length=height in the
95                             # direction of the culvert
96
97        gap = (1 + self.enquiry_gap)*h
98
99        self.inlet_polygons = []
100        self.inlet_equiry_points = []
101
102        # Build exchange polygon and enquiry point
103        for i in [0, 1]:
104            i0 = (2*i-1)
105            p0 = self.end_points[i] + w
106            p1 = self.end_points[i] - w
107            p2 = p1 + i0*h
108            p3 = p0 + i0*h
109            ep = self.end_points[i] + i0*gap
110
111            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
112            self.inlet_equiry_points.append(ep)
113
114        # Check that enquiry points are outside inlet polygons
115        for i in [0,1]:
116            polygon = self.inlet_polygons[i]
117            ep = self.inlet_equiry_points[i]
118           
119            area = polygon_area(polygon)
120           
121            msg = 'Polygon %s ' %(polygon)
122            msg += ' has area = %f' % area
123            assert area > 0.0, msg
124
125            msg = 'Enquiry point falls inside an exchange polygon.'
126            assert not inside_polygon(ep, polygon), msg
127   
128           
129        #print '   outflow volume ',outflow.get_total_water_volume()
130       
131
132    def print_stats(self):
133
134        print '====================================='
135        print 'Generic Culvert Operator'
136        print '====================================='
137
138        print 'Culvert'
139        print self.culvert
140
141        print 'Culvert Routine'
142        print self.routine
143       
144        for i, inlet in enumerate(self.inlets):
145            print '-------------------------------------'
146            print 'Inlet %i' % i
147            print '-------------------------------------'
148
149            print 'inlet triangle indices and centres'
150            print inlet.triangle_indices
151            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
152       
153            print 'polygon'
154            print inlet.polygon
155
156        print '====================================='
157
[7995]158
159    def structure_statistics(self):
160
161        message = '---------------------------\n'
[7996]162        message += 'Structure report for structure %s:\n' % self.description
[7995]163        message += '--------------------------\n'
164        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
165        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
[7996]166        message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
167        message += 'delta total energy %.2f\n' % self.delta_total_energy
[7995]168#        message += 'Net boundary flow by tags [m^3/s]\n'
169#        for tag in boundary_flows:
170#            message += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
171#
172#        message += 'Total net boundary flow [m^3/s]: %.2f\n' % \
173#                    (total_boundary_inflow + total_boundary_outflow)
174#        message += 'Total volume in domain [m^3]: %.2f\n' % \
175#                    self.compute_total_volume()
176#
177#        # The go through explicit forcing update and record the rate of change
178#        # for stage and
179#        # record into forcing_inflow and forcing_outflow. Finally compute
180#        # integral of depth to obtain total volume of domain.
181#
182        # FIXME(Ole): This part is not yet done.
183
184        return message
185
[7993]186    def get_inlets(self):
187       
188        return self.inlets
189       
190       
191    def get_culvert_length(self):
192       
193        return self.culvert_length
194       
195       
196    def get_culvert_width(self):
197       
198        return self.width
199       
200       
201    def get_culvert_height(self):
202   
203        return self.height
204
205
206    def get_culvert_apron(self):
207
208        return self.apron
Note: See TracBrowser for help on using the repository browser.