Ignore:
Timestamp:
Aug 26, 2010, 5:17:24 PM (12 years ago)
Author:
habili
Message:

Refactoring of routines.

File:
1 edited

Legend:

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

    r7976 r7977  
    22from anuga.config import g
    33import anuga.utilities.log as log
    4 import inlet
    5 import numpy as num
    6 import math
     4import box_culvert
    75
    8 class Below_interval(Exception): pass
    9 class Above_interval(Exception): pass
    10 
    11 
    12 class Generic_box_culvert:
     6class Culvert_operator:
    137    """Culvert flow - transfer water from one rectangular box to another.
    148    Sets up the geometry of problem
     
    2923                 verbose=False):
    3024       
    31         # Input check
    32        
    3325        self.domain = domain
    34 
    3526        self.domain.set_fractional_step_operator(self)
    36 
    37         self.end_points = [end_point0, end_point1]
     27        end_points = [end_point0, end_point1]
    3828       
    3929        if height is None:
     
    4232        self.width = width
    4333        self.height = height
    44        
    45         self.verbose=verbose
    4634       
    47         # Create the fundamental culvert polygons and create inlet objects
    48         self.create_culvert_polygons()
    49 
    50         #FIXME (SR) Put this into a foe loop to deal with more inlets
    51         self.inlets = []
    52         polygon0 = self.inlet_polygons[0]
    53         inlet0_vector = self.culvert_vector
    54         self.inlets.append(inlet.Inlet(self.domain, polygon0))
    55 
    56         polygon1 = self.inlet_polygons[1]
    57         inlet1_vector = - self.culvert_vector
    58         self.inlets.append(inlet.Inlet(self.domain, polygon1))
     35        self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height)
     36        self.inlets = self.culvert.get_inlets()
    5937   
    6038        self.print_stats()
     
    6341    def __call__(self):
    6442
    65         # Time stuff
    66         time     = self.domain.get_time()
    6743        timestep = self.domain.get_timestep()
    68 
    69 
    7044
    7145        inflow  = self.inlets[0]
     
    8155
    8256        delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
    83         culvert_slope = delta_z/self.culvert_length
     57        culvert_slope = delta_z/self.culvert.get_culvert_length()
    8458
    8559        # Determine controlling energy (driving head) for culvert
     
    9165            driving_head = inflow.get_average_specific_energy()
    9266           
    93 
    9467        # Transfer
    95         from culvert_routines import boyd_generalised_culvert_model
     68        from culvert_routines import boyd_box, boyd_circle
    9669        Q, barrel_velocity, culvert_outlet_depth =\
    97                     boyd_generalised_culvert_model(inflow.get_average_height(),
     70                              boyd_circle(inflow.get_average_height(),
    9871                                         outflow.get_average_height(),
    9972                                         inflow.get_average_speed(),
     
    10174                                         inflow.get_average_specific_energy(),
    10275                                         delta_total_energy,
    103                                          g,
    104                                          culvert_length=self.culvert_length,
     76                                         culvert_length=self.culvert.get_culvert_length(),
    10577                                         culvert_width=self.width,
    10678                                         culvert_height=self.height,
    107                                          culvert_type='box',
    10879                                         manning=0.01)
    10980
    11081        transfer_water = Q*timestep
    111 
    11282
    11383        inflow.set_heights(inflow.get_average_height() - transfer_water)
     
    11585        inflow.set_ymoms(0.0)
    11686
    117 
    11887        outflow.set_heights(outflow.get_average_height() + transfer_water)
    11988        outflow.set_xmoms(0.0)
    12089        outflow.set_ymoms(0.0)
    121 
    122 
    12390
    12491    def print_stats(self):
     
    141108
    142109        print '====================================='
    143 
    144 
    145 
    146 
    147 
    148     def create_culvert_polygons(self):
    149 
    150         """Create polygons at the end of a culvert inlet and outlet.
    151         At either end two polygons will be created; one for the actual flow to pass through and one a little further away
    152         for enquiring the total energy at both ends of the culvert and transferring flow.
    153         """
    154 
    155         # Calculate geometry
    156         x0, y0 = self.end_points[0]
    157         x1, y1 = self.end_points[1]
    158 
    159         dx = x1 - x0
    160         dy = y1 - y0
    161 
    162         self.culvert_vector = num.array([dx, dy])
    163         self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
    164         assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
    165 
    166         # Unit direction vector and normal
    167         self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
    168         self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
    169 
    170         # Short hands
    171         w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
    172         h = self.height*self.culvert_vector    # Vector of length=height in the
    173                              # direction of the culvert
    174 
    175         self.inlet_polygons = []
    176 
    177         # Build exchange polygon and enquiry points 0 and 1
    178         for i in [0, 1]:
    179             i0 = (2*i-1)
    180             p0 = self.end_points[i] + w
    181             p1 = self.end_points[i] - w
    182             p2 = p1 + i0*h
    183             p3 = p0 + i0*h
    184             self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
    185 
    186         # Check that enquiry points are outside inlet polygons
    187         for i in [0,1]:
    188             polygon = self.inlet_polygons[i]
    189             # FIXME (SR) Probably should calculate the area of all the triangles
    190             # associated with this polygon, as there is likely to be some
    191             # inconsistency between triangles and ploygon
    192             area = polygon_area(polygon)
    193            
    194             msg = 'Polygon %s ' %(polygon)
    195             msg += ' has area = %f' % area
    196             assert area > 0.0, msg
    197    
    198110
    199111                       
Note: See TracChangeset for help on using the changeset viewer.