source: anuga_core/source/anuga/culvert_flows/culvert_polygons.py @ 5585

Last change on this file since 5585 was 5585, checked in by ole, 14 years ago

Multibarrel culvert functionality.

File size: 3.5 KB
Line 
1"""Functions for geometries related to culvert flows
2"""
3
4# Import necessary modules
5from math import sqrt
6from Numeric import array, sum
7
8def create_culvert_polygons(end_point0,
9                            end_point1, 
10                            width, height=None,
11                            enquiry_gap_factor=1.0,
12                            enquiry_shape_factor=2.0,
13                            number_of_barrels=1):
14    """Create polygons at the end of a culvert inlet and outlet.
15    At either end two polygons will be created; one for the actual flow to pass through and one a little further away
16    for enquiring the total energy at both ends of the culvert and transferring flow.
17
18    Input (mandatory):
19        end_point0 - one end of the culvert (x,y)
20        end_point1 - other end of the culvert (x,y)       
21        width - culvert width
22
23    Input (optional):       
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)
29        number_of_barrels - number of identical pipes.
30       
31    Output:
32
33        Dictionary of four polygons. The dictionary keys are:
34            'exchange_polygon0' - polygon defining the flow area at end_point0
35            '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           
38    """   
39
40
41    # Input check
42    if height is None:
43        height = width
44
45    # Dictionary for calculated polygons
46    culvert_polygons = {}
47   
48
49    # Calculate geometry
50    x0, y0 = end_point0
51    x1, y1 = end_point1
52
53    dx = x1-x0
54    dy = y1-y0
55
56    vector = array([dx, dy])
57    length = sqrt(sum(vector**2))
58
59    # Adjust polygon width to number of barrels in this culvert
60    width *= number_of_barrels
61   
62   
63    # Unit direction vector and normal
64    vector /= length
65    normal = array([-dy, dx])/length
66   
67
68    # Short hands
69    w = width/2*normal # Perpendicular vector of 1/2 width
70    h = height*vector  # Vector of length=height in the
71                       # direction of the culvert
72
73    # Build exchange polygon 0
74    p0 = end_point0 + w
75    p1 = end_point0 - w
76    p2 = p1 - h
77    p3 = p0 - h
78    culvert_polygons['exchange_polygon0'] = array([p0,p1,p2,p3])
79
80   
81    # Build exchange polygon 1
82    p0 = end_point1 + w
83    p1 = end_point1 - w
84    p2 = p1 + h
85    p3 = p0 + h
86    culvert_polygons['exchange_polygon1'] = array([p0,p1,p2,p3])
87
88
89
90    # Redefine shorthands for enquiry polygons
91    w = w*enquiry_shape_factor
92    h = h/enquiry_shape_factor
93    gap = (enquiry_gap_factor + h)*vector
94   
95    # Build enquiry polygon 0
96    p0 = end_point0 + w - gap
97    p1 = end_point0 - w - gap
98    p2 = p1 - h
99    p3 = p0 - h
100    culvert_polygons['enquiry_polygon0'] = array([p0,p1,p2,p3])
101
102    # Build enquiry polygon 1
103    p0 = end_point1 + w + gap
104    p1 = end_point1 - w + gap
105    p2 = p1 + h
106    p3 = p0 + h
107    culvert_polygons['enquiry_polygon1'] = array([p0,p1,p2,p3])   
108
109    # Return results
110    culvert_polygons['vector'] = vector
111    culvert_polygons['length'] = length
112    culvert_polygons['normal'] = normal   
113    return culvert_polygons
Note: See TracBrowser for help on using the repository browser.