Python Domain.get_vertex_coordinates方法代码示例

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import get_vertex_coordinates [as 别名]
    def setUp(self):
        # print "****set up****"
        # Create an sww file

        # Set up an sww that has a geo ref.
        # have it cover an area in Australia.  'gong maybe
        # Don't have many triangles though!

        # Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
        # Zone:   56
        # Easting:  222908.705  Northing: 6233785.284
        # Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
        # Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660

        # geo-ref
        # Zone:   56
        # Easting:  220000  Northing: 6230000

        # have  a big area covered.
        mesh_file = tempfile.mktemp(".tsh")
        points_lat_long = [[-33, 152], [-35, 152], [-35, 150], [-33, 150]]
        spat = Geospatial_data(data_points=points_lat_long, points_are_lats_longs=True)
        points_ab = spat.get_data_points(absolute=True)

        geo = Geo_reference(56, 400000, 6000000)
        spat.set_geo_reference(geo)
        m = Mesh()
        m.add_vertices(spat)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)

        # Create shallow water domain
        domain = Domain(mesh_file)

        os.remove(mesh_file)

        domain.default_order = 2
        # Set some field values
        # domain.set_quantity('stage', 1.0)
        domain.set_quantity("elevation", -0.5)
        domain.set_quantity("friction", 0.03)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary({"exterior": 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.set_quantity("xmomentum", stage * 22.0)
        domain.set_quantity("ymomentum", stage * 55.0)

        domain.distribute_to_vertices_and_edges()

        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

        # sww_file = tempfile.mktemp("")
        self.domain.set_name("tid_P0")
        self.domain.format = "sww"
        self.domain.smooth = True
        self.domain.reduction = mean

        sww = SWW_file(self.domain)
        sww.store_connectivity()
        sww.store_timestep()
        self.domain.time = 2.0
        sww.store_timestep()
        self.sww = sww  # so it can be deleted

        # Create another sww file
        mesh_file = tempfile.mktemp(".tsh")
        points_lat_long = [[-35, 152], [-36, 152], [-36, 150], [-35, 150]]
        spat = Geospatial_data(data_points=points_lat_long, points_are_lats_longs=True)
        points_ab = spat.get_data_points(absolute=True)

        geo = Geo_reference(56, 400000, 6000000)
        spat.set_geo_reference(geo)
        m = Mesh()
        m.add_vertices(spat)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import get_vertex_coordinates [as 别名]
    def test_read_sww(self):
        """
        Save to an sww file and then read back the info.
        Here we store the info "uniquely"
        """

        # ---------------------------------------------------------------------
        # Import necessary modules
        # ---------------------------------------------------------------------
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
        from anuga.shallow_water.shallow_water_domain import Domain
        from anuga import Reflective_boundary
        from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary, Time_boundary

        # ---------------------------------------------------------------------
        # Setup computational domain
        # ---------------------------------------------------------------------
        length = 8.0
        width = 4.0
        dx = dy = 2  # Resolution: Length of subdivisions on both axes

        inc = 0.05  # Elevation increment

        points, vertices, boundary = rectangular_cross(int(length / dx), int(width / dy), len1=length, len2=width)
        domain = Domain(points, vertices, boundary)
        domain.set_name("read_sww_test" + str(domain.processor))  # Output name
        domain.set_quantities_to_be_stored({"elevation": 2, "stage": 2, "xmomentum": 2, "ymomentum": 2, "friction": 1})

        domain.set_store_vertices_uniquely(True)

        # ---------------------------------------------------------------------
        # Setup initial conditions
        # ---------------------------------------------------------------------
        domain.set_quantity("elevation", 0.0)  # Flat bed initially
        domain.set_quantity("friction", 0.01)  # Constant friction
        domain.set_quantity("stage", 0.0)  # Dry initial condition

        # ------------------------------------------------------------------
        # Setup boundary conditions
        # ------------------------------------------------------------------
        Bi = Dirichlet_boundary([0.4, 0, 0])  # Inflow
        Br = Reflective_boundary(domain)  # Solid reflective wall
        Bo = Dirichlet_boundary([-5, 0, 0])  # Outflow

        domain.set_boundary({"left": Bi, "right": Bo, "top": Br, "bottom": Br})

        # -------------------------------------------------------------------
        # Evolve system through time
        # -------------------------------------------------------------------

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

        # Check that quantities have been stored correctly
        source = domain.get_name() + ".sww"

        # x = fid.variables['x'][:]
        # y = fid.variables['y'][:]
        # stage = fid.variables['stage'][:]
        # elevation = fid.variables['elevation'][:]
        # fid.close()

        # assert len(stage.shape) == 2
        # assert len(elevation.shape) == 2

        # M, N = stage.shape

        sww_file = sww.Read_sww(source)

        # print 'last frame number',sww_file.get_last_frame_number()

        assert num.allclose(sww_file.x, domain.get_vertex_coordinates()[:, 0])
        assert num.allclose(sww_file.y, domain.get_vertex_coordinates()[:, 1])

        assert num.allclose(sww_file.time, [0.0, 1.0, 2.0, 3.0, 4.0])

        M = domain.get_number_of_triangles()

        assert num.allclose(num.reshape(num.arange(3 * M), (M, 3)), sww_file.vertices)

        last_frame_number = sww_file.get_last_frame_number()
        assert last_frame_number == 4

        assert num.allclose(sww_file.get_bounds(), [0.0, length, 0.0, width])

        assert "stage" in sww_file.quantities.keys()
        assert "friction" in sww_file.quantities.keys()
        assert "elevation" in sww_file.quantities.keys()
        assert "xmomentum" in sww_file.quantities.keys()
        assert "ymomentum" in sww_file.quantities.keys()

        for qname, q in sww_file.read_quantities(last_frame_number).items():

            # print qname
            # print num.linalg.norm(num.abs((domain.get_quantity(qname).get_values()-q).flatten()), ord=1)

            assert num.allclose(domain.get_quantity(qname).get_values(), q)

        # -----------------------------------------
        # Start the evolution off again at frame 3
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import get_vertex_coordinates [as 别名]
    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],
#.........这里部分代码省略.........