- Timestamp:
- Jan 7, 2009, 4:34:45 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/culvert_flows/culvert_polygons.py
r6098 r6121 5 5 from math import sqrt 6 6 from Numeric import array, sum 7 from anuga.utilities.polygon import inside_polygon, polygon_area 7 8 8 9 def create_culvert_polygons(end_point0, 9 10 end_point1, 10 11 width, height=None, 11 enquiry_gap_factor=1.0, 12 enquiry_shape_factor=2.0, 12 enquiry_gap_factor=0.2, 13 13 number_of_barrels=1): 14 14 """Create polygons at the end of a culvert inlet and outlet. … … 23 23 Input (optional): 24 24 height - culvert height, defaults to width making a square culvert 25 enquiry_gap_factor - sets the distance to the enquiry polygon 26 enquiry_shape_factor - sets the shape of the enquiry polygon 27 (large value widens polygon but reduces height 28 to preserve same area as exchange polygon) 25 enquiry_gap_factor - sets the distance to the enquiry point as fraction of the height 29 26 number_of_barrels - number of identical pipes. 30 27 … … 34 31 'exchange_polygon0' - polygon defining the flow area at end_point0 35 32 'exchange_polygon1' - polygon defining the flow area at end_point1 36 'enquiry_polygon0' - polygon defining the enquiry field near end_point0 37 'enquiry_polygon1' - polygon defining the enquiry field near end_point1 33 'enquiry_point0' - point beyond exchange_polygon0 34 'enquiry_point1' - point beyond exchange_polygon1 35 'vector' 36 'length' 37 'normal' 38 38 """ 39 39 … … 62 62 63 63 # Unit direction vector and normal 64 vector /= length 65 normal = array([-dy, dx])/length 64 vector /= length # Unit vector in culvert direction 65 normal = array([-dy, dx])/length # Normal vector 66 66 67 culvert_polygons['vector'] = vector 68 culvert_polygons['length'] = length 69 culvert_polygons['normal'] = normal 67 70 68 71 # Short hands 69 72 w = 0.5*width*normal # Perpendicular vector of 1/2 width 70 h = height*vector # Vector of length=height in the 71 # direction of the culvert 73 h = height*vector # Vector of length=height in the 74 # direction of the culvert 75 gap = (1 + enquiry_gap_factor)*h 76 72 77 73 # Build exchange polygon 078 # Build exchange polygon and enquiry point for opening 0 74 79 p0 = end_point0 + w 75 80 p1 = end_point0 - w … … 77 82 p3 = p0 - h 78 83 culvert_polygons['exchange_polygon0'] = array([p0,p1,p2,p3]) 84 culvert_polygons['enquiry_point0'] = end_point0 - gap 85 79 86 80 # Build exchange polygon 187 # Build exchange polygon and enquiry point for opening 1 81 88 p0 = end_point1 + w 82 89 p1 = end_point1 - w … … 84 91 p3 = p0 + h 85 92 culvert_polygons['exchange_polygon1'] = array([p0,p1,p2,p3]) 93 culvert_polygons['enquiry_point1'] = end_point1 + gap 86 94 95 # Check that enquiry polygons are outside exchange polygons 96 for key1 in ['exchange_polygon0', 97 'exchange_polygon1']: 98 polygon = culvert_polygons[key1] 99 area = polygon_area(polygon) 100 101 msg = 'Polygon %s ' %(polygon) 102 msg += ' has area = %f' % area 103 assert area > 0.0, msg 87 104 88 89 # Redefine shorthands for enquiry polygons 90 w = w*enquiry_shape_factor 91 h = h/enquiry_shape_factor 92 gap = (enquiry_gap_factor + h)*vector 93 94 # Build enquiry polygon 0 95 p0 = end_point0 + w - gap 96 p1 = end_point0 - w - gap 97 p2 = p1 - h 98 p3 = p0 - h 99 culvert_polygons['enquiry_polygon0'] = array([p0,p1,p2,p3]) 100 101 # Build enquiry polygon 1 102 p0 = end_point1 + w + gap 103 p1 = end_point1 - w + gap 104 p2 = p1 + h 105 p3 = p0 + h 106 culvert_polygons['enquiry_polygon1'] = array([p0,p1,p2,p3]) 105 for key2 in ['enquiry_point0', 'enquiry_point1']: 106 point = culvert_polygons[key2] 107 msg = 'Enquiry point falls inside an enquiry point.' 108 msg += 'Email Ole.Nielsen@ga.gov.au' 109 assert not inside_polygon(point, polygon), msg 107 110 108 111 # Return results 109 culvert_polygons['vector'] = vector110 culvert_polygons['length'] = length111 culvert_polygons['normal'] = normal112 112 return culvert_polygons
Note: See TracChangeset
for help on using the changeset viewer.