Changeset 5301


Ignore:
Timestamp:
May 9, 2008, 4:09:56 PM (17 years ago)
Author:
ole
Message:

Refactored generation of culvert polygons to a robust calculation.

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  
    33
    44# 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
     5from math import sqrt
     6from Numeric import array, sum
    107
    11 
    12 def create_culvert_polygons(center1, center2,
     8def create_culvert_polygons(end_point0,
     9                            end_point1,
    1310                            width, height=None,
    1411                            enquiry_gap_factor=1.0,
    … …  
    1916
    2017    Input (mandatory):
    21         center1 - center point at one end of the culvert (x,y)
    22         center2 - center point at other 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)       
    2320        width - culvert width
    2421
    … …  
    2724        enquiry_gap_factor - sets the distance to the enquiry polygon
    2825        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)
    2928       
    3029    Output:
    3130
    3231        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           
    3836    """   
    3937
    … …  
    4341        height = width
    4442
     43    # Dictionary for calculated polygons
     44    culvert_polygons = {}
     45   
    4546
    4647    # Calculate geometry
    47     x1, y1 = center1
    48     x2, y2 = center2
     48    x0, y0 = end_point0
     49    x1, y1 = end_point1
    4950
    50     dx=x2-x1
    51     dy=y2-y1
     51    dx = x1-x0
     52    dy = y1-y0
    5253
     54    vector = array([dx, dy])
     55    length = sqrt(sum(vector**2))
    5356
    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   
    5561
    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
    6166
     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])
    6273
    63     #  Intercept    b=Y1-slope*X1
    64 
    65 
    66     # Calculate basic trigonometric values
    6774   
    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/2-angle_atan
    88     else:
    89         if slope > 0:
    90             angle_rad=pi/2-angle_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])
    9881
    9982
    10083
    101     # From the 2 X,Y co-ordinates 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])
    10595
    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])   
    107102
    108 
    109     # Exchange polygon 1
    110     pt1x1=x1
    111     pt1y1=y1
    112     pt2x1=pt1x1-angle_sin*width/2
    113     pt2y1=pt1y1+angle_cos*width/2
    114     pt3x1=pt2x1-angle_per_sin*height
    115     pt3y1=pt2y1-angle_per_cos*height
    116     pt4x1=pt3x1+angle_sin*width
    117     pt4y1=pt3y1-angle_cos*width
    118     pt5x1=pt4x1+angle_per_sin*height
    119     pt5y1=pt4y1+angle_per_cos*height
    120     pt6x1=pt5x1-angle_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=pt1y2-angle_cos*width/2
    133     pt3x2=pt2x2+angle_per_sin*height
    134     pt3y2=pt2y2+angle_per_cos*height
    135     pt4x2=pt3x2-angle_sin*width
    136     pt4y2=pt3y2+angle_cos*width
    137     pt5x2=pt4x2-angle_per_sin*height
    138     pt5y2=pt4y2-angle_per_cos*height
    139     pt6x2=pt5x2+angle_sin*width/2
    140     pt6y2=pt5y2-angle_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=x1-angle_per_sin*enq_dist
    158     pt1qy1=y1-angle_per_cos*enq_dist
    159     pt2qx1=pt1qx1-angle_sin*enq_width/2
    160     pt2qy1=pt1qy1+angle_cos*enq_width/2
    161     pt3qx1=pt2qx1-angle_per_sin*enq_height
    162     pt3qy1=pt2qy1-angle_per_cos*enq_height
    163     pt4qx1=pt3qx1+angle_sin*enq_width
    164     pt4qy1=pt3qy1-angle_cos*enq_width
    165     pt5qx1=pt4qx1+angle_per_sin*enq_height
    166     pt5qy1=pt4qy1+angle_per_cos*enq_height
    167     pt6qx1=pt5qx1-angle_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=pt1qy2-angle_cos*enq_width/2
    181     pt3qx2=pt2qx2+angle_per_sin*enq_height
    182     pt3qy2=pt2qy2+angle_per_cos*enq_height
    183     pt4qx2=pt3qx2-angle_sin*enq_width
    184     pt4qy2=pt3qy2+angle_cos*enq_width
    185     pt5qx2=pt4qx2-angle_per_sin*enq_height
    186     pt5qy2=pt4qy2-angle_per_cos*enq_height
    187     pt6qx2=pt5qx2+angle_sin*enq_width/2
    188     pt6qy2=pt5qy2-angle_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   
    198107    return culvert_polygons

  • anuga_core/source/anuga/culvert_flows/culvert_polygons_example.py

    r5299 r5301  
    88"""
    99
    10 from anuga.utilities.polygon import read_polygon, plot_polygons
     10from anuga.utilities.polygon import plot_polygons
    1111from culvert_polygons import create_culvert_polygons
    1212
    13 # Read these values or passed from some where
    14 culvert_description= '2.4mx1.2m Box Culvert'
    15 
    1613# Culvert location
    17 x1= 5.0; y1= 5.0  # One end
    18 x2=15.0; y2=10.0  # Other end
     14x0 = 5.0; y0 = 5.0  # One end
     15x1 = 15.0; y1 = 20.0  # Other end
    1916
    2017# Culvert Size
    … …  
    2421
    2522# Call function
    26 P = create_culvert_polygons(center1=[x1, y1],
    27                             center2=[x2, y2],
     23P = create_culvert_polygons(end_point0=[x0, y0],
     24                            end_point1=[x1, y1],
    2825                            width=culvert_width,
    2926                            height=culvert_height)
    … …  
    3128
    3229
    33 culv_polygon = [[x1,y1], [x2,y2]]
     30culv_polygon = [[x0,y0], [x1,y1]]
    3431plot_polygons([culv_polygon,
     32               P['exchange_polygon0'],
    3533               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')
    4037
    4138
Note: See TracChangeset for help on using the changeset viewer.