 Timestamp:
 Aug 26, 2010, 5:17:24 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7976 r7977 2 2 from anuga.config import g 3 3 import anuga.utilities.log as log 4 import inlet 5 import numpy as num 6 import math 4 import box_culvert 7 5 8 class Below_interval(Exception): pass 9 class Above_interval(Exception): pass 10 11 12 class Generic_box_culvert: 6 class Culvert_operator: 13 7 """Culvert flow  transfer water from one rectangular box to another. 14 8 Sets up the geometry of problem … … 29 23 verbose=False): 30 24 31 # Input check32 33 25 self.domain = domain 34 35 26 self.domain.set_fractional_step_operator(self) 36 37 self.end_points = [end_point0, end_point1] 27 end_points = [end_point0, end_point1] 38 28 39 29 if height is None: … … 42 32 self.width = width 43 33 self.height = height 44 45 self.verbose=verbose46 34 47 # Create the fundamental culvert polygons and create inlet objects 48 self.create_culvert_polygons() 49 50 #FIXME (SR) Put this into a foe loop to deal with more inlets 51 self.inlets = [] 52 polygon0 = self.inlet_polygons[0] 53 inlet0_vector = self.culvert_vector 54 self.inlets.append(inlet.Inlet(self.domain, polygon0)) 55 56 polygon1 = self.inlet_polygons[1] 57 inlet1_vector =  self.culvert_vector 58 self.inlets.append(inlet.Inlet(self.domain, polygon1)) 35 self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height) 36 self.inlets = self.culvert.get_inlets() 59 37 60 38 self.print_stats() … … 63 41 def __call__(self): 64 42 65 # Time stuff66 time = self.domain.get_time()67 43 timestep = self.domain.get_timestep() 68 69 70 44 71 45 inflow = self.inlets[0] … … 81 55 82 56 delta_z = inflow.get_average_elevation()  outflow.get_average_elevation() 83 culvert_slope = delta_z/self.culvert _length57 culvert_slope = delta_z/self.culvert.get_culvert_length() 84 58 85 59 # Determine controlling energy (driving head) for culvert … … 91 65 driving_head = inflow.get_average_specific_energy() 92 66 93 94 67 # Transfer 95 from culvert_routines import boyd_ generalised_culvert_model68 from culvert_routines import boyd_box, boyd_circle 96 69 Q, barrel_velocity, culvert_outlet_depth =\ 97 boyd_generalised_culvert_model(inflow.get_average_height(),70 boyd_circle(inflow.get_average_height(), 98 71 outflow.get_average_height(), 99 72 inflow.get_average_speed(), … … 101 74 inflow.get_average_specific_energy(), 102 75 delta_total_energy, 103 g, 104 culvert_length=self.culvert_length, 76 culvert_length=self.culvert.get_culvert_length(), 105 77 culvert_width=self.width, 106 78 culvert_height=self.height, 107 culvert_type='box',108 79 manning=0.01) 109 80 110 81 transfer_water = Q*timestep 111 112 82 113 83 inflow.set_heights(inflow.get_average_height()  transfer_water) … … 115 85 inflow.set_ymoms(0.0) 116 86 117 118 87 outflow.set_heights(outflow.get_average_height() + transfer_water) 119 88 outflow.set_xmoms(0.0) 120 89 outflow.set_ymoms(0.0) 121 122 123 90 124 91 def print_stats(self): … … 141 108 142 109 print '=====================================' 143 144 145 146 147 148 def create_culvert_polygons(self):149 150 """Create polygons at the end of a culvert inlet and outlet.151 At either end two polygons will be created; one for the actual flow to pass through and one a little further away152 for enquiring the total energy at both ends of the culvert and transferring flow.153 """154 155 # Calculate geometry156 x0, y0 = self.end_points[0]157 x1, y1 = self.end_points[1]158 159 dx = x1  x0160 dy = y1  y0161 162 self.culvert_vector = num.array([dx, dy])163 self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))164 assert self.culvert_length > 0.0, 'The length of culvert is less than 0'165 166 # Unit direction vector and normal167 self.culvert_vector /= self.culvert_length # Unit vector in culvert direction168 self.culvert_normal = num.array([dy, dx])/self.culvert_length # Normal vector169 170 # Short hands171 w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width172 h = self.height*self.culvert_vector # Vector of length=height in the173 # direction of the culvert174 175 self.inlet_polygons = []176 177 # Build exchange polygon and enquiry points 0 and 1178 for i in [0, 1]:179 i0 = (2*i1)180 p0 = self.end_points[i] + w181 p1 = self.end_points[i]  w182 p2 = p1 + i0*h183 p3 = p0 + i0*h184 self.inlet_polygons.append(num.array([p0, p1, p2, p3]))185 186 # Check that enquiry points are outside inlet polygons187 for i in [0,1]:188 polygon = self.inlet_polygons[i]189 # FIXME (SR) Probably should calculate the area of all the triangles190 # associated with this polygon, as there is likely to be some191 # inconsistency between triangles and ploygon192 area = polygon_area(polygon)193 194 msg = 'Polygon %s ' %(polygon)195 msg += ' has area = %f' % area196 assert area > 0.0, msg197 198 110 199 111
Note: See TracChangeset
for help on using the changeset viewer.