Changeset 5301
 Timestamp:
 May 9, 2008, 4:09:56 PM (16 years ago)
 Location:
 anuga_core/source/anuga/culvert_flows
 Files:

 1 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/culvert_flows/culvert_polygons.py
r5298 r5301 3 3 4 4 # Import necessary modules 5 from math import atan 6 from math import sin 7 from math import cos 8 from math import pi 9 from math import degrees 5 from math import sqrt 6 from Numeric import array, sum 10 7 11 12 def create_culvert_polygons(center1, center2, 8 def create_culvert_polygons(end_point0, 9 end_point1, 13 10 width, height=None, 14 11 enquiry_gap_factor=1.0, … … 19 16 20 17 Input (mandatory): 21 center1  center point atone end of the culvert (x,y)22 center2  center point atother end of the culvert (x,y)18 end_point0  one end of the culvert (x,y) 19 end_point1  other end of the culvert (x,y) 23 20 width  culvert width 24 21 … … 27 24 enquiry_gap_factor  sets the distance to the enquiry polygon 28 25 enquiry_shape_factor  sets the shape of the enquiry polygon 26 (large value widens polygon but reduces height 27 to preserve same area as exchange polygon) 29 28 30 29 Output: 31 30 32 31 Dictionary of four polygons. The dictionary keys are: 33 'exchange_polygon1'  polygon defining the flow area at center1 34 'exchange_polygon2'  polygon defining the flow area at center2 35 'enquiry_polygon1'  polygon defining the enquiry field near center1 36 'enquiry_polygon2'  polygon defining the enquiry field near center2 37 32 'exchange_polygon0'  polygon defining the flow area at end_point0 33 'exchange_polygon1'  polygon defining the flow area at end_point1 34 'enquiry_polygon0'  polygon defining the enquiry field near end_point0 35 'enquiry_polygon1'  polygon defining the enquiry field near end_point1 38 36 """ 39 37 … … 43 41 height = width 44 42 43 # Dictionary for calculated polygons 44 culvert_polygons = {} 45 45 46 46 47 # Calculate geometry 47 x 1, y1 = center148 x 2, y2 = center248 x0, y0 = end_point0 49 x1, y1 = end_point1 49 50 50 dx =x2x151 dy =y2y151 dx = x1x0 52 dy = y1y0 52 53 54 vector = array([dx, dy]) 55 length = sqrt(sum(vector**2)) 53 56 54 # Calculate slope a in y=a*x+b to project points along the direction of the culvert 57 # Unit direction vector and normal 58 vector /= length 59 normal = array([dy, dx])/length 60 55 61 56 if dx==0.0: 57 # vertical 58 slope=0.0 59 else: 60 slope=dy/dx 62 # Short hands 63 w = width/2*normal # Perpendicular vector of 1/2 width 64 h = height*vector # Vector of length=height in the 65 # direction of the culvert 61 66 67 # Build exchange polygon 0 68 p0 = end_point0 + w 69 p1 = end_point0  w 70 p2 = p1  h 71 p3 = p0  h 72 culvert_polygons['exchange_polygon0'] = array([p0,p1,p2,p3]) 62 73 63 # Intercept b=Y1slope*X164 65 66 # Calculate basic trigonometric values67 74 68 # FIXME (Ole): Suggest using 69 # from anuga.utilities.numerical_tools import angle 70 # theta = angle(slope) 71 # sin_theta=sin(theta) 72 # cos_theta=cos(theta) 73 74 75 angle_atan=atan(slope) 76 angle_sin=sin(angle_atan) 77 angle_cos=cos(angle_atan) 78 79 # If you need to establish orientation USE THIS other Wise SKIP !!! 80 if dx == 0.0: 81 if dy > 0.0: # vertical up 82 angle_rad = 0.0 83 else: # dy < 0 vertical down 84 angle_rad = pi # 180 Degrees 85 elif dx<0: 86 if slope > 0: # Pt 2 less than Pt 1 87 angle_rad=3*pi/2angle_atan 88 else: 89 if slope > 0: 90 angle_rad=pi/2angle_atan 91 92 # Calculate derived perpendicular trigonometric values 93 angle_per=atan(angle_rad) 94 angle_per_sin=sin(angle_rad) 95 angle_per_cos=cos(angle_rad) 96 97 angle_deg = degrees(angle_rad) # Convert to degrees 75 # Build exchange polygon 1 76 p0 = end_point1 + w 77 p1 = end_point1  w 78 p2 = p1 + h 79 p3 = p0 + h 80 culvert_polygons['exchange_polygon1'] = array([p0,p1,p2,p3]) 98 81 99 82 100 83 101 # From the 2 X,Y coordinates at each end of the Culvert 102 # establish Polygons 103 #  to abstract flow and/or inject flow, 104 #  and set up enquiry area for total energy 84 # Redefine shorthands for enquiry polygons 85 w = w*enquiry_shape_factor 86 h = h/enquiry_shape_factor 87 gap = (enquiry_gap_factor + h)*vector 88 89 # Build enquiry polygon 0 90 p0 = end_point0 + w  gap 91 p1 = end_point0  w  gap 92 p2 = p1  h 93 p3 = p0  h 94 culvert_polygons['enquiry_polygon0'] = array([p0,p1,p2,p3]) 105 95 106 culvert_polygons = {} 96 # Build enquiry polygon 1 97 p0 = end_point1 + w + gap 98 p1 = end_point1  w + gap 99 p2 = p1 + h 100 p3 = p0 + h 101 culvert_polygons['enquiry_polygon1'] = array([p0,p1,p2,p3]) 107 102 108 109 # Exchange polygon 1 110 pt1x1=x1 111 pt1y1=y1 112 pt2x1=pt1x1angle_sin*width/2 113 pt2y1=pt1y1+angle_cos*width/2 114 pt3x1=pt2x1angle_per_sin*height 115 pt3y1=pt2y1angle_per_cos*height 116 pt4x1=pt3x1+angle_sin*width 117 pt4y1=pt3y1angle_cos*width 118 pt5x1=pt4x1+angle_per_sin*height 119 pt5y1=pt4y1+angle_per_cos*height 120 pt6x1=pt5x1angle_sin*width/2 121 pt6y1=pt5y1+angle_cos*width/2 122 123 culvert_polygons['exchange_polygon1'] = [[pt1x1,pt1y1], [pt2x1,pt2y1], 124 [pt3x1,pt3y1], [pt4x1,pt4y1], 125 [pt5x1,pt5y1]] 126 127 128 # Exchange polygon 2 129 pt1x2=x2 130 pt1y2=y2 131 pt2x2=pt1x2+angle_sin*width/2 132 pt2y2=pt1y2angle_cos*width/2 133 pt3x2=pt2x2+angle_per_sin*height 134 pt3y2=pt2y2+angle_per_cos*height 135 pt4x2=pt3x2angle_sin*width 136 pt4y2=pt3y2+angle_cos*width 137 pt5x2=pt4x2angle_per_sin*height 138 pt5y2=pt4y2angle_per_cos*height 139 pt6x2=pt5x2+angle_sin*width/2 140 pt6y2=pt5y2angle_cos*width/2 141 142 culvert_polygons['exchange_polygon2'] = [[pt1x2,pt1y2], [pt2x2,pt2y2], 143 [pt3x2,pt3y2], [pt4x2,pt4y2], 144 [pt5x2,pt5y2]] 145 146 147 # Calculate the the energy enquiry fields 148 # using a couple of distances first using factors from above 149 150 enq_dist=height*(1+enquiry_gap_factor) 151 enq_width=width*enquiry_shape_factor 152 enq_height=height/enquiry_shape_factor 153 154 155 156 # Enquiry polygon 1 157 pt1qx1=x1angle_per_sin*enq_dist 158 pt1qy1=y1angle_per_cos*enq_dist 159 pt2qx1=pt1qx1angle_sin*enq_width/2 160 pt2qy1=pt1qy1+angle_cos*enq_width/2 161 pt3qx1=pt2qx1angle_per_sin*enq_height 162 pt3qy1=pt2qy1angle_per_cos*enq_height 163 pt4qx1=pt3qx1+angle_sin*enq_width 164 pt4qy1=pt3qy1angle_cos*enq_width 165 pt5qx1=pt4qx1+angle_per_sin*enq_height 166 pt5qy1=pt4qy1+angle_per_cos*enq_height 167 pt6qx1=pt5qx1angle_sin*enq_width/2 168 pt6qy1=pt5qy1+angle_cos*enq_width/2 169 170 171 culvert_polygons['enquiry_polygon1'] = [[pt1qx1,pt1qy1], [pt2qx1,pt2qy1], 172 [pt3qx1,pt3qy1], [pt4qx1,pt4qy1], 173 [pt5qx1,pt5qy1]] 174 175 176 # Enquiry polygon 2 177 pt1qx2=x2+angle_per_sin*enq_dist 178 pt1qy2=y2+angle_per_cos*enq_dist 179 pt2qx2=pt1qx2+angle_sin*enq_width/2 180 pt2qy2=pt1qy2angle_cos*enq_width/2 181 pt3qx2=pt2qx2+angle_per_sin*enq_height 182 pt3qy2=pt2qy2+angle_per_cos*enq_height 183 pt4qx2=pt3qx2angle_sin*enq_width 184 pt4qy2=pt3qy2+angle_cos*enq_width 185 pt5qx2=pt4qx2angle_per_sin*enq_height 186 pt5qy2=pt4qy2angle_per_cos*enq_height 187 pt6qx2=pt5qx2+angle_sin*enq_width/2 188 pt6qy2=pt5qy2angle_cos*enq_width/2 189 190 191 culvert_polygons['enquiry_polygon2'] = [[pt1qx2,pt1qy2], [pt2qx2,pt2qy2], 192 [pt3qx2,pt3qy2], [pt4qx2,pt4qy2], 193 [pt5qx2,pt5qy2]] 194 195 196 197 # Return dictionary of all four polygons 103 # Return results 104 culvert_polygons['vector'] = vector 105 culvert_polygons['length'] = length 106 culvert_polygons['normal'] = normal 198 107 return culvert_polygons 
anuga_core/source/anuga/culvert_flows/culvert_polygons_example.py
r5299 r5301 8 8 """ 9 9 10 from anuga.utilities.polygon import read_polygon,plot_polygons10 from anuga.utilities.polygon import plot_polygons 11 11 from culvert_polygons import create_culvert_polygons 12 12 13 # Read these values or passed from some where14 culvert_description= '2.4mx1.2m Box Culvert'15 16 13 # Culvert location 17 x 1= 5.0; y1= 5.0 # One end18 x 2=15.0; y2=10.0 # Other end14 x0 = 5.0; y0 = 5.0 # One end 15 x1 = 15.0; y1 = 20.0 # Other end 19 16 20 17 # Culvert Size … … 24 21 25 22 # Call function 26 P = create_culvert_polygons( center1=[x1, y1],27 center2=[x2, y2],23 P = create_culvert_polygons(end_point0=[x0, y0], 24 end_point1=[x1, y1], 28 25 width=culvert_width, 29 26 height=culvert_height) … … 31 28 32 29 33 culv_polygon = [[x 1,y1], [x2,y2]]30 culv_polygon = [[x0,y0], [x1,y1]] 34 31 plot_polygons([culv_polygon, 32 P['exchange_polygon0'], 35 33 P['exchange_polygon1'], 36 P['exchange_polygon2'], 37 P['enquiry_polygon1'], 38 P['enquiry_polygon2']], 39 figname='Culvert Check') 34 P['enquiry_polygon0'], 35 P['enquiry_polygon1']], 36 figname='culvert_polygon_example') 40 37 41 38
Note: See TracChangeset
for help on using the changeset viewer.