| 1 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect |
|---|
| 2 | from anuga.config import velocity_protection, g |
|---|
| 3 | import math |
|---|
| 4 | |
|---|
| 5 | import numpy as num |
|---|
| 6 | |
|---|
| 7 | import inlet |
|---|
| 8 | |
|---|
| 9 | class Inlet_enquiry(inlet.Inlet): |
|---|
| 10 | """Contains information associated with each inlet plus an enquiry point |
|---|
| 11 | """ |
|---|
| 12 | |
|---|
| 13 | def __init__(self, |
|---|
| 14 | domain, |
|---|
| 15 | poly, |
|---|
| 16 | enquiry_pt, |
|---|
| 17 | invert_elevation = None, |
|---|
| 18 | outward_culvert_vector=None, |
|---|
| 19 | verbose=False): |
|---|
| 20 | """ |
|---|
| 21 | poly can be a line or a polygon |
|---|
| 22 | """ |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | inlet.Inlet.__init__(self, domain, poly, verbose) |
|---|
| 26 | |
|---|
| 27 | self.enquiry_pt = enquiry_pt |
|---|
| 28 | self.invert_elevation = invert_elevation |
|---|
| 29 | self.outward_culvert_vector = outward_culvert_vector |
|---|
| 30 | |
|---|
| 31 | self.compute_enquiry_index() |
|---|
| 32 | |
|---|
| 33 | |
|---|
| 34 | def compute_enquiry_index(self): |
|---|
| 35 | |
|---|
| 36 | # Get boundary (in absolute coordinates) |
|---|
| 37 | bounding_polygon = self.domain_bounding_polygon |
|---|
| 38 | #domain_centroids = self.domain.get_centroid_coordinates(absolute=True) |
|---|
| 39 | #vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True) |
|---|
| 40 | |
|---|
| 41 | |
|---|
| 42 | point = self.enquiry_pt |
|---|
| 43 | msg = 'Enquiry Point %s ' % str(point) |
|---|
| 44 | msg += ' did not fall within the domain boundary.' |
|---|
| 45 | assert is_inside_polygon(point, bounding_polygon), msg |
|---|
| 46 | |
|---|
| 47 | try: |
|---|
| 48 | self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt) |
|---|
| 49 | except: |
|---|
| 50 | msg = "Enquiry point %s doesn't intersect mesh, maybe inside a building, try reducing enquiry_gap" % str(self.enquiry_pt) |
|---|
| 51 | raise Exception(msg) |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | if self.enquiry_index in self.triangle_indices: |
|---|
| 55 | msg = 'Enquiry point %s' % (self.enquiry_pt) |
|---|
| 56 | msg += 'is in an inlet triangle' |
|---|
| 57 | import warnings |
|---|
| 58 | warnings.warn(msg) |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | def get_enquiry_position(self): |
|---|
| 62 | |
|---|
| 63 | return self.domain.get_centroid_coordinates(absolute=True)[self.enquiry_index] |
|---|
| 64 | |
|---|
| 65 | def get_enquiry_stage(self): |
|---|
| 66 | |
|---|
| 67 | return self.domain.quantities['stage'].centroid_values[self.enquiry_index] |
|---|
| 68 | |
|---|
| 69 | |
|---|
| 70 | def get_enquiry_xmom(self): |
|---|
| 71 | |
|---|
| 72 | return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index] |
|---|
| 73 | |
|---|
| 74 | def get_enquiry_ymom(self): |
|---|
| 75 | |
|---|
| 76 | return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index] |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | def get_enquiry_elevation(self): |
|---|
| 80 | |
|---|
| 81 | return self.domain.quantities['elevation'].centroid_values[self.enquiry_index] |
|---|
| 82 | |
|---|
| 83 | def get_enquiry_depth(self): |
|---|
| 84 | |
|---|
| 85 | return max(self.get_enquiry_stage() - self.get_enquiry_invert_elevation(), 0.0) |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | def get_enquiry_water_depth(self): |
|---|
| 89 | |
|---|
| 90 | return self.get_enquiry_stage() - self.get_enquiry_elevation() |
|---|
| 91 | |
|---|
| 92 | |
|---|
| 93 | def get_enquiry_invert_elevation(self): |
|---|
| 94 | |
|---|
| 95 | if self.invert_elevation == None: |
|---|
| 96 | return self.get_enquiry_elevation() |
|---|
| 97 | else: |
|---|
| 98 | return self.invert_elevation |
|---|
| 99 | |
|---|
| 100 | |
|---|
| 101 | def get_enquiry_velocity(self): |
|---|
| 102 | |
|---|
| 103 | depth = self.get_enquiry_water_depth() |
|---|
| 104 | u = depth*self.get_enquiry_xmom()/(depth**2 + velocity_protection) |
|---|
| 105 | v = depth*self.get_enquiry_ymom()/(depth**2 + velocity_protection) |
|---|
| 106 | |
|---|
| 107 | return u, v |
|---|
| 108 | |
|---|
| 109 | |
|---|
| 110 | def get_enquiry_xvelocity(self): |
|---|
| 111 | |
|---|
| 112 | depth = self.get_enquiry_water_depth() |
|---|
| 113 | return depth*self.get_enquiry_xmom()/(depth**2 + velocity_protection) |
|---|
| 114 | |
|---|
| 115 | def get_enquiry_yvelocity(self): |
|---|
| 116 | |
|---|
| 117 | depth = self.get_enquiry_water_depth() |
|---|
| 118 | return depth*self.get_enquiry_ymom()/(depth**2 + velocity_protection) |
|---|
| 119 | |
|---|
| 120 | |
|---|
| 121 | def get_enquiry_speed(self): |
|---|
| 122 | |
|---|
| 123 | u, v = self.get_enquiry_velocity() |
|---|
| 124 | |
|---|
| 125 | return math.sqrt(u**2 + v**2) |
|---|
| 126 | |
|---|
| 127 | |
|---|
| 128 | def get_enquiry_velocity_head(self): |
|---|
| 129 | |
|---|
| 130 | return 0.5*self.get_enquiry_speed()**2/g |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | def get_enquiry_total_energy(self): |
|---|
| 134 | |
|---|
| 135 | return self.get_enquiry_velocity_head() + self.get_enquiry_stage() |
|---|
| 136 | |
|---|
| 137 | |
|---|
| 138 | def get_enquiry_specific_energy(self): |
|---|
| 139 | |
|---|
| 140 | return self.get_enquiry_velocity_head() + self.get_enquiry_depth() |
|---|
| 141 | |
|---|
| 142 | |
|---|
