source: anuga_work/development/demos/cornell_room.py @ 5162

Last change on this file since 5162 was 4539, checked in by ole, 18 years ago

Moved files from examples to anuga_work

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.store = True    #Store for visualisation purposes
66domain.format = 'sww'   #Native netcdf visualisation format
67import sys, os
68root, ext = os.path.splitext(sys.argv[0])
69if domain.smooth is True:
70    s = 'smooth'
71else:
72    s = 'nonsmooth'       
73domain.set_name(root + '_' + s)
74
75#Set bed-slope and friction
76manning = 0.0
77
78print 'Field values'
79domain.set_quantity('elevation', mesh_values)
80domain.set_quantity('friction', manning)
81
82
83######################
84# Boundary conditions
85#
86print 'Boundaries'
87Br = Reflective_boundary(domain)
88domain.set_boundary({'exterior': Br})
89
90                 
91
92######################
93#Initial condition
94#
95print 'Initial condition'
96
97#Define water height as a lump in one corner
98def height(x, y):
99    from Numeric import zeros, Float
100   
101    N = len(x)
102    assert N == len(y)   
103
104    xmin = min(x); xmax = max(x)
105    ymin = min(y); ymax = max(y)
106
107    xrange = xmax - xmin
108    yrange = ymax - ymin   
109
110    z = zeros(N, Float)
111    for i in range(N):
112        if x[i] <= xmin + 0.25*xrange and y[i] <= ymin + 0.25*yrange:
113            z[i] = 300
114
115    return z
116
117domain.set_quantity('stage', height)
118
119E = domain.quantities['elevation'].vertex_values
120L = domain.quantities['stage'].vertex_values
121domain.set_quantity('stage', E+L)
122
123#Evolve
124for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
125    domain.write_time()
126   
127
128   
Note: See TracBrowser for help on using the repository browser.