[7972] | 1 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon |
---|
[7975] | 2 | from anuga.config import velocity_protection, g |
---|
[7973] | 3 | import math |
---|
[7972] | 4 | |
---|
| 5 | import numpy as num |
---|
| 6 | |
---|
| 7 | class Inlet: |
---|
| 8 | """Contains information associated with each inlet |
---|
| 9 | """ |
---|
| 10 | |
---|
[7984] | 11 | def __init__(self, domain, polygon, enquiry_pt, outward_culvert_vector=None, verbose=False): |
---|
[7972] | 12 | |
---|
| 13 | self.domain = domain |
---|
| 14 | self.domain_bounding_polygon = self.domain.get_boundary_polygon() |
---|
| 15 | self.polygon = polygon |
---|
[7984] | 16 | self.enquiry_pt = enquiry_pt |
---|
[7980] | 17 | self.outward_culvert_vector = outward_culvert_vector |
---|
[7984] | 18 | self.verbose = verbose |
---|
[7972] | 19 | |
---|
[7984] | 20 | self.compute_indices() |
---|
[7972] | 21 | self.compute_area() |
---|
| 22 | |
---|
| 23 | |
---|
[7984] | 24 | def compute_indices(self): |
---|
[7972] | 25 | |
---|
| 26 | # Get boundary (in absolute coordinates) |
---|
| 27 | bounding_polygon = self.domain_bounding_polygon |
---|
| 28 | domain_centroids = self.domain.get_centroid_coordinates(absolute=True) |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | # Check that polygon lies within the mesh. |
---|
[8003] | 32 | for point in self.polygon: |
---|
[7984] | 33 | msg = 'Point %s ' % str(point) |
---|
[7972] | 34 | msg += ' did not fall within the domain boundary.' |
---|
| 35 | assert is_inside_polygon(point, bounding_polygon), msg |
---|
[8003] | 36 | |
---|
| 37 | point = self.enquiry_pt |
---|
| 38 | msg = 'Enquiry Point %s ' % str(point) |
---|
| 39 | msg += ' did not fall within the domain boundary.' |
---|
| 40 | assert is_inside_polygon(point, bounding_polygon), msg |
---|
[7972] | 41 | |
---|
| 42 | |
---|
[7984] | 43 | self.triangle_indices = inside_polygon(domain_centroids, self.polygon, verbose=self.verbose) |
---|
[7980] | 44 | |
---|
[7972] | 45 | if len(self.triangle_indices) == 0: |
---|
[7984] | 46 | region = 'Inlet polygon=%s' % (self.polygon) |
---|
| 47 | msg = 'No triangles have been identified in region ' |
---|
[7972] | 48 | raise Exception, msg |
---|
| 49 | |
---|
[7984] | 50 | self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt) |
---|
[7972] | 51 | |
---|
[7980] | 52 | |
---|
| 53 | |
---|
[7972] | 54 | def compute_area(self): |
---|
| 55 | |
---|
| 56 | # Compute inlet area as the sum of areas of triangles identified |
---|
| 57 | # by polygon. Must be called after compute_inlet_triangle_indices(). |
---|
| 58 | if len(self.triangle_indices) == 0: |
---|
| 59 | region = 'Inlet polygon=%s' % (self.inlet_polygon) |
---|
[7984] | 60 | msg = 'No triangles have been identified in region ' |
---|
[7972] | 61 | raise Exception, msg |
---|
| 62 | |
---|
| 63 | self.area = 0.0 |
---|
| 64 | for j in self.triangle_indices: |
---|
| 65 | self.area += self.domain.areas[j] |
---|
| 66 | |
---|
| 67 | msg = 'Inlet exchange area has area = %f' % self.area |
---|
| 68 | assert self.area > 0.0 |
---|
| 69 | |
---|
| 70 | |
---|
[7980] | 71 | def get_area(self): |
---|
| 72 | |
---|
| 73 | return self.area |
---|
| 74 | |
---|
| 75 | |
---|
[7972] | 76 | def get_areas(self): |
---|
| 77 | |
---|
| 78 | # Must be called after compute_inlet_triangle_indices(). |
---|
| 79 | return self.domain.areas.take(self.triangle_indices) |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | def get_stages(self): |
---|
| 83 | |
---|
| 84 | return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices) |
---|
| 85 | |
---|
[7973] | 86 | |
---|
[7972] | 87 | def get_average_stage(self): |
---|
[7980] | 88 | |
---|
| 89 | return num.sum(self.get_stages()*self.get_areas())/self.area |
---|
[7972] | 90 | |
---|
| 91 | def get_elevations(self): |
---|
| 92 | |
---|
| 93 | return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices) |
---|
| 94 | |
---|
| 95 | def get_average_elevation(self): |
---|
| 96 | |
---|
[7980] | 97 | |
---|
| 98 | return num.sum(self.get_elevations()*self.get_areas())/self.area |
---|
[7972] | 99 | |
---|
| 100 | |
---|
| 101 | def get_xmoms(self): |
---|
| 102 | |
---|
| 103 | return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices) |
---|
| 104 | |
---|
[7973] | 105 | |
---|
[7972] | 106 | def get_average_xmom(self): |
---|
[7980] | 107 | |
---|
| 108 | return num.sum(self.get_xmoms()*self.get_areas())/self.area |
---|
[7972] | 109 | |
---|
| 110 | |
---|
| 111 | def get_ymoms(self): |
---|
| 112 | |
---|
| 113 | return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices) |
---|
| 114 | |
---|
| 115 | |
---|
| 116 | def get_average_ymom(self): |
---|
| 117 | |
---|
[7980] | 118 | return num.sum(self.get_ymoms()*self.get_areas())/self.area |
---|
[7972] | 119 | |
---|
[7984] | 120 | |
---|
[7972] | 121 | def get_heights(self): |
---|
| 122 | |
---|
| 123 | return self.get_stages() - self.get_elevations() |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | def get_total_water_volume(self): |
---|
| 127 | |
---|
[7975] | 128 | return num.sum(self.get_heights()*self.get_areas()) |
---|
[7980] | 129 | |
---|
| 130 | |
---|
[7972] | 131 | def get_average_height(self): |
---|
| 132 | |
---|
| 133 | return self.get_total_water_volume()/self.area |
---|
| 134 | |
---|
| 135 | |
---|
| 136 | def get_velocities(self): |
---|
| 137 | |
---|
[7980] | 138 | heights = self.get_heights() |
---|
| 139 | u = self.get_xmoms()/(heights + velocity_protection/heights) |
---|
| 140 | v = self.get_ymoms()/(heights + velocity_protection/heights) |
---|
[7972] | 141 | |
---|
| 142 | return u, v |
---|
[7980] | 143 | |
---|
| 144 | |
---|
| 145 | def get_xvelocities(self): |
---|
| 146 | |
---|
| 147 | heights = self.get_heights() |
---|
| 148 | return self.get_xmoms()/(heights + velocity_protection/heights) |
---|
| 149 | |
---|
| 150 | def get_yvelocities(self): |
---|
| 151 | |
---|
| 152 | heights = self.get_heights() |
---|
| 153 | return self.get_ymoms()/(heights + velocity_protection/heights) |
---|
[7972] | 154 | |
---|
| 155 | |
---|
[7975] | 156 | def get_average_speed(self): |
---|
[7972] | 157 | |
---|
[7973] | 158 | u, v = self.get_velocities() |
---|
[7972] | 159 | |
---|
[7980] | 160 | average_u = num.sum(u*self.get_areas())/self.area |
---|
| 161 | average_v = num.sum(v*self.get_areas())/self.area |
---|
[7973] | 162 | |
---|
[7975] | 163 | return math.sqrt(average_u**2 + average_v**2) |
---|
[7973] | 164 | |
---|
| 165 | |
---|
[7975] | 166 | def get_average_velocity_head(self): |
---|
| 167 | |
---|
| 168 | return 0.5*self.get_average_speed()**2/g |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | def get_average_total_energy(self): |
---|
[7973] | 172 | |
---|
[7975] | 173 | return self.get_average_velocity_head() + self.get_average_stage() |
---|
[7973] | 174 | |
---|
[7975] | 175 | |
---|
| 176 | def get_average_specific_energy(self): |
---|
[7973] | 177 | |
---|
[7975] | 178 | return self.get_average_velocity_head() + self.get_average_height() |
---|
[7973] | 179 | |
---|
[7975] | 180 | |
---|
[7984] | 181 | def get_enquiry_stage(self): |
---|
| 182 | |
---|
| 183 | return self.domain.quantities['stage'].centroid_values[self.enquiry_index] |
---|
| 184 | |
---|
| 185 | |
---|
| 186 | def get_enquiry_xmom(self): |
---|
| 187 | |
---|
| 188 | return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index] |
---|
| 189 | |
---|
| 190 | def get_enquiry_ymom(self): |
---|
| 191 | |
---|
| 192 | return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index] |
---|
| 193 | |
---|
| 194 | |
---|
| 195 | def get_enquiry_elevation(self): |
---|
| 196 | |
---|
| 197 | return self.domain.quantities['elevation'].centroid_values[self.enquiry_index] |
---|
| 198 | |
---|
| 199 | def get_enquiry_height(self): |
---|
| 200 | |
---|
| 201 | return self.get_enquiry_stage() - self.get_enquiry_elevation() |
---|
| 202 | |
---|
| 203 | |
---|
| 204 | def get_enquiry_velocity(self): |
---|
| 205 | |
---|
| 206 | height = self.get_enquiry_height() |
---|
| 207 | u = self.get_enquiry_xmom()/(height + velocity_protection/height) |
---|
| 208 | v = self.get_enquiry_ymom()/(height + velocity_protection/height) |
---|
| 209 | |
---|
| 210 | return u, v |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | def get_enquiry_xvelocity(self): |
---|
| 214 | |
---|
| 215 | height = self.get_enquiry_height() |
---|
| 216 | return self.get_enquiry_xmom()/(height + velocity_protection/height) |
---|
| 217 | |
---|
| 218 | def get_enquiry_yvelocity(self): |
---|
| 219 | |
---|
| 220 | height = self.get_enquiry_height() |
---|
| 221 | return self.get_enquiry_ymom()/(height + velocity_protection/height) |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | def get_enquiry_speed(self): |
---|
| 225 | |
---|
| 226 | u, v = self.get_enquiry_velocity() |
---|
| 227 | |
---|
| 228 | return math.sqrt(u**2 + v**2) |
---|
| 229 | |
---|
| 230 | |
---|
| 231 | def get_enquiry_velocity_head(self): |
---|
| 232 | |
---|
| 233 | return 0.5*self.get_enquiry_speed()**2/g |
---|
| 234 | |
---|
| 235 | |
---|
| 236 | def get_enquiry_total_energy(self): |
---|
| 237 | |
---|
| 238 | return self.get_enquiry_velocity_head() + self.get_enquiry_stage() |
---|
| 239 | |
---|
| 240 | |
---|
| 241 | def get_enquiry_specific_energy(self): |
---|
| 242 | |
---|
| 243 | return self.get_enquiry_velocity_head() + self.get_enquiry_height() |
---|
| 244 | |
---|
| 245 | |
---|
[7975] | 246 | def set_heights(self,height): |
---|
| 247 | |
---|
| 248 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height) |
---|
| 249 | |
---|
| 250 | |
---|
| 251 | def set_stages(self,stage): |
---|
| 252 | |
---|
| 253 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage) |
---|
| 254 | |
---|
[7977] | 255 | |
---|
[7975] | 256 | def set_xmoms(self,xmom): |
---|
| 257 | |
---|
[7984] | 258 | self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom) |
---|
[7975] | 259 | |
---|
| 260 | |
---|
| 261 | def set_ymoms(self,ymom): |
---|
| 262 | |
---|
[7984] | 263 | self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom) |
---|
[7975] | 264 | |
---|
| 265 | |
---|
| 266 | def set_elevations(self,elevation): |
---|
| 267 | |
---|
[7976] | 268 | self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation) |
---|