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