Index: /inundation/pyvolution-1d/domain.py
===================================================================
--- /inundation/pyvolution-1d/domain.py	(revision 1722)
+++ /inundation/pyvolution-1d/domain.py	(revision 1722)
@@ -0,0 +1,107 @@
+"""Class Domain - 1D domains for finite-volume computations of
+   the shallow water wave equation
+
+
+   Copyright 2004
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia
+"""
+
+class Domain:
+
+    def __init__(self, coordinates, conserved_quantities = None,
+                 other_quantities = None):
+        """
+        Build 1D elements from x coordinates
+        """
+
+        from Numeric import array, zeros, Float, Int
+
+        self.beta = 1.0
+        #Store Points
+        self.coordinates = array(coordinates)
+
+        #Register number of Elements
+        self.number_of_elements = N = len(self.coordinates)-1
+        
+
+        #Allocate space for geometric quantities
+        self.vertices  = zeros((N, 2), Float)
+        self.centroids = zeros(N, Float)
+        self.areas     = zeros(N, Float)
+        for i in range(N):
+            xl = self.coordinates[i]
+            xr = self.coordinates[i+1]
+            self.vertices[i,0] = xl
+            self.vertices[i,1] = xr
+            
+            centroid = (xl+xr)/2
+            self.centroids[i] = centroid
+
+            msg = 'Coordinates should be ordered, smallest to largest'
+            assert xr>xl, msg
+
+            self.areas[i] = (xr-xl)
+
+##         print 'N', N
+##         print 'Centroid', self.centroids
+##         print 'Areas', self.areas
+##         print 'Vertex_Coordinates', self.vertices
+        
+    def get_centroids(self):
+        """Return all coordinates of centroids
+        Return x coordinate of centroid for each element as a N array
+        """
+        
+        return self.centroids
+
+    def get_vertices(self):
+        """Return all coordinates of centroids
+        Return x coordinate of centroid for each element as a N array
+        """
+
+        return self.vertices 
+
+    def get_coordinate(self, elem_id, vertex=None):
+        """Return coordinate of centroid,
+        or left or right vertex. 
+        Left vertex (vertex=0). Right vertex (vertex=1)
+        """
+
+        if vertex is None:
+            return self.centroids[elem_id]
+        else:
+            return self.vertices[elem_id,vertex]
+
+    def get_area(self, elem_id):
+        """Return area of element id
+        """
+
+        return self.areas[elem_id]
+    
+
+if __name__ == "__main__":
+
+    points1 = [0.0, 1.0, 2.0, 3.0]
+    D1 = Domain(points1)
+
+    print D1.get_coordinate(0)
+    print D1.get_coordinate(0,1)
+
+    try:
+        print D1.get_coordinate(3)
+    except:
+        pass
+    else:
+        msg =  'Should have raised an out of bounds exception'
+    raise msg
+    
+    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
+    #D2 = Domain(points2)
+
+    
+
+
+
+
+
Index: /inundation/pyvolution-1d/quantity.py
===================================================================
--- /inundation/pyvolution-1d/quantity.py	(revision 1722)
+++ /inundation/pyvolution-1d/quantity.py	(revision 1722)
@@ -0,0 +1,464 @@
+"""Class Quantity - Implements values at each 1d element
+
+To create:
+
+   Quantity(domain, vertex_values)
+
+   domain: Associated domain structure. Required.
+   
+   vertex_values: N x 2 array of values at each vertex for each element.
+                  Default None
+
+   If vertex_values are None Create array of zeros compatible with domain.
+   Otherwise check that it is compatible with dimenions of domain.
+   Otherwise raise an exception
+"""
+
+
+class Quantity:
+
+    def __init__(self, domain, vertex_values=None):
+
+        from domain import Domain
+        from Numeric import array, zeros, Float
+        
+        msg = 'First argument in Quantity.__init__ '
+        msg += 'must be of class Domain (or a subclass thereof)'
+        assert isinstance(domain, Domain), msg
+
+        if vertex_values is None:
+            N = domain.number_of_elements            
+            self.vertex_values = zeros((N, 2), Float)
+        else:    
+            self.vertex_values = array(vertex_values, Float)
+
+            N, V = self.vertex_values.shape
+            assert V == 2,\
+                   'Two vertex values per element must be specified'
+            
+
+            msg = 'Number of vertex values (%d) must be consistent with'\
+                  %N
+            msg += 'number of elements in specified domain (%d).'\
+                   %domain.number_of_elements
+            
+            assert N == domain.number_of_elements, msg
+
+        self.domain = domain    
+
+        #Allocate space for other quantities
+        self.centroid_values = zeros(N, Float)
+
+        #Intialise centroid values
+        self.interpolate()
+        
+    
+    def interpolate(self):
+        """Compute interpolated values at centroid
+        Pre-condition: vertex_values have been set
+        """
+
+        N = self.vertex_values.shape[0]        
+        for i in range(N):
+            v0 = self.vertex_values[i, 0]
+            v1 = self.vertex_values[i, 1]
+            
+            self.centroid_values[i] = (v0 + v1)/2
+
+    def set_values(self, X, location='vertices'):
+        """Set values for quantity
+
+        X: Compatible list, Numeric array (see below), constant or function
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If values are described a function, it will be evaluated at specified points
+
+        If selected location is vertices, values for centroid and edges
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+        """
+
+        if location not in ['vertices', 'centroids']:
+            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
+            raise msg
+
+        if X is None:
+            msg = 'Given values are None'
+            raise msg            
+        
+        import types
+        
+        if callable(X):
+            #Use function specific method
+            self.set_function_values(X, location)            
+        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
+            if location == 'centroids':
+                self.centroid_values[:] = X
+            else:
+                self.vertex_values[:] = X                
+
+        else:
+            #Use array specific method
+            self.set_array_values(X, location)
+
+        if location == 'vertices':
+            #Intialise centroid
+            self.interpolate()
+            
+
+
+
+    def set_function_values(self, f, location='vertices'):
+        """Set values for quantity using specified function
+
+        f: x -> z Function where x and z are arrays
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid
+                  Default is "vertices"
+        """
+
+        if location == 'centroids':
+            P = self.domain.centroids
+            self.set_values(f(P), location)
+        else:
+            #Vertices
+            P = self.domain.get_vertices()
+            for i in range(2):
+                self.vertex_values[:,i] = f(P[:,i])
+
+
+    def set_array_values(self, values, location='vertices'):
+        """Set values for quantity
+
+        values: Numeric array
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx2
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If selected location is vertices, values for centroid
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+        """
+
+        from Numeric import array, Float
+
+        values = array(values).astype(Float)
+
+        N = self.centroid_values.shape[0]
+        
+        msg = 'Number of values must match number of elements'
+        assert values.shape[0] == N, msg
+        
+        if location == 'centroids':
+            assert len(values.shape) == 1, 'Values array must be 1d'
+            self.centroid_values = values                        
+        else:
+            assert len(values.shape) == 2, 'Values array must be 2d'
+            msg = 'Array must be N x 2'            
+            assert values.shape[1] == 2, msg
+            
+            self.vertex_values = values
+
+
+
+
+class Conserved_quantity(Quantity):
+    """Class conserved quantity adds to Quantity:
+
+    storage and method for updating, and
+    methods for extrapolation from centropid to vertices inluding
+    gradients and limiters
+    """
+
+    def __init__(self, domain, vertex_values=None):
+        Quantity.__init__(self, domain, vertex_values)
+        
+        from Numeric import zeros, Float
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+
+        N = domain.number_of_elements
+        self.explicit_update = zeros(N, Float )
+        self.semi_implicit_update = zeros(N, Float )
+
+        self.gradients = zeros(N, Float)
+        self.qmax = zeros(self.centroid_values.shape, Float)
+        self.qmin = zeros(self.centroid_values.shape, Float)              
+                
+
+    def update(self, timestep):
+        """Update centroid values based on values stored in
+        explicit_update and semi_implicit_update as well as given timestep
+        """
+
+        from Numeric import sum, equal, ones, Float
+        
+        N = self.centroid_values.shape[0]
+        
+        #Explicit updates
+        self.centroid_values += timestep*self.explicit_update
+            
+        #Semi implicit updates
+        denominator = ones(N, Float)-timestep*self.semi_implicit_update
+
+        if sum(equal(denominator, 0.0)) > 0.0:
+            msg = 'Zero division in semi implicit update. Call Stephen :-)'
+            raise msg
+        else:
+            #Update conserved_quantities from semi implicit updates
+            self.centroid_values /= denominator
+
+
+    def compute_gradients(self):
+        """Compute gradients of piecewise linear function defined by centroids of
+        neighbouring volumes.
+        """
+
+        
+        from Numeric import array, zeros, Float
+
+        N = self.centroid_values.shape[0]
+
+
+        G = self.gradients
+        Q = self.centroid_values
+        X = self.domain.centroids
+        
+        for k in range(N):
+
+    	    # first and last elements have boundaries
+
+            if k == 0:
+
+                #Get data
+                k0 = k
+                k1 = k+1
+                
+                q0 = Q[k0]
+                q1 = Q[k1]
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k1] #V1 centroid        
+
+                #Gradient
+                G[k] = (q1 - q0)/(x1 - x0)
+
+            elif k == N-1:
+                
+                #Get data
+                k0 = k
+                k1 = k-1
+                
+                q0 = Q[k0]
+                q1 = Q[k1]
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k1] #V1 centroid        
+
+                #Gradient
+                G[k] = (q1 - q0)/(x1 - x0)
+                
+            else:
+                #Interior Volume (2 neighbours)
+
+                #Get data
+                k0 = k-1
+                k2 = k+1
+
+                q0 = Q[k0]
+                q1 = Q[k]
+                q2 = Q[k2]                    
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k] #V1 centroid (Self)
+                x2 = X[k2] #V2 centroid
+
+                #Gradient
+                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+        
+        return
+        
+    def extrapolate_first_order(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using
+        first order scheme.
+        """
+        
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        for i in range(2):
+            qv[:,i] = qc
+            
+            
+    def extrapolate_second_order(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using
+        second order scheme.
+        """
+        
+        self.compute_gradients()
+
+        G = self.gradients
+        V = self.domain.vertices
+        Qc = self.centroid_values
+        Qv = self.vertex_values
+        
+        #Check each triangle
+        for k in range(self.domain.number_of_elements):
+            #Centroid coordinates            
+            x = self.domain.centroids[k]
+
+            #vertex coordinates
+            x0, x1 = V[k,:]
+
+            #Extrapolate
+            Qv[k,0] = Qc[k] + G[k]*(x0-x)
+            Qv[k,1] = Qc[k] + G[k]*(x1-x)
+
+
+            
+
+    def limit(self):        
+        """Limit slopes for each volume to eliminate artificial variance
+        introduced by e.g. second order extrapolator
+
+        This is an unsophisticated limiter as it does not take into
+        account dependencies among quantities.
+
+        precondition:
+        vertex values are estimated from gradient
+        postcondition:
+        vertex values are updated
+        """
+
+        from Numeric import zeros, Float
+
+        N = self.domain.number_of_elements
+        beta = self.domain.beta
+
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        #Find min and max of this and neighbour's centroid values
+        qmax = self.qmax
+        qmin = self.qmin
+
+        for k in range(N):
+            qmax[k] = qmin[k] = qc[k]
+            
+            for i in [-1,1]:
+                n = k+i
+                if (n >= 0) & (n <= N-1):
+                    qn = qc[n] #Neighbour's centroid value
+
+                    qmin[k] = min(qmin[k], qn)
+                    qmax[k] = max(qmax[k], qn)
+
+        #Phi limiter    
+        for k in range(N):
+
+            #Diffences between centroids and maxima/minima
+            dqmax = qmax[k] - qc[k] 
+            dqmin = qmin[k] - qc[k]
+
+            #Deltas between vertex and centroid values
+            dq = [0.0, 0.0]
+            for i in range(2):
+                dq[i] = qv[k,i] - qc[k]            
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(2):
+                r = 1.0
+                if (dq[i] > 0): r = dqmax/dq[i]
+                if (dq[i] < 0): r = dqmin/dq[i]
+
+                phi = min( min(r*beta, 1), phi )    
+
+            #Then update using phi limiter
+            for i in range(2):
+                qv[k,i] = qc[k] + phi*dq[i]
+
+
+
+if __name__ == "__main__":
+    from domain import Domain
+
+    points1 = [0.0, 1.0, 2.0, 3.0]
+    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+    
+    D1 = Domain(points1)
+
+    Q1 = Conserved_quantity(D1, vertex_values)
+
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
+
+    Q1.set_values(new_vertex_values)
+    
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    new_centroid_values = [20,30,40]
+    Q1.set_values(new_centroid_values,'centroids')
+            
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    def fun(x):
+        return x**2
+
+    Q1.set_values(fun,'vertices')    
+            
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    Xc = Q1.domain.vertices
+    Qc = Q1.vertex_values
+    print Xc
+    print Qc
+
+    Qc[1,0] = 3
+
+    Q1.extrapolate_second_order()
+    Q1.limit()
+    
+    from pylab import *
+    plot(Xc,Qc)
+
+    show()
+
+    
+
+
+
+    
+    
+
+    
+
+
+
+
Index: /inundation/pyvolution-1d/test_domain.py
===================================================================
--- /inundation/pyvolution-1d/test_domain.py	(revision 1722)
+++ /inundation/pyvolution-1d/test_domain.py	(revision 1722)
@@ -0,0 +1,107 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt
+
+from domain import *
+from Numeric import allclose, array
+
+
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+
+
+    def tearDown(self):
+        pass
+
+
+    def test_construct(self):
+
+        a = 0.0
+        b = 2.0
+        c = 2.5
+        d = 3.0
+        e = 10.0
+        f = 12.0
+
+        points = [a, b, c, d, e, f]
+
+        domain = Domain(points) 
+        assert allclose(domain.get_centroids(),[(a+b)*0.5, (b+c)*0.5, (c+d)*0.5, (d+e)*0.5, (e+f)*0.5])
+
+
+    def test_area(self):
+
+        a = 0.0
+        b = 2.0
+        c = 2.5
+        d = 3.0
+        e = 10.0
+        f = 12.0
+
+        points = [a, b, c, d, e, f]
+
+        domain = Domain(points) 
+        assert allclose(domain.get_area(0),b-a)
+        assert allclose(domain.get_area(1),c-b)
+        assert allclose(domain.get_area(2),d-c)
+        assert allclose(domain.get_area(3),e-d)
+        assert allclose(domain.get_area(4),f-e)
+        
+    def test_in_order(self):
+        a = 0.0
+        b = 2.0
+        c = 2.5
+        d = 3.0
+        e = 10.0
+        f = 12.0
+
+        points = [a, c, b, d, e, f]
+
+        try:
+            domain = Domain(points)
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+
+    def test_point(self):
+        a = 0.0
+        b = 2.0
+        c = 2.5
+        d = 3.0
+
+        points = [a, b, c, d]
+
+        domain = Domain(points)
+
+        assert allclose(domain.get_coordinate(0), 0.5*(a+b))
+        assert allclose(domain.get_coordinate(1), 0.5*(b+c))
+        assert allclose(domain.get_coordinate(2), 0.5*(c+d))
+        
+        try:
+            allclose(domain.get_coordinate(3), 0.5*(a+b))
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+        
+        assert allclose(domain.get_coordinate(0,0), a)
+        assert allclose(domain.get_coordinate(1,0), b)
+        assert allclose(domain.get_coordinate(2,0), c)
+
+        assert allclose(domain.get_coordinate(0,1), b)
+        assert allclose(domain.get_coordinate(1,1), c)
+        assert allclose(domain.get_coordinate(2,1), d)
+
+
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /inundation/pyvolution-1d/test_quantity.py
===================================================================
--- /inundation/pyvolution-1d/test_quantity.py	(revision 1722)
+++ /inundation/pyvolution-1d/test_quantity.py	(revision 1722)
@@ -0,0 +1,391 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from quantity import *
+#from config import epsilon
+from Numeric import allclose, array, zeros, Float
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        from domain import Domain
+        
+        a = 0.0
+        b = 1.0
+        c = 2.0
+        d = 2.5
+        e = 3.1
+        f = 4.0
+
+        self.points          = [a, b, c, d, e, f]
+        self.centroids       = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2]
+        self.vertex_values   = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]]
+        self.centroid_values = [[1.5,      2.5,      3.5,      4.75,   5.55]]
+
+        self.domain1 = Domain(self.points[0:2])
+        self.domain5 = Domain(self.points)
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        
+        quantity = Quantity(self.domain5, self.vertex_values)
+        assert allclose(quantity.centroid_values, self.centroid_values)
+
+        try:
+            quantity = Quantity()
+        except:
+            pass
+        else:
+            raise 'Should have raised empty quantity exception'
+
+
+        try:
+            quantity = Quantity([1,2])
+        except AssertionError:
+            pass
+        except:
+            raise 'Should have raised "missing domain object" error'        
+        
+
+    def test_creation_zeros(self):
+       
+        quantity = Quantity(self.domain1)
+        assert allclose(quantity.vertex_values, [[0.,0.]])
+
+
+        quantity = Quantity(self.domain5)
+        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
+                                                 [0.,0.], [0.,0.],
+                                                 [0.,0.]])        
+
+        
+    def test_interpolation(self):
+        quantity = Quantity(self.domain1, [[1,2]])
+        assert allclose(quantity.centroid_values, 1.5) #Centroid
+
+
+    def test_interpolation2(self):
+        quantity = Quantity(self.domain5,
+                            [[1,2], [5,5], [0,9], [-6, 3], [3,4]])
+        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid
+
+    def test_set_values(self):
+        quantity = Quantity(self.domain5)
+
+
+        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
+                            location = 'vertices')
+        assert allclose(quantity.vertex_values,
+                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
+        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
+
+        #Test default
+        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
+        assert allclose(quantity.vertex_values,
+                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])        
+        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
+
+        #Test centroids
+        quantity.set_values([1,2,3,4,5], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
+
+        #Test exceptions
+        try:
+            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
+                                location = 'bas kamel tuba')
+        except:
+            pass
+
+
+        try:
+            quantity.set_values([[1,2], [0,0]])
+        except AssertionError:
+            pass
+        except:
+            raise 'should have raised AssertionError'
+
+    def test_set_values_const(self):
+        quantity = Quantity(self.domain5)
+
+        quantity.set_values(1.0, location = 'vertices')
+        assert allclose(quantity.vertex_values,
+                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
+        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
+
+
+        quantity.set_values(2.0, location = 'centroids')
+        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
+
+
+    def test_set_values_func(self):
+        quantity = Quantity(self.domain5)
+
+        def f(x):
+            return x**2
+
+        quantity.set_values(f, location = 'vertices')
+        assert allclose(quantity.vertex_values,
+                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])        
+        assert allclose(quantity.centroid_values,
+                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
+        
+        quantity.set_values(f, location = 'centroids')
+        assert allclose(quantity.centroid_values,
+                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])        
+
+
+    def test_gradient(self):
+        
+        quantity = Conserved_quantity(self.domain5)
+
+        #Set up for a gradient of (3,0) at mid triangle 
+        quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0],
+                            location = 'centroids')
+        
+        #Gradients
+        quantity.compute_gradients()
+        N = quantity.gradients.shape[0]
+        
+        G = quantity.gradients
+        G0 = zeros(N, Float)
+        Q = quantity.centroid_values
+        X = quantity.domain.centroids
+
+        def grad0(x0,x1,x2,q0,q1,q2):
+            return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+
+        def grad1(x0,x1,q0,q1):
+            return  (q1 - q0)/(x1 - x0)
+
+        #Check Gradients
+        G0[0] = grad1(X[0],X[1],Q[0],Q[1])
+        for i in range(1,4):
+            G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1])
+        G0[4] = grad1(X[3],X[4],Q[3],Q[4])
+
+        assert allclose(G,G0)
+
+        quantity.extrapolate_first_order()
+        
+        assert allclose(quantity.vertex_values, [[2., 2.],
+                                                 [4., 4.],
+                                                 [8., 8.],
+                                                 [2., 2.],
+                                                 [5., 5.]])
+
+        quantity.extrapolate_second_order()
+
+        V = quantity.domain.vertices
+        for k in range(quantity.domain.number_of_elements):
+            G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5
+
+        assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]],
+                                                 [4.-G0[1], 4.+G0[1]],
+                                                 [8.-G0[2], 8.+G0[2]],
+                                                 [2.-G0[3], 2.+G0[3]],
+                                                 [5.-G0[4], 5.+G0[4]]])
+
+
+
+    def test_second_order_extrapolation2(self):
+        quantity = Conserved_quantity(self.domain5)        
+
+        #Set up for a gradient of (2), f(x) = 2x
+        quantity.set_values([1., 3., 4.5, 5.6, 7.1],
+                            location = 'centroids')
+        
+        #Gradients
+        quantity.compute_gradients()
+        G = quantity.gradients
+        
+        allclose(G, 2.0)
+
+        quantity.extrapolate_second_order()
+
+        V = quantity.domain.vertices
+        Q = quantity.vertex_values
+        
+        assert allclose(Q[:,0], 2.0*V[:,0])
+        assert allclose(Q[:,1], 2.0*V[:,1])        
+
+
+
+
+
+        
+
+
+##     def test_first_order_extrapolator(self):
+##         quantity = Conserved_quantity(self.mesh4)                
+
+##         #Test centroids
+##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+##         #Extrapolate
+##         quantity.extrapolate_first_order()
+
+##         #Check vertices but not edge values
+##         assert allclose(quantity.vertex_values,
+##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+
+
+    def test_second_order_limit_extrapolator(self):
+        quantity = Conserved_quantity(self.domain5)                        
+
+        #Set up for a gradient of (3,0) at mid triangle 
+        quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0],
+                            location = 'centroids')
+        
+        quantity.extrapolate_second_order()
+        quantity.limit()
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.vertex_values[0,0] >= 2.0
+        assert quantity.vertex_values[0,0] <= 4.0
+        assert quantity.vertex_values[0,1] >= 2.0
+        assert quantity.vertex_values[0,1] <= 4.0
+        
+        assert quantity.vertex_values[1,0] >= 2.0
+        assert quantity.vertex_values[1,0] <= 8.0
+        assert quantity.vertex_values[1,1] >= 2.0
+        assert quantity.vertex_values[1,1] <= 8.0
+
+        assert quantity.vertex_values[2,0] >= 2.0
+        assert quantity.vertex_values[2,0] <= 8.0
+        assert quantity.vertex_values[2,1] >= 2.0
+        assert quantity.vertex_values[2,1] <= 8.0
+
+        assert quantity.vertex_values[3,0] >= 0.0
+        assert quantity.vertex_values[3,0] <= 8.0
+        assert quantity.vertex_values[3,1] >= 0.0
+        assert quantity.vertex_values[3,1] <= 8.0
+
+        assert quantity.vertex_values[4,0] >= 0.0
+        assert quantity.vertex_values[4,0] <= 2.0
+        assert quantity.vertex_values[4,1] >= 0.0
+        assert quantity.vertex_values[4,1] <= 2.0
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/2)
+        
+        
+        
+
+    def test_limiter(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[3,4], [5,5], [0,0], [-6, 3], [-2,4]],
+                            location = 'vertices')        
+
+        #Limit
+        quantity.limit()
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.vertex_values[0,1] >= quantity.centroid_values[0]
+        assert quantity.vertex_values[1,0] <= quantity.centroid_values[1]
+        
+        assert quantity.vertex_values[1,1] >= quantity.centroid_values[1]
+
+        
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/2)
+        
+
+
+##     def test_distribute_first_order(self):
+##         quantity = Conserved_quantity(self.mesh4)        
+
+##         #Test centroids
+##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+
+##         #Extrapolate
+##         quantity.extrapolate_first_order()
+
+##         #Interpolate
+##         quantity.interpolate_from_vertices_to_edges()        
+
+##         assert allclose(quantity.vertex_values,
+##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+##         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
+##                                                [3,3,3], [4, 4, 4]])
+
+
+
+    def test_update_explicit(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_update_semi_implicit(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_both_updates(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [4.,3.,2.,1.,0.0] )
+        
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] )
+        x /= array( [.9,.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)        
+
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner(verbosity=2)
+    runner.run(suite)
+
+
