Index: anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py	(revision 7855)
@@ -20,6 +20,7 @@
     from numpy import sign, abs, minimum, where
 
-    return where( (sign(a)*sign(b) > 0.0) & (sign(a)*sign(c)>0.0), 
-        sign(a)*minimum(abs(a),abs(b),abs(c)), 0.0 )
+    return where( (sign(a)*sign(b) > 0.0) & (sign(a)*sign(c)>0.0),
+        sign(a)*minimum(minimum(abs(a),abs(b)),abs(c)), 0.0 )
+
 
 
Index: anuga_work/development/2010-projects/anuga_1d/base/quantity.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/base/quantity.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/base/quantity.py	(revision 7855)
@@ -26,8 +26,8 @@
         #Initialise Quantity using optional vertex values.
         
-        from generic_domain import Generic_domain
+        from anuga_1d.base.generic_domain import Generic_domain
 
         msg = 'First argument in Quantity.__init__ '
-        msg += 'must be of class Domain (or a subclass thereof)'
+        msg += 'must be of class Generic_domain (or a subclass thereof)'
         assert isinstance(domain, Generic_domain), msg
 
@@ -78,5 +78,14 @@
         self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
 
-        self.beta = domain.beta        
+        self.beta = domain.beta
+
+
+        self.beta_p = numpy.zeros(N,numpy.float)
+        self.beta_m = numpy.zeros(N,numpy.float)
+        self.beta_x = numpy.zeros(N,numpy.float)
+
+
+        self.dx  = numpy.zeros((N,2), numpy.float)
+        self.phi = numpy.zeros(N, numpy.float)
 
 
@@ -602,19 +611,20 @@
 
 
-
-        N = self.centroid_values.shape[0]
-
         #Explicit updates
         self.centroid_values += timestep*self.explicit_update
         
         #Semi implicit updates
-        denominator = numpy.ones(N, numpy.float)-timestep*self.semi_implicit_update
-
-        if sum(numpy.equal(denominator, 0.0)) > 0.0:
-            msg = 'Zero division in semi implicit update. Call Stephen :-)'
-            raise Exception(msg)
-        else:
-            #Update conserved_quantities from semi implicit updates
-            self.centroid_values /= denominator
+        denominator = 1.0-timestep*self.semi_implicit_update
+
+#        if sum(numpy.equal(denominator, 0.0)) > 0.0:
+#            msg = 'Zero division in semi implicit update. Call Stephen :-)'
+#            raise Exception(msg)
+#        else:
+#            #Update conserved_quantities from semi implicit updates
+#            self.centroid_values /= denominator
+#            
+
+        #Update conserved_quantities from semi implicit updates
+        self.centroid_values /= denominator
 
 
@@ -625,6 +635,4 @@
 
 
-
-
         N = self.centroid_values.shape[0]
 
@@ -634,77 +642,65 @@
         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
-                k2 = k+2
-
-                q0 = Q[k0]
-                q1 = Q[k1]
-                q2 = Q[k2]
-
-                x0 = X[k0] #V0 centroid
-                x1 = X[k1] #V1 centroid
-                x2 = X[k2]
-
-                #Gradient
-                #G[k] = (q1 - q0)/(x1 - x0)
-                
-                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
-                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
-
-            elif k == N-1:
-
-                #Get data
-                k0 = k
-                k1 = k-1
-                k2 = k-2
-
-                q0 = Q[k0]
-                q1 = Q[k1]
-                q2 = Q[k2]
-
-                x0 = X[k0] #V0 centroid
-                x1 = X[k1] #V1 centroid
-                x2 = X[k2]
-
-                #Gradient
-                #G[k] = (q1 - q0)/(x1 - x0)
-                
-                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
-                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
-
-##                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] = (q2-q0)/(x2-x0)
-                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+    	# first element
+
+        k = 0
+
+        #Get data
+        k0 = k
+        k1 = k+1
+        k2 = k+2
+
+        q0 = Q[k0]
+        q1 = Q[k1]
+        q2 = Q[k2]
+
+        x0 = X[k0] #V0 centroid
+        x1 = X[k1] #V1 centroid
+        x2 = X[k2]
+
+        #Gradient
+        #G[k] = (q1 - q0)/(x1 - x0)
+
+        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+        #last element
+        k = N-1
+
+
+        k0 = k
+        k1 = k-1
+        k2 = k-2
+
+        q0 = Q[k0]
+        q1 = Q[k1]
+        q2 = Q[k2]
+
+        x0 = X[k0] #V0 centroid
+        x1 = X[k1] #V1 centroid
+        x2 = X[k2]
+
+        #Gradient
+        #G[k] = (q1 - q0)/(x1 - x0)
+
+        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+
+        #Interior Volume (2 neighbours)
+
+
+        q0 = Q[0:-2]
+        q1 = Q[1:-1]
+        q2 = Q[2:]
+
+        x0 = X[0:-2] #V0 centroid
+        x1 = X[1:-1]  #V1 centroid (Self)
+        x2 = X[2:] #V2 centroid
+
+        #Gradient
+        #G[k] = (q2-q0)/(x2-x0)
+        G[1:-1] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+
 
 
@@ -719,5 +715,6 @@
         
         def xmin(a,b):
-            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
+            from numpy import sign, minimum
+            return 0.5*(sign(a)+sign(b))*minimum(abs(a),abs(b))
 
         def xmic(t,a,b):
@@ -733,89 +730,71 @@
         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
-                k2 = k+2
-
-                q0 = Q[k0]
-                q1 = Q[k1]
-                q2 = Q[k2]
-
-                x0 = X[k0] #V0 centroid
-                x1 = X[k1] #V1 centroid
-                x2 = X[k2]
-
-                #Gradient
-                #G[k] = (q1 - q0)/(x1 - x0)
-                
-                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
-                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
-
-            elif k == N-1:
-
-                #Get data
-                k0 = k
-                k1 = k-1
-                k2 = k-2
-
-                q0 = Q[k0]
-                q1 = Q[k1]
-                q2 = Q[k2]
-
-                x0 = X[k0] #V0 centroid
-                x1 = X[k1] #V1 centroid
-                x2 = X[k2]
-
-                #Gradient
-                #G[k] = (q1 - q0)/(x1 - x0)
-                
-                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
-                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
-
-##                #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 (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
-                G[k] = 0.0
-
-            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
-
-                # assuming uniform grid
-                d1 = (q1 - q0)/(x1-x0)
-                d2 = (q2 - q1)/(x2-x1)
-
-                #Gradient
-                #G[k] = (d1+d2)*0.5
-                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)                
-                G[k] = xmic( self.domain.beta, d1, d2 )
+        #-----------------
+        #first element
+        #-----------------
+        k = 0
+
+        k0 = k
+        k1 = k+1
+        k2 = k+2
+
+        q0 = Q[k0]
+        q1 = Q[k1]
+        q2 = Q[k2]
+
+        x0 = X[k0] #V0 centroid
+        x1 = X[k1] #V1 centroid
+        x2 = X[k2]
+
+        #Gradient
+        #G[k] = (q1 - q0)/(x1 - x0)
+
+        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+        #-------------------
+        # Last element
+        #-------------------
+        k = N-1
+
+        k0 = k
+        k1 = k-1
+        k2 = k-2
+
+        q0 = Q[k0]
+        q1 = Q[k1]
+        q2 = Q[k2]
+
+        x0 = X[k0] #V0 centroid
+        x1 = X[k1] #V1 centroid
+        x2 = X[k2]
+
+        #Gradient
+        #G[k] = (q1 - q0)/(x1 - x0)
+
+        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
+        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
+
+
+
+        #------------------------------
+        #Interior Volume (2 neighbours)
+        #------------------------------
+
+        q0 = Q[0:-2]
+        q1 = Q[1:-1]
+        q2 = Q[2:]
+
+        x0 = X[0:-2] #V0 centroid
+        x1 = X[1:-1]  #V1 centroid (Self)
+        x2 = X[2:] #V2 centroid
+
+        # assuming uniform grid
+        d1 = (q1 - q0)/(x1-x0)
+        d2 = (q2 - q1)/(x2-x1)
+
+        #Gradient
+        G[1:-1] = xmic( self.domain.beta, d1, d2 )      
+
         
 
@@ -855,4 +834,5 @@
         """ Find min and max of this and neighbour's centroid values"""
 
+        from numpy import maximum, minimum
 
         qmax = self.qmax
@@ -860,14 +840,27 @@
 
         qc = self.centroid_values
-
-        for k in range(N):
-            qmax[k] = qmin[k] = qc[k]
-            for i in range(2):
-                n = self.domain.neighbours[k,i]
-                if n >= 0:
-                    qn = qc[n] #Neighbour's centroid value
-
-                    qmin[k] = min(qmin[k], qn)
-                    qmax[k] = max(qmax[k], qn)
+        
+        qmax[:] = qc
+        qmin[:] = qc
+        
+        # check left
+        qmax[1:] = maximum(qmax[1:], qc[0:-1])
+        qmin[1:] = minimum(qmin[1:], qc[0:-1])
+        
+        # check right
+        qmax[0:-1] = maximum(qmax[0:-1], qc[1:])
+        qmin[0:-1] = minimum(qmin[0:-1], qc[1:])
+
+
+
+#        for k in range(N):
+#            qmax[k] = qmin[k] = qc[k]
+#            for i in range(2):
+#                n = self.domain.neighbours[k,i]
+#                if n >= 0:
+#                    qn = qc[n] #Neighbour's centroid value
+#
+#                    qmin[k] = min(qmin[k], qn)
+#                    qmax[k] = max(qmax[k], qn)
 
 
@@ -990,5 +983,5 @@
         import sys
 
-        from limiters_python import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
+        from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
 
         limiter = self.domain.limiter
@@ -1003,23 +996,9 @@
         x1 = self.domain.vertices[:,1]
 
-        beta_p = numpy.zeros(N,numpy.float)
-        beta_m = numpy.zeros(N,numpy.float)
-        beta_x = numpy.zeros(N,numpy.float)
-        
-#        for k in range(N):
-#
-#            n0 = self.domain.neighbours[k,0]
-#            n1 = self.domain.neighbours[k,1]
-#
-#            if ( n0 >= 0) & (n1 >= 0):
-#                #SLOPE DERIVATIVE LIMIT
-#                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
-#                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
-#                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
-
-
-
-        n0 = self.domain.neighbours[:,0]
-        n1 = self.domain.neighbours[:,1]
+        beta_p = self.beta_p
+        beta_m = self.beta_m
+        beta_x = self.beta_x
+        phi = self.phi
+        dx  = self.dx
 
 
@@ -1027,12 +1006,7 @@
         beta_m[:-1] = beta_p[1:]
         beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2])
-                
-        dx = numpy.zeros(qv.shape, numpy.float)
 
         dx[:,0] = x0 - xc
         dx[:,1] = x1 - xc
-
-
-        phi = numpy.zeros(qc.shape, numpy.float)
 
         phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0])
@@ -1053,9 +1027,7 @@
 
         elif limiter == "superbee":
-
             slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1])
             slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1])
             phi[1:-1] = maxmod(slope1,slope2)
-
 
         else:
Index: anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py	(revision 7855)
@@ -3,8 +3,8 @@
 import unittest
 from math import sqrt, pi
-
+import numpy
 
 from anuga_1d.base.generic_domain import *
-import numpy
+
 
 
@@ -47,4 +47,5 @@
 
 
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
Index: anuga_work/development/2010-projects/anuga_1d/base/test_quantity.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/base/test_quantity.py	(revision 7852)
+++ 	(revision )
@@ -1,604 +1,0 @@
-import quantity
-#!/usr/bin/env python
-
-import unittest
-from math import sqrt, pi
-
-
-from generic_domain import Generic_domain as Domain
-#from shallow_water_domain import flux_function as domain_flux_function
-
-from quantity import *
-
-
-
-from numpy import allclose, array, ones, zeros
-import numpy
-
-
-class Test_Quantity(unittest.TestCase):
-    def setUp(self):
-        self.points3        = [0.0, 1.0, 2.0, 3.0]
-        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
-        self.domain3        = Domain(self.points3)
-
-
-        
-        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
-        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
-        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
-        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
-        self.domain4          = Domain(self.points4,self.boundary4)
-
-        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
-        self.domain10 = Domain(self.points10)
-
-        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
-        self.domain6 = Domain(self.points6)
-
-
-    def tearDown(self):
-        pass
-        #print "  Tearing down"
-
-
-    def test_creat_with_boundary(self):
-
-        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
-        
-    def test_creation(self):
-
-        quantity = Quantity(self.domain3)
-        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
-
-        
-        try:
-            quantity = Quantity()
-        except:
-            pass
-        else:
-            raise 'Should have raised empty quantity exception'
-
-
-        try:
-            quantity = Quantity([1,2,3])
-        except AssertionError:
-            pass
-        except:
-            raise 'Should have raised "mising domain object" error'
-
-
-    def test_creation_zeros(self):
-
-        quantity = Quantity(self.domain3)
-        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
-
-
-        quantity = Quantity(self.domain4)
-        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
-                                                 [0.,0.], [0.,0.]])
-
-
-    def test_interpolation(self):
-        quantity = Quantity(self.domain4, self.vertex_values4)
-        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
-
-
-
-    def test_interpolation2(self):
-        quantity = Quantity(self.domain4, self.vertex_values4)
-        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
-
-        quantity.extrapolate_second_order()
-
-        #print quantity.vertex_values
-        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
-                                                [ 4.5,  4.5],
-                                                [ 0.5,  0.5],
-                                                [ 1.3,  7.7]])
-
- 
-
- 
-    def test_boundary_allocation(self):
-        quantity = Quantity(self.domain4,
-                            [[1,2], [5,5], [0,9], [-6, 3]])
-
-
-        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
-
-
-    def test_set_values(self):
-        quantity = Quantity(self.domain4)
-
-        # Test vertices
-        quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices')
-        assert allclose(quantity.vertex_values,
-                        [[1,2], [5,5], [0,9], [-6, 3]])
-        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid
-
-
-
-        # Test unique_vertices
-        quantity.set_values([1,2,4,-5,6], location='unique_vertices')
-        assert allclose(quantity.vertex_values,
-                        [[1,2], [2,4], [4,-5], [-5,6]])
-        assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid
-
-
-        # Test centroids
-        quantity.set_values([1,2,3,4], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
-
-        # Test exceptions
-        try:
-            quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]],
-                                location = 'bas kamel tuba')
-        except:
-            pass
-
-
-        try:
-            quantity.set_values([[1,3], [0,9]])
-        except AssertionError:
-            pass
-        except:
-            raise 'should have raised Assertionerror'
-
-
-
-    def test_set_values_const(self):
-        quantity = Quantity(self.domain4)
-
-        quantity.set_values(1.0, location = 'vertices')
-
-        assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]])
-        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
-
-
-        quantity.set_values(2.0, location = 'centroids')
-
-        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
-
-
-    def test_set_values_func(self):
-        quantity = Quantity(self.domain4)
-
-        def f(x):
-            return x*x
-
-
-
-        quantity.set_values(f, location = 'vertices')
-
-
-        assert allclose(quantity.vertex_values,
-                        [[0,1], [1,6.25], [6.25,9], [9,25]])
-
-        assert allclose(quantity.centroid_values,
-                        [0.5, 3.625, 7.625, 34*0.5])
-
-
-        quantity.set_values(f, location = 'centroids')
-
-        
-        assert allclose(quantity.centroid_values,
-                        [0.25, 3.0625, 7.5625, 16.0])
-
-
-    def test_integral(self):
-        quantity = Quantity(self.domain4)
-
-        #Try constants first
-        const = 5
-        quantity.set_values(const, location = 'vertices')
-        #print 'Q', quantity.get_integral()
-
-        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
-
-        #Try with a linear function
-        def f(x):
-            return x
-
-        quantity.set_values(f, location = 'vertices')
-
- 
-        assert allclose (quantity.centroid_values,
-        [ 0.5,   1.75,  2.75,  4.  ])
-
-        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
-                                                  [ 1.,   2.5],
-                                                  [ 2.5,  3. ],
-                                                  [ 3.,   5. ]])
-
-
-        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
-
-        assert allclose (quantity.get_integral(), ref_integral)
-
-
-
-
-
-    def test_set_values_from_quantity(self):
-
-        quantity1 = Quantity(self.domain4)
-        quantity1.set_values([0,1,2,3,4], location='unique_vertices')
-
-        assert allclose(quantity1.vertex_values,
-                        [[0,1], [1,2], [2,3], [3,4]])
-
-
-        quantity2 = Quantity(self.domain4)
-        quantity2.set_values(quantity1)
-        assert allclose(quantity2.vertex_values,
-                        [[0,1], [1,2], [2,3], [3,4]])
-
-        quantity2.set_values(2*quantity1)
-
-        assert allclose(quantity2.vertex_values,
-                        [[0,2], [2,4], [4,6], [6,8]])
-
-        quantity2.set_values(2*quantity1 + 3)
-        assert allclose(quantity2.vertex_values,
-                        [[3,5], [5,7], [7,9], [9,11]])
-
-
-
-    def test_overloading(self):
-
-        quantity1 = Quantity(self.domain4)
-        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
-                            location = 'vertices')
-
-        assert allclose(quantity1.vertex_values,
-                        [[0,1], [1,2], [2,3], [3,4]])
-
-
-        quantity2 = Quantity(self.domain4)
-        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
-                             location = 'vertices')
-
-
-
-        quantity3 = Quantity(self.domain4)
-        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
-                             location = 'vertices')
-
-
-        # Negation
-        Q = -quantity1
-        assert allclose(Q.vertex_values, -quantity1.vertex_values)
-        assert allclose(Q.centroid_values, -quantity1.centroid_values)
-
-
-        # Addition
-        Q = quantity1 + 7
-        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
-        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-
-        Q = 7 + quantity1
-        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
-        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-
-        Q = quantity1 + quantity2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values + quantity2.vertex_values)
-        assert allclose(Q.centroid_values,
-                        quantity1.centroid_values + quantity2.centroid_values)
-
-        Q = quantity1 + quantity2 - 3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values + quantity2.vertex_values - 3)
-
-        Q = quantity1 - quantity2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values - quantity2.vertex_values)
-
-        #Scaling
-        Q = quantity1*3
-        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
-        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
-
-        Q = 3*quantity1
-        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
-
-        #Multiplication
-        Q = quantity1 * quantity2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values)
-
-        #Linear combinations
-        Q = 4*quantity1 + 2
-        assert allclose(Q.vertex_values,
-                        4*quantity1.vertex_values + 2)
-
-        Q = quantity1*quantity2 + 2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values + 2)
-
-        Q = quantity1*quantity2 + quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        quantity3.vertex_values)
-        Q = quantity1*quantity2 + 3*quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        3*quantity3.vertex_values)
-        Q = quantity1*quantity2 + 3*quantity3 + 5.0
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        3*quantity3.vertex_values + 5)
-
-        Q = quantity1*quantity2 - quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values -
-                        quantity3.vertex_values)
-        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
-        assert allclose(Q.vertex_values,
-                        1.5*quantity1.vertex_values * quantity2.vertex_values -
-                        3*quantity3.vertex_values + 5)
-
-        #Try combining quantities and arrays and scalars
-        Q = 1.5*quantity1*quantity2.vertex_values -\
-            3*quantity3.vertex_values + 5.0
-        assert allclose(Q.vertex_values,
-                        1.5*quantity1.vertex_values * quantity2.vertex_values -
-                        3*quantity3.vertex_values + 5)
-
-
-        #Powers
-        Q = quantity1**2
-        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
-
-        Q = quantity1**2 +quantity2**2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values**2 + \
-                        quantity2.vertex_values**2)
-
-        Q = (quantity1**2 +quantity2**2)**0.5
-        assert allclose(Q.vertex_values,
-                        (quantity1.vertex_values**2 + \
-                        quantity2.vertex_values**2)**0.5)
-
-    def test_compute_gradient(self):
-        quantity = Quantity(self.domain6)
-
-        #Set up for a gradient of (2,0) at mid triangle
-        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
-                            location = 'centroids')
-
-
-        #Gradients
-        quantity.compute_gradients()
-
-        a = quantity.gradients
-
-        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
-        
-        quantity.extrapolate_second_order()
-
-
-        assert allclose(quantity.vertex_values, [[ 1., 3. ],
-                                                 [ 4., 4. ],
-                                                 [ 4., 4. ],
-                                                 [ 4., 6.],
-                                                 [ 8.25, 11.75],
-                                                 [ 11.,  13. ]])
-
-                                                 
-
-    def test_second_order_extrapolation2(self):
-        quantity = Quantity(self.domain4)
-
-        #Set up for a gradient of (3,1), f(x) = 3x+y
-        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
-                            location = 'centroids')
-
-        #Gradients
-        quantity.compute_gradients()
-
-        a = quantity.gradients
-        
-
-        assert allclose(a[1], 2.8)
-
-        #Work out the others
-
-        quantity.extrapolate_second_order()
-        
-
-        assert allclose(quantity.vertex_values[1,0], 3.33333333)
-        assert allclose(quantity.vertex_values[1,1], 7.33333333)
-
-        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
-
-
-
-
-    def test_backup_saxpy_centroid_values(self):
-        quantity = Quantity(self.domain4)
-
-        #Set up for a gradient of (3,1), f(x) = 3x+y
-        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
-        d_values = array([1.0, 2.0, 3.0, 4.0])
-        quantity.set_values(c_values, location = 'centroids')
-
-        #Backup
-        quantity.backup_centroid_values()
-
-        #print quantity.vertex_values
-        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
-
-
-        quantity.set_values(d_values, location = 'centroids')
-
-        quantity.saxpy_centroid_values(2.0, 3.0)
-
-        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
-
-
-
-    def test_first_order_extrapolator(self):
-        quantity = Quantity(self.domain4)
-
-        centroid_values = array([1.,2.,3.,4.])
-        quantity.set_values(centroid_values, location = 'centroids')
-        assert allclose(quantity.centroid_values, centroid_values) #Centroid
-
-        #Extrapolate
-        quantity.extrapolate_first_order()
-
-        #Check that gradient is zero
-        a = quantity.gradients
-        assert allclose(a, [0,0,0,0])
-       
-
-        #Check vertices but not edge values
-        assert allclose(quantity.vertex_values,
-                        [[1,1], [2,2], [3,3], [4,4]])
-
-
-    def test_second_order_extrapolator(self):
-        quantity = Quantity(self.domain4)
-
-        #Set up for a gradient of (3,0) at mid triangle
-        quantity.set_values([2.0, 4.0, 8.0, 2.0],
-                            location = 'centroids')
-
-
-
-        quantity.extrapolate_second_order()
-
-
-        #Assert that quantities are conserved
-        for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             numpy.sum(quantity.vertex_values[k,:])/2.0)
-
-
-    def test_limit(self):
-        quantity = Quantity(self.domain4)
-
-        #Create a deliberate overshoot (e.g. from gradient computation)
-        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
-
-        #Limit
-        quantity.limit_minmod()
-
-       
-        #cells 1 and 2 are local max and min
-        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
-        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
-
-        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
-        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
-
-
-
-    def test_distribute_first_order(self):
-        quantity = Quantity(self.domain4)
-
-        #Test centroids
-        centroid_values = array([1.,2.,3.,4.])
-        quantity.set_values(centroid_values, location = 'centroids')
-        assert allclose(quantity.centroid_values, centroid_values) #Centroid
-
-
-        #Extrapolate from centroid to vertices and edges
-        quantity.extrapolate_first_order()
-        
-        assert allclose(quantity.vertex_values,[[ 1.,  1.],
-                                                [ 2.,  2.],
-                                                [ 3.,  3.],
-                                                [ 4.,  4.]])
- 
-
-
-    def test_distribute_second_order(self):
-        quantity = Quantity(self.domain4)
-
-        #Test centroids
-        centroid_values = array([2.,4.,8.,2.])
-        quantity.set_values(centroid_values, location = 'centroids')
-        assert allclose(quantity.centroid_values, centroid_values) #Centroid
-
-
-        #Extrapolate
-        quantity.extrapolate_second_order()
-
-        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
-                                                 [ 2.,  6. ],
-                                                 [ 8.,   8. ],
-                                                 [ 6.8, -2.8]])
-
-
-    def test_update_explicit(self):
-        quantity = Quantity(self.domain4)
-
-        #Test centroids
-        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-        #Set explicit_update
-        quantity.explicit_update = array( [1.,1.,1.,1.] )
-
-        #Update with given timestep
-        quantity.update(0.1)
-
-        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
-        assert allclose( quantity.centroid_values, x)
-
-    def test_update_semi_implicit(self):
-        quantity = Quantity(self.domain4)
-
-        #Test centroids
-        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-        #Set semi implicit update
-        quantity.semi_implicit_update = array([1.,1.,1.,1.])
-
-        #Update with given timestep
-        timestep = 0.1
-        quantity.update(timestep)
-
-        sem = array([1.,1.,1.,1.])
-        denom = ones(4, numpy.float)-timestep*sem
-
-        x = array([1, 2, 3, 4])/denom
-        assert allclose( quantity.centroid_values, x)
-
-
-    def test_both_updates(self):
-        quantity = Quantity(self.domain4)
-
-        #Test centroids
-        centroid_values = array( [1, 2, 3, 4] )
-        quantity.set_values(centroid_values, location = 'centroids')
-        assert allclose(quantity.centroid_values, centroid_values) #Centroid
-
-        #Set explicit_update
-        explicit_update = array( [4.,3.,2.,1.] )
-        quantity.explicit_update[:,] = explicit_update
-
-        #Set semi implicit update
-        semi_implicit_update = array( [1.,1.,1.,1.] )
-        quantity.semi_implicit_update[:,] = semi_implicit_update
-
-        #Update with given timestep
-        timestep = 0.1
-        quantity.update(0.1)
-
-        denom = 1.0-timestep*semi_implicit_update
-        x = (centroid_values + timestep*explicit_update)/denom
- 
-        assert allclose( quantity.centroid_values, x)
-
-
- 
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
Index: anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py	(revision 7855)
+++ anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py	(revision 7855)
@@ -0,0 +1,620 @@
+import quantity
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from anuga_1d.base.generic_domain import Generic_domain as Domain
+#from shallow_water_domain import flux_function as domain_flux_function
+
+from anuga_1d.base.quantity import *
+
+
+
+from numpy import allclose, array, ones, zeros
+import numpy
+
+
+class Test_Quantity(unittest.TestCase):
+    def setUp(self):
+        self.points3        = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.domain3        = Domain(self.points3)
+
+
+        
+        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
+        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
+        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
+        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
+        self.domain4          = Domain(self.points4,self.boundary4)
+
+        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        self.domain10 = Domain(self.points10)
+
+        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+        self.domain6 = Domain(self.points6)
+
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creat_with_boundary(self):
+
+        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
+
+        
+    def test_creation(self):
+
+        quantity = Quantity(self.domain3)
+        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
+
+        
+        try:
+            quantity = Quantity()
+        except:
+            pass
+        else:
+            raise 'Should have raised empty quantity exception'
+
+
+        try:
+            quantity = Quantity([1,2,3])
+        except AssertionError:
+            pass
+        except:
+            raise 'Should have raised "mising domain object" error'
+
+
+    def test_creation_zeros(self):
+
+        quantity = Quantity(self.domain3)
+        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
+
+
+        quantity = Quantity(self.domain4)
+        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
+                                                 [0.,0.], [0.,0.]])
+
+
+    def test_interpolation(self):
+        quantity = Quantity(self.domain4, self.vertex_values4)
+        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
+
+
+
+    def test_interpolation2(self):
+        quantity = Quantity(self.domain4, self.vertex_values4)
+        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
+
+        quantity.extrapolate_second_order()
+
+        #print quantity.vertex_values
+        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
+                                                [ 4.5,  4.5],
+                                                [ 0.5,  0.5],
+                                                [ 1.3,  7.7]])
+
+ 
+
+ 
+    def test_boundary_allocation(self):
+        quantity = Quantity(self.domain4,
+                            [[1,2], [5,5], [0,9], [-6, 3]])
+
+
+        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
+
+
+    def test_set_values(self):
+        quantity = Quantity(self.domain4)
+
+        # Test vertices
+        quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices')
+        assert allclose(quantity.vertex_values,
+                        [[1,2], [5,5], [0,9], [-6, 3]])
+        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid
+
+
+
+        # Test unique_vertices
+        quantity.set_values([1,2,4,-5,6], location='unique_vertices')
+        assert allclose(quantity.vertex_values,
+                        [[1,2], [2,4], [4,-5], [-5,6]])
+        assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid
+
+
+        # Test centroids
+        quantity.set_values([1,2,3,4], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
+
+        # Test exceptions
+        try:
+            quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]],
+                                location = 'bas kamel tuba')
+        except:
+            pass
+
+
+        try:
+            quantity.set_values([[1,3], [0,9]])
+        except AssertionError:
+            pass
+        except:
+            raise 'should have raised Assertionerror'
+
+
+
+    def test_set_values_const(self):
+        quantity = Quantity(self.domain4)
+
+        quantity.set_values(1.0, location = 'vertices')
+
+        assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]])
+        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
+
+
+        quantity.set_values(2.0, location = 'centroids')
+
+        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
+
+
+    def test_set_values_func(self):
+        quantity = Quantity(self.domain4)
+
+        def f(x):
+            return x*x
+
+
+
+        quantity.set_values(f, location = 'vertices')
+
+
+        assert allclose(quantity.vertex_values,
+                        [[0,1], [1,6.25], [6.25,9], [9,25]])
+
+        assert allclose(quantity.centroid_values,
+                        [0.5, 3.625, 7.625, 34*0.5])
+
+
+        quantity.set_values(f, location = 'centroids')
+
+        
+        assert allclose(quantity.centroid_values,
+                        [0.25, 3.0625, 7.5625, 16.0])
+
+
+    def test_integral(self):
+        quantity = Quantity(self.domain4)
+
+        #Try constants first
+        const = 5
+        quantity.set_values(const, location = 'vertices')
+        #print 'Q', quantity.get_integral()
+
+        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
+
+        #Try with a linear function
+        def f(x):
+            return x
+
+        quantity.set_values(f, location = 'vertices')
+
+ 
+        assert allclose (quantity.centroid_values,
+        [ 0.5,   1.75,  2.75,  4.  ])
+
+        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
+                                                  [ 1.,   2.5],
+                                                  [ 2.5,  3. ],
+                                                  [ 3.,   5. ]])
+
+
+        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
+
+        assert allclose (quantity.get_integral(), ref_integral)
+
+
+
+
+
+    def test_set_values_from_quantity(self):
+
+        quantity1 = Quantity(self.domain4)
+        quantity1.set_values([0,1,2,3,4], location='unique_vertices')
+
+        assert allclose(quantity1.vertex_values,
+                        [[0,1], [1,2], [2,3], [3,4]])
+
+
+        quantity2 = Quantity(self.domain4)
+        quantity2.set_values(quantity1)
+        assert allclose(quantity2.vertex_values,
+                        [[0,1], [1,2], [2,3], [3,4]])
+
+        quantity2.set_values(2*quantity1)
+
+        assert allclose(quantity2.vertex_values,
+                        [[0,2], [2,4], [4,6], [6,8]])
+
+        quantity2.set_values(2*quantity1 + 3)
+        assert allclose(quantity2.vertex_values,
+                        [[3,5], [5,7], [7,9], [9,11]])
+
+
+
+    def test_overloading(self):
+
+        quantity1 = Quantity(self.domain4)
+        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
+                            location = 'vertices')
+
+        assert allclose(quantity1.vertex_values,
+                        [[0,1], [1,2], [2,3], [3,4]])
+
+
+        quantity2 = Quantity(self.domain4)
+        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
+                             location = 'vertices')
+
+
+
+        quantity3 = Quantity(self.domain4)
+        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
+                             location = 'vertices')
+
+
+        # Negation
+        Q = -quantity1
+        assert allclose(Q.vertex_values, -quantity1.vertex_values)
+        assert allclose(Q.centroid_values, -quantity1.centroid_values)
+
+
+        # Addition
+        Q = quantity1 + 7
+        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
+        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
+
+        Q = 7 + quantity1
+        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
+        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
+
+        Q = quantity1 + quantity2
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values + quantity2.vertex_values)
+        assert allclose(Q.centroid_values,
+                        quantity1.centroid_values + quantity2.centroid_values)
+
+        Q = quantity1 + quantity2 - 3
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values + quantity2.vertex_values - 3)
+
+        Q = quantity1 - quantity2
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values - quantity2.vertex_values)
+
+        #Scaling
+        Q = quantity1*3
+        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
+        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
+
+        Q = 3*quantity1
+        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
+
+        #Multiplication
+        Q = quantity1 * quantity2
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values)
+
+        #Linear combinations
+        Q = 4*quantity1 + 2
+        assert allclose(Q.vertex_values,
+                        4*quantity1.vertex_values + 2)
+
+        Q = quantity1*quantity2 + 2
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values + 2)
+
+        Q = quantity1*quantity2 + quantity3
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values +
+                        quantity3.vertex_values)
+        Q = quantity1*quantity2 + 3*quantity3
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values +
+                        3*quantity3.vertex_values)
+        Q = quantity1*quantity2 + 3*quantity3 + 5.0
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values +
+                        3*quantity3.vertex_values + 5)
+
+        Q = quantity1*quantity2 - quantity3
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values * quantity2.vertex_values -
+                        quantity3.vertex_values)
+        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
+        assert allclose(Q.vertex_values,
+                        1.5*quantity1.vertex_values * quantity2.vertex_values -
+                        3*quantity3.vertex_values + 5)
+
+        #Try combining quantities and arrays and scalars
+        Q = 1.5*quantity1*quantity2.vertex_values -\
+            3*quantity3.vertex_values + 5.0
+        assert allclose(Q.vertex_values,
+                        1.5*quantity1.vertex_values * quantity2.vertex_values -
+                        3*quantity3.vertex_values + 5)
+
+
+        #Powers
+        Q = quantity1**2
+        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
+
+        Q = quantity1**2 +quantity2**2
+        assert allclose(Q.vertex_values,
+                        quantity1.vertex_values**2 + \
+                        quantity2.vertex_values**2)
+
+        Q = (quantity1**2 +quantity2**2)**0.5
+        assert allclose(Q.vertex_values,
+                        (quantity1.vertex_values**2 + \
+                        quantity2.vertex_values**2)**0.5)
+
+    def test_compute_gradient(self):
+        quantity = Quantity(self.domain6)
+
+        #Set up for a gradient of (2,0) at mid triangle
+        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
+                            location = 'centroids')
+
+
+        #Gradients
+        quantity.compute_gradients()
+
+        a = quantity.gradients
+
+        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
+        
+        quantity.extrapolate_second_order()
+
+
+        assert allclose(quantity.vertex_values, [[ 1., 3. ],
+                                                 [ 4., 4. ],
+                                                 [ 4., 4. ],
+                                                 [ 4., 6.],
+                                                 [ 8.25, 11.75],
+                                                 [ 11.,  13. ]])
+
+                                                 
+
+    def test_second_order_extrapolation2(self):
+        quantity = Quantity(self.domain4)
+
+        #Set up for a gradient of (3,1), f(x) = 3x+y
+        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
+                            location = 'centroids')
+
+        #Gradients
+        quantity.compute_gradients()
+
+        a = quantity.gradients
+        
+
+        assert allclose(a[1], 2.8)
+
+        #Work out the others
+
+        quantity.extrapolate_second_order()
+        
+
+        assert allclose(quantity.vertex_values[1,0], 3.33333333)
+        assert allclose(quantity.vertex_values[1,1], 7.33333333)
+
+        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
+
+
+
+
+    def test_backup_saxpy_centroid_values(self):
+        quantity = Quantity(self.domain4)
+
+        #Set up for a gradient of (3,1), f(x) = 3x+y
+        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
+        d_values = array([1.0, 2.0, 3.0, 4.0])
+        quantity.set_values(c_values, location = 'centroids')
+
+        #Backup
+        quantity.backup_centroid_values()
+
+        #print quantity.vertex_values
+        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
+
+
+        quantity.set_values(d_values, location = 'centroids')
+
+        quantity.saxpy_centroid_values(2.0, 3.0)
+
+        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
+
+
+
+    def test_first_order_extrapolator(self):
+        quantity = Quantity(self.domain4)
+
+        centroid_values = array([1.,2.,3.,4.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
+
+        #Extrapolate
+        quantity.extrapolate_first_order()
+
+        #Check that gradient is zero
+        a = quantity.gradients
+        assert allclose(a, [0,0,0,0])
+       
+
+        #Check vertices but not edge values
+        assert allclose(quantity.vertex_values,
+                        [[1,1], [2,2], [3,3], [4,4]])
+
+
+    def test_second_order_extrapolator(self):
+        quantity = Quantity(self.domain4)
+
+        #Set up for a gradient of (3,0) at mid triangle
+        quantity.set_values([2.0, 4.0, 8.0, 2.0],
+                            location = 'centroids')
+
+
+
+        quantity.extrapolate_second_order()
+
+
+        #Assert that quantities are conserved
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             numpy.sum(quantity.vertex_values[k,:])/2.0)
+
+
+    def test_limit(self):
+        quantity = Quantity(self.domain4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
+
+        #Limit
+        quantity.limit_minmod()
+
+       
+        #cells 1 and 2 are local max and min
+        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
+        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
+
+        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
+        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
+
+
+
+    def test_find_qmax_qmin(self):
+        quantity = Quantity(self.domain4)
+
+
+        quantity.set_values([1.0, 4.0, 8.0, 2.0],
+                            location = 'centroids')
+
+
+
+        quantity.find_qmax_qmin()
+
+
+        assert allclose(quantity.qmax, [4.0, 8.0, 8.0, 8.0])
+        assert allclose(quantity.qmin, [1.0, 1.0, 2.0, 2.0])
+
+
+    def test_distribute_first_order(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        centroid_values = array([1.,2.,3.,4.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
+
+
+        #Extrapolate from centroid to vertices and edges
+        quantity.extrapolate_first_order()
+        
+        assert allclose(quantity.vertex_values,[[ 1.,  1.],
+                                                [ 2.,  2.],
+                                                [ 3.,  3.],
+                                                [ 4.,  4.]])
+ 
+
+
+    def test_distribute_second_order(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        centroid_values = array([2.,4.,8.,2.])
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
+
+
+        #Extrapolate
+        quantity.extrapolate_second_order()
+
+        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
+                                                 [ 2.,  6. ],
+                                                 [ 8.,   8. ],
+                                                 [ 6.8, -2.8]])
+
+
+    def test_update_explicit(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_update_semi_implicit(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+        #Set semi implicit update
+        quantity.semi_implicit_update = array([1.,1.,1.,1.])
+
+        #Update with given timestep
+        timestep = 0.1
+        quantity.update(timestep)
+
+        sem = array([1.,1.,1.,1.])
+        denom = ones(4, numpy.float)-timestep*sem
+
+        x = array([1, 2, 3, 4])/denom
+        assert allclose( quantity.centroid_values, x)
+
+
+    def test_both_updates(self):
+        quantity = Quantity(self.domain4)
+
+        #Test centroids
+        centroid_values = array( [1, 2, 3, 4] )
+        quantity.set_values(centroid_values, location = 'centroids')
+        assert allclose(quantity.centroid_values, centroid_values) #Centroid
+
+        #Set explicit_update
+        explicit_update = array( [4.,3.,2.,1.] )
+        quantity.explicit_update[:,] = explicit_update
+
+        #Set semi implicit update
+        semi_implicit_update = array( [1.,1.,1.,1.] )
+        quantity.semi_implicit_update[:,] = semi_implicit_update
+
+        #Update with given timestep
+        timestep = 0.1
+        quantity.update(0.1)
+
+        denom = 1.0-timestep*semi_implicit_update
+        x = (centroid_values + timestep*explicit_update)/denom
+ 
+        assert allclose( quantity.centroid_values, x)
+
+
+ 
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Quantity, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py	(revision 7855)
@@ -202,8 +202,7 @@
 
     if domain.discontinousb:
-        domain.quantities['width'].extrapolate_second_order()
-        
-
-    # FIXME (SGR): Replace with numpy.where to speed up code
+        width.extrapolate_second_order()
+        
+
     h0 = 1.0e-12
 
@@ -217,25 +216,4 @@
     w_C[:] = h_C + z_C
     
-#    for i in range(N):
-#
-#        if a_C[i] <= h0:
-#	    a_C[i] = 0.0
-#	    h_C[i] = 0.0
-#            d_C[i] = 0.0
-#            u_C[i] = 0.0
-#            w_C[i] = z_C[i]
-#        else:
-#
-#            if b_C[i]<=h0:
-#                a_C[i] = 0.0
-#                h_C[i] = 0.0
-#                d_C[i] = 0.0
-#                u_C[i] = 0.0
-#                w_C[i] = z_C[i]
-#
-#            else:
-#                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
-#                w_C[i] = h_C[i]+z_C[i]
-#                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
 		
 
@@ -260,7 +238,4 @@
     b_V  = width.vertex_values
 
-   
-    #FIXME (SGR): replace with numpy.where
-
 
     h_V[:] = w_V-z_V
@@ -271,12 +246,4 @@
     d_V[:] = u_V*a_V
 
-#    for i in range(len(h_C)):
-#        for j in range(2):
-#            h_V[i,j] = w_V[i,j]-z_V[i,j]
-#            if h_V[i,j]<h0:
-#                h_V[i,j]=0
-#                w_V[i,j]=z_V[i,j]
-#            a_V[i,j] = b_V[i,j]*h_V[i,j]
-#            d_V[i,j]=u_V[i,j]*a_V[i,j]
     
     return
Index: anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py	(revision 7855)
@@ -92,5 +92,7 @@
 #p.strip_dirs().sort_stats(-1).print_stats(20)
 
-p.sort_stats('cumulative').print_stats(10)
+p.sort_stats('cumulative').print_stats(30)
+
+print p
 
 
Index: anuga_work/development/2010-projects/anuga_1d/test_all.py
===================================================================
--- anuga_work/development/2010-projects/anuga_1d/test_all.py	(revision 7852)
+++ anuga_work/development/2010-projects/anuga_1d/test_all.py	(revision 7855)
@@ -23,7 +23,6 @@
 
 # Directories that should not be searched for test files.
-exclude_dirs = ['channel',    # Special requirements
-                '.svn',       # subversion
-                'props', 'wcprops', 'prop-base', 'text-base', 'tmp']
+exclude_dirs = ['channel',  '.svn'] 
+                
 
 
@@ -178,10 +177,10 @@
 # @brief Check that the environment is sane.
 # @note Stops here if there is an error.
-def check_anuga_import():
+def check_anuga_1d_import():
     try:
         # importing something that loads quickly
-        import anuga.anuga_exceptions
+        import anuga_1d.config
     except ImportError:
-        print "Python cannot import ANUGA module."
+        print "Python cannot import ANUGA_1D module."
         print "Check you have followed all steps of its installation."
         import sys
@@ -190,5 +189,7 @@
 
 if __name__ == '__main__':
-    check_anuga_import()
+    check_anuga_1d_import()
+
+    print 'test-all'
 
     if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
