Changeset 5570 for anuga_core/source/anuga/culvert_flows/culvert_class.py
 Timestamp:
 Jul 25, 2008, 1:35:45 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/culvert_flows/culvert_class.py
r5453 r5570 2 2 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons 3 3 from anuga.utilities.system_tools import log_to_file 4 from anuga.utilities.polygon import inside_polygon 5 from anuga.utilities.polygon import is_inside_polygon 6 4 7 5 8 class Culvert_flow: … … 125 128 # P['enquiry_polygon1']], 126 129 # figname='culvert_polygon_output') 127 130 131 132 # Check that all polygons lie within the mesh. 133 bounding_polygon = domain.get_boundary_polygon() 134 for key in P.keys(): 135 print 'Key', key 136 for point in P[key]: 137 138 print 'Passing in:', point 139 msg = 'Point %s in polygon %s for culvert %s did not'\ 140 %(str(point), key, self.label) 141 msg += 'fall within the domain boundary.' 142 assert is_inside_polygon(point, bounding_polygon), msg 143 144 145 # Create inflow object at each end of the culvert. 128 146 self.openings = [] 129 147 self.openings.append(Inflow(domain, … … 188 206 def __call__(self, domain): 189 207 from anuga.utilities.numerical_tools import mean 190 from anuga.utilities.polygon import inside_polygon208 191 209 from anuga.config import g, epsilon 192 210 from Numeric import take, sqrt … … 209 227 openings = self.openings # There are two Opening [0] and [1] 210 228 for i, opening in enumerate(openings): 211 stage = domain.quantities['stage'].get_values(location='centroids', 212 indices=opening.exchange_indices) 213 elevation = domain.quantities['elevation'].get_values(location='centroids', 214 indices=opening.exchange_indices) 229 dq = domain.quantities 230 231 stage = dq['stage'].get_values(location='centroids', 232 indices=opening.exchange_indices) 233 elevation = dq['elevation'].get_values(location='centroids', 234 indices=opening.exchange_indices) 215 235 216 236 # Indices corresponding to energy enquiry field for this opening 217 237 coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y) 218 enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 238 enquiry_indices = inside_polygon(coordinates, 239 self.enquiry_polygons[i]) 219 240 220 241 if len(enquiry_indices) == 0: … … 223 244 224 245 # Get model values for points in enquiry polygon for this opening 225 dq = domain.quantities 226 stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices) 227 xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices) 228 ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices) 229 elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices) 246 stage = dq['stage'].get_values(location='centroids', 247 indices=enquiry_indices) 248 xmomentum = dq['xmomentum'].get_values(location='centroids', 249 indices=enquiry_indices) 250 ymomentum = dq['ymomentum'].get_values(location='centroids', 251 indices=enquiry_indices) 252 elevation = dq['elevation'].get_values(location='centroids', 253 indices=enquiry_indices) 230 254 depth = stage  elevation 231 255 … … 234 258 ux = xmomentum/(depth+velocity_protection/depth) # Velocity (xdirection) 235 259 uy = ymomentum/(depth+velocity_protection/depth) # Velocity (ydirection) 260 print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum 236 261 v = mean(sqrt(ux**2+uy**2)) # Average velocity 237 262 w = mean(stage) # Average stage
Note: See TracChangeset
for help on using the changeset viewer.