[7993] | 1 | from anuga.geometry.polygon import inside_polygon, polygon_area |
---|
| 2 | from anuga.config import g |
---|
| 3 | import numpy as num |
---|
| 4 | import math |
---|
| 5 | import inlet |
---|
| 6 | |
---|
| 7 | class Structure_operator: |
---|
[7995] | 8 | """Structure Operator - transfer water from one rectangular box to another. |
---|
[7993] | 9 | Sets up the geometry of problem |
---|
| 10 | |
---|
| 11 | This is the base class for culverts. Inherit from this class (and overwrite |
---|
| 12 | compute_discharge method for specific subclasses) |
---|
| 13 | |
---|
| 14 | Input: Two points, pipe_size (either diameter or width, height), |
---|
| 15 | mannings_rougness, |
---|
| 16 | """ |
---|
| 17 | |
---|
| 18 | def __init__(self, |
---|
| 19 | domain, |
---|
| 20 | end_point0, |
---|
| 21 | end_point1, |
---|
| 22 | width, |
---|
| 23 | height, |
---|
| 24 | apron, |
---|
| 25 | manning, |
---|
| 26 | enquiry_gap, |
---|
[7996] | 27 | description, |
---|
[7993] | 28 | verbose): |
---|
| 29 | |
---|
| 30 | self.domain = domain |
---|
| 31 | self.domain.set_fractional_step_operator(self) |
---|
| 32 | self.end_points = [end_point0, end_point1] |
---|
| 33 | |
---|
| 34 | if height is None: |
---|
| 35 | height = width |
---|
| 36 | |
---|
| 37 | if apron is None: |
---|
| 38 | apron = width |
---|
| 39 | |
---|
| 40 | self.width = width |
---|
| 41 | self.height = height |
---|
| 42 | self.apron = apron |
---|
| 43 | self.manning = manning |
---|
| 44 | self.enquiry_gap = enquiry_gap |
---|
[7996] | 45 | self.description = description |
---|
[7993] | 46 | self.verbose = verbose |
---|
[7995] | 47 | |
---|
| 48 | self.discharge = 0.0 |
---|
| 49 | self.velocity = 0.0 |
---|
[7996] | 50 | self.delta_total_energy = 0.0 |
---|
| 51 | self.driving_energy = 0.0 |
---|
[7993] | 52 | |
---|
| 53 | self.__create_exchange_polygons() |
---|
| 54 | |
---|
| 55 | self.inlets = [] |
---|
| 56 | polygon0 = self.inlet_polygons[0] |
---|
| 57 | enquiry_point0 = self.inlet_equiry_points[0] |
---|
| 58 | outward_vector0 = self.culvert_vector |
---|
| 59 | self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_point0, outward_vector0)) |
---|
| 60 | |
---|
| 61 | polygon1 = self.inlet_polygons[1] |
---|
| 62 | exchange_polygon1 = self.inlet_equiry_points[1] |
---|
| 63 | outward_vector1 = - self.culvert_vector |
---|
| 64 | self.inlets.append(inlet.Inlet(self.domain, polygon1, exchange_polygon1, outward_vector1)) |
---|
| 65 | |
---|
| 66 | def __call__(self): |
---|
| 67 | |
---|
| 68 | pass |
---|
| 69 | |
---|
| 70 | def __create_exchange_polygons(self): |
---|
| 71 | |
---|
| 72 | """Create polygons at the end of a culvert inlet and outlet. |
---|
| 73 | At either end two polygons will be created; one for the actual flow to pass through and one a little further away |
---|
| 74 | for enquiring the total energy at both ends of the culvert and transferring flow. |
---|
| 75 | """ |
---|
| 76 | |
---|
| 77 | # Calculate geometry |
---|
| 78 | x0, y0 = self.end_points[0] |
---|
| 79 | x1, y1 = self.end_points[1] |
---|
| 80 | |
---|
| 81 | dx = x1 - x0 |
---|
| 82 | dy = y1 - y0 |
---|
| 83 | |
---|
| 84 | self.culvert_vector = num.array([dx, dy]) |
---|
| 85 | self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2)) |
---|
| 86 | assert self.culvert_length > 0.0, 'The length of culvert is less than 0' |
---|
| 87 | |
---|
| 88 | # Unit direction vector and normal |
---|
| 89 | self.culvert_vector /= self.culvert_length # Unit vector in culvert direction |
---|
| 90 | self.culvert_normal = num.array([-dy, dx])/self.culvert_length # Normal vector |
---|
| 91 | |
---|
| 92 | # Short hands |
---|
| 93 | w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width |
---|
| 94 | h = self.apron*self.culvert_vector # Vector of length=height in the |
---|
| 95 | # direction of the culvert |
---|
| 96 | |
---|
| 97 | gap = (1 + self.enquiry_gap)*h |
---|
| 98 | |
---|
| 99 | self.inlet_polygons = [] |
---|
| 100 | self.inlet_equiry_points = [] |
---|
| 101 | |
---|
| 102 | # Build exchange polygon and enquiry point |
---|
| 103 | for i in [0, 1]: |
---|
| 104 | i0 = (2*i-1) |
---|
| 105 | p0 = self.end_points[i] + w |
---|
| 106 | p1 = self.end_points[i] - w |
---|
| 107 | p2 = p1 + i0*h |
---|
| 108 | p3 = p0 + i0*h |
---|
| 109 | ep = self.end_points[i] + i0*gap |
---|
| 110 | |
---|
| 111 | self.inlet_polygons.append(num.array([p0, p1, p2, p3])) |
---|
| 112 | self.inlet_equiry_points.append(ep) |
---|
| 113 | |
---|
| 114 | # Check that enquiry points are outside inlet polygons |
---|
| 115 | for i in [0,1]: |
---|
| 116 | polygon = self.inlet_polygons[i] |
---|
| 117 | ep = self.inlet_equiry_points[i] |
---|
| 118 | |
---|
| 119 | area = polygon_area(polygon) |
---|
| 120 | |
---|
| 121 | msg = 'Polygon %s ' %(polygon) |
---|
| 122 | msg += ' has area = %f' % area |
---|
| 123 | assert area > 0.0, msg |
---|
| 124 | |
---|
| 125 | msg = 'Enquiry point falls inside an exchange polygon.' |
---|
| 126 | assert not inside_polygon(ep, polygon), msg |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | #print ' outflow volume ',outflow.get_total_water_volume() |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | def print_stats(self): |
---|
| 133 | |
---|
| 134 | print '=====================================' |
---|
| 135 | print 'Generic Culvert Operator' |
---|
| 136 | print '=====================================' |
---|
| 137 | |
---|
| 138 | print 'Culvert' |
---|
| 139 | print self.culvert |
---|
| 140 | |
---|
| 141 | print 'Culvert Routine' |
---|
| 142 | print self.routine |
---|
| 143 | |
---|
| 144 | for i, inlet in enumerate(self.inlets): |
---|
| 145 | print '-------------------------------------' |
---|
| 146 | print 'Inlet %i' % i |
---|
| 147 | print '-------------------------------------' |
---|
| 148 | |
---|
| 149 | print 'inlet triangle indices and centres' |
---|
| 150 | print inlet.triangle_indices |
---|
| 151 | print self.domain.get_centroid_coordinates()[inlet.triangle_indices] |
---|
| 152 | |
---|
| 153 | print 'polygon' |
---|
| 154 | print inlet.polygon |
---|
| 155 | |
---|
| 156 | print '=====================================' |
---|
| 157 | |
---|
[7995] | 158 | |
---|
| 159 | def structure_statistics(self): |
---|
| 160 | |
---|
| 161 | message = '---------------------------\n' |
---|
[7996] | 162 | message += 'Structure report for structure %s:\n' % self.description |
---|
[7995] | 163 | message += '--------------------------\n' |
---|
| 164 | message += 'Discharge [m^3/s]: %.2f\n' % self.discharge |
---|
| 165 | message += 'Velocity [m/s]: %.2f\n' % self.velocity |
---|
[7996] | 166 | message += 'Inlet Driving Energy %.2f\n' % self.driving_energy |
---|
| 167 | message += 'delta total energy %.2f\n' % self.delta_total_energy |
---|
[7995] | 168 | # message += 'Net boundary flow by tags [m^3/s]\n' |
---|
| 169 | # for tag in boundary_flows: |
---|
| 170 | # message += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag]) |
---|
| 171 | # |
---|
| 172 | # message += 'Total net boundary flow [m^3/s]: %.2f\n' % \ |
---|
| 173 | # (total_boundary_inflow + total_boundary_outflow) |
---|
| 174 | # message += 'Total volume in domain [m^3]: %.2f\n' % \ |
---|
| 175 | # self.compute_total_volume() |
---|
| 176 | # |
---|
| 177 | # # The go through explicit forcing update and record the rate of change |
---|
| 178 | # # for stage and |
---|
| 179 | # # record into forcing_inflow and forcing_outflow. Finally compute |
---|
| 180 | # # integral of depth to obtain total volume of domain. |
---|
| 181 | # |
---|
| 182 | # FIXME(Ole): This part is not yet done. |
---|
| 183 | |
---|
| 184 | return message |
---|
| 185 | |
---|
[7993] | 186 | def get_inlets(self): |
---|
| 187 | |
---|
| 188 | return self.inlets |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | def get_culvert_length(self): |
---|
| 192 | |
---|
| 193 | return self.culvert_length |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | def get_culvert_width(self): |
---|
| 197 | |
---|
| 198 | return self.width |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | def get_culvert_height(self): |
---|
| 202 | |
---|
| 203 | return self.height |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | def get_culvert_apron(self): |
---|
| 207 | |
---|
| 208 | return self.apron |
---|