Python shallow_water_domain.Domain类代码示例

    def test_unique_vertices_average_loc_unique_vert(self):
        """
        get values based on triangle lists.
        """

        #Create basic mesh
        points, vertices, boundary = rectangular(1, 3)
        #Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.build_tagged_elements_dictionary({'bottom':[0,1],
                                                 'top':[4,5],
                                                 'not_bottom':[2,3,4,5]})

        #Set friction
        domain.set_quantity('friction', add_x_y)
        av_bottom = 2.0/3.0
        add = 60.0
        calc_frict = av_bottom + add
        domain.set_tag_region(Add_value_to_region('bottom', 'friction', add,
                          initial_quantity='friction',
                           location='unique vertices',
                           average=True
                          ))

        #print domain.quantities['friction'].get_values()
        frict_points = domain.quantities['friction'].get_values()
        assert num.allclose(frict_points[0],\
                            [ calc_frict, calc_frict, calc_frict])
        assert num.allclose(frict_points[1],\
                            [ calc_frict, calc_frict, calc_frict])
        assert num.allclose(frict_points[2],\
                            [ calc_frict, 1.0 + 2.0/3.0, calc_frict])
        assert num.allclose(frict_points[3],\
                            [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0])

    def test_unique_verticesII(self):
        """
        get values based on triangle lists.
        """

        #Create basic mesh
        points, vertices, boundary = rectangular(1, 3)

        #Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.build_tagged_elements_dictionary({'bottom':[0,1],
                                                 'top':[4,5],
                                                 'all':[0,1,2,3,4,5]})

        #Set friction
        manning = 0.07
        domain.set_quantity('friction', manning)

        domain.set_tag_region(Add_value_to_region('bottom', 'friction', 1.0,initial_quantity='friction', location = 'unique vertices'))

        #print domain.quantities['friction'].get_values()
        assert num.allclose(domain.quantities['friction'].get_values(),\
                            [[ 1.07,  1.07,  1.07],
                             [ 1.07,  1.07,  1.07],
                             [ 1.07,  0.07,  1.07],
                             [ 0.07,  1.07,  0.07],
                             [ 0.07,  0.07,  0.07],
                             [ 0.07,  0.07,  0.07]])

def load_pts_as_polygon(points_file, minimum_triangle_angle=3.0):
    """
    WARNING: This function is not fully working.

    Function to return a polygon returned from alpha shape, given a points file.

    WARNING: Alpha shape returns multiple polygons, but this function only
             returns one polygon.
    """

    from anuga.pmesh.mesh import importMeshFromFile
    from anuga.shallow_water.shallow_water_domain import Domain

    mesh = importMeshFromFile(points_file)
    mesh.auto_segment()
    mesh.exportASCIIsegmentoutlinefile("outline.tsh")
    mesh2 = importMeshFromFile("outline.tsh")
    mesh2.generate_mesh(maximum_triangle_area=1000000000,
                        minimum_triangle_angle=minimum_triangle_angle,
                        verbose=False)
    mesh2.export_mesh_file('outline_meshed.tsh')
    domain = Domain("outline_meshed.tsh", use_cache = False)
    polygon =  domain.get_boundary_polygon()
    return polygon

    def test_get_maximum_inundation_de0(self):
        """Test that sww information can be converted correctly to maximum
        runup elevation and location (without and with georeferencing)

        This test creates a slope and a runup which is maximal (~11m) at around 10s
        and levels out to the boundary condition (1m) at about 30s.
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        #Setup

        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (100m x 100m)
        points, vertices, boundary = rectangular(20, 5, 100, 50)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2
        domain.set_minimum_storable_height(0.01)

        filename = 'runup_test_3'
        domain.set_name(filename)
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        # FIXME (Ole): Backwards compatibility
        # Look at sww file and see what happens when
        # domain.tight_slope_limiters = 1
        domain.tight_slope_limiters = 0
        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)        
        
        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([1.0,0,0])


        #---------- First run without geo referencing
        
        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
        domain.set_quantity('stage', -6)
        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = 50):
            pass


        # Check maximal runup
        runup = get_maximum_inundation_elevation(swwfile)
        location = get_maximum_inundation_location(swwfile)
        #print 'Runup, location', runup, location
        assert num.allclose(runup, 4.66666666667)
        assert num.allclose(location[0], 46.666668) 
               
        # Check final runup
        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
        #print 'Runup, location:',runup, location

        assert num.allclose(runup, 3.81481488546)
        assert num.allclose(location[0], 51.666668)

        # Check runup restricted to a polygon
        p = [[50,1], [99,1], [99,49], [50,49]]
        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
        location = get_maximum_inundation_location(swwfile, polygon=p)
        #print runup, location

        assert num.allclose(runup, 3.81481488546) 
        assert num.allclose(location[0], 51.6666666)                

        # Check that mimimum_storable_height works
        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
        
        stage = fid.variables['stage_c'][:]
        z = fid.variables['elevation_c'][:]
        xmomentum = fid.variables['xmomentum_c'][:]
        ymomentum = fid.variables['ymomentum_c'][:]
        
        for i in range(stage.shape[0]):
            h = stage[i]-z # depth vector at time step i
            
            # Check every node location
            for j in range(stage.shape[1]):
                # Depth being either exactly zero implies
                # momentum being zero.
                # Or else depth must be greater than or equal to
                # the minimal storable height
                if h[j] == 0.0:
                    assert xmomentum[i,j] == 0.0
                    assert ymomentum[i,j] == 0.0                
                else:
                    assert h[j] >= 0.0
        
        fid.close()

#.........这里部分代码省略.........

    def test_merge_swwfiles(self):
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular, \
                                                                    rectangular_cross
        from anuga.shallow_water.shallow_water_domain import Domain
        from anuga.file.sww import SWW_file
        from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import \
            Dirichlet_boundary

        Bd = Dirichlet_boundary([0.5, 0., 0.])

        # Create shallow water domain
        domain = Domain(*rectangular_cross(2, 2))
        domain.set_name('test1')
        domain.set_quantity('elevation', 2)
        domain.set_quantity('stage', 5)
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
        for t in domain.evolve(yieldstep=0.5, finaltime=1):
            pass
            
        domain = Domain(*rectangular(3, 3))
        domain.set_name('test2')
        domain.set_quantity('elevation', 3)
        domain.set_quantity('stage', 50)
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
        for t in domain.evolve(yieldstep=0.5, finaltime=1):
            pass
                
        outfile = 'test_out.sww'
        _sww_merge(['test1.sww', 'test2.sww'], outfile)
        self.assertTrue(os.access(outfile, os.F_OK))  
        
        # remove temp files
        if not sys.platform == 'win32':		
			os.remove('test1.sww')
			os.remove('test2.sww')
			os.remove(outfile)      

    def test_sww2pts_centroids_de0(self):
        """Test that sww information can be converted correctly to pts data at specified coordinates
        - in this case, the centroids.
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile
        # Used for points that lie outside mesh
        NODATA_value = 1758323

        # Setup
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create shallow water domain
        domain = Domain(*rectangular(2, 2))

        B = Transmissive_boundary(domain)
        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})

        domain.set_name('datatest_de0')

        ptsfile = domain.get_name() + '_elevation.pts'
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.set_quantity('elevation', lambda x,y: -x-y)

        domain.geo_reference = Geo_reference(56,308500,6189000)

        sww = SWW_file(domain)
        sww.store_connectivity()
        sww.store_timestep()

        #self.domain.tight_slope_limiters = 1
        domain.evolve_to_end(finaltime = 0.01)
        sww.store_timestep()

        # Check contents in NetCDF
        fid = NetCDFFile(sww.filename, netcdf_mode_r)

        # Get the variables
        x = fid.variables['x'][:]
        y = fid.variables['y'][:]
        elevation = fid.variables['elevation'][:]
        time = fid.variables['time'][:]
        stage = fid.variables['stage'][:]

        volumes = fid.variables['volumes'][:]


        # Invoke interpolation for vertex points       
        points = num.concatenate( (x[:,num.newaxis],y[:,num.newaxis]), axis=1 )
        points = num.ascontiguousarray(points)
        sww2pts(domain.get_name() + '.sww',
                quantity = 'elevation',
                data_points = points,
                NODATA_value = NODATA_value)
        ref_point_values = elevation
        point_values = Geospatial_data(ptsfile).get_attributes()
        #print 'P', point_values
        #print 'Ref', ref_point_values        
        assert num.allclose(point_values, ref_point_values)        



        # Invoke interpolation for centroids
        points = domain.get_centroid_coordinates()
        #print points
        sww2pts(domain.get_name() + '.sww',
                quantity = 'elevation',
                data_points = points,
                NODATA_value = NODATA_value)
        #ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids

        ref_point_values = [-0.77777777, -0.77777777, -0.99999998, -0.99999998, 
                             -0.99999998, -0.99999998, -1.22222221, -1.22222221]
        point_values = Geospatial_data(ptsfile).get_attributes()
        #print 'P', point_values
        #print 'Ref', ref_point_values        
        assert num.allclose(point_values, ref_point_values)        

        fid.close()

        #Cleanup
        os.remove(sww.filename)
        os.remove(ptsfile)

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
print 'Setting up domain'

length = 200. #x-Dir
width = 200.  #y-dir

dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
#dx = dy = .5           # Resolution: Length of subdivisions on both axes
#dx = dy = .5           # Resolution: Length of subdivisions on both axes
#dx = dy = .1           # Resolution: Length of subdivisions on both axes

points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                                    len1=length, len2=width)
domain = Domain(points, vertices, boundary)   
domain.set_name('Test_WIDE_BRIDGE')                 # Output name
#domain.set_default_order(2)
#omain.H0 = 0.01
#domain.tight_slope_limiters = 1

domain.set_flow_algorithm('2_0')

print 'Size', len(domain)

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

def topography(x, y):
    """Set up a weir

    def test_get_mesh_and_quantities_from_de0_sww_file(self):
        """test_get_mesh_and_quantities_from_sww_file(self):
        """     
        
        # Generate a test sww file with non trivial georeference
        
        import time, os

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (100m x 5m)
        width = 5
        length = 50
        t_end = 10
        points, vertices, boundary = rectangular(length, width, 50, 5)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary,
                        geo_reference = Geo_reference(56,308500,6189000))

        domain.set_name('test_get_mesh_and_quantities_from_sww_file')
        swwfile = domain.get_name() + '.sww'
        domain.set_datadir('.')
        domain.set_flow_algorithm('DE0')

        Br = Reflective_boundary(domain)    # Side walls
        Bd = Dirichlet_boundary([1, 0, 0])  # inflow

        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = t_end):
            pass

        
        # Read it

        # Get mesh and quantities from sww file
        X = get_mesh_and_quantities_from_file(swwfile,
                                              quantities=['elevation',
                                                          'stage',
                                                          'xmomentum',
                                                          'ymomentum'], 
                                              verbose=False)
        mesh, quantities, time = X
        

        
        # Check that mesh has been recovered
        assert num.alltrue(mesh.triangles == domain.get_triangles())
        assert num.allclose(mesh.nodes, domain.get_nodes())

        # Check that time has been recovered
        assert num.allclose(time, range(t_end+1))

        # Check that quantities have been recovered
        # (sww files use single precision)
        z=domain.get_quantity('elevation').get_values(location='unique vertices')
        assert num.allclose(quantities['elevation'], z)

        for q in ['stage', 'xmomentum', 'ymomentum']:
            # Get quantity at last timestep
            q_ref=domain.get_quantity(q).get_values(location='unique vertices')

            #print q,quantities[q]
            q_sww=quantities[q][-1,:]
            
            msg = 'Quantity %s failed to be recovered' %q
            assert num.allclose(q_ref, q_sww, atol=1.0e-2), msg

    def test_get_maximum_inundation_from_sww(self):
        """test_get_maximum_inundation_from_sww(self)

        Test of get_maximum_inundation_elevation()
        and get_maximum_inundation_location().
   
        This is based on test_get_maximum_inundation_3(self) but works with the
        stored results instead of with the internal data structure.

        This test uses the underlying get_maximum_inundation_data for tests
        """

        verbose = False
        from anuga.config import minimum_storable_height
        
        initial_runup_height = -0.4
        final_runup_height = -0.3
        filename = 'runup_test_2'

        #--------------------------------------------------------------
        # Setup computational domain
        #--------------------------------------------------------------
        N = 10
        points, vertices, boundary = rectangular_cross(N, N)
        domain = Domain(points, vertices, boundary)
        domain.set_name(filename)
        domain.set_maximum_allowed_speed(1.0)
        #domain.set_minimum_storable_height(1.0e-5)
        domain.set_store_vertices_uniquely()

        # FIXME: This works better with old limiters so far
        domain.tight_slope_limiters = 0

        #--------------------------------------------------------------
        # Setup initial conditions
        #--------------------------------------------------------------
        def topography(x, y):
            return -x/2                             # linear bed slope

        # Use function for elevation
        domain.set_quantity('elevation', topography)
        domain.set_quantity('friction', 0.)                # Zero friction
        # Constant negative initial stage
        domain.set_quantity('stage', initial_runup_height)

        #--------------------------------------------------------------
        # Setup boundary conditions
        #--------------------------------------------------------------
        Br = Reflective_boundary(domain)                       # Reflective wall
        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow

        # All reflective to begin with (still water)
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

        #--------------------------------------------------------------
        # Test initial inundation height
        #--------------------------------------------------------------
        indices = domain.get_wet_elements()
        z = domain.get_quantity('elevation').\
                get_values(location='centroids', indices=indices)
        assert num.alltrue(z < initial_runup_height)

        q_ref = domain.get_maximum_inundation_elevation(minimum_height=minimum_storable_height)
        # First order accuracy
        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)

        #--------------------------------------------------------------
        # Let triangles adjust
        #--------------------------------------------------------------
        q_max = None
        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
            q = domain.get_maximum_inundation_elevation(minimum_height=minimum_storable_height)

            if verbose:
                domain.write_time()
                print q
                
            if q > q_max:
                q_max = q

        #--------------------------------------------------------------
        # Test inundation height again
        #--------------------------------------------------------------
        #q_ref = domain.get_maximum_inundation_elevation()
        q = get_maximum_inundation_elevation(filename+'.sww')
        msg = 'We got %f, should have been %f' % (q, q_max)
        assert num.allclose(q, q_max, rtol=2.0/N), msg

        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg

        # Test error condition if time interval is out
        try:
            q = get_maximum_inundation_elevation(filename+'.sww',
                                                 time_interval=[2.0, 3.0])
        except ValueError:
            pass
        else:
            msg = 'should have caught wrong time interval'
            raise Exception, msg
#.........这里部分代码省略.........

    def test_region_tags(self):
        """get values based on triangle lists."""

        #Create basic mesh
        points, vertices, boundary = rectangular(1, 3)

        #Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.build_tagged_elements_dictionary({'bottom': [0,1],
                                                 'top': [4,5],
                                                 'all': [0,1,2,3,4,5]})

        #Set friction
        manning = 0.07
        domain.set_quantity('friction', manning)

        a = Set_tag_region('bottom', 'friction', 0.09)
        b = Set_tag_region('top', 'friction', 1.0)
        domain.set_tag_region([a, b])

        expected = [[ 0.09,  0.09,  0.09],
                    [ 0.09,  0.09,  0.09],
                    [ 0.07,  0.07,  0.07],
                    [ 0.07,  0.07,  0.07],
                    [ 1.0,   1.0,   1.0],
                    [ 1.0,   1.0,   1.0]]
        msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s"
               % (str(domain.quantities['friction'].get_values()),
                  str(expected)))
        assert num.allclose(domain.quantities['friction'].get_values(),
                            expected), msg

        #c = Add_Value_To_region('all', 'friction', 10.0)
        domain.set_tag_region(Add_value_to_region('all', 'friction', 10.0))
        #print domain.quantities['friction'].get_values()
        assert num.allclose(domain.quantities['friction'].get_values(),
                            [[ 10.09, 10.09, 10.09],
                             [ 10.09, 10.09, 10.09],
                             [ 10.07, 10.07, 10.07],
                             [ 10.07, 10.07, 10.07],
                             [ 11.0,  11.0,  11.0],
                             [ 11.0,  11.0,  11.0]])

        # trying a function
        domain.set_tag_region(Set_tag_region('top', 'friction', add_x_y))
        #print domain.quantities['friction'].get_values()
        assert num.allclose(domain.quantities['friction'].get_values(),
                            [[ 10.09, 10.09, 10.09],
                             [ 10.09, 10.09, 10.09],
                             [ 10.07, 10.07, 10.07],
                             [ 10.07, 10.07, 10.07],
                             [ 5./3,  2.0,  2./3],
                             [ 1.0,  2./3,  2.0]])

        domain.set_quantity('elevation', 10.0)
        domain.set_quantity('stage', 10.0)
        domain.set_tag_region(Add_value_to_region('top', 'stage', 1.0,initial_quantity='elevation'))
        #print domain.quantities['stage'].get_values()
        assert num.allclose(domain.quantities['stage'].get_values(),
                            [[ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 11.0,  11.0,  11.0],
                             [ 11.0,  11.0,  11.0]])

        
        domain.set_quantity('elevation', 10.0)
        domain.set_quantity('stage', give_me_23)
        #this works as well, (is cleaner, but doesn't work for regions)
        #domain.set_quantity('stage',
        #                    domain.quantities['stage'].vertex_values+ \
        #                    domain.quantities['elevation'].vertex_values)
        domain.set_tag_region(Add_quantities('top', 'elevation','stage'))
        #print domain.quantities['stage'].get_values()
        assert num.allclose(domain.quantities['elevation'].get_values(),
                            [[ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 10., 10., 10.],
                             [ 33.,  33.0,  33.],
                             [ 33.0,  33.,  33.]])

    def _create_domain(self,d_length,
                            d_width,
                            dx,
                            dy,
                            elevation_0,
                            elevation_1,
                            stage_0,
                            stage_1):
        
        points, vertices, boundary = rectangular_cross(int(d_length/dx), int(d_width/dy),
                                                        len1=d_length, len2=d_width)
        domain = Domain(points, vertices, boundary)   
        domain.set_name('Test_Outlet_Inlet')                 # Output name
        domain.set_store()
        domain.set_default_order(2)
        domain.H0 = 0.01
        domain.tight_slope_limiters = 1

        #print 'Size', len(domain)

        #------------------------------------------------------------------------------
        # Setup initial conditions
        #------------------------------------------------------------------------------

        def elevation(x, y):
            """Set up a elevation
            """
            
            z = numpy.zeros(x.shape,dtype='d')
            z[:] = elevation_0
            
            numpy.putmask(z, x > d_length/2, elevation_1)
    
            return z
            
        def stage(x,y):
            """Set up stage
            """
            z = numpy.zeros(x.shape,dtype='d')
            z[:] = stage_0
            
            numpy.putmask(z, x > d_length/2, stage_1)

            return z
            
        #print 'Setting Quantities....'
        domain.set_quantity('elevation', elevation)  # Use function for elevation
        domain.set_quantity('stage',  stage)   # Use function for elevation

        Br = anuga.Reflective_boundary(domain)
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
        
        return domain

    def test_earthquake_tsunami(self):
        from os import sep, getenv
        import sys
        from anuga.abstract_2d_finite_volumes.mesh_factory \
                        import rectangular_cross
        from anuga.abstract_2d_finite_volumes.quantity import Quantity
        from anuga.utilities.system_tools import get_pathname_from_package
        """
        Pick the test you want to do; T= 0 test a point source,
        T= 1  test single rectangular source, T= 2 test multiple
        rectangular sources
        """

        # Get path where this test is run
        path= get_pathname_from_package('anuga.tsunami_source')
        
        # Choose what test to proceed
        T=1

        if T==0:
            # Fortran output file
            filename = path+sep+'tests'+sep+'data'+sep+'fullokada_SP.txt'
            
            # Initial condition of earthquake for multiple source
            x0 = 7000.0
            y0 = 10000.0
            length = 0
            width =0
            strike = 0.0
            depth = 15.0
            slip = 10.0
            dip =15.0
            rake =90.0
            ns=1
            NSMAX=1
        elif T==1:
            # Fortran output file
            filename = path+sep+'tests'+sep+'data'+sep+'fullokada_SS.txt'
            
            # Initial condition of earthquake for multiple source
            x0 = 7000.0
            y0 = 10000.0
            length = 10.0
            width =6.0
            strike = 0.0
            depth = 15.0
            slip = 10.0
            dip =15.0
            rake =90.0
            ns=1
            NSMAX=1
            
        elif T==2:

            # Fortran output file
            filename = path+sep+'tests'+sep+'data'+sep+'fullokada_MS.txt'
            
            # Initial condition of earthquake for multiple source
            x0 = [7000.0,10000.0]
            y0 = [10000.0,7000.0]
            length = [10.0,10.0]
            width =[6.0,6.0]
            strike = [0.0,0.0]
            depth = [15.0,15.0]
            slip = [10.0,10.0]
            dip = [15.0,15.0]
            rake = [90.0,90.0]
            ns=2
            NSMAX=2



        # Get output file from original okada fortran script.
        # Vertical displacement is listed under tmp.
        polyline_file=open(filename,'r')
        lines=polyline_file.readlines()
        polyline_file.close()
        tmp=[]
        stage=[]
        for line in lines [0:]:
            field = line.split('    ')
            z=float(field[2])
            tmp.append(z)

         
        # Create domain 
        dx = dy = 4000
        l=20000
        w=20000
        
        # Create topography
        def topography(x,y):
            el=-1000
            return el
        
        points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
                                               len1=l, len2=w)
        domain = Domain(points, vertices, boundary)   
        domain.set_name('test')
        domain.set_quantity('elevation',topography)
#.........这里部分代码省略.........

    def test_get_flow_through_cross_section_stored_uniquely(self):
        """test_get_flow_through_cross_section_stored_uniquely(self):

        Test that the total flow through a cross section can be
        correctly obtained from an sww file.
        
        This test creates a flat bed with a known flow through it and tests
        that the function correctly returns the expected flow.

        The specifics are
        u = 2 m/s
        h = 1 m
        w = 3 m (width of channel)

        q = u*h*w = 6 m^3/s
       
        
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (20m x 3m)
        width = 3
        length = 20
        t_end = 3
        points, vertices, boundary = rectangular(length, width,
                                                 length, width)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2
        domain.set_minimum_storable_height(0.01)

        domain.set_name('flowtest_uniquely')
        swwfile = domain.get_name() + '.sww'

        domain.set_store_vertices_uniquely()
        
        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        h = 1.0
        u = 2.0
        uh = u*h

        Br = Reflective_boundary(domain)     # Side walls
        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet: 


        
        domain.set_quantity('elevation', 0.0)
        domain.set_quantity('stage', h)
        domain.set_quantity('xmomentum', uh)
        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = t_end):
            pass

        # Check that momentum is as it should be in the interior

        I = [[0, width/2.],
             [length/2., width/2.],
             [length, width/2.]]
        
        f = file_function(swwfile,
                          quantities=['stage', 'xmomentum', 'ymomentum'],
                          interpolation_points=I,
                          verbose=False)
        for t in range(t_end+1):
            for i in range(3):
                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
            

        # Check flows through the middle
        for i in range(5):
            x = length/2. + i*0.23674563 # Arbitrary
            cross_section = [[x, 0], [x, width]]
            time, Q = get_flow_through_cross_section(swwfile,
                                                     cross_section,
                                                     verbose=False)

            assert num.allclose(Q, uh*width)


       
        # Try the same with partial lines
        x = length/2.
        for i in range(5):
            start_point = [length/2., i*width/5.]
            #print start_point
                            
            cross_section = [start_point, [length/2., width]]
            time, Q = get_flow_through_cross_section(swwfile,
                                                     cross_section,
                                                     verbose=False)
#.........这里部分代码省略.........

    def test_get_energy_through_cross_section(self):
        """test_get_energy_through_cross_section(self):

        Test that the specific and total energy through a cross section can be
        correctly obtained from an sww file.
        
        This test creates a flat bed with a known flow through it and tests
        that the function correctly returns the expected energies.

        The specifics are
        u = 2 m/s
        h = 1 m
        w = 3 m (width of channel)

        q = u*h*w = 6 m^3/s
        Es = h + 0.5*v*v/g  # Specific energy head [m]
        Et = w + 0.5*v*v/g  # Total energy head [m]        


        This test uses georeferencing
        
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (20m x 3m)
        width = 3
        length = 20
        t_end = 1
        points, vertices, boundary = rectangular(length, width,
                                                 length, width)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary,
                        geo_reference = Geo_reference(56,308500,6189000))

        domain.default_order = 2
        domain.set_minimum_storable_height(0.01)

        domain.set_name('flowtest')
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        e = -1.0
        w = 1.0
        h = w-e
        u = 2.0
        uh = u*h

        Br = Reflective_boundary(domain)     # Side walls
        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 

        
        domain.set_quantity('elevation', e)
        domain.set_quantity('stage', w)
        domain.set_quantity('xmomentum', uh)
        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = t_end):
            pass

        # Check that momentum is as it should be in the interior

        I = [[0, width/2.],
             [length/2., width/2.],
             [length, width/2.]]
        
        I = domain.geo_reference.get_absolute(I)
        f = file_function(swwfile,
                          quantities=['stage', 'xmomentum', 'ymomentum'],
                          interpolation_points=I,
                          verbose=False)

        for t in range(t_end+1):
            for i in range(3):
                #print i, t, f(t, i)
                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
            

        # Check energies through the middle
        for i in range(5):
            x = length/2. + i*0.23674563 # Arbitrary
            cross_section = [[x, 0], [x, width]]

            cross_section = domain.geo_reference.get_absolute(cross_section)            
            
            time, Es = get_energy_through_cross_section(swwfile,
                                                       cross_section,
                                                       kind='specific',
                                                       verbose=False)
            assert num.allclose(Es, h + 0.5*u*u/g)
            
            time, Et = get_energy_through_cross_section(swwfile,
#.........这里部分代码省略.........

    def setUp(self):
        import time
        
        self.verbose = Test_File_Conversion.verbose
        # Create basic mesh
        points, vertices, boundary = rectangular(2, 2)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2

        # Set some field values
        domain.set_quantity('elevation', lambda x,y: -x)
        domain.set_quantity('friction', 0.03)


        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})


        ######################
        #Initial condition - with jumps
        bed = domain.quantities['elevation'].vertex_values
        stage = num.zeros(bed.shape, num.float)

        h = 0.3
        for i in range(stage.shape[0]):
            if i % 2 == 0:
                stage[i,:] = bed[i,:] + h
            else:
                stage[i,:] = bed[i,:]

        domain.set_quantity('stage', stage)


        domain.distribute_to_vertices_and_edges()               
        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)


        self.domain = domain

        C = domain.get_vertex_coordinates()
        self.X = C[:,0:6:2].copy()
        self.Y = C[:,1:6:2].copy()

        self.F = bed

        #Write A testfile (not realistic. Values aren't realistic)
        self.test_MOST_file = 'most_small'

        longitudes = [150.66667, 150.83334, 151., 151.16667]
        latitudes = [-34.5, -34.33333, -34.16667, -34]

        long_name = 'LON'
        lat_name = 'LAT'

        nx = 4
        ny = 4
        six = 6


        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)

            fid.createDimension(long_name,nx)
            fid.createVariable(long_name,netcdf_float,(long_name,))
            fid.variables[long_name].point_spacing='uneven'
            fid.variables[long_name].units='degrees_east'
            fid.variables[long_name][:] = longitudes

            fid.createDimension(lat_name,ny)
            fid.createVariable(lat_name,netcdf_float,(lat_name,))
            fid.variables[lat_name].point_spacing='uneven'
            fid.variables[lat_name].units='degrees_north'
            fid.variables[lat_name][:] = latitudes

            fid.createDimension('TIME',six)
            fid.createVariable('TIME',netcdf_float,('TIME',))
            fid.variables['TIME'].point_spacing='uneven'
            fid.variables['TIME'].units='seconds'
            fid.variables['TIME'][:] = [0.0, 0.1, 0.6, 1.1, 1.6, 2.1]


            name = ext[1:3].upper()
            if name == 'E.': name = 'ELEVATION'
            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
            fid.variables[name].units='CENTIMETERS'
            fid.variables[name].missing_value=-1.e+034

            fid.variables[name][:] = [[[0.3400644, 0, -46.63519, -6.50198],
                                              [-0.1214216, 0, 0, 0],
                                              [0, 0, 0, 0],
                                              [0, 0, 0, 0]],
                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
                                              [0, 2.291054e-005, 2.291054e-005, 0]],
                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
#.........这里部分代码省略.........