source: inundation-numpy-branch/debug/mesh_error_reporting/show_mesh_bug.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.4 KB
Line 
1"""Example of shallow water wave equation.
2
3This script sets up LWRU2 benchmark with initial condition stated
4
5See also
6
7http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
8
9Depth at western boundary is d = 13.5 cm
10"""
11
12import sys
13from os import sep
14sys.path.append('..'+sep+'..'+sep)
15
16def prepare_timeboundary(filename):
17    """Converting benchmark 2 time series to NetCDF sww file.
18    """
19
20
21    print 'Preparing time boundary from %s' %filename
22    from Numeric import array
23
24    fid = open(filename)
25
26    #Skip first line
27    line = fid.readline()
28
29    lines = fid.readlines()
30    fid.close()
31
32    N = len(lines)
33
34    T = zeros(N, Float)  #Time
35    Q = zeros(N, Float)  #Values
36
37    for i, line in enumerate(lines):
38        fields = line.split()
39
40        T[i] = float(fields[0])
41        Q[i] = float(fields[1])
42
43    #Create sww file
44    from Scientific.IO.NetCDF import NetCDFFile
45
46    outfile = filename[:-4] + '.sww'
47    print 'Writing to', outfile
48    fid = NetCDFFile(outfile, 'w')
49
50    fid.institution = 'Geoscience Australia'
51    fid.description = 'Input wave for Benchmark 2'
52    fid.starttime = 0.0
53    fid.createDimension('number_of_timesteps', len(T))
54    fid.createVariable('time', Float, ('number_of_timesteps',))
55    fid.variables['time'][:] = T
56
57    fid.createVariable('stage', Float, ('number_of_timesteps',))
58    fid.variables['stage'][:] = Q[:]
59
60    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
61    fid.variables['xmomentum'][:] = 0.0
62
63    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
64    fid.variables['ymomentum'][:] = 0.0
65
66    fid.close()
67
68
69# Module imports
70from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
71            File_boundary, Transmissive_Momentum_Set_Stage_boundary
72from anuga.pyvolution.mesh_factory import rectangular
73from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
74from Numeric import array, zeros, Float, allclose
75import filenames
76
77picklefile = 'domain.pck'
78from cPickle import dump, load
79
80N = 50
81
82#Make mesh
83execfile('create_mesh.py')
84
85
86#Preparing points
87from anuga.pyvolution.data_manager import xya2pts
88xya2pts(filenames.bathymetry_filename, verbose = True,
89        stride = 16,
90        z_func = lambda z: -z)
91
92
93#######################
94# Domain
95
96print 'Creating domain from', filenames.mesh_filename
97domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)
98
99print "Number of triangles = ", len(domain)
100print 'The extent is ', domain.get_extent()
101
102
103domain.check_integrity()
104domain.default_order = 2
105
106print "Number of triangles = ", len(domain)
107domain.store = True    #Store for visualisation purposes
108
109import sys, os
110base = os.path.basename(sys.argv[0])
111domain.filename, _ = os.path.splitext(base)
112
113#LS code to be included in set_quantity
114from anuga.pyvolution import util, least_squares
115points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
116
117#Fit attributes to mesh
118vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
119                                              domain.triangles,
120                                              points,
121                                              attributes['elevation'],
122                                              alpha = 0.01,
123                                              verbose=True)
124
125
126
127
128
Note: See TracBrowser for help on using the repository browser.