Index: trunk/anuga_core/source/anuga/__init__.py
===================================================================
--- trunk/anuga_core/source/anuga/__init__.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/__init__.py	(revision 7841)
@@ -30,7 +30,9 @@
 from anuga.shallow_water.shallow_water_domain import Domain
 
-from anuga.abstract_2d_finite_volumes.util import file_function, sww2timeseries
-
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.util import file_function, \
+                                        sww2timeseries, sww2csv_gauges
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross, \
+                                                    rectangular
 
 from anuga.file.csv_file import load_csv_as_building_polygons,  \
@@ -47,5 +49,5 @@
 
 from anuga.utilities.system_tools import file_length
-
+from anuga.utilities.sww_merge import sww_merge
 from anuga.utilities.file_utils import copy_code_files
 
@@ -189,5 +191,7 @@
     throw an error
     
-    Interior_holes is a list of ploygons for each hole.
+    interior_holes is a list of ploygons for each hole. These polygons do not
+    need to be closed, but their points must be specified in a counter-clockwise
+    order.
 
     This function does not allow segments to share points - use underlying
Index: trunk/anuga_core/source/anuga/config.py
===================================================================
--- trunk/anuga_core/source/anuga/config.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/config.py	(revision 7841)
@@ -34,5 +34,5 @@
 # Major revision number for use with create_distribution
 # and update_anuga_user_guide
-major_revision = '1.1beta'
+major_revision = '1.2'
 
 ################################################################################
Index: trunk/anuga_core/source/anuga/file/sww.py
===================================================================
--- trunk/anuga_core/source/anuga/file/sww.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/file/sww.py	(revision 7841)
@@ -47,146 +47,4 @@
 
         #FIXME: Should we have a general set_precision function?
-
-
-##
-# @brief Class for handling checkpoints data
-# @note This is not operational at the moment
-class CPT_file(Data_format):
-    """Interface to native NetCDF format (.cpt) to be 
-    used for checkpointing (one day)
-    """
-
-    ##
-    # @brief Initialize this instantiation.
-    # @param domain ??
-    # @param mode Mode of underlying data file (default WRITE).
-    def __init__(self, domain, mode=netcdf_mode_w):
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        self.precision = netcdf_float #Use full precision
-
-        Data_format.__init__(self, domain, 'sww', mode)
-
-        # NetCDF file definition
-        fid = NetCDFFile(self.filename, mode)
-        if mode[0] == 'w':
-            # Create new file
-            fid.institution = 'Geoscience Australia'
-            fid.description = 'Checkpoint data'
-            #fid.smooth = domain.smooth
-            fid.order = domain.default_order
-
-            # Dimension definitions
-            fid.createDimension('number_of_volumes', self.number_of_volumes)
-            fid.createDimension('number_of_vertices', 3)
-
-            # Store info at all vertices (no smoothing)
-            fid.createDimension('number_of_points', 3*self.number_of_volumes)
-            fid.createDimension('number_of_timesteps', None) #extensible
-
-            # Variable definitions
-
-            # Mesh
-            fid.createVariable('x', self.precision, ('number_of_points',))
-            fid.createVariable('y', self.precision, ('number_of_points',))
-
-
-            fid.createVariable('volumes', netcdf_int, ('number_of_volumes',
-                                                       'number_of_vertices'))
-
-            fid.createVariable('time', self.precision, ('number_of_timesteps',))
-
-            #Allocate space for all quantities
-            for name in domain.quantities.keys():
-                fid.createVariable(name, self.precision,
-                                   ('number_of_timesteps',
-                                    'number_of_points'))
-
-        fid.close()
-
-    ##
-    # @brief Store connectivity data to underlying data file.
-    def store_checkpoint(self):
-        """Write x,y coordinates of triangles.
-        Write connectivity (
-        constituting
-        the bed elevation.
-        """
-
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        domain = self.domain
-
-        #Get NetCDF
-        fid = NetCDFFile(self.filename, netcdf_mode_a)
-
-        # Get the variables
-        x = fid.variables['x']
-        y = fid.variables['y']
-
-        volumes = fid.variables['volumes']
-
-        # Get X, Y and bed elevation Z
-        Q = domain.quantities['elevation']
-        X,Y,Z,V = Q.get_vertex_values(xy=True, precision=self.precision)
-
-        x[:] = X.astype(self.precision)
-        y[:] = Y.astype(self.precision)
-        z[:] = Z.astype(self.precision)
-
-        volumes[:] = V
-
-        fid.close()
-
-    ##
-    # @brief Store time and named quantities to underlying data file.
-    # @param name 
-    def store_timestep(self, name):
-        """Store time and named quantity to file
-        """
-
-        from Scientific.IO.NetCDF import NetCDFFile
-        from time import sleep
-
-        #Get NetCDF
-        retries = 0
-        file_open = False
-        while not file_open and retries < 10:
-            try:
-                fid = NetCDFFile(self.filename, netcdf_mode_a)
-            except IOError:
-                #This could happen if someone was reading the file.
-                #In that case, wait a while and try again
-                msg = 'Warning (store_timestep): File %s could not be opened' \
-                      ' - trying again' % self.filename
-                log.critical(msg)
-                retries += 1
-                sleep(1)
-            else:
-                file_open = True
-
-        if not file_open:
-            msg = 'File %s could not be opened for append' % self.filename
-            raise DataFileNotOpenError, msg
-
-        domain = self.domain
-
-        # Get the variables
-        time = fid.variables['time']
-        stage = fid.variables['stage']
-        i = len(time)
-
-        #Store stage
-        time[i] = self.domain.time
-
-        # Get quantity
-        Q = domain.quantities[name]
-        A,V = Q.get_vertex_values(xy=False, precision=self.precision)
-
-        stage[i,:] = A.astype(self.precision)
-
-        #Flush and close
-        fid.sync()
-        fid.close()
 
 
@@ -536,4 +394,6 @@
 
     def read_mesh(self):
+        """ Read and store the mesh data contained within this sww file.
+        """
         fin = NetCDFFile(self.source, 'r')
 
@@ -551,9 +411,11 @@
 
 
-
         fin.close()
         
     def read_quantities(self, frame_number=0):
-
+        """
+        Read the quantities contained in this file.
+        frame_number is the time index to load.
+        """
         assert frame_number >= 0 and frame_number <= self.last_frame_number
 
@@ -575,21 +437,32 @@
 
     def get_bounds(self):
+        """
+            Get the bounding rect around the mesh.
+        """
         return [self.xmin, self.xmax, self.ymin, self.ymax]
 
     def get_last_frame_number(self):
+        """
+            Return the last loaded frame index.
+        """
         return self.last_frame_number
 
     def get_time(self):
+        """
+            Get time at the current frame num, in secs.
+        """
         return self.time[self.frame_number]
 
 
-# @brief A class to write an SWW file.
 class Write_sww:
-    
+    """
+        A class to write an SWW file.
+        
+        It is domain agnostic, and requires all the data to be fed in
+        manually.
+    """
     RANGE = '_range'
     EXTREMA = ':extrema'
 
-    ##
-    # brief Instantiate the SWW writer class.
     def __init__(self, static_quantities, dynamic_quantities):
         """Initialise Write_sww with two list af quantity names: 
@@ -607,16 +480,4 @@
 
 
-    ##
-    # @brief Store a header in the SWW file.
-    # @param outfile Open handle to the file that will be written.
-    # @param times A list of time slices *or* a start time.
-    # @param number_of_volumes The number of triangles.
-    # @param number_of_points The number of points.
-    # @param description The internal file description string.
-    # @param smoothing True if smoothing is to be used.
-    # @param order 
-    # @param sww_precision Data type of the quantity written (netcdf constant)
-    # @param verbose True if this function is to be verbose.
-    # @note If 'times' is a list, the info will be made relative.
     def store_header(self,
                      outfile,
@@ -631,8 +492,11 @@
         """Write an SWW file header.
 
+        Writes the first section of the .sww file.
+
         outfile - the open file that will be written
         times - A list of the time slice times OR a start time
         Note, if a list is given the info will be made relative.
         number_of_volumes - the number of triangles
+        number_of_points - the number of vertices in the mesh
         """
 
@@ -746,13 +610,5 @@
                          % (num.min(times), num.max(times), len(times.flat)))
 
-    ##
-    # @brief Store triangulation data in the underlying file.
-    # @param outfile Open handle to underlying file.
-    # @param points_utm List or array of points in UTM.
-    # @param volumes 
-    # @param zone 
-    # @param new_origin georeference that the points can be set to.
-    # @param points_georeference The georeference of the points_utm.
-    # @param verbose True if this function is to be verbose.
+
     def store_triangulation(self,
                             outfile,
@@ -764,4 +620,16 @@
                             verbose=False):
         """
+        Store triangulation data in the underlying file.
+        
+        Stores the points and triangle indices in the sww file
+        
+        outfile Open handle to underlying file.
+
+        new_origin georeference that the points can be set to.
+
+        points_georeference The georeference of the points_utm.
+
+        verbose True if this function is to be verbose.
+
         new_origin - qa georeference that the points can be set to. (Maybe
         do this before calling this function.)
@@ -840,10 +708,4 @@
 
 
-
-    # @brief Write the static quantity data to the underlying file.
-    # @param outfile Handle to open underlying file.
-    # @param sww_precision Format of quantity data to write (default Float32).
-    # @param verbose True if this function is to be verbose.
-    # @param **quant
     def store_static_quantities(self, 
                                 outfile, 
@@ -896,12 +758,4 @@
         
         
-    ##
-    # @brief Write the quantity data to the underlying file.
-    # @param outfile Handle to open underlying file.
-    # @param sww_precision Format of quantity data to write (default Float32).
-    # @param slice_index
-    # @param time
-    # @param verbose True if this function is to be verbose.
-    # @param **quant
     def store_quantities(self, 
                          outfile, 
@@ -986,10 +840,6 @@
 
 
-##
-# @brief Get the extents of a NetCDF data file.
-# @param file_name The path to the SWW file.
-# @return A list of x, y, z and stage limits (min, max).
 def extent_sww(file_name):
-    """Read in an sww file.
+    """Read in an sww file, then get its extents
 
     Input:
@@ -1015,18 +865,10 @@
 
 
-##
-# @brief 
-# @param filename
-# @param boundary
-# @param t
-# @param fail_if_NaN
-# @param NaN_filler
-# @param verbose
-# @param very_verbose
-# @return 
 def load_sww_as_domain(filename, boundary=None, t=None,
                fail_if_NaN=True, NaN_filler=0,
                verbose=False, very_verbose=False):
     """
+    Load an sww file into a domain.
+
     Usage: domain = load_sww_as_domain('file.sww',
                         t=time (default = last time in file))
@@ -1194,12 +1036,4 @@
 
 
-##
-# @brief Get mesh and quantity data from an SWW file.
-# @param filename Path to data file to read.
-# @param quantities UNUSED!
-# @param verbose True if this function is to be verbose.
-# @return (mesh, quantities, time) where mesh is the mesh data, quantities is
-#         a dictionary of {name: value}, and time is the time vector.
-# @note Quantities extracted: 'elevation', 'stage', 'xmomentum' and 'ymomentum'
 def get_mesh_and_quantities_from_file(filename,
                                       quantities=None,
@@ -1280,13 +1114,10 @@
 
 
-
-##
-# @brief 
-# @parm time 
-# @param t 
-# @return An (index, ration) tuple.
 def get_time_interp(time, t=None):
-    #Finds the ratio and index for time interpolation.
-    #It is borrowed from previous abstract_2d_finite_volumes code.
+    """Finds the ratio and index for time interpolation.
+        time is an array of time steps
+        t is the sample time.
+        returns a tuple containing index into time, and ratio
+    """
     if t is None:
         t=time[-1]
@@ -1337,10 +1168,7 @@
 
 
-##
-# @brief 
-# @param coordinates 
-# @param volumes 
-# @param boundary 
 def weed(coordinates, volumes, boundary=None):
+    """ Excise all duplicate points.
+    """
     if isinstance(coordinates, num.ndarray):
         coordinates = coordinates.tolist()
Index: trunk/anuga_core/source/anuga/file_conversion/dem2pts.py
===================================================================
--- trunk/anuga_core/source/anuga/file_conversion/dem2pts.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/file_conversion/dem2pts.py	(revision 7841)
@@ -6,16 +6,8 @@
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
                             netcdf_float
+
+from asc2dem import asc2dem
                             
-##
-# @brief Convert DEM data  to PTS data.
-# @param basename_in Stem of input filename.
-# @param basename_out Stem of output filename.
-# @param easting_min 
-# @param easting_max 
-# @param northing_min 
-# @param northing_max 
-# @param use_cache 
-# @param verbose 
-# @return 
+
 def dem2pts(name_in, name_out=None,
             easting_min=None, easting_max=None,
@@ -33,4 +25,6 @@
     NODATA_value  -9999
     138.3698 137.4194 136.5062 135.5558 ..........
+
+    name_in may be a .asc or .dem file to be converted.
 
     Convert to NetCDF pts format which is
@@ -81,11 +75,18 @@
     from Scientific.IO.NetCDF import NetCDFFile
 
-    if name_in[-4:] != '.dem':
-        raise IOError('Input file %s should be of type .dem.' % name_in)
+    root = name_in[:-4]
+
+    if name_in[-4:] == '.asc':
+        intermediate = root + '.dem'
+        if verbose:
+            log.critical('Preconvert %s from asc to %s' % \
+                                    (name_in, intermediate))
+        asc2dem(name_in)
+        name_in = intermediate
+    elif name_in[-4:] != '.dem':
+        raise IOError('Input file %s should be of type .asc or .dem.' % name_in)
 
     if name_out != None and basename_out[-4:] != '.pts':
         raise IOError('Input file %s should be of type .pts.' % name_out)
-
-    root = name_in[:-4]
 
     # Get NetCDF
Index: trunk/anuga_core/source/anuga/file_conversion/sww2dem.py
===================================================================
--- trunk/anuga_core/source/anuga/file_conversion/sww2dem.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/file_conversion/sww2dem.py	(revision 7841)
@@ -124,4 +124,6 @@
     if block_size is None:
         block_size = DEFAULT_BLOCK_SIZE
+
+    assert(isinstance(block_size, (int, long, float)))
 
     # Read sww file
@@ -217,4 +219,9 @@
     result = num.zeros(number_of_points, num.float)
 
+    if verbose:
+        msg = 'Slicing sww file, num points: ' + str(number_of_points)
+        msg += ', block size: ' + str(block_size)
+        log.critical(msg)
+
     for start_slice in xrange(0, number_of_points, block_size):
         # Limit slice size to array end if at last block
@@ -357,5 +364,6 @@
 
         #Write
-        if verbose: log.critical('Writing %s' % demfile)
+        if verbose:
+            log.critical('Writing %s' % name_out)
 
         import ermapper_grids
@@ -437,7 +445,9 @@
                 format='ers'):
     """Wrapper for sww2dem.
-    See sww2dem to find out what most of the parameters do.
+    See sww2dem to find out what most of the parameters do. Note that since this
+    is a batch command, the normal filename naming conventions do not apply.
 
     basename_in is a path to sww file/s, without the .sww extension.
+    extra_name_out is a postfix to add to the output filename.
 
     Quantities is a list of quantities.  Each quantity will be
@@ -474,6 +484,12 @@
                 basename_out = sww_file + '_' + quantity + '_' + extra_name_out
 
-            file_out = sww2dem(dir+os.sep+sww_file+'.sww',
-                               dir+os.sep+basename_out,
+            swwin = dir+os.sep+sww_file+'.sww'
+            demout = dir+os.sep+basename_out+'.'+format
+
+            if verbose:
+                log.critical('sww2dem: %s => %s' % (swwin, demout))
+
+            file_out = sww2dem(swwin,
+                               demout,
                                quantity,
                                reduction,
@@ -487,6 +503,6 @@
                                verbose,
                                origin,
-                               datum,
-                               format)
+                               datum)
+                               
             files_out.append(file_out)
     return files_out
Index: trunk/anuga_core/source/anuga/file_conversion/test_dem2pts.py
===================================================================
--- trunk/anuga_core/source/anuga/file_conversion/test_dem2pts.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/file_conversion/test_dem2pts.py	(revision 7841)
@@ -79,6 +79,5 @@
 
         #Convert to NetCDF pts
-        asc2dem(filename)
-        dem2pts(root+'.dem', easting_min=2002.0, easting_max=2007.0,
+        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
                 northing_min=3003.0, northing_max=3006.0,
                 verbose=False)
@@ -199,6 +198,5 @@
 
         #Convert to NetCDF pts
-        asc2dem(filename)
-        dem2pts(root+'.dem', easting_min=2002.0, easting_max=2007.0,
+        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
                 northing_min=3003.0, northing_max=3006.0)
 
@@ -328,6 +326,5 @@
 
         #Convert to NetCDF pts
-        asc2dem(filename)
-        dem2pts(root+'.dem', easting_min=2002.0, easting_max=2007.0,
+        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
                 northing_min=3003.0, northing_max=3006.0)
 
Index: trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py
===================================================================
--- trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py	(revision 7841)
@@ -20,5 +20,5 @@
 
 # local modules
-from sww2dem import sww2dem
+from sww2dem import sww2dem, sww2dem_batch
 
 class Test_Sww2Dem(unittest.TestCase):
@@ -1392,4 +1392,133 @@
         os.remove(self.domain.get_name() + '_elevation')
         os.remove(self.domain.get_name() + '_elevation.ers')
+        
+    def test_export_grid_parallel(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        base_name = 'tegp'
+        #Setup
+        self.domain.set_name(base_name+'_P0_8')
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.0001)
+        #Setup
+        self.domain.set_name(base_name+'_P1_8')
+        swwfile2 = self.domain.get_name() + '.sww'
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.0002)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        fid.close()
+
+        #Export to ascii/prj files
+        extra_name_out = 'yeah'
+        sww2dem_batch(base_name,
+                    quantities = ['elevation', 'depth'],
+                    extra_name_out = extra_name_out,
+                    cellsize = cellsize,
+                    verbose = self.verbose,
+                    format = 'asc')
+
+        prjfile = base_name + '_P0_8_elevation_yeah.prj'
+        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)               
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        prjfile = base_name + '_P1_8_elevation_yeah.prj'
+        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)               
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+        #Check asc file
+        ascfile = base_name + '_P0_8_depth_yeah.asc'
+        prjfile = base_name + '_P0_8_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        #Check asc file
+        ascfile = base_name + '_P1_8_depth_yeah.asc'
+        prjfile = base_name + '_P1_8_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile2)
+        
 
 #################################################################################
Index: trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 7841)
@@ -1,6 +1,13 @@
 """Least squares interpolation.
 
-   Implements a least-squares interpolation.
-   Putting mesh data onto points.
+   These functions and classes calculate a value at a particular point on
+   the given mesh. It interpolates the values stored at the vertices of the
+   mesh.
+   
+   For example, if you want to get the height of a terrain mesh at particular
+   point, you pass the point to an Interpolate class. The point will intersect
+   one of the triangles on the mesh, and the interpolated height will be an
+   intermediate value between the three vertices of that triangle.
+   This value is returned by the class.
 
    Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
@@ -896,16 +903,13 @@
                     if verbose is True:
                         import sys
-                        if sys.platform == 'win32':
-                            # FIXME (Ole): Why only Windoze?
-                            from anuga.geometry.polygon import plot_polygons
-                            title = ('Interpolation points fall '
-                                     'outside specified mesh')
-                            plot_polygons([mesh_boundary_polygon,
-                                           interpolation_points,
-                                           out_interp_pts],
-                                          ['line', 'point', 'outside'],
-                                          figname='points_boundary_out',
-                                          label=title,
-                                          verbose=verbose)
+                        from anuga.geometry.polygon import plot_polygons
+                        title = ('Interpolation points fall '
+                                 'outside specified mesh')
+                        plot_polygons([mesh_boundary_polygon,
+                                       interpolation_points,
+                                       out_interp_pts],
+                                      ['line', 'point', 'outside'],
+                                      figname='points_boundary_out',
+                                      label=title)
 
                     # Joaquim Luis suggested this as an Exception, so
Index: trunk/anuga_core/source/anuga/geometry/__init__.py
===================================================================
--- trunk/anuga_core/source/anuga/geometry/__init__.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/geometry/__init__.py	(revision 7841)
@@ -2,5 +2,4 @@
 """
 
-pass
 
 #Add path of package to PYTHONPATH to allow C-extensions to be loaded
Index: trunk/anuga_core/source/anuga/geometry/polygon.py
===================================================================
--- trunk/anuga_core/source/anuga/geometry/polygon.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/geometry/polygon.py	(revision 7841)
@@ -5,9 +5,7 @@
 import numpy as num
 
-from math import sqrt
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.geospatial_data.geospatial_data import ensure_absolute, \
                                                     Geospatial_data
-from anuga.config import netcdf_float
 import anuga.utilities.log as log
 
@@ -129,9 +127,9 @@
     line1 = ensure_numeric(line1, num.float)
 
-    x0 = line0[0,0]; y0 = line0[0,1]
-    x1 = line0[1,0]; y1 = line0[1,1]
-
-    x2 = line1[0,0]; y2 = line1[0,1]
-    x3 = line1[1,0]; y3 = line1[1,1]
+    x0 = line0[0, 0]; y0 = line0[0, 1]
+    x1 = line0[1, 0]; y1 = line0[1, 1]
+
+    x2 = line1[0, 0]; y2 = line1[0, 1]
+    x3 = line1[1, 0]; y3 = line1[1, 1]
 
     denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
@@ -208,8 +206,8 @@
     line1 = ensure_numeric(line1, num.float)
 
-    status, value = _intersection(line0[0,0], line0[0,1],
-                                  line0[1,0], line0[1,1],
-                                  line1[0,0], line1[0,1],
-                                  line1[1,0], line1[1,1])
+    status, value = _intersection(line0[0, 0], line0[0, 1],
+                                  line0[1, 0], line0[1, 1],
+                                  line1[0, 0], line1[0, 1],
+                                  line1[1, 0], line1[1, 1])
 
     return status, value
@@ -219,6 +217,5 @@
                        rtol=1.0e-12,
                        atol=1.0e-12,                      
-                       check_inputs=True, 
-                       verbose=False):
+                       check_inputs=True):
     """Determine if one point is inside a triangle
     
@@ -278,8 +275,12 @@
     
     # Quickly reject points that are clearly outside
-    if point[0] < min(triangle[:, 0]): return False 
-    if point[0] > max(triangle[:, 0]): return False    
-    if point[1] < min(triangle[:, 1]): return False
-    if point[1] > max(triangle[:, 1]): return False        
+    if point[0] < min(triangle[:, 0]):
+        return False 
+    if point[0] > max(triangle[:, 0]): 
+        return False    
+    if point[1] < min(triangle[:, 1]):
+        return False
+    if point[1] > max(triangle[:, 1]):
+        return False        
 
 
@@ -314,5 +315,5 @@
         # Check if point lies on one of the edges
         
-        for X, Y in [[A,B], [B,C], [C,A]]:
+        for X, Y in [[A, B], [B, C], [C, A]]:
             res = _point_on_line(point[0], point[1],
                                  X[0], X[1],
@@ -341,4 +342,5 @@
     
     def segments_joined(seg0, seg1):
+        """ See if there are identical segments in the 2 lists. """
         for i in seg0:
             for j in seg1:    
@@ -381,28 +383,4 @@
     return False
     
-    
-def is_inside_polygon_quick(point, polygon, closed=True, verbose=False):
-    """Determine if one point is inside a polygon
-    Both point and polygon are assumed to be numeric arrays or lists and
-    no georeferencing etc or other checks will take place.
-    
-    As such it is faster than is_inside_polygon 
-    """
-
-    # FIXME(Ole): This function isn't being used
-    polygon = ensure_numeric(polygon, num.float)
-    points = ensure_numeric(point, num.float) # Convert point to array of points
-    points = num.ascontiguousarray(points[num.newaxis, :])
-    msg = ('is_inside_polygon() must be invoked with one point only.\n'
-           'I got %s and converted to %s' % (str(point), str(points.shape)))
-    assert points.shape[0] == 1 and points.shape[1] == 2, msg
-    
-    indices = num.zeros(1, num.int)
-
-    count = _separate_points_by_polygon(points, polygon, indices,
-                                        int(closed), int(verbose))
-
-    return count > 0
-
 
 def is_inside_polygon(point, polygon, closed=True, verbose=False):
@@ -420,5 +398,5 @@
     else:
         msg = 'is_inside_polygon must be invoked with one point only'
-        raise msg
+        raise Exception(msg)
 
 ##
@@ -743,6 +721,5 @@
                   figname=None,
                   label=None,
-                  alpha=None,
-                  verbose=False):
+                  alpha=None):
     """ Take list of polygons and plot.
 
@@ -769,5 +746,5 @@
     """
 
-    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, \
+    from pylab import ion, hold, plot, savefig, xlabel, \
                       ylabel, title, close, title, fill
 
@@ -793,8 +770,5 @@
             alpha = None
         else:
-            if alpha < 0.0:
-                alpha = 0.0
-            if alpha > 1.0:
-                alpha = 1.0
+            alpha = max(0.0, min(1.0, alpha))
 
     n = len(polygons_points)
@@ -818,9 +792,13 @@
     for i, item in enumerate(polygons_points):
         x, y = poly_xy(item)
-        if min(x) < minx: minx = min(x)
-        if max(x) > maxx: maxx = max(x)
-        if min(y) < miny: miny = min(y)
-        if max(y) > maxy: maxy = max(y)
-        plot(x,y,colour[i])
+        if min(x) < minx:
+            minx = min(x)
+        if max(x) > maxx:
+            maxx = max(x)
+        if min(y) < miny:
+            miny = min(y)
+        if max(y) > maxy:
+            maxy = max(y)
+        plot(x, y, colour[i])
         if alpha:
             fill(x, y, colour[i], alpha=alpha)
@@ -840,5 +818,5 @@
 
 
-def poly_xy(polygon, verbose=False):
+def poly_xy(polygon):
     """ this is used within plot_polygons so need to duplicate
         the first point so can have closed polygon in plot
@@ -858,8 +836,8 @@
         raise Exception, msg
 
-    x = polygon[:,0]
-    y = polygon[:,1]
-    x = num.concatenate((x, [polygon[0,0]]), axis = 0)
-    y = num.concatenate((y, [polygon[0,1]]), axis = 0)
+    x = polygon[:, 0]
+    y = polygon[:, 1]
+    x = num.concatenate((x, [polygon[0, 0]]), axis = 0)
+    y = num.concatenate((y, [polygon[0, 1]]), axis = 0)
 
     return x, y
@@ -960,7 +938,7 @@
         assert y.shape[0] == N
 
-        points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
-                                                        y[:,num.newaxis]),
-                                                       axis=1 ))
+        points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis],
+                                                        y[:, num.newaxis]),
+                                                       axis = 1 ))
 
         if callable(self.default):
@@ -1042,11 +1020,4 @@
         fid.write('%f, %f\n' % point)
     fid.close()
-
-##
-# @brief Unimplemented.
-def read_tagged_polygons(filename):
-    """
-    """
-    pass
 
 
@@ -1097,9 +1068,9 @@
             if exclude is not None:
                 for ex_poly in exclude:
-                    if is_inside_polygon([x,y], ex_poly):
+                    if is_inside_polygon([x, y], ex_poly):
                         append = False
 
         if append is True:
-            points.append([x,y])
+            points.append([x, y])
 
     return points
@@ -1122,5 +1093,6 @@
     import exceptions
 
-    class Found(exceptions.Exception): pass
+    class Found(exceptions.Exception):
+        pass
 
     polygon = ensure_numeric(polygon)
@@ -1232,8 +1204,8 @@
     # Find outer extent of polygon
     num_polygon = ensure_numeric(polygon)
-    max_x = max(num_polygon[:,0])
-    max_y = max(num_polygon[:,1])
-    min_x = min(num_polygon[:,0])
-    min_y = min(num_polygon[:,1])
+    max_x = max(num_polygon[:, 0])
+    max_y = max(num_polygon[:, 1])
+    min_x = min(num_polygon[:, 0])
+    min_y = min(num_polygon[:, 1])
 
     # Keep only some points making sure extrema are kept
@@ -1257,6 +1229,5 @@
                          interpolation_points=None,
                          rtol=1.0e-6,
-                         atol=1.0e-8,
-                         verbose=False):
+                         atol=1.0e-8):
     """Interpolate linearly between values data on polyline nodes
     of a polyline to list of interpolation points.
@@ -1289,20 +1260,21 @@
     gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
 
-    n = polyline_nodes.shape[0]    # Number of nodes in polyline
+    num_nodes = polyline_nodes.shape[0]    # Number of nodes in polyline
 
     # Input sanity check
-    msg = 'interpolation_points are not given (interpolate.py)'
-    assert interpolation_points is not None, msg
-
-    msg = 'function value must be specified at every interpolation node'
-    assert data.shape[0] == polyline_nodes.shape[0], msg
-
-    msg = 'Must define function value at one or more nodes'
-    assert data.shape[0] > 0, msg
-
-    if n == 1:
-        msg = 'Polyline contained only one point. I need more. ' + str(data)
-        raise Exception, msg
-    elif n > 1:
+    assert_msg = 'interpolation_points are not given (interpolate.py)'
+    assert interpolation_points is not None, assert_msg
+
+    assert_msg = 'function value must be specified at every interpolation node'
+    assert data.shape[0] == polyline_nodes.shape[0], assert_msg
+
+    assert_msg = 'Must define function value at one or more nodes'
+    assert data.shape[0] > 0, assert_msg
+
+    if num_nodes == 1:
+        assert_msg = 'Polyline contained only one point. I need more. '
+        assert_msg += str(data)
+        raise Exception, assert_msg
+    elif num_nodes > 1:
         _interpolate_polyline(data,
                               polyline_nodes,
@@ -1346,8 +1318,8 @@
 
 else:
-    msg = 'C implementations could not be accessed by %s.\n ' %__file__
-    msg += 'Make sure compile_all.py has been run as described in '
-    msg += 'the ANUGA installation guide.'
-    raise Exception, msg
+    error_msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    error_msg += 'Make sure compile_all.py has been run as described in '
+    error_msg += 'the ANUGA installation guide.'
+    raise Exception(error_msg)
 
 
Index: trunk/anuga_core/source/anuga/geometry/test_polygon.py
===================================================================
--- trunk/anuga_core/source/anuga/geometry/test_polygon.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/geometry/test_polygon.py	(revision 7841)
@@ -178,5 +178,5 @@
 
     def test_inside_polygon_main(self):
-        """test_is_inside_polygon_quick
+        """test_is_inside_polygon
         
         Test fast version of of is_inside_polygon
@@ -186,20 +186,20 @@
         polygon = [[0,0], [1,0], [1,1], [0,1]]
 
-        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
-        assert not is_inside_polygon_quick( (-0.5, 0.5), polygon )
-        assert not is_inside_polygon_quick( (1.5, 0.5), polygon )
+        assert is_inside_polygon( (0.5, 0.5), polygon )
+        assert not is_inside_polygon( (0.5, 1.5), polygon )
+        assert not is_inside_polygon( (0.5, -0.5), polygon )
+        assert not is_inside_polygon( (-0.5, 0.5), polygon )
+        assert not is_inside_polygon( (1.5, 0.5), polygon )
 
         # Try point on borders
-        assert is_inside_polygon_quick( (1., 0.5), polygon, closed=True)
-        assert is_inside_polygon_quick( (0.5, 1), polygon, closed=True)
-        assert is_inside_polygon_quick( (0., 0.5), polygon, closed=True)
-        assert is_inside_polygon_quick( (0.5, 0.), polygon, closed=True)
-
-        assert not is_inside_polygon_quick( (0.5, 1), polygon, closed=False)
-        assert not is_inside_polygon_quick( (0., 0.5), polygon, closed=False)
-        assert not is_inside_polygon_quick( (0.5, 0.), polygon, closed=False)
-        assert not is_inside_polygon_quick( (1., 0.5), polygon, closed=False)
+        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
+        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
+        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
+        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
+
+        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
+        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
+        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
+        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
 
 
@@ -229,10 +229,10 @@
         assert not is_inside_polygon( (0.5, -0.5), polygon )
 
-        assert is_inside_polygon_quick( (0.5, 0.5), polygon )
-        assert is_inside_polygon_quick( (1, -0.5), polygon )
-        assert is_inside_polygon_quick( (1.5, 0), polygon )
-
-        assert not is_inside_polygon_quick( (0.5, 1.5), polygon )
-        assert not is_inside_polygon_quick( (0.5, -0.5), polygon )
+        assert is_inside_polygon( (0.5, 0.5), polygon )
+        assert is_inside_polygon( (1, -0.5), polygon )
+        assert is_inside_polygon( (1.5, 0), polygon )
+
+        assert not is_inside_polygon( (0.5, 1.5), polygon )
+        assert not is_inside_polygon( (0.5, -0.5), polygon )
 
         # Very convoluted polygon
@@ -449,5 +449,5 @@
             assert is_inside_polygon(point, polygon)
 
-            assert is_inside_polygon_quick(point, polygon)
+            assert is_inside_polygon(point, polygon)
 
 
@@ -458,5 +458,5 @@
         for point in points:
             assert is_inside_polygon(point, polygon)
-            assert is_inside_polygon_quick(point, polygon)
+            assert is_inside_polygon(point, polygon)
 
 
@@ -1598,6 +1598,6 @@
         assert y[4] == 6
 
-    # Disabled
-    def xtest_plot_polygons(self):
+
+    def test_plot_polygons(self):
         import os
 
@@ -1605,5 +1605,5 @@
         polygon1 = [[0,0], [1,0], [1,1], [0,1]]
         polygon2 = [[1,1], [2,1], [3,2], [2,2]]
-        v = plot_polygons([polygon1, polygon2], 'test1')
+        v = plot_polygons([polygon1, polygon2], figname='test1')
         assert len(v) == 4
         assert v[0] == 0
@@ -1614,5 +1614,5 @@
         # Another case
         polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
-        v = plot_polygons([polygon2,polygon3], 'test2')
+        v = plot_polygons([polygon2,polygon3], figname='test2')
         assert len(v) == 4
         assert v[0] == 1
Index: trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7841)
@@ -1121,367 +1121,4 @@
             self.failUnless(0 ==1,  'Bad input did not throw exception error!')
 
-    def test_export_grid_parallel(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        base_name = 'tegp'
-        #Setup
-        self.domain.set_name(base_name+'_P0_8')
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 1.0)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        self.domain.evolve_to_end(finaltime = 0.0001)
-        #Setup
-        self.domain.set_name(base_name+'_P1_8')
-        swwfile2 = self.domain.get_name() + '.sww'
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        self.domain.evolve_to_end(finaltime = 0.0002)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-        fid.close()
-
-        #Export to ascii/prj files
-        extra_name_out = 'yeah'
-        sww2dem_batch(base_name,
-                    quantities = ['elevation', 'depth'],
-                    extra_name_out = extra_name_out,
-                    cellsize = cellsize,
-                    verbose = self.verbose,
-                    format = 'asc')
-
-        prjfile = base_name + '_P0_8_elevation_yeah.prj'
-        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                #print " -i*cellsize - y",  -i*cellsize - y
-                #print "float(L[i])", float(L[i])
-                assert num.allclose(float(L[i]), -i*cellsize - y)               
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-
-        prjfile = base_name + '_P1_8_elevation_yeah.prj'
-        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                #print " -i*cellsize - y",  -i*cellsize - y
-                #print "float(L[i])", float(L[i])
-                assert num.allclose(float(L[i]), -i*cellsize - y)               
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-        #Check asc file
-        ascfile = base_name + '_P0_8_depth_yeah.asc'
-        prjfile = base_name + '_P0_8_depth_yeah.prj'
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-
-        #Check asc file
-        ascfile = base_name + '_P1_8_depth_yeah.asc'
-        prjfile = base_name + '_P1_8_depth_yeah.prj'
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile2)
-
-
-
-    def DISABLEDtest_sww2domain2(self):
-        ##################################################################
-        #Same as previous test, but this checks how NaNs are handled.
-        ##################################################################
-
-        #FIXME: See ticket 223
-
-        from mesh_factory import rectangular
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(2,2)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.smooth = False
-        domain.store = True
-        domain.set_name('test_file')
-        domain.set_datadir('.')
-        domain.default_order=2
-
-        domain.set_quantity('elevation', lambda x,y: -x/3)
-        domain.set_quantity('friction', 0.1)
-
-        from math import sin, pi
-        Br = Reflective_boundary(domain)
-        Bt = Transmissive_boundary(domain)
-        Bd = Dirichlet_boundary([0.2,0.,0.])
-        Bw = Time_boundary(domain=domain,
-                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
-
-        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
-
-        h = 0.05
-        elevation = domain.quantities['elevation'].vertex_values
-        domain.set_quantity('stage', elevation + h)
-
-        domain.check_integrity()
-
-        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
-            pass
-            #domain.write_time()
-
-
-        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
-
-        # Fail because NaNs are present
-        #domain2 = sww2domain(filename,
-        #                     boundary,
-        #                     fail_if_NaN=True,
-        #                     verbose=self.verbose)        
-        try:
-            domain2 = load_sww_as_domain(filename,
-                                 boundary,
-                                 fail_if_NaN=True,
-                                 verbose=self.verbose)
-        except DataDomainError:
-            # Now import it, filling NaNs to be -9999
-            filler = -9999
-            domain2 = load_sww_as_domain(filename,
-                                 None,
-                                 fail_if_NaN=False,
-                                 NaN_filler=filler,
-                                 verbose=self.verbose)
-        else:
-            raise Exception, 'should have failed'
-
-            
-        # Now import it, filling NaNs to be 0
-        filler = -9999
-        domain2 = load_sww_as_domain(filename,
-                             None,
-                             fail_if_NaN=False,
-                             NaN_filler=filler,
-                             verbose=self.verbose)            
-                             
-        import sys; sys.exit() 
-            
-        # Clean up
-        os.remove(filename)
-
-
-        bits = ['geo_reference.get_xllcorner()',
-                'geo_reference.get_yllcorner()',
-                'vertex_coordinates']
-
-        for quantity in domain.quantities_to_be_stored:
-            bits.append('get_quantity("%s").get_integral()' %quantity)
-            bits.append('get_quantity("%s").get_values()' %quantity)
-
-        for bit in bits:
-        #    print 'testing that domain.'+bit+' has been restored'
-            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
-
-        print 
-        print
-        print domain2.get_quantity('xmomentum').get_values()
-        print
-        print domain2.get_quantity('stage').get_values()
-        print
-             
-        print 'filler', filler
-        print max(domain2.get_quantity('xmomentum').get_values().flat)
-        
-        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
-        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
-        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
-        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
-
-
-
-    #FIXME This fails - smooth makes the comparism too hard for allclose
-    def ztest_sww2domain3(self):
-        ################################################
-        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
-        ################################################
-        from mesh_factory import rectangular
-        #Create basic mesh
-
-        yiel=0.01
-        points, vertices, boundary = rectangular(10,10)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.geo_reference = Geo_reference(56,11,11)
-        domain.smooth = True
-        domain.store = True
-        domain.set_name('bedslope')
-        domain.default_order=2
-        #Bed-slope and friction
-        domain.set_quantity('elevation', lambda x,y: -x/3)
-        domain.set_quantity('friction', 0.1)
-        # Boundary conditions
-        from math import sin, pi
-        Br = Reflective_boundary(domain)
-        Bt = Transmissive_boundary(domain)
-        Bd = Dirichlet_boundary([0.2,0.,0.])
-        Bw = Time_boundary(domain=domain,
-                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
-
-        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
-
-        domain.quantities_to_be_stored['xmomentum'] = 2
-        domain.quantities_to_be_stored['ymomentum'] = 2        
-        #Initial condition
-        h = 0.05
-        elevation = domain.quantities['elevation'].vertex_values
-        domain.set_quantity('stage', elevation + h)
-
-
-        domain.check_integrity()
-        #Evolution
-        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
-        #    domain.write_time()
-            pass
-
-
-        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
-        domain2 = load_sww_as_domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
-        #points, vertices, boundary = rectangular(15,15)
-        #domain2.boundary = boundary
-        ###################
-        ##NOW TEST IT!!!
-        ###################
-
-        os.remove(domain.get_name() + '.sww')
-
-        #FIXME smooth domain so that they can be compared
-
-
-        bits = []
-        for quantity in domain.quantities_to_be_stored:
-            bits.append('quantities["%s"].get_integral()'%quantity)
-
-
-        for bit in bits:
-            #print 'testing that domain.'+bit+' has been restored'
-            #print bit
-            #print 'done'
-            #print ('domain.'+bit), eval('domain.'+bit)
-            #print ('domain2.'+bit), eval('domain2.'+bit)
-            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
-            pass
-
-        ######################################
-        #Now evolve them both, just to be sure
-        ######################################x
-        domain.time = 0.
-        from time import sleep
-
-        final = .5
-        domain.set_quantity('friction', 0.1)
-        domain.store = False
-        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep = yiel, finaltime = final):
-            #domain.write_time()
-            pass
-
-        domain2.smooth = True
-        domain2.store = False
-        domain2.default_order=2
-        domain2.set_quantity('friction', 0.1)
-        #Bed-slope and friction
-        # Boundary conditions
-        Bd2=Dirichlet_boundary([0.2,0.,0.])
-        Br2 = Reflective_boundary(domain2)
-        domain2.boundary = domain.boundary
-        #print 'domain2.boundary'
-        #print domain2.boundary
-        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
-        #domain2.boundary = domain.boundary
-        #domain2.set_boundary({'exterior': Bd})
-
-        domain2.check_integrity()
-
-        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
-            #domain2.write_time()
-            pass
-
-        ###################
-        ##NOW TEST IT!!!
-        ##################
-
-        print '><><><><>>'
-        bits = [ 'vertex_coordinates']
-
-        for quantity in ['elevation','xmomentum','ymomentum']:
-            #bits.append('quantities["%s"].get_integral()'%quantity)
-            bits.append('get_quantity("%s").get_values()' %quantity)
-
-        for bit in bits:
-            print bit
-            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
 
     def test_file_boundary_stsIV_sinewave_ordering(self):
Index: trunk/anuga_core/source/anuga/shallow_water/test_loadsave.py
===================================================================
--- trunk/anuga_core/source/anuga/shallow_water/test_loadsave.py	(revision 7841)
+++ trunk/anuga_core/source/anuga/shallow_water/test_loadsave.py	(revision 7841)
@@ -0,0 +1,207 @@
+#!/usr/bin/env python
+
+import unittest, os, time
+import os.path
+from math import pi, sqrt
+import tempfile
+
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.file.sww import extent_sww
+
+from anuga.config import g, epsilon
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.utilities.numerical_tools import mean, ensure_numeric
+from anuga.geometry.polygon import is_inside_polygon
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.geospatial_data.geospatial_data import Geospatial_data
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross, \
+                                            rectangular
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+from anuga.shallow_water.forcing import Inflow, Cross_section
+from anuga.geospatial_data.geospatial_data import ensure_geospatial
+
+from anuga.utilities.system_tools import get_pathname_from_package
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
+        import Dirichlet_boundary
+from anuga.shallow_water.forcing import Rainfall, Wind_stress
+from anuga.shallow_water.forcing import Inflow, Cross_section
+from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
+
+from shallow_water_domain import Domain
+
+# boundary functions
+from boundaries import Reflective_boundary, \
+            Field_boundary, Transmissive_momentum_set_stage_boundary, \
+            Transmissive_stage_zero_momentum_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary, Dirichlet_boundary, \
+            Time_boundary, File_boundary, AWI_boundary
+
+import numpy as num
+from anuga.config import g
+
+# Get gateway to C implementation of flux function for direct testing
+from shallow_water_ext import flux_function_central as flux_function
+from shallow_water_ext import rotate
+
+
+def set_bottom_friction(tag, elements, domain):
+    if tag == "bottom":
+        domain.set_quantity('friction', 0.09, indices = elements)
+
+def set_top_friction(tag, elements, domain):
+    if tag == "top":
+        domain.set_quantity('friction', 1., indices = elements)
+
+
+def set_all_friction(tag, elements, domain):
+    if tag == 'all':
+        new_values = domain.get_quantity('friction').get_values(indices = elements) + 10.0
+
+        domain.set_quantity('friction', new_values, indices = elements)
+
+
+# For test_fitting_using_shallow_water_domain example
+def linear_function(point):
+    point = num.array(point)
+    return point[:,0]+point[:,1]
+
+
+
+def scalar_func(t, x, y):
+    """Function that returns a scalar.
+
+    Used to test error message when numeric array is expected
+    """
+
+    return 17.7
+
+
+
+class Test_LoadSave(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+        
+    def test_get_flow_through_cross_section_with_geo(self):
+        """test_get_flow_through_cross_section(self):
+
+        Test that the total flow through a cross section can be
+        correctly obtained at run-time from the ANUGA domain.
+
+        This test creates a flat bed with a known flow through it and tests
+        that the function correctly returns the expected flow.
+
+        The specifics are
+        e = -1 m
+        u = 2 m/s
+        h = 2 m
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
+
+        This run tries it with georeferencing and with elevation = -1
+        """
+
+        # Create basic mesh (20m x 3m)
+        width = 3
+        length = 20
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width, length, width)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference=Geo_reference(56, 308500, 6189000))
+
+        domain.default_order = 2
+        domain.set_quantities_to_be_stored(None)
+
+        e = -1.0
+        w = 1.0
+        h = w-e
+        u = 2.0
+        uh = u*h
+
+        Br = Reflective_boundary(domain)     # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
+
+
+        # Initial conditions
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
+        domain.set_quantity('xmomentum', uh)
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        # Interpolation points down the middle
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        interpolation_points = domain.geo_reference.get_absolute(I)
+
+        for t in domain.evolve(yieldstep=0.1, finaltime=0.5):
+            # Shortcuts to quantites
+            stage = domain.get_quantity('stage')
+            xmomentum = domain.get_quantity('xmomentum')
+            ymomentum = domain.get_quantity('ymomentum')
+
+            # Check that quantities are they should be in the interior
+            w_t = stage.get_values(interpolation_points)
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)
+
+            assert num.allclose(w_t, w)
+            assert num.allclose(uh_t, uh)
+            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
+
+            # Check flows through the middle
+            for i in range(5):
+                x = length/2. + i*0.23674563    # Arbitrary
+                cross_section = [[x, 0], [x, width]]
+
+                cross_section = domain.geo_reference.get_absolute(cross_section)
+                Q = domain.get_flow_through_cross_section(cross_section,
+                                                          verbose=False)
+
+                assert num.allclose(Q, uh*width)
+
+        import cPickle        
+        cPickle.dump(domain, open('domain_pickle.txt', 'w'))
+        domain_restored = cPickle.load(open('domain_pickle.txt'))
+
+        for t in domain_restored.evolve(yieldstep=0.1, finaltime=1.0):
+            # Shortcuts to quantites
+            stage = domain_restored.get_quantity('stage')
+            xmomentum = domain_restored.get_quantity('xmomentum')
+            ymomentum = domain_restored.get_quantity('ymomentum')
+            
+            # Check that quantities are they should be in the interior
+            w_t = stage.get_values(interpolation_points)
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)
+
+            assert num.allclose(w_t, w)
+            assert num.allclose(uh_t, uh)
+            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
+
+            # Check flows through the middle
+            for i in range(5):
+                x = length/2. + i*0.23674563    # Arbitrary
+                cross_section = [[x, 0], [x, width]]
+
+                cross_section = domain_restored.geo_reference.get_absolute(cross_section)
+                Q = domain_restored.get_flow_through_cross_section(cross_section,
+                                                          verbose=False)
+
+                assert num.allclose(Q, uh*width)
+                
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_LoadSave, 'test')
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)    
Index: trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7841)
@@ -1339,4 +1339,5 @@
 
                 assert num.allclose(Q, uh*width)
+
 
     def test_get_energy_through_cross_section_with_geo(self):
Index: trunk/anuga_core/source/anuga/shallow_water/tsh2sww.py
===================================================================
--- trunk/anuga_core/source/anuga/shallow_water/tsh2sww.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/shallow_water/tsh2sww.py	(revision 7841)
@@ -13,5 +13,5 @@
 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
 import time, os 
-from anuga.pyvolution.sww_file import SWW_file
+from anuga.file.sww import SWW_file
 from anuga.utilities.numerical_tools import mean
 import anuga.utilities.log as log
Index: trunk/anuga_core/source/anuga/utilities/file_utils.py
===================================================================
--- trunk/anuga_core/source/anuga/utilities/file_utils.py	(revision 7838)
+++ trunk/anuga_core/source/anuga/utilities/file_utils.py	(revision 7841)
@@ -8,4 +8,5 @@
 import numpy as num
 import shutil
+import log
 
 def make_filename(s):
@@ -333,5 +334,5 @@
         if verbose:
             log.critical('Make directory %s' % dir_name)
-        mkdir(dir_name, 0777)
+        os.mkdir(dir_name, 0777)
 
     if verbose:
