1 | """Functions for geometries related to culvert flows |
---|
2 | """ |
---|
3 | |
---|
4 | # Import necessary modules |
---|
5 | from math import sqrt |
---|
6 | from anuga.geometry.polygon import inside_polygon, polygon_area |
---|
7 | |
---|
8 | import numpy as num |
---|
9 | |
---|
10 | |
---|
11 | def 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 |
---|