Index: /anuga_work/development/pipeflow/generic_domain.py
===================================================================
--- /anuga_work/development/pipeflow/generic_domain.py	(revision 7822)
+++ /anuga_work/development/pipeflow/generic_domain.py	(revision 7823)
@@ -12,5 +12,5 @@
 
 
-class Generic_Domain:
+class Generic_domain:
 
     def __init__(self,
@@ -587,9 +587,13 @@
             return self.vertices[elem_id,vertex]
 
-    def get_area(self, elem_id):
+    def get_area(self, elem_id=None):
         """Return area of element id
         """
 
-        return self.areas[elem_id]
+        if elem_id is None:
+            return sum(self.areas)
+        else:
+            return self.areas[elem_id]
+
 
     def get_quantity(self, name, location='vertices', indices = None):
Index: /anuga_work/development/pipeflow/old_quantity.py
===================================================================
--- /anuga_work/development/pipeflow/old_quantity.py	(revision 7823)
+++ /anuga_work/development/pipeflow/old_quantity.py	(revision 7823)
@@ -0,0 +1,1003 @@
+"""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
+"""
+
+import numpy
+
+class Quantity:
+
+    
+    def __init__(self, domain, vertex_values=None):
+        #Initialise Quantity using optional vertex values.
+        
+        from generic_domain import Generic_Domain
+
+        msg = 'First argument in Quantity.__init__ '
+        msg += 'must be of class Domain (or a subclass thereof)'
+        assert isinstance(domain, Generic_Domain), msg
+
+        if vertex_values is None:
+            N = domain.number_of_elements
+            self.vertex_values = numpy.zeros((N, 2), numpy.float)
+        else:
+            self.vertex_values = numpy.array(vertex_values, numpy.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 = numpy.zeros(N, numpy.float)
+        self.centroid_backup_values = numpy.zeros(N, numpy.float)
+        #self.edge_values = zeros((N, 2), numpy.float)
+        #edge values are values of the ends of each interval
+      
+        #Intialise centroid values
+        self.interpolate()
+
+        
+        #Allocate space for boundary values
+        #L = len(domain.boundary)
+        self.boundary_values = numpy.zeros(2, numpy.float) #assumes no parrellism
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+        
+        N = domain.number_of_elements
+        self.explicit_update = numpy.zeros(N, numpy.float )
+        self.semi_implicit_update = numpy.zeros(N, numpy.float )
+
+        self.gradients = numpy.zeros(N, numpy.float)
+        self.qmax = numpy.zeros(self.centroid_values.shape, numpy.float)
+        self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
+
+        self.beta = domain.beta        
+
+
+    def __len__(self):
+        """
+        Returns number of intervals.
+        """
+        return self.centroid_values.shape[0]
+    
+    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.0
+
+    def set_values(self, X, location='vertices'):
+        """Set values for quantity
+
+        X: Compatible list, numpy 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 numpyal 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 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 and edge values
+            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: numpy array
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid, edges
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a numpyal 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.
+        """
+
+
+        values = numpy.array(values, numpy.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
+        #elif location == 'edges':
+        #    assert len(values.shape) == 2, 'Values array must be 2d'
+        #    msg = 'Array must be N x 2'
+        #    self.edge_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
+
+
+    def get_values(self, location='vertices', indices = None):
+        """get values for quantity
+
+        return X, Compatible list, numpy array (see below)
+        location: Where values are to be stored.
+                  Permissible options are: vertices, centroid
+                  and unique vertices. Default is 'vertices'
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a numpyal array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3
+
+        The returned values with be a list the length of indices
+        (N if indices = None).  Each value will be a list of the three
+        vertex values for this quantity.
+
+        Indices is the set of element ids that the operation applies to.
+
+        """
+        from numpy import take
+
+        if location not in ['vertices', 'centroids', 'unique vertices']:
+            msg = 'Invalid location: %s' %location
+            raise msg
+
+        import types
+        assert type(indices) in [types.ListType, types.NoneType,
+                                 numpy.ArrayType],\
+                                 'Indices must be a list or None'
+
+        if location == 'centroids':
+            if (indices ==  None):
+                indices = range(len(self))
+            return take(self.centroid_values,indices)
+        elif location == 'unique vertices':
+            if (indices ==  None):
+                indices=range(self.domain.coordinates.shape[0])
+            vert_values = []
+            #Go through list of unique vertices
+            for unique_vert_id in indices:
+                cells = self.domain.vertexlist[unique_vert_id]
+
+                #In case there are unused points
+                if cells is None:
+                    msg = 'Unique vertex not associated with cells'
+                    raise msg
+
+                # Go through all cells, vertex pairs
+                # Average the values
+                sum = 0
+                for cell_id, vertex_id in cells:
+                    sum += self.vertex_values[cell_id, vertex_id]
+                vert_values.append(sum/len(cells))
+            return numpy.array(vert_values)
+        else:
+            if (indices ==  None):
+                indices = range(len(self))
+            return numpy.take(self.vertex_values,indices)
+
+
+    def get_vertex_values(self,
+                          x=True,
+                          smooth = None,
+                          precision = None,
+                          reduction = None):
+        """Return vertex values like an OBJ format
+
+        The vertex values are returned as one sequence in the 1D numpy.float array A.
+        If requested the coordinates will be returned in 1D arrays X.
+
+        The connectivity is represented as an integer array, V, of dimension
+        M x 2, where M is the number of volumes. Each row has two indices
+        into the X, A arrays defining the element.
+
+        if smooth is True, vertex values corresponding to one common
+        coordinate set will be smoothed according to the given
+        reduction operator. In this case vertex coordinates will be
+        de-duplicated.
+
+        If no smoothings is required, vertex coordinates and values will
+        be aggregated as a concatenation of values at
+        vertices 0, vertices 1
+
+
+        Calling convention
+        if x is True:
+           X,A,V = get_vertex_values
+        else:
+           A,V = get_vertex_values
+
+        """
+
+
+        if smooth is None:
+            smooth = self.domain.smooth
+
+        if precision is None:
+            precision = numpy.float
+
+        if reduction is None:
+            reduction = self.domain.reduction
+
+        #Create connectivity
+
+        if smooth == True:
+
+            V = self.domain.get_vertices()
+            N = len(self.domain.vertexlist)
+            #N = len(self.domain.vertices)
+            A = numpy.zeros(N, precision)
+
+            #Smoothing loop
+            for k in range(N):
+                L = self.domain.vertexlist[k]
+                #L = self.domain.vertices[k]
+
+                #Go through all triangle, vertex pairs
+                #contributing to vertex k and register vertex value
+
+                if L is None: continue #In case there are unused points
+
+                contributions = []
+                for volume_id, vertex_id in L:
+                    v = self.vertex_values[volume_id, vertex_id]
+                    contributions.append(v)
+
+                A[k] = reduction(contributions)
+
+            if x is True:
+                 #X = self.domain.coordinates[:,0].astype(precision)
+                 X = self.domain.coordinates[:].astype(precision)
+                 #Y = self.domain.coordinates[:,1].astype(precision)
+
+                 #return X, Y, A, V
+                 return X, A, V
+            
+            #else:
+            return A, V
+        else:
+            #Don't smooth
+            #obj machinery moved to general_mesh
+
+            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
+            # These vert_id's will relate to the verts created below
+            #m = len(self.domain)  #Number of volumes
+            #M = 3*m        #Total number of unique vertices
+            #V = reshape(array(range(M)).astype(Int), (m,3))
+
+            #V = self.domain.get_triangles(obj=True)
+            V = self.domain.get_vertices
+            #FIXME use get_vertices, when ready
+
+            A = self.vertex_values.flat
+
+            #Do vertex coordinates
+            if x is True:
+                X = self.domain.get_vertex_coordinates()
+
+                #X = C[:,0:6:2].copy()
+                #Y = C[:,1:6:2].copy()
+
+                return X.flat, A, V
+            else:
+                return A, V
+
+    def get_integral(self):
+        """Compute the integral of quantity across entire domain
+        """
+        integral = 0
+        for k in range(self.domain.number_of_elements):
+            area = self.domain.areas[k]
+            qc = self.centroid_values[k]
+            integral += qc*area
+
+        return integral
+
+
+    def update(self, timestep):
+        """Update centroid values based on values stored in
+        explicit_update and semi_implicit_update as well as given timestep
+        """
+
+        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 numpy.sum(numpy.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.
+        """
+
+        #print 'compute_gradient'
+
+        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
+                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)
+
+
+    def compute_minmod_gradients(self):
+        """Compute gradients of piecewise linear function defined by centroids of
+        neighbouring volumes.
+        """
+
+        #print 'compute_minmod_gradients'
+        
+        def xmin(a,b):
+            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
+
+        def xmic(t,a,b):
+            return xmin(t*xmin(a,b), 0.50*(a+b) )
+
+        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
+                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 )
+        
+
+    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.
+        """
+        if self.domain.limiter == "pyvolution":
+            #Z = self.gradients
+            #print "gradients 1",Z
+            self.compute_gradients()
+            #print "gradients 2",Z
+
+            #Z = self.gradients
+            #print "gradients 1",Z
+            #self.compute_minmod_gradients()
+            #print "gradients 2", Z
+
+            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)
+            self.limit_pyvolution()
+        elif self.domain.limiter == "minmod_steve":
+            self.limit_minmod()
+        else:
+            self.limit_range()
+        
+        
+
+    def limit_minmod(self):
+        #Z = self.gradients
+        #print "gradients 1",Z
+        self.compute_minmod_gradients()
+        #print "gradients 2", Z
+
+        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_pyvolution(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
+        """
+
+        N = self.domain.number_of_elements
+        beta = self.domain.beta
+        #beta = 0.8
+
+        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 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)
+
+
+        #Diffences between centroids and maxima/minima
+        dqmax = qmax - qc
+        dqmin = qmin - qc
+
+        #Deltas between vertex and centroid values
+        dq = zeros(qv.shape, numpy.float)
+        for i in range(2):
+            dq[:,i] = qv[:,i] - qc
+
+        #Phi limiter
+        for k in range(N):
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(2):
+                r = 1.0
+                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
+                if (dq[k,i] < 0): r = dqmin[k]/dq[k,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[k,i]
+
+    def limit_range(self):
+        import sys
+
+        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
+        limiter = self.domain.limiter
+        #print limiter
+        
+        #print 'limit_range'
+        N  = self.domain.number_of_elements
+        qc = self.centroid_values
+        qv = self.vertex_values
+        C  = self.domain.centroids
+        X  = self.domain.vertices
+        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])
+                
+        dq = zeros(qv.shape, numpy.float)
+        for i in range(2):
+            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
+            
+        #Phi limiter
+        for k in range(N):
+            n0 = self.domain.neighbours[k,0]
+            n1 = self.domain.neighbours[k,1]
+            if n0 < 0:
+                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
+            elif n1 < 0:
+                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
+            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
+            #    phi = 0.0
+            else:
+                if limiter == "minmod":
+                    phi = minmod(beta_p[k],beta_m[k])
+
+                elif limiter == "minmod_kurganov":#Change this
+                    # Also known as monotonized central difference limiter
+                    # if theta = 2.0
+                    theta = 2.0 
+                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
+                elif limiter == "superbee":
+                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
+                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
+                    phi = maxmod(slope1,slope2)
+
+                elif limiter == "vanleer":
+                    phi = vanleer(beta_p[k],beta_m[k])
+
+                elif limiter == "vanalbada":
+                    phi = vanalbada(beta_m[k],beta_p[k])
+            
+            for i in range(2):
+                qv[k,i] = qc[k] + phi*dq[k,i]
+
+    def limit_steve_slope(self):    
+
+        import sys
+        from util import minmod, minmod_kurganov, maxmod, vanleer
+
+        N = self.domain.number_of_elements
+        limiter = self.domain.limiter
+        limiter_type = self.domain.limiter_type
+            
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        #Find min and max of this and neighbour's centroid values
+        beta_p = numpy.zeros(N,numpy.float)
+        beta_m = numpy.zeros(N,numpy.float)
+        beta_x = numpy.zeros(N,numpy.float)
+        C = self.domain.centroids
+        X = self.domain.vertices
+
+        for k in range(N):
+        
+            n0 = self.domain.neighbours[k,0]
+            n1 = self.domain.neighbours[k,1]
+            
+            if (n0 >= 0) & (n1 >= 0):
+                # Check denominator not zero
+                if (qc[k+1]-qc[k]) == 0.0:
+                    beta_p[k] = numpy.float(sys.maxint)
+                    beta_m[k] = numpy.float(sys.maxint)
+                else:
+                    #STEVE LIMIT
+                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
+                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
+
+        #Deltas between vertex and centroid values
+        dq = zeros(qv.shape, numpy.float)
+        for i in range(2):
+            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
+            
+        #Phi limiter
+        for k in range(N):
+                
+            phi = 0.0
+            if limiter == "flux_minmod":
+                #FLUX MINMOD
+                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
+            elif limiter == "flux_superbee":
+                #FLUX SUPERBEE
+                phi = max(0.0,min(1.0,2.0*beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0*beta_p[k]),min(2.0,beta_p[k]))-1.0
+            elif limiter == "flux_muscl":
+                #FLUX MUSCL
+                phi = max(0.0,min(2.0,2.0*beta_m[k],2.0*beta_p[k],0.5*(beta_m[k]+beta_p[k])))
+            elif limiter == "flux_vanleer":
+                #FLUX VAN LEER
+                phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))+(beta_p[k]+abs(beta_p[k]))/(1.0+abs(beta_p[k]))-1.0
+                
+                #Then update using phi limiter
+                n = self.domain.neighbours[k,1]
+                if n>=0:
+                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
+                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
+                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
+                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
+                else:
+                    qv[k,i] = qc[k]
+
+    def backup_centroid_values(self):
+        # Call correct module function
+        # (either from this module or C-extension)
+        #backup_centroid_values(self)
+
+        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
+
+    def saxpy_centroid_values(self,a,b):
+        # Call correct module function
+        # (either from this module or C-extension)
+        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
+        
+
+
+def newLinePlot(title='Simple Plot'):
+    import pylab as g
+    g.ion()
+    g.hold(False)
+    g.title(title)
+    g.xlabel('x')
+    g.ylabel('y')
+    
+
+def linePlot(x,y):
+    import pylab as g
+    g.plot(x.flat,y.flat)
+
+
+def closePlots():
+    import pylab as g
+    g.close('all')
+    
+if __name__ == "__main__":
+    #from domain import Domain
+    from shallow_water_domain import Domain     
+    from numpy import arange
+    
+    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 = 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
+
+    class FunClass:
+        def __init__(self,value):
+            self.value = value
+
+        def __call__(self,x):
+            return self.value*(x**2)
+
+
+    fun = FunClass(1.0)
+    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_minmod()
+
+    newLinePlot('plots')
+    linePlot(Xc,Qc)
+    raw_input('press return')
+
+    points2 = arange(10)
+    D2 = Domain(points2)
+
+    Q2 = Quantity(D2)
+    Q2.set_values(fun,'vertices')
+    Xc = Q2.domain.vertices
+    Qc = Q2.vertex_values
+    linePlot(Xc,Qc)
+    raw_input('press return')
+
+    
+    Q2.extrapolate_second_order()
+    #Q2.limit_minmod()
+    Xc = Q2.domain.vertices
+    Qc = Q2.vertex_values
+    print Q2.centroid_values
+    print Qc
+    linePlot(Xc,Qc)
+    raw_input('press return')
+
+
+    for i in range(10):
+        import pylab as g
+        g.hold(True)
+        fun = FunClass(i/10.0)
+        Q2.set_values(fun,'centroids')
+        Q2.extrapolate_second_order()
+        #Q2.limit_minmod()
+        Qc = Q2.vertex_values
+        linePlot(Xc,Qc)
+        raw_input('press return')
+
+    raw_input('press return to quit')
+closePlots()
Index: /anuga_work/development/pipeflow/old_sww_comp_flux_ext.c
===================================================================
--- /anuga_work/development/pipeflow/old_sww_comp_flux_ext.c	(revision 7823)
+++ /anuga_work/development/pipeflow/old_sww_comp_flux_ext.c	(revision 7823)
@@ -0,0 +1,312 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+const double pi = 3.14159265358979;
+
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+//Innermost flux function (using w=z+h)
+int _flux_function(double *q_left, double *q_right,
+		   double normals, double g, double epsilon, double h0,
+		   double *edgeflux, double *max_speed) {
+		
+		int i;
+		double flux_left[2], flux_right[2];
+		double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
+		double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
+		double z, s_max, s_min, denom;
+		
+		//printf("h0 = %f \n",h0);
+		w_left  = q_left[0];
+		uh_left = q_left[1]*normals;
+		h_left  = q_left[2];
+		z_left  = q_left[3];
+		u_left  = q_left[4]*normals;
+	
+		w_right  = q_right[0];
+		uh_right = q_right[1]*normals;
+		h_right  = q_right[2];
+		z_right  = q_right[3];
+		u_right  = q_right[4]*normals;		
+	  
+		z = (z_left+z_right)/2.0;		   
+				   
+		soundspeed_left = sqrt(g*h_left);
+		soundspeed_right = sqrt(g*h_right);
+		
+		s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+		if (s_max < 0.0) s_max = 0.0;
+	
+		s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+		if (s_min > 0.0) s_min = 0.0;
+		
+		
+		// Flux formulas
+		flux_left[0] = u_left*h_left;
+		flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
+
+		flux_right[0] = u_right*h_right;
+		flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
+
+		// Flux computation
+		denom = s_max-s_min;
+		if (denom < epsilon) {
+			for (i=0; i<2; i++) edgeflux[i] = 0.0;
+			*max_speed = 0.0;
+		} else {
+			edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
+			edgeflux[0] += s_max*s_min*(q_right[0]-q_left[0]);
+			edgeflux[0] /= denom;
+			edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
+			edgeflux[1] += s_max*s_min*(q_right[1]-q_left[1]);
+			edgeflux[1] /= denom;
+			edgeflux[1] *= normals;
+    
+			// Maximal wavespeed
+			*max_speed = max(fabs(s_max), fabs(s_min));
+		}  
+		return 0;		
+}
+		
+		
+	
+	
+// Computational function for flux computation
+double _compute_fluxes_ext(
+			   double cfl,
+			   double timestep,
+			   double epsilon,
+			   double g,
+			   double h0,
+			   long* neighbours,
+			   long* neighbour_vertices,
+			   double* normals,
+			   double* areas,
+			   double* stage_edge_values,
+			   double* xmom_edge_values,
+			   double* bed_edge_values,
+			   double* height_edge_values,
+			   double* velocity_edge_values,
+			   double* stage_boundary_values,
+			   double* xmom_boundary_values,
+			   double* bed_boundary_values,
+			   double* height_boundary_values,
+			   double* velocity_boundary_values,
+			   double* stage_explicit_update,
+			   double* xmom_explicit_update,
+			   int number_of_elements,
+			   double* max_speed_array) {
+		
+                double flux[2], ql[5], qr[5], edgeflux[2];
+		double max_speed, normal;
+		int k, i, ki, n, m, nm=0;
+		
+		
+		for (k=0; k<number_of_elements; k++) {
+			flux[0] = 0.0;
+			flux[1] = 0.0;
+			
+			for (i=0; i<2; i++) {
+				ki = k*2+i;
+				
+				ql[0] = stage_edge_values[ki];
+				ql[1] = xmom_edge_values[ki];
+				ql[2] = bed_edge_values[ki];
+				ql[3] = height_edge_values[ki];
+				ql[4] = velocity_edge_values[ki];
+				
+				n = neighbours[ki];
+				if (n<0) {
+					m = -n-1;
+					qr[0] = stage_boundary_values[m];
+					qr[1] = xmom_boundary_values[m];
+					qr[2] = bed_edge_values[m];
+					qr[3] = height_edge_values[m];
+					qr[4] = velocity_edge_values[m];
+				} else {
+					m = neighbour_vertices[ki];
+					nm = n*2+m;
+					qr[0] = stage_edge_values[nm];
+					qr[1] = xmom_edge_values[nm];
+					qr[2] = bed_edge_values[nm];
+					qr[3] = height_edge_values[nm];
+					qr[4] = velocity_edge_values[nm];
+				}
+				
+				normal = normals[ki];
+				_flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
+				flux[0] -= edgeflux[0];
+				flux[1] -= edgeflux[1];
+				
+				// Update timestep based on edge i and possibly neighbour n
+				if (max_speed > epsilon) {
+				    // Original CFL calculation
+				    
+				    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
+				    if (n>=0) {
+					timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
+				    }
+				}
+			} // End edge i (and neighbour n)
+			flux[0] /= areas[k];
+			stage_explicit_update[k] = flux[0];
+			flux[1] /= areas[k];
+			xmom_explicit_update[k] = flux[1];
+			
+			//Keep track of maximal speeds
+			max_speed_array[k]=max_speed;
+		}
+		return timestep;	
+	}
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
+
+    PyObject
+	*domain,
+	*stage,
+	*xmom,
+	*bed,
+	*height,
+	*velocity;
+
+    PyArrayObject
+	*neighbours,
+	*neighbour_vertices,
+	*normals,
+	*areas,
+	*stage_vertex_values,
+	*xmom_vertex_values,
+	*bed_vertex_values,
+	*height_vertex_values,
+	*velocity_vertex_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*bed_boundary_values,
+	*height_boundary_values,
+	*velocity_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+
+  double timestep, epsilon, g, h0, cfl;
+  int number_of_elements;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOOOO",
+			&timestep,
+			&domain,
+			&stage,
+			&xmom,
+			&bed,
+			&height,
+			&velocity)) {
+      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
+      return NULL;
+  }
+
+
+    epsilon           = get_python_double(domain,"epsilon");
+    g                 = get_python_double(domain,"g");
+    h0                = get_python_double(domain,"h0");
+    cfl               = get_python_double(domain,"cfl");
+
+
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
+    normals           = get_consecutive_array(domain, "normals");
+    areas             = get_consecutive_array(domain, "areas");
+    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
+
+    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");
+    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");
+    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
+    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
+    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
+
+    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");
+    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");
+    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
+    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
+    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");
+    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");
+
+    number_of_elements = stage_vertex_values -> dimensions[0];
+
+    // Call underlying flux computation routine and update
+    // the explicit update arrays
+    timestep = _compute_fluxes_ext(
+	cfl,
+	timestep,
+	epsilon,
+	g,
+	h0,
+	(long*) neighbours -> data,
+	(long*) neighbour_vertices -> data,
+	(double*) normals -> data,
+	(double*) areas -> data,
+	(double*) stage_vertex_values -> data,
+	(double*) xmom_vertex_values -> data,
+	(double*) bed_vertex_values -> data,
+	(double*) height_vertex_values -> data,
+	(double*) velocity_vertex_values -> data,
+	(double*) stage_boundary_values -> data,
+	(double*) xmom_boundary_values -> data,
+	(double*) bed_boundary_values -> data,
+	(double*) height_boundary_values -> data,
+	(double*) velocity_boundary_values -> data,
+	(double*) stage_explicit_update -> data,
+	(double*) xmom_explicit_update -> data,
+	number_of_elements,
+	(double*) max_speed_array -> data);
+
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_vertices);
+  Py_DECREF(normals);
+  Py_DECREF(areas);
+  Py_DECREF(stage_vertex_values);
+  Py_DECREF(xmom_vertex_values);
+  Py_DECREF(bed_vertex_values);
+  Py_DECREF(height_vertex_values);
+  Py_DECREF(velocity_vertex_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(bed_boundary_values);
+  Py_DECREF(height_boundary_values);
+  Py_DECREF(velocity_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(max_speed_array);
+
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+
+static struct PyMethodDef MethodTable[] = {
+  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initsww_comp_flux_ext(void){
+  Py_InitModule("sww_comp_flux_ext", MethodTable);
+  import_array();
+}
Index: /anuga_work/development/pipeflow/old_test_quantity.py
===================================================================
--- /anuga_work/development/pipeflow/old_test_quantity.py	(revision 7823)
+++ /anuga_work/development/pipeflow/old_test_quantity.py	(revision 7823)
@@ -0,0 +1,447 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from shallow_water_domain import *
+from shallow_water_domain import flux_function as domain_flux_function
+
+
+
+import numpy
+
+
+class Test_Shallow_Water(unittest.TestCase):
+    def setUp(self):
+        self.points = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        domain = Domain(self.points)
+        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])
+
+    def test_compute_fluxes(self):
+        """
+        Compare shallow_water_domain flux calculation against a previous
+        Python implementation (defined in this file)
+        """
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        stage_ud, xmom_ud = local_compute_fluxes(domain)
+
+        domain.distribute_to_vertices_and_edges()
+        domain.compute_fluxes()
+
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
+
+        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
+        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = numpy.array([1.0, 2.0],numpy.float)
+        qr = numpy.array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        #This assumes h0 = 1.0e-3!!
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        #print maxspeed
+        #print edgeflux
+        
+        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
+        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+        normal = -1.0
+        ql = numpy.array([1.0, 2.0],numpy.float)
+        qr = numpy.array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+
+        #print maxspeed
+        #print edgeflux        
+        
+        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
+        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+    def test_domain_flux_function(self):
+        normal = 1.0
+        ql = numpy.array([1.0, 2.0],numpy.float)
+        qr = numpy.array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+        #print edgeflux
+
+
+
+        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
+
+        #print domainedgeflux
+
+        assert numpy.allclose(domainedgeflux, edgeflux)
+        assert numpy.allclose(domainmaxspeed, maxspeed)
+
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+
+        def slope_one(x):
+            return x
+
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
+
+
+    def xtest_evolve_first_order(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+        
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 1
+        domain.set_timestepping_method('euler')
+        
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            #pass
+        
+            print domain.quantities['stage'].vertex_values
+            print domain.quantities['elevation'].vertex_values
+            print domain.quantities['xmomentum'].vertex_values
+          
+        
+            print domain.quantities['stage'].centroid_values  
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values    
+
+        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+    def xtest_evolve_euler_second_order_space(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.set_timestepping_method('rk2')
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+    def xtest_evolve_second_order_space_time(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.set_timestepping_method('rk3')
+
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+
+#==============================================================================
+
+def local_compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    tmp0 = zeros(N,Float)
+    tmp1 = zeros(N,Float)
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    #    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
+    #   ymom =  Ymom.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    #print 'stage edge values', stage
+    #print 'xmom edge values', xmom
+    #print 'bed values', bed
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    #print 'stage_bdry',stage_bdry
+    #print 'xmom_bdry', xmom_bdry
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, Float) #Work numpy.array for summing up fluxes
+    flux = zeros(2, Float) #Work numpy.array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work numpy.array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            #print 'ql',ql
+            #print 'qr',qr
+            
+
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+
+        #Stage.explicit_update[k] = flux[0]
+        tmp0[k] = flux[0]
+        tmp1[k] = flux[1]
+
+
+    return tmp0, tmp1
+
+
+def local_flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/(h_left +  h0/h_left)
+
+
+    uh_left = u_left*h_left
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/(h_right + h0/h_right)
+
+    uh_right = u_right*h_right
+    
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = numpy.array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = numpy.array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = numpy.array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = numpy.array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = numpy.array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
+
+
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /anuga_work/development/pipeflow/old_test_sww_domain_vel.py
===================================================================
--- /anuga_work/development/pipeflow/old_test_sww_domain_vel.py	(revision 7823)
+++ /anuga_work/development/pipeflow/old_test_sww_domain_vel.py	(revision 7823)
@@ -0,0 +1,517 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from shallow_water_vel_domain import *
+
+
+from numpy import allclose, array, ones, numpy.float, maximum, zeros
+
+
+class Test_Shallow_Water(unittest.TestCase):
+    def setUp(self):
+        self.points = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        domain = Domain(self.points)
+        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
+
+    def test_reflective_boundary(self):
+        """
+        Test setting reflective boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test reflective'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [-6.0, -6.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [0.0, 0.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [-3.0, -3.0])
+
+        
+
+    def test_dirichlet_boundary(self):
+        """
+        Test setting dirichlet boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Dirichlet_boundary([3.0, 8.0, 1.0, 2.0, 4.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test dirichlet'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [3.0, 3.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [8.0, 8.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [1.0, 1.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [4.0, 4.0])
+
+
+    def test_compute_fluxes(self):
+        """
+        Compare shallow_water_domain flux calculation against a previous
+        Python implementation (defined in this file)
+        """
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_boundary({'exterior' : Dirichlet_boundary([0.0, 0.0, 0.0, 0.0, 0.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        stage_ud, xmom_ud = local_compute_fluxes(domain)
+
+        domain.compute_fluxes()
+        
+        print domain.quantities['stage'].vertex_values
+        print domain.quantities['xmomentum'].vertex_values
+        print domain.quantities['elevation'].vertex_values
+        print domain.quantities['height'].vertex_values
+        print domain.quantities['velocity'].vertex_values
+
+        print domain.quantities['stage'].boundary_values
+        print domain.quantities['xmomentum'].boundary_values
+        print domain.quantities['elevation'].boundary_values
+        print domain.quantities['height'].boundary_values
+        print domain.quantities['velocity'].boundary_values
+        
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
+
+        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
+        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],numpy.float)
+        qr = array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        #This assumes h0 = 1.0e-3!!
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        #print maxspeed
+        #print edgeflux
+        
+        assert allclose(array([2.0, 8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],numpy.float)
+        qr = array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+
+        #print maxspeed
+        #print edgeflux        
+        
+        assert allclose(array([-2.0, -8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+
+        def slope_one(x):
+            return x
+
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert allclose( array([-34.3, -24.5, -14.7], numpy.float), domain.quantities['xmomentum'].explicit_update )
+
+
+    def test_evolve_first_order(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+        
+        def slope_square(x):
+            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 1
+        domain.set_timestepping_method('euler')
+        
+        yieldstep=0.25
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        ##	print domain.quantities['stage'].vertex_values 
+        ##	print domain.quantities['elevation'].vertex_values 
+        ##        print domain.quantities['xmomentum'].vertex_values
+        ##  
+        ##
+        ##        print domain.quantities['stage'].centroid_values  
+        ##	print domain.quantities['elevation'].centroid_values    
+        ##        print domain.quantities['xmomentum'].centroid_values    
+
+        #assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        #assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+##     def test_evolve_euler_second_order_space(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk2')
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+##     def test_evolve_second_order_space_time(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk3')
+
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+
+#==============================================================================
+
+def local_compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from numpy import zeros, numpy.float
+
+    N = domain.number_of_elements
+
+    tmp0 = zeros(N,numpy.float)
+    tmp1 = zeros(N,numpy.float)
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    #    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
+    #   ymom =  Ymom.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    #print 'stage edge values', stage
+    #print 'xmom edge values', xmom
+    #print 'bed values', bed
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    #print 'stage_bdry',stage_bdry
+    #print 'xmom_bdry', xmom_bdry
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, numpy.float) #Work array for summing up fluxes
+    flux = zeros(2, numpy.float) #Work array for summing up fluxes
+    ql = zeros(2, numpy.float)
+    qr = zeros(2, numpy.float)
+
+    #Loop
+    timestep = numpy.float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            #print 'ql',ql
+            #print 'qr',qr
+            
+
+            edgeflux, max_speed = local_flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+
+        #Stage.explicit_update[k] = flux[0]
+        tmp0[k] = flux[0]
+        tmp1[k] = flux[1]
+
+
+    return tmp0, tmp1
+
+
+def local_flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numpyal vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+    from numpy import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/(h_left +  h0/h_left)
+
+
+    uh_left = u_left*h_left
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/(h_right + h0/h_right)
+
+    uh_right = u_right*h_right
+    
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')    
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
+
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
+    #print "restricted test"
+    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /anuga_work/development/pipeflow/quantity.py
===================================================================
--- /anuga_work/development/pipeflow/quantity.py	(revision 7822)
+++ /anuga_work/development/pipeflow/quantity.py	(revision 7823)
@@ -1,3 +1,4 @@
-"""Class Quantity - Implements values at each 1d element
+"""
+Class Quantity - Implements values at each 1d element
 
 To create:
@@ -17,4 +18,6 @@
 import numpy
 
+
+
 class Quantity:
 
@@ -23,9 +26,9 @@
         #Initialise Quantity using optional vertex values.
         
-        from generic_domain import Generic_Domain
+        from generic_domain import Generic_domain
 
         msg = 'First argument in Quantity.__init__ '
         msg += 'must be of class Domain (or a subclass thereof)'
-        assert isinstance(domain, Generic_Domain), msg
+        assert isinstance(domain, Generic_domain), msg
 
         if vertex_values is None:
@@ -52,5 +55,5 @@
         self.centroid_values = numpy.zeros(N, numpy.float)
         self.centroid_backup_values = numpy.zeros(N, numpy.float)
-        #self.edge_values = zeros((N, 2), numpy.float)
+        #self.edge_values = numpy.zeros((N, 2), numpy.float)
         #edge values are values of the ends of each interval
       
@@ -58,4 +61,5 @@
         self.interpolate()
 
+
         
         #Allocate space for boundary values
@@ -82,5 +86,176 @@
         """
         return self.centroid_values.shape[0]
-    
+
+    def __neg__(self):
+        """Negate all values in this quantity giving meaning to the
+        expression -Q where Q is an instance of class Quantity
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(-self.vertex_values)
+        return Q
+
+
+
+    def __add__(self, other):
+        """Add to self anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
+        return result
+
+    def __radd__(self, other):
+        """Handle cases like 7+Q, where Q is an instance of class Quantity
+        """
+
+        return self + other
+
+    def __sub__(self, other):
+        return self + -other            # Invoke self.__neg__()
+
+    def __mul__(self, other):
+        """Multiply self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The product of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be products from q1 and q2
+        result.vertex_values = self.vertex_values * Q.vertex_values
+        result.centroid_values = self.centroid_values * Q.centroid_values
+
+        return result
+
+    def __rmul__(self, other):
+        """Handle cases like 3*Q, where Q is an instance of class Quantity
+        """
+
+        return self * other
+
+    def __div__(self, other):
+        """Divide self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+
+        Zero division is dealt with by adding an epsilon to the divisore
+        FIXME (Ole): Replace this with native INF once we migrate to NumPy
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The quotient of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be quotient of q1 and q2
+        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
+        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
+
+        return result
+
+    def __rdiv__(self, other):
+        """Handle cases like 3/Q, where Q is an instance of class Quantity
+        """
+
+        return self / other
+
+    def __pow__(self, other):
+        """Raise quantity to (numerical) power
+
+        As with __mul__ vertex values are processed entry by entry
+        while centroid and edge values are re-interpolated.
+
+        Example using __pow__:
+          Q = (Q1**2 + Q2**2)**0.5
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The power of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be correct
+        result.vertex_values = self.vertex_values ** other
+        result.centroid_values = self.centroid_values ** other
+
+        return result
+
+    def set_values_from_numeric(self, numeric):
+        
+        x = numpy.array([1.0,2.0])
+        y = [1.0,2.0]
+
+        if type(numeric) == type(y):
+            self.set_values_from_array(numeric)
+        elif type(numeric) == type(x):
+            self.set_values_from_array(numeric)
+        elif callable(numeric):
+            self.set_values_from_function(numeric)
+        elif isinstance(numeric, Quantity):
+            self.set_values_from_quantity(numeric)
+        else:   # see if it's coercible to a float (float, int or long, etc)
+            try:
+                numeric = float(numeric)
+            except ValueError:
+                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
+                raise Exception(msg)
+            self.set_values_from_constant(numeric)
+
+    def set_values_from_constant(self,numeric):
+
+        self.vertex_values[:,:]   = numeric
+        self.centroid_values[:,] = numeric
+
+
+    def set_values_from_array(self,numeric):
+
+        self.vertex_values[:,:]   = numeric
+        self.interpolate()
+
+
+    def set_values_from_quantity(self,quantity):
+
+        self.vertex_values[:,:]   = quantity.vertex_values
+        self.centroid_values[:,] = quantity.centroid_values
+
+    def set_values_from_function(self,function):
+
+        self.vertex_values[:,:]   = map(function, self.domain.vertices)
+        self.centroid_values[:,] = map(function, self.domain.centroids)
+
+
     def interpolate(self):
         """
@@ -96,8 +271,11 @@
             self.centroid_values[i] = (v0 + v1)/2.0
 
+
+
+
     def set_values(self, X, location='vertices'):
         """Set values for quantity
 
-        X: Compatible list, numpy array (see below), constant or function
+        X: Compatible list, Numeric array (see below), constant or function
         location: Where values are to be stored.
                   Permissible options are: vertices, centroid
@@ -105,5 +283,5 @@
 
         In case of location == 'centroid' the dimension values must
-        be a list of a numpyal array of length N, N being the number
+        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
 
@@ -119,11 +297,11 @@
         """
 
-        if location not in ['vertices', 'centroids']:
-            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
-            raise msg
+        if location not in ['vertices', 'centroids', 'unique_vertices']:
+            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
+            raise Exception(msg)
 
         if X is None:
             msg = 'Given values are None'
-            raise msg
+            raise Exception(msg)
 
         import types
@@ -134,8 +312,10 @@
           
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
-            if location == 'centroids':
-                self.centroid_values[:] = X
-            else:
-                self.vertex_values[:] = X
+            
+            self.centroid_values[:,] = float(X)
+            self.vertex_values[:,:] = float(X)
+
+        elif isinstance(X, Quantity):
+            self.set_array_values(X.vertex_values, location)
 
         else:
@@ -143,7 +323,10 @@
             self.set_array_values(X, location)
 
-        if location == 'vertices':
-            #Intialise centroid and edge values
+        if location == 'vertices' or location == 'unique_vertices':
+            #Intialise centroid
             self.interpolate()
+
+        if location == 'centroid':
+            self.extrapolate_first_order()
 
 
@@ -176,12 +359,14 @@
         """Set values for quantity
 
-        values: numpy array
+        values: Numeric array
         location: Where values are to be stored.
-                  Permissible options are: vertices, centroid, edges
+                  Permissible options are: vertices, centroid, unique_vertices
                   Default is "vertices"
 
         In case of location == 'centroid' the dimension values must
-        be a list of a numpyal array of length N, N being the number
-        of elements in the mesh. Otherwise it must be of dimension Nx2
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. If location == 'unique_vertices' the
+        dimension values must be a list or a Numeric array of length N+1.
+        Otherwise it must be of dimension Nx2
 
         The values will be stored in elements following their
@@ -195,24 +380,39 @@
 
 
-        values = numpy.array(values, numpy.float)
+        values = numpy.array(values).astype(numpy.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':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             assert len(values.shape) == 1, 'Values array must be 1d'
-            self.centroid_values = values
-        #elif location == 'edges':
-        #    assert len(values.shape) == 2, 'Values array must be 2d'
-        #    msg = 'Array must be N x 2'
-        #    self.edge_values = values
-        else:
+
+            for i in range(N):
+                self.centroid_values[i] = values[i]
+
+            self.vertex_values[:,0] = values
+            self.vertex_values[:,1] = values
+ 
+        if location == 'vertices':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             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
+            self.vertex_values[:,:] = values[:,:]
+
+        if location == 'unique_vertices':
+            msg = 'Number of values must match number of elements +1'
+            assert values.shape[0] == N+1, msg
+            assert len(values.shape) == 1, 'Values array must be 1d'
+
+            self.vertex_values[:,0] = values[0:-1]
+            self.vertex_values[:,1] = values[1:N+1]  
+
+
+            
 
 
@@ -220,5 +420,5 @@
         """get values for quantity
 
-        return X, Compatible list, numpy array (see below)
+        return X, Compatible list, Numeric array (see below)
         location: Where values are to be stored.
                   Permissible options are: vertices, centroid
@@ -226,5 +426,5 @@
 
         In case of location == 'centroids' the dimension values must
-        be a list of a numpyal array of length N, N being the number
+        be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3
 
@@ -236,11 +436,11 @@
 
         """
-        from numpy import take
+
 
         if location not in ['vertices', 'centroids', 'unique vertices']:
             msg = 'Invalid location: %s' %location
-            raise msg
-
-        import types
+            raise Exception(msg)
+
+        import types, numpy
         assert type(indices) in [types.ListType, types.NoneType,
                                  numpy.ArrayType],\
@@ -262,5 +462,5 @@
                 if cells is None:
                     msg = 'Unique vertex not associated with cells'
-                    raise msg
+                    raise Exception(msg)
 
                 # Go through all cells, vertex pairs
@@ -274,5 +474,5 @@
             if (indices ==  None):
                 indices = range(len(self))
-            return numpy.take(self.vertex_values,indices)
+            return take(self.vertex_values,indices)
 
 
@@ -284,5 +484,5 @@
         """Return vertex values like an OBJ format
 
-        The vertex values are returned as one sequence in the 1D numpy.float array A.
+        The vertex values are returned as one sequence in the 1D float array A.
         If requested the coordinates will be returned in 1D arrays X.
 
@@ -308,4 +508,6 @@
 
         """
+
+
 
 
@@ -399,4 +601,6 @@
         """
 
+
+
         N = self.centroid_values.shape[0]
 
@@ -407,7 +611,7 @@
         denominator = numpy.ones(N, numpy.float)-timestep*self.semi_implicit_update
 
-        if numpy.sum(numpy.equal(denominator, 0.0)) > 0.0:
+        if sum(numpy.equal(denominator, 0.0)) > 0.0:
             msg = 'Zero division in semi implicit update. Call Stephen :-)'
-            raise msg
+            raise Exception(msg)
         else:
             #Update conserved_quantities from semi implicit updates
@@ -420,5 +624,6 @@
         """
 
-        #print 'compute_gradient'
+
+
 
         N = self.centroid_values.shape[0]
@@ -510,4 +715,6 @@
 
         #print 'compute_minmod_gradients'
+        from numpy import sign
+
         
         def xmin(a,b):
@@ -516,4 +723,6 @@
         def xmic(t,a,b):
             return xmin(t*xmin(a,b), 0.50*(a+b) )
+
+
 
         N = self.centroid_values.shape[0]
@@ -702,4 +911,5 @@
         """
 
+
         N = self.domain.number_of_elements
         beta = self.domain.beta
@@ -729,5 +939,5 @@
 
         #Deltas between vertex and centroid values
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] = qv[:,i] - qc
@@ -753,13 +963,14 @@
 
         from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
+
         limiter = self.domain.limiter
         #print limiter
         
         #print 'limit_range'
-        N  = self.domain.number_of_elements
+        N = self.domain.number_of_elements
         qc = self.centroid_values
         qv = self.vertex_values
-        C  = self.domain.centroids
-        X  = self.domain.vertices
+        C = self.domain.centroids
+        X = self.domain.vertices
         beta_p = numpy.zeros(N,numpy.float)
         beta_m = numpy.zeros(N,numpy.float)
@@ -777,5 +988,5 @@
                 beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
                 
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
@@ -817,4 +1028,5 @@
 
         import sys
+
         from util import minmod, minmod_kurganov, maxmod, vanleer
 
@@ -841,6 +1053,6 @@
                 # Check denominator not zero
                 if (qc[k+1]-qc[k]) == 0.0:
-                    beta_p[k] = numpy.float(sys.maxint)
-                    beta_m[k] = numpy.float(sys.maxint)
+                    beta_p[k] = float(sys.maxint)
+                    beta_m[k] = float(sys.maxint)
                 else:
                     #STEVE LIMIT
@@ -849,5 +1061,5 @@
 
         #Deltas between vertex and centroid values
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
@@ -885,10 +1097,10 @@
         #backup_centroid_values(self)
 
-        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
+        self.centroid_backup_values[:,] = (self.centroid_values).astype('f')
 
     def saxpy_centroid_values(self,a,b):
         # Call correct module function
         # (either from this module or C-extension)
-        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
+        self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
         
 
@@ -914,6 +1126,5 @@
 if __name__ == "__main__":
     #from domain import Domain
-    from shallow_water_domain import Domain     
-    from numpy import arange
+    from sww_domain import Domain
     
     points1 = [0.0, 1.0, 2.0, 3.0]
@@ -968,5 +1179,5 @@
     raw_input('press return')
 
-    points2 = arange(10)
+    points2 = numpy.arange(10)
     D2 = Domain(points2)
 
Index: /anuga_work/development/pipeflow/sww_comp_flux_ext.c
===================================================================
--- /anuga_work/development/pipeflow/sww_comp_flux_ext.c	(revision 7822)
+++ /anuga_work/development/pipeflow/sww_comp_flux_ext.c	(revision 7823)
@@ -10,30 +10,84 @@
 
 
+/* double max(double a, double b) { */
+/* 	double z; */
+/* 	z=(a>b)?a:b; */
+/* 	return z;} */
+	
+/* double min(double a, double b) { */
+/* 	double z; */
+/* 	z=(a<b)?a:b; */
+/* 	return z;} */
+
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+		      double *h, 
+		      double epsilon, 
+		      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+     u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);    
+  }
+  
+
+  // Adjust momentum to be consistent with speed
+  *uh = u * *h;
+  
+  return u;
+}
+
+
+
+
 //Innermost flux function (using w=z+h)
 int _flux_function(double *q_left, double *q_right,
+        double z_left, double z_right,
 		   double normals, double g, double epsilon, double h0,
-		   double *edgeflux, double *max_speed) {
+		double *edgeflux, double *max_speed) {
 		
 		int i;
-		double flux_left[2], flux_right[2];
-		double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
-		double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
-		double z, s_max, s_min, denom;
+		double ql[2], qr[2], flux_left[2], flux_right[2];
+		double z, w_left, h_left, uh_left, soundspeed_left, u_left;
+		double w_right, h_right, uh_right, soundspeed_right, u_right;
+		double s_max, s_min, denom;
 		
 		//printf("h0 = %f \n",h0);
-		w_left  = q_left[0];
-		uh_left = q_left[1]*normals;
-		h_left  = q_left[2];
-		z_left  = q_left[3];
-		u_left  = q_left[4]*normals;
-	
-		w_right  = q_right[0];
-		uh_right = q_right[1]*normals;
-		h_right  = q_right[2];
-		z_right  = q_right[3];
-		u_right  = q_right[4]*normals;		
+		ql[0] = q_left[0];
+		ql[1] = q_left[1];
+		ql[1] = ql[1]*normals;
+	
+		qr[0] = q_right[0];
+		qr[1] = q_right[1];
+		qr[1] = qr[1]*normals;
 	  
 		z = (z_left+z_right)/2.0;		   
 				   
+		//w_left = ql[0];
+		//h_left = w_left-z;
+		//uh_left = ql[1];
+		
+
+
+		// Compute speeds in x-direction
+		w_left = ql[0];          
+		h_left = w_left-z;
+		uh_left = ql[1];
+
+		u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
+
+		w_right = qr[0];
+		h_right = w_right-z;
+		uh_right = qr[1];
+
+		u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
+  
 		soundspeed_left = sqrt(g*h_left);
 		soundspeed_right = sqrt(g*h_right);
@@ -48,8 +102,8 @@
 		// Flux formulas
 		flux_left[0] = u_left*h_left;
-		flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
+		flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
 
 		flux_right[0] = u_right*h_right;
-		flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
+		flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
 
 		// Flux computation
@@ -60,16 +114,16 @@
 		} else {
 			edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
-			edgeflux[0] += s_max*s_min*(q_right[0]-q_left[0]);
+			edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
 			edgeflux[0] /= denom;
 			edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
-			edgeflux[1] += s_max*s_min*(q_right[1]-q_left[1]);
+			edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
 			edgeflux[1] /= denom;
 			edgeflux[1] *= normals;
     
-			// Maximal wavespeed
-			*max_speed = max(fabs(s_max), fabs(s_min));
+		// Maximal wavespeed
+        *max_speed = max(fabs(s_max), fabs(s_min));
 		}  
-		return 0;		
-}
+    return 0;		
+	}
 		
 		
@@ -90,11 +144,6 @@
 			   double* xmom_edge_values,
 			   double* bed_edge_values,
-			   double* height_edge_values,
-			   double* velocity_edge_values,
 			   double* stage_boundary_values,
 			   double* xmom_boundary_values,
-			   double* bed_boundary_values,
-			   double* height_boundary_values,
-			   double* velocity_boundary_values,
 			   double* stage_explicit_update,
 			   double* xmom_explicit_update,
@@ -102,6 +151,6 @@
 			   double* max_speed_array) {
 		
-                double flux[2], ql[5], qr[5], edgeflux[2];
-		double max_speed, normal;
+                double flux[2], ql[2], qr[2], edgeflux[2];
+		double zl, zr, max_speed, normal;
 		int k, i, ki, n, m, nm=0;
 		
@@ -116,7 +165,5 @@
 				ql[0] = stage_edge_values[ki];
 				ql[1] = xmom_edge_values[ki];
-				ql[2] = bed_edge_values[ki];
-				ql[3] = height_edge_values[ki];
-				ql[4] = velocity_edge_values[ki];
+				zl = bed_edge_values[ki];
 				
 				n = neighbours[ki];
@@ -125,7 +172,5 @@
 					qr[0] = stage_boundary_values[m];
 					qr[1] = xmom_boundary_values[m];
-					qr[2] = bed_edge_values[m];
-					qr[3] = height_edge_values[m];
-					qr[4] = velocity_edge_values[m];
+					zr = zl;
 				} else {
 					m = neighbour_vertices[ki];
@@ -133,11 +178,9 @@
 					qr[0] = stage_edge_values[nm];
 					qr[1] = xmom_edge_values[nm];
-					qr[2] = bed_edge_values[nm];
-					qr[3] = height_edge_values[nm];
-					qr[4] = velocity_edge_values[nm];
+					zr = bed_edge_values[nm];				
 				}
 				
 				normal = normals[ki];
-				_flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
+				_flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
 				flux[0] -= edgeflux[0];
 				flux[1] -= edgeflux[1];
@@ -152,5 +195,5 @@
 				    }
 				}
-			} // End edge i (and neighbour n)
+            } // End edge i (and neighbour n)
 			flux[0] /= areas[k];
 			stage_explicit_update[k] = flux[0];
@@ -163,4 +206,11 @@
 		return timestep;	
 	}
+	
+	
+	
+	
+	
+
+
 
 //=========================================================================
@@ -168,45 +218,107 @@
 //=========================================================================
 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-
-    PyObject
+  
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_vertices,
+	*normals, 
+	*areas,
+	*stage_edge_values,
+	*xmom_edge_values,
+	*bed_edge_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+    
+  double timestep, epsilon, g, h0, cfl;
+  int number_of_elements;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddddOOOOOOOOOOOiO",
+			&cfl,
+			&timestep,
+			&epsilon,
+			&g,
+			&h0,
+			&neighbours,
+			&neighbour_vertices,
+			&normals,
+			&areas,
+			&stage_edge_values,
+			&xmom_edge_values,
+			&bed_edge_values,
+			&stage_boundary_values,
+			&xmom_boundary_values,
+			&stage_explicit_update,
+			&xmom_explicit_update,
+			&number_of_elements,
+			&max_speed_array)) {
+    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
+    return NULL;
+  }
+
+ 
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_ext(cfl,
+				 timestep,
+				 epsilon,
+				 g,
+				 h0,
+				 (long*) neighbours -> data,
+				 (long*) neighbour_vertices -> data,
+				 (double*) normals -> data,
+				 (double*) areas -> data,
+				 (double*) stage_edge_values -> data,
+				 (double*) xmom_edge_values -> data,
+				 (double*) bed_edge_values -> data,
+				 (double*) stage_boundary_values -> data,
+				 (double*) xmom_boundary_values -> data,
+				 (double*) stage_explicit_update -> data,
+				 (double*) xmom_explicit_update -> data,
+				 number_of_elements,
+				 (double*) max_speed_array -> data);
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+//-----------------------
+PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
+  
+   PyObject 
 	*domain,
-	*stage,
-	*xmom,
-	*bed,
-	*height,
-	*velocity;
-
-    PyArrayObject
-	*neighbours,
+	*stage, 
+	*xmom, 
+	*bed;
+
+    PyArrayObject 
+	*neighbours, 
 	*neighbour_vertices,
-	*normals,
+	*normals, 
 	*areas,
 	*stage_vertex_values,
 	*xmom_vertex_values,
 	*bed_vertex_values,
-	*height_vertex_values,
-	*velocity_vertex_values,
 	*stage_boundary_values,
 	*xmom_boundary_values,
-	*bed_boundary_values,
-	*height_boundary_values,
-	*velocity_boundary_values,
 	*stage_explicit_update,
 	*xmom_explicit_update,
 	*max_speed_array;
-
+    
   double timestep, epsilon, g, h0, cfl;
   int number_of_elements;
 
-
+    
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dOOOOOO",
+  if (!PyArg_ParseTuple(args, "dOOOO",
 			&timestep,
 			&domain,
 			&stage,
 			&xmom,
-			&bed,
-			&height,
-			&velocity)) {
+			&bed)) {
       PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
       return NULL;
@@ -217,33 +329,32 @@
     g                 = get_python_double(domain,"g");
     h0                = get_python_double(domain,"h0");
-    cfl               = get_python_double(domain,"cfl");
-
-
+    cfl               = get_python_double(domain,"CFL");
+ 
+    
     neighbours        = get_consecutive_array(domain, "neighbours");
-    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
     normals           = get_consecutive_array(domain, "normals");
-    areas             = get_consecutive_array(domain, "areas");
+    areas             = get_consecutive_array(domain, "areas");    
     max_speed_array   = get_consecutive_array(domain, "max_speed_array");
-
-    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");
-    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");
-    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
-    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
-    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
-
-    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");
-    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");
-    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
-    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
-    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
-
-
-    stage_explicit_update = get_consecutive_array(stage, "explicit_update");
-    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");
+    
+    stage_vertex_values = get_consecutive_array(stage, "vertex_values");    
+    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");    
+    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+
+
 
     number_of_elements = stage_vertex_values -> dimensions[0];
 
-    // Call underlying flux computation routine and update
-    // the explicit update arrays
+
+ 
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
     timestep = _compute_fluxes_ext(
 	cfl,
@@ -259,11 +370,6 @@
 	(double*) xmom_vertex_values -> data,
 	(double*) bed_vertex_values -> data,
-	(double*) height_vertex_values -> data,
-	(double*) velocity_vertex_values -> data,
 	(double*) stage_boundary_values -> data,
 	(double*) xmom_boundary_values -> data,
-	(double*) bed_boundary_values -> data,
-	(double*) height_boundary_values -> data,
-	(double*) velocity_boundary_values -> data,
 	(double*) stage_explicit_update -> data,
 	(double*) xmom_explicit_update -> data,
@@ -279,11 +385,6 @@
   Py_DECREF(xmom_vertex_values);
   Py_DECREF(bed_vertex_values);
-  Py_DECREF(height_vertex_values);
-  Py_DECREF(velocity_vertex_values);
   Py_DECREF(stage_boundary_values);
   Py_DECREF(xmom_boundary_values);
-  Py_DECREF(bed_boundary_values);
-  Py_DECREF(height_boundary_values);
-  Py_DECREF(velocity_boundary_values);
   Py_DECREF(stage_explicit_update);
   Py_DECREF(xmom_explicit_update);
@@ -291,8 +392,11 @@
 
 
+
+
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
 
+ 
 
 
@@ -303,4 +407,5 @@
 static struct PyMethodDef MethodTable[] = {
   {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
   {NULL, NULL}
 };
Index: /anuga_work/development/pipeflow/sww_domain.py
===================================================================
--- /anuga_work/development/pipeflow/sww_domain.py	(revision 7822)
+++ /anuga_work/development/pipeflow/sww_domain.py	(revision 7823)
@@ -50,5 +50,5 @@
 
 #Shallow water domain
-class Domain(Generic_Domain):
+class Domain(Generic_domain):
 
     def __init__(self, coordinates, boundary = None, tagged_elements = None):
@@ -57,5 +57,5 @@
         evolved_quantities   = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
         other_quantities     = ['friction']
-        Generic_Domain.__init__(self,
+        Generic_domain.__init__(self,
                                 coordinates          = coordinates,
                                 boundary             = boundary,
@@ -143,5 +143,5 @@
 
 
-        timestep = numpy.float(sys.maxint)
+        timestep = float(sys.maxint)
 
         stage    = self.quantities['stage']
@@ -152,7 +152,11 @@
 
 
-        from sww_comp_flux_ext import compute_fluxes_ext
-
-        self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
+        from sww_comp_flux_ext import compute_fluxes_ext_short
+
+        #self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
+
+    
+ 	self.flux_timestep = compute_fluxes_ext_short(timestep,self,stage,xmom,bed)
+
 
 
@@ -168,5 +172,5 @@
 
         #Call basic machinery from parent class
-        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
+        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
                                        skip_initial_step):
 
@@ -266,5 +270,5 @@
 ##	else:
 
-    u_C[:]  = uh_C/(h_C + h0/h_C)
+    u_C[:,]  = uh_C/(h_C + h0/h_C)
 	
     for name in [ 'velocity', 'stage' ]:
@@ -286,5 +290,5 @@
     xmom_V  = domain.quantities['xmomentum'].vertex_values
 
-    h_V[:]    = stage_V - bed_V
+    h_V[:,:]    = stage_V - bed_V
     for i in range(len(h_C)):
         for j in range(2):
@@ -297,5 +301,5 @@
                 stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
                 
-    xmom_V[:] = u_V * h_V
+    xmom_V[:,:] = u_V * h_V
     
     return
@@ -314,6 +318,6 @@
 
     #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
+    wc    = domain.quantities['stage'].centroid_values
+    zc    = domain.quantities['elevation'].centroid_values
     xmomc = domain.quantities['xmomentum'].centroid_values
     hc = wc - zc  #Water depths at centroids
Index: /anuga_work/development/pipeflow/sww_python.py
===================================================================
--- /anuga_work/development/pipeflow/sww_python.py	(revision 7822)
+++ /anuga_work/development/pipeflow/sww_python.py	(revision 7823)
@@ -58,5 +58,5 @@
     for k in range(N):
 
-        flux[:] = 0.  #Reset work numpy.array
+        flux[:,] = 0.  #Reset work numpy.array
         #for i in range(3):
         for i in range(2):
Index: /anuga_work/development/pipeflow/test_quantity.py
===================================================================
--- /anuga_work/development/pipeflow/test_quantity.py	(revision 7822)
+++ /anuga_work/development/pipeflow/test_quantity.py	(revision 7823)
@@ -1,2 +1,3 @@
+import quantity
 #!/usr/bin/env python
 
@@ -5,18 +6,36 @@
 
 
-from shallow_water_domain import *
-from shallow_water_domain import flux_function as domain_flux_function
-
-
-
+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_Shallow_Water(unittest.TestCase):
+class Test_Quantity(unittest.TestCase):
     def setUp(self):
-        self.points = [0.0, 1.0, 2.0, 3.0]
-        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
-        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
-        
+        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
@@ -24,424 +43,563 @@
 
 
+    def test_creat_with_boundary(self):
+
+        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
+        
     def test_creation(self):
-        domain = Domain(self.points)
-        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])
-
-    def test_compute_fluxes(self):
-        """
-        Compare shallow_water_domain flux calculation against a previous
-        Python implementation (defined in this file)
-        """
-        domain = Domain(self.points)
-        domain.set_quantity('stage',2.0)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        stage_ud, xmom_ud = local_compute_fluxes(domain)
-
-        domain.distribute_to_vertices_and_edges()
-        domain.compute_fluxes()
-
-        print domain.quantities['stage'].explicit_update
-        print domain.quantities['xmomentum'].explicit_update
-        print stage_ud
-        print xmom_ud
-
-        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
-        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
-
-
-    def test_local_flux_function(self):
-        normal = 1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        #This assumes h0 = 1.0e-3!!
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-        #print maxspeed
-        #print edgeflux
-        
-        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
-        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-        normal = -1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-
-
-        #print maxspeed
-        #print edgeflux        
-        
-        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
-        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-    def test_domain_flux_function(self):
-        normal = 1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-
-        #print edgeflux
-
-
-
-        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
-
-        #print domainedgeflux
-
-        assert numpy.allclose(domainedgeflux, edgeflux)
-        assert numpy.allclose(domainmaxspeed, maxspeed)
-
-    def test_gravity(self):
-        """
-        Compare shallow_water_domain gravity calculation
-        """
-
-        def slope_one(x):
+
+        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
 
-        domain = Domain(self.points)
-        domain.set_quantity('stage',4.0)
-        domain.set_quantity('elevation',slope_one)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        gravity(domain)
-        
-        #print domain.quantities['stage'].vertex_values 
-        #print domain.quantities['elevation'].vertex_values 
-        #print domain.quantities['xmomentum'].explicit_update       
-
-        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
-
-
-    def xtest_evolve_first_order(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-        
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 1
-        domain.set_timestepping_method('euler')
-        
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            #pass
-        
-            print domain.quantities['stage'].vertex_values
-            print domain.quantities['elevation'].vertex_values
-            print domain.quantities['xmomentum'].vertex_values
-          
-        
-            print domain.quantities['stage'].centroid_values  
-            print domain.quantities['elevation'].centroid_values
-            print domain.quantities['xmomentum'].centroid_values    
-
-        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-    def xtest_evolve_euler_second_order_space(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 2
-        domain.set_timestepping_method('rk2')
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            pass
-        
-        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-    def xtest_evolve_second_order_space_time(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 2
-        domain.set_timestepping_method('rk3')
-
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            pass
-        
-        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-
-#==============================================================================
-
-def local_compute_fluxes(domain):
-    """Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    """
-
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-    tmp0 = zeros(N,Float)
-    tmp1 = zeros(N,Float)
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    #    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    #Arrays
-    #stage = Stage.edge_values
-    #xmom =  Xmom.edge_values
-    #   ymom =  Ymom.edge_values
-    #bed =   Bed.edge_values
-    
-    stage = Stage.vertex_values
-    xmom =  Xmom.vertex_values
-    bed =   Bed.vertex_values
-    
-    #print 'stage edge values', stage
-    #print 'xmom edge values', xmom
-    #print 'bed values', bed
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-    #print 'stage_bdry',stage_bdry
-    #print 'xmom_bdry', xmom_bdry
-    #    ymom_bdry =  Ymom.boundary_values
-
-    #    flux = zeros(3, Float) #Work numpy.array for summing up fluxes
-    flux = zeros(2, Float) #Work numpy.array for summing up fluxes
-    ql = zeros(2, Float)
-    qr = zeros(2, Float)
-
-    #Loop
-    timestep = float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work numpy.array
-        #for i in range(3):
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            #ql[0] = stage[k, i]
-            #ql[1] = xmom[k, i]
-            ql = [stage[k, i], xmom[k, i]]
-            zl = bed[k, i]
-
-            #Quantities at neighbour on nearest face
-            n = domain.neighbours[k,i]
-            if n < 0:
-                m = -n-1 #Convert negative flag to index
-                qr[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                #print i, ' ' , m
-                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                qr[0] = stage[n, m]
-                qr[1] = xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-        
-            #Flux computation using provided function
-            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            #print 'ql',ql
-            #print 'qr',qr
-            
-
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-            #print 'edgeflux', edgeflux
-
-            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
-            # flux = edgefluxleft - edgefluxright
-            flux -= edgeflux #* domain.edgelengths[k,i]
-            #Update optimal_timestep
-            try:
-                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
-            except ZeroDivisionError:
-                pass
-
-        #Normalise by area and store for when all conserved
-        #quantities get updated
-        flux /= domain.areas[k]
-
-        #Stage.explicit_update[k] = flux[0]
-        tmp0[k] = flux[0]
-        tmp1[k] = flux[1]
-
-
-    return tmp0, tmp1
-
-
-def local_flux_function(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numerical vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from config import g, epsilon, h0
-    from math import sqrt
-
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #Take average of field values
-
-    w_left  = q_left[0]   #w=h+z
-    h_left  = w_left-z
-    uh_left = q_left[1]
-
-    if h_left < epsilon:
-        u_left = 0.0  #Could have been negative
-        h_left = 0.0
-    else:
-        u_left  = uh_left/(h_left +  h0/h_left)
-
-
-    uh_left = u_left*h_left
-
-    w_right  = q_right[0]  #w=h+z
-    h_right  = w_right-z
-    uh_right = q_right[1]
-
-
-    if h_right < epsilon:
-        u_right = 0.0  #Could have been negative
-        h_right = 0.0
-    else:
-        u_right  = uh_right/(h_right + h0/h_right)
-
-    uh_right = u_right*h_right
-    
-    #vh_left  = q_left[2]
-    #vh_right = q_right[2]
-
-    #print h_right
-    #print u_right
-    #print h_left
-    #print u_right
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    
-    #Flux computation
-
-    #flux_left  = numpy.array([u_left*h_left,
-    #                    u_left*uh_left + 0.5*g*h_left**2])
-    #flux_right = numpy.array([u_right*h_right,
-    #                    u_right*uh_right + 0.5*g*h_right**2])
-    flux_left  = numpy.array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left*h_left])
-    flux_right = numpy.array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right*h_right])
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = numpy.array([0.0, 0.0])
-        max_speed = 0.0
-    else:
-        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
-        edgeflux += s_max*s_min*(q_right-q_left)/denom
-        
-        edgeflux[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
+        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
+        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.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_Shallow_Water, 'test')
-    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
-
-
+    suite = unittest.makeSuite(Test_Quantity, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /anuga_work/development/pipeflow/test_sww_domain.py
===================================================================
--- /anuga_work/development/pipeflow/test_sww_domain.py	(revision 7822)
+++ /anuga_work/development/pipeflow/test_sww_domain.py	(revision 7823)
@@ -43,4 +43,5 @@
         print domain.quantities['stage'].explicit_update
         print domain.quantities['xmomentum'].explicit_update
+
         print stage_ud
         print xmom_ud
@@ -50,5 +51,5 @@
 
 
-    def test_local_flux_function(self):
+    def test_sww_python_flux_function(self):
         normal = 1.0
         ql = numpy.array([1.0, 2.0],numpy.float)
@@ -58,5 +59,5 @@
 
         #This assumes h0 = 1.0e-3!!
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
         #print maxspeed
         #print edgeflux
@@ -71,5 +72,5 @@
         zr = 0.0
 
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
 
 
@@ -104,5 +105,5 @@
 
 
-    def xtest_evolve_first_order(self):
+    def test_evolve_first_order(self):
         """
         Compare still lake solution for various versions of shallow_water_domain
@@ -139,5 +140,5 @@
 
 
-    def xtest_evolve_euler_second_order_space(self):
+    def test_evolve_euler_second_order_space(self):
         """
         Compare still lake solution for various versions of shallow_water_domain
@@ -163,5 +164,5 @@
         assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
 
-    def xtest_evolve_second_order_space_time(self):
+    def test_evolve_second_order_space_time(self):
         """
         Compare still lake solution for various versions of shallow_water_domain
@@ -192,232 +193,4 @@
 #==============================================================================
 
-def local_compute_fluxes(domain):
-    """Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    """
-
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-    tmp0 = zeros(N,Float)
-    tmp1 = zeros(N,Float)
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    #    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    #Arrays
-    #stage = Stage.edge_values
-    #xmom =  Xmom.edge_values
-    #   ymom =  Ymom.edge_values
-    #bed =   Bed.edge_values
-    
-    stage = Stage.vertex_values
-    xmom =  Xmom.vertex_values
-    bed =   Bed.vertex_values
-    
-    #print 'stage edge values', stage
-    #print 'xmom edge values', xmom
-    #print 'bed values', bed
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-    #print 'stage_bdry',stage_bdry
-    #print 'xmom_bdry', xmom_bdry
-    #    ymom_bdry =  Ymom.boundary_values
-
-    #    flux = zeros(3, Float) #Work numpy.array for summing up fluxes
-    flux = zeros(2, Float) #Work numpy.array for summing up fluxes
-    ql = zeros(2, Float)
-    qr = zeros(2, Float)
-
-    #Loop
-    timestep = float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work numpy.array
-        #for i in range(3):
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            #ql[0] = stage[k, i]
-            #ql[1] = xmom[k, i]
-            ql = [stage[k, i], xmom[k, i]]
-            zl = bed[k, i]
-
-            #Quantities at neighbour on nearest face
-            n = domain.neighbours[k,i]
-            if n < 0:
-                m = -n-1 #Convert negative flag to index
-                qr[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                #print i, ' ' , m
-                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                qr[0] = stage[n, m]
-                qr[1] = xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-        
-            #Flux computation using provided function
-            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            #print 'ql',ql
-            #print 'qr',qr
-            
-
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-            #print 'edgeflux', edgeflux
-
-            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
-            # flux = edgefluxleft - edgefluxright
-            flux -= edgeflux #* domain.edgelengths[k,i]
-            #Update optimal_timestep
-            try:
-                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
-            except ZeroDivisionError:
-                pass
-
-        #Normalise by area and store for when all conserved
-        #quantities get updated
-        flux /= domain.areas[k]
-
-        #Stage.explicit_update[k] = flux[0]
-        tmp0[k] = flux[0]
-        tmp1[k] = flux[1]
-
-
-    return tmp0, tmp1
-
-
-def local_flux_function(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numerical vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from config import g, epsilon, h0
-    from math import sqrt
-
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #Take average of field values
-
-    w_left  = q_left[0]   #w=h+z
-    h_left  = w_left-z
-    uh_left = q_left[1]
-
-    if h_left < epsilon:
-        u_left = 0.0  #Could have been negative
-        h_left = 0.0
-    else:
-        u_left  = uh_left/(h_left +  h0/h_left)
-
-
-    uh_left = u_left*h_left
-
-    w_right  = q_right[0]  #w=h+z
-    h_right  = w_right-z
-    uh_right = q_right[1]
-
-
-    if h_right < epsilon:
-        u_right = 0.0  #Could have been negative
-        h_right = 0.0
-    else:
-        u_right  = uh_right/(h_right + h0/h_right)
-
-    uh_right = u_right*h_right
-    
-    #vh_left  = q_left[2]
-    #vh_right = q_right[2]
-
-    #print h_right
-    #print u_right
-    #print h_left
-    #print u_right
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    
-    #Flux computation
-
-    #flux_left  = numpy.array([u_left*h_left,
-    #                    u_left*uh_left + 0.5*g*h_left**2])
-    #flux_right = numpy.array([u_right*h_right,
-    #                    u_right*uh_right + 0.5*g*h_right**2])
-    flux_left  = numpy.array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left*h_left])
-    flux_right = numpy.array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right*h_right])
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = numpy.array([0.0, 0.0])
-        max_speed = 0.0
-    else:
-        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
-        edgeflux += s_max*s_min*(q_right-q_left)/denom
-        
-        edgeflux[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
-#-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Shallow_Water, 'test')
Index: /anuga_work/development/pipeflow/test_sww_domain_vel.py
===================================================================
--- /anuga_work/development/pipeflow/test_sww_domain_vel.py	(revision 7822)
+++ /anuga_work/development/pipeflow/test_sww_domain_vel.py	(revision 7823)
@@ -8,5 +8,5 @@
 
 
-from numpy import allclose, array, ones, numpy.float, maximum, zeros
+from Numeric import allclose, array, ones, Float, maximum, zeros
 
 
@@ -114,18 +114,19 @@
         domain.compute_fluxes()
         
-        print domain.quantities['stage'].vertex_values
-        print domain.quantities['xmomentum'].vertex_values
-        print domain.quantities['elevation'].vertex_values
-        print domain.quantities['height'].vertex_values
-        print domain.quantities['velocity'].vertex_values
-
-        print domain.quantities['stage'].boundary_values
-        print domain.quantities['xmomentum'].boundary_values
-        print domain.quantities['elevation'].boundary_values
-        print domain.quantities['height'].boundary_values
-        print domain.quantities['velocity'].boundary_values
+#        print domain.quantities['stage'].vertex_values
+#        print domain.quantities['xmomentum'].vertex_values
+#        print domain.quantities['elevation'].vertex_values
+#        print domain.quantities['height'].vertex_values
+#        print domain.quantities['velocity'].vertex_values
+#
+#        print domain.quantities['stage'].boundary_values
+#        print domain.quantities['xmomentum'].boundary_values
+#        print domain.quantities['elevation'].boundary_values
+#        print domain.quantities['height'].boundary_values
+#        print domain.quantities['velocity'].boundary_values
         
         print domain.quantities['stage'].explicit_update
         print domain.quantities['xmomentum'].explicit_update
+
         print stage_ud
         print xmom_ud
@@ -137,6 +138,6 @@
     def test_local_flux_function(self):
         normal = 1.0
-        ql = array([1.0, 2.0],numpy.float)
-        qr = array([1.0, 2.0],numpy.float)
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
         zl = 0.0
         zr = 0.0
@@ -147,10 +148,10 @@
         #print edgeflux
         
-        assert allclose(array([2.0, 8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(array([2.0, 8.9],Float), edgeflux, rtol=1.0e-005)
         assert allclose(5.1305, maxspeed, rtol=1.0e-005)
 
         normal = -1.0
-        ql = array([1.0, 2.0],numpy.float)
-        qr = array([1.0, 2.0],numpy.float)
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
         zl = 0.0
         zr = 0.0
@@ -162,5 +163,5 @@
         #print edgeflux        
         
-        assert allclose(array([-2.0, -8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(array([-2.0, -8.9],Float), edgeflux, rtol=1.0e-005)
         assert allclose(5.1305, maxspeed, rtol=1.0e-005)
 
@@ -188,5 +189,5 @@
         #print domain.quantities['xmomentum'].explicit_update       
 
-        assert allclose( array([-34.3, -24.5, -14.7], numpy.float), domain.quantities['xmomentum'].explicit_update )
+        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
 
 
@@ -300,10 +301,10 @@
 
     import sys
-    from numpy import zeros, numpy.float
+    from Numeric import zeros, Float
 
     N = domain.number_of_elements
 
-    tmp0 = zeros(N,numpy.float)
-    tmp1 = zeros(N,numpy.float)
+    tmp0 = zeros(N,Float)
+    tmp1 = zeros(N,Float)
 
     #Shortcuts
@@ -333,15 +334,15 @@
     #    ymom_bdry =  Ymom.boundary_values
 
-    #    flux = zeros(3, numpy.float) #Work array for summing up fluxes
-    flux = zeros(2, numpy.float) #Work array for summing up fluxes
-    ql = zeros(2, numpy.float)
-    qr = zeros(2, numpy.float)
+    #    flux = zeros(3, Float) #Work array for summing up fluxes
+    flux = zeros(2, Float) #Work array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
 
     #Loop
-    timestep = numpy.float(sys.maxint)
+    timestep = float(sys.maxint)
     enter = True
     for k in range(N):
 
-        flux[:] = 0.  #Reset work array
+        flux[:,] = 0.  #Reset work array
         #for i in range(3):
         for i in range(2):
@@ -415,5 +416,5 @@
 
     Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numpyal vectors ql an qr.
+    in the numerical vectors ql an qr.
 
     Bed elevations zl and zr.
@@ -422,5 +423,5 @@
     from config import g, epsilon, h0
     from math import sqrt
-    from numpy import array
+    from Numeric import array
 
     #print 'ql',ql
