Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1018)
@@ -5,5 +5,5 @@
   2: Variable data: Conserved quantities. Stored once per timestep.
 
-All data is assumed to reside at vertex locations.  
+All data is assumed to reside at vertex locations.
 
 
@@ -38,5 +38,5 @@
 PTS + TSH -> TSH with elevation: Least squares fit
 
-TSH -> SWW:  Conversion of TSH to sww viewable using Swollen 
+TSH -> SWW:  Conversion of TSH to sww viewable using Swollen
 
 TSH + Boundary SWW -> SWW: SImluation using pyvolution
@@ -45,5 +45,5 @@
 """
 
-from Numeric import concatenate        
+from Numeric import concatenate
 
 
@@ -56,6 +56,6 @@
     s = s.replace('(', '')
     s = s.replace(')', '')
-    s = s.replace('__', '_')                
-        
+    s = s.replace('__', '_')
+
     return s
 
@@ -63,5 +63,5 @@
 def check_dir(path, verbose=None):
     """Check that specified path exists.
-    If path does not exist it will be created if possible    
+    If path does not exist it will be created if possible
 
     USAGE:
@@ -74,5 +74,5 @@
     RETURN VALUE:
         Verified path including trailing separator
-        
+
     """
 
@@ -86,5 +86,5 @@
 
 
-    if path[-1] != os.sep: 
+    if path[-1] != os.sep:
         path = path + os.sep  # Add separator for directories
 
@@ -100,7 +100,7 @@
             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.'
@@ -109,10 +109,10 @@
             else:
                 path = 'C:'
-                
+
             print 'Using directory %s instead' %path
 
     return(path)
 
-    
+
 
 def del_dir(path):
@@ -128,5 +128,5 @@
 
             if os.path.isdir(X) and not os.path.islink(X):
-                del_dir(X) 
+                del_dir(X)
             else:
                 try:
@@ -138,21 +138,21 @@
 
 
-	
+
 def create_filename(datadir, filename, format, size=None, time=None):
 
     import os
     #from config import data_dir
-           
+
     FN = check_dir(datadir) + filename
-    
-    if size is not None:	   
+
+    if size is not None:
     	FN += '_size%d' %size
 
     if time is not None:
         FN += '_time%.2f' %time
-        
+
     FN += '.' + format
     return FN
-    
+
 
 def get_files(datadir, filename, format, size):
@@ -164,5 +164,5 @@
     import os
     #from config import data_dir
-           
+
     dir = check_dir(datadir)
 
@@ -177,18 +177,18 @@
     """
 
-        
+
     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)'''    
+                                        '''   'a' (append)'''
 
         #Create filename
-        #self.filename = create_filename(domain.get_datadir(), 
-	#	domain.get_name(), extension, len(domain))
-
-	
-        self.filename = create_filename(domain.get_datadir(), 
-			domain.get_name(), extension)	
+        #self.filename = create_filename(domain.get_datadir(),
+        #domain.get_name(), extension, len(domain))
+
+
+        self.filename = create_filename(domain.get_datadir(),
+			domain.get_name(), extension)
 
         #print 'F', self.filename
@@ -196,8 +196,8 @@
         self.number_of_volumes = len(domain)
         self.domain = domain
-    
-    
-    #FIXME: Should we have a general set_precision function?
-    
+
+
+        #FIXME: Should we have a general set_precision function?
+
 
 
@@ -207,11 +207,11 @@
     """
 
-        
+
     def __init__(self, domain, mode = 'w'):
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import Int, Float, Float32 
+        from Numeric import Int, Float, Float32
 
         self.precision = Float32 #Use single precision
-        
+
         Data_format.__init__(self, domain, 'sww', mode)
 
@@ -219,21 +219,21 @@
         # NetCDF file definition
         fid = NetCDFFile(self.filename, mode)
-        
+
         if mode == 'w':
-            
+
             #Create new file
             fid.institution = 'Geoscience Australia'
             fid.description = 'Output from pyvolution suitable for plotting'
-            
+
             if domain.smooth:
                 fid.smoothing = 'Yes'
-            else:    
+            else:
                 fid.smoothing = 'No'
 
-            fid.order = domain.default_order    
-
-	    #Reference point
-	    #Start time in seconds since the epoch (midnight 1/1/1970)
-	    fid.starttime = domain.starttime
+            fid.order = domain.default_order
+
+            #Reference point
+            #Start time in seconds since the epoch (midnight 1/1/1970)
+            fid.starttime = domain.starttime
             fid.xllcorner = domain.xllcorner
             fid.yllcorner = domain.yllcorner
@@ -241,12 +241,12 @@
 
             # dimension definitions
-            fid.createDimension('number_of_volumes', self.number_of_volumes)  
+            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:    
+            else:
                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
-                
+
             fid.createDimension('number_of_timesteps', None) #extensible
 
@@ -255,9 +255,9 @@
             fid.createVariable('y', self.precision, ('number_of_points',))
             fid.createVariable('elevation', self.precision, ('number_of_points',))
-            
+
 	    #FIXME: Backwards compatibility
 	    fid.createVariable('z', self.precision, ('number_of_points',))
 	    #################################
-        
+
             fid.createVariable('volumes', Int, ('number_of_volumes',
                                                 'number_of_vertices'))
@@ -265,5 +265,5 @@
             fid.createVariable('time', self.precision,
                                ('number_of_timesteps',))
-        
+
             fid.createVariable('stage', self.precision,
                                ('number_of_timesteps',
@@ -279,5 +279,5 @@
 
         #Close
-        fid.close()                
+        fid.close()
 
 
@@ -297,10 +297,10 @@
         #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['elevation']
-        
+
         volumes = fid.variables['volumes']
 
@@ -309,5 +309,5 @@
         X,Y,Z,V = Q.get_vertex_values(xy=True,
 				      precision = self.precision)
-				      
+
 
 
@@ -315,12 +315,12 @@
         y[:] = Y.astype(self.precision)
         z[:] = Z.astype(self.precision)
-	
+
 	#FIXME: Backwards compatibility
         z = fid.variables['z']
-        z[:] = Z.astype(self.precision)	
+        z[:] = Z.astype(self.precision)
 	################################
 
         volumes[:] = V.astype(volumes.typecode())
-        
+
         #Close
         fid.close()
@@ -334,10 +334,10 @@
 	from time import sleep
 
-		
+
         #Get NetCDF
 	retries = 0
 	file_open = False
 	while not file_open and retries < 10:
-  	    try:	
+  	    try:
                 fid = NetCDFFile(self.filename, 'a') #Open existing file
 	    except IOError:
@@ -352,17 +352,17 @@
 	    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']
         xmomentum = fid.variables['xmomentum']
-        ymomentum = fid.variables['ymomentum']                        
+        ymomentum = fid.variables['ymomentum']
         i = len(time)
 
@@ -374,5 +374,5 @@
             names = [names]
 
-        for name in names:    
+        for name in names:
             # Get quantity
             Q = domain.quantities[name]
@@ -386,11 +386,11 @@
                 xmomentum[i,:] = A.astype(self.precision)
             elif name == 'ymomentum':
-                ymomentum[i,:] = A.astype(self.precision)   
-		                     
+                ymomentum[i,:] = A.astype(self.precision)
+
             #As in....
 	    #eval( name + '[i,:] = A.astype(self.precision)' )
 	    #FIXME: But we need a UNIT test for that before refactoring
-	    
-	    
+
+
 
         #Flush and close
@@ -404,11 +404,11 @@
     """
 
-        
+
     def __init__(self, domain, mode = 'w'):
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import Int, Float, Float 
+        from Numeric import Int, Float, Float
 
         self.precision = Float #Use full precision
-        
+
         Data_format.__init__(self, domain, 'sww', mode)
 
@@ -416,14 +416,14 @@
         # NetCDF file definition
         fid = NetCDFFile(self.filename, mode)
-        
+
         if mode == 'w':
             #Create new file
             fid.institution = 'Geoscience Australia'
             fid.description = 'Checkpoint data'
-            #fid.smooth = domain.smooth                
-            fid.order = domain.default_order    
+            #fid.smooth = domain.smooth
+            fid.order = domain.default_order
 
             # dimension definitions
-            fid.createDimension('number_of_volumes', self.number_of_volumes)  
+            fid.createDimension('number_of_volumes', self.number_of_volumes)
             fid.createDimension('number_of_vertices', 3)
 
@@ -438,5 +438,5 @@
             fid.createVariable('y', self.precision, ('number_of_points',))
 
-        
+
             fid.createVariable('volumes', Int, ('number_of_volumes',
                                                 'number_of_vertices'))
@@ -452,5 +452,5 @@
 
         #Close
-        fid.close()                
+        fid.close()
 
 
@@ -471,9 +471,9 @@
         #Get NetCDF
         fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
-        
+
         # Get the variables
         x = fid.variables['x']
         y = fid.variables['y']
-        
+
         volumes = fid.variables['volumes']
 
@@ -482,5 +482,5 @@
         X,Y,Z,V = Q.get_vertex_values(xy=True,
 				      precision = self.precision)
-				      
+
 
 
@@ -490,5 +490,5 @@
 
         volumes[:] = V
-        
+
         #Close
         fid.close()
@@ -500,10 +500,10 @@
         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:	
+  	    try:
                 fid = NetCDFFile(self.filename, 'a') #Open existing file
 	    except IOError:
@@ -518,15 +518,15 @@
 	    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']        
+        stage = fid.variables['stage']
         i = len(time)
 
@@ -538,5 +538,5 @@
         A,V = Q.get_vertex_values(xy=False,
 				  precision = self.precision)
-        
+
         stage[i,:] = A.astype(self.precision)
 
@@ -560,13 +560,13 @@
     def __init__(self, domain, mode = 'w'):
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import Int, Float, Float32 
+        from Numeric import Int, Float, Float32
 
         self.precision = Float32 #Use single precision
-        
+
         Data_format.__init__(self, domain, 'xya', mode)
-        
-    
-    
-    #FIXME -This is the old xya format  
+
+
+
+    #FIXME -This is the old xya format
     def store_all(self):
         """Specialisation of store all for xya format
@@ -579,5 +579,5 @@
 
         domain = self.domain
-        
+
         fd = open(self.filename, 'w')
 
@@ -585,5 +585,5 @@
         if domain.smooth is True:
             number_of_points =  len(domain.vertexlist)
-        else:    
+        else:
             number_of_points = 3*self.number_of_volumes
 
@@ -613,5 +613,5 @@
             s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
             fd.write(s)
-            
+
         #close
         fd.close()
@@ -622,5 +622,5 @@
         """
         pass
-    
+
 
 #Auxiliary
@@ -634,5 +634,5 @@
 
        e.g. the x coordinate of vertex 1 of element i is in x[i,1]
-       
+
     """
     #print 'Writing obj to %s' % filename
@@ -645,5 +645,5 @@
     else:
         FN = filename + '.obj'
-        
+
 
     outfile = open(FN, 'wb')
@@ -669,5 +669,5 @@
     """Convert netcdf based data output to obj
     """
-    from Scientific.IO.NetCDF import NetCDFFile        
+    from Scientific.IO.NetCDF import NetCDFFile
 
     from Numeric import Float, zeros
@@ -678,12 +678,12 @@
     fid = NetCDFFile(FN, 'r')  #Open existing file for read
 
-        
+
     # Get the variables
     x = fid.variables['x']
     y = fid.variables['y']
-    z = fid.variables['elevation']        
+    z = fid.variables['elevation']
     time = fid.variables['time']
     stage = fid.variables['stage']
-        
+
     M = size  #Number of lines
     xx = zeros((M,3), Float)
@@ -696,7 +696,7 @@
             yy[i,j] = y[i+j*M]
             zz[i,j] = z[i+j*M]
-        
+
     #Write obj for bathymetry
-    FN = create_filename('.', basefilename, 'obj', size)        
+    FN = create_filename('.', basefilename, 'obj', size)
     write_obj(FN,xx,yy,zz)
 
@@ -716,5 +716,5 @@
         #Write obj for variable data
         #FN = create_filename(basefilename, 'obj', size, time=t)
-        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)   
+        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
         write_obj(FN,xx,yy,zz)
 
@@ -727,9 +727,9 @@
     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
 
@@ -742,11 +742,11 @@
     for i, line in enumerate(lines):
         tokens = line.split()
-        values = map(float,tokens)        
+        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]        
-        
+            z[i,j] = values[j*3+2]
+
         ##i += 1
 
@@ -774,7 +774,7 @@
             #Skip bathymetry file
             continue
-        
-        i0 += 6  #Position where time starts           
-        i1 = filename.find('.dat')	    
+
+        i0 += 6  #Position where time starts
+        i1 = filename.find('.dat')
 
         if i1 > i0:
@@ -782,15 +782,15 @@
         else:
             raise 'Hmmmm'
-                
-        
-        
+
+
+
         ##i = 0
         for i, line in enumerate(lines):
             tokens = line.split()
-            values = map(float,tokens)        
+            values = map(float,tokens)
 
             for j in range(3):
-                z[i,j] = values[j]        
-        
+                z[i,j] = values[j]
+
             ##i += 1
 
@@ -803,10 +803,10 @@
     nettcdf file filename2
     """
-    from Scientific.IO.NetCDF import NetCDFFile        
-    
+    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
@@ -817,6 +817,6 @@
         var = infile.variables[name]
         outfile.createVariable(name, var.typecode(), var.dimensions)
-        
-    
+
+
     #Copy the static variables
     for name in infile.variables:
@@ -825,12 +825,12 @@
         else:
             #Copy
-            outfile.variables[name][:] = infile.variables[name][:]            
-
-    #Copy selected timesteps 
+            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']    
+    newstage = outfile.variables['stage']
 
     if last is None:
@@ -841,6 +841,6 @@
         print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
         newtime[i] = time[j]
-        newstage[i,:] = stage[j,:]        
-        
+        newstage[i,:] = stage[j,:]
+
     #Close
     infile.close()
@@ -873,9 +873,9 @@
     points:  (Nx2) Float array
     elevation: N Float array
-    """ 
+    """
 
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, arrayrange, concatenate    
+    from Numeric import Float, arrayrange, concatenate
 
     root = basename_in
@@ -885,5 +885,5 @@
 
     if verbose: print 'Reading DEM from %s' %(root + '.dem')
-    
+
     ncols = infile.ncols[0]
     nrows = infile.nrows[0]
@@ -902,5 +902,5 @@
     datum = infile.datum
     units = infile.units
-    
+
 
     #Get output file
@@ -909,9 +909,9 @@
     else:
         ptsname = basename_out + '.pts'
-        
+
     if verbose: print 'Store to NetCDF file %s' %ptsname
     # NetCDF file definition
     outfile = NetCDFFile(ptsname, 'w')
-        
+
     #Create new file
     outfile.institution = 'Geoscience Australia'
@@ -925,5 +925,5 @@
     outfile.false_easting = false_easting
     outfile.false_northing =false_northing
-    
+
     outfile.projection = projection
     outfile.datum = datum
@@ -935,9 +935,9 @@
     outfile.nrows = nrows
 
-    
+
     # dimension definitions
-    outfile.createDimension('number_of_points', nrows*ncols)    
+    outfile.createDimension('number_of_points', nrows*ncols)
     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-    
+
     # variable definitions
     outfile.createVariable('points', Float, ('number_of_points',
@@ -953,17 +953,17 @@
     for i in range(nrows):
         if verbose: print 'Processing row %d of %d' %(i, nrows)
-        
-        y = (nrows-i)*cellsize            
+
+        y = (nrows-i)*cellsize
         for j in range(ncols):
             index = i*ncols + j
-            
+
             x = j*cellsize
             points[index, :] = [x,y]
-            elevation[index] = dem_elevation[i, j]            
-
-    infile.close()            
+            elevation[index] = dem_elevation[i, j]
+
+    infile.close()
     outfile.close()
 
-				      
+
 def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
                                   verbose=False):
@@ -988,5 +988,5 @@
 
     The prj format is assumed to be as
-    
+
     Projection    UTM
     Zone          56
@@ -998,5 +998,5 @@
     Yshift        10000000.0000000000
     Parameters
-    """ 
+    """
 
     import os
@@ -1010,5 +1010,5 @@
     # Read Meta data
     if verbose: print 'Reading METADATA from %s' %root + '.prj'
-    metadatafile = open(root + '.prj')    
+    metadatafile = open(root + '.prj')
     metalines = metadatafile.readlines()
     metadatafile.close()
@@ -1020,11 +1020,11 @@
     L = metalines[1].strip().split()
     assert L[0].strip().lower() == 'zone'
-    zone = int(L[1].strip())    
+    zone = int(L[1].strip())
 
     L = metalines[2].strip().split()
     assert L[0].strip().lower() == 'datum'
-    datum = L[1].strip()                        #TEXT 
-
-    L = metalines[3].strip().split()      
+    datum = L[1].strip()                        #TEXT
+
+    L = metalines[3].strip().split()
     assert L[0].strip().lower() == 'zunits'     #IGNORE
     zunits = L[1].strip()                       #TEXT
@@ -1039,10 +1039,10 @@
 
     L = metalines[6].strip().split()
-    assert L[0].strip().lower() == 'xshift'     
+    assert L[0].strip().lower() == 'xshift'
     false_easting = float(L[1].strip())
 
     L = metalines[7].strip().split()
     assert L[0].strip().lower() == 'yshift'
-    false_northing = float(L[1].strip())    
+    false_northing = float(L[1].strip())
 
     #print false_easting, false_northing, zone, datum
@@ -1051,5 +1051,5 @@
     ###########################################
     #Read DEM data
-    
+
     datafile = open(basename_in + '.asc')
 
@@ -1059,5 +1059,5 @@
 
     if verbose: print 'Got', len(lines), ' lines'
-    
+
     ncols = int(lines[0].split()[1].strip())
     nrows = int(lines[1].split()[1].strip())
@@ -1076,10 +1076,10 @@
         netcdfname = root + '.dem'
     else:
-        netcdfname = basename_out + '.dem'        
-        
+        netcdfname = basename_out + '.dem'
+
     if verbose: print 'Store to NetCDF file %s' %netcdfname
     # NetCDF file definition
     fid = NetCDFFile(netcdfname, 'w')
-        
+
     #Create new file
     fid.institution = 'Geoscience Australia'
@@ -1096,13 +1096,13 @@
     fid.zone = zone
     fid.false_easting = false_easting
-    fid.false_northing = false_northing    
+    fid.false_northing = false_northing
     fid.projection = projection
     fid.datum = datum
     fid.units = units
-    
+
 
     # dimension definitions
     fid.createDimension('number_of_rows', nrows)
-    fid.createDimension('number_of_columns', ncols)            
+    fid.createDimension('number_of_columns', ncols)
 
     # variable definitions
@@ -1122,5 +1122,5 @@
     fid.close()
 
-    
+
 
 def ferret2sww(basename_in, basename_out = None,
@@ -1153,10 +1153,10 @@
     ferret2sww uses grid points as vertices in a triangular grid
     counting vertices from lower left corner upwards, then right
-    """ 
+    """
 
     import os
     from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import Float, Int, Int32, searchsorted, zeros, array
-    precision = Float 
+    precision = Float
 
 
@@ -1165,10 +1165,10 @@
     file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
     file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
-    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
+    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
 
     if basename_out is None:
         swwname = basename_in + '.sww'
     else:
-        swwname = basename_out + '.sww'        
+        swwname = basename_out + '.sww'
 
     times = file_h.variables['TIME']
@@ -1178,12 +1178,12 @@
     if mint == None:
         jmin = 0
-    else:    
+    else:
         jmin = searchsorted(times, mint)
-    
+
     if maxt == None:
         jmax=len(times)
-    else:    
+    else:
         jmax = searchsorted(times, maxt)
-        
+
     if minlat == None:
         kmin=0
@@ -1198,13 +1198,13 @@
     if minlon == None:
         lmin=0
-    else:    
+    else:
         lmin = searchsorted(longitudes, minlon)
-    
+
     if maxlon == None:
         lmax = len(longitudes)
     else:
-        lmax = searchsorted(longitudes, maxlon)            
-
-    times = times[jmin:jmax]        
+        lmax = searchsorted(longitudes, maxlon)
+
+    times = times[jmin:jmax]
     latitudes = latitudes[kmin:kmax]
     longitudes = longitudes[lmin:lmax]
@@ -1214,5 +1214,5 @@
     amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
     xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax]
-    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]    
+    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]
 
     number_of_times = times.shape[0]
@@ -1222,5 +1222,5 @@
     assert amplitudes.shape[0] == number_of_times
     assert amplitudes.shape[1] == number_of_latitudes
-    assert amplitudes.shape[2] == number_of_longitudes          
+    assert amplitudes.shape[2] == number_of_longitudes
 
     #print times
@@ -1229,16 +1229,16 @@
 
     #print 'MIN', min(min(min(amplitudes)))
-    #print 'MAX', max(max(max(amplitudes)))    
+    #print 'MAX', max(max(max(amplitudes)))
 
     #print number_of_latitudes, number_of_longitudes
     number_of_points = number_of_latitudes*number_of_longitudes
     number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
-    
+
     #print file_h.dimensions.keys()
-    #print file_h.variables.keys()    
+    #print file_h.variables.keys()
 
     file_h.close()
     file_u.close()
-    file_v.close()    
+    file_v.close()
 
 
@@ -1246,5 +1246,5 @@
     # NetCDF file definition
     outfile = NetCDFFile(swwname, 'w')
-        
+
     #Create new file
     outfile.institution = 'Geoscience Australia'
@@ -1256,10 +1256,10 @@
 
     #For sww compatibility
-    outfile.smoothing = 'Yes' 
+    outfile.smoothing = 'Yes'
     outfile.order = 1
-	    
+
     #Start time in seconds since the epoch (midnight 1/1/1970)
     outfile.starttime = times[0]
-    
+
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
@@ -1267,6 +1267,6 @@
     outfile.createDimension('number_of_vertices', 3)
     outfile.createDimension('number_of_points', number_of_points)
-                            
-                
+
+
     #outfile.createDimension('number_of_timesteps', len(times))
     outfile.createDimension('number_of_timesteps', len(times))
@@ -1276,15 +1276,15 @@
     outfile.createVariable('y', precision, ('number_of_points',))
     outfile.createVariable('elevation', precision, ('number_of_points',))
-            
+
     #FIXME: Backwards compatibility
     outfile.createVariable('z', precision, ('number_of_points',))
     #################################
-        
+
     outfile.createVariable('volumes', Int, ('number_of_volumes',
                                             'number_of_vertices'))
-    
+
     outfile.createVariable('time', precision,
                            ('number_of_timesteps',))
-        
+
     outfile.createVariable('stage', precision,
                            ('number_of_timesteps',
@@ -1294,5 +1294,5 @@
                            ('number_of_timesteps',
                             'number_of_points'))
-    
+
     outfile.createVariable('ymomentum', precision,
                            ('number_of_timesteps',
@@ -1308,13 +1308,13 @@
 
     #Check zone boundaries
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])    
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
 
     vertices = {}
-    for k, lat in enumerate(latitudes):    
-        for l, lon in enumerate(longitudes):    
+    for k, lat in enumerate(latitudes):
+        for l, lon in enumerate(longitudes):
 
             vertices[l,k] = i
 
-            zone, easting, northing = redfearn(lat,lon)                
+            zone, easting, northing = redfearn(lat,lon)
 
             msg = 'Zone boundary crossed at longitude =', lon
@@ -1331,15 +1331,15 @@
         for k in range(number_of_latitudes-1):
             v1 = vertices[l,k+1]
-            v2 = vertices[l,k]            
-            v3 = vertices[l+1,k+1]            
-            v4 = vertices[l+1,k]            
+            v2 = vertices[l,k]
+            v3 = vertices[l+1,k+1]
+            v4 = vertices[l+1,k]
 
             volumes.append([v1,v2,v3]) #Upper element
-            volumes.append([v4,v3,v2]) #Lower element            
-
-    volumes = array(volumes)        
+            volumes.append([v4,v3,v2]) #Lower element
+
+    volumes = array(volumes)
 
     if origin == None:
-        zone = refzone        
+        zone = refzone
         xllcorner = min(x)
         yllcorner = min(y)
@@ -1347,7 +1347,7 @@
         zone = origin[0]
         xllcorner = origin[1]
-        yllcorner = origin[2]                
-
-    
+        yllcorner = origin[2]
+
+
     outfile.xllcorner = xllcorner
     outfile.yllcorner = yllcorner
@@ -1358,7 +1358,7 @@
     outfile.variables['z'][:] = 0.0
     outfile.variables['elevation'][:] = 0.0
-    outfile.variables['time'][:] = times        
+    outfile.variables['time'][:] = times
     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
-    
+
 
 
@@ -1366,26 +1366,26 @@
     stage = outfile.variables['stage']
     xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum']        
+    ymomentum = outfile.variables['ymomentum']
 
     for j in range(len(times)):
         i = 0
         for k in range(number_of_latitudes):
-            for l in range(number_of_longitudes):                    
+            for l in range(number_of_longitudes):
                 h = zscale*amplitudes[j,k,l]/100 + mean_stage
                 stage[j,i] = h
                 xmomentum[j,i] = xspeed[j,k,l]/100*h
-                ymomentum[j,i] = yspeed[j,k,l]/100*h                
+                ymomentum[j,i] = yspeed[j,k,l]/100*h
                 i += 1
-        
-    outfile.close()                
-    
- 
+
+    outfile.close()
+
+
 def extent_sww(file_name):
     """
     Read in an sww file.
-    
+
     Input;
     file_name - the sww file
-    
+
     Output;
     z - Vector of bed elevation
@@ -1395,12 +1395,12 @@
     stage - array with respect to time and vertices (x,y)
     """
-        
-    
-    from Scientific.IO.NetCDF import NetCDFFile        	
-    
+
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
     #Check contents
     #Get NetCDF
-    fid = NetCDFFile(file_name, 'r') 
-    
+    fid = NetCDFFile(file_name, 'r')
+
     # Get the variables
     x = fid.variables['x'][:]
@@ -1413,5 +1413,5 @@
 
     fid.close()
-    
+
     return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
 
@@ -1426,17 +1426,17 @@
 
     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])                        
+            fid.write('%.16e ' %v0[i,j])
         for j in range(N):
-            fid.write('%.16e ' %v2[i,j])                        
-            
+            fid.write('%.16e ' %v1[i,j])
+        for j in range(N):
+            fid.write('%.16e ' %v2[i,j])
+
         fid.write('\n')
     fid.close()
@@ -1448,8 +1448,8 @@
 
     M, N = v0.shape
-	
+
     FN = create_filename(filename, 'cpt', M, t)
     #print 'Reading from %s' %FN
-    
+
     fid = open(FN)
 
@@ -1462,5 +1462,5 @@
             v1[i,j] = float(values[3+j])
             v2[i,j] = float(values[6+j])
-    
+
     fid.close()
 
@@ -1474,25 +1474,25 @@
 
     print X0
-    import sys; sys.exit() 
+    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):        
+        for j in range(2):
             fid.write('%.16e ' %X0[i,j])   #x, y
-        for j in range(N):                    
+        for j in range(N):
             fid.write('%.16e ' %v0[i,j])       #z,z,z,
-            
-        for j in range(2):                    
+
+        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):                                
+        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]) 
-            
+        for j in range(N):
+            fid.write('%.16e ' %v2[i,j])
+
         fid.write('\n')
     fid.close()
@@ -1508,22 +1508,22 @@
 
     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):        
+        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):                    
+
+        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):                                
+
+        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()
@@ -1536,15 +1536,14 @@
 
     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('%.4f ' %v2[i,0])
 
         fid.write('\n')
     fid.close()
-
Index: /inundation/ga/storm_surge/pyvolution/test_advection.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_advection.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_advection.py	(revision 1018)
@@ -8,9 +8,9 @@
 from advection import *
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_Advection(unittest.TestCase):
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
@@ -25,10 +25,10 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
-        
-        domain = Domain(points, vertices)        
+
+        domain = Domain(points, vertices)
         domain.check_integrity()
-        
+
         assert domain.quantities.has_key('stage')
 
@@ -38,15 +38,15 @@
     def test_flux_1_triangle0(self):
         a = [0.0, 0.5]
-        b = [0.0, 0.0]                
+        b = [0.0, 0.0]
         c = [0.5, 0.5]
 
         points = [a, b, c]
         vertices = [ [0,1,2] ]
-        domain = Domain(points, vertices)        
+        domain = Domain(points, vertices)
         domain.check_integrity()
 
 
         #Populate boundary array with dirichlet conditions.
- 	domain.neighbours = array([[-1,-2,-3]])	        
+ 	domain.neighbours = array([[-1,-2,-3]])
         domain.quantities['stage'].boundary_values[:] = 1.0
 
@@ -54,5 +54,5 @@
 
 	domain.distribute_to_vertices_and_edges() #Use first order default
-        
+
         domain.check_integrity()
 
@@ -61,16 +61,16 @@
         R = -0.5/domain.areas[0]
 
-        assert U==R, '%s %s' %(U, R) 
+        assert U==R, '%s %s' %(U, R)
 
 
     def test_flux_1_triangle1(self):
- 
+
         a = [0.0, 0.5]
-        b = [0.0, 0.0]                
+        b = [0.0, 0.0]
         c = [0.5, 0.5]
 
         points = [a, b, c]
         vertices = [ [0,1,2] ]
-        domain = Domain(points, vertices)        
+        domain = Domain(points, vertices)
         domain.check_integrity()
 
@@ -79,32 +79,32 @@
 	domain.distribute_to_vertices_and_edges()
         domain.check_integrity()
-        
-        
+
+
 	domain.compute_fluxes()
         U = -domain.quantities['stage'].explicit_update
         R = 0.5/domain.areas[0]
 
-        assert U==R, '%s %s' %(U, R) 
+        assert U==R, '%s %s' %(U, R)
 
 
 
     def test_flux_1_triangle2(self):
-        
+
         a = [0.0, 0.5]
-        b = [0.0, 0.0]                
+        b = [0.0, 0.0]
         c = [0.5, 0.5]
 
         points = [a, b, c]
         vertices = [ [0,1,2] ]
-        domain = Domain(points, vertices)        
+        domain = Domain(points, vertices)
         domain.check_integrity()
 
-		
+
         #Populate boundary array with dirichlet conditions.
- 	domain.neighbours = array([[-1,-2,-3]])	                
+ 	domain.neighbours = array([[-1,-2,-3]])
         domain.quantities['stage'].boundary_values[0] = 1.0
-        
+
 	domain.distribute_to_vertices_and_edges() #Use first order default
-        
+
         domain.check_integrity()
 
@@ -113,5 +113,5 @@
         assert allclose(U, 0)
 
-        
+
 
 
@@ -122,11 +122,11 @@
 
         a = [0.0, 0.5]
-        b = [0.0, 0.0]                
+        b = [0.0, 0.0]
         c = [0.5, 0.5]
-        d = [0.5, 0.0]        
+        d = [0.5, 0.0]
 
         points = [a, b, c, d]
         vertices = [ [0,1,2], [3,2,1] ]
-        domain = Domain(points, vertices)        
+        domain = Domain(points, vertices)
         domain.check_integrity()
 
@@ -137,26 +137,26 @@
 	domain.distribute_to_vertices_and_edges()
 
-	domain.compute_fluxes()	        
+	domain.compute_fluxes()
 
         X = domain.quantities['stage'].explicit_update
         assert X[0] == -X[1]
-			
+
 
     def test_advection_example(self):
         #Test that system can evolve
-        
+
         from mesh_factory import rectangular
 
         points, vertices, boundary = rectangular(6, 6)
 
-        #Create advection domain with direction (1,-1)        
+        #Create advection domain with direction (1,-1)
         domain = Domain(points, vertices, boundary, velocity=[1.0, -1.0])
 
         # Initial condition is zero by default
-        
+
         #Boundaries
         T = Transmissive_boundary(domain)
         D = Dirichlet_boundary(array([3.1415]))
-        
+
         domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
         domain.check_integrity()
@@ -166,11 +166,11 @@
             if allclose(domain.quantities['stage'].centroid_values, 3.1415):
                 break
-            
+
         assert allclose(domain.quantities['stage'].centroid_values, 3.1415)
 
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Advection,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_all.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_all.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_all.py	(revision 1018)
@@ -23,5 +23,5 @@
 
     import sys
-    
+
     files = os.listdir(path)
 
@@ -38,10 +38,10 @@
         else:
             pass
-    return test_files    
+    return test_files
 
 
 
 def regressionTest():
-    import sys, os, re, unittest    
+    import sys, os, re, unittest
     path = os.path.split(sys.argv[0])[0] or os.getcwd()
 
@@ -53,17 +53,17 @@
     #test = re.compile('^test_[\w]*.py$', re.IGNORECASE)
     #files = filter(test.search, files)
-    
+
 
     try:
         files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
     except:
-        files.remove('test_all.py') 
+        files.remove('test_all.py')
 
 
-    if globals().has_key('exclude'):    
+    if globals().has_key('exclude'):
         for file in exclude:
             files.remove(file)
             print 'WARNING: File '+ file + ' excluded from testing'
-        
+
 
     filenameToModuleName = lambda f: os.path.splitext(f)[0]
@@ -76,5 +76,5 @@
 
     os.system('python compile.py') #Attempt to compile all extensions
-    
+
+    #print regressionTest()
     unittest.main(defaultTest='regressionTest')
-    
Index: /inundation/ga/storm_surge/pyvolution/test_cg_solve.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_cg_solve.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_cg_solve.py	(revision 1018)
@@ -9,14 +9,14 @@
 
 
-class TestCase(unittest.TestCase):
+class Test_CG_Solve(unittest.TestCase):
 
     def test_sparse_solve(self):
         """Small Sparse Matrix"""
-        
+
         A = [[2.0, -1.0, 0.0, 0.0 ],
              [-1.0, 2.0, -1.0, 0.0],
              [0.0, -1.0, 2.0, -1.0],
              [0.0,0.0, -1.0, 2.0]]
-        
+
         A = Sparse(A)
 
@@ -35,5 +35,5 @@
         n = 50
         A = Sparse(n,n)
-        
+
         for i in arange(0,n):
             A[i,i] = 1.0
@@ -42,5 +42,5 @@
             if i < n-1 :
                 A[i,i+1] = -0.5
-                
+
         xe = ones( (n,), Float)
 
@@ -52,5 +52,5 @@
     def test_solve_large_2d(self):
         """Standard 2d laplacian"""
-        
+
         n = 20
         m = 10
@@ -70,5 +70,5 @@
                 if j < m-1 :
                     A[I,I+1] = -1.0
-                
+
         xe = ones( (n*m,), Float)
 
@@ -81,5 +81,5 @@
         """Standard 2d laplacian with csr format
         """
-        
+
         n = 100
         m = 100
@@ -99,5 +99,5 @@
                 if j < m-1 :
                     A[I,I+1] = -1.0
-                
+
         xe = ones( (n*m,), Float)
 
@@ -114,5 +114,5 @@
     def test_solve_large_2d_with_default_guess(self):
         """Standard 2d laplacian using default first guess"""
-        
+
         n = 20
         m = 10
@@ -132,5 +132,5 @@
                 if j < m-1 :
                     A[I,I+1] = -1.0
-                
+
         xe = ones( (n*m,), Float)
 
@@ -143,10 +143,10 @@
     def test_vector_shape_error(self):
         """Raise VectorShapeError"""
-        
+
         A = [[2.0, -1.0, 0.0, 0.0 ],
              [-1.0, 2.0, -1.0, 0.0],
              [0.0, -1.0, 2.0, -1.0],
              [0.0,0.0, -1.0, 2.0]]
-        
+
         A = Sparse(A)
 
@@ -161,13 +161,11 @@
             raise msg
 
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-     suite = unittest.makeSuite(TestCase,'test')
+     suite = unittest.makeSuite(Test_CG_Solve,'test')
      runner = unittest.TextTestRunner() #(verbosity=2)
      runner.run(suite)
 
-    
-    
 
 
Index: /inundation/ga/storm_surge/pyvolution/test_data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1018)
@@ -11,5 +11,5 @@
 from config import epsilon
 
-class dataTestCase(unittest.TestCase):
+class Test_Data_Manager(unittest.TestCase):
     def setUp(self):
         import time
@@ -24,17 +24,17 @@
         domain.default_order=2
 
-        
+
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: -x)        
+        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
 
@@ -45,13 +45,13 @@
         h = 0.3
         for i in range(stage.shape[0]):
-            if i % 2 == 0:            
+            if i % 2 == 0:
                 stage[i,:] = bed[i,:] + h
             else:
                 stage[i,:] = bed[i,:]
-                
+
         domain.set_quantity('stage', stage)
         self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
 
-        domain.distribute_to_vertices_and_edges()    
+        domain.distribute_to_vertices_and_edges()
 
 
@@ -61,13 +61,13 @@
         self.X = C[:,0:6:2].copy()
         self.Y = C[:,1:6:2].copy()
-        
+
         self.F = bed
 
-        
+
     def tearDown(self):
         pass
 
 
-        
+
 
 #     def test_xya(self):
@@ -79,9 +79,9 @@
 
 #         domain = self.domain
-        
+
 #         domain.filename = 'datatest' + str(time.time())
 #         domain.format = 'xya'
 #         domain.smooth = True
-        
+
 #         xya = get_dataobject(self.domain)
 #         xya.store_all()
@@ -99,6 +99,6 @@
 #         if domain.smooth:
 #             self.failUnless(lFile[0] == '9 3 # <vertex #> <x> <y> [attributes]')
-#         else:    
-#             self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]')            
+#         else:
+#             self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]')
 
 #         #Get smoothed field values with X and Y
@@ -108,8 +108,8 @@
 
 #         Q,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
-#                                            indices = (0,), precision = Float)        
-
-
-        
+#                                            indices = (0,), precision = Float)
+
+
+
 #         for i, line in enumerate(lFile[1:]):
 #             fields = line.split()
@@ -121,5 +121,5 @@
 #             assert allclose(float(fields[2]), F[i,0])
 #             assert allclose(float(fields[3]), Q[i,0])
-#             assert allclose(float(fields[4]), F[i,1])            
+#             assert allclose(float(fields[4]), F[i,1])
 
 
@@ -133,10 +133,10 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
+        from Scientific.IO.NetCDF import NetCDFFile
 
         self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = False
-      
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -145,5 +145,5 @@
         #Get NetCDF
         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
-      
+
         # Get the variables
         x = fid.variables['x']
@@ -151,5 +151,5 @@
         z = fid.variables['elevation']
 
-        volumes = fid.variables['volumes']        
+        volumes = fid.variables['volumes']
 
 
@@ -164,9 +164,9 @@
             assert V[k, 0] == 3*k
             assert V[k, 1] == 3*k+1
-            assert V[k, 2] == 3*k+2            
-      
-      
-        fid.close()
-      
+            assert V[k, 2] == 3*k+2
+
+
+        fid.close()
+
         #Cleanup
         os.remove(sww.filename)
@@ -180,17 +180,17 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
-
-        self.domain.filename = 'datatest' + str(id(self))        
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = True
-      
+
         sww = get_dataobject(self.domain)
-        sww.store_connectivity()	
+        sww.store_connectivity()
 
         #Check contents
         #Get NetCDF
         fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
-      
+
         # Get the variables
         x = fid.variables['x']
@@ -198,9 +198,9 @@
         z = fid.variables['elevation']
 
-        volumes = fid.variables['volumes']        
+        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]))
@@ -209,5 +209,5 @@
         assert allclose([X[4], Y[4]], array([0.5, 0.5]))
 
-        assert allclose([X[7], Y[7]], array([1.0, 0.5]))        
+        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
 
         Z = z[:]
@@ -221,9 +221,9 @@
         assert V[4,0] == 6
         assert V[4,1] == 7
-        assert V[4,2] == 3                      
-      
-
-        fid.close()
-      
+        assert V[4,2] == 3
+
+
+        fid.close()
+
         #Cleanup
         os.remove(sww.filename)
@@ -237,11 +237,11 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
+        from Scientific.IO.NetCDF import NetCDFFile
 
         self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = True
-        self.domain.reduction = mean             
-      
+        self.domain.reduction = mean
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -251,63 +251,5 @@
         #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['elevation']        
-        time = fid.variables['time']
-        stage = fid.variables['stage']
-      
-
-	Q = self.domain.quantities['stage']	
-        Q0 = Q.vertex_values[:,0]
-        Q1 = Q.vertex_values[:,1]
-        Q2 = Q.vertex_values[:,2]
-
-        A = stage[0,:]
-        #print A[0], (Q2[0,0] + Q1[1,0])/2
-        assert allclose(A[0], (Q2[0] + Q1[1])/2)  
-        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
-        assert allclose(A[2], Q0[3])
-        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
-
-        #Center point
-        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
-                                 Q0[5] + Q2[6] + Q1[7])/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(id(self))        
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-
-        self.domain.reduction = mean     
-      
-        sww = get_dataobject(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep('stage')      
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep('stage')
-
-
-        #Check contents
-        #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+
 
         # Get the variables
@@ -316,14 +258,15 @@
         z = fid.variables['elevation']
         time = fid.variables['time']
-        stage = fid.variables['stage']        
-
-        #Check values
-	Q = self.domain.quantities['stage']	
+        stage = fid.variables['stage']
+
+
+	Q = self.domain.quantities['stage']
         Q0 = Q.vertex_values[:,0]
         Q1 = Q.vertex_values[:,1]
         Q2 = Q.vertex_values[:,2]
 
-        A = stage[1,:]
-        assert allclose(A[0], (Q2[0] + Q1[1])/2)  
+        A = stage[0,:]
+        #print A[0], (Q2[0,0] + Q1[1,0])/2
+        assert allclose(A[0], (Q2[0] + Q1[1])/2)
         assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
         assert allclose(A[2], Q0[3])
@@ -333,17 +276,16 @@
         assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
                                  Q0[5] + Q2[6] + Q1[7])/6)
-                               
-      
-      
-          
-      
-        fid.close()
-      
+
+
+
+        fid.close()
+
         #Cleanup
         os.remove(sww.filename)
-        
-    def test_sww_variable3(self):
+
+
+    def test_sww_variable2(self):
         """Test that sww information can be written correctly
-        multiple timesteps using a different reduction operator (min)
+        multiple timesteps. Use average as reduction operator
         """
 
@@ -352,22 +294,20 @@
         from Scientific.IO.NetCDF import NetCDFFile
 
-        self.domain.filename = 'datatest' + str(id(self))        
+        self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = True
-        self.domain.reduction = min
-      
+
+        self.domain.reduction = mean
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
         sww.store_timestep('stage')
-      
         self.domain.evolve_to_end(finaltime = 0.01)
         sww.store_timestep('stage')
 
 
-
         #Check contents
         #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')
-      
+        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
 
         # Get the variables
@@ -376,9 +316,8 @@
         z = fid.variables['elevation']
         time = fid.variables['time']
-        stage = fid.variables['stage']        
+        stage = fid.variables['stage']
 
         #Check values
-        #Check values
-	Q = self.domain.quantities['stage']	
+	Q = self.domain.quantities['stage']
         Q0 = Q.vertex_values[:,0]
         Q1 = Q.vertex_values[:,1]
@@ -386,5 +325,66 @@
 
         A = stage[1,:]
-        assert allclose(A[0], min(Q2[0], Q1[1]))  
+        assert allclose(A[0], (Q2[0] + Q1[1])/2)
+        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
+        assert allclose(A[2], Q0[3])
+        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
+
+        #Center point
+        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
+                                 Q0[5] + Q2[6] + Q1[7])/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(id(self))
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.reduction = min
+
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+
+
+        #Check contents
+        #Get NetCDF
+        fid = NetCDFFile(sww.filename, 'r')
+
+
+        # Get the variables
+        x = fid.variables['x']
+        y = fid.variables['y']
+        z = fid.variables['elevation']
+        time = fid.variables['time']
+        stage = fid.variables['stage']
+
+        #Check values
+        #Check values
+	Q = self.domain.quantities['stage']
+        Q0 = Q.vertex_values[:,0]
+        Q1 = Q.vertex_values[:,1]
+        Q2 = Q.vertex_values[:,2]
+
+        A = stage[1,:]
+        assert allclose(A[0], min(Q2[0], Q1[1]))
         assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
         assert allclose(A[2], Q0[3])
@@ -394,8 +394,8 @@
         assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
                                   Q0[5], Q2[6], Q1[7]))
-                               
-      
-        fid.close()
-      
+
+
+        fid.close()
+
         #Cleanup
         os.remove(sww.filename)
@@ -408,5 +408,5 @@
         import time, os, config
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
+        from Scientific.IO.NetCDF import NetCDFFile
 
         self.domain.filename = 'synctest'
@@ -426,8 +426,8 @@
             if t == 0.0:
                 assert allclose(stage, self.initial_stage)
-                assert allclose(stage_file[:], stage.flat) 
+                assert allclose(stage_file[:], stage.flat)
             else:
                 assert not allclose(stage, self.initial_stage)
-                assert not allclose(stage_file[:], stage.flat)                 
+                assert not allclose(stage_file[:], stage.flat)
 
             fid.close()
@@ -442,11 +442,11 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
-
-        self.domain.filename = 'datatest' + str(id(self))                
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = True
-        self.domain.reduction = mean             
-      
+        self.domain.reduction = mean
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -455,6 +455,6 @@
         #Check contents
         #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')  
-      
+        fid = NetCDFFile(sww.filename, 'r')
+
         # Get the variables
         x = fid.variables['x']
@@ -462,9 +462,9 @@
         z = fid.variables['elevation']
 
-        volumes = fid.variables['volumes']   
+        volumes = fid.variables['volumes']
         time = fid.variables['time']
 
-        # 2D 
-        stage = fid.variables['stage']       
+        # 2D
+        stage = fid.variables['stage']
 
         X = x[:]
@@ -474,5 +474,5 @@
         T = time[:]
         S = stage[:,:]
-      
+
 #         print "****************************"
 #         print "X ",X
@@ -489,7 +489,7 @@
 #         print "****************************"
 
-        
-        fid.close()
-      
+
+        fid.close()
+
         #Cleanup
         os.remove(sww.filename)
@@ -522,5 +522,5 @@
         ref_elevation = []
         for i in range(6):
-            y = (6-i)*25.0            
+            y = (6-i)*25.0
             for j in range(5):
                 x = j*25.0
@@ -530,5 +530,5 @@
                 ref_elevation.append(z)
                 fid.write('%f ' %z)
-            fid.write('\n')    
+            fid.write('\n')
 
         fid.close()
@@ -550,5 +550,5 @@
 """)
         fid.close()
-        
+
         #Convert to NetCDF xya
         convert_dem_from_ascii2netcdf(root)
@@ -558,5 +558,5 @@
         #Get NetCDF
         fid = NetCDFFile(root+'.pts', 'r')
-      
+
         # Get the variables
         #print fid.variables.keys()
@@ -572,14 +572,14 @@
         #print attributes[:]
         #print ref_elevation
-        assert allclose(elevation, ref_elevation)        
-      
-        #Cleanup
-        fid.close()
-        
+        assert allclose(elevation, ref_elevation)
+
+        #Cleanup
+        fid.close()
+
 
         os.remove(root + '.pts')
-        os.remove(root + '.dem')               
-        os.remove(root + '.asc')        
-        os.remove(root + '.prj')                
+        os.remove(root + '.dem')
+        os.remove(root + '.asc')
+        os.remove(root + '.prj')
 
 
@@ -590,5 +590,5 @@
         from Scientific.IO.NetCDF import NetCDFFile
 
-        #The test file has 
+        #The test file has
         # LON = 150.66667, 150.83334, 151, 151.16667
         # LAT = -34.5, -34.33333, -34.16667, -34 ;
@@ -597,5 +597,5 @@
         # First value (index=0) in small_ha.nc is 0.3400644 cm,
         # Fourth value (index==3) is -6.50198 cm
-        
+
 
         from coordinate_transforms.redfearn import redfearn
@@ -605,19 +605,19 @@
         first_value = fid.variables['HA'][:][0,0,0]
         fourth_value = fid.variables['HA'][:][0,0,3]
-        
+
         #Call conversion (with zero origin)
         ferret2sww('small', verbose=False,
                    origin = (56, 0, 0))
- 
+
 
         #Work out the UTM coordinates for first point
         zone, e, n = redfearn(-34.5, 150.66667)
         #print zone, e, n
-        
+
         #Read output file 'small.sww'
         fid = NetCDFFile('small.sww')
 
         x = fid.variables['x'][:]
-        y = fid.variables['y'][:]        
+        y = fid.variables['y'][:]
 
         #Check that first coordinate is correctly represented
@@ -631,5 +631,5 @@
 
         #Check fourth value
-        assert allclose(stage[0,3], fourth_value/100)  #Meters        
+        assert allclose(stage[0,3], fourth_value/100)  #Meters
 
         fid.close()
@@ -639,5 +639,5 @@
         os.remove('small.sww')
 
-        
+
 
     def test_ferret2sww_2(self):
@@ -647,5 +647,5 @@
         from Scientific.IO.NetCDF import NetCDFFile
 
-        #The test file has 
+        #The test file has
         # LON = 150.66667, 150.83334, 151, 151.16667
         # LAT = -34.5, -34.33333, -34.16667, -34 ;
@@ -654,5 +654,5 @@
         # First value (index=0) in small_ha.nc is 0.3400644 cm,
         # Fourth value (index==3) is -6.50198 cm
-        
+
 
         from coordinate_transforms.redfearn import redfearn
@@ -665,26 +665,26 @@
         lat_index = 0
         lon_index = 2
-        
+
         test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
         test_time = fid.variables['TIME'][:][time_index]
         test_lat = fid.variables['LAT'][:][lat_index]
-        test_lon = fid.variables['LON'][:][lon_index]        
-        
-        linear_point_index = lat_index*4 + lon_index 
-        fid.close()
-        
+        test_lon = fid.variables['LON'][:][lon_index]
+
+        linear_point_index = lat_index*4 + lon_index
+        fid.close()
+
         #Call conversion (with zero origin)
         ferret2sww('small', verbose=False,
                    origin = (56, 0, 0))
- 
+
 
         #Work out the UTM coordinates for test point
         zone, e, n = redfearn(test_lat, test_lon)
-        
+
         #Read output file 'small.sww'
         fid = NetCDFFile('small.sww')
 
         x = fid.variables['x'][:]
-        y = fid.variables['y'][:]        
+        y = fid.variables['y'][:]
 
         #Check that test coordinate is correctly represented
@@ -705,5 +705,5 @@
 
 
-        
+
     def test_sww_extent(self):
         """Not a test, rather a look at the sww format
@@ -712,13 +712,13 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
-
-        self.domain.filename = 'datatest' + str(id(self))                        
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        self.domain.filename = 'datatest' + str(id(self))
         self.domain.format = 'sww'
         self.domain.smooth = True
-        self.domain.reduction = mean 
+        self.domain.reduction = mean
         self.domain.set_datadir('.')
 
-      
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -736,10 +736,10 @@
         [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
                extent_sww(file_and_extension_name )
-        
-        assert allclose(xmin, 0.0) 
-        assert allclose(xmax, 1.0) 
-        assert allclose(ymin, 0.0) 
+
+        assert allclose(xmin, 0.0)
+        assert allclose(xmax, 1.0)
+        assert allclose(ymin, 0.0)
         assert allclose(ymax, 1.0)
-        assert allclose(stagemin, -0.85) 
+        assert allclose(stagemin, -0.85)
         assert allclose(stagemax, 0.15)
 
@@ -747,9 +747,9 @@
         #Cleanup
         os.remove(sww.filename)
-        
-        
+
+
     def test_ferret2sww_nz_origin(self):
         from Scientific.IO.NetCDF import NetCDFFile
-        from coordinate_transforms.redfearn import redfearn        
+        from coordinate_transforms.redfearn import redfearn
 
         #Call conversion (with nonzero origin)
@@ -765,5 +765,5 @@
 
         x = fid.variables['x'][:]
-        y = fid.variables['y'][:]        
+        y = fid.variables['y'][:]
 
         #Check that first coordinate is correctly represented
@@ -777,8 +777,8 @@
         os.remove('small.sww')
 
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(dataTestCase,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 1018)
@@ -18,27 +18,27 @@
         #print 'bottom - indexes',elements
         domain.set_quantity('friction', 0.09, indexes = elements)
- 
+
 def set_top_friction(tag, elements, domain):
     if tag == "top":
         #print 'top - indexes',elements
         domain.set_quantity('friction', 1., indexes = elements)
-        
- 
+
+
 def set_all_friction(tag, elements, domain):
     if tag == "all":
         new_values = domain.get_quantity('friction', indexes = elements) + 10.0
-        
+
         domain.set_quantity('friction', new_values, indexes = elements)
-                  
-      
-class TestCase(unittest.TestCase):
+
+
+class Test_Domain(unittest.TestCase):
     def setUp(self):
         pass
 
-        
+
     def tearDown(self):
         pass
 
-            
+
     def test_simple(self):
         a = [0.0, 0.0]
@@ -50,12 +50,12 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
 
         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
         other_quantities = ['elevation', 'friction']
-        
+
         domain = Domain(points, vertices, None,
-                        conserved_quantities, other_quantities)        
+                        conserved_quantities, other_quantities)
         domain.check_integrity()
 
@@ -64,18 +64,18 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.    
+        assert domain.get_conserved_quantities(0, edge=1) == 0.
 
 
     def test_conserved_quantities(self):
-        
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0,0.0]
-
-        points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
 
@@ -87,9 +87,9 @@
         domain.set_quantity('stage', [[1,2,3], [5,5,5],
                                       [0,0,9], [-6, 3, 3],
-                                      [0,0,0], [0,0,0]])        
+                                      [0,0,0], [0,0,0]])
 
         domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
                                           [0,0,9], [-6, 3, 3],
-                                          [0,0,0], [0,0,0]])        
+                                          [0,0,0], [0,0,0]])
 
         domain.check_integrity()
@@ -100,5 +100,5 @@
 
         q = domain.get_conserved_quantities(1)
-        assert allclose(q, [5., 5., 0.])         
+        assert allclose(q, [5., 5., 0.])
 
         q = domain.get_conserved_quantities(2)
@@ -106,6 +106,6 @@
 
         q = domain.get_conserved_quantities(3)
-        assert allclose(q, [0., 0., 0.])                         
-        
+        assert allclose(q, [0., 0., 0.])
+
 
         #Edges
@@ -135,5 +135,5 @@
         assert allclose(q, [-1.5, -1.5, 0.])
         q = domain.get_conserved_quantities(3, edge=2)
-        assert allclose(q, [-1.5, -1.5, 0.])                    
+        assert allclose(q, [-1.5, -1.5, 0.])
 
 
@@ -142,7 +142,7 @@
         from config import default_boundary_tag
 
-        
+
         a = [0.0, 0.5]
-        b = [0.0, 0.0]                
+        b = [0.0, 0.0]
         c = [0.5, 0.5]
 
@@ -153,13 +153,13 @@
         domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} )
 
-        
-        domain.check_integrity()
-
- 	assert allclose(domain.neighbours, [[-1,-2,-3]])	
-		
+
+        domain.check_integrity()
+
+ 	assert allclose(domain.neighbours, [[-1,-2,-3]])
+
 
 
     def test_boundary_conditions(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -177,10 +177,10 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Second'}                                         
-                     
+                     (3, 2): 'Second'}
+
 
         domain = Domain(points, vertices, boundary,
                         conserved_quantities =\
-                        ['stage', 'xmomentum', 'ymomentum'])        
+                        ['stage', 'xmomentum', 'ymomentum'])
         domain.check_integrity()
 
@@ -205,5 +205,5 @@
                domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
         assert domain.quantities['stage'].boundary_values[5] ==\
-               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)         
+               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)
 
         #Check enumeration
@@ -217,5 +217,5 @@
         """Domain implements a default first order gradient limiter
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -233,10 +233,10 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Third'}                                         
-                     
+                     (3, 2): 'Third'}
+
 
         domain = Domain(points, vertices, boundary,
                         conserved_quantities =\
-                        ['stage', 'xmomentum', 'ymomentum'])                
+                        ['stage', 'xmomentum', 'ymomentum'])
         domain.check_integrity()
 
@@ -247,10 +247,10 @@
         assert allclose( domain.quantities['stage'].centroid_values,
                          [2,5,3,0] )
-                         
+
         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
                                           [3,3,3], [4, 4, 4]])
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -283,10 +283,10 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Third'}                                         
-                     
+                     (3, 2): 'Third'}
+
 
         domain = Domain(points, vertices, boundary,
                         conserved_quantities =\
-                        ['stage', 'xmomentum', 'ymomentum'])                
+                        ['stage', 'xmomentum', 'ymomentum'])
         domain.check_integrity()
 
@@ -295,5 +295,5 @@
         domain.set_quantity('xmomentum', [1,2,3,4], 'centroids')
         domain.set_quantity('ymomentum', [1,2,3,4], 'centroids')
-                                          
+
 
         #Assign some values to update vectors
@@ -303,5 +303,5 @@
             domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
             domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
-            
+
 
         #Update with given timestep (assuming no other forcing terms)
@@ -316,5 +316,5 @@
 
         for name in domain.conserved_quantities:
-            assert allclose(domain.quantities[name].centroid_values, x)   
+            assert allclose(domain.quantities[name].centroid_values, x)
 
 
@@ -322,5 +322,5 @@
         """Set quantities for sub region
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -338,9 +338,9 @@
                      (2, 1): 'Second',
                      (3, 1): 'Second',
-                     (3, 2): 'Third'}             
+                     (3, 2): 'Third'}
 
         domain = Domain(points, vertices, boundary,
                         conserved_quantities =\
-                        ['stage', 'xmomentum', 'ymomentum'])                
+                        ['stage', 'xmomentum', 'ymomentum'])
         domain.check_integrity()
 
@@ -350,10 +350,10 @@
         assert allclose( domain.quantities['stage'].centroid_values,
                          [2,5,3,0] )
-                         
+
         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
                                           [3,3,3], [4, 4, 4]])
 
         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
+                                          [30,30,30], [40, 40, 40]])
 
 
@@ -369,8 +369,8 @@
         domain.build_tagged_elements_dictionary({'mound':[0,1]})
         domain.set_region([add_to_verts])
-        
+
         self.failUnless(domain.test == "Mound",
                         'set region failed')
-        
+
 
 
@@ -382,5 +382,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -392,5 +392,5 @@
                                                  'all':[0,1,2,3,4,5]})
 
-        
+
         #Set friction
         manning = 0.07
@@ -406,5 +406,5 @@
                          [ 1.0,  1.0,  1.0],
                          [ 1.0,  1.0,  1.0]])
-       
+
         domain.set_region([set_all_friction])
         #print domain.quantities['friction'].get_values()
@@ -416,7 +416,9 @@
                          [ 11.0,  11.0,  11.0],
                          [ 11.0,  11.0,  11.0]])
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Domain,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_general_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_general_mesh.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_general_mesh.py	(revision 1018)
@@ -9,9 +9,9 @@
 from Numeric import allclose, array, ones, Float
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_General_Mesh(unittest.TestCase):
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
@@ -24,12 +24,12 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
         domain = Domain(points, vertices, boundary)
-        
+
         value = [7]
         indexes = [1]
-        assert  domain.get_vertices() == domain.triangles 
+        assert  domain.get_vertices() == domain.triangles
         assert domain.get_vertices([0,4]) == [domain.triangles[0],
                                               domain.triangles[4]]
@@ -38,10 +38,10 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
         domain = Domain(points, vertices, boundary)
 
-	assert domain.get_area() == 1.0	        
+	assert domain.get_area() == 1.0
 
 
@@ -53,23 +53,22 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
         domain = Domain(points, vertices, boundary)
-        
+
         assert  domain.get_unique_vertices() == [0,1,2,3,4,5,6,7]
         unique_vertices = domain.get_unique_vertices([0,1,4])
         unique_vertices.sort()
         assert unique_vertices == [0,1,2,4,5,6,7]
-        
+
         unique_vertices = domain.get_unique_vertices([0,4])
         unique_vertices.sort()
         assert unique_vertices == [0,2,4,5,6,7]
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_General_Mesh,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
 
-
Index: /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 1018)
@@ -8,10 +8,10 @@
 from Numeric import allclose, array
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_Generic_Boundary_Conditions(unittest.TestCase):
     def setUp(self):
         pass
         #print "  Setting up"
-        
+
     def tearDown(self):
         pass
@@ -28,5 +28,5 @@
         else:
             raise 'Should have raised exception'
-        
+
 
     def test_dirichlet_empty(self):
@@ -50,5 +50,5 @@
         from domain import Domain
         from quantity import Quantity
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -62,9 +62,9 @@
         #bac, bce, ecf, dbe
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-        
+
         domain = Domain(points, elements)
         domain.check_integrity()
 
-        domain.conserved_quantities = ['stage', 'ymomentum']        
+        domain.conserved_quantities = ['stage', 'ymomentum']
         domain.quantities['stage'] =\
                                    Quantity(domain, [[1,2,3], [5,5,5],
@@ -110,5 +110,5 @@
         from math import sin, pi
         from config import time_format
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -122,5 +122,5 @@
         #bac, bce, ecf, dbe
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-        
+
         domain = Domain(points, elements)
         domain.conserved_quantities = ['stage', 'ymomentum']
@@ -134,5 +134,5 @@
 
         domain.check_integrity()
-        
+
 
         #Write file
@@ -145,17 +145,17 @@
             t_string = time.strftime(time_format, time.gmtime(t))
 
-            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 
+            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10)))
         fid.close()
-                      
-        
+
+
         F = File_boundary(filename, domain)
         os.remove(filename)
-        
+
 
         #Check that midpoint coordinates at boundary are correctly computed
-	assert allclose( F.midpoint_coordinates, 
-	                 [[0.0, 3.0], [1.0, 3.0], [0.0, 1.0], 
+	assert allclose( F.midpoint_coordinates,
+	                 [[0.0, 3.0], [1.0, 3.0], [0.0, 1.0],
 			 [1.0, 0.0], [3.0, 0.0], [3.0, 1.0]])
-	
+
         #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0])
         #assert allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0])
@@ -163,5 +163,5 @@
         #assert allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0])
         #assert allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0])
-        #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0])          
+        #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0])
 
 
@@ -177,5 +177,5 @@
         domain.time = 2.5*5*60  #Half way between steps 2 and 3
         q = F.evaluate()
-        assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])        
+        assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])
 
 
@@ -191,5 +191,5 @@
         from math import sin, pi
         from config import time_format
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -203,5 +203,5 @@
         #bac, bce, ecf, dbe
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-        
+
         domain = Domain(points, elements)
         domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
@@ -228,7 +228,7 @@
             t_string = time.strftime(time_format, time.gmtime(t))
 
-            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 
+            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10)))
         fid.close()
-                      
+
 
         try:
@@ -242,11 +242,11 @@
 
 
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
 
 
-        
+
Index: /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py	(revision 1018)
@@ -16,5 +16,5 @@
 
 
-class dataTestCase(unittest.TestCase):
+class Test_Interpolate_sww(unittest.TestCase):
     def setUp(self):
         import time
@@ -30,17 +30,17 @@
         domain.beta_h = 0
 
-        
+
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: -x)        
+        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
 
@@ -50,12 +50,12 @@
         h = 0.3
         for i in range(stage.shape[0]):
-            if i % 2 == 0:            
+            if i % 2 == 0:
                 stage[i,:] = bed[i,:] + h
             else:
                 stage[i,:] = bed[i,:]
-                
+
         domain.set_quantity('stage', stage)
 
-        domain.distribute_to_vertices_and_edges()    
+        domain.distribute_to_vertices_and_edges()
 
 
@@ -65,8 +65,8 @@
         self.X = C[:,0:6:2].copy()
         self.Y = C[:,1:6:2].copy()
-        
+
         self.F = bed
 
-        
+
     def tearDown(self):
         pass
@@ -78,21 +78,21 @@
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
+        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             
-      
+        self.domain.reduction = mean
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
         sww.store_timestep('stage')
-        self.domain.time = 2. 
+        self.domain.time = 2.
         sww.store_timestep('stage')
 
         #Check contents
         #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')  
-      
+        fid = NetCDFFile(sww.filename, 'r')
+
         # Get the variables
         x = fid.variables['x']
@@ -100,9 +100,9 @@
         z = fid.variables['elevation']
 
-        volumes = fid.variables['volumes']   
+        volumes = fid.variables['volumes']
         time = fid.variables['time']
 
-        # 2D 
-        stage = fid.variables['stage']       
+        # 2D
+        stage = fid.variables['stage']
 
         X = x[:]
@@ -127,15 +127,15 @@
             print "Stage ",S
             print "****************************"
-            
-
-        
-        
+
+
+
+
         fid.close()
-      
+
         #Cleanup
         os.remove(sww.filename)
- 
+
     def test_interpolate_sww(self):
-        """Not reaa unit test, rather a system test for 
+        """Not reaa unit test, rather a system test for
         """
 
@@ -143,5 +143,5 @@
         from Numeric import array, zeros, allclose, Float, concatenate, \
              transpose
-        from Scientific.IO.NetCDF import NetCDFFile    
+        from Scientific.IO.NetCDF import NetCDFFile
         import tempfile
         from load_mesh.loadASCII import  load_points_file, \
@@ -151,15 +151,15 @@
         self.domain.format = 'sww'
         self.domain.smooth = True
-        self.domain.reduction = mean             
-      
+        self.domain.reduction = mean
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
         sww.store_timestep('stage')
-        self.domain.time = 2. 
+        self.domain.time = 2.
         sww.store_timestep('stage')
 
         #print "self.domain.filename",self.domain.filename
         interp = Interpolate_sww(sww.filename, 'depth')
-        
+
         assert allclose(interp.time,[0.0,2.0])
 
@@ -186,9 +186,9 @@
         point_file_out = tempfile.mktemp(".xya")
         interp.write_depth_xya(point_file_out)
-        
+
         #check the output file
         xya_dict = load_points_file(point_file_out)
         #print "xya_dict", xya_dict
-        title_list,point_attributes = concatinate_attributelist(xya_dict['attributelist'])  
+        title_list,point_attributes = concatinate_attributelist(xya_dict['attributelist'])
         assert allclose(interp.point_coordinates, xya_dict['pointlist'])
         assert allclose(interp.interpolated_quantity,
@@ -201,5 +201,5 @@
         #time_list = []
 
-        # Try another quantity      
+        # Try another quantity
         interp = Interpolate_sww(sww.filename, 'stage')
         interp.interpolate_xya(point_file)
@@ -209,5 +209,5 @@
         assert allclose(interp.interpolated_quantity,answer)
 
-        
+
         # look at error catching
         try:
@@ -218,5 +218,5 @@
             self.failUnless(0==1,
                         'bad key did not raise an error!')
-        
+
         # look at error catching
         try:
@@ -227,13 +227,13 @@
             self.failUnless(0==1,
                         'bad key did not raise an error!')
-                
+
         #Cleanup
         os.remove(sww.filename)
         os.remove(point_file_out)
         os.remove(point_file)
- 
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(dataTestCase,'test')
+    suite = unittest.makeSuite(Test_Interpolate_sww,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_least_squares.py	(revision 1018)
@@ -11,16 +11,16 @@
 
 def distance(x, y):
-    return sqrt( sum( (array(x)-array(y))**2 ))		
+    return sqrt( sum( (array(x)-array(y))**2 ))
 
 def linear_function(point):
     point = array(point)
     return point[:,0]+point[:,1]
-    
-        
-class TestCase(unittest.TestCase):
+
+
+class Test_Least_Squares(unittest.TestCase):
 
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
@@ -33,9 +33,9 @@
         vertices = [ [1,0,2] ]   #bac
 
-        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point        
-        
+        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
+
         interp = Interpolation(points, vertices, data)
-        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 
-        
+        assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
+
 
     def test_quad_tree(self):
@@ -50,5 +50,5 @@
         p8 = [-66.0, 20.0]
         p9 = [10.0, -66.0]
-        
+
         points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
         triangles = [ [0, 1, 2],
@@ -60,7 +60,7 @@
                       [7, 4, 5],
                       [8, 0, 2],
-                      [0, 9, 1]]   
-
-        data = [ [4,4] ] 
+                      [0, 9, 1]]
+
+        data = [ [4,4] ]
         interp = Interpolation(points, triangles, data, alpha = 0.0,
                                max_points_per_cell = 4)
@@ -75,5 +75,5 @@
         assert allclose(interp.get_A(), answer)
 
-        
+
         #point outside of quad tree root cell
         interp.set_point_coordinates([[-70, -70]])
@@ -82,5 +82,5 @@
                       0., 0. , 0., 0., 0., 0.]]
         assert allclose(interp.get_A(), answer)
-     
+
     def test_expand_search(self):
         p0 = [-10.0, -10.0]
@@ -94,5 +94,5 @@
         p8 = [-66.0, 20.0]
         p9 = [10.0, -66.0]
-        
+
         points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
         triangles = [ [0, 1, 2],
@@ -104,5 +104,5 @@
                       [7, 4, 5],
                       [8, 0, 2],
-                      [0, 9, 1]]   
+                      [0, 9, 1]]
 
         data = [ [4,4],
@@ -134,5 +134,5 @@
                  [50,50],
                  [30,0],
-                 [-20,-20]] 
+                 [-20,-20]]
         point_attributes = [ -400000,
                      10,
@@ -163,6 +163,6 @@
                      10,
                    10,
-                   99] 
-        
+                   99]
+
         interp = Interpolation(points, triangles, data,
                                alpha=0.0, expand_search=False, #verbose = True, #False,
@@ -170,5 +170,5 @@
         calc = interp.fit_points(point_attributes, )
         #print "calc",calc
-        
+
         # the point at 4,4 is ignored.  An expanded search has to be done
         # to fine which triangel it's in.
@@ -194,5 +194,5 @@
         p4 = [10.0, 60.0]
         p5 = [60.0, 60.0]
-        
+
         points = [p0, p1, p2, p3, p4, p5]
         triangles = [ [0, 1, 2],
@@ -201,5 +201,5 @@
                       [0, 2, 4],
                       [4, 2, 5],
-                      [5, 2, 3]]   
+                      [5, 2, 3]]
 
         data = [ [-26.0,-26.0] ]
@@ -216,5 +216,5 @@
         assert allclose(interp.get_A(), answer)
 
-        
+
         #point outside of quad tree root cell
         interp.set_point_coordinates([[-70, -70]])
@@ -223,10 +223,10 @@
                       0., 0. ]]
         assert allclose(interp.get_A(), answer)
-           
+
 
     def test_datapoints_at_vertices(self):
         """Test that data points coinciding with vertices yield a diagonal matrix
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -235,11 +235,11 @@
         vertices = [ [1,0,2] ]   #bac
 
-        data = points #Use data at vertices        
-        
+        data = points #Use data at vertices
+
         interp = Interpolation(points, vertices, data)
         assert allclose(interp.get_A(), [[1., 0., 0.],
                                    [0., 1., 0.],
-                                   [0., 0., 1.]]) 
-        
+                                   [0., 0., 1.]])
+
 
 
@@ -248,5 +248,5 @@
         each point should affect two matrix entries equally
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -256,8 +256,8 @@
 
         data = [ [0., 1.], [1., 0.], [1., 1.] ]
-        
+
         interp = Interpolation(points, vertices, data)
-        
-        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0   
+
+        assert allclose(interp.get_A(), [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0
                                    [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
                                    [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
@@ -268,5 +268,5 @@
         each point should affect two matrix entries in proportion
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -276,7 +276,7 @@
 
         data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
-        
+
         interp = Interpolation(points, vertices, data)
-        
+
         assert allclose(interp.get_A(), [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0
                                    [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
@@ -288,5 +288,5 @@
 
         from Numeric import sum
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -296,9 +296,9 @@
 
         data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
-        
+
         interp = Interpolation(points, vertices, data)
         #print "interp.get_A()", interp.get_A()
         assert allclose(sum(interp.get_A(), axis=1), 1.0)
-        
+
     def test_arbitrary_datapoints_some_outside(self):
         """Try arbitrary datapoints one outside the triangle.
@@ -307,5 +307,5 @@
 
         from Numeric import sum
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -322,10 +322,10 @@
         interp = Interpolation(points, vertices, data, precrop = False)
         assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
-        
-        
-
-    # this causes a memory error in scipy.sparse       
+
+
+
+    # this causes a memory error in scipy.sparse
     def test_more_triangles(self):
-        
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -340,7 +340,7 @@
         #Data points
         data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
-        interp = Interpolation(points, triangles, data_points) 
-
-        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d     
+        interp = Interpolation(points, triangles, data_points)
+
+        answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],    #Affects point d
                   [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],    #Affects points a and d
                   [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
@@ -354,5 +354,5 @@
             for j in range(A.shape[1]):
                 if not allclose(A[i,j], answer[i][j]):
-                    print i,j,':',A[i,j], answer[i][j] 
+                    print i,j,':',A[i,j], answer[i][j]
 
 
@@ -368,5 +368,5 @@
         points = [a, b, c]
         triangles = [ [1,0,2] ]   #bac
-         
+
         d1 = [1.0, 1.0]
         d2 = [1.0, 3.0]
@@ -375,6 +375,6 @@
         z2 = 4
         z3 = 4
-        data_coords = [d1, d2, d3] 
-        
+        data_coords = [d1, d2, d3]
+
         interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
         z = [z1, z2, z3]
@@ -384,10 +384,10 @@
         #print "f\n",f
         #print "answer\n",answer
-        
+
         assert allclose(f, answer, atol=1e-7)
-        
-        
+
+
     def test_smooth_att_to_meshII(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -395,18 +395,18 @@
         points = [a, b, c]
         triangles = [ [1,0,2] ]   #bac
-         
+
         d1 = [1.0, 1.0]
         d2 = [1.0, 2.0]
         d3 = [3.0,1.0]
-        data_coords = [d1, d2, d3]       
+        data_coords = [d1, d2, d3]
         z = linear_function(data_coords)
         interp = Interpolation(points, triangles, data_coords, alpha=0.0)
         f = interp.fit(z)
         answer = linear_function(points)
-        
+
         assert allclose(f, answer)
-    
+
     def test_smooth_attributes_to_meshIII(self):
-        
+
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -441,10 +441,10 @@
         #print "answer\n",answer
         assert allclose(f, answer)
-        
+
 
     def test_smooth_attributes_to_meshIV(self):
         """ Testing 2 attributes smoothed to the mesh
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -452,5 +452,5 @@
         points = [a, b, c]
         triangles = [ [1,0,2] ]   #bac
-         
+
         d1 = [1.0, 1.0]
         d2 = [1.0, 3.0]
@@ -459,6 +459,6 @@
         z2 = [4, 8]
         z3 = [4, 8]
-        data_coords = [d1, d2, d3] 
-        
+        data_coords = [d1, d2, d3]
+
         interp = Interpolation(points, triangles, data_coords, alpha=0.0)
         z = [z1, z2, z3]
@@ -466,27 +466,27 @@
         answer = [[0,0], [5., 10.], [5., 10.]]
         assert allclose(f, answer)
-        
+
     def test_interpolate_attributes_to_points(self):
         v0 = [0.0, 0.0]
         v1 = [0.0, 5.0]
         v2 = [5.0, 0.0]
-        
+
         vertices = [v0, v1, v2]
         triangles = [ [1,0,2] ]   #bac
-         
+
         d0 = [1.0, 1.0]
         d1 = [1.0, 2.0]
         d2 = [3.0, 1.0]
         point_coords = [ d0, d1, d2]
-        
+
         interp = Interpolation(vertices, triangles, point_coords)
         f = linear_function(vertices)
         z = interp.interpolate(f)
         answer = linear_function(point_coords)
-        
+
 
         assert allclose(z, answer)
 
-       
+
     def test_interpolate_attributes_to_pointsII(self):
         a = [-1.0, 0.0]
@@ -500,5 +500,5 @@
         triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
-         
+
         point_coords = [[-2.0, 2.0],
                         [-1.0, 1.0],
@@ -513,5 +513,5 @@
                         [0.5, -1.9],
                         [3.0, 1.0]]
-        
+
         interp = Interpolation(vertices, triangles, point_coords)
         f = linear_function(vertices)
@@ -521,5 +521,5 @@
         #print "answer",answer
         assert allclose(z, answer)
-    
+
     def test_interpolate_attributes_to_pointsIII(self):
         """Test linear interpolation of known values at vertices to
@@ -530,5 +530,5 @@
         c = [5.0, 0.0]
         d = [5.0, 5.0]
-        
+
         vertices = [a, b, c, d]
         triangles = [ [1,0,2], [2,3,0] ]   #bac, cdb
@@ -545,20 +545,20 @@
         d4 = [2.5, 2.5]
         d5 = [4.0, 1.0]
-         
+
         #Point on common vertex
         d6 = [0., 5.]
 
-        
+
         point_coords = [d0, d1, d2, d3, d4, d5, d6]
-        
+
         interp = Interpolation(vertices, triangles, point_coords)
 
         #Known values at vertices
-        #Functions are x+y, x+2y, 2x+y, x-y-5        
-        f = [ [0., 0., 0., -5.],        # (0,0)    
+        #Functions are x+y, x+2y, 2x+y, x-y-5
+        f = [ [0., 0., 0., -5.],        # (0,0)
               [5., 10., 5., -10.],      # (0,5)
               [5., 5., 10.0, 0.],       # (5,0)
               [10., 15., 15., -5.]]     # (5,5)
-         
+
         z = interp.interpolate(f)
         answer = [ [2., 3., 3., -5.],   # (1,1)
@@ -576,5 +576,5 @@
 
         assert allclose(z, answer)
-   
+
     def test_interpolate_attributes_to_pointsIV(self):
         a = [-1.0, 0.0]
@@ -588,5 +588,5 @@
         triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
-         
+
         point_coords = [[-2.0, 2.0],
                         [-1.0, 1.0],
@@ -601,5 +601,5 @@
                         [0.5, -1.9],
                         [3.0, 1.0]]
-        
+
         interp = Interpolation(vertices, triangles, point_coords)
         f = array([linear_function(vertices),2*linear_function(vertices) ])
@@ -607,5 +607,5 @@
         #print "f",f
         z = interp.interpolate(f)
-        answer = [linear_function(point_coords),             
+        answer = [linear_function(point_coords),
                   2*linear_function(point_coords) ]
         answer = transpose(answer)
@@ -613,9 +613,9 @@
         #print "answer",answer
         assert allclose(z, answer)
-        
+
     def test_smooth_attributes_to_mesh_function(self):
         """ Testing 2 attributes smoothed to the mesh
         """
-        
+
         a = [0.0, 0.0]
         b = [0.0, 5.0]
@@ -623,5 +623,5 @@
         points = [a, b, c]
         triangles = [ [1,0,2] ]   #bac
-         
+
         d1 = [1.0, 1.0]
         d2 = [1.0, 3.0]
@@ -632,12 +632,12 @@
         data_coords = [d1, d2, d3]
         z = [z1, z2, z3]
-        
+
         f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
         answer = [[0, 0], [5., 10.], [5., 10.]]
-        
+
         assert allclose(f, answer)
 
 
-	
+
     def test_pts2rectangular(self):
 
@@ -645,18 +645,18 @@
         FN = 'xyatest' + str(time.time()) + '.xya'
     	fid = open(FN, 'w')
-	fid.write(' %s \n' %('elevation'))        
+	fid.write(' %s \n' %('elevation'))
 	fid.write('%f %f %f\n' %(1,1,2) )
 	fid.write('%f %f %f\n' %(1,3,4) )
-	fid.write('%f %f %f\n' %(3,1,4) )		
+	fid.write('%f %f %f\n' %(3,1,4) )
 	fid.close()
-	
+
 	points, triangles, boundary, attributes =\
                 pts2rectangular(FN, 4, 4, format = 'asc')
-	
+
 
         data_coords = [ [1,1], [1,3], [3,1] ]
         z = [2, 4, 4]
 
-        ref = fit_to_mesh(points, triangles, data_coords, z)		
+        ref = fit_to_mesh(points, triangles, data_coords, z)
 
         #print attributes
@@ -665,7 +665,7 @@
 
         os.remove(FN)
-	
-
-    #Tests of smoothing matrix    
+
+
+    #Tests of smoothing matrix
     def test_smoothing_matrix_one_triangle(self):
         from Numeric import dot
@@ -674,5 +674,5 @@
         c = [2.0,0.0]
         points = [a, b, c]
-        
+
         vertices = [ [1,0,2] ]   #bac
 
@@ -689,5 +689,5 @@
         #           int 1 dx dy = area = 2
         assert dot(dot(f, interp.get_D()), f) == 2
-        
+
         #Define f(x,y) = y
         f = array([0,2,0])  #Value at global vertex 1
@@ -702,11 +702,11 @@
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
         #           int 2 dx dy = 2*area = 4
-        assert dot(dot(f, interp.get_D()), f) == 4        
-        
+        assert dot(dot(f, interp.get_D()), f) == 4
+
 
 
     def test_smoothing_matrix_more_triangles(self):
         from Numeric import dot
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -717,5 +717,5 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
@@ -728,12 +728,12 @@
 
         #Define f(x,y) = x
-        f = array([0,0,2,0,2,4]) #f evaluated at points a-f 
+        f = array([0,0,2,0,2,4]) #f evaluated at points a-f
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
         #           int 1 dx dy = total area = 8
         assert dot(dot(f, interp.get_D()), f) == 8
-        
+
         #Define f(x,y) = y
-        f = array([0,2,0,4,2,0]) #f evaluated at points a-f 
+        f = array([0,2,0,4,2,0]) #f evaluated at points a-f
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -742,14 +742,14 @@
 
         #Define f(x,y) = x+y
-        f = array([0,2,2,4,4,4])  #f evaluated at points a-f   
+        f = array([0,2,2,4,4,4])  #f evaluated at points a-f
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
         #           int 2 dx dy = 2*area = 16
-        assert dot(dot(f, interp.get_D()), f) == 16        
+        assert dot(dot(f, interp.get_D()), f) == 16
 
 
     def test_fit_and_interpolation(self):
         from mesh import Mesh
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -760,5 +760,5 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
@@ -773,9 +773,9 @@
                        [ 1.0, 1.0],
                        [ 1.0, 2.0],
-                       [ 1.0, 3.0],                                          
+                       [ 1.0, 3.0],
                        [ 2.0, 1.0],
                        [ 3.0, 0.0],
-                       [ 3.0, 1.0]]                                          
-        
+                       [ 3.0, 1.0]]
+
         interp = Interpolation(points, triangles, data_points, alpha=0.0)
 
@@ -797,5 +797,5 @@
 
     def test_smoothing_and_interpolation(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -806,5 +806,5 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
@@ -816,9 +816,9 @@
 
         z = linear_function(data_points)
-        answer = linear_function(points)        
+        answer = linear_function(points)
 
         #Make interpolator with too few data points and no smoothing
         interp = Interpolation(points, triangles, data_points, alpha=0.0)
-        #Must raise an exception        
+        #Must raise an exception
         try:
             f = interp.fit(z)
@@ -828,5 +828,5 @@
         #Now try with smoothing parameter
         interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
-        
+
         f = interp.fit(z)
         #f will be different from answerr due to smoothing
@@ -855,5 +855,5 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
@@ -869,5 +869,5 @@
                         [ 15, -17],   #Outside mesh
                         [ 1.0, 2.0],
-                        [ 1.0, 3.0],                                          
+                        [ 1.0, 3.0],
                         [ 2.0, 1.0],
                         [ 3.0, 0.0],
@@ -889,6 +889,6 @@
                         [ 1.0, 0.5],
                         [ 2.0, 0.4],
-                        [ 2.8, 1.2]]                                          
-        
+                        [ 2.8, 1.2]]
+
 
 
@@ -907,5 +907,5 @@
         import Numeric
         z1 = Numeric.take(z1, [0,1,2,4,5,6])
-        answer = Numeric.take(answer, [0,1,2,4,5,6])                
+        answer = Numeric.take(answer, [0,1,2,4,5,6])
         assert allclose(z1, answer)
 
@@ -942,5 +942,5 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
@@ -955,5 +955,5 @@
                         [ 1.0, 1.0],
                         [ 1.0, 2.0],
-                        [ 1.0, 3.0],                                          
+                        [ 1.0, 3.0],
                         [ 2.0, 1.0],
                         [ 3.0, 0.0],
@@ -963,5 +963,5 @@
         #First check that things are OK when using same origin
         mesh_origin = (56, 290000, 618000) #zone, easting, northing
-        data_origin = (56, 290000, 618000) #zone, easting, northing 
+        data_origin = (56, 290000, 618000) #zone, easting, northing
 
 
@@ -971,5 +971,5 @@
                                data_origin = data_origin,
                                mesh_origin = mesh_origin)
-        
+
         z = linear_function(data_points1) #Example z-values
         f = interp.fit(z)                 #Fitted values at vertices
@@ -982,6 +982,6 @@
                         [ 1.0, 0.5],
                         [ 2.0, 0.4],
-                        [ 2.8, 1.2]]                                          
-        
+                        [ 2.8, 1.2]]
+
 
         #Build new A matrix based on new points
@@ -1000,7 +1000,7 @@
         #Then check situtaion where points are relative to a different
         #origin (same zone, though, until we figure that out (FIXME))
-        
+
         mesh_origin = (56, 290000, 618000) #zone, easting, northing
-        data_origin = (56, 10000, 10000) #zone, easting, northing 
+        data_origin = (56, 10000, 10000) #zone, easting, northing
 
         #Shift datapoints according to new origin
@@ -1012,8 +1012,8 @@
         for k in range(len(data_points2)):
             data_points2[k][0] += mesh_origin[1] - data_origin[1]
-            data_points2[k][1] += mesh_origin[2] - data_origin[2]            
-
-
-        
+            data_points2[k][1] += mesh_origin[2] - data_origin[2]
+
+
+
         #Fit surface to mesh
         interp = Interpolation(points, triangles, data_points1,
@@ -1021,5 +1021,5 @@
                                data_origin = data_origin,
                                mesh_origin = mesh_origin)
-        
+
         f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
 
@@ -1055,8 +1055,8 @@
     def test_fit_to_mesh_file(self):
         from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
-             export_mesh_file    
+             export_mesh_file
         import tempfile
         import os
-        
+
         # create a .tsh file, no user outline
         mesh_dic = {}
@@ -1072,8 +1072,8 @@
         mesh_dic['segment_tags'] = ['external',
                                                   'external',
-                                                  'external']       
+                                                  'external']
         mesh_file = tempfile.mktemp(".tsh")
         export_mesh_file(mesh_file,mesh_dic)
-        
+
         # create an .xya file
         point_file = tempfile.mktemp(".xya")
@@ -1081,33 +1081,33 @@
         fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
         fd.close()
-        
+
         mesh_output_file = "new_trianlge.tsh"
         fit_to_mesh_file(mesh_file,
                          point_file,
                          mesh_output_file,
-                         alpha = 0.0)  
+                         alpha = 0.0)
         # load in the .tsh file we just wrote
         mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
-        #print "mesh_dic",mesh_dic 
+        #print "mesh_dic",mesh_dic
         ans =[[0.0, 0.0],
               [5.0, 10.0],
               [5.0,10.0]]
         assert allclose(mesh_dic['vertex_attributes'],ans)
-        
+
         self.failUnless(mesh_dic['vertex_attribute_titles']  ==
                         ['elevation','stage'],
                         'test_fit_to_mesh_file failed')
-        
+
         #clean up
         os.remove(mesh_file)
         os.remove(point_file)
-        os.remove(mesh_output_file)        
-        
+        os.remove(mesh_output_file)
+
     def test_fit_to_mesh_fileII(self):
         from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
-             export_mesh_file    
+             export_mesh_file
         import tempfile
         import os
-        
+
         # create a .tsh file, no user outline
         mesh_dic = {}
@@ -1123,8 +1123,8 @@
         mesh_dic['segment_tags'] = ['external',
                                                   'external',
-                                                  'external']       
+                                                  'external']
         mesh_file = tempfile.mktemp(".tsh")
         export_mesh_file(mesh_file,mesh_dic)
-        
+
         # create an .xya file
         point_file = tempfile.mktemp(".xya")
@@ -1132,10 +1132,10 @@
         fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
         fd.close()
-        
+
         mesh_output_file = "new_triangle.tsh"
         fit_to_mesh_file(mesh_file,
                          point_file,
                          mesh_output_file,
-                         alpha = 0.0)  
+                         alpha = 0.0)
         # load in the .tsh file we just wrote
         mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
@@ -1145,19 +1145,19 @@
                          [1.0, 2.0,5.0, 10.0],
                          [1.0, 2.0,5.0,10.0]])
-        
+
         self.failUnless(mesh_dic['vertex_attribute_titles']  ==
                         ['density', 'temp','elevation','stage'],
                         'test_fit_to_mesh_file failed')
-        
+
         #clean up
         os.remove(mesh_file)
-        os.remove(mesh_output_file)        
+        os.remove(mesh_output_file)
         os.remove(point_file)
-        
+
     def test_fit_to_msh_netcdf_fileII(self):
-        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary,export_mesh_file    
+        from load_mesh.loadASCII import mesh_file_to_mesh_dictionary,export_mesh_file
         import tempfile
         import os
-        
+
         # create a .tsh file, no user outline
         mesh_dic = {}
@@ -1173,8 +1173,8 @@
         mesh_dic['segment_tags'] = ['external',
                                                   'external',
-                                                  'external']       
+                                                  'external']
         mesh_file = tempfile.mktemp(".msh")
         export_mesh_file(mesh_file,mesh_dic)
-        
+
         # create an .xya file
         point_file = tempfile.mktemp(".xya")
@@ -1182,10 +1182,10 @@
         fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
         fd.close()
-        
+
         mesh_output_file = "new_triangle.msh"
         fit_to_mesh_file(mesh_file,
                          point_file,
                          mesh_output_file,
-                         alpha = 0.0)  
+                         alpha = 0.0)
         # load in the .tsh file we just wrote
         mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
@@ -1195,26 +1195,25 @@
                          [1.0, 2.0,5.0, 10.0],
                          [1.0, 2.0,5.0,10.0]])
-        
+
         self.failUnless(mesh_dic['vertex_attribute_titles']  ==
                         ['density', 'temp','elevation','stage'],
                         'test_fit_to_mesh_file failed')
-        
+
         #clean up
         os.remove(mesh_file)
-        os.remove(mesh_output_file)        
+        os.remove(mesh_output_file)
         os.remove(point_file)
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
-
-    #suite = unittest.makeSuite(TestCase,'test_fit_to_msh_netcdf_fileII')
-    #suite = unittest.makeSuite(TestCase,'test_fit_to_mesh_fileII')
+    suite = unittest.makeSuite(Test_Least_Squares,'test')
+
+    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
+    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)
 
-    
-    
-
-
-
+
+
+
+
Index: /inundation/ga/storm_surge/pyvolution/test_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 1018)
@@ -12,10 +12,10 @@
 
 def distance(x, y):
-    return sqrt( sum( (array(x)-array(y))**2 ))		
-        
-class TestCase(unittest.TestCase):
+    return sqrt( sum( (array(x)-array(y))**2 ))
+
+class Test_Mesh(unittest.TestCase):
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
@@ -32,6 +32,6 @@
             msg = 'Should have raised exception'
             raise msg
-        
-        
+
+
     def test_basic_triangle(self):
 
@@ -39,8 +39,8 @@
         b = [4.0, 0.0]
         c = [0.0, 3.0]
-        
+
         points = [a, b, c]
         vertices = [[0,1,2]]
-        mesh = Mesh(points, vertices) 
+        mesh = Mesh(points, vertices)
 
         #Centroid
@@ -57,10 +57,10 @@
         assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
         assert allclose(normals[0, 2:4], [-1.0, 0.0])
-        assert allclose(normals[0, 4:6], [0.0, -1.0])        
+        assert allclose(normals[0, 4:6], [0.0, -1.0])
 
         assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
         assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
-        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])                
-        
+        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])
+
         #Edge lengths
         assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
@@ -78,5 +78,5 @@
 
         V2 = mesh.get_vertex_coordinate(0, 2)
-        assert allclose(V2, [0.0, 3.0])                
+        assert allclose(V2, [0.0, 3.0])
 
 
@@ -86,10 +86,10 @@
         mesh.check_integrity()
 
-        
+
         #Test that the centroid is located 2/3 of the way
         #from each vertex to the midpoint of the opposite side
 
         V = mesh.get_vertex_coordinates()
-        
+
         x0 = V[0,0]
         y0 = V[0,1]
@@ -98,5 +98,5 @@
         x2 = V[0,4]
         y2 = V[0,5]
-        
+
         m0 = [(x1 + x2)/2, (y1 + y2)/2]
         m1 = [(x0 + x2)/2, (y0 + y2)/2]
@@ -108,21 +108,21 @@
         #
         d0 = distance(centroid, [x1, y1])
-        d1 = distance(m1, [x1, y1])        
+        d1 = distance(m1, [x1, y1])
         assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
 
         d0 = distance(centroid, [x2, y2])
         d1 = distance(m2, [x2, y2])
-        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)        
-        
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
+
         #Radius
         d0 = distance(centroid, m0)
         assert d0 == 5.0/6
-        
+
         d1 = distance(centroid, m1)
-        assert d1 == sqrt(73.0/36)        
-
-        d2 = distance(centroid, m2)        
-        assert d2 == sqrt(13.0/9)        
-        
+        assert d1 == sqrt(73.0/36)
+
+        d2 = distance(centroid, m2)
+        assert d2 == sqrt(13.0/9)
+
         assert mesh.radii[0] == min(d0, d1, d2)
         assert mesh.radii[0] == 5.0/6
@@ -134,10 +134,10 @@
         vertices = [[0,3,2], [2,3,1], [1,3,0]]
         new_mesh = Mesh(points, vertices)
-        
+
         assert new_mesh.areas[0] == new_mesh.areas[1]
-        assert new_mesh.areas[1] == new_mesh.areas[2]        
-        assert new_mesh.areas[1] == new_mesh.areas[2]        
-
-        assert new_mesh.areas[1] == mesh.areas[0]/3        
+        assert new_mesh.areas[1] == new_mesh.areas[2]
+        assert new_mesh.areas[1] == new_mesh.areas[2]
+
+        assert new_mesh.areas[1] == mesh.areas[0]/3
 
 
@@ -151,13 +151,13 @@
         vertices = [[0,1,2]]
 
-        mesh = Mesh(points, vertices) 
+        mesh = Mesh(points, vertices)
         centroid = mesh.centroid_coordinates[0]
 
-        
+
         #Test that the centroid is located 2/3 of the way
         #from each vertex to the midpoint of the opposite side
 
         V = mesh.get_vertex_coordinates()
-        
+
         x0 = V[0,0]
         y0 = V[0,1]
@@ -166,5 +166,5 @@
         x2 = V[0,4]
         y2 = V[0,5]
-        
+
         m0 = [(x1 + x2)/2, (y1 + y2)/2]
         m1 = [(x0 + x2)/2, (y0 + y2)/2]
@@ -176,15 +176,15 @@
         #
         d0 = distance(centroid, [x1, y1])
-        d1 = distance(m1, [x1, y1])        
+        d1 = distance(m1, [x1, y1])
         assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
 
         d0 = distance(centroid, [x2, y2])
         d1 = distance(m2, [x2, y2])
-        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)        
-        
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
+
         #Radius
         d0 = distance(centroid, m0)
         d1 = distance(centroid, m1)
-        d2 = distance(centroid, m2)        
+        d2 = distance(centroid, m2)
         assert mesh.radii[0] == min(d0, d1, d2)
 
@@ -197,11 +197,11 @@
         vertices = [[0,3,2], [2,3,1], [1,3,0]]
         new_mesh = Mesh(points, vertices)
-        
+
         assert new_mesh.areas[0] == new_mesh.areas[1]
-        assert new_mesh.areas[1] == new_mesh.areas[2]        
-        assert new_mesh.areas[1] == new_mesh.areas[2]        
-
-        assert new_mesh.areas[1] == mesh.areas[0]/3        
-        
+        assert new_mesh.areas[1] == new_mesh.areas[2]
+        assert new_mesh.areas[1] == new_mesh.areas[2]
+
+        assert new_mesh.areas[1] == mesh.areas[0]/3
+
 
         #Test that points are arranged in a counter clock wise order
@@ -217,19 +217,19 @@
         e = [2.0, 2.0]
         points = [a, b, c, e]
-        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
         mesh = Mesh(points, vertices)
-        
+
         assert mesh.areas[0] == 2.0
 
-        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])        
+        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
 
 
         #Test that points are arranged in a counter clock wise order
-        mesh.check_integrity()        
-
-
-            
+        mesh.check_integrity()
+
+
+
     def test_more_triangles(self):
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -240,10 +240,10 @@
 
         points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
-        mesh = Mesh(points, vertices)        
+        mesh = Mesh(points, vertices)
 
         #Test that points are arranged in a counter clock wise order
-        mesh.check_integrity()        
+        mesh.check_integrity()
 
         assert mesh.areas[0] == 2.0
@@ -252,6 +252,6 @@
         assert mesh.areas[3] == 2.0
         assert mesh.areas[4] == 8.0
-        assert mesh.areas[5] == 4.0                                
-        
+        assert mesh.areas[5] == 4.0
+
         assert mesh.edgelengths[1,0] == 2.0
         assert mesh.edgelengths[1,1] == 2.0
@@ -259,7 +259,7 @@
 
         assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
-        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 
+        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
         assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
-        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])         
+        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
 
 
@@ -278,5 +278,5 @@
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
         mesh = Mesh(points, vertices)
-        
+
         mesh.check_integrity()
 
@@ -284,11 +284,11 @@
         T = mesh
         tid = 0
-        assert T.number_of_boundaries[tid] == 2    
+        assert T.number_of_boundaries[tid] == 2
         assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
         assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
-        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)        
+        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)
 
         tid = 1
-        assert T.number_of_boundaries[tid] == 0            
+        assert T.number_of_boundaries[tid] == 0
         assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
         assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
@@ -296,11 +296,11 @@
 
         tid = 2
-        assert T.number_of_boundaries[tid] == 2            
+        assert T.number_of_boundaries[tid] == 2
         assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
         assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
-        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)        
+        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)
 
         tid = 3
-        assert T.number_of_boundaries[tid] == 2            
+        assert T.number_of_boundaries[tid] == 2
         assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
         assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
@@ -311,22 +311,22 @@
         assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
         assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
-        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)        
-
-        tid = 1        
+        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)
+
+        tid = 1
         assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
         assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
         assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
 
-        tid = 2        
+        tid = 2
         assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
         assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
-        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)        
-
-        tid = 3        
+        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
+
+        tid = 3
         assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
         assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
         assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
 
-        
+
     def test_rectangular_mesh_basic(self):
         M=1
@@ -336,5 +336,5 @@
 	mesh = Mesh(points, vertices, boundary)
 
-        #Test that points are arranged in a counter clock wise order        
+        #Test that points are arranged in a counter clock wise order
         mesh.check_integrity()
 
@@ -343,21 +343,21 @@
         points, vertices, boundary = rectangular(M, N)
 	mesh = Mesh(points, vertices, boundary)
-	
-        #Test that points are arranged in a counter clock wise order        
+
+        #Test that points are arranged in a counter clock wise order
         mesh.check_integrity()
 
         #assert mesh.boundary[(7,1)] == 2 # top
         assert mesh.boundary[(7,1)] == 'top' # top
-        assert mesh.boundary[(3,1)] == 'top' # top                
-
-        
+        assert mesh.boundary[(3,1)] == 'top' # top
+
+
     def test_boundary_tags(self):
-        
+
 
         points, vertices, boundary = rectangular(4, 4)
 	mesh = Mesh(points, vertices, boundary)
-	
-
-        #Test that points are arranged in a counter clock wise order        
+
+
+        #Test that points are arranged in a counter clock wise order
         mesh.check_integrity()
 
@@ -369,6 +369,6 @@
 
         for k in [24, 26, 28, 30]:
-            assert mesh.boundary[(k,2)] == 'right'            
-            
+            assert mesh.boundary[(k,2)] == 'right'
+
         for k in [7, 15, 23, 31]:
             assert mesh.boundary[(k,1)] == 'top'
@@ -377,5 +377,5 @@
 
 
-        
+
     def test_rectangular_mesh(self):
         M=4
@@ -383,9 +383,9 @@
         len1 = 100.0
         len2 = 17.0
-	
+
         points, vertices, boundary = rectangular(M, N, len1, len2)
 	mesh = Mesh(points, vertices, boundary)
 
-        assert len(mesh) == 2*M*N        
+        assert len(mesh) == 2*M*N
 
         for i in range(len(mesh)):
@@ -396,7 +396,7 @@
             assert mesh.edgelengths[i, 1] == len1/M #x direction
             assert mesh.edgelengths[i, 2] == len2/N #y direction
-            
-        #Test that points are arranged in a counter clock wise order      
-        mesh.check_integrity()    
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()
 
 
@@ -404,5 +404,5 @@
         #Check that integers don't cause trouble
         N = 16
-	
+
         points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
 	mesh = Mesh(points, vertices, boundary)
@@ -422,5 +422,5 @@
         #bac, bce, ecf, dbe
         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
-        mesh = Mesh(points, vertices)        
+        mesh = Mesh(points, vertices)
         mesh.check_integrity()
 
@@ -428,11 +428,11 @@
         T = mesh
         tid = 0
-        assert T.number_of_boundaries[tid] == 2    
-        assert T.surrogate_neighbours[tid, 0] == tid 
-        assert T.surrogate_neighbours[tid, 1] == 1   
-        assert T.surrogate_neighbours[tid, 2] == tid 
+        assert T.number_of_boundaries[tid] == 2
+        assert T.surrogate_neighbours[tid, 0] == tid
+        assert T.surrogate_neighbours[tid, 1] == 1
+        assert T.surrogate_neighbours[tid, 2] == tid
 
         tid = 1
-        assert T.number_of_boundaries[tid] == 0            
+        assert T.number_of_boundaries[tid] == 0
         assert T.surrogate_neighbours[tid, 0] == 2
         assert T.surrogate_neighbours[tid, 1] == 3
@@ -440,12 +440,12 @@
 
         tid = 2
-        assert T.number_of_boundaries[tid] == 2            
-        assert T.surrogate_neighbours[tid, 0] == tid 
-        assert T.surrogate_neighbours[tid, 1] == tid 
+        assert T.number_of_boundaries[tid] == 2
+        assert T.surrogate_neighbours[tid, 0] == tid
+        assert T.surrogate_neighbours[tid, 1] == tid
         assert T.surrogate_neighbours[tid, 2] == 1
 
         tid = 3
-        assert T.number_of_boundaries[tid] == 2            
-        assert T.surrogate_neighbours[tid, 0] == 1 
+        assert T.number_of_boundaries[tid] == 2
+        assert T.surrogate_neighbours[tid, 0] == 1
         assert T.surrogate_neighbours[tid, 1] == tid
         assert T.surrogate_neighbours[tid, 2] == tid
@@ -470,13 +470,13 @@
                      (2, 1): 'Fourth',
                      (3, 1): 'Fifth',
-                     (3, 2): 'Sixth'}                                         
-                     
-        
-        mesh = Mesh(points, vertices, boundary)        
+                     (3, 2): 'Sixth'}
+
+
+        mesh = Mesh(points, vertices, boundary)
         mesh.check_integrity()
 
 
         #Check enumeration
-        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):        
+        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
         #    b = -k-1
         #    assert mesh.neighbours[vol_id, edge_id] == b
@@ -502,8 +502,8 @@
                      (2, 1): 'Fourth',
                      #(3, 1): 'Fifth',  #Skip this
-                     (3, 2): 'Sixth'}                                         
-                     
-        
-        mesh = Mesh(points, vertices, boundary)        
+                     (3, 2): 'Sixth'}
+
+
+        mesh = Mesh(points, vertices, boundary)
         mesh.check_integrity()
 
@@ -513,5 +513,5 @@
 
         #Check enumeration
-        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):        
+        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
         #    b = -k-1
         #    assert mesh.neighbours[vol_id, edge_id] == b
@@ -535,8 +535,8 @@
                      (2, 1): 'Fourth',
                      #(3, 1): 'Fifth',  #Skip this
-                     (3, 2): 'Sixth'}                                         
-                     
-        
-        mesh = Mesh(points, vertices) #, boundary)        
+                     (3, 2): 'Sixth'}
+
+
+        mesh = Mesh(points, vertices) #, boundary)
         mesh.check_integrity()
 
@@ -545,11 +545,11 @@
         assert mesh.boundary[ (0, 2) ] == default_boundary_tag
         assert mesh.boundary[ (2, 0) ] == default_boundary_tag
-        assert mesh.boundary[ (2, 1) ] == default_boundary_tag                
+        assert mesh.boundary[ (2, 1) ] == default_boundary_tag
         assert mesh.boundary[ (3, 1) ] == default_boundary_tag
-        assert mesh.boundary[ (3, 2) ] == default_boundary_tag                
+        assert mesh.boundary[ (3, 2) ] == default_boundary_tag
 
 
         #Check enumeration
-        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):        
+        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
         #    b = -k-1
         #    assert mesh.neighbours[vol_id, edge_id] == b
@@ -575,5 +575,5 @@
         #Too few points
         try:
-            mesh = Mesh([points[0]], vertices)        
+            mesh = Mesh([points[0]], vertices)
         except AssertionError:
             pass
@@ -583,13 +583,5 @@
         #Too few points - 1 element
         try:
-            mesh = Mesh([points[0]], [vertices[0]])        
-        except AssertionError:
-            pass
-        else:
-            raise 'Should have raised an exception'        
-        
-        #Wrong dimension of vertices
-        try:
-            mesh = Mesh(points, vertices[0])        
+            mesh = Mesh([points[0]], [vertices[0]])
         except AssertionError:
             pass
@@ -597,7 +589,15 @@
             raise 'Should have raised an exception'
 
+        #Wrong dimension of vertices
+        try:
+            mesh = Mesh(points, vertices[0])
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
         #Unsubscriptable coordinates object raises exception
         try:
-            mesh = Mesh(points[0], [vertices[0]])        
+            mesh = Mesh(points[0], [vertices[0]])
         except AssertionError:
             pass
@@ -610,5 +610,5 @@
         #Not specifying all boundary tags
         #try:
-        #    mesh = Mesh(points, vertices, {(3,0): 'x'})        
+        #    mesh = Mesh(points, vertices, {(3,0): 'x'})
         #except AssertionError:
         #    pass
@@ -618,5 +618,5 @@
         #Specifying wrong non existing segment
         try:
-            mesh = Mesh(points, vertices, {(5,0): 'x'})        
+            mesh = Mesh(points, vertices, {(5,0): 'x'})
         except AssertionError:
             pass
@@ -634,5 +634,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -648,5 +648,5 @@
                                                  'all':[0,1,2,3,4,5]})
 
-        
+
     def test_boundary_polygon(self):
         from mesh_factory import rectangular
@@ -654,9 +654,9 @@
         from Numeric import zeros, Float
         from util import inside_polygon
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(2, 2)
         mesh = Mesh(points, vertices, boundary)
-        
+
 
         P = mesh.get_boundary_polygon()
@@ -668,5 +668,5 @@
         for p in points:
             #print p, P
-            assert inside_polygon(p, P) 
+            assert inside_polygon(p, P)
 
 
@@ -675,5 +675,5 @@
         from Numeric import zeros, Float
         from util import inside_polygon
-        
+
         #Points
         a = [0.0, 0.0] #0
@@ -695,7 +695,7 @@
 
         mesh = Mesh(points, vertices)
-        
-        mesh.check_integrity()
-        
+
+        mesh.check_integrity()
+
         P = mesh.get_boundary_polygon()
 
@@ -705,5 +705,5 @@
         for p in points:
             #print p, P
-            assert inside_polygon(p, P) 
+            assert inside_polygon(p, P)
 
 
@@ -711,9 +711,9 @@
         """Same as II but vertices ordered differently
         """
-        
+
         from mesh import Mesh
         from Numeric import zeros, Float
         from util import inside_polygon
-        
+
         #Points
         a = [0.0, 0.0] #0
@@ -738,5 +738,5 @@
         mesh.check_integrity()
 
-        
+
         P = mesh.get_boundary_polygon()
 
@@ -745,17 +745,16 @@
 
         for p in points:
-            assert inside_polygon(p, P) 
-
-
-        
+            assert inside_polygon(p, P)
+
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
 
-    
-    
-
-
-
+
+
+
+
Index: /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1018)
@@ -23,17 +23,17 @@
      Transmissive_boundary
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_pmesh2domain(unittest.TestCase):
 
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
-               
+
     def test_pmesh2Domain(self):
          import os
          import tempfile
-       
+
          fileName = tempfile.mktemp(".tsh")
          file = open(fileName,"w")
@@ -66,5 +66,5 @@
 0 # <x> <y> <attribute>...Mesh Regions... \n")
          file.close()
-       
+
          tags = {}
          b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
@@ -74,15 +74,15 @@
          tags["2"] = b2
          tags["3"] = b3
-         
+
          domain = pmesh_to_domain_instance(fileName, Domain)
          os.remove(fileName)
 
-         ## check the quantities       
+         ## check the quantities
          #print domain.quantities['elevation'].vertex_values
          answer = [[0., 8., 0.],
                    [0., 10., 8.]]
          assert allclose(domain.quantities['elevation'].vertex_values,
-                        answer) 
-         
+                        answer)
+
          #print domain.quantities['stage'].vertex_values
          answer = [[0., 12., 10.],
@@ -90,5 +90,5 @@
          assert allclose(domain.quantities['stage'].vertex_values,
                         answer)
-         
+
          #print domain.quantities['friction'].vertex_values
          answer = [[0.01, 0.04, 0.03],
@@ -100,5 +100,5 @@
          assert allclose(domain.tagged_elements['dsg'][0],0)
          assert allclose(domain.tagged_elements['ole nielsen'][0],1)
-         
+
          self.failUnless( domain.boundary[(1, 0)]  == '1',
                           "test_tags_to_boundaries  failed. Single boundary wasn't added.")
@@ -123,7 +123,7 @@
          meshDict['triangle_tags'] = [6.6,6.6]
          meshDict['triangle_neighbors'] = [[-1,-1,1],[-1,0,-1]]
-         meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]] 
+         meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]]
          meshDict['segment_tags'] = [2,3,1,1]
-       
+
          domain = Domain.pmesh_dictionary_to_domain(meshDict)
 
@@ -143,12 +143,12 @@
          boundary_value = Boundary_value.instances[id]
          boundary_obj = boundary_value.boundary_object
-         
+
          self.failUnless( boundary_obj  == b1,
                           "test_tags_to_boundaries  failed. Single boundary wasn't added.")
-         
+
          inverted_id = Volume.instances[0].neighbours[1]
          id = -(inverted_id+1)
          boundary_value = Boundary_value.instances[id]
-         boundary_obj = boundary_value.boundary_object   
+         boundary_obj = boundary_value.boundary_object
          self.failUnless( boundary_obj  == b3,
                           "test_tags_to_boundaries  failed. Single boundary wasn't added.")
@@ -156,11 +156,10 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase, 'test')
+    suite = unittest.makeSuite(Test_pmesh2domain, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
-    
-
-    
 
 
 
+
+
Index: /inundation/ga/storm_surge/pyvolution/test_quad.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_quad.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_quad.py	(revision 1018)
@@ -9,9 +9,9 @@
 #-------------------------------------------------------------
 
-class TestCase(unittest.TestCase):
+class Test_Quad(unittest.TestCase):
 
     def setUp(self):
         self.cell = Cell(100, 140, 0, 40, 'cell')
-        
+
         a = [3, 107]
         b = [5, 107]
@@ -24,11 +24,11 @@
 
         points = [a, b, c, d, e, f, g, h]
-        #bac, bce, ecf, dbe, daf, dae 
-        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]] 
+        #bac, bce, ecf, dbe, daf, dae
+        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
 
 	mesh = Mesh(points, vertices)
         self.mesh = mesh
         Cell.initialise(mesh)
-        
+
     def tearDown(self):
         pass
@@ -37,16 +37,16 @@
         self.cell.insert(0)
         self.cell.insert(1)
-        
+
         result = self.cell.retrieve()
         assert type(result) in [types.ListType,types.TupleType],\
-                                'should be a list'        
+                                'should be a list'
         self.assertEqual(len(result),2)
-        
+
     def test_add_points_2_cellII(self):
         self.cell.insert([0,1,2,3,4,5,6,7])
-        
+
         result = self.cell.retrieve()
         assert type(result) in [types.ListType,types.TupleType],\
-	                        'should be a list'        
+	                        'should be a list'
         self.assertEqual(len(result),8)
 
@@ -55,35 +55,35 @@
         self.cell.insert([0,1,2,3,4,5,6,7])
 	self.cell.split(4)
-        
+
         result =  self.cell.search(x = 1, y = 101)
         assert type(result) in [types.ListType,types.TupleType],\
-	                        'should be a list'        
+	                        'should be a list'
         self.assertEqual(result, [0,1,2,3])
-	
-	
+
+
     def test_clear_1(self):
         self.cell.insert([0,1,2,3,4,5,6,7])
-	assert self.cell.count() == 8	
+	assert self.cell.count() == 8
 	self.cell.clear()
-	
+
         #This one actually revealed a bug :-)
-	assert self.cell.count() == 0	
-		
+	assert self.cell.count() == 0
+
     def test_clear_2(self):
         self.cell.insert([0,1,2,3,4,5,6,7])
-	assert self.cell.count() == 8	
+	assert self.cell.count() == 8
 	self.cell.split(2)
-	assert self.cell.count() == 8		
-		
+	assert self.cell.count() == 8
+
 	self.cell.clear()
-	assert self.cell.count() == 0	
-	
-	
-	
+	assert self.cell.count() == 0
+
+
+
     def test_split(self):
-        self.cell.insert([0,1,2,3,4,5,6,7], split = False)    
-        
+        self.cell.insert([0,1,2,3,4,5,6,7], split = False)
+
 	#No children yet
-	assert self.cell.children is None 
+	assert self.cell.children is None
 	assert self.cell.count() == 8
 
@@ -92,28 +92,28 @@
 	#self.cell.show()
 	#self.cell.show_all()
-	
-	
+
+
 	#Now there are children
-	assert self.cell.children is not None 
-	assert self.cell.count() == 8	
-	
-	
-		
+	assert self.cell.children is not None
+	assert self.cell.count() == 8
+
+
+
     def test_collapse(self):
-        self.cell.insert([0,1,2,3,4,5,6,7], split = False)    
-        
+        self.cell.insert([0,1,2,3,4,5,6,7], split = False)
+
         #Split maximally
 	self.cell.split(1)
-	
+
 	#Now there are children
-	assert self.cell.children is not None 
-	assert self.cell.count() == 8	
-	
+	assert self.cell.children is not None
+	assert self.cell.count() == 8
+
 	#Collapse
 	self.cell.collapse(8)
-	
+
 	#No children
-	assert self.cell.children is None 
-	assert self.cell.count() == 8	
+	assert self.cell.children is None
+	assert self.cell.count() == 8
 
     def test_build_quadtree(self):
@@ -123,5 +123,5 @@
 
         #Q.show()
-        
+
         result = Q.search(3, 105)
         assert type(result) in [types.ListType,types.TupleType],\
@@ -130,9 +130,9 @@
         self.assertEqual(result, [0,1,2,3])
 
-     
+
     def test_build_quadtreeII(self):
 
         self.cell = Cell(100, 140, 0, 40, 'cell')
-        
+
         p0 = [34.6292076111,-7999.92529297]
         p1 = [8000.0, 7999.0]
@@ -141,5 +141,5 @@
 
         points = [p0,p1,p2, p3]
-        #bac, bce, ecf, dbe, daf, dae 
+        #bac, bce, ecf, dbe, daf, dae
         vertices = [[0,1,2],[0,2,3]]
 
@@ -148,10 +148,9 @@
         #This was causing round off error
         Q = build_quadtree(mesh)
-        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
 
-    mysuite = unittest.makeSuite(TestCase,'test')
+    mysuite = unittest.makeSuite(Test_Sparse,'test')
     runner = unittest.TextTestRunner()
     runner.run(mysuite)
-
Index: /inundation/ga/storm_surge/pyvolution/test_quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 1018)
@@ -9,9 +9,9 @@
 from Numeric import allclose, array, ones, Float
 
-        
-class TestCase(unittest.TestCase):
+
+class Test_Quantity(unittest.TestCase):
     def setUp(self):
         from domain import Domain
-        
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -26,10 +26,10 @@
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
 
-        self.mesh1 = Domain(points[:3], [elements[0]])        
+        self.mesh1 = Domain(points[:3], [elements[0]])
         self.mesh1.check_integrity()
-        
-        self.mesh4 = Domain(points, elements)        
+
+        self.mesh4 = Domain(points, elements)
         self.mesh4.check_integrity()
-        
+
     def tearDown(self):
         pass
@@ -38,5 +38,5 @@
 
     def test_creation(self):
-        
+
         quantity = Quantity(self.mesh1, [[1,2,3]])
         assert allclose(quantity.vertex_values, [[1.,2.,3.]])
@@ -55,9 +55,9 @@
             pass
         except:
-            raise 'Should have raised "mising mesh object" error'        
-        
+            raise 'Should have raised "mising mesh object" error'
+
 
     def test_creation_zeros(self):
-       
+
         quantity = Quantity(self.mesh1)
         assert allclose(quantity.vertex_values, [[0.,0.,0.]])
@@ -66,12 +66,12 @@
         quantity = Quantity(self.mesh4)
         assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
-                                                 [0.,0.,0.], [0.,0.,0.]])        
-
-        
+                                                 [0.,0.,0.], [0.,0.,0.]])
+
+
     def test_interpolation(self):
         quantity = Quantity(self.mesh1, [[1,2,3]])
         assert allclose(quantity.centroid_values, [2.0]) #Centroid
 
-        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 
+        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
 
 
@@ -88,6 +88,6 @@
                                                [7.5, 0.5, 1.],
                                                [-5, -2.5, 7.5]])
-        
-        #print quantity.edge_values        
+
+        #print quantity.edge_values
         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
                                                [5., 5., 5.],
@@ -120,5 +120,5 @@
         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
         assert allclose(quantity.vertex_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])        
+                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
@@ -135,5 +135,5 @@
                             location = 'edges')
         assert allclose(quantity.edge_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])        
+                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
 
         #Test exceptions
@@ -174,5 +174,5 @@
                                                [3, 3, 3],
                                                [3, 3, 3],
-                                               [3, 3, 3]])        
+                                               [3, 3, 3]])
 
 
@@ -186,11 +186,11 @@
         #print "quantity.vertex_values",quantity.vertex_values
         assert allclose(quantity.vertex_values,
-                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])        
+                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
         assert allclose(quantity.centroid_values,
                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
         assert allclose(quantity.edge_values,
-                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])                
-
-        
+                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
+
+
         quantity.set_values(f, location = 'centroids')
         assert allclose(quantity.centroid_values,
@@ -201,12 +201,12 @@
         quantity = Quantity(self.mesh4)
 
-	#Try constants first	
+	#Try constants first
 	const = 5
-        quantity.set_values(const, location = 'vertices')	
+        quantity.set_values(const, location = 'vertices')
 	#print 'Q', quantity.get_integral()
-		
-	assert allclose(quantity.get_integral(), 
+
+	assert allclose(quantity.get_integral(),
 			self.mesh4.get_area() * const)
-	
+
 	#Try with a linear function
         def f(x, y):
@@ -214,9 +214,9 @@
 
         quantity.set_values(f, location = 'vertices')
-	
-	
+
+
 	ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
-        
-	assert allclose (quantity.get_integral(), ref_integral)	
+
+	assert allclose (quantity.get_integral(), ref_integral)
 
 
@@ -256,5 +256,5 @@
         #Centroid
         assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
-        
+
         assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
                                                [3., 2.5, 1.5],
@@ -268,5 +268,5 @@
         quantity = Conserved_quantity(self.mesh4)
 
-        #Set up for a gradient of (3,0) at mid triangle 
+        #Set up for a gradient of (3,0) at mid triangle
         quantity.set_values([2.0, 4.0, 8.0, 2.0],
                             location = 'centroids')
@@ -294,12 +294,12 @@
                                                  [0., 0.,  6.]])
 
-		
+
     def test_second_order_extrapolation2(self):
-        quantity = Conserved_quantity(self.mesh4)        
+        quantity = Conserved_quantity(self.mesh4)
 
         #Set up for a gradient of (3,1), f(x) = 3x+y
         quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
                             location = 'centroids')
-        
+
         #Gradients
         a, b = quantity.compute_gradients()
@@ -314,8 +314,8 @@
 
         quantity.extrapolate_second_order()
-        
+
         #print quantity.vertex_values
         assert allclose(quantity.vertex_values[1,0], 2.0)
-        assert allclose(quantity.vertex_values[1,1], 6.0)        
+        assert allclose(quantity.vertex_values[1,1], 6.0)
         assert allclose(quantity.vertex_values[1,2], 8.0)
 
@@ -323,5 +323,5 @@
 
     def test_first_order_extrapolator(self):
-        quantity = Conserved_quantity(self.mesh4)                
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids
@@ -338,10 +338,10 @@
 
     def test_second_order_extrapolator(self):
-        quantity = Conserved_quantity(self.mesh4)                        
-
-        #Set up for a gradient of (3,0) at mid triangle 
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Set up for a gradient of (3,0) at mid triangle
         quantity.set_values([2.0, 4.0, 8.0, 2.0],
                             location = 'centroids')
-        
+
 
 
@@ -353,12 +353,12 @@
         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
         assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
-        
+
         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
-        
+
         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
         assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
 
-        
+
         #Assert that quantities are conserved
         from Numeric import sum
@@ -366,8 +366,8 @@
             assert allclose (quantity.centroid_values[k],
                              sum(quantity.vertex_values[k,:])/3)
-        
-        
-
-        
+
+
+
+
 
     def test_limiter(self):
@@ -384,13 +384,13 @@
         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
-        
+
         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
-        
+
         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
         assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
 
 
-        
+
         #Assert that quantities are conserved
         from Numeric import sum
@@ -398,5 +398,5 @@
             assert allclose (quantity.centroid_values[k],
                              sum(quantity.vertex_values[k,:])/3)
-        
+
 
     def test_limiter2(self):
@@ -418,7 +418,7 @@
         quantity.limit()
 
-        
+
         assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
-        
+
 
         #Assert that quantities are conserved
@@ -429,9 +429,9 @@
 
 
-        
+
 
 
     def test_distribute_first_order(self):
-        quantity = Conserved_quantity(self.mesh4)        
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids
@@ -444,5 +444,5 @@
 
         #Interpolate
-        quantity.interpolate_from_vertices_to_edges()        
+        quantity.interpolate_from_vertices_to_edges()
 
         assert allclose(quantity.vertex_values,
@@ -453,5 +453,5 @@
 
     def test_distribute_second_order(self):
-        quantity = Conserved_quantity(self.mesh4)        
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids
@@ -498,5 +498,5 @@
         sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
         denom = ones(4, Float)-timestep*sem
-        
+
         x = array([1, 2, 3, 4])/denom
         assert allclose( quantity.centroid_values, x)
@@ -512,10 +512,10 @@
         #Set explicit_update
         quantity.explicit_update = array( [4.,3.,2.,1.] )
-       
+
         #Set semi implicit update
         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
 
         #Update with given timestep
-        timestep = 0.1        
+        timestep = 0.1
         quantity.update(0.1)
 
@@ -525,10 +525,10 @@
         x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
         x /= denom
-        assert allclose( quantity.centroid_values, x)        
-
-        
-
-
-    #Test smoothing    
+        assert allclose( quantity.centroid_values, x)
+
+
+
+
+    #Test smoothing
     def test_smoothing(self):
 
@@ -546,17 +546,17 @@
         domain.reduction = mean
 
-        
+
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: x)        
+        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
 
@@ -566,9 +566,9 @@
         h = 0.03
         for i in range(stage.shape[0]):
-            if i % 2 == 0:            
+            if i % 2 == 0:
                 stage[i,:] = bed[i,:] + h
             else:
                 stage[i,:] = bed[i,:]
-                
+
         domain.set_quantity('stage', stage)
 
@@ -579,11 +579,11 @@
         Q = stage.vertex_values
 
-        
-	assert A.shape[0] == 9	
-	assert V.shape[0] == 8	
-	assert V.shape[1] == 3			
-	
-        #First four points 
-        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)  
+
+	assert A.shape[0] == 9
+	assert V.shape[0] == 8
+	assert V.shape[1] == 3
+
+        #First four points
+        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
         assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
         assert allclose(A[2], Q[3,0])
@@ -593,15 +593,15 @@
         assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
                                Q[5,0] + Q[6,2] + Q[7,1])/6)
-                                 
+
 
         #Check V
         assert allclose(V[0,:], [3,4,0])
         assert allclose(V[1,:], [1,0,4])
-        assert allclose(V[2,:], [4,5,1])                        
+        assert allclose(V[2,:], [4,5,1])
         assert allclose(V[3,:], [2,1,5])
-        assert allclose(V[4,:], [6,7,3])        
+        assert allclose(V[4,:], [6,7,3])
         assert allclose(V[5,:], [4,3,7])
         assert allclose(V[6,:], [7,8,4])
-        assert allclose(V[7,:], [5,4,8])                        
+        assert allclose(V[7,:], [5,4,8])
 
         #Get smoothed stage with XY
@@ -609,5 +609,5 @@
 
         assert allclose(A, A1)
-        assert allclose(V, V1)        
+        assert allclose(V, V1)
 
         #Check XY
@@ -616,5 +616,5 @@
 
         assert allclose(X[7], 1.0)
-        assert allclose(Y[7], 0.5)                 
+        assert allclose(Y[7], 0.5)
 
 
@@ -628,5 +628,5 @@
         from util import mean
 
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(2, 2)
@@ -637,11 +637,11 @@
         domain.reduction = mean
 
-        
+
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: x)        
+        domain.set_quantity('elevation', lambda x,y: x)
         domain.set_quantity('friction', 0.03)
 
 
-        ###################### 
+        ######################
         #Initial condition - with jumps
 
@@ -651,13 +651,13 @@
         h = 0.03
         for i in range(stage.shape[0]):
-            if i % 2 == 0:            
+            if i % 2 == 0:
                 stage[i,:] = bed[i,:] + h
             else:
                 stage[i,:] = bed[i,:]
-                
+
         domain.set_quantity('stage', stage)
 
         #Get stage
-        stage = domain.quantities['stage']        
+        stage = domain.quantities['stage']
         A, V = stage.get_vertex_values(xy=False, smooth=False)
         Q = stage.vertex_values.flat
@@ -666,17 +666,17 @@
             assert allclose(A[k], Q[k])
 
-            
+
         for k in range(8):
             assert V[k, 0] == 3*k
             assert V[k, 1] == 3*k+1
-            assert V[k, 2] == 3*k+2            
-            
+            assert V[k, 2] == 3*k+2
+
 
 
         X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
 
-        
+
         assert allclose(A, A1)
-        assert allclose(V, V1)        
+        assert allclose(V, V1)
 
         #Check XY
@@ -686,8 +686,8 @@
         assert allclose(Y[4], 0.0)
         assert allclose(X[12], 1.0)
-        assert allclose(Y[12], 0.0)                
-        
-        
- 
+        assert allclose(Y[12], 0.0)
+
+
+
     def set_array_values_by_index(self):
 
@@ -695,5 +695,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 1)
@@ -701,5 +701,5 @@
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements 
+        #print "domain.number_of_elements ",domain.number_of_elements
         quantity = Quantity(domain,[[1,1,1],[2,2,2]])
         value = [7]
@@ -709,5 +709,5 @@
                                            indexes = indexes)
         #print "quantity.centroid_values",quantity.centroid_values
-        
+
         assert allclose(quantity.centroid_values, [1,7])
 
@@ -717,5 +717,5 @@
         quantity.set_array_values([15,20,25], indexes = indexes)
         assert allclose(quantity.centroid_values, [1,20])
-        
+
     def test_setting_some_vertex_values(self):
         """
@@ -725,5 +725,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -731,5 +731,5 @@
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements 
+        #print "domain.number_of_elements ",domain.number_of_elements
         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5],[6,6,6]])
@@ -741,5 +741,5 @@
         #print "quantity.centroid_values",quantity.centroid_values
         assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
-        
+
         value = [[15,20,25]]
         quantity.set_values(value, indexes = indexes)
@@ -758,5 +758,5 @@
         assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
 
-        
+
         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5],[6,6,6]])
@@ -764,5 +764,5 @@
         #this will be per unique vertex, indexing the vertices
         #print "quantity.vertex_values",quantity.vertex_values
-        quantity.set_values(values, indexes = [0,1,5]) 
+        quantity.set_values(values, indexes = [0,1,5])
         #print "quantity.vertex_values",quantity.vertex_values
         assert allclose(quantity.vertex_values[0], [1,50,10])
@@ -780,9 +780,9 @@
                                     [4,4,4],[5,5,5],[6,6,6]]
         quantity.set_values(values)
-        
+
         # testing the standard set values by vertex
         # indexed by vertex_id in general_mesh.coordinates
         values = [0,1,2,3,4,5,6,7]
-        
+
         quantity.set_values(values)
         #print "1 quantity.vertex_values",quantity.vertex_values
@@ -793,5 +793,5 @@
                                                 [ 6.,  7.,  2.],
                                                 [ 3.,  2.,  7.]])
-    
+
     def test_setting_unique_vertex_values(self):
         """
@@ -801,5 +801,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -807,5 +807,5 @@
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements 
+        #print "domain.number_of_elements ",domain.number_of_elements
         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5]])
@@ -819,6 +819,6 @@
         assert allclose(quantity.vertex_values[1], [7,1,7])
         assert allclose(quantity.vertex_values[2], [7,2,7])
-        
-    
+
+
     def test_get_values(self):
         """
@@ -828,5 +828,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -838,8 +838,8 @@
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements 
+        #print "domain.number_of_elements ",domain.number_of_elements
         quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5]])
-        
+
         #print "quantity.get_values(location = 'unique vertices')", \
         #      quantity.get_values(location = 'unique vertices')
@@ -848,9 +848,9 @@
         #      quantity.get_values(indexes=[0,1,2,3,4,5,6,7], \
         #                          location = 'unique vertices')
-        
+
         answer = [0.5,2,4,5,0,1,3,4.5]
         assert allclose(answer,
                         quantity.get_values(location = 'unique vertices'))
-        
+
         indexes = [0,5,3]
         answer = [0.5,1,5]
@@ -861,5 +861,5 @@
         #print "quantity.get_values(location = 'centroids') ",\
         #      quantity.get_values(location = 'centroids')
-        
+
     def test_getting_some_vertex_values(self):
         """
@@ -869,5 +869,5 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
@@ -879,5 +879,5 @@
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
-        #print "domain.number_of_elements ",domain.number_of_elements 
+        #print "domain.number_of_elements ",domain.number_of_elements
         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
                                     [4,4,4],[5,5,5],[6,6,6]])
@@ -893,5 +893,5 @@
                         quantity.get_values(location = 'centroids'))
 
-        
+
         value = [[15,20,25]]
         quantity.set_values(value, indexes = indexes)
@@ -902,5 +902,5 @@
                         quantity.get_values(location = 'edges'))
 
-        # get a subset of elements 
+        # get a subset of elements
         subset = quantity.get_values(location='centroids', indexes=[0,5])
         answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
@@ -920,13 +920,12 @@
         assert allclose(subset, answer)
 
-        
-        
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Quantity,'test')
     #print "restricted test"
-    #suite = unittest.makeSuite(TestCase,'test_set_vertex_values_subset')
+    #suite = unittest.makeSuite(Test_Quantity,'test_set_vertex_values_subset')
     runner = unittest.TextTestRunner()
     runner.run(suite)
 
-
Index: /inundation/ga/storm_surge/pyvolution/test_region.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_region.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_region.py	(revision 1018)
@@ -12,10 +12,10 @@
 def add_x_y(x, y):
     return x+y
-        
-class TestCase(unittest.TestCase):
+
+class Test_Region(unittest.TestCase):
     def setUp(self):
         pass
 
-        
+
     def tearDown(self):
         pass
@@ -29,8 +29,8 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -39,5 +39,5 @@
                                                  'all':[0,1,2,3,4,5]})
 
-        
+
         #Set friction
         manning = 0.07
@@ -66,6 +66,6 @@
                          [ 11.0,  11.0,  11.0],
                          [ 11.0,  11.0,  11.0]])
-                         
-        # trying a function                 
+
+        # trying a function
         domain.set_region(Set_region('top', 'friction', add_x_y))
         #print domain.quantities['friction'].get_values()
@@ -77,5 +77,5 @@
                          [ 5./3,  2.0,  2./3],
                          [ 1.0,  2./3,  2.0]])
-                         
+
         domain.set_quantity('elevation', 10.0)
         domain.set_quantity('stage', 10.0)
@@ -89,5 +89,5 @@
                          [ 11.0,  11.0,  11.0],
                          [ 11.0,  11.0,  11.0]])
-    
+
     def test_unique_vertices(self):
         """
@@ -97,8 +97,8 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -106,5 +106,5 @@
                                                  'top':[4,5],
                                                  'all':[0,1,2,3,4,5]})
-      
+
         #Set friction
         manning = 0.07
@@ -122,5 +122,5 @@
                          [ 0.07,  0.07,  0.07]])
 
-        
+
     def test_unique_verticesII(self):
         """
@@ -130,8 +130,8 @@
         from shallow_water import Domain
         from Numeric import zeros, Float
-        
+
         #Create basic mesh
         points, vertices, boundary = rectangular(1, 3)
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -139,5 +139,5 @@
                                                  'top':[4,5],
                                                  'all':[0,1,2,3,4,5]})
-      
+
         #Set friction
         manning = 0.07
@@ -145,5 +145,5 @@
 
         domain.set_region(Add_value_to_region('bottom', 'friction', 1.0,initial_quantity='friction', location = 'unique vertices'))
-        
+
         #print domain.quantities['friction'].get_values()
         assert allclose(domain.quantities['friction'].get_values(),\
@@ -156,5 +156,5 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Region,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_sparse.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_sparse.py	(revision 1018)
@@ -5,11 +5,11 @@
 
 from sparse import *
-from Numeric import allclose, array, transpose, Float 
-        
-class TestCase(unittest.TestCase):
+from Numeric import allclose, array, transpose, Float
+
+class Test_Sparse(unittest.TestCase):
 
     def setUp(self):
         pass
-        
+
     def tearDown(self):
         pass
@@ -17,5 +17,5 @@
     def test_init1(Self):
         """Test initialisation from dimensions
-        """        
+        """
         A = Sparse(3,3)
         A[1,1] = 4
@@ -26,10 +26,10 @@
                     assert A[i,j] == 4.0
                 else:
-                    assert A[i,j] == 0.0                    
-                    
+                    assert A[i,j] == 0.0
+
     def test_init2(Self):
         """Test initialisation from dense matrix
         """
-        
+
         A = Sparse(4,3)
         A[1,1] = 4
@@ -37,17 +37,17 @@
         A[0,1] = 13
         A[0,2] = -6
-        A[2,2] = 1                        
+        A[2,2] = 1
 
         B = A.todense()
 
         C = Sparse(B)
-        
+
         assert allclose(C.todense(), B)
 
-                    
+
     def test_dense(self):
         A = Sparse(4,3)
-        A[1,1] = 4        
-                     
+        A[1,1] = 4
+
         assert allclose(A.todense(), [[0,0,0], [0,4,0], [0,0,0], [0,0,0]])
 
@@ -59,5 +59,5 @@
 
         A = Sparse(3,3)
-        A[1,1] = 4        
+        A[1,1] = 4
         A[1,1] = 0
 
@@ -68,9 +68,9 @@
         A[1,2] = 0
         assert len(A) == 0
-        assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])        
+        assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
 
     def test_sparse_multiplication_vector(self):
         A = Sparse(3,3)
-        
+
         A[0,0] = 3
         A[1,1] = 2
@@ -91,10 +91,10 @@
 
         u = A*v[:,1]
-        assert allclose(u, [12,16,4])        
+        assert allclose(u, [12,16,4])
 
 
     def test_sparse_multiplication_matrix(self):
         A = Sparse(3,3)
-        
+
         A[0,0] = 3
         A[1,1] = 2
@@ -106,11 +106,11 @@
 
         u = A*v
-        assert allclose(u, [[6,12], [14,16], [4,4]])  
-
-
-
-    def test_sparse_transpose_multiplication(self):        
-        A = Sparse(3,3)
-        
+        assert allclose(u, [[6,12], [14,16], [4,4]])
+
+
+
+    def test_sparse_transpose_multiplication(self):
+        A = Sparse(3,3)
+
         A[0,0] = 3
         A[1,1] = 2
@@ -128,5 +128,5 @@
         """Test method __rmul__
         """
-        
+
         A = Sparse(3,3)
 
@@ -140,5 +140,5 @@
 
         B = A*3
-        assert allclose(B.todense(), 3*A.todense())        
+        assert allclose(B.todense(), 3*A.todense())
 
         try:
@@ -152,5 +152,5 @@
         """ Test sparse addition with dok format
         """
-        
+
         A = Sparse(3,3)
 
@@ -171,5 +171,5 @@
         """ Test conversion to csr format
         """
-        
+
         A = Sparse(4,3)
 
@@ -181,5 +181,5 @@
         A[2,0] = 5
 
-        #print ' ' 
+        #print ' '
         #print A.todense()
 
@@ -197,11 +197,10 @@
 
         assert allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
-        
+
 
 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase, 'test')
+    suite = unittest.makeSuite(Test_Sparse, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
-
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 1017)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 1018)
@@ -13,5 +13,5 @@
     return x+y
 
-class TestCase(unittest.TestCase):
+class Test_Util(unittest.TestCase):
     def setUp(self):
         pass
@@ -25,5 +25,5 @@
         x1 = 1.0; y1 = 0.0; z1 = -1.0
         x2 = 0.0; y2 = 1.0; z2 = 0.0
-        
+
         zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
 
@@ -44,11 +44,11 @@
         q1 = 8.0+2.0/3
         q2 = 2.0+8.0/3
-        
-        #Gradient of fitted pwl surface    
+
+        #Gradient of fitted pwl surface
         a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        
+
         assert abs(a - 3.0) < epsilon
         assert abs(b - 1.0) < epsilon
-    
+
 
 
@@ -67,5 +67,5 @@
                 import util_ext
 
-            
+
     def test_gradient_C_extension(self):
         from util_ext import gradient as gradient_c
@@ -78,15 +78,15 @@
         q1 = 8.0+2.0/3
         q2 = 2.0+8.0/3
-        
-        #Gradient of fitted pwl surface    
+
+        #Gradient of fitted pwl surface
         a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        
+
         assert abs(a - 3.0) < epsilon
         assert abs(b - 1.0) < epsilon
 
-        
+
     def test_gradient_C_extension3(self):
         from util_ext import gradient as gradient_c
-        
+
         from RandomArray import uniform, seed
         seed(17, 53)
@@ -96,10 +96,10 @@
         q0 = uniform(0.0, 10.0, 4)
         q1 = uniform(1.0, 3.0, 4)
-        q2 = uniform(7.0, 20.0, 4)        
-        
+        q2 = uniform(7.0, 20.0, 4)
+
 
         for i in range(4):
             #Gradient of fitted pwl surface
-            from util import gradient_python	        
+            from util import gradient_python
             a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
                                     q0[i], q1[i], q2[i])
@@ -112,8 +112,8 @@
             assert abs(a - a_ref) < epsilon
             assert abs(b - b_ref) < epsilon
-        
-
-
-    #Geometric        
+
+
+
+    #Geometric
     #def test_distance(self):
     #    from util import distance#
@@ -127,5 +127,5 @@
     #
     #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
-    #                    'distance is wrong!')        
+    #                    'distance is wrong!')
 
     def test_angle(self):
@@ -136,5 +136,5 @@
 
     def test_anglediff(self):
-        from util import anglediff        
+        from util import anglediff
 
         assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
@@ -163,5 +163,5 @@
             fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
             t += dt
-    
+
         fid.close()
 
@@ -174,13 +174,13 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t) 
+            assert allclose(q[0], 2*t)
             if i%6 == 0:
                 assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))                
+                assert allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
 
         t = 90 #Halfway between 60 and 120
-        q = F(t) 
+        q = F(t)
         assert allclose( (120**2 + 60**2)/2, q[1] )
         assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
@@ -188,10 +188,10 @@
 
         t = 100 #Two thirds of the way between between 60 and 120
-        q = F(t) 
+        q = F(t)
         assert allclose( 2*120**2/3 + 60**2/3, q[1] )
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
-        
+
         os.remove(filename)
-        
+
 
 
@@ -215,5 +215,5 @@
 
         #Nice test that may render some of the others redundant.
-        
+
         import os, time
         from config import time_format
@@ -231,5 +231,5 @@
                 rectangular(4, 4, 15, 30, origin = (0, -20))
 
-        
+
         #print 'Number of elements', len(vertices)
         domain = Domain(points, vertices, boundary)
@@ -246,5 +246,5 @@
         domain.starttime = start
 
-        
+
         #Store structure
         domain.initialise_storage()
@@ -262,5 +262,5 @@
 
             f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
-            domain.set_quantity('ymomentum', f3)                        
+            domain.set_quantity('ymomentum', f3)
 
             #Store and advance time
@@ -275,5 +275,5 @@
 
         #Set domain.starttime to too early
-        domain.starttime = start - 1                
+        domain.starttime = start - 1
 
         #Create file function
@@ -286,38 +286,38 @@
 
         #Check that domain.starttime isn't updated if later
-        domain.starttime = start + 1                
+        domain.starttime = start + 1
         F = file_function(filename + '.sww', domain,
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+1)                
+        assert allclose(domain.starttime, start+1)
 
         domain.starttime = start
 
 
-        
+
         #Check linear interpolation in time
         #for id in range(len(interpolation_points)):
         for id in [1]:
             x = interpolation_points[id][0]
-            y = interpolation_points[id][1]            
-            
+            y = interpolation_points[id][1]
+
             for i in range(20):
                 t = i*10
                 k = i%6
-            
+
                 if k == 0:
-                    q0 = F(t, point_id=id)                                
-                    q1 = F(t+60, point_id=id)                                
-            
-            
+                    q0 = F(t, point_id=id)
+                    q1 = F(t+60, point_id=id)
+
+
                 q = F(t, point_id=id)
                 assert allclose(q, (k*q1 + (6-k)*q0)/6)
-            
- 
+
+
         #Another check of linear interpolation in time
         for id in range(len(interpolation_points)):
             q60 = F(60, point_id=id)
-            q120 = F(120, point_id=id)        
-        
+            q120 = F(120, point_id=id)
+
             t = 90 #Halfway between 60 and 120
             q = F(t,point_id=id)
@@ -333,32 +333,32 @@
         #than file end time
 	delta = 23
-        domain.starttime = start + delta                
+        domain.starttime = start + delta
         F = file_function(filename + '.sww', domain,
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+delta)                
-
-        
-	
-       
+        assert allclose(domain.starttime, start+delta)
+
+
+
+
         #Now try interpolation with delta offset
         for id in [1]:
             x = interpolation_points[id][0]
-            y = interpolation_points[id][1]            
-            
+            y = interpolation_points[id][1]
+
             for i in range(20):
                 t = i*10
                 k = i%6
-            
+
                 if k == 0:
-                    q0 = F(t-delta, point_id=id)      
-                    q1 = F(t+60-delta, point_id=id)                            
-            
+                    q0 = F(t-delta, point_id=id)
+                    q1 = F(t+60-delta, point_id=id)
+
                 q = F(t-delta, point_id=id)
                 assert allclose(q, (k*q1 + (6-k)*q0)/6)
-            
-        
+
+
         os.remove(filename + '.sww')
-        
+
 
 
@@ -385,5 +385,5 @@
             fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
             t += dt
-   
+
         fid.close()
 
@@ -391,26 +391,26 @@
         b = [4.0, 0.0]
         c = [0.0, 3.0]
-       
+
         points = [a, b, c]
         vertices = [[0,1,2]]
-        domain = Domain(points, vertices) 
+        domain = Domain(points, vertices)
 
         #Check that domain.starttime is updated if non-existing
         F = file_function(filename, domain)
-        assert allclose(domain.starttime, start)        
+        assert allclose(domain.starttime, start)
 
         #Check that domain.starttime is updated if too early
-        domain.starttime = start - 1                
+        domain.starttime = start - 1
         F = file_function(filename, domain)
         assert allclose(domain.starttime, start)
 
         #Check that domain.starttime isn't updated if later
-        domain.starttime = start + 1                
+        domain.starttime = start + 1
         F = file_function(filename, domain)
-        assert allclose(domain.starttime, start+1)                
+        assert allclose(domain.starttime, start+1)
 
         domain.starttime = start
-        
-       
+
+
         #Now try interpolation
         for i in range(20):
@@ -419,13 +419,13 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t) 
+            assert allclose(q[0], 2*t)
             if i%6 == 0:
                 assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))                
+                assert allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
 
         t = 90 #Halfway between 60 and 120
-        q = F(t) 
+        q = F(t)
         assert allclose( (120**2 + 60**2)/2, q[1] )
         assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
@@ -433,10 +433,10 @@
 
         t = 100 #Two thirds of the way between between 60 and 120
-        q = F(t) 
+        q = F(t)
         assert allclose( 2*120**2/3 + 60**2/3, q[1] )
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
-        
+
         os.remove(filename)
-        
+
 
     def test_file_function_time_with_domain_different_start(self):
@@ -464,5 +464,5 @@
             fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
             t += dt
-   
+
         fid.close()
 
@@ -470,20 +470,20 @@
         b = [4.0, 0.0]
         c = [0.0, 3.0]
-       
+
         points = [a, b, c]
         vertices = [[0,1,2]]
-        domain = Domain(points, vertices) 
+        domain = Domain(points, vertices)
 
         #Check that domain.starttime isn't updated if later than file starttime but earlier
         #than file end time
 	delta = 23
-        domain.starttime = start + delta                
+        domain.starttime = start + delta
         F = file_function(filename, domain)
-        assert allclose(domain.starttime, start+delta)                
-
-        
-	
-       
-        #Now try interpolation with delta offset 
+        assert allclose(domain.starttime, start+delta)
+
+
+
+
+        #Now try interpolation with delta offset
         for i in range(20):
             t = i*10
@@ -491,13 +491,13 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t) 
+            assert allclose(q[0], 2*t)
             if i%6 == 0:
                 assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))                
+                assert allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
 
         t = 90 #Halfway between 60 and 120
-        q = F(t-delta) 
+        q = F(t-delta)
         assert allclose( (120**2 + 60**2)/2, q[1] )
         assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
@@ -505,8 +505,8 @@
 
         t = 100 #Two thirds of the way between between 60 and 120
-        q = F(t-delta) 
+        q = F(t-delta)
         assert allclose( 2*120**2/3 + 60**2/3, q[1] )
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
-        
+
         os.remove(filename)
 
@@ -517,34 +517,34 @@
         """
         import time
-       
-	#Create sww file of simple propagation from left to right 
+
+	#Create sww file of simple propagation from left to right
 	#through rectangular domain
 	from shallow_water import Domain, Dirichlet_boundary
         from mesh_factory import rectangular
-        from Numeric import take, concatenate, reshape        
+        from Numeric import take, concatenate, reshape
 
         #Create basic mesh and shallow water domain
         points, vertices, boundary = rectangular(3, 3)
         domain1 = Domain(points, vertices, boundary)
-	
+
         from util import mean
-        domain1.reduction = mean 
+        domain1.reduction = mean
         domain1.smooth = True #NOTE: Mimic sww output where each vertex has
                               # only one value.
 
         domain1.default_order = 2
-	domain1.store = True	
+	domain1.store = True
         domain1.set_datadir('.')
         domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
-		
+
         #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
         domain1.set_quantity('elevation', 0)
         domain1.set_quantity('friction', 0)
-        domain1.set_quantity('stage', 0)        
+        domain1.set_quantity('stage', 0)
 
         # Boundary conditions
-        B0 = Dirichlet_boundary([0,0,0])	        
-        B6 = Dirichlet_boundary([0.6,0,0])	
+        B0 = Dirichlet_boundary([0,0,0])
+        B6 = Dirichlet_boundary([0.6,0,0])
         domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
         domain1.check_integrity()
@@ -560,8 +560,8 @@
         from Scientific.IO.NetCDF import NetCDFFile
         filename = domain1.get_name() + '.' + domain1.format
-        fid = NetCDFFile(filename)  
+        fid = NetCDFFile(filename)
 
         x = fid.variables['x'][:]
-        y = fid.variables['y'][:]        
+        y = fid.variables['y'][:]
         stage = fid.variables['stage'][:]
         xmomentum = fid.variables['xmomentum'][:]
@@ -581,10 +581,10 @@
         #this timestep are
         r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2         
-        r2 = (D[2] + D[3])/2         
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
 
         #And the midpoints are found now
         Dx = take(reshape(x, (16,1)), [0,5,10,15])
-        Dy = take(reshape(y, (16,1)), [0,5,10,15])                
+        Dy = take(reshape(y, (16,1)), [0,5,10,15])
 
         diag = concatenate( (Dx, Dy), axis=1)
@@ -597,5 +597,5 @@
 
         q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)        
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
         q = f(timestep/10., point_id=2); assert allclose(r2, q)
 
@@ -613,9 +613,9 @@
         #this timestep are
         r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2         
-        r2 = (D[2] + D[3])/2         
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
 
         #Let us see if the file function can find the correct
-        #values 
+        #values
         q = f(0, point_id=0); assert allclose(r0, q)
         q = f(0, point_id=1); assert allclose(r1, q)
@@ -626,5 +626,5 @@
         #Now do it again for a timestep in the middle
 
-        timestep = 33 
+        timestep = 33
         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
@@ -635,9 +635,9 @@
         #this timestep are
         r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2         
-        r2 = (D[2] + D[3])/2         
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
 
         q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)        
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
         q = f(timestep/10., point_id=2); assert allclose(r2, q)
 
@@ -656,6 +656,6 @@
         #this timestep are
         r0_0 = (D[0] + D[1])/2
-        r1_0 = (D[1] + D[2])/2         
-        r2_0 = (D[2] + D[3])/2         
+        r1_0 = (D[1] + D[2])/2
+        r2_0 = (D[2] + D[3])/2
 
         #
@@ -669,14 +669,14 @@
         #this timestep are
         r0_1 = (D[0] + D[1])/2
-        r1_1 = (D[1] + D[2])/2         
-        r2_1 = (D[2] + D[3])/2         
+        r1_1 = (D[1] + D[2])/2
+        r2_1 = (D[2] + D[3])/2
 
         # The reference values are
         r0 = (r0_0 + r0_1)/2
         r1 = (r1_0 + r1_1)/2
-        r2 = (r2_0 + r2_1)/2        
+        r2 = (r2_0 + r2_1)/2
 
         q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)        
+        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
         q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
 
@@ -688,15 +688,15 @@
         r0 = (r0_0 + 2*r0_1)/3
         r1 = (r1_0 + 2*r1_1)/3
-        r2 = (r2_0 + 2*r2_1)/3        
+        r2 = (r2_0 + 2*r2_1)/3
 
         #And the file function gives
         q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)   
+        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
         q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
-        
+
         fid.close()
         import os
         os.remove(filename)
-        
+
 
     def test_xya_ascii(self):
@@ -704,10 +704,10 @@
         FN = 'xyatest' + str(time.time()) + '.xya'
     	fid = open(FN, 'w')
-	fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )        
+	fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
 	fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
 	fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
-	fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )		
+	fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
 	fid.close()
-	
+
 	points, attributes = read_xya(FN, format = 'asc')
 
@@ -715,18 +715,18 @@
 	assert allclose(attributes['a1'], [10,30,40.2])
 	assert allclose(attributes['a2'], [20,20,40.3])
-	assert allclose(attributes['a3'], [30,10,40.4])        
-	
+	assert allclose(attributes['a3'], [30,10,40.4])
+
 	os.remove(FN)
-	
+
     def test_xya_ascii_w_names(self):
     	import time, os
         FN = 'xyatest' + str(time.time()) + '.xya'
     	fid = open(FN, 'w')
-	fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )        
+	fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
 	fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
 	fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
-	fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )		
+	fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
 	fid.close()
-	
+
 	points, attributes = read_xya(FN, format = 'asc')
 
@@ -735,7 +735,7 @@
 	assert allclose(attributes['a1'], [10,30,40.2])
 	assert allclose(attributes['a2'], [20,20,40.3])
-	assert allclose(attributes['a3'], [30,10,40.4])        
-
-	
+	assert allclose(attributes['a3'], [30,10,40.4])
+
+
 	os.remove(FN)
 
@@ -743,178 +743,178 @@
 
 
-    #Polygon stuff	
+    #Polygon stuff
     def test_polygon_function_constants(self):
         p1 = [[0,0], [10,0], [10,10], [0,10]]
         p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-	
-	f = Polygon_function( [(p1, 1.0)] )	
+
+	f = Polygon_function( [(p1, 1.0)] )
 	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
 	assert allclose(z, [1,1,0,0])
-		
-	
+
+
 	f = Polygon_function( [(p2, 2.0)] )
 	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
-	assert allclose(z, [2,0,0,2])	
-		
-	
-	#Combined	
+	assert allclose(z, [2,0,0,2])
+
+
+	#Combined
 	f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
-	z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
-	assert allclose(z, [2,1,0,2])	
-	
-	
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,1,0,2])
+
+
     def test_polygon_function_callable(self):
-        """Check that values passed into Polygon_function can be callable 
+        """Check that values passed into Polygon_function can be callable
 	themselves.
 	"""
         p1 = [[0,0], [10,0], [10,10], [0,10]]
         p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-	
-	f = Polygon_function( [(p1, test_function)] )	
+
+	f = Polygon_function( [(p1, test_function)] )
 	z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
 	assert allclose(z, [10,14,0,0])
-		
-	#Combined	
+
+	#Combined
 	f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
-	z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
-	assert allclose(z, [2,14,0,2])	
-	
-	
-	#Combined w default	
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,0,2])
+
+
+	#Combined w default
 	f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
-	z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
-	assert allclose(z, [2,14,3.14,2])		
-	
-	
-	#Combined w default func	
-	f = Polygon_function( [(p1, test_function), (p2, 2.0)], 
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,3.14,2])
+
+
+	#Combined w default func
+	f = Polygon_function( [(p1, test_function), (p2, 2.0)],
 			      default = test_function)
-	z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
-	assert allclose(z, [2,14,35,2])			
-	
-	
+	z = f([5, 5, 27, 35], [5, 9, 8, -5])
+	assert allclose(z, [2,14,35,2])
+
+
     def test_point_on_line(self):
-            	
-	#Endpoints first	
-	assert point_on_line( 0, 0, 0,0, 1,0 ) 	
-	assert point_on_line( 1, 0, 0,0, 1,0 ) 		
-
-	#Then points on line	
-	assert point_on_line( 0.5, 0, 0,0, 1,0 ) 
-	assert point_on_line( 0, 0.5, 0,1, 0,0 ) 		
-	assert point_on_line( 1, 0.5, 1,1, 1,0 ) 	
-	assert point_on_line( 0.5, 0.5, 0,0, 1,1 ) 		
-	
-	#Then points not on line		
-	assert not point_on_line( 0.5, 0, 0,1, 1,1 ) 	
-	assert not point_on_line( 0, 0.5, 0,0, 1,1 ) 
-	
+
+	#Endpoints first
+	assert point_on_line( 0, 0, 0,0, 1,0 )
+	assert point_on_line( 1, 0, 0,0, 1,0 )
+
+	#Then points on line
+	assert point_on_line( 0.5, 0, 0,0, 1,0 )
+	assert point_on_line( 0, 0.5, 0,1, 0,0 )
+	assert point_on_line( 1, 0.5, 1,1, 1,0 )
+	assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
+
+	#Then points not on line
+	assert not point_on_line( 0.5, 0, 0,1, 1,1 )
+	assert not point_on_line( 0, 0.5, 0,0, 1,1 )
+
 	#From real example that failed
 	assert not point_on_line( 40,50, 40,20, 40,40 )
-	
-	
+
+
 	#From real example that failed
-	assert not point_on_line( 40,19, 40,20, 40,40 )	
-	 		
-	
-	
-	
+	assert not point_on_line( 40,19, 40,20, 40,40 )
+
+
+
+
     def test_inside_polygon_main(self):
 
-        
+
         #Simplest case: Polygon is the unit square
         polygon = [[0,0], [1,0], [1,1], [0,1]]
 
 	assert inside_polygon( (0.5, 0.5), polygon )
-	assert not inside_polygon( (0.5, 1.5), polygon )   
-	assert not inside_polygon( (0.5, -0.5), polygon )   	
-	assert not inside_polygon( (-0.5, 0.5), polygon )   		
-	assert not inside_polygon( (1.5, 0.5), polygon )   	
-	
+	assert not inside_polygon( (0.5, 1.5), polygon )
+	assert not inside_polygon( (0.5, -0.5), polygon )
+	assert not inside_polygon( (-0.5, 0.5), polygon )
+	assert not inside_polygon( (1.5, 0.5), polygon )
+
 	#Try point on borders
-	assert inside_polygon( (1., 0.5), polygon, closed=True)	
-	assert inside_polygon( (0.5, 1), polygon, closed=True)	
-	assert inside_polygon( (0., 0.5), polygon, closed=True)	
-	assert inside_polygon( (0.5, 0.), polygon, closed=True)	
-	
-	assert not inside_polygon( (0.5, 1), polygon, closed=False)	
-	assert not inside_polygon( (0., 0.5), polygon, closed=False)	
-	assert not inside_polygon( (0.5, 0.), polygon, closed=False)	
-	assert not inside_polygon( (1., 0.5), polygon, closed=False)	
-
-
-    
+	assert inside_polygon( (1., 0.5), polygon, closed=True)
+	assert inside_polygon( (0.5, 1), polygon, closed=True)
+	assert inside_polygon( (0., 0.5), polygon, closed=True)
+	assert inside_polygon( (0.5, 0.), polygon, closed=True)
+
+	assert not inside_polygon( (0.5, 1), polygon, closed=False)
+	assert not inside_polygon( (0., 0.5), polygon, closed=False)
+	assert not inside_polygon( (0.5, 0.), polygon, closed=False)
+	assert not inside_polygon( (1., 0.5), polygon, closed=False)
+
+
+
         #From real example (that failed)
-	polygon = [[20,20], [40,20], [40,40], [20,40]]         
+	polygon = [[20,20], [40,20], [40,40], [20,40]]
 	points = [ [40, 50] ]
 	res = inside_polygon(points, polygon)
 	assert len(res) == 0
-	
-	polygon = [[20,20], [40,20], [40,40], [20,40]]         
+
+	polygon = [[20,20], [40,20], [40,40], [20,40]]
         points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
 	res = inside_polygon(points, polygon)
 	assert len(res) == 2
 	assert allclose(res, [0,1])
-        
-		
-	      
+
+
+
 	#More convoluted and non convex polygon
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
 	assert inside_polygon( (0.5, 0.5), polygon )
-	assert inside_polygon( (1, -0.5), polygon ) 
+	assert inside_polygon( (1, -0.5), polygon )
 	assert inside_polygon( (1.5, 0), polygon )
-	
-	assert not inside_polygon( (0.5, 1.5), polygon )	
-	assert not inside_polygon( (0.5, -0.5), polygon )		
-	
-	
+
+	assert not inside_polygon( (0.5, 1.5), polygon )
+	assert not inside_polygon( (0.5, -0.5), polygon )
+
+
 	#Very convoluted polygon
         polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-	assert inside_polygon( (5, 5), polygon )	
-	assert inside_polygon( (17, 7), polygon )		
-	assert inside_polygon( (27, 2), polygon )			
-	assert inside_polygon( (35, -5), polygon )			
+	assert inside_polygon( (5, 5), polygon )
+	assert inside_polygon( (17, 7), polygon )
+	assert inside_polygon( (27, 2), polygon )
+	assert inside_polygon( (35, -5), polygon )
 	assert not inside_polygon( (15, 7), polygon )
 	assert not inside_polygon( (24, 3), polygon )
-	assert not inside_polygon( (25, -10), polygon )	
-	
-	
-	
+	assert not inside_polygon( (25, -10), polygon )
+
+
+
 	#Another combination (that failed)
         polygon = [[0,0], [10,0], [10,10], [0,10]]
-	assert inside_polygon( (5, 5), polygon )	
-	assert inside_polygon( (7, 7), polygon )	
-	assert not inside_polygon( (-17, 7), polygon )		      	
-	assert not inside_polygon( (7, 17), polygon )			
-	assert not inside_polygon( (17, 7), polygon )			
-	assert not inside_polygon( (27, 8), polygon )			
-	assert not inside_polygon( (35, -5), polygon )			
-	
-	
-			
-	
-	#Multiple polygons	
-	
+	assert inside_polygon( (5, 5), polygon )
+	assert inside_polygon( (7, 7), polygon )
+	assert not inside_polygon( (-17, 7), polygon )
+	assert not inside_polygon( (7, 17), polygon )
+	assert not inside_polygon( (17, 7), polygon )
+	assert not inside_polygon( (27, 8), polygon )
+	assert not inside_polygon( (35, -5), polygon )
+
+
+
+
+	#Multiple polygons
+
         polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
 		   [10,10], [11,10], [11,11], [10,11], [10,10]]
-        assert inside_polygon( (0.5, 0.5), polygon )		   
-        assert inside_polygon( (10.5, 10.5), polygon )		   	
-	
+        assert inside_polygon( (0.5, 0.5), polygon )
+        assert inside_polygon( (10.5, 10.5), polygon )
+
 	#FIXME: Fails if point is 5.5, 5.5
-        assert not inside_polygon( (0, 5.5), polygon )		   		
-
-	#Polygon with a hole	
+        assert not inside_polygon( (0, 5.5), polygon )
+
+	#Polygon with a hole
         polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
 	           [0,0], [1,0], [1,1], [0,1], [0,0]]
-		   	
+
         assert inside_polygon( (0, -0.5), polygon )
-        assert not inside_polygon( (0.5, 0.5), polygon )	
-
-    def test_inside_polygon_vector_version(self):	
+        assert not inside_polygon( (0.5, 0.5), polygon )
+
+    def test_inside_polygon_vector_version(self):
 	#Now try the vector formulation returning indices
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]        
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
 	points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
 	res = inside_polygon( points, polygon )
-	
+
 	assert allclose( res, [0,1,2] )
 
@@ -938,16 +938,15 @@
         for point in points:
             assert inside_polygon(point, polygon)
-            
-
-		      
+
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
     #suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Util,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
-    
-    
-
-
-
+
+
+
+
