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

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

using "import anuga"

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