# source:trunk/anuga_core/source/anuga/structures/culvert_polygons.py@7954

Last change on this file since 7954 was 7954, checked in by habili, 14 years ago

Minor refactoring

File size: 3.6 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,
14                            height=None,
15                            enquiry_gap_factor=0.2,
16                            number_of_barrels=1):
17    """Create polygons at the end of a culvert inlet and outlet.
18    At either end two polygons will be created; one for the actual flow to pass through and one a little further away
19    for enquiring the total energy at both ends of the culvert and transferring flow.
20
21    Input (mandatory):
22        end_point0 - one end of the culvert (x,y)
23        end_point1 - other end of the culvert (x,y)
24        width - culvert width
25
26    Input (optional):
27        height - culvert height, defaults to width making a square culvert
28        enquiry_gap_factor - sets the distance to the enquiry point as fraction of the height
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_point0' - point beyond exchange_polygon0
37            'enquiry_point1' - point beyond exchange_polygon1
38            'vector'
39            'length'
40            'normal'
41    """
42
43
44    # Input check
45    if height is None:
46        height = width
47
48    # Dictionary for calculated polygons
49    culvert_polygons = {}
50
51
52    # Calculate geometry
53    x0, y0 = end_point0
54    x1, y1 = end_point1
55
56    dx = x1-x0
57    dy = y1-y0
58
59    dxdy = num.array([dx, dy])
60    length = sqrt(num.sum(dxdy**2))
61
62    # Adjust polygon width to number of barrels in this culvert
63    width *= number_of_barrels
64
65
66    # Unit direction vector and normal
67    dxdy /= length                 # Unit vector in culvert direction
68    normal = num.array([-dy, dx])/length # Normal vector
69
70    culvert_polygons['vector'] = dxdy
71    culvert_polygons['length'] = length
72    culvert_polygons['normal'] = normal
73
74    # Short hands
75    w = 0.5*width*normal # Perpendicular vector of 1/2 width
76    h = height*dxdy # Vector of length=height in the
77                    # direction of the culvert
78    gap = (1 + enquiry_gap_factor)*h
79
80
81    # Build exchange polygon and enquiry point for opening 0
82    p0 = end_point0 + w
83    p1 = end_point0 - w
84    p2 = p1 - h
85    p3 = p0 - h
86    culvert_polygons['exchange_polygon0'] = num.array([p0,p1,p2,p3])
87    culvert_polygons['enquiry_point0'] = end_point0 - gap
88
89
90    # Build exchange polygon and enquiry point for opening 1
91    p0 = end_point1 + w
92    p1 = end_point1 - w
93    p2 = p1 + h
94    p3 = p0 + h
95    culvert_polygons['exchange_polygon1'] = num.array([p0,p1,p2,p3])
96    culvert_polygons['enquiry_point1'] = end_point1 + gap
97
98    # Check that enquiry polygons are outside exchange polygons
99    for key1 in ['exchange_polygon0', '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            assert not inside_polygon(point, polygon), msg
111
112    # Return results
113    return culvert_polygons
Note: See TracBrowser for help on using the repository browser.