source: trunk/anuga_core/source/anuga/structures/culvert_polygons.py @ 7939

Last change on this file since 7939 was 7939, checked in by steve, 14 years ago

Adding in culvert routines into structures directory

File size: 3.7 KB
Line 
1"""Functions for geometries related to culvert flows
2"""
3
4# Import necessary modules
5from math import sqrt
6from anuga.geometry.polygon import inside_polygon, polygon_area
7
8import numpy as num
9
10
11def create_culvert_polygons(end_point0,
12                            end_point1, 
13                            width, height=None,
14                            enquiry_gap_factor=0.2,
15                            number_of_barrels=1):
16    """Create polygons at the end of a culvert inlet and outlet.
17    At either end two polygons will be created; one for the actual flow to pass through and one a little further away
18    for enquiring the total energy at both ends of the culvert and transferring flow.
19
20    Input (mandatory):
21        end_point0 - one end of the culvert (x,y)
22        end_point1 - other end of the culvert (x,y)       
23        width - culvert width
24
25    Input (optional):       
26        height - culvert height, defaults to width making a square culvert
27        enquiry_gap_factor - sets the distance to the enquiry point as fraction of the height
28        number_of_barrels - number of identical pipes.
29       
30    Output:
31
32        Dictionary of four polygons. The dictionary keys are:
33            'exchange_polygon0' - polygon defining the flow area at end_point0
34            'exchange_polygon1' - polygon defining the flow area at end_point1
35            'enquiry_point0' - point beyond exchange_polygon0
36            'enquiry_point1' - point beyond exchange_polygon1           
37            'vector'
38            'length'
39            'normal'
40    """   
41
42
43    # Input check
44    if height is None:
45        height = width
46
47    # Dictionary for calculated polygons
48    culvert_polygons = {}
49   
50
51    # Calculate geometry
52    x0, y0 = end_point0
53    x1, y1 = end_point1
54
55    dx = x1-x0
56    dy = y1-y0
57
58    vector = num.array([dx, dy])
59    length = sqrt(num.sum(vector**2))
60
61    # Adjust polygon width to number of barrels in this culvert
62    width *= number_of_barrels
63   
64   
65    # Unit direction vector and normal
66    vector /= length                 # Unit vector in culvert direction
67    normal = num.array([-dy, dx])/length # Normal vector
68   
69    culvert_polygons['vector'] = vector
70    culvert_polygons['length'] = length
71    culvert_polygons['normal'] = normal   
72
73    # Short hands
74    w = 0.5*width*normal # Perpendicular vector of 1/2 width
75    h = height*vector    # Vector of length=height in the
76                         # direction of the culvert
77    gap = (1 + enquiry_gap_factor)*h
78                         
79
80    # Build exchange polygon and enquiry point for opening 0
81    p0 = end_point0 + w
82    p1 = end_point0 - w
83    p2 = p1 - h
84    p3 = p0 - h
85    culvert_polygons['exchange_polygon0'] = num.array([p0,p1,p2,p3])
86    culvert_polygons['enquiry_point0'] = end_point0 - gap
87   
88
89    # Build exchange polygon and enquiry point for opening 1
90    p0 = end_point1 + w
91    p1 = end_point1 - w
92    p2 = p1 + h
93    p3 = p0 + h
94    culvert_polygons['exchange_polygon1'] = num.array([p0,p1,p2,p3])
95    culvert_polygons['enquiry_point1'] = end_point1 + gap 
96
97    # Check that enquiry polygons are outside exchange polygons
98    for key1 in ['exchange_polygon0',
99                 'exchange_polygon1']:
100        polygon = culvert_polygons[key1]
101        area = polygon_area(polygon)
102       
103        msg = 'Polygon %s ' %(polygon)
104        msg += ' has area = %f' % area
105        assert area > 0.0, msg
106
107        for key2 in ['enquiry_point0', 'enquiry_point1']:
108            point = culvert_polygons[key2]
109            msg = 'Enquiry point falls inside an enquiry point.'
110            msg += 'Email Ole.Nielsen@ga.gov.au'
111            assert not inside_polygon(point, polygon), msg
112
113    # Return results
114    return culvert_polygons
Note: See TracBrowser for help on using the repository browser.