source: anuga_core/source/anuga/pyvolution/cornell_room.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 3.1 KB
Line 
1"""Example of shallow water wave equation.
2
3This is called Netherlands because it shows a dam with a gap in it and
4stylised housed behind it and below the water surface.
5
6"""
7
8######################
9# Module imports
10#
11from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
12     Transmissive_boundary, Time_boundary, Constant_height
13
14from mesh_factory import from_polyfile, rectangular
15from Numeric import array
16from math import sqrt
17#from least_squares import Interpolation
18from anuga.fit_interpolate.interpolate import Interpolate
19   
20
21print 'Creating domain'
22#data_points, _, data_values = from_polyfile('cornell_room_medres')
23#points, triangles, values = from_polyfile('hires2')
24data_points, _, data_values = from_polyfile('hires2')
25
26
27#Regrid onto numerically stable mesh
28#
29#Compute regular mesh based on resolution and extent of data
30data_points = array(data_points)
31pmax = max(data_points)
32pmin = min(data_points)
33
34M = len(data_points)
35
36N = int(0.8*sqrt(M))
37
38#print N
39
40mesh_points, vertices, boundary = rectangular(N, N, 
41                                              len1=pmax[0]-pmin[0],
42                                              len2=pmax[1]-pmin[1],
43                                              origin = pmin)
44                                             
45
46#Compute smooth surface on new mesh based on values from old (regrid)
47print 'Interp'
48interp = Interpolate(mesh_points, vertices, alpha=10)
49mesh_values = interp.fit( data_points, data_values) # this has not been tested
50print 'Len mesh values', len(mesh_values)
51print 'Len mesh points', len(mesh_points)
52
53   
54#Create shallow water domain
55print 'Creating domain'
56domain = Domain(mesh_points, vertices) #, boundary)
57
58domain.check_integrity()
59domain.default_order = 2
60domain.smooth = True
61domain.reduction = min  #Looks a lot better on top of steep slopes
62
63print "Number of triangles = ", len(domain)
64
65domain.visualise = False
66domain.checkpoint = False
67domain.store = True    #Store for visualisation purposes
68domain.format = 'sww'   #Native netcdf visualisation format
69import sys, os
70root, ext = os.path.splitext(sys.argv[0])
71if domain.smooth is True:
72    s = 'smooth'
73else:
74    s = 'nonsmooth'       
75domain.filename = root + '_' + s
76
77#Set bed-slope and friction
78manning = 0.0
79
80print 'Field values'
81domain.set_quantity('elevation', mesh_values)
82domain.set_quantity('friction', manning)
83
84
85######################
86# Boundary conditions
87#
88print 'Boundaries'
89Br = Reflective_boundary(domain)
90domain.set_boundary({'exterior': Br})
91
92                 
93
94######################
95#Initial condition
96#
97print 'Initial condition'
98
99#Define water height as a lump in one corner
100def height(x, y):
101    from Numeric import zeros, Float
102   
103    N = len(x)
104    assert N == len(y)   
105
106    xmin = min(x); xmax = max(x)
107    ymin = min(y); ymax = max(y)
108
109    xrange = xmax - xmin
110    yrange = ymax - ymin   
111
112    z = zeros(N, Float)
113    for i in range(N):
114        if x[i] <= xmin + 0.25*xrange and y[i] <= ymin + 0.25*yrange:
115            z[i] = 300
116
117    return z
118
119domain.set_quantity('stage', height)
120
121E = domain.quantities['elevation'].vertex_values
122L = domain.quantities['stage'].vertex_values
123domain.set_quantity('stage', E+L)
124
125#Evolve
126for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
127    domain.write_time()
128   
129
130   
Note: See TracBrowser for help on using the repository browser.