| 1 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon |
|---|
| 2 | from anuga.config import velocity_protection, g |
|---|
| 3 | import math |
|---|
| 4 | |
|---|
| 5 | import numpy as num |
|---|
| 6 | |
|---|
| 7 | class Inlet: |
|---|
| 8 | """Contains information associated with each inlet |
|---|
| 9 | """ |
|---|
| 10 | |
|---|
| 11 | def __init__(self, domain, polygon, outward_culvert_vector=None): |
|---|
| 12 | |
|---|
| 13 | self.domain = domain |
|---|
| 14 | self.domain_bounding_polygon = self.domain.get_boundary_polygon() |
|---|
| 15 | self.polygon = polygon |
|---|
| 16 | self.outward_culvert_vector = outward_culvert_vector |
|---|
| 17 | |
|---|
| 18 | # FIXME (SR) Using get_triangle_containing_point which needs to be sped up |
|---|
| 19 | |
|---|
| 20 | self.compute_triangle_indices() |
|---|
| 21 | self.compute_area() |
|---|
| 22 | |
|---|
| 23 | |
|---|
| 24 | def compute_triangle_indices(self): |
|---|
| 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 | self.inlet_polygon = self.polygon |
|---|
| 31 | |
|---|
| 32 | # Check that polygon lies within the mesh. |
|---|
| 33 | for point in self.inlet_polygon: |
|---|
| 34 | msg = 'Point %s in polygon for forcing term' % str(point) |
|---|
| 35 | msg += ' did not fall within the domain boundary.' |
|---|
| 36 | assert is_inside_polygon(point, bounding_polygon), msg |
|---|
| 37 | |
|---|
| 38 | |
|---|
| 39 | self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon, verbose=True) |
|---|
| 40 | |
|---|
| 41 | if len(self.triangle_indices) == 0: |
|---|
| 42 | region = 'Inlet polygon=%s' % (self.inlet_polygon) |
|---|
| 43 | msg = 'No triangles have been identified in ' |
|---|
| 44 | msg += 'specified region: %s' % region |
|---|
| 45 | raise Exception, msg |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | |
|---|
| 50 | def compute_area(self): |
|---|
| 51 | |
|---|
| 52 | # Compute inlet area as the sum of areas of triangles identified |
|---|
| 53 | # by polygon. Must be called after compute_inlet_triangle_indices(). |
|---|
| 54 | if len(self.triangle_indices) == 0: |
|---|
| 55 | region = 'Inlet polygon=%s' % (self.inlet_polygon) |
|---|
| 56 | msg = 'No triangles have been identified in ' |
|---|
| 57 | msg += 'specified region: %s' % region |
|---|
| 58 | raise Exception, msg |
|---|
| 59 | |
|---|
| 60 | self.area = 0.0 |
|---|
| 61 | for j in self.triangle_indices: |
|---|
| 62 | self.area += self.domain.areas[j] |
|---|
| 63 | |
|---|
| 64 | msg = 'Inlet exchange area has area = %f' % self.area |
|---|
| 65 | assert self.area > 0.0 |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | def get_area(self): |
|---|
| 69 | |
|---|
| 70 | return self.area |
|---|
| 71 | |
|---|
| 72 | |
|---|
| 73 | def get_areas(self): |
|---|
| 74 | |
|---|
| 75 | # Must be called after compute_inlet_triangle_indices(). |
|---|
| 76 | return self.domain.areas.take(self.triangle_indices) |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | def get_stages(self): |
|---|
| 80 | |
|---|
| 81 | return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices) |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | def get_average_stage(self): |
|---|
| 85 | |
|---|
| 86 | return num.sum(self.get_stages()*self.get_areas())/self.area |
|---|
| 87 | |
|---|
| 88 | def get_elevations(self): |
|---|
| 89 | |
|---|
| 90 | return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices) |
|---|
| 91 | |
|---|
| 92 | def get_average_elevation(self): |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | return num.sum(self.get_elevations()*self.get_areas())/self.area |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | def get_xmoms(self): |
|---|
| 99 | |
|---|
| 100 | return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices) |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | def get_average_xmom(self): |
|---|
| 104 | |
|---|
| 105 | return num.sum(self.get_xmoms()*self.get_areas())/self.area |
|---|
| 106 | |
|---|
| 107 | |
|---|
| 108 | def get_ymoms(self): |
|---|
| 109 | |
|---|
| 110 | return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices) |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | def get_average_ymom(self): |
|---|
| 114 | |
|---|
| 115 | return num.sum(self.get_ymoms()*self.get_areas())/self.area |
|---|
| 116 | |
|---|
| 117 | def get_heights(self): |
|---|
| 118 | |
|---|
| 119 | return self.get_stages() - self.get_elevations() |
|---|
| 120 | |
|---|
| 121 | |
|---|
| 122 | def get_total_water_volume(self): |
|---|
| 123 | |
|---|
| 124 | return num.sum(self.get_heights()*self.get_areas()) |
|---|
| 125 | |
|---|
| 126 | |
|---|
| 127 | def get_average_height(self): |
|---|
| 128 | |
|---|
| 129 | return self.get_total_water_volume()/self.area |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | def get_velocities(self): |
|---|
| 133 | |
|---|
| 134 | heights = self.get_heights() |
|---|
| 135 | u = self.get_xmoms()/(heights + velocity_protection/heights) |
|---|
| 136 | v = self.get_ymoms()/(heights + velocity_protection/heights) |
|---|
| 137 | |
|---|
| 138 | return u, v |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | def get_xvelocities(self): |
|---|
| 142 | |
|---|
| 143 | heights = self.get_heights() |
|---|
| 144 | return self.get_xmoms()/(heights + velocity_protection/heights) |
|---|
| 145 | |
|---|
| 146 | def get_yvelocities(self): |
|---|
| 147 | |
|---|
| 148 | heights = self.get_heights() |
|---|
| 149 | return self.get_ymoms()/(heights + velocity_protection/heights) |
|---|
| 150 | |
|---|
| 151 | |
|---|
| 152 | def get_average_speed(self): |
|---|
| 153 | |
|---|
| 154 | u, v = self.get_velocities() |
|---|
| 155 | |
|---|
| 156 | average_u = num.sum(u*self.get_areas())/self.area |
|---|
| 157 | average_v = num.sum(v*self.get_areas())/self.area |
|---|
| 158 | |
|---|
| 159 | return math.sqrt(average_u**2 + average_v**2) |
|---|
| 160 | |
|---|
| 161 | |
|---|
| 162 | def get_average_velocity_head(self): |
|---|
| 163 | |
|---|
| 164 | return 0.5*self.get_average_speed()**2/g |
|---|
| 165 | |
|---|
| 166 | |
|---|
| 167 | def get_average_total_energy(self): |
|---|
| 168 | |
|---|
| 169 | return self.get_average_velocity_head() + self.get_average_stage() |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | def get_average_specific_energy(self): |
|---|
| 173 | |
|---|
| 174 | return self.get_average_velocity_head() + self.get_average_height() |
|---|
| 175 | |
|---|
| 176 | |
|---|
| 177 | def set_heights(self,height): |
|---|
| 178 | |
|---|
| 179 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height) |
|---|
| 180 | |
|---|
| 181 | |
|---|
| 182 | def set_stages(self,stage): |
|---|
| 183 | |
|---|
| 184 | self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage) |
|---|
| 185 | |
|---|
| 186 | |
|---|
| 187 | def set_xmoms(self,xmom): |
|---|
| 188 | |
|---|
| 189 | self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom) |
|---|
| 190 | |
|---|
| 191 | |
|---|
| 192 | def set_ymoms(self,ymom): |
|---|
| 193 | |
|---|
| 194 | self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom) |
|---|
| 195 | |
|---|
| 196 | |
|---|
| 197 | def set_elevations(self,elevation): |
|---|
| 198 | |
|---|
| 199 | self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation) |
|---|
