source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py @ 5442

#!/usr/bin/env python


import unittest
from Numeric import zeros, array, allclose, Float
from math import sqrt, pi
import tempfile, os
from os import access, F_OK,sep, removedirs,remove,mkdir,getcwd

from anuga.abstract_2d_finite_volumes.util import *
from anuga.config import epsilon
from anuga.shallow_water.data_manager import timefile2netcdf,del_dir

from anuga.utilities.numerical_tools import NAN

from sys import platform 

from anuga.pmesh.mesh import Mesh
from anuga.shallow_water import Domain, Transmissive_boundary
from anuga.shallow_water.data_manager import get_dataobject
from csv import reader,writer
import time
import string

def test_function(x, y):
    return x+y

class Test_Util(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass




    #Geometric
    #def test_distance(self):
    #    from anuga.abstract_2d_finite_volumes.util import distance#
    #
    #    self.failUnless( distance([4,2],[7,6]) == 5.0,
    #                     'Distance is wrong!')
    #    self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
    #                    'distance is wrong!')
    #    self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
    #                    'distance is wrong!')
    #
    #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
    #                    'distance is wrong!')


    def test_file_function_time1(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        """

        #Write file
        import os, time
        from anuga.config import time_format
        from math import sin, pi

        #Typical ASCII file
        finaltime = 1200
        filename = 'test_file_function'
        fid = open(filename + '.txt', 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt

        fid.close()

        #Convert ASCII file to NetCDF (Which is what we really like!)
        timefile2netcdf(filename)


        #Create file function from time series
        F = file_function(filename + '.tms',
                          quantities = ['Attribute0',
                                        'Attribute1',
                                        'Attribute2'])
        
        #Now try interpolation
        for i in range(20):
            t = i*10
            q = F(t)

            #Exact linear intpolation
            assert allclose(q[0], 2*t)
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        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] )


        t = 100 #Two thirds of the way between between 60 and 120
        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 + '.txt')
        os.remove(filename + '.tms')        


        
    def test_spatio_temporal_file_function_basic(self):
        """Test that spatio temporal file function performs the correct
        interpolations in both time and space
        NetCDF version (x,y,t dependency)        
        """
        import time

        #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

        #Create basic mesh and shallow water domain
        points, vertices, boundary = rectangular(3, 3)
        domain1 = Domain(points, vertices, boundary)

        from anuga.utilities.numerical_tools import 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.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)

        # Boundary conditions
        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()

        finaltime = 8
        #Evolution
        t0 = -1
        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
            #print 'Timesteps: %.16f, %.16f' %(t0, t)
            #if t == t0:
            #    msg = 'Duplicate timestep found: %f, %f' %(t0, t)
            #   raise msg
            t0 = t
             
            #domain1.write_time()


        #Now read data from sww and check
        from Scientific.IO.NetCDF import NetCDFFile
        filename = domain1.get_name() + '.' + domain1.format
        fid = NetCDFFile(filename)

        x = fid.variables['x'][:]
        y = fid.variables['y'][:]
        stage = fid.variables['stage'][:]
        xmomentum = fid.variables['xmomentum'][:]
        ymomentum = fid.variables['ymomentum'][:]
        time = fid.variables['time'][:]

        #Take stage vertex values at last timestep on diagonal
        #Diagonal is identified by vertices: 0, 5, 10, 15

        last_time_index = len(time)-1 #Last last_time_index
        d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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])

        diag = concatenate( (Dx, Dy), axis=1)
        d_midpoints = (diag[1:] + diag[:-1])/2

        #Let us see if the file function can find the correct
        #values at the midpoints at the last timestep:
        f = file_function(filename, domain1,
                          interpolation_points = d_midpoints)

        T = f.get_time()
        msg = 'duplicate timesteps: %.16f and %.16f' %(T[-1], T[-2])
        assert not T[-1] == T[-2], msg
        t = time[last_time_index]
        q = f(t, point_id=0); assert allclose(r0, q)
        q = f(t, point_id=1); assert allclose(r1, q)
        q = f(t, point_id=2); assert allclose(r2, q)


        ##################
        #Now do the same for the first timestep

        timestep = 0 #First timestep
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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
        q = f(0, point_id=0); assert allclose(r0, q)
        q = f(0, point_id=1); assert allclose(r1, q)
        q = f(0, point_id=2); assert allclose(r2, q)


        ##################
        #Now do it again for a timestep in the middle

        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))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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=2); assert allclose(r2, q)


        ##################
        #Now check temporal interpolation
        #Halfway between timestep 15 and 16

        timestep = 15
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_0 = (D[0] + D[1])/2
        r1_0 = (D[1] + D[2])/2
        r2_0 = (D[2] + D[3])/2

        #
        timestep = 16
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_1 = (D[0] + D[1])/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

        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=2); assert allclose(r2, q)

        ##################
        #Finally check interpolation 2 thirds of the way
        #between timestep 15 and 16

        # The reference values are
        r0 = (r0_0 + 2*r0_1)/3
        r1 = (r1_0 + 2*r1_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=2); assert allclose(r2, q)

        fid.close()
        import os
        os.remove(filename)



    def test_spatio_temporal_file_function_different_origin(self):
        """Test that spatio temporal file function performs the correct
        interpolations in both time and space where space is offset by
        xllcorner and yllcorner
        NetCDF version (x,y,t dependency)        
        """
        import time

        #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 anuga.coordinate_transforms.geo_reference import Geo_reference
        xllcorner = 2048
        yllcorner = 11000
        zone = 2

        #Create basic mesh and shallow water domain
        points, vertices, boundary = rectangular(3, 3)
        domain1 = Domain(points, vertices, boundary,
                         geo_reference = Geo_reference(xllcorner = xllcorner,
                                                       yllcorner = yllcorner))
        

        from anuga.utilities.numerical_tools import 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.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)

        # Boundary conditions
        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()

        finaltime = 8
        #Evolution
        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
            pass
            #domain1.write_time()


        #Now read data from sww and check
        from Scientific.IO.NetCDF import NetCDFFile
        filename = domain1.get_name() + '.' + domain1.format
        fid = NetCDFFile(filename)

        x = fid.variables['x'][:]
        y = fid.variables['y'][:]
        stage = fid.variables['stage'][:]
        xmomentum = fid.variables['xmomentum'][:]
        ymomentum = fid.variables['ymomentum'][:]
        time = fid.variables['time'][:]

        #Take stage vertex values at last timestep on diagonal
        #Diagonal is identified by vertices: 0, 5, 10, 15

        last_time_index = len(time)-1 #Last last_time_index     
        d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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])

        diag = concatenate( (Dx, Dy), axis=1)
        d_midpoints = (diag[1:] + diag[:-1])/2


        #Adjust for georef - make interpolation points absolute
        d_midpoints[:,0] += xllcorner
        d_midpoints[:,1] += yllcorner                

        #Let us see if the file function can find the correct
        #values at the midpoints at the last timestep:
        f = file_function(filename, domain1,
                          interpolation_points = d_midpoints)

        t = time[last_time_index]                         
        q = f(t, point_id=0); assert allclose(r0, q)
        q = f(t, point_id=1); assert allclose(r1, q)
        q = f(t, point_id=2); assert allclose(r2, q)


        ##################
        #Now do the same for the first timestep

        timestep = 0 #First timestep
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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
        q = f(0, point_id=0); assert allclose(r0, q)
        q = f(0, point_id=1); assert allclose(r1, q)
        q = f(0, point_id=2); assert allclose(r2, q)


        ##################
        #Now do it again for a timestep in the middle

        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))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/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=2); assert allclose(r2, q)


        ##################
        #Now check temporal interpolation
        #Halfway between timestep 15 and 16

        timestep = 15
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_0 = (D[0] + D[1])/2
        r1_0 = (D[1] + D[2])/2
        r2_0 = (D[2] + D[3])/2

        #
        timestep = 16
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_1 = (D[0] + D[1])/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

        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=2); assert allclose(r2, q)

        ##################
        #Finally check interpolation 2 thirds of the way
        #between timestep 15 and 16

        # The reference values are
        r0 = (r0_0 + 2*r0_1)/3
        r1 = (r1_0 + 2*r1_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=2); assert allclose(r2, q)

        fid.close()
        import os
        os.remove(filename)

        


    def qtest_spatio_temporal_file_function_time(self):
        """Test that File function interpolates correctly
        between given times.
        NetCDF version (x,y,t dependency)
        """

        #Create NetCDF (sww) file to be read
        # x: 0, 5, 10, 15
        # y: -20, -10, 0, 10
        # t: 0, 60, 120, ...., 1200
        #
        # test quantities (arbitrary but non-trivial expressions):
        #
        #   stage     = 3*x - y**2 + 2*t
        #   xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2
        #   ymomentum = x**2 + y**2 * sin(t*pi/600)

        #NOTE: Nice test that may render some of the others redundant.

        import os, time
        from anuga.config import time_format
        from Numeric import sin, pi, exp
        from mesh_factory import rectangular
        from shallow_water import Domain
        import anuga.shallow_water.data_manager

        finaltime = 1200
        filename = 'test_file_function'

        #Create a domain to hold test grid
        #(0:15, -20:10)
        points, vertices, boundary =\
                rectangular(4, 4, 15, 30, origin = (0, -20))
        print "points", points

        #print 'Number of elements', len(vertices)
        domain = Domain(points, vertices, boundary)
        domain.smooth = False
        domain.default_order = 2
        domain.set_datadir('.')
        domain.set_name(filename)
        domain.store = True
        domain.format = 'sww'   #Native netcdf visualisation format

        #print points
        start = time.mktime(time.strptime('2000', '%Y'))
        domain.starttime = start


        #Store structure
        domain.initialise_storage()

        #Compute artificial time steps and store
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            #Compute quantities
            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
            domain.set_quantity('stage', f1)

            f2 = lambda x,y: x+y+t**2
            domain.set_quantity('xmomentum', f2)

            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
            domain.set_quantity('ymomentum', f3)

            #Store and advance time
            domain.time = t
            domain.store_timestep(domain.conserved_quantities)
            t += dt


        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
      
        #Deliberately set domain.starttime to too early
        domain.starttime = start - 1

        #Create file function
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)

        #Check that FF updates fixes domain starttime
        assert allclose(domain.starttime, start)

        #Check that domain.starttime isn't updated if later
        domain.starttime = start + 1
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)
        assert allclose(domain.starttime, start+1)
        domain.starttime = start


        #Check linear interpolation in time
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)                
        for id in range(len(interpolation_points)):
            x = interpolation_points[id][0]
            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)

                if q0 == NAN:
                    actual = q0
                else:
                    actual = (k*q1 + (6-k)*q0)/6
                q = F(t, point_id=id)
                #print i, k, t, q
                #print ' ', q0
                #print ' ', q1
                print "q",q
                print "actual", actual
                #print
                if q0 == NAN:
                     self.failUnless( q == actual, 'Fail!')
                else:
                    assert allclose(q, actual)


        #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)

            t = 90 #Halfway between 60 and 120
            q = F(t, point_id=id)
            assert allclose( (q120+q60)/2, q )

            t = 100 #Two thirds of the way between between 60 and 120
            q = F(t, point_id=id)
            assert allclose(q60/3 + 2*q120/3, q)



        #Check that domain.starttime isn't updated if later than file starttime but earlier
        #than file end time
        delta = 23
        domain.starttime = start + delta
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)
        assert allclose(domain.starttime, start+delta)




        #Now try interpolation with delta offset
        for id in range(len(interpolation_points)):            
            x = interpolation_points[id][0]
            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)

                q = F(t-delta, point_id=id)
                assert allclose(q, (k*q1 + (6-k)*q0)/6)


        os.remove(filename + '.sww')



    def NOtest_spatio_temporal_file_function_time(self):
        # FIXME: This passes but needs some TLC
        # Test that File function interpolates correctly
        # When some points are outside the mesh

        import os, time
        from anuga.config import time_format
        from Numeric import sin, pi, exp
        from mesh_factory import rectangular
        from shallow_water import Domain
        import anuga.shallow_water.data_manager 
        from anuga.pmesh.mesh_interface import create_mesh_from_regions
        finaltime = 1200
        
        filename = tempfile.mktemp()
        #print "filename",filename 
        filename = 'test_file_function'

        meshfilename = tempfile.mktemp(".tsh")

        boundary_tags = {'walls':[0,1],'bom':[2]}
        
        polygon_absolute = [[0,-20],[10,-20],[10,15],[-20,15]]
        
        create_mesh_from_regions(polygon_absolute,
                                 boundary_tags,
                                 10000000,
                                 filename=meshfilename)
        domain = Domain(mesh_filename=meshfilename)
        domain.smooth = False
        domain.default_order = 2
        domain.set_datadir('.')
        domain.set_name(filename)
        domain.store = True
        domain.format = 'sww'   #Native netcdf visualisation format

        #print points
        start = time.mktime(time.strptime('2000', '%Y'))
        domain.starttime = start
        

        #Store structure
        domain.initialise_storage()

        #Compute artificial time steps and store
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            #Compute quantities
            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
            domain.set_quantity('stage', f1)

            f2 = lambda x,y: x+y+t**2
            domain.set_quantity('xmomentum', f2)

            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
            domain.set_quantity('ymomentum', f3)

            #Store and advance time
            domain.time = t
            domain.store_timestep(domain.conserved_quantities)
            t += dt

        interpolation_points = [[1,0]]
        interpolation_points = [[100,1000]]
        
        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5],
                                [78787,78787],[7878,3432]]
            
        #Deliberately set domain.starttime to too early
        domain.starttime = start - 1

        #Create file function
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)

        #Check that FF updates fixes domain starttime
        assert allclose(domain.starttime, start)

        #Check that domain.starttime isn't updated if later
        domain.starttime = start + 1
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)
        assert allclose(domain.starttime, start+1)
        domain.starttime = start


        #Check linear interpolation in time
        # checking points inside and outside the mesh
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)
        
        for id in range(len(interpolation_points)):
            x = interpolation_points[id][0]
            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)

                if q0 == NAN:
                    actual = q0
                else:
                    actual = (k*q1 + (6-k)*q0)/6
                q = F(t, point_id=id)
                #print i, k, t, q
                #print ' ', q0
                #print ' ', q1
                #print "q",q
                #print "actual", actual
                #print
                if q0 == NAN:
                     self.failUnless( q == actual, 'Fail!')
                else:
                    assert allclose(q, actual)

        # now lets check points inside the mesh
        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14]] #, [10,-12.5]] - this point doesn't work WHY?
        interpolation_points = [[10,-12.5]]
            
        print "len(interpolation_points)",len(interpolation_points) 
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)

        domain.starttime = start


        #Check linear interpolation in time
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)                
        for id in range(len(interpolation_points)):
            x = interpolation_points[id][0]
            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)

                if q0 == NAN:
                    actual = q0
                else:
                    actual = (k*q1 + (6-k)*q0)/6
                q = F(t, point_id=id)
                print "############"
                print "id, x, y ", id, x, y #k, t, q
                print "t", t
                #print ' ', q0
                #print ' ', q1
                print "q",q
                print "actual", actual
                #print
                if q0 == NAN:
                     self.failUnless( q == actual, 'Fail!')
                else:
                    assert allclose(q, actual)


        #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)

            t = 90 #Halfway between 60 and 120
            q = F(t, point_id=id)
            assert allclose( (q120+q60)/2, q )

            t = 100 #Two thirds of the way between between 60 and 120
            q = F(t, point_id=id)
            assert allclose(q60/3 + 2*q120/3, q)



        #Check that domain.starttime isn't updated if later than file starttime but earlier
        #than file end time
        delta = 23
        domain.starttime = start + delta
        F = file_function(filename + '.sww', domain,
                          quantities = domain.conserved_quantities,
                          interpolation_points = interpolation_points)
        assert allclose(domain.starttime, start+delta)




        #Now try interpolation with delta offset
        for id in range(len(interpolation_points)):            
            x = interpolation_points[id][0]
            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)

                q = F(t-delta, point_id=id)
                assert allclose(q, (k*q1 + (6-k)*q0)/6)


        os.remove(filename + '.sww')

    def test_file_function_time_with_domain(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        Use domain with starttime
        """

        #Write file
        import os, time, calendar
        from anuga.config import time_format
        from math import sin, pi
        from domain import Domain

        finaltime = 1200
        filename = 'test_file_function'
        fid = open(filename + '.txt', 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt

        fid.close()


        #Convert ASCII file to NetCDF (Which is what we really like!)
        timefile2netcdf(filename)



        a = [0.0, 0.0]
        b = [4.0, 0.0]
        c = [0.0, 3.0]

        points = [a, b, c]
        vertices = [[0,1,2]]
        domain = Domain(points, vertices)

        # Check that domain.starttime is updated if non-existing
        F = file_function(filename + '.tms',
                          domain,
                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])  
        assert allclose(domain.starttime, start)

        # Check that domain.starttime is updated if too early
        domain.starttime = start - 1
        F = file_function(filename + '.tms',
                          domain,
                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
        assert allclose(domain.starttime, start)

        # Check that domain.starttime isn't updated if later
        domain.starttime = start + 1
        F = file_function(filename + '.tms',
                          domain,
                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
        assert allclose(domain.starttime, start+1)

        domain.starttime = start
        F = file_function(filename + '.tms',
                          domain,
                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'],
                          use_cache=True)
        

        #print F.precomputed_values
        #print 'F(60)', F(60)
        
        #Now try interpolation
        for i in range(20):
            t = i*10
            q = F(t)

            #Exact linear intpolation
            assert allclose(q[0], 2*t)
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        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] )


        t = 100 #Two thirds of the way between between 60 and 120
        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 + '.tms')
        os.remove(filename + '.txt')        

    def test_file_function_time_with_domain_different_start(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        Use domain with a starttime later than that of file

        ASCII version
        """

        #Write file
        import os, time, calendar
        from anuga.config import time_format
        from math import sin, pi
        from domain import Domain

        finaltime = 1200
        filename = 'test_file_function'
        fid = open(filename + '.txt', 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt

        fid.close()

        #Convert ASCII file to NetCDF (Which is what we really like!)
        timefile2netcdf(filename)        

        a = [0.0, 0.0]
        b = [4.0, 0.0]
        c = [0.0, 3.0]

        points = [a, b, c]
        vertices = [[0,1,2]]
        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
        F = file_function(filename + '.tms', domain,
                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])        
        assert allclose(domain.starttime, start+delta)




        #Now try interpolation with delta offset
        for i in range(20):
            t = i*10
            q = F(t-delta)

            #Exact linear intpolation
            assert allclose(q[0], 2*t)
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        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] )


        t = 100 #Two thirds of the way between between 60 and 120
        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 + '.tms')
        os.remove(filename + '.txt')                



    def test_apply_expression_to_dictionary(self):

        #FIXME: Division is not expected to work for integers.
        #This must be caught.
        foo = array([[1,2,3],
                     [4,5,6]], Float)

        bar = array([[-1,0,5],
                     [6,1,1]], Float)                  

        D = {'X': foo, 'Y': bar}

        Z = apply_expression_to_dictionary('X+Y', D)        
        assert allclose(Z, foo+bar)

        Z = apply_expression_to_dictionary('X*Y', D)        
        assert allclose(Z, foo*bar)        

        Z = apply_expression_to_dictionary('4*X+Y', D)        
        assert allclose(Z, 4*foo+bar)        

        # test zero division is OK
        Z = apply_expression_to_dictionary('X/Y', D)
        assert allclose(1/Z, 1/(foo/bar)) # can't compare inf to inf

        # make an error for zero on zero
        # this is really an error in Numeric, SciPy core can handle it
        # Z = apply_expression_to_dictionary('0/Y', D)

        #Check exceptions
        try:
            #Wrong name
            Z = apply_expression_to_dictionary('4*X+A', D)        
        except NameError:
            pass
        else:
            msg = 'Should have raised a NameError Exception'
            raise msg


        try:
            #Wrong order
            Z = apply_expression_to_dictionary(D, '4*X+A')        
        except AssertionError:
            pass
        else:
            msg = 'Should have raised a AssertionError Exception'
            raise msg        
        

    def test_multiple_replace(self):
        """Hard test that checks a true word-by-word simultaneous replace
        """
        
        D = {'x': 'xi', 'y': 'eta', 'xi':'lam'}
        exp = '3*x+y + xi'
        
        new = multiple_replace(exp, D)
        
        assert new == '3*xi+eta + lam'
                          


    def test_point_on_line_obsolete(self):
        """Test that obsolete call issues appropriate warning"""

        #Turn warning into an exception
        import warnings
        warnings.filterwarnings('error')

        try:
            assert point_on_line( 0, 0.5, 0,1, 0,0 )
        except DeprecationWarning:
            pass
        else:
            msg = 'point_on_line should have issued a DeprecationWarning'
            raise Exception(msg)    

        warnings.resetwarnings()
    
    def test_get_revision_number(self):
        """test_get_revision_number(self):

        Test that revision number can be retrieved.
        """
        if os.environ.has_key('USER') and os.environ['USER'] == 'dgray':
            # I have a known snv incompatability issue,
            # so I'm skipping this test.
            # FIXME when SVN is upgraded on our clusters
            pass
        else:    
            n = get_revision_number()
            assert n>=0


        
    def test_add_directories(self):
        
        import tempfile
        root_dir = tempfile.mkdtemp('_test_util', 'test_util_')
        directories = ['ja','ne','ke']
        kens_dir = add_directories(root_dir, directories)
        assert kens_dir == root_dir + sep + 'ja' + sep + 'ne' + \
               sep + 'ke'
        assert access(root_dir,F_OK)

        add_directories(root_dir, directories)
        assert access(root_dir,F_OK)
        
        #clean up!
        os.rmdir(kens_dir)
        os.rmdir(root_dir + sep + 'ja' + sep + 'ne')
        os.rmdir(root_dir + sep + 'ja')
        os.rmdir(root_dir)

    def test_add_directories_bad(self):
        
        import tempfile
        root_dir = tempfile.mkdtemp('_test_util', 'test_util_')
        directories = ['/\/!@#@#$%^%&*((*:*:','ne','ke']
        
        try:
            kens_dir = add_directories(root_dir, directories)
        except OSError:
            pass
        else:
            msg = 'bad dir name should give OSError'
            raise Exception(msg)    
            
        #clean up!
        os.rmdir(root_dir)

    def test_check_list(self):

        check_list(['stage','xmomentum'])

        
    def test_add_directories(self):
        
        import tempfile
        root_dir = tempfile.mkdtemp('_test_util', 'test_util_')
        directories = ['ja','ne','ke']
        kens_dir = add_directories(root_dir, directories)
        assert kens_dir == root_dir + sep + 'ja' + sep + 'ne' + \
               sep + 'ke'
        assert access(root_dir,F_OK)

        add_directories(root_dir, directories)
        assert access(root_dir,F_OK)
        
        #clean up!
        os.rmdir(kens_dir)
        os.rmdir(root_dir + sep + 'ja' + sep + 'ne')
        os.rmdir(root_dir + sep + 'ja')
        os.rmdir(root_dir)

    def test_add_directories_bad(self):
        
        import tempfile
        root_dir = tempfile.mkdtemp('_test_util', 'test_util_')
        directories = ['/\/!@#@#$%^%&*((*:*:','ne','ke']
        
        try:
            kens_dir = add_directories(root_dir, directories)
        except OSError:
            pass
        else:
            msg = 'bad dir name should give OSError'
            raise Exception(msg)    
            
        #clean up!
        os.rmdir(root_dir)

    def test_check_list(self):

        check_list(['stage','xmomentum'])
        
    def test_remove_lone_verts_d(self):
        verts = [[0,0],[1,0],[0,1]]
        tris = [[0,1,2]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        assert new_verts == verts
        assert new_tris == tris
     

    def test_remove_lone_verts_e(self):
        verts = [[0,0],[1,0],[0,1],[99,99]]
        tris = [[0,1,2]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        assert new_verts == verts[0:3]
        assert new_tris == tris
        
    def test_remove_lone_verts_a(self):
        verts = [[99,99],[0,0],[1,0],[99,99],[0,1],[99,99]]
        tris = [[1,2,4]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        #print "new_verts", new_verts
        assert new_verts == [[0,0],[1,0],[0,1]]
        assert new_tris == [[0,1,2]]
     
    def test_remove_lone_verts_c(self):
        verts = [[0,0],[1,0],[99,99],[0,1]]
        tris = [[0,1,3]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        #print "new_verts", new_verts
        assert new_verts == [[0,0],[1,0],[0,1]]
        assert new_tris == [[0,1,2]]
        
    def test_remove_lone_verts_b(self):
        verts = [[0,0],[1,0],[0,1],[99,99],[99,99],[99,99]]
        tris = [[0,1,2]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        assert new_verts == verts[0:3]
        assert new_tris == tris
     

    def test_remove_lone_verts_e(self):
        verts = [[0,0],[1,0],[0,1],[99,99]]
        tris = [[0,1,2]]
        new_verts, new_tris = remove_lone_verts(verts, tris)
        assert new_verts == verts[0:3]
        assert new_tris == tris
        
    def test_get_min_max_values(self):
        
        list=[8,9,6,1,4]
        min1, max1 = get_min_max_values(list)
        
        assert min1==1 
        assert max1==9
        
    def test_get_min_max_values1(self):
        
        list=[-8,-9,-6,-1,-4]
        min1, max1 = get_min_max_values(list)
        
#        print 'min1,max1',min1,max1
        assert min1==-9 
        assert max1==-1

#    def test_get_min_max_values2(self):
#        '''
#        The min and max supplied are greater than the ones in the 
#        list and therefore are the ones returned
#        '''
#        list=[-8,-9,-6,-1,-4]
#        min1, max1 = get_min_max_values(list,-10,10)
#        
##        print 'min1,max1',min1,max1
#        assert min1==-10 
#        assert max1==10
        
    def test_make_plots_from_csv_files(self):
        
        #if sys.platform == 'win32':  #Windows
            try: 
                import pylab
            except ImportError:
                #ANUGA don't need pylab to work so the system doesn't 
                #rely on pylab being installed 
                return
            
        
            current_dir=getcwd()+sep+'abstract_2d_finite_volumes'
            temp_dir = tempfile.mkdtemp('','figures')
    #        print 'temp_dir',temp_dir
            fileName = temp_dir+sep+'time_series_3.csv'
            file = open(fileName,"w")
            file.write("time,stage,speed,momentum,elevation\n\
1.0, 0, 0, 0, 10 \n\
2.0, 5, 2, 4, 10 \n\
3.0, 3, 3, 5, 10 \n")
            file.close()
    
            fileName1 = temp_dir+sep+'time_series_4.csv'
            file1 = open(fileName1,"w")
            file1.write("time,stage,speed,momentum,elevation\n\
1.0, 0, 0, 0, 5 \n\
2.0, -5, -2, -4, 5 \n\
3.0, -4, -3, -5, 5 \n")
            file1.close()
    
            fileName2 = temp_dir+sep+'time_series_5.csv'
            file2 = open(fileName2,"w")
            file2.write("time,stage,speed,momentum,elevation\n\
1.0, 0, 0, 0, 7 \n\
2.0, 4, -0.45, 57, 7 \n\
3.0, 6, -0.5, 56, 7 \n")
            file2.close()
            
            dir, name=os.path.split(fileName)
            csv2timeseries_graphs(directories_dic={dir:['gauge', 0, 0]},
                                  output_dir=temp_dir,
                                  base_name='time_series_',
                                  plot_numbers=['3-5'],
                                  quantities=['speed','stage','momentum'],
                                  assess_all_csv_files=True,
                                  extra_plot_name='test')
            
            #print dir+sep+name[:-4]+'_stage_test.png'
            assert(access(dir+sep+name[:-4]+'_stage_test.png',F_OK)==True)
            assert(access(dir+sep+name[:-4]+'_speed_test.png',F_OK)==True)
            assert(access(dir+sep+name[:-4]+'_momentum_test.png',F_OK)==True)
    
            dir1, name1=os.path.split(fileName1)
            assert(access(dir+sep+name1[:-4]+'_stage_test.png',F_OK)==True)
            assert(access(dir+sep+name1[:-4]+'_speed_test.png',F_OK)==True)
            assert(access(dir+sep+name1[:-4]+'_momentum_test.png',F_OK)==True)
    
    
            dir2, name2=os.path.split(fileName2)
            assert(access(dir+sep+name2[:-4]+'_stage_test.png',F_OK)==True)
            assert(access(dir+sep+name2[:-4]+'_speed_test.png',F_OK)==True)
            assert(access(dir+sep+name2[:-4]+'_momentum_test.png',F_OK)==True)
    
            del_dir(temp_dir)
        

    def test_sww2csv_gauges(self):

        def elevation_function(x, y):
            return -x
        
        """Most of this test was copied from test_interpolate
        test_interpole_sww2csv
        
        This is testing the gauge_sww2csv function, by creating a sww file and
        then exporting the gauges and checking the results.
        """
        
        # Create mesh
        mesh_file = tempfile.mktemp(".tsh")    
        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
        m = Mesh()
        m.add_vertices(points)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)
        
        # Create shallow water domain
        domain = Domain(mesh_file)
        os.remove(mesh_file)
        
        domain.default_order=2
        
        # This test was made before tight_slope_limiters were introduced
        # Since were are testing interpolation values this is OK
        domain.tight_slope_limiters = 0 
        

        # Set some field values
        domain.set_quantity('elevation', elevation_function)
        domain.set_quantity('friction', 0.03)
        domain.set_quantity('xmomentum', 3.0)
        domain.set_quantity('ymomentum', 4.0)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary( {'exterior': B})

        # This call mangles the stage values.
        domain.distribute_to_vertices_and_edges()
        domain.set_quantity('stage', 1.0)


        domain.set_name('datatest' + str(time.time()))
        domain.format = 'sww'
        domain.smooth = True
        domain.reduction = mean


        sww = get_dataobject(domain)
        sww.store_connectivity()
        sww.store_timestep(['stage', 'xmomentum', 'ymomentum','elevation'])
        domain.set_quantity('stage', 10.0) # This is automatically limited
        # so it will not be less than the elevation
        domain.time = 2.
        sww.store_timestep(['stage','elevation', 'xmomentum', 'ymomentum'])

        # test the function
        points = [[5.0,1.],[0.5,2.]]

        points_file = tempfile.mktemp(".csv")
#        points_file = 'test_point.csv'
        file_id = open(points_file,"w")
        file_id.write("name, easting, northing, elevation \n\
point1, 5.0, 1.0, 3.0\n\
point2, 0.5, 2.0, 9.0\n")
        file_id.close()

        
        sww2csv_gauges(sww.filename, 
                       points_file,
                       verbose=False,
                       use_cache=False)

#        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
        point1_answers_array = [[0.0,1.0,6.0,-5.0,3.0,4.0], [2.0,10.0,15.0,-5.0,3.0,4.0]]
        point1_filename = 'gauge_point1.csv'
        point1_handle = file(point1_filename)
        point1_reader = reader(point1_handle)
        point1_reader.next()

        line=[]
        for i,row in enumerate(point1_reader):
            #print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
            #print 'assert line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point1_answers_array[i])

        point2_answers_array = [[0.0,1.0,1.5,-0.5,3.0,4.0], [2.0,10.0,10.5,-0.5,3.0,4.0]]
        point2_filename = 'gauge_point2.csv' 
        point2_handle = file(point2_filename)
        point2_reader = reader(point2_handle)
        point2_reader.next()
                        
        line=[]
        for i,row in enumerate(point2_reader):
            #print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
            #print 'assert line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point2_answers_array[i])
                         
        # clean up
        point1_handle.close()
        point2_handle.close()
        #print "sww.filename",sww.filename 
        os.remove(sww.filename)
        os.remove(points_file)
        os.remove(point1_filename)
        os.remove(point2_filename)



    def test_sww2csv_gauges1(self):
        from anuga.pmesh.mesh import Mesh
        from anuga.shallow_water import Domain, Transmissive_boundary
        from anuga.shallow_water.data_manager import get_dataobject
        from csv import reader,writer
        import time
        import string

        def elevation_function(x, y):
            return -x
        
        """Most of this test was copied from test_interpolate
        test_interpole_sww2csv
        
        This is testing the gauge_sww2csv function, by creating a sww file and
        then exporting the gauges and checking the results.
        
        This tests the ablity not to have elevation in the points file and 
        not store xmomentum and ymomentum
        """
        
        # Create mesh
        mesh_file = tempfile.mktemp(".tsh")    
        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
        m = Mesh()
        m.add_vertices(points)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)
        
        # Create shallow water domain
        domain = Domain(mesh_file)
        os.remove(mesh_file)
        
        domain.default_order=2

        # Set some field values
        domain.set_quantity('elevation', elevation_function)
        domain.set_quantity('friction', 0.03)
        domain.set_quantity('xmomentum', 3.0)
        domain.set_quantity('ymomentum', 4.0)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary( {'exterior': B})

        # This call mangles the stage values.
        domain.distribute_to_vertices_and_edges()
        domain.set_quantity('stage', 1.0)


        domain.set_name('datatest' + str(time.time()))
        domain.format = 'sww'
        domain.smooth = True
        domain.reduction = mean

        sww = get_dataobject(domain)
        sww.store_connectivity()
        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
        domain.set_quantity('stage', 10.0) # This is automatically limited
        # so it will not be less than the elevation
        domain.time = 2.
        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])

        # test the function
        points = [[5.0,1.],[0.5,2.]]

        points_file = tempfile.mktemp(".csv")
#        points_file = 'test_point.csv'
        file_id = open(points_file,"w")
        file_id.write("name,easting,northing \n\
point1, 5.0, 1.0\n\
point2, 0.5, 2.0\n")
        file_id.close()

        sww2csv_gauges(sww.filename, 
                            points_file,
                            quantities=['stage', 'elevation'],
                            use_cache=False,
                            verbose=False)

        point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
        point1_filename = 'gauge_point1.csv'
        point1_handle = file(point1_filename)
        point1_reader = reader(point1_handle)
        point1_reader.next()

        line=[]
        for i,row in enumerate(point1_reader):
#            print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2])])
            #print 'line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point1_answers_array[i])

        point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
        point2_filename = 'gauge_point2.csv' 
        point2_handle = file(point2_filename)
        point2_reader = reader(point2_handle)
        point2_reader.next()
                        
        line=[]
        for i,row in enumerate(point2_reader):
#            print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2])])
#            print 'line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point2_answers_array[i])
                         
        # clean up
        point1_handle.close()
        point2_handle.close()
        #print "sww.filename",sww.filename 
        os.remove(sww.filename)
        os.remove(points_file)
        os.remove(point1_filename)
        os.remove(point2_filename)


    def test_sww2csv_gauges2(self):

        def elevation_function(x, y):
            return -x
        
        """Most of this test was copied from test_interpolate
        test_interpole_sww2csv
        
        This is testing the gauge_sww2csv function, by creating a sww file and
        then exporting the gauges and checking the results.
        
        This is the same as sww2csv_gauges except set domain.set_starttime to 5.
        Therefore testing the storing of the absolute time in the csv files
        """
        
        # Create mesh
        mesh_file = tempfile.mktemp(".tsh")    
        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
        m = Mesh()
        m.add_vertices(points)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)
        
        # Create shallow water domain
        domain = Domain(mesh_file)
        os.remove(mesh_file)
        
        domain.default_order=2

        # This test was made before tight_slope_limiters were introduced
        # Since were are testing interpolation values this is OK
        domain.tight_slope_limiters = 0         

        # Set some field values
        domain.set_quantity('elevation', elevation_function)
        domain.set_quantity('friction', 0.03)
        domain.set_quantity('xmomentum', 3.0)
        domain.set_quantity('ymomentum', 4.0)
        domain.set_starttime(5)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary( {'exterior': B})

        # This call mangles the stage values.
        domain.distribute_to_vertices_and_edges()
        domain.set_quantity('stage', 1.0)
        


        domain.set_name('datatest' + str(time.time()))
        domain.format = 'sww'
        domain.smooth = True
        domain.reduction = mean

        sww = get_dataobject(domain)
        sww.store_connectivity()
        sww.store_timestep(['stage', 'xmomentum', 'ymomentum','elevation'])
        domain.set_quantity('stage', 10.0) # This is automatically limited
        # so it will not be less than the elevation
        domain.time = 2.
        sww.store_timestep(['stage','elevation', 'xmomentum', 'ymomentum'])

        # test the function
        points = [[5.0,1.],[0.5,2.]]

        points_file = tempfile.mktemp(".csv")
#        points_file = 'test_point.csv'
        file_id = open(points_file,"w")
        file_id.write("name, easting, northing, elevation \n\
point1, 5.0, 1.0, 3.0\n\
point2, 0.5, 2.0, 9.0\n")
        file_id.close()

        
        sww2csv_gauges(sww.filename, 
                            points_file,
                            verbose=False,
                            use_cache=False)

#        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
        point1_answers_array = [[5.0,1.0,6.0,-5.0,3.0,4.0], [7.0,10.0,15.0,-5.0,3.0,4.0]]
        point1_filename = 'gauge_point1.csv'
        point1_handle = file(point1_filename)
        point1_reader = reader(point1_handle)
        point1_reader.next()

        line=[]
        for i,row in enumerate(point1_reader):
            #print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
            #print 'assert line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point1_answers_array[i])

        point2_answers_array = [[5.0,1.0,1.5,-0.5,3.0,4.0], [7.0,10.0,10.5,-0.5,3.0,4.0]]
        point2_filename = 'gauge_point2.csv' 
        point2_handle = file(point2_filename)
        point2_reader = reader(point2_handle)
        point2_reader.next()
                        
        line=[]
        for i,row in enumerate(point2_reader):
            #print 'i',i,'row',row
            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
            #print 'assert line',line[i],'point1',point1_answers_array[i]
            assert allclose(line[i], point2_answers_array[i])
                         
        # clean up
        point1_handle.close()
        point2_handle.close()
        #print "sww.filename",sww.filename 
        os.remove(sww.filename)
        os.remove(points_file)
        os.remove(point1_filename)
        os.remove(point2_filename)


    def test_greens_law(self):

        from math import sqrt
        
        d1 = 80.0
        d2 = 20.0
        h1 = 1.0
        h2 = greens_law(d1,d2,h1)

        assert h2==sqrt(2.0)
        
    def test_calc_bearings(self):
 
        from math import atan, degrees
        #Test East
        uh = 1
        vh = 1.e-15
        angle = calc_bearing(uh, vh)
        if 89 < angle < 91: v=1
        assert v==1
        #Test West
        uh = -1
        vh = 1.e-15
        angle = calc_bearing(uh, vh)
        if 269 < angle < 271: v=1
        assert v==1
        #Test North
        uh = 1.e-15
        vh = 1
        angle = calc_bearing(uh, vh)
        if -1 < angle < 1: v=1
        assert v==1
        #Test South
        uh = 1.e-15
        vh = -1
        angle = calc_bearing(uh, vh)
        if 179 < angle < 181: v=1
        assert v==1
        #Test South-East
        uh = 1
        vh = -1
        angle = calc_bearing(uh, vh)
        if 134 < angle < 136: v=1
        assert v==1
        #Test North-East
        uh = 1
        vh = 1
        angle = calc_bearing(uh, vh)
        if 44 < angle < 46: v=1
        assert v==1
        #Test South-West
        uh = -1
        vh = -1
        angle = calc_bearing(uh, vh)
        if 224 < angle < 226: v=1
        assert v==1
        #Test North-West
        uh = -1
        vh = 1
        angle = calc_bearing(uh, vh)
        if 314 < angle < 316: v=1
        assert v==1

 


       

#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Util,'test')
#    suite = unittest.makeSuite(Test_Util,'test_sww2csv')
#    runner = unittest.TextTestRunner(verbosity=2)
    runner = unittest.TextTestRunner(verbosity=1)
    runner.run(suite)