Changeset 7958


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

Corrections made to the code

File:
1 edited

Legend:

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

    r7957 r7958  
    22
    33from anuga.shallow_water.forcing import Inflow, General_forcing
    4 from anuga.structures.culvert_polygons import create_culvert_polygons
    54from anuga.utilities.system_tools import log_to_file
    65from anuga.geometry.polygon import inside_polygon, is_inside_polygon
     
    1716import numpy as num
    1817from math import sqrt
    19 from math import sqrt
    2018
    2119class Below_interval(Exception): pass
     
    2523
    2624class Generic_box_culvert:
    27     """Culvert flow - transfer water from one rectngular box to another.
     25    """Culvert flow - transfer water from one rectangular box to another.
    2826    Sets up the geometry of problem
    2927   
     
    4644        # Input check
    4745       
     46        self.domain = domain
     47
     48        self.end_points= [end_point0, end_point1]
     49        self.enquiry_gap_factor = enquiry_gap_factor
     50       
    4851        if height is None:
    4952            height = width
    5053
     54        self.width = width
    5155        self.height = height
    52         self.width = width
    53        
    54         self.domain = domain
    55         self.end_points= [end_point0, end_point1]
    56         self.enquiry_gap_factor=enquiry_gap_factor
    57 
     56       
    5857        self.verbose=verbose
    5958        self.filename = None
    6059       
    61 
    6260        # Create the fundamental culvert polygons from polygon
    6361        self.create_culvert_polygons()
     
    6563        self.check_culvert_inside_domain()
    6664
    67 
    6865        # Establish initial values at each enquiry point
    6966        self.enquiry_quantity_values = []
    7067        dq = domain.quantities
    71         for i in [0,1]:
     68        for i in [0, 1]:
    7269            idx = self.enquiry_indices[i]
    7370            elevation = dq['elevation'].get_values(location='centroids',
     
    7572            stage = dq['stage'].get_values(location='centroids',
    7673                                           indices=[idx])[0]
    77             depth = stage-elevation
     74            depth = stage - elevation
    7875           
    7976            quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth }
    8077            self.enquiry_quantity_values.append(quantity_values)
    81        
    82        
    83         assert self.culvert_length > 0.0
    84 
    85 
    86     def set_store_hydrograph_discharge(self,filename=None):
     78
     79    def set_store_hydrograph_discharge(self, filename=None):
    8780
    8881        if filename is None:
     
    9790        fid.write('time, discharge\n')
    9891        fid.close()
    99 
    100 
    101 
    10292
    10393    def create_culvert_polygons(self):
     
    111101        x1, y1 = self.end_points[1]
    112102
    113         dx = x1-x0
    114         dy = y1-y0
     103        dx = x1 - x0
     104        dy = y1 - y0
    115105
    116106        self.culvert_vector = num.array([dx, dy])
    117107        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
    118 
     108        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
    119109
    120110        # Unit direction vector and normal
     
    123113
    124114        # Short hands
    125         w = 0.5*width*normal # Perpendicular vector of 1/2 width
    126         h = height*vector    # Vector of length=height in the
     115        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
     116        h = self.height*self.culvert_vector    # Vector of length=height in the
    127117                             # direction of the culvert
    128         gap = (1 + enquiry_gap_factor)*h
    129 
     118        gap = (1 + self.enquiry_gap_factor)*h
    130119
    131120        self.exchange_polygons = []
     
    133122
    134123        # Build exchange polygon and enquiry points 0 and 1
    135         for i in [0,1]:
     124        for i in [0, 1]:
    136125            p0 = self.end_points[i] + w
    137             p1 = self.end_point[i] - w
     126            p1 = self.end_points[i] - w
    138127            p2 = p1 - h
    139128            p3 = p0 - h
    140             self.exchange_polygons.append(num.array([p0,p1,p2,p3]))
    141             self.enquiry_points.append(end_point0 - gap)
     129            self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
     130            self.enquiry_points.append(self.end_points[i] + (2*i-1)*gap)
    142131
    143132        self.polygon_areas = []
     
    239228        """
    240229        bounding_polygon = self.domain.get_boundary_polygon()
    241         P = self.culvert_polygons
    242         for key in P.keys():
    243             if key in ['exchange_polygon0',
    244                        'exchange_polygon1']:
    245                 for point in list(P[key]) + self.enquiry_points:
    246                     msg = 'Point %s in polygon %s for culvert %s did not'\
    247                         %(str(point), key, self.label)
    248                     msg += 'fall within the domain boundary.'
    249                     assert is_inside_polygon(point, bounding_polygon), msg
     230        for i in [0, 1]:
     231            for point in list(self.exchange_polygons[i]) + self.enquiry_points:
     232                msg = 'Point %s did not '\
     233                    %(str(point))
     234                msg += 'fall within the domain boundary.'
     235                assert is_inside_polygon(point, bounding_polygon), msg
    250236           
    251237
     
    719705    return label, type, width, height, length, number_of_barrels, description, rating_curve
    720706
    721            
    722 
    723 Culvert_flow = Culvert_flow_general       
Note: See TracChangeset for help on using the changeset viewer.