# Changeset 7958

Ignore:
Timestamp:
Aug 20, 2010, 2:02:43 PM (12 years ago)
Message:

Corrections made to the code

File:
1 edited

Unmodified
Added
Removed
• ## trunk/anuga_core/source/anuga/structures/culvert_operator.py

 r7957 from anuga.shallow_water.forcing import Inflow, General_forcing from anuga.structures.culvert_polygons import create_culvert_polygons from anuga.utilities.system_tools import log_to_file from anuga.geometry.polygon import inside_polygon, is_inside_polygon import numpy as num from math import sqrt from math import sqrt class Below_interval(Exception): pass class Generic_box_culvert: """Culvert flow - transfer water from one rectngular box to another. """Culvert flow - transfer water from one rectangular box to another. Sets up the geometry of problem # Input check self.domain = domain self.end_points= [end_point0, end_point1] self.enquiry_gap_factor = enquiry_gap_factor if height is None: height = width self.width = width self.height = height self.width = width self.domain = domain self.end_points= [end_point0, end_point1] self.enquiry_gap_factor=enquiry_gap_factor self.verbose=verbose self.filename = None # Create the fundamental culvert polygons from polygon self.create_culvert_polygons() self.check_culvert_inside_domain() # Establish initial values at each enquiry point self.enquiry_quantity_values = [] dq = domain.quantities for i in [0,1]: for i in [0, 1]: idx = self.enquiry_indices[i] elevation = dq['elevation'].get_values(location='centroids', stage = dq['stage'].get_values(location='centroids', indices=[idx])[0] depth = stage-elevation depth = stage - elevation quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth } self.enquiry_quantity_values.append(quantity_values) assert self.culvert_length > 0.0 def set_store_hydrograph_discharge(self,filename=None): def set_store_hydrograph_discharge(self, filename=None): if filename is None: fid.write('time, discharge\n') fid.close() def create_culvert_polygons(self): x1, y1 = self.end_points[1] dx = x1-x0 dy = y1-y0 dx = x1 - x0 dy = y1 - y0 self.culvert_vector = num.array([dx, dy]) self.culvert_length = sqrt(num.sum(self.culvert_vector**2)) assert self.culvert_length > 0.0, 'The length of culvert is less than 0' # Unit direction vector and normal # Short hands w = 0.5*width*normal # Perpendicular vector of 1/2 width h = height*vector    # Vector of length=height in the w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width h = self.height*self.culvert_vector    # Vector of length=height in the # direction of the culvert gap = (1 + enquiry_gap_factor)*h gap = (1 + self.enquiry_gap_factor)*h self.exchange_polygons = [] # Build exchange polygon and enquiry points 0 and 1 for i in [0,1]: for i in [0, 1]: p0 = self.end_points[i] + w p1 = self.end_point[i] - w p1 = self.end_points[i] - w p2 = p1 - h p3 = p0 - h self.exchange_polygons.append(num.array([p0,p1,p2,p3])) self.enquiry_points.append(end_point0 - gap) self.exchange_polygons.append(num.array([p0, p1, p2, p3])) self.enquiry_points.append(self.end_points[i] + (2*i-1)*gap) self.polygon_areas = [] """ bounding_polygon = self.domain.get_boundary_polygon() P = self.culvert_polygons for key in P.keys(): if key in ['exchange_polygon0', 'exchange_polygon1']: for point in list(P[key]) + self.enquiry_points: msg = 'Point %s in polygon %s for culvert %s did not'\ %(str(point), key, self.label) msg += 'fall within the domain boundary.' assert is_inside_polygon(point, bounding_polygon), msg for i in [0, 1]: for point in list(self.exchange_polygons[i]) + self.enquiry_points: msg = 'Point %s did not '\ %(str(point)) msg += 'fall within the domain boundary.' assert is_inside_polygon(point, bounding_polygon), msg return label, type, width, height, length, number_of_barrels, description, rating_curve Culvert_flow = Culvert_flow_general
Note: See TracChangeset for help on using the changeset viewer.