[7939] | 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, |
---|
[7954] | 13 | width, |
---|
| 14 | height=None, |
---|
[7939] | 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 | |
---|
[7954] | 59 | dxdy = num.array([dx, dy]) |
---|
| 60 | length = sqrt(num.sum(dxdy**2)) |
---|
[7939] | 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 |
---|
[7954] | 67 | dxdy /= length # Unit vector in culvert direction |
---|
[7939] | 68 | normal = num.array([-dy, dx])/length # Normal vector |
---|
| 69 | |
---|
[7954] | 70 | culvert_polygons['vector'] = dxdy |
---|
[7939] | 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 |
---|
[7954] | 76 | h = height*dxdy # Vector of length=height in the |
---|
| 77 | # direction of the culvert |
---|
[7939] | 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 |
---|
[7954] | 99 | for key1 in ['exchange_polygon0', 'exchange_polygon1']: |
---|
[7939] | 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 |
---|