Changeset 6121


Ignore:
Timestamp:
Jan 7, 2009, 4:34:45 PM (15 years ago)
Author:
ole
Message:

Changed culvert polygons to compute enquiry points and verify that they are always outside exchange area.
Worked on volumetric conservation of culverts - not succeeding. Test is disabled.

Location:
anuga_core/source/anuga
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6111 r6121  
    66from anuga.utilities.polygon import plot_polygons
    77
     8from anuga.utilities.numerical_tools import mean
    89from anuga.utilities.numerical_tools import ensure_numeric
     10       
     11from anuga.config import g, epsilon
     12from Numeric import take, sqrt
     13from anuga.config import velocity_protection       
     14
     15       
     16
     17
    918from Numeric import allclose
     19from Numeric import sqrt, sum
     20
     21
    1022
    1123import sys
    … …  
    1628
    1729
    18 # FIXME(Ole): Write in C and reuse this function by similar code in interpolate.py
     30# FIXME(Ole): Write in C and reuse this function by similar code
     31# in interpolate.py
    1932def interpolate_linearly(x, xvec, yvec):
     33
     34    msg = 'Input to function interpolate_linearly could not be converted '
     35    msg += 'to numerical scalar: x = %s' % str(x)
     36    try:
     37        x = float(x)
     38    except:
     39        raise Exception, msg
     40
    2041
    2142    # Check bounds
    2243    if x < xvec[0]:
    23         msg = 'Value provided = %.2f, interpolation minimum = %.2f.' %(x, xvec[0])
     44        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
     45            % (x, xvec[0])
    2446        raise Below_interval, msg
    2547       
    2648    if x > xvec[-1]:
    27         msg = 'Value provided = %.2f, interpolation maximum = %.2f.' %(x, xvec[-1])
     49        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
     50            %(x, xvec[-1])
    2851        raise Above_interval, msg       
    2952       
    … …  
    6184            barrel_velocity = float(fields[2].strip())
    6285           
    63             rating_curve.append( [head_difference, flow_rate, barrel_velocity] )
     86            rating_curve.append([head_difference, flow_rate, barrel_velocity])
    6487       
    6588        if i == 0:
    … …  
    107130                 end_point0=None,
    108131                 end_point1=None,
     132                 enquiry_point0=None,
     133                 enquiry_point1=None,                                           
    109134                 update_interval=None,
    110135                 log_file=False,
    … …  
    112137                 verbose=False):
    113138       
    114         from Numeric import sqrt, sum
    115139
    116140       
    … …  
    155179                                    number_of_barrels=number_of_barrels)
    156180       
     181        # Select enquiry points
     182        if enquiry_point0 is None:
     183            enquiry_point0 = P['enquiry_point0']
     184           
     185        if enquiry_point1 is None:
     186            enquiry_point1 = P['enquiry_point1']           
     187           
    157188        if verbose is True:
    158189            pass
    … …  
    160191            #               P['exchange_polygon0'],
    161192            #               P['exchange_polygon1'],
    162             #               P['enquiry_polygon0'],
    163             #               P['enquiry_polygon1']],
     193            #               [enquiry_point0, 1.005*enquiry_point0],
     194            #               [enquiry_point1, 1.005*enquiry_point1]],                           
    164195            #              figname='culvert_polygon_output')
    165196
    166                          
    167         # Compute the average point for enquiry
    168         enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
    169         enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
    170        
     197           
     198           
    171199        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
     200
    172201        self.enquiry_indices = []                 
    173202        for point in self.enquiry_points:
    … …  
    225254        dq = domain.quantities                                           
    226255        for i, opening in enumerate(self.openings):                           
    227             #elevation = dq['elevation'].get_values(location='centroids',
    228             #                                      interpolation_points=[self.enquiry_points[i]])
    229            
    230256            elevation = dq['elevation'].get_values(location='centroids',
    231257                                                   indices=[self.enquiry_indices[i]])           
    … …  
    258284        self.verbose = verbose
    259285        self.last_update = 0.0 # For use with update_interval       
     286        self.last_time = 0.0               
    260287        self.update_interval = update_interval
    261288       
    … …  
    274301       
    275302    def __call__(self, domain):
    276         from anuga.utilities.numerical_tools import mean
    277        
    278         from anuga.config import g, epsilon
    279         from Numeric import take, sqrt
    280         from anuga.config import velocity_protection       
    281 
    282303
    283304        # Time stuff
    … …  
    288309        if self.update_interval is None:
    289310            update = True
     311            delta_t = domain.timestep # Next timestep has been computed in domain.py
    290312        else:   
    291313            if time - self.last_update > self.update_interval or time == 0.0:
    292314                update = True
    293 
     315            delta_t = self.update_interval
     316           
     317        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
     318        if hasattr(self, 'log_filename'):           
     319            log_to_file(self.log_filename, s)
     320               
    294321                               
    295322        if update is True:
    296323            self.last_update = time
     324       
    297325            dq = domain.quantities
    298326                       
    … …  
    303331                # Compute mean values of selected quantitites in the
    304332                # enquiry area in front of the culvert
    305                 # Stage and velocity comes from enquiry area
    306                 # and elevation from exchange area
    307333               
    308334                stage = dq['stage'].get_values(location='centroids',
    … …  
    346372                        log_to_file(self.log_filename, msg)
    347373                except Above_interval, e:
    348                     Q = self.rating_curve[-1,1]             
     374                    Q = self.rating_curve[-1,1]         
    349375                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
    350376                    msg += str(e)
    … …  
    360386            Q *= self.number_of_barrels
    361387           
    362            
     388
     389            # Adjust Q downwards depending on available water at inlet
     390            stage = self.inlet.get_quantity_values(quantity_name='stage')
     391            elevation = self.inlet.get_quantity_values(quantity_name='elevation')
     392            depth = stage-elevation
     393           
     394           
     395            V = 0
     396            for i, d in enumerate(depth):
     397                V += d * domain.areas[i]
     398           
     399            #Vsimple = mean(depth)*self.inlet.exchange_area # Current volume in exchange area 
     400            #print 'Q', Q, 'dt', delta_t, 'Q*dt', Q*delta_t, 'V', V, 'Vsimple', Vsimple
     401
     402            dt = delta_t           
     403            if Q*dt > V:
     404           
     405                Q_reduced = 0.9*V/dt # Reduce with safety factor
     406               
     407                msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
     408                msg += 'greater than current volume (V) at inlet.\n'
     409                msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
     410               
     411                #print msg
     412               
     413                if self.verbose is True:
     414                    print msg
     415                if hasattr(self, 'log_filename'):                   
     416                    log_to_file(self.log_filename, msg)                                       
     417               
     418                Q = Q_reduced
     419       
    363420            self.inlet.rate = -Q
    364421            self.outlet.rate = Q
    … …  
    373430                fid.close()
    374431
    375 
     432            # Store value of time
     433            self.last_time = time
    376434
    377435           
    … …  
    445503                 culvert_routine=None,
    446504                 number_of_barrels=1,
     505                 enquiry_point0=None,
     506                 enquiry_point1=None,                                           
    447507                 update_interval=None,
    448508                 verbose=False):
    … …  
    525585                                    number_of_barrels=number_of_barrels)
    526586       
     587        # Select enquiry points
     588        if enquiry_point0 is None:
     589            enquiry_point0 = P['enquiry_point0']
     590           
     591        if enquiry_point1 is None:
     592            enquiry_point1 = P['enquiry_point1']           
     593           
    527594        if verbose is True:
    528595            pass
    … …  
    530597            #               P['exchange_polygon0'],
    531598            #               P['exchange_polygon1'],
    532             #               P['enquiry_polygon0'],
    533             #               P['enquiry_polygon1']],
     599            #               [enquiry_point0, 1.005*enquiry_point0],
     600            #               [enquiry_point1, 1.005*enquiry_point1]],
    534601            #              figname='culvert_polygon_output')
    535             #import sys; sys.exit()                           
    536 
    537            
    538         # Compute the average point for enquiry
    539         enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
    540         enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
    541        
     602
     603
    542604        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
     605       
     606       
    543607        self.enquiry_indices = []                 
    544608        for point in self.enquiry_points:
    … …  
    579643        for key in P.keys():
    580644            if key in ['exchange_polygon0',
    581                        'exchange_polygon1',
    582                        'enquiry_polygon0',
    583                        'enquiry_polygon1']:
     645                       'exchange_polygon1']:
    584646                for point in P[key]:
    585647               
    … …  
    606668        self.invert_levels = [invert_level0, invert_level1]               
    607669        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
    608         self.enquiry_polylines = [P['enquiry_polygon0'][:2],
    609                                   P['enquiry_polygon1'][:2]]
     670        #self.enquiry_polylines = [P['enquiry_polygon0'][:2],
     671        #                          P['enquiry_polygon1'][:2]]
    610672        self.vector = P['vector']
    611673        self.length = P['length']; assert self.length > 0.0
    … …  
    657719       
    658720    def __call__(self, domain):
    659         from anuga.utilities.numerical_tools import mean
    660        
    661         from anuga.config import g, epsilon
    662         from Numeric import take, sqrt
    663         from anuga.config import velocity_protection       
    664 
    665721
    666722        log_filename = self.log_filename
    … …  
    669725        time = domain.get_time()
    670726       
    671        
     727        # Short hand
     728        dq = domain.quantities
     729               
     730
    672731        update = False
    673732        if self.update_interval is None:
    674733            update = True
     734            delta_t = domain.timestep # Next timestep has been computed in domain.py
    675735        else:   
    676736            if time - self.last_update > self.update_interval or time == 0.0:
    677737                update = True
    678 
    679         #print 'call', time, time - self.last_update
     738            delta_t = self.update_interval
     739           
     740        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
     741        if hasattr(self, 'log_filename'):           
     742            log_to_file(log_filename, s)
    680743               
    681744                               
    682745        if update is True:
    683             #print 'Updating', time, time - self.last_update
    684746            self.last_update = time
    685        
    686             delta_t = time-self.last_time
    687             s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    688             log_to_file(log_filename, s)
    689        
     747                       
    690748            msg = 'Time did not advance'
    691749            if time > 0.0: assert delta_t > 0.0, msg
    … …  
    695753            openings = self.openings   # There are two Opening [0] and [1]
    696754            for i, opening in enumerate(openings):
    697                 dq = domain.quantities
    698755               
    699756                # Compute mean values of selected quantitites in the
    700757                # exchange area in front of the culvert
    701                 # Stage and velocity comes from enquiry area
    702                 # and elevation from exchange area
    703                
    704                 stage = dq['stage'].get_values(location='centroids',
    705                                                indices=opening.exchange_indices)           
     758                     
     759                stage = opening.get_quantity_values(quantity_name='stage')
    706760                w = mean(stage) # Average stage
    707761
    … …  
    711765                    z = invert_level
    712766                else:
    713                     elevation = dq['elevation'].get_values(location='centroids',
    714                                                            indices=opening.exchange_indices)
     767                    elevation = opening.get_quantity_values(quantity_name='elevation')
    715768                    z = mean(elevation) # Average elevation
    716769
    … …  
    734787                   
    735788                # Average measures of energy in front of this opening
    736                 #polyline = self.enquiry_polylines[i]
    737                 #opening.total_energy = domain.get_energy_through_cross_section(polyline,
    738                 #                                                               kind='total')           
    739789               
    740790                id = [self.enquiry_indices[i]]
    … …  
    742792                                               indices=id)
    743793                elevation = dq['elevation'].get_values(location='centroids',
    744                                                indices=id)                                               
     794                                                       indices=id)                                               
    745795                xmomentum = dq['xmomentum'].get_values(location='centroids',
    746                                                indices=id)                                               
     796                                                       indices=id)                                               
    747797                ymomentum = dq['xmomentum'].get_values(location='centroids',
    748798                                                       indices=id)                                                                                             
    … …  
    840890
    841891            # Update momentum       
    842             delta_t = time - self.last_time
     892
    843893            if delta_t > 0.0:
    844894                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    … …  
    873923                log_to_file(log_filename, s)
    874924           
     925            # Store value of time
     926            self.last_time = time
    875927               
    876928
    … …  
    886938        self.outlet.momentum[1](domain)       
    887939           
    888         # Store value of time #FIXME(Ole): Maybe only every time we update   
    889         self.last_time = time
    890940
    891941

  • anuga_core/source/anuga/culvert_flows/culvert_polygons.py

    r6098 r6121  
    55from math import sqrt
    66from Numeric import array, sum
     7from anuga.utilities.polygon import inside_polygon, polygon_area
    78
    89def create_culvert_polygons(end_point0,
    910                            end_point1,
    1011                            width, height=None,
    11                             enquiry_gap_factor=1.0,
    12                             enquiry_shape_factor=2.0,
     12                            enquiry_gap_factor=0.2,
    1313                            number_of_barrels=1):
    1414    """Create polygons at the end of a culvert inlet and outlet.
    … …  
    2323    Input (optional):       
    2424        height - culvert height, defaults to width making a square culvert
    25         enquiry_gap_factor - sets the distance to the enquiry polygon
    26         enquiry_shape_factor - sets the shape of the enquiry polygon
    27                                (large value widens polygon but reduces height
    28                                to preserve same area as exchange polygon)
     25        enquiry_gap_factor - sets the distance to the enquiry point as fraction of the height
    2926        number_of_barrels - number of identical pipes.
    3027       
    … …  
    3431            'exchange_polygon0' - polygon defining the flow area at end_point0
    3532            'exchange_polygon1' - polygon defining the flow area at end_point1
    36             'enquiry_polygon0' - polygon defining the enquiry field near end_point0
    37             'enquiry_polygon1' - polygon defining the enquiry field near end_point1           
     33            'enquiry_point0' - point beyond exchange_polygon0
     34            'enquiry_point1' - point beyond exchange_polygon1           
     35            'vector'
     36            'length'
     37            'normal'
    3838    """   
    3939
    … …  
    6262   
    6363    # Unit direction vector and normal
    64     vector /= length
    65     normal = array([-dy, dx])/length
     64    vector /= length                 # Unit vector in culvert direction
     65    normal = array([-dy, dx])/length # Normal vector
    6666   
     67    culvert_polygons['vector'] = vector
     68    culvert_polygons['length'] = length
     69    culvert_polygons['normal'] = normal   
    6770
    6871    # Short hands
    6972    w = 0.5*width*normal # Perpendicular vector of 1/2 width
    70     h = height*vector  # Vector of length=height in the
    71                        # direction of the culvert
     73    h = height*vector    # Vector of length=height in the
     74                         # direction of the culvert
     75    gap = (1 + enquiry_gap_factor)*h
     76                         
    7277
    73     # Build exchange polygon 0
     78    # Build exchange polygon and enquiry point for opening 0
    7479    p0 = end_point0 + w
    7580    p1 = end_point0 - w
    … …  
    7782    p3 = p0 - h
    7883    culvert_polygons['exchange_polygon0'] = array([p0,p1,p2,p3])
     84    culvert_polygons['enquiry_point0'] = end_point0 - gap
     85   
    7986
    80     # Build exchange polygon 1
     87    # Build exchange polygon and enquiry point for opening 1
    8188    p0 = end_point1 + w
    8289    p1 = end_point1 - w
    … …  
    8491    p3 = p0 + h
    8592    culvert_polygons['exchange_polygon1'] = array([p0,p1,p2,p3])
     93    culvert_polygons['enquiry_point1'] = end_point1 + gap 
    8694
     95    # Check that enquiry polygons are outside exchange polygons
     96    for key1 in ['exchange_polygon0',
     97                 'exchange_polygon1']:
     98        polygon = culvert_polygons[key1]
     99        area = polygon_area(polygon)
     100       
     101        msg = 'Polygon %s ' %(polygon)
     102        msg += ' has area = %f' % area
     103        assert area > 0.0, msg
    87104
    88 
    89     # Redefine shorthands for enquiry polygons
    90     w = w*enquiry_shape_factor
    91     h = h/enquiry_shape_factor
    92     gap = (enquiry_gap_factor + h)*vector
    93    
    94     # Build enquiry polygon 0
    95     p0 = end_point0 + w - gap
    96     p1 = end_point0 - w - gap
    97     p2 = p1 - h
    98     p3 = p0 - h
    99     culvert_polygons['enquiry_polygon0'] = array([p0,p1,p2,p3])
    100 
    101     # Build enquiry polygon 1
    102     p0 = end_point1 + w + gap
    103     p1 = end_point1 - w + gap
    104     p2 = p1 + h
    105     p3 = p0 + h
    106     culvert_polygons['enquiry_polygon1'] = array([p0,p1,p2,p3])   
     105        for key2 in ['enquiry_point0', 'enquiry_point1']:
     106            point = culvert_polygons[key2]
     107            msg = 'Enquiry point falls inside an enquiry point.'
     108            msg += 'Email Ole.Nielsen@ga.gov.au'
     109            assert not inside_polygon(point, polygon), msg
    107110
    108111    # Return results
    109     culvert_polygons['vector'] = vector
    110     culvert_polygons['length'] = length
    111     culvert_polygons['normal'] = normal   
    112112    return culvert_polygons

  • anuga_core/source/anuga/culvert_flows/culvert_routines.py

    r5862 r6121  
    66
    77"""
     8
     9#NOTE:
     10# Inlet control:  Delta_Et > Es at the inlet
     11# Outlet control: Delta_Et < Es at the inlet
     12# where Et is total energy (w + 0.5*v^2/g) and
     13# Es is the specific energy (h + 0.5*v^2/g)
     14
     15
    816
    917#   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
    … …  
    6169    if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01:
    6270        # Calculate driving energy
     71        # FIXME(Ole): Should this be specific energy?
    6372        E = inlet.total_energy
    6473

  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r6111 r6121  
    143143       
    144144
    145     def test_that_culvert_dry_bed_rating(self):
     145    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
    146146        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
    147147       
    … …  
    243243   
    244244
    245 
    246     def test_that_culvert_dry_bed_boyd(self):
    247         """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
     245   
     246    def Xtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
     247        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
     248       
     249        Test that culvert on a sloping dry bed limits flows when very little water
     250        is present at inlet
     251
     252        This one is using the rating curve variant
     253        """
     254
     255        path = get_pathname_from_package('anuga.culvert_flows')   
     256       
     257        length = 40.
     258        width = 5.
     259
     260        dx = dy = 1           # Resolution: Length of subdivisions on both axes
     261
     262        points, vertices, boundary = rectangular_cross(int(length/dx),
     263                                                       int(width/dy),
     264                                                       len1=length,
     265                                                       len2=width)
     266        domain = Domain(points, vertices, boundary)   
     267        domain.set_name('Test_culvert_shallow') # Output name
     268        domain.set_default_order(2)
     269
     270
     271        #----------------------------------------------------------------------
     272        # Setup initial conditions
     273        #----------------------------------------------------------------------
     274
     275        def topography(x, y):
     276            """Set up a weir
     277           
     278            A culvert will connect either side
     279            """
     280            # General Slope of Topography
     281            z=-x/1000
     282           
     283            N = len(x)
     284            for i in range(N):
     285
     286               # Sloping Embankment Across Channel
     287                if 5.0 < x[i] < 10.1:
     288                    # Cut Out Segment for Culvert FACE               
     289                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
     290                       z[i]=z[i]
     291                    else:
     292                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
     293                if 10.0 < x[i] < 12.1:
     294                   z[i] +=  2.5                    # Flat Crest of Embankment
     295                if 12.0 < x[i] < 14.5:
     296                    # Cut Out Segment for Culvert FACE               
     297                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
     298                       z[i]=z[i]
     299                    else:
     300                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
     301                                   
     302                       
     303            return z
     304
     305
     306        domain.set_quantity('elevation', topography)
     307        domain.set_quantity('friction', 0.01)         # Constant friction
     308        domain.set_quantity('stage',
     309                            expression='elevation + 0.2') # Shallow initial condition
     310                           
     311        # NOTE: Shallow values may cause this test to fail regardless of the
     312        # culvert due to initial adjustments. A good value is 0.2
     313
     314
     315        filename = os.path.join(path, 'example_rating_curve.csv')
     316        culvert = Culvert_flow_rating(domain,
     317                                      culvert_description_filename=filename,
     318                                      end_point0=[9.0, 2.5],
     319                                      end_point1=[13.0, 2.5],
     320                                      verbose=False)
     321
     322        domain.forcing_terms.append(culvert)
     323       
     324
     325        #-----------------------------------------------------------------------
     326        # Setup boundary conditions
     327        #-----------------------------------------------------------------------
     328
     329        # Inflow based on Flow Depth and Approaching Momentum
     330
     331        Br = Reflective_boundary(domain)              # Solid reflective wall
     332        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     333
     334
     335        #-----------------------------------------------------------------------
     336        # Evolve system through time
     337        #-----------------------------------------------------------------------
     338
     339        ref_volume = domain.get_quantity('stage').get_integral()
     340        for t in domain.evolve(yieldstep = 1, finaltime = 25):
     341            print domain.timestepping_statistics()
     342            new_volume = domain.get_quantity('stage').get_integral()
     343           
     344            msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
     345                % (new_volume, ref_volume)
     346            assert allclose(new_volume, ref_volume), msg
     347   
     348   
     349   
     350    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
     351        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
    248352       
    249353        Test that culvert on a sloping dry bed doesn't produce flows
    … …  
    436540#-------------------------------------------------------------
    437541if __name__ == "__main__":
    438     #suite = unittest.makeSuite(Test_Culvert, 'test_predicted_boyd_flow')
    439     suite = unittest.makeSuite(Test_Culvert, 'test')
     542    suite = unittest.makeSuite(Test_Culvert, 'test_that_culvert_rating_limits_flow_in_shallow_inlet_condition')
     543    #suite = unittest.makeSuite(Test_Culvert, 'test')
    440544    runner = unittest.TextTestRunner()
    441545    runner.run(suite)

  • anuga_core/source/anuga/culvert_flows/test_culvert_polygons.py

    r6098 r6121  
    2727
    2828        P = create_culvert_polygons(end_point0,
    29                                 end_point1,
    30                                 width=width,   
    31                                 height=height,
    32                                 number_of_barrels=number_of_barrels)
     29                                    end_point1,
     30                                    width=width,   
     31                                    height=height,
     32                                    number_of_barrels=number_of_barrels)
    3333       
    3434        # Compute area and check that it is greater than 0
    35         for key in ['exchange_polygon0',
    36                     'exchange_polygon1',
    37                     'enquiry_polygon0',
    38                     'enquiry_polygon1']:
    39             polygon = P[key]
     35        for key1 in ['exchange_polygon0',
     36                     'exchange_polygon1']:
     37            polygon = P[key1]
    4038            area = polygon_area(polygon)
    4139           
    … …  
    4442            assert area > 0.0, msg
    4543
     44            for key2 in ['enquiry_point0', 'enquiry_point1']:
     45                point = P[key2]
     46                assert not inside_polygon(point, polygon)
     47       
    4648
    4749    def test_2(self):
    … …  
    5557
    5658        P = create_culvert_polygons(end_point0,
    57                                 end_point1,
    58                                 width=width,   
    59                                 height=height,
    60                                 number_of_barrels=number_of_barrels)
     59                                    end_point1,
     60                                    width=width,   
     61                                    height=height,
     62                                    number_of_barrels=number_of_barrels)
    6163       
    6264        # Compute area and check that it is greater than 0
    63         for key in ['exchange_polygon0',
    64                     'exchange_polygon1',
    65                     'enquiry_polygon0',
    66                     'enquiry_polygon1']:
    67             polygon = P[key]
     65        for key1 in ['exchange_polygon0',
     66                    'exchange_polygon1']:
     67            polygon = P[key1]
    6868            area = polygon_area(polygon)
    6969           
    … …  
    7171            msg += ' has area = %f' % area
    7272            assert area > 0.0, msg
    73            
     73
     74            for key2 in ['enquiry_point0', 'enquiry_point1']:
     75                point = P[key2]
     76                assert not inside_polygon(point, polygon)                       
    7477
    7578   

  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6086 r6121  
    19291929
    19301930
    1931     def get_quantity_values(self):
     1931    def get_quantity_values(self, quantity_name=None):
    19321932        """Return values for specified quantity restricted to opening
    1933         """
    1934        
    1935         q = self.domain.quantities[self.quantity_name]
    1936         return q.get_values(indices=self.exchange_indices)
     1933       
     1934        Optionally a quantity name can be specified if values from another
     1935        quantity is sought
     1936        """
     1937       
     1938        if quantity_name is None:
     1939            quantity_name = self.quantity_name
     1940           
     1941        q = self.domain.quantities[quantity_name]
     1942        return q.get_values(location='centroids',
     1943                            indices=self.exchange_indices)
    19371944   
    19381945
    1939     def set_quantity_values(self, val):
     1946    def set_quantity_values(self, val, quantity_name=None):
    19401947        """Set values for specified quantity restricted to opening
    1941         """
    1942 
     1948       
     1949        Optionally a quantity name can be specified if values from another
     1950        quantity is sought       
     1951        """
     1952
     1953        if quantity_name is None:
     1954            quantity_name = self.quantity_name
     1955                   
    19431956        q = self.domain.quantities[self.quantity_name]               
    1944         q.set_values(val, indices=self.exchange_indices)   
     1957        q.set_values(val,
     1958                     location='centroids',
     1959                     indices=self.exchange_indices)   
    19451960
    19461961
Note: See TracChangeset for help on using the changeset viewer.