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

Last change on this file since 5301 was 5301, checked in by ole, 15 years ago

Refactored generation of culvert polygons to a robust calculation.

File size: 3.3 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    """Create polygons at the end of a culvert inlet and outlet.
14    At either end two polygons will be created; one for the actual flow to pass through and one a little further away
15    for enquiring the total energy at both ends of the culvert and transferring flow.
16
17    Input (mandatory):
18        end_point0 - one end of the culvert (x,y)
19        end_point1 - other end of the culvert (x,y)       
20        width - culvert width
21
22    Input (optional):       
23        height - culvert height, defaults to width making a square culvert
24        enquiry_gap_factor - sets the distance to the enquiry polygon
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)
28       
29    Output:
30
31        Dictionary of four polygons. The dictionary keys are:
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           
36    """   
37
38
39    # Input check
40    if height is None:
41        height = width
42
43    # Dictionary for calculated polygons
44    culvert_polygons = {}
45   
46
47    # Calculate geometry
48    x0, y0 = end_point0
49    x1, y1 = end_point1
50
51    dx = x1-x0
52    dy = y1-y0
53
54    vector = array([dx, dy])
55    length = sqrt(sum(vector**2))
56
57    # Unit direction vector and normal
58    vector /= length
59    normal = array([-dy, dx])/length
60   
61
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
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])
73
74   
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])
81
82
83
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])
95
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])   
102
103    # Return results
104    culvert_polygons['vector'] = vector
105    culvert_polygons['length'] = length
106    culvert_polygons['normal'] = normal   
107    return culvert_polygons
Note: See TracBrowser for help on using the repository browser.