Changeset 8056


Ignore:
Timestamp:
Oct 29, 2010, 11:08:43 AM (12 years ago)
Author:
habili
Message:

The variables end_points, exchange_lines and enquiry_points are now all user defined. Either end_points or exchange_lines must be defined. Use example in test_boyd_box_operator.py

Location:
trunk/anuga_core/source/anuga/structures
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/structures/boyd_box_operator.py

    r8035 r8056  
    1717    def __init__(self,
    1818                 domain,
    19                  end_point0,
    20                  end_point1,
    2119                 losses,
    2220                 width,
    2321                 height=None,
    24                  apron=None,
     22                 end_points=None,
     23                 exchange_lines=None,
     24                 enquiry_points=None,
     25                 apron=0.1,
    2526                 manning=0.013,
    2627                 enquiry_gap=0.0,
    … …  
    3233                 logging=False,
    3334                 verbose=False):
    34 
    3535                     
    3636        anuga.Structure_operator.__init__(self,
    3737                                          domain,
    38                                           end_point0,
    39                                           end_point1,
    4038                                          width,
    4139                                          height,
     40                                          end_points,
     41                                          exchange_lines,
     42                                          enquiry_points,
    4243                                          apron,
    4344                                          manning,

  • trunk/anuga_core/source/anuga/structures/boyd_pipe_operator.py

    r8034 r8056  
    1616    def __init__(self,
    1717                 domain,
    18                  end_point0,
    19                  end_point1,
    2018                 losses,
     19                 end_points=None,
     20                 exchange_lines=None,
     21                 enquiry_points=None,
    2122                 diameter=None,
    22                  apron=None,
     23                 apron=0.1,
    2324                 manning=0.013,
    2425                 enquiry_gap=0.2,
    … …  
    3031                 logging=False,
    3132                 verbose=False):
    32 
    33 
    3433                     
    3534        anuga.Structure_operator.__init__(self,
    3635                                          domain,
    37                                           end_point0,
    38                                           end_point1,
     36                                          end_points,
     37                                          exchange_lines,
     38                                          enquiry_points,
    3939                                          width=diameter,
    4040                                          height=None,

  • trunk/anuga_core/source/anuga/structures/inlet.py

    r8050 r8056  
    1 from anuga.geometry.polygon import inside_polygon, is_inside_polygon, polyline_overlap
     1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
    22from anuga.config import velocity_protection, g
    33import math
    … …  
    99    """
    1010
    11     def __init__(self, domain, polyline, verbose=False):
     11    def __init__(self, domain, line, verbose=False):
    1212
    1313        self.domain = domain
    1414        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
    15         self.polyline = polyline
     15        self.line = line
    1616        self.verbose = verbose
    1717
    … …  
    2727        vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
    2828
    29         # Check that polyline lies within the mesh.
    30         for point in self.polyline: 
     29        # Check that line lies within the mesh.
     30        for point in self.line: 
    3131                msg = 'Point %s ' %  str(point)
    3232                msg += ' did not fall within the domain boundary.'
    … …  
    3535
    3636
    37         self.triangle_indices = polyline_overlap(vertex_coordinates, self.polyline)
     37        self.triangle_indices = line_intersect(vertex_coordinates, self.line)
    3838
    3939        if len(self.triangle_indices) == 0:
    40             msg = 'Inlet polyline=%s ' % (self.polyline)
    41             msg += 'No triangles intersecting polyline '
     40            msg = 'Inlet line=%s ' % (self.line)
     41            msg += 'No triangles intersecting line '
    4242            raise Exception, msg
    4343
    … …  
    4747       
    4848        # Compute inlet area as the sum of areas of triangles identified
    49         # by polyline. Must be called after compute_inlet_triangle_indices().
     49        # by line. Must be called after compute_inlet_triangle_indices().
    5050        if len(self.triangle_indices) == 0:
    51             region = 'Inlet polyline=%s' % (self.inlet_polyline)
     51            region = 'Inlet line=%s' % (self.inlet_line)
    5252            msg = 'No triangles have been identified in region '
    5353            raise Exception, msg

  • trunk/anuga_core/source/anuga/structures/inlet_enquiry.py

    r8050 r8056  
    1 from anuga.geometry.polygon import inside_polygon, is_inside_polygon, polyline_overlap
     1from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
    22from anuga.config import velocity_protection, g
    33import math

  • trunk/anuga_core/source/anuga/structures/structure_operator.py

    r8050 r8056  
    2222    def __init__(self,
    2323                 domain,
    24                  end_point0,
    25                  end_point1,
    2624                 width,
    2725                 height,
     26                 end_points,
     27                 exchange_lines,
     28                 enquiry_points,
    2829                 apron,
    2930                 manning,
    … …  
    3738        self.domain = domain
    3839        self.domain.set_fractional_step_operator(self)
    39         self.end_points = [end_point0, end_point1]
    40 
    41        
     40        self.end_points = end_points
     41        self.exchange_lines = exchange_lines
     42        self.enquiry_points = enquiry_points
    4243       
    4344        if height is None:
    … …  
    5859            self.description = description
    5960       
    60 
    6161        if label == None:
    6262            self.label = "structure_%g" % Structure_operator.counter
    … …  
    6464            self.label = label + '_%g' % Structure_operator.counter
    6565
    66 
    6766        if structure_type == None:
    6867            self.structure_type = 'generic structure'
    … …  
    7069            self.structure_type = structure_type
    7170           
    72         self.verbose = verbose
    73 
    74        
     71        self.verbose = verbose       
    7572       
    7673        # Keep count of structures
    … …  
    8380        self.driving_energy = 0.0
    8481       
    85         self.__create_exchange_polylines()
    86 
     82        if exchange_lines is not None:
     83            self.__process_skew_culvert()
     84        elif end_points is not None:
     85            self.__process_non_skew_culvert()
     86        else:
     87            raise Exception, 'Define either exchange_lines or end_points'
     88       
    8789        self.inlets = []
    88         polyline0 = self.inlet_polylines[0]
    89         enquiry_point0 = self.inlet_equiry_points[0]
     90        line0 = self.exchange_lines[0] #self.inlet_lines[0]
     91        enquiry_point0 = self.enquiry_points[0]
    9092        outward_vector0 = self.culvert_vector
    91         self.inlets.append(inlet_enquiry.Inlet_enquiry(self.domain, polyline0,
     93        self.inlets.append(inlet_enquiry.Inlet_enquiry(self.domain, line0,
    9294                           enquiry_point0, outward_vector0, self.verbose))
    9395
    94         polyline1 = self.inlet_polylines[1]
    95         enquiry_point1 = self.inlet_equiry_points[1]
     96        line1 = self.exchange_lines[1]
     97        enquiry_point1 = self.enquiry_points[1]
    9698        outward_vector1  = - self.culvert_vector
    97         self.inlets.append(inlet_enquiry.Inlet_enquiry(self.domain, polyline1,
     99        self.inlets.append(inlet_enquiry.Inlet_enquiry(self.domain, line1,
    98100                           enquiry_point1, outward_vector1, self.verbose))
    99101
    … …  
    185187            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
    186188
    187 
    188189        self.inflow  = self.inlets[0]
    189190        self.outflow = self.inlets[1]
    190        
    191191
    192192        if self.delta_total_energy < 0:
    … …  
    196196
    197197
    198     def __create_exchange_polylines(self):
    199 
    200         """Create polylines at the end of a culvert inlet and outlet.
    201         At either end two polylines will be created; one for the actual flow to pass through and one a little further away
     198    def __process_non_skew_culvert(self):
     199
     200        """Create lines at the end of a culvert inlet and outlet.
     201        At either end two lines will be created; one for the actual flow to pass through and one a little further away
    202202        for enquiring the total energy at both ends of the culvert and transferring flow.
    203203        """
    204 
    205         # Calculate geometry
    206         x0, y0 = self.end_points[0]
    207         x1, y1 = self.end_points[1]
    208 
    209         dx = x1 - x0
    210         dy = y1 - y0
    211 
    212         self.culvert_vector = num.array([dx, dy])
    213         self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
     204       
     205        self.culvert_vector = self.end_points[1] - self.end_points[0]
     206        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))   
    214207        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
    215 
    216         # Unit direction vector and normal
    217         self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
    218         self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
    219 
    220         # Short hands
    221         w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
    222         #h = self.apron*self.culvert_vector    # Vector of length=height in the
    223                              # direction of the culvert
    224 
    225         #gap = 1.5*h + self.enquiry_gap
    226         gap = (self.apron+ self.enquiry_gap)*self.culvert_vector
    227 
    228         self.inlet_polylines = []
    229         self.inlet_equiry_points = []
     208       
     209        self.culvert_vector /= self.culvert_length
     210       
     211        culvert_normal = num.array([-self.culvert_vector[1], self.culvert_vector[0]])  # Normal vector
     212        w = 0.5*self.width*culvert_normal # Perpendicular vector of 1/2 width
     213
     214        self.exchange_lines = []
    230215
    231216        # Build exchange polyline and enquiry point
    232         for i in [0, 1]:
    233             i0 = (2*i-1) #i0 determines the sign of the points
    234             p0 = self.end_points[i] + w
    235             p1 = self.end_points[i] - w
    236             ep = self.end_points[i] + i0*gap
    237 
    238             self.inlet_polylines.append(num.array([p0, p1]))
    239             self.inlet_equiry_points.append(ep)           
    240    
     217        if self.enquiry_points is None:
     218           
     219            gap = (self.apron + self.enquiry_gap)*self.culvert_vector
     220            self.enquiry_points = []
     221           
     222            for i in [0, 1]:
     223                p0 = self.end_points[i] + w
     224                p1 = self.end_points[i] - w
     225                self.exchange_lines.append(num.array([p0, p1]))
     226                ep = self.end_points[i] + (2*i - 1)*gap #(2*i - 1) determines the sign of the points
     227                self.enquiry_points.append(ep)
     228           
     229        else:           
     230            for i in [0, 1]:
     231                p0 = self.end_points[i] + w
     232                p1 = self.end_points[i] - w
     233                self.exchange_lines.append(num.array([p0, p1]))
     234           
     235 
     236    def __process_skew_culvert(self):   
     237       
     238        """Compute skew culvert.
     239        If exchange lines are given, the enquiry points are determined. This is for enquiring
     240        the total energy at both ends of the culvert and transferring flow.
     241        """
     242           
     243        centre_point0 = 0.5*(self.exchange_lines[0][0] + self.exchange_lines[0][1])
     244        centre_point1 = 0.5*(self.exchange_lines[1][0] + self.exchange_lines[1][1])
     245       
     246        if self.end_points is None:
     247            self.culvert_vector = centre_point1 - centre_point0
     248        else:
     249            self.culvert_vector = self.end_points[1] - end_points[0]
     250       
     251        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))   
     252        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
     253       
     254        if self.enquiry_points is None:
     255       
     256            self.culvert_vector /= self.culvert_length
     257            gap = (self.apron + self.enquiry_gap)*self.culvert_vector
     258       
     259            self.enquiry_points = []
     260
     261            self.enquiry_points.append(centre_point0 - gap)
     262            self.enquiry_points.append(centre_point1 + gap)
     263           
     264       
    241265    def discharge_routine(self):
    242266       
    … …  
    269293            message += '\n'
    270294
    271             message += 'polyline\n'
    272             message += '%s' % inlet.polyline
     295            message += 'line\n'
     296            message += '%s' % inlet.line
    273297            message += '\n'
    274298

  • trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py

    r8036 r8056  
    2626
    2727
    28     def test_boyd_1(self):
    29         """test_boyd_1
     28    def test_boyd_non_skew(self):
     29        """test_boyd_non_skew
    3030       
    3131        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson,
    … …  
    4444        culvert_mannings = 0.013
    4545       
    46         culvert_apron = 0.001
    47         enquiry_gap = 1.0
     46        culvert_apron = 0.0
     47        enquiry_gap = 10.0
    4848       
    4949        expected_Q = 4.55
    … …  
    5656        domain_length = 200.  #x-Dir
    5757        domain_width  = 200.  #y-dir
    58         dx = dy = 10.0          # Resolution: Length of subdivisions on both axes
     58        dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
    5959
    6060
    … …  
    101101        print 'Defining Structures'
    102102       
     103        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
     104        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
     105       
     106       
    103107        culvert = Boyd_box_operator(domain,
    104             label='3.6x3.6RCBC',
    105             end_point0=[domain_length/2-culvert_length/2, 100.0],
    106             end_point1=[domain_length/2+culvert_length/2, 100.0],
    107             losses=culvert_losses,
    108             width=culvert_width,
    109             height=culvert_height,
    110             apron=culvert_apron,
    111             enquiry_gap=enquiry_gap,
    112             use_momentum_jet=False,
    113             use_velocity_head=False,
    114             manning=culvert_mannings,
    115             verbose=False)
    116 
     108                                    losses=culvert_losses,
     109                                    width=culvert_width,
     110                                    end_points=[ep0, ep1],
     111                                    height=culvert_height,
     112                                    apron=culvert_apron,
     113                                    enquiry_gap=enquiry_gap,
     114                                    use_momentum_jet=False,
     115                                    use_velocity_head=False,
     116                                    manning=culvert_mannings,
     117                                    label='3.6x3.6RCBC',
     118                                    verbose=False)
    117119
    118120        culvert.determine_inflow_outflow()
    119121       
    120122        ( Q, v, d ) = culvert.discharge_routine()
     123       
     124        print 'test_boyd_non_skew'
     125        print 'Q: ', Q, 'expected_Q: ', expected_Q
    121126       
    122127
    … …  
    125130        assert numpy.allclose(d, expected_d, rtol=1.0e-2) #depth at outlet used to calc v
    126131       
     132       
     133    def test_boyd_skew(self):
     134        """test_boyd_skew
     135       
     136        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson,
     137        calculation code by MJ Boyd
     138        """
     139
     140        stage_0 = 11.0
     141        stage_1 = 10.0
     142        elevation_0 = 10.0
     143        elevation_1 = 10.0
     144
     145        culvert_length = 20.0
     146        culvert_width = 3.66
     147        culvert_height = 3.66
     148        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
     149        culvert_mannings = 0.013
     150       
     151        culvert_apron = 0.0
     152        enquiry_gap = 10.0
     153       
     154        expected_Q = 4.55
     155        expected_v = 2.3
     156        expected_d = 0.54
     157       
     158
     159        # Probably no need to change below here
     160       
     161        domain_length = 200.  #x-Dir
     162        domain_width  = 200.  #y-dir
     163        dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
     164
     165
     166        points, vertices, boundary = rectangular_cross(int(domain_length/dx), int(domain_width/dy),
     167                                                        len1=domain_length, len2=domain_width)
     168        domain = Domain(points, vertices, boundary)   
     169        domain.set_name('Test_Outlet_Inlet')                 # Output name
     170        domain.set_default_order(2)
     171        domain.H0 = 0.01
     172        domain.tight_slope_limiters = 1
     173
     174        print 'Size', len(domain)
     175
     176        #------------------------------------------------------------------------------
     177        # Setup initial conditions
     178        #------------------------------------------------------------------------------
     179
     180        def elevation(x, y):
     181            """Set up a elevation
     182            """
     183           
     184            z = numpy.zeros(x.shape,dtype='d')
     185            z[:] = elevation_0
     186           
     187            numpy.putmask(z, x > domain_length/2, elevation_1)
     188   
     189            return z
     190           
     191        def stage(x,y):
     192            """Set up stage
     193            """
     194            z = numpy.zeros(x.shape,dtype='d')
     195            z[:] = stage_0
     196           
     197            numpy.putmask(z, x > domain_length/2, stage_1)
     198
     199            return z
     200           
     201        print 'Setting Quantities....'
     202        domain.set_quantity('elevation', elevation)  # Use function for elevation
     203        domain.set_quantity('stage',  stage)   # Use function for elevation
     204
     205
     206        print 'Defining Structures'
     207       
     208        a = domain_length/2 - culvert_length/2
     209        b = domain_length/2 + culvert_length/2
     210       
     211        el0 = numpy.array([[a, 100.0 - culvert_width/2], [a, 100.0 + culvert_width/2]])
     212        el1 = numpy.array([[b, 100.0 - culvert_width/2], [b, 100.0 + culvert_width/2]])
     213       
     214        culvert = Boyd_box_operator(domain,
     215                                    losses=culvert_losses,
     216                                    width=culvert_width,
     217                                    exchange_lines=[el0, el1],
     218                                    height=culvert_height,
     219                                    apron=culvert_apron,
     220                                    enquiry_gap=enquiry_gap,
     221                                    use_momentum_jet=False,
     222                                    use_velocity_head=False,
     223                                    manning=culvert_mannings,
     224                                    label='3.6x3.6RCBC',
     225                                    verbose=False)
     226
     227        culvert.determine_inflow_outflow()
     228       
     229        ( Q, v, d ) = culvert.discharge_routine()
     230       
     231        print 'test_boyd_skew'
     232        print 'Q: ', Q, 'expected_Q: ', expected_Q
     233
     234        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow
     235        assert numpy.allclose(v, expected_v, rtol=1.0e-2) #outflow velocity
     236        assert numpy.allclose(d, expected_d, rtol=1.0e-2) #depth at outlet used to calc v         
     237
     238
     239    def test_boyd_non_skew_enquiry_points(self):
     240        """test_boyd_skew_enquiry_points
     241       
     242        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson,
     243        calculation code by MJ Boyd
     244        """
     245
     246        stage_0 = 11.0
     247        stage_1 = 10.0
     248        elevation_0 = 10.0
     249        elevation_1 = 10.0
     250
     251        culvert_length = 20.0
     252        culvert_width = 3.66
     253        culvert_height = 3.66
     254        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
     255        culvert_mannings = 0.013
     256       
     257        culvert_apron = 0.0
     258        enquiry_gap = 10.0
     259       
     260        expected_Q = 4.55
     261        expected_v = 2.3
     262        expected_d = 0.54
     263       
     264
     265        # Probably no need to change below here
     266       
     267        domain_length = 200.  #x-Dir
     268        domain_width  = 200.  #y-dir
     269        dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
     270
     271
     272        points, vertices, boundary = rectangular_cross(int(domain_length/dx), int(domain_width/dy),
     273                                                        len1=domain_length, len2=domain_width)
     274        domain = Domain(points, vertices, boundary)   
     275        domain.set_name('Test_Outlet_Inlet')                 # Output name
     276        domain.set_default_order(2)
     277        domain.H0 = 0.01
     278        domain.tight_slope_limiters = 1
     279
     280        print 'Size', len(domain)
     281
     282        #------------------------------------------------------------------------------
     283        # Setup initial conditions
     284        #------------------------------------------------------------------------------
     285
     286        def elevation(x, y):
     287            """Set up a elevation
     288            """
     289           
     290            z = numpy.zeros(x.shape,dtype='d')
     291            z[:] = elevation_0
     292           
     293            numpy.putmask(z, x > domain_length/2, elevation_1)
     294   
     295            return z
     296           
     297        def stage(x,y):
     298            """Set up stage
     299            """
     300            z = numpy.zeros(x.shape,dtype='d')
     301            z[:] = stage_0
     302           
     303            numpy.putmask(z, x > domain_length/2, stage_1)
     304
     305            return z
     306           
     307        print 'Setting Quantities....'
     308        domain.set_quantity('elevation', elevation)  # Use function for elevation
     309        domain.set_quantity('stage',  stage)   # Use function for elevation
     310
     311
     312        print 'Defining Structures'
     313       
     314        a = domain_length/2 - culvert_length/2
     315        b = domain_length/2 + culvert_length/2
     316       
     317        el0 = numpy.array([[a, 100.0 - culvert_width/2], [a, 100.0 + culvert_width/2]])
     318        el1 = numpy.array([[b, 100.0 - culvert_width/2], [b, 100.0 + culvert_width/2]])
     319       
     320        enquiry_points = (numpy.array([85, 100]), numpy.array([115, 100]))
     321       
     322        culvert = Boyd_box_operator(domain,
     323                                    losses=culvert_losses,
     324                                    width=culvert_width,
     325                                    exchange_lines=[el0, el1],
     326                                    enquiry_points=enquiry_points,
     327                                    height=culvert_height,
     328                                    apron=culvert_apron,
     329                                    enquiry_gap=enquiry_gap,
     330                                    use_momentum_jet=False,
     331                                    use_velocity_head=False,
     332                                    manning=culvert_mannings,
     333                                    label='3.6x3.6RCBC',
     334                                    verbose=False)
     335
     336        culvert.determine_inflow_outflow()
     337       
     338        ( Q, v, d ) = culvert.discharge_routine()
     339       
     340        print test_boyd_non_skew_enquiry_points
     341        print 'Q: ', Q, 'expected_Q: ', expected_Q
     342
     343        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow
     344        assert numpy.allclose(v, expected_v, rtol=1.0e-2) #outflow velocity
     345        assert numpy.allclose(d, expected_d, rtol=1.0e-2) #depth at outlet used to calc v         
     346
    127347
    128348# =========================================================================
Note: See TracChangeset for help on using the changeset viewer.