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