Index: inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 276)
+++ inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 276)
@@ -0,0 +1,790 @@
+"""Functions to store and retrieve data for the Shallow Water Wave equation.
+There are two kinds of data
+
+  1: Constant data: Vertex coordinates and field values
+  2: Variable data: Conserved quantities
+
+All data is assumed to reside at vertex locations.  
+
+"""
+
+from Numeric import concatenate        
+
+
+def make_filename(s):
+    """Transform argument string into a suitable filename
+    """
+
+    s = s.strip()
+    s = s.replace(' ', '_')
+    #s = s.replace('(', '{')
+    #s = s.replace(')', '}')
+    s = s.replace('(', '')
+    s = s.replace(')', '')
+    s = s.replace('__', '_')                
+        
+    return s
+
+
+def check_dir(path, verbose=None):
+    """Check that specified path exists.
+    If path does not exist it will be created if possible    
+
+    USAGE:
+       checkdir(path, verbose):
+
+    ARGUMENTS:
+        path -- Directory
+        verbose -- Flag verbose output (default: None)
+
+    RETURN VALUE:
+        Verified path including trailing separator
+        
+    """
+
+    import os, sys
+    import os.path
+
+    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
+        unix = 0
+    else:
+        unix = 1
+
+
+    if path[-1] != os.sep: 
+        path = path + os.sep  # Add separator for directories
+
+    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
+    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
+        try:
+            exitcode=os.mkdir(path)
+
+            # Change access rights if possible
+            #
+            if unix:
+                exitcode=os.system('chmod 775 '+path)
+            else:
+                pass  # FIXME: What about acces rights under Windows?
+            
+            if verbose: print 'MESSAGE: Directory', path, 'created.'
+          
+        except:
+            print 'WARNING: Directory', path, 'could not be created.'
+            if unix:
+                path = '/tmp/'
+            else:
+                path = 'C:'
+                
+            print 'Using directory %s instead' %path
+
+    return(path)
+
+    
+
+def del_dir(path):
+    """Recursively delete directory path and all its contents
+    """
+
+    import os
+
+    if os.path.isdir(path):
+        for file in os.listdir(path):
+            X = os.path.join(path, file)
+
+
+            if os.path.isdir(X) and not os.path.islink(X):
+                del_dir(X) 
+            else:
+                try:
+                    os.remove(X)
+                except:
+                    print "Could not remove file %s" %X
+
+        os.rmdir(path)
+
+
+	
+def create_filename(filename, format, size, time=None):
+
+    import os
+    from config import data_dir
+
+           
+    FN = check_dir(data_dir) + filename + '_size%d' %size
+	
+    if time is not None:
+        FN += '_time%.2f' %time
+        
+    FN += '.' + format
+        
+    return FN
+    
+
+def get_files(filename, format, size):
+    """Get all file (names) with gven name, size and format
+    """
+
+    import glob
+
+    import os
+    from config import data_dir
+
+           
+    dir = check_dir(data_dir)
+
+    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
+    return glob.glob(pattern)
+
+#Native checkpoint format.
+#Information needed to recreate a state is preserved
+def cpt_variable_writer(filename, t, v0, v1, v2):
+    """Store all conserved quantities to file
+    """
+
+    M, N = v0.shape
+	
+    FN = create_filename(filename, 'cpt', M, t)
+    #print 'Writing to %s' %FN
+    
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(N):        
+            fid.write('%.16e ' %v0[i,j])            
+        for j in range(N):
+            fid.write('%.16e ' %v1[i,j])                        
+        for j in range(N):
+            fid.write('%.16e ' %v2[i,j])                        
+            
+        fid.write('\n')
+    fid.close()
+
+
+def cpt_variable_reader(filename, t, v0, v1, v2):
+    """Store all conserved quantities to file
+    """
+
+    M, N = v0.shape
+	
+    FN = create_filename(filename, 'cpt', M, t)
+    #print 'Reading from %s' %FN
+    
+    fid = open(FN)
+
+
+    for i in range(M):
+        values = fid.readline().split() #Get one line
+
+        for j in range(N):
+            v0[i,j] = float(values[j])
+            v1[i,j] = float(values[3+j])
+            v2[i,j] = float(values[6+j])
+    
+    fid.close()
+
+def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
+    """Writes x,y,z,z,z coordinates of triangles constituting the bed
+    elevation.
+    Not in use pt
+    """
+
+    M, N = v0.shape
+
+    print X0
+    import sys; sys.exit() 
+    FN = create_filename(filename, 'cpt', M)
+    print 'Writing to %s' %FN
+    
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(2):        
+            fid.write('%.16e ' %X0[i,j])   #x, y
+        for j in range(N):                    
+            fid.write('%.16e ' %v0[i,j])       #z,z,z,
+            
+        for j in range(2):                    
+            fid.write('%.16e ' %X1[i,j])   #x, y
+        for j in range(N):                                
+            fid.write('%.16e ' %v1[i,j])  
+        
+        for j in range(2):                                
+            fid.write('%.16e ' %X2[i,j])   #x, y
+        for j in range(N):                                            
+            fid.write('%.16e ' %v2[i,j]) 
+            
+        fid.write('\n')
+    fid.close()
+
+
+
+#Function for storing out to e.g. visualisation
+def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
+    """Writes x,y,z coordinates of triangles constituting the bed elevation.
+    """
+
+    M, N = v0.shape
+	
+    FN = create_filename(filename, 'dat', M)
+    #print 'Writing to %s' %FN
+    
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(2):        
+            fid.write('%f ' %X0[i,j])   #x, y
+        fid.write('%f ' %v0[i,0])       #z
+            
+        for j in range(2):                    
+            fid.write('%f ' %X1[i,j])   #x, y
+        fid.write('%f ' %v1[i,0])       #z
+        
+        for j in range(2):                                
+            fid.write('%f ' %X2[i,j])   #x, y
+        fid.write('%f ' %v2[i,0])       #z
+            
+        fid.write('\n')
+    fid.close()
+
+
+
+def dat_variable_writer(filename, t, v0, v1, v2):
+    """Store water height to file
+    """
+
+    M, N = v0.shape
+	
+    FN = create_filename(filename, 'dat', M, t)
+    #print 'Writing to %s' %FN
+    
+    fid = open(FN, 'w')
+    for i in range(M):
+        fid.write('%.4f ' %v0[i,0])
+        fid.write('%.4f ' %v1[i,0])
+        fid.write('%.4f ' %v2[i,0])        
+
+        fid.write('\n')
+    fid.close()
+
+
+#Function for storing out to e.g. visualisation
+class Data_format:
+    """Generic interface to data formats
+    """
+
+        
+    def __init__(self, domain, extension, mode = 'w'):
+        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
+                                        '''   'w' (write)'''+\
+                                        '''   'r' (read)''' +\
+                                        '''   'a' (append)'''    
+
+        #Create filename
+        #FIXME: This may change if separate timesteps are stored in separate files.
+        self.filename = create_filename(domain.get_name(), extension, len(domain))
+        self.timestep = 0
+        self.number_of_volumes = len(domain)
+        self.domain = domain
+    
+    
+
+    def store_all(self):
+        """Writes x,y,z coordinates of triangles constituting
+        the bed elevation.
+        """
+
+        msg = 'Method should be overridden in concrete class'
+        raise msg
+
+    
+    def store_timestep(self, t, V0, V1, V2):
+        """Store time, water heights (and momentums) to file
+        """
+
+        msg = 'Method should be overridden in concrete class'
+        raise msg        
+    
+    #FIXME: Should we have a general set_precision function?
+    
+
+#Function for storing out to e.g. visualisation
+class Data_format_sww(Data_format):
+    """Interface to native NetCDF format (.sww)
+    """
+
+        
+    def __init__(self, domain, mode = 'w'):
+        from Scientific.IO.NetCDF import NetCDFFile
+        from Numeric import Int, Float, Float32 
+
+        self.precision = Float32 #Use single precision
+        
+        Data_format.__init__(self, domain, 'sww', mode)
+
+
+        # NetCDF file definition
+        fid = NetCDFFile(self.filename, mode)
+        
+        if mode == 'w':
+            
+            #Create new file
+            fid.institution = 'Geoscience Australia'
+            if domain.smooth:
+                fid.smoothing = 'Yes'
+            else:    
+                fid.smoothing = 'No'
+
+            fid.order = domain.default_order    
+
+            # dimension definitions
+            fid.createDimension('number_of_volumes', self.number_of_volumes)  
+            fid.createDimension('number_of_vertices', 3)
+
+            if domain.smooth is True:
+                fid.createDimension('number_of_points', len(domain.vertexlist))
+            else:    
+                fid.createDimension('number_of_points', 3*self.number_of_volumes)
+                
+            fid.createDimension('number_of_timesteps', None) #extensible
+
+            # variable definitions
+            fid.createVariable('x', self.precision, ('number_of_points',))
+            fid.createVariable('y', self.precision, ('number_of_points',))
+            fid.createVariable('z', self.precision, ('number_of_points',))
+        
+            fid.createVariable('volumes', Int, ('number_of_volumes',
+                                                'number_of_vertices'))
+
+            fid.createVariable('time', self.precision,
+                               ('number_of_timesteps',))
+        
+            fid.createVariable('stage', self.precision,
+                               ('number_of_timesteps',
+                                'number_of_points'))
+
+
+        #Close
+        fid.close()                
+
+    def store_all(self):
+        """Specialisation of store all for net CDF format
+
+        Writes x,y,z coordinates of triangles constituting
+        the bed elevation.
+        """
+
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        from Numeric import concatenate
+
+        domain = self.domain
+
+        #Get NetCDF
+        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
+        
+        # Get the variables
+        x = fid.variables['x']
+        y = fid.variables['y']
+        z = fid.variables['z']
+        
+        volumes = fid.variables['volumes']
+
+        # Get X, Y and bed elevation (index=0)
+        X,Y,A,V = domain.get_vertex_values(xy=True,
+                                           value_array='field_values',
+                                           indices = (0,),
+                                           precision = self.precision)
+
+
+        x[:] = X
+        y[:] = Y
+        z[:] = A[:,0]
+
+        volumes[:] = V
+        
+        #Close
+        fid.close()
+
+
+    def store_timestep(self):
+        """Store time, water heights (and momentums) to file
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+	from time import sleep
+		
+        #Get NetCDF
+	retries = 0
+	file_open = False
+	while not file_open and retries < 10:
+  	    try:	
+                fid = NetCDFFile(self.filename, 'a') #Open existing file
+	    except IOError:
+	        #This could happen if someone was reading the file.
+		#In that case, wait a while and try again
+		msg = 'Warning (store_timestep): File %s could not be opened'\
+		%self.filename
+		msg += ' - trying again'
+		print msg
+		retries += 1
+		sleep(1)
+	    else:
+	       file_open = True
+	       
+	if not file_open:
+	    msg = 'File %s could not be opened for append' %self.filename
+	    raise msg
+	    	
+	    
+        domain = self.domain
+        
+        # Get the variables
+        time = fid.variables['time']
+        stage = fid.variables['stage']        
+        i = len(time)
+
+        #Store stage
+        time[i] = self.domain.time
+
+        # Get stage (index = 0)
+        A,V = domain.get_vertex_values(xy=False,
+                                       value_array='conserved_quantities',
+                                       indices = (0,),
+                                       precision = self.precision)
+
+        
+        stage[i,:] = A[:,0]  #Only store stage
+
+        #Flush and close
+        fid.sync()
+        fid.close()
+
+
+#Function for storing out to e.g. visualisation
+class Data_format_xya(Data_format):
+    """Generic interface to data formats
+    """
+
+
+    def __init__(self, domain, mode = 'w'):
+        from Scientific.IO.NetCDF import NetCDFFile
+        from Numeric import Int, Float, Float32 
+
+        self.precision = Float32 #Use single precision
+        
+        Data_format.__init__(self, domain, 'xya', mode)
+        
+    
+    
+
+    def store_all(self):
+        """Specialisation of store all for xya format
+
+        Writes x,y,z coordinates of triangles constituting
+        the bed elevation.
+        """
+
+        from Numeric import concatenate
+
+        domain = self.domain
+        
+        fd = open(self.filename, 'w')
+
+
+        if domain.smooth is True:
+            number_of_points =  len(domain.vertexlist)
+        else:    
+            number_of_points = 3*self.number_of_volumes
+
+        numVertAttrib = 3 #Three attributes is what is assumed by the xya format
+
+        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
+                 " # <vertex #> <x> <y> [attributes]" + "\n")
+
+
+        # Get X, Y, bed elevation and friction (index=0,1)
+        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
+                                           indices = (0,1), precision = self.precision)
+
+        bed_eles = A[:,0]
+        fricts = A[:,1]
+
+        # Get stage (index=0)
+        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
+                                       indices = (0,), precision = self.precision)
+
+        stages = B[:,0]
+
+        #<vertex #> <x> <y> [attributes]
+        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
+                                               stages, fricts):
+
+            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
+            fd.write(s)
+            
+        #close
+        fd.close()
+
+
+        
+
+
+    def store_timestep(self, t, V0, V1, V2):
+        """Store time, water heights (and momentums) to file
+        """
+        pass
+    
+
+#Auxiliary
+def write_obj(filename,x,y,z):
+    """Store x,y,z vectors into filename (obj format)
+       Vectors are assumed to have dimension (M,3) where
+       M corresponds to the number elements.
+       triangles are assumed to be disconnected
+
+       The three numbers in each vector correspond to three vertices,
+
+       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
+       
+    """
+    #print 'Writing obj to %s' % filename
+
+    import os.path
+
+    root, ext = os.path.splitext(filename)
+    if ext == '.obj':
+        FN = filename
+    else:
+        FN = filename + '.obj'
+        
+
+    outfile = open(FN, 'wb')
+    outfile.write("# Triangulation as an obj file\n")
+
+    M, N = x.shape
+    assert N==3  #Assuming three vertices per element
+
+    for i in range(M):
+        for j in range(N):
+            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))
+
+    for i in range(M):
+        base = i*N
+        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))
+
+    outfile.close()
+
+
+
+#Conversion routines
+def sww2obj(basefilename, size):
+    """Convert netcdf based data output to obj
+    """
+    from Scientific.IO.NetCDF import NetCDFFile        
+
+    from Numeric import Float, zeros
+
+    #Get NetCDF
+    FN = create_filename(basefilename, 'sww', size)
+    print 'Reading from ', FN
+    fid = NetCDFFile(FN, 'r')  #Open existing file for read
+
+        
+    # Get the variables
+    x = fid.variables['x']
+    y = fid.variables['y']
+    z = fid.variables['z']        
+    time = fid.variables['time']
+    stage = fid.variables['stage']
+        
+    M = size  #Number of lines
+    xx = zeros((M,3), Float)
+    yy = zeros((M,3), Float)
+    zz = zeros((M,3), Float)
+
+    for i in range(M):
+        for j in range(3):
+            xx[i,j] = x[i+j*M]
+            yy[i,j] = y[i+j*M]
+            zz[i,j] = z[i+j*M]
+        
+    #Write obj for bathymetry
+    FN = create_filename(basefilename, 'obj', size)        
+    write_obj(FN,xx,yy,zz)
+
+
+    #Now read all the data with variable information, combine with
+    #x,y info and store as obj
+
+    for k in range(len(time)):
+        t = time[k]
+        print 'Processing timestep %f' %t
+
+        for i in range(M):
+            for j in range(3):
+                zz[i,j] = stage[k,i+j*M]
+
+
+        #Write obj for variable data
+        #FN = create_filename(basefilename, 'obj', size, time=t)
+        FN = create_filename(basefilename[:5], 'obj', size, time=t)   
+        write_obj(FN,xx,yy,zz)
+
+
+def dat2obj(basefilename):
+    """Convert line based data output to obj
+    """
+
+    import glob, os
+    from config import data_dir
+   
+    
+    #Get bathymetry and x,y's
+    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
+    
+    from Numeric import zeros, Float
+
+    M = len(lines)  #Number of lines
+    x = zeros((M,3), Float)
+    y = zeros((M,3), Float)
+    z = zeros((M,3), Float)
+
+    ##i = 0
+    for i, line in enumerate(lines):
+        tokens = line.split()
+        values = map(float,tokens)        
+
+        for j in range(3):
+            x[i,j] = values[j*3]
+            y[i,j] = values[j*3+1]
+            z[i,j] = values[j*3+2]        
+        
+        ##i += 1
+
+
+    #Write obj for bathymetry
+    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)
+
+
+    #Now read all the data files with variable information, combine with
+    #x,y info
+    #and store as obj
+
+    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')
+
+    for filename in files:
+        print 'Processing %s' % filename
+
+        lines = open(data_dir+os.sep+filename,'r').readlines()
+        assert len(lines) == M
+        root, ext = os.path.splitext(filename)
+
+        #Get time from filename
+        i0 = filename.find('_time=')
+        if i0 == -1:
+            #Skip bathymetry file
+            continue
+        
+        i0 += 6  #Position where time starts           
+        i1 = filename.find('.dat')	    
+
+        if i1 > i0:
+            t = float(filename[i0:i1])
+        else:
+            raise 'Hmmmm'
+                
+        
+        
+        ##i = 0
+        for i, line in enumerate(lines):
+            tokens = line.split()
+            values = map(float,tokens)        
+
+            for j in range(3):
+                z[i,j] = values[j]        
+        
+            ##i += 1
+
+        #Write obj for variable data
+        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)
+
+
+def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
+    """Read netcdf filename1, pick timesteps first:step:last and save to
+    nettcdf file filename2
+    """
+    from Scientific.IO.NetCDF import NetCDFFile        
+    
+    #Get NetCDF
+    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
+    outfile = NetCDFFile(filename2, 'w')  #Open new file
+    
+
+    #Copy dimensions
+    for d in infile.dimensions:
+        outfile.createDimension(d, infile.dimensions[d])
+
+    for name in infile.variables:
+        var = infile.variables[name]
+        outfile.createVariable(name, var.typecode(), var.dimensions)
+        
+    
+    #Copy the static variables
+    for name in infile.variables:
+        if name == 'time' or name == 'stage':
+            pass
+        else:
+            #Copy
+            outfile.variables[name][:] = infile.variables[name][:]            
+
+    #Copy selected timesteps 
+    time = infile.variables['time']
+    stage = infile.variables['stage']
+
+    newtime = outfile.variables['time']
+    newstage = outfile.variables['stage']    
+
+    if last is None:
+        last = len(time)
+
+    selection = range(first, last, step)
+    for i, j in enumerate(selection):
+        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
+        newtime[i] = time[j]
+        newstage[i,:] = stage[j,:]        
+        
+    #Close
+    infile.close()
+    outfile.close()
+
+
+#Get data objects
+def get_dataobject(domain, mode='w'):
+    """Return instance of class of given format using filename
+    """
+
+    cls = eval('Data_format_%s' %domain.format)
+    return cls(domain, mode)
+
+
+
+
+
+#Mappings from format to reader (FIXME: Soon Obsoleted)
+def get_variable_reader(format):
+    
+    func = eval('%s_variable_reader' %format)
+    return func
+
+
+def get_variable_writer(format):
+    
+    func = eval('%s_variable_writer' %format)
+    return func    
+
+def get_constant_reader(format):
+    
+    func = eval('%s_constant_reader' %format)
+    return func
+
+def get_constant_writer(format):
+    
+    func = eval('%s_constant_writer' %format)
+    return func    
+
Index: inundation/ga/storm_surge/pyvolution/test_data_manager.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 276)
+++ inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 276)
@@ -0,0 +1,426 @@
+#!/usr/bin/env python
+#
+
+import unittest
+from Numeric import zeros, array, allclose, Float
+from util import mean
+
+from data_manager import *
+from shallow_water import *
+from config import epsilon
+
+class dataTestCase(unittest.TestCase):
+    def setUp(self):
+        import time
+        from mesh_factory import rectangular
+
+
+        #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
+        level = zeros(bed.shape, Float)
+
+        h = 0.3
+        for i in range(level.shape[0]):
+            if i % 2 == 0:            
+                level[i,:] = bed[i,:] + h
+            else:
+                level[i,:] = bed[i,:]
+                
+        domain.set_quantity('level', level)
+
+        domain.distribute_to_vertices_and_edges()    
+
+
+        self.domain = domain
+
+        self.X = domain.get_vertex_coordinates()
+        self.X0 = self.X[:,0]
+        self.Y0 = self.X[:,1]
+        self.X1 = self.X[:,2]
+        self.Y1 = self.X[:,3]
+        self.X2 = self.X[:,4]
+        self.Y2 = self.X[:,5]                
+        
+        self.F = bed
+
+        
+    def tearDown(self):
+        pass
+
+
+        
+
+#     def test_xya(self):
+#         import os
+#         from Numeric import concatenate
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+
+#         domain = self.domain
+        
+#         domain.filename = 'datatest' + str(time.time())
+#         domain.format = 'xya'
+#         domain.smooth = True
+        
+#         xya = get_dataobject(self.domain)
+#         xya.store_all()
+
+
+#         #Read back
+#         file = open(xya.filename)
+#         lFile = file.read().split('\n')
+#         lFile = lFile[:-1]
+
+#         file.close()
+#         os.remove(xya.filename)
+
+#         #Check contents
+#         if domain.smooth:
+#             self.failUnless(lFile[0] == '9 3 # <vertex #> <x> <y> [attributes]')
+#         else:    
+#             self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]')            
+
+#         #Get smoothed field values with X and Y
+#         X,Y,F,V = domain.get_vertex_values(xy=True, value_array='field_values',
+#                                            indices = (0,1), precision = Float)
+
+
+#         Q,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
+#                                            indices = (0,), precision = Float)        
+
+
+        
+#         for i, line in enumerate(lFile[1:]):
+#             fields = line.split()
+
+#             assert len(fields) == 5
+
+#             assert allclose(float(fields[0]), X[i])
+#             assert allclose(float(fields[1]), Y[i])
+#             assert allclose(float(fields[2]), F[i,0])
+#             assert allclose(float(fields[3]), Q[i,0])
+#             assert allclose(float(fields[4]), F[i,1])            
+
+
+
+
+#     def test_sww_constant(self):
+#         """Test that constant sww information can be written correctly
+#         (non smooth)
+#         """
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+#         from Scientific.IO.NetCDF import NetCDFFile        
+
+#         self.domain.filename = 'datatest' + str(time.time())
+#         self.domain.format = 'sww'
+#         self.domain.smooth = False
+        
+#         sww = get_dataobject(self.domain)
+#         sww.store_all()
+
+#         #Check contents
+#         #Get NetCDF
+#         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        
+#         # Get the variables
+#         x = fid.variables['x']
+#         y = fid.variables['y']
+#         z = fid.variables['z']
+
+#         volumes = fid.variables['volumes']        
+
+#         assert allclose (x[:], concatenate( (self.X0, self.X1, self.X2) ))
+#         assert allclose (y[:], concatenate( (self.Y0, self.Y1, self.Y2) ))
+#         assert allclose (z[:], concatenate( (self.F0[:,0],
+#                                              self.F1[:,0],
+#                                              self.F2[:,0]) ))
+
+#         V = volumes
+
+#         P = len(self.domain)
+#         for k in range(P):
+#              assert V[k,0] == V[k,1] - P
+#              assert V[k,1] == V[k,2] - P
+        
+
+
+        
+#         fid.close()
+        
+#         #Cleanup
+#         os.remove(sww.filename)
+
+
+#     def test_sww_constant_smooth(self):
+#         """Test that constant sww information can be written correctly
+#         (non smooth)
+#         """
+
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+#         from Scientific.IO.NetCDF import NetCDFFile        
+
+#         self.domain.filename = 'datatest' + str(time.time())
+#         self.domain.format = 'sww'
+#         self.domain.smooth = True
+        
+#         sww = get_dataobject(self.domain)
+#         sww.store_all()
+
+#         #Check contents
+#         #Get NetCDF
+#         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        
+#         # Get the variables
+#         x = fid.variables['x']
+#         y = fid.variables['y']
+#         z = fid.variables['z']
+
+#         volumes = fid.variables['volumes']        
+
+#         X = x[:]
+#         Y = y[:]
+        
+#         assert allclose([X[0], Y[0]], array([0.0, 0.0]))
+#         assert allclose([X[1], Y[1]], array([0.0, 0.5]))
+#         assert allclose([X[2], Y[2]], array([0.0, 1.0]))
+
+#         assert allclose([X[4], Y[4]], array([0.5, 0.5]))
+
+#         assert allclose([X[7], Y[7]], array([1.0, 0.5]))        
+
+#         Z = z[:]
+#         assert Z[4] == -0.5
+
+
+
+#         V = volumes
+#         assert V[2,0] == 4
+#         assert V[2,1] == 5
+#         assert V[2,2] == 1
+
+#         assert V[4,0] == 6
+#         assert V[4,1] == 7
+#         assert V[4,2] == 3                      
+        
+        
+
+
+        
+#         fid.close()
+        
+#         #Cleanup
+#         os.remove(sww.filename)
+
+
+
+#     def test_sww_variable(self):
+#         """Test that sww information can be written correctly
+#         """
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+#         from Scientific.IO.NetCDF import NetCDFFile        
+
+#         self.domain.filename = 'datatest' + str(time.time())
+#         self.domain.format = 'sww'
+#         self.domain.smooth = True
+#         self.domain.reduction = mean             
+        
+#         sww = get_dataobject(self.domain)
+#         sww.store_all()
+#         sww.store_timestep()
+
+#         #Check contents
+#         #Get NetCDF
+#         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        
+
+#         # Get the variables
+#         x = fid.variables['x']
+#         y = fid.variables['y']
+#         z = fid.variables['z']        
+#         time = fid.variables['time']
+#         stage = fid.variables['stage']
+
+
+#         Q0 = self.domain.volume_class.conserved_quantities_vertex0
+#         Q1 = self.domain.volume_class.conserved_quantities_vertex1
+#         Q2 = self.domain.volume_class.conserved_quantities_vertex2        
+
+#         A = stage[0,:]
+#         #print A[0], (Q2[0,0] + Q1[1,0])/2
+#         assert allclose(A[0], (Q2[0,0] + Q1[1,0])/2)  
+#         assert allclose(A[1], (Q0[1,0] + Q1[3,0] + Q2[2,0])/3)
+#         assert allclose(A[2], Q0[3,0])
+#         assert allclose(A[3], (Q0[0,0] + Q1[5,0] + Q2[4,0])/3)
+
+#         #Center point
+#         assert allclose(A[4], (Q1[0,0] + Q2[1,0] + Q0[2,0] +\
+#                                  Q0[5,0] + Q2[6,0] + Q1[7,0])/6)
+                                 
+        
+        
+#         fid.close()
+        
+#         #Cleanup
+#         os.remove(sww.filename)
+        
+        
+#     def test_sww_variable2(self):
+#         """Test that sww information can be written correctly
+#         multiple timesteps. Use average as reduction operator
+#         """
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+#         from Scientific.IO.NetCDF import NetCDFFile
+
+#         self.domain.filename = 'datatest' + str(time.time())
+#         self.domain.format = 'sww'
+#         self.domain.smooth = True
+
+#         self.domain.reduction = mean     
+        
+#         sww = get_dataobject(self.domain)
+#         sww.store_all()
+#         sww.store_timestep()
+        
+#         self.domain.evolve_to_end(max_timestep = 0.01, finaltime = 0.01)
+#         sww.store_timestep()
+
+
+
+#         #Check contents
+#         #Get NetCDF
+#         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        
+
+#         # Get the variables
+#         x = fid.variables['x']
+#         y = fid.variables['y']
+#         z = fid.variables['z']
+#         time = fid.variables['time']
+#         stage = fid.variables['stage']        
+
+#         #Check values
+
+#         Q0 = self.domain.volume_class.conserved_quantities_vertex0
+#         Q1 = self.domain.volume_class.conserved_quantities_vertex1
+#         Q2 = self.domain.volume_class.conserved_quantities_vertex2        
+
+#         A = stage[1,:]
+#         assert allclose(A[0], (Q2[0,0] + Q1[1,0])/2)  
+#         assert allclose(A[1], (Q0[1,0] + Q1[3,0] + Q2[2,0])/3)
+#         assert allclose(A[2], Q0[3,0])
+#         assert allclose(A[3], (Q0[0,0] + Q1[5,0] + Q2[4,0])/3)
+
+#         #Center point
+#         assert allclose(A[4], (Q1[0,0] + Q2[1,0] + Q0[2,0] +\
+#                                  Q0[5,0] + Q2[6,0] + Q1[7,0])/6)
+                                 
+        
+        
+            
+        
+#         fid.close()
+        
+#         #Cleanup
+#         os.remove(sww.filename)
+        
+#     def test_sww_variable3(self):
+#         """Test that sww information can be written correctly
+#         multiple timesteps using a different reduction operator (min)
+#         """
+
+#         import time, os
+#         from Numeric import array, zeros, allclose, Float, concatenate
+#         from Scientific.IO.NetCDF import NetCDFFile
+
+#         self.domain.filename = 'datatest' + str(time.time())
+#         self.domain.format = 'sww'
+#         self.domain.smooth = True
+#         self.domain.reduction = min
+        
+#         sww = get_dataobject(self.domain)
+#         sww.store_all()
+#         sww.store_timestep()
+        
+#         self.domain.evolve_to_end(max_timestep = 0.01, finaltime = 0.01)
+#         sww.store_timestep()
+
+
+
+#         #Check contents
+#         #Get NetCDF
+#         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        
+
+#         # Get the variables
+#         x = fid.variables['x']
+#         y = fid.variables['y']
+#         z = fid.variables['z']
+#         time = fid.variables['time']
+#         stage = fid.variables['stage']        
+
+#         #Check values
+
+#         Q0 = self.domain.volume_class.conserved_quantities_vertex0
+#         Q1 = self.domain.volume_class.conserved_quantities_vertex1
+#         Q2 = self.domain.volume_class.conserved_quantities_vertex2        
+
+#         A = stage[1,:]
+#         assert allclose(A[0], min(Q2[0,0], Q1[1,0]) )  
+#         assert allclose(A[1], min(Q0[1,0], Q1[3,0], Q2[2,0]) )
+#         assert allclose(A[2], Q0[3,0])
+#         assert allclose(A[3], min(Q0[0,0], Q1[5,0], Q2[4,0]) )
+
+#         #Center point
+#         assert allclose(A[4], min(Q1[0,0], Q2[1,0], Q0[2,0],
+#                                   Q0[5,0], Q2[6,0], Q1[7,0]) )
+                                 
+        
+        
+            
+        
+#         fid.close()
+        
+#         #Cleanup
+#         os.remove(sww.filename)
+
+
+    def test_dummy(self):
+        pass
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    #suite = unittest.makeSuite(dataTestCase,'test_sww_constant')
+    suite = unittest.makeSuite(dataTestCase,'test_dummy')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
