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

Last change on this file since 7995 was 7995, checked in by steve, 13 years ago

Added stats

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