Changeset 7958
 Timestamp:
 Aug 20, 2010, 2:02:43 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7957 r7958 2 2 3 3 from anuga.shallow_water.forcing import Inflow, General_forcing 4 from anuga.structures.culvert_polygons import create_culvert_polygons5 4 from anuga.utilities.system_tools import log_to_file 6 5 from anuga.geometry.polygon import inside_polygon, is_inside_polygon … … 17 16 import numpy as num 18 17 from math import sqrt 19 from math import sqrt20 18 21 19 class Below_interval(Exception): pass … … 25 23 26 24 class Generic_box_culvert: 27 """Culvert flow  transfer water from one rect ngular box to another.25 """Culvert flow  transfer water from one rectangular box to another. 28 26 Sets up the geometry of problem 29 27 … … 46 44 # Input check 47 45 46 self.domain = domain 47 48 self.end_points= [end_point0, end_point1] 49 self.enquiry_gap_factor = enquiry_gap_factor 50 48 51 if height is None: 49 52 height = width 50 53 54 self.width = width 51 55 self.height = height 52 self.width = width 53 54 self.domain = domain 55 self.end_points= [end_point0, end_point1] 56 self.enquiry_gap_factor=enquiry_gap_factor 57 56 58 57 self.verbose=verbose 59 58 self.filename = None 60 59 61 62 60 # Create the fundamental culvert polygons from polygon 63 61 self.create_culvert_polygons() … … 65 63 self.check_culvert_inside_domain() 66 64 67 68 65 # Establish initial values at each enquiry point 69 66 self.enquiry_quantity_values = [] 70 67 dq = domain.quantities 71 for i in [0, 1]:68 for i in [0, 1]: 72 69 idx = self.enquiry_indices[i] 73 70 elevation = dq['elevation'].get_values(location='centroids', … … 75 72 stage = dq['stage'].get_values(location='centroids', 76 73 indices=[idx])[0] 77 depth = stage elevation74 depth = stage  elevation 78 75 79 76 quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth } 80 77 self.enquiry_quantity_values.append(quantity_values) 81 82 83 assert self.culvert_length > 0.0 84 85 86 def set_store_hydrograph_discharge(self,filename=None): 78 79 def set_store_hydrograph_discharge(self, filename=None): 87 80 88 81 if filename is None: … … 97 90 fid.write('time, discharge\n') 98 91 fid.close() 99 100 101 102 92 103 93 def create_culvert_polygons(self): … … 111 101 x1, y1 = self.end_points[1] 112 102 113 dx = x1 x0114 dy = y1 y0103 dx = x1  x0 104 dy = y1  y0 115 105 116 106 self.culvert_vector = num.array([dx, dy]) 117 107 self.culvert_length = sqrt(num.sum(self.culvert_vector**2)) 118 108 assert self.culvert_length > 0.0, 'The length of culvert is less than 0' 119 109 120 110 # Unit direction vector and normal … … 123 113 124 114 # Short hands 125 w = 0.5* width*normal # Perpendicular vector of 1/2 width126 h = height*vector # Vector of length=height in the115 w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width 116 h = self.height*self.culvert_vector # Vector of length=height in the 127 117 # direction of the culvert 128 gap = (1 + enquiry_gap_factor)*h 129 118 gap = (1 + self.enquiry_gap_factor)*h 130 119 131 120 self.exchange_polygons = [] … … 133 122 134 123 # Build exchange polygon and enquiry points 0 and 1 135 for i in [0, 1]:124 for i in [0, 1]: 136 125 p0 = self.end_points[i] + w 137 p1 = self.end_point [i]  w126 p1 = self.end_points[i]  w 138 127 p2 = p1  h 139 128 p3 = p0  h 140 self.exchange_polygons.append(num.array([p0, p1,p2,p3]))141 self.enquiry_points.append( end_point0 gap)129 self.exchange_polygons.append(num.array([p0, p1, p2, p3])) 130 self.enquiry_points.append(self.end_points[i] + (2*i1)*gap) 142 131 143 132 self.polygon_areas = [] … … 239 228 """ 240 229 bounding_polygon = self.domain.get_boundary_polygon() 241 P = self.culvert_polygons 242 for key in P.keys(): 243 if key in ['exchange_polygon0', 244 'exchange_polygon1']: 245 for point in list(P[key]) + self.enquiry_points: 246 msg = 'Point %s in polygon %s for culvert %s did not'\ 247 %(str(point), key, self.label) 248 msg += 'fall within the domain boundary.' 249 assert is_inside_polygon(point, bounding_polygon), msg 230 for i in [0, 1]: 231 for point in list(self.exchange_polygons[i]) + self.enquiry_points: 232 msg = 'Point %s did not '\ 233 %(str(point)) 234 msg += 'fall within the domain boundary.' 235 assert is_inside_polygon(point, bounding_polygon), msg 250 236 251 237 … … 719 705 return label, type, width, height, length, number_of_barrels, description, rating_curve 720 706 721 722 723 Culvert_flow = Culvert_flow_general
Note: See TracChangeset
for help on using the changeset viewer.