source: trunk/anuga_core/documentation/user_manual/demos/buildings.py @ 7838

Last change on this file since 7838 was 7838, checked in by hudson, 14 years ago

Got all demos running with new API.

File size: 3.7 KB
Line 
1"""Simple buildings example using ANUGA
2
3    A simple example using internal boundaries to represent buildings, using new
4    functionality found in ANUGA 1.2.
5   
6    written by James Hudson
7"""
8
9#------------------------------------------------------------------------------
10# Import necessary modules
11#------------------------------------------------------------------------------
12import anuga
13
14def poly_from_box(W, E, N, S):
15    """
16        Helper function to create a counter-clockwise polygon
17        from the given box.
18       
19        N, E, S, W are the extents of the box.
20    """
21    return [[W, N], [W, S], [E, S], [E, N]]
22
23#------------------------------------------------------------------------------
24# Setup computational domain - a long skinny channel with obstructing boxes
25#------------------------------------------------------------------------------
26length = 50
27width = 10
28resolution = 1.0 # make this number smaller to make the simulation more accurate
29
30# Create the "world" as a long, skinny channel
31boundary_poly = poly_from_box(0, length, 0, width)
32
33# Place 3 buildings downstream
34building_polys = [  poly_from_box(10, 15, 2.5, 7.5),        # upstream box
35                    poly_from_box(22.5, 27.5, 1.5, 6.5),    # middle box
36                    poly_from_box(35, 40, 3.5, 8.5)]        # downstream box
37
38# create a domain mesh, with 3 building holes in it
39domain = anuga.create_domain_from_regions(boundary_poly,
40                                            boundary_tags={'left': [0],
41                                                   'bottom': [1],
42                                                   'right': [2],
43                                                   'top': [3]},
44                                            maximum_triangle_area = resolution,
45                                            mesh_filename = 'building.msh',
46                                            interior_holes = building_polys,
47                                            use_cache=True, # to speed it up
48                                            verbose=True)   # log output on
49 
50domain.set_name('buildings')                 # Output name
51
52
53#------------------------------------------------------------------------------
54# Setup initial conditions
55#------------------------------------------------------------------------------
56def topography(x,y):
57    return -x/15                             # gentle linear bed slope
58
59domain.set_quantity('elevation', topography) # Use function for elevation
60domain.set_quantity('friction', 0.03)        # Constant friction
61domain.set_quantity('stage',
62                    expression='elevation')  # Dry initial condition
63
64#------------------------------------------------------------------------------
65# Setup boundary conditions
66#------------------------------------------------------------------------------
67Bi = anuga.Dirichlet_boundary([1.0, 0, 0])   # Inflow
68Br = anuga.Reflective_boundary(domain)       # Solid reflective wall
69Bo = anuga.Dirichlet_boundary([-5, 0, 0])    # Outflow
70
71domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br,
72                        'exterior': Br # this is the internal building boundary
73                    })
74
75import cPickle
76
77#------------------------------------------------------------------------------
78# Evolve system through time
79#------------------------------------------------------------------------------
80#for t in domain.evolve(yieldstep=0.2, finaltime=15.0):
81#    print domain.timestepping_statistics()
82
83#cPickle.dump(domain, open('domain_pickle.txt', 'w'))
84 
85domain = cPickle.load(open('domain_pickle.txt'))
86
87   
88# now turn off the tap
89domain.set_boundary({'left': Bo})
90
91for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
92    print domain.timestepping_statistics()       
Note: See TracBrowser for help on using the repository browser.