Index: inundation/ga/storm_surge/examples/advection_dirichlet.py
===================================================================
--- inundation/ga/storm_surge/examples/advection_dirichlet.py	(revision 195)
+++ inundation/ga/storm_surge/examples/advection_dirichlet.py	(revision 195)
@@ -0,0 +1,37 @@
+"""Simple example using the module for solving the advection equation
+
+Initial condition is zero, boudary conditions are dirichlet at the left edge
+and transmissive everywhere else
+"""
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+from mesh_factory import rectangular
+from advection import Domain, Transmissive_boundary, Dirichlet_boundary
+
+from Numeric import array
+
+#Create basic mesh
+points, vertices, boundary = rectangular(8, 8)
+
+#Create advection domain with direction (1,-1)
+# - Initial condition is zero by default
+domain = Domain(points, vertices, boundary, velocity=[1.0, -1.0])
+domain.smooth = False
+domain.visualise = True
+
+#Boundaries
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary(array([0.1]))
+
+domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
+domain.check_integrity()
+
+###################
+# Evolution
+for t in domain.evolve(yieldstep = 0.02, finaltime = 2):
+    domain.write_time()
+
+
+
Index: inundation/ga/storm_surge/examples/advection_waves.py
===================================================================
--- inundation/ga/storm_surge/examples/advection_waves.py	(revision 195)
+++ inundation/ga/storm_surge/examples/advection_waves.py	(revision 195)
@@ -0,0 +1,49 @@
+"""Example using the module for solving the advection equation
+
+
+Initial condition is zero, boundary conditions are transmissive everywhere except
+for three segments where three different time dependent conditions are applied.
+"""
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+from mesh_factory import rectangular
+from advection import Domain, Transmissive_boundary, Dirichlet_boundary, Time_boundary
+
+#Create basic mesh
+points, vertices, boundary = rectangular(20, 4)
+
+#Create advection domain
+# - Initial condition is zero by default
+domain = Domain(points, vertices, boundary, velocity=[1.0, -0.8])
+#domain.smooth = False
+domain.visualise = True
+
+#Boundaries
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([0.1])
+
+from math import sin, pi
+T1 = Time_boundary(domain, f=lambda x: [0.2*(sin(2*x*pi)+1)/2])
+T2 = Time_boundary(domain, f=lambda x: [0.2*(sin(4*x*pi)+1)/3])
+T3 = Time_boundary(domain, f=lambda x: [0.2*(sin(6*x*pi)+1)/4])
+
+#Modify boundary tags
+domain.boundary[(7, 1)] = 't1'
+domain.boundary[(5, 2)] = 't2'
+domain.boundary[(3, 2)] = 't3'
+
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T,
+                      't1': T1, 't2': T2, 't3': T3} )
+domain.check_integrity()
+
+
+###################
+# Evolution
+for t in domain.evolve(yieldstep = 0.02, finaltime = None):
+    domain.write_time()
+
+
+
+
Index: inundation/ga/storm_surge/pyvolution/advection.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/advection.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/advection.py	(revision 195)
@@ -0,0 +1,161 @@
+"""Class Domain -
+2D triangular domains for finite-volume computations of
+the advection equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the advection equantion
+
+The equation is
+
+  u_t + (v_1 u)_x + (v_2 u)_y = 0
+
+There is only one conserved quantity, the level u
+
+The advection equation is a very simple specialisation of the generic
+domain and may serve as an instructive example or a test of other
+components such as visualisation.
+
+Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2004   
+"""
+
+from domain import *
+Generic_domain = Domain #Rename
+
+class Domain(Generic_domain):
+
+    def __init__(self, coordinates, vertices, boundary = None, velocity = None):
+
+        
+        Generic_domain.__init__(self, coordinates, vertices, boundary,
+                                ['level'])
+
+        if velocity is None:
+            self.velocity = [1,0]
+        else:
+            self.velocity = velocity
+
+        #Only first is implemented for advection    
+        self.default_order = self.order = 1  
+
+
+    def check_integrity(self):
+        Generic_domain.check_integrity(self)
+
+        msg = 'Conserved quantity must be "level"'
+        assert self.conserved_quantities[0] == 'level', msg
+        
+
+    def flux_function(self, normal, ql, qr, zl=None, zr=None):
+        """Compute outward flux as inner product between velocity
+        vector v=(v_1, v_2) and normal vector n.
+            
+        if <n,v> > 0 flux direction is outward bound and its magnitude is
+        determined by the quantity inside volume: ql.
+        Otherwise it is inbound and magnitude is determined by the
+        quantity outside the volume: qr.
+        """
+    
+        v1 = self.velocity[0]
+        v2 = self.velocity[1]
+
+
+        normal_velocity = v1*normal[0] + v2*normal[1]
+
+        if normal_velocity < 0:
+            flux = qr * normal_velocity
+        else:
+            flux = ql * normal_velocity
+        
+        max_speed = abs(normal_velocity)
+        return flux, max_speed        
+
+        
+    def compute_fluxes(self):
+        """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
+        domain.explicit_update
+
+        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
+        from config import max_timestep
+
+        N = self.number_of_elements
+    
+        neighbours = self.neighbours
+        neighbour_edges = self.neighbour_edges
+        normals = self.normals
+
+        areas = self.areas
+        radii = self.radii
+        edgelengths = self.edgelengths
+            
+        timestep = max_timestep #FIXME: Get rid of this
+
+        #Shortcuts
+        Level = self.quantities['level']
+
+        #Arrays
+        level = Level.edge_values
+
+        level_bdry = Level.boundary_values
+    
+        flux = zeros(1, Float) #Work array for summing up fluxes
+
+        #Loop
+        for k in range(N):
+            optimal_timestep = float(sys.maxint)
+
+            flux[:] = 0.  #Reset work array
+            for i in range(3):
+                #Quantities inside volume facing neighbour i
+                ql = level[k, i]
+
+                #Quantities at neighbour on nearest face
+                n = neighbours[k,i] 
+                if n < 0:
+                    m = -n-1 #Convert neg flag to index
+                    qr = level_bdry[m]
+                else:    
+                    m = neighbour_edges[k,i]
+                    qr = level[n, m]
+
+                    
+                #Outward pointing normal vector   
+                normal = normals[k, 2*i:2*i+2]
+
+                #Flux computation using provided function
+                edgeflux, max_speed = self.flux_function(normal, ql, qr)
+                flux -= edgeflux * edgelengths[k,i]
+                
+                #Update optimal_timestep
+                try:
+                    optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
+                except ZeroDivisionError:
+                    pass
+
+            #Normalise by area and store for when all conserved
+            #quantities get updated
+            flux /= areas[k]
+            Level.explicit_update[k] = flux[0]
+            
+            timestep = min(timestep, optimal_timestep)
+
+        self.timestep = timestep    
+	
+
Index: inundation/ga/storm_surge/pyvolution/config.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/config.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/config.py	(revision 195)
@@ -0,0 +1,68 @@
+"""Module where global model parameters are set
+"""
+
+epsilon = 1.0e-13
+
+default_boundary_tag = 'exterior'
+
+
+min_timestep = 1.0e-6 #Should be computed based on geometry
+max_timestep = 1000
+#This is how:
+#Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)
+#Then work out minimal internal distance in mesh r_min and set
+#min_timestep = r_min/v_max 
+# 
+#Max speeds are calculated in the flux function as
+#
+#lambda = v +/- sqrt(gh)
+#
+# so with 500 m/s, h ~ 500^2/g = 2500 m well out of the domain of the
+# shallow water wave equation
+#
+#The actual soundspeed can be as high as 1530m
+#(see http://staff.washington.edu/aganse/public.projects/clustering/clustering.html), but that would only happen with h>225000m in this equation. Why ?
+#The maximal speed we specify is really related to the max speed
+#of surface pertubation
+
+
+
+v_max = 100 #For use in domain_ext.c
+sound_speed = 500
+
+
+max_smallsteps = 10  #Max number of degenerate steps allowed b4 trying first order
+
+g = 9.8        #Gravity
+manning = 0.3 #Manning's friction coefficient
+
+
+
+beta = 0.9  #]0;1] (e.g. 0.9). 1 allows the steepest gradients,
+            #lower values are the more conservative
+            #the limiter being a first order method.
+
+pmesh_filename = '.\\pmesh'
+
+
+import os, sys
+
+if sys.platform == 'win32':
+    data_dir = 'C:\grohm_output'
+else:    
+    data_dir = os.path.expanduser('~'+os.sep+'grohm_output')
+
+
+
+
+use_extensions = True   #Try to use C-extensions
+#use_extensions = False   #Do not use C-extensions
+
+use_psyco = True  #Use psyco optimisations
+#use_psyco = False  #Do not use psyco optimisations
+
+
+#Specific to shallow water W.E.
+minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0
+number_of_conserved_quantities = 3 #Should really be set dynamically - but hey!
+number_of_field_values = 2
Index: inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/domain.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/domain.py	(revision 195)
@@ -0,0 +1,483 @@
+"""Class Domain - 2D triangular domains for finite-volume computations of
+   the shallow water wave equation
+
+
+   Copyright 2004
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia
+"""
+
+from mesh import Mesh
+from generic_boundary_conditions import *
+
+class Domain(Mesh):
+
+    def __init__(self, coordinates, vertices, boundary = None,
+                 conserved_quantities = None, other_quantities = None):
+
+        Mesh.__init__(self, coordinates, vertices, boundary)
+
+        from Numeric import zeros, Float 
+        from quantity import Quantity
+
+        #List of quantity names entering
+        #the conservation equations
+        #(Must be a subset of quantities)
+        if conserved_quantities is None:            
+            self.conserved_quantities = []
+        else:
+            self.conserved_quantities = conserved_quantities            
+
+        if other_quantities is None:
+            other_quantities = []
+
+        #Build dictionary of Quantity instances keyed by quantity names
+        self.quantities = {}
+        for name in self.conserved_quantities + other_quantities:
+            self.quantities[name] = Quantity(self)
+
+        #Defaults
+        from config import max_smallsteps, beta
+        self.beta = beta
+        self.default_order = 1
+        self.order = self.default_order
+        self.smallsteps = 0
+        self.max_smallsteps = max_smallsteps        
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0        
+
+        #Model time    
+        self.time = 0.0
+        self.finaltime = None
+        self.min_timestep = self.max_timestep = 0.0        
+	
+        #Checkpointing		
+	self.filename = 'domain'
+	self.checkpoint = False
+
+        #Realtime visualisation
+        self.visualise = False
+        
+        #Stored output
+        self.store=False
+	self.format = 'dat'    
+        #self.smooth = False
+        self.smooth = True
+
+        #Reduction operation for get_vertex_values              
+        #from pytools.stats import mean
+        #self.reduction = mean 
+        self.reduction = min  #Looks better near steep slopes
+        
+
+    #Public interface to Domain        
+    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
+        """Get conserved quantities at volume vol_id
+
+        If vertex is specified use it as index for vertex values
+        If edge is specified use it as index for edge values
+        If neither are specified use centroid values
+        If both are specified an exeception is raised
+
+        Return value: Vector of length == number_of_conserved quantities 
+        
+        """
+
+        from Numeric import zeros, Float
+        
+        if not (vertex is None or edge is None):
+            msg = 'Values for both vertex and edge was specified.'
+            msg += 'Only one (or none) is allowed.' 
+            raise msg
+
+        q = zeros( len(self.conserved_quantities), Float)
+        
+        for i, name in enumerate(self.conserved_quantities):
+            Q = self.quantities[name]
+            if vertex is not None:
+                q[i] = Q.vertex_values[vol_id, vertex]
+            elif edge is not None:
+                q[i] = Q.edge_values[vol_id, edge]
+            else:
+                q[i] = Q.centroid_values[vol_id]                
+
+        return q 
+
+    def set_quantity(self, name, values, location='vertices'):
+        """Set values for named quantity
+
+        name: Name of quantity
+        values: Compatible list or Numeric array (see below)
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+        """
+
+        self.quantities[name].set_values(values, location)
+        
+
+    def set_boundary(self, boundary_map):
+        """Associate boundary objects with tagged boundary segments.
+
+        Input boundary_map is a dictionary of boundary objects keyed
+        by symbolic tags to matched against tags in the internal dictionary
+        self.boundary.
+
+        As result one pointer to a boundary object is stored for each vertex
+        in the list self.boundary_objects.
+        More entries may point to the same boundary object
+        
+        Schematically the mapping is:
+
+        self.boundary_segments: k: (vol_id, edge_id)
+        self.boundary:          (vol_id, edge_id): tag
+        boundary_map (input):   tag: boundary_object
+        ----------------------------------------------
+        self.boundary_objects:  k: boundary_object
+        
+
+        Pre-condition:
+          self.boundary and self.boundary_segments have been built.
+
+        Post-condition:
+          self.boundary_objects is built
+
+        If a tag from the domain doesn't appear in the input dictionary an
+        exception is raised.
+        However, if a tag is not used to the domain, no error is thrown.
+        FIXME: This would lead to implementation of a
+        default boundary condition
+        """
+
+        self.boundary_objects = []
+        for k, (vol_id, edge_id) in enumerate(self.boundary_segments):
+            tag = self.boundary[ (vol_id, edge_id) ]
+
+            if boundary_map.has_key(tag):
+                B = boundary_map[tag]
+                self.boundary_objects.append(B)                
+
+            else:
+                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
+                msg += 'bound to a boundary object.\n'
+                msg += 'All boundary tags defined in domain must appear '
+                msg += 'in the supplied dictionary.\n'
+                msg += 'The tags are: %s' %self.get_boundary_tags()
+                raise msg
+
+
+
+    #MISC        
+    def check_integrity(self):
+        Mesh.check_integrity(self)
+
+        for quantity in self.conserved_quantities:
+            msg = 'Conserved quantities must be a subset of all quantities'
+            assert quantity in self.quantities, msg
+    
+    def write_time(self):
+        if self.min_timestep == self.max_timestep:
+            print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+                  %(self.time, self.min_timestep, self.number_of_steps,
+                    self.number_of_first_order_steps)
+        else:    
+            print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+                  %(self.time, self.min_timestep,
+                    self.max_timestep, self.number_of_steps,
+                    self.number_of_first_order_steps)
+
+
+
+
+    ###########################
+    #Main components of evolve
+
+    def evolve(self, yieldstep = None, finaltime = None):
+        """Evolve model from time=0.0 to finaltime yielding results
+        every yieldstep.
+        
+        Internally, smaller timesteps may be taken.
+
+        Evolve is implemented as a generator and is to be called as such, e.g.
+        
+        for t in domain.evolve(timestep, yieldstep, finaltime):
+            <Do something with domain and t>
+        
+        """
+
+        #import data_manager
+        from config import min_timestep, max_timestep, epsilon
+
+        if yieldstep is None:
+            yieldstep = max_timestep
+
+
+        self.order = self.default_order
+        
+            
+        self.yieldtime = 0.0 #Time between 'yields'
+
+        #Initialise interval of timestep sizes (for reporting only)
+        self.min_timestep = max_timestep
+        self.max_timestep = min_timestep
+        self.finaltime = finaltime
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0                
+
+        #Initial update of vertex and edge values
+        self.distribute_to_vertices_and_edges()
+	
+	
+        #Or maybe restore from latest checkpoint
+	if self.checkpoint is True:
+            self.goto_latest_checkpoint()
+            
+
+        #Store model data, e.g. for visualisation    
+        if self.store is True and self.time == 0.0:
+            self.store_bathymetry()        
+            self.store_conserved_quantities()	
+
+        if self.visualise is True and self.time == 0.0:
+            import realtime_visualisation as visualise
+            visualise.create_surface(self)
+        
+	
+        yield(self.time)  #Yield initial values
+	
+        while True:
+            #Update boundary values
+            self.update_boundary()            
+
+            #Compute all fluxes and timestep suitable for all volumes
+            self.compute_fluxes()
+
+            #Update timestep to fit yieldstep and finaltime
+            self.update_timestep(yieldstep, finaltime)
+
+            #Update conserved quantities
+            self.update_conserved_quantities()            
+
+            #Update vertex and edge values
+            self.distribute_to_vertices_and_edges()
+                
+            #Update time    
+            self.time += self.timestep
+            self.yieldtime += self.timestep
+            self.number_of_steps += 1
+            if self.order == 1:
+                self.number_of_first_order_steps += 1        
+
+            #Yield results
+            if finaltime is not None and abs(self.time - finaltime) < epsilon:
+                # Yield final time and allow inspection of domain
+                yield(self.time)
+                break
+                
+            if abs(self.yieldtime - yieldstep) < epsilon: 
+                # Yield (intermediate) time and allow inspection of domain
+
+                if self.checkpoint is True:
+                    self.store_checkpoint()
+                    self.delete_old_checkpoints()
+                    
+                #Store model data, e.g. for subsequent visualisation    
+                if self.store is True:
+                    self.store_conserved_quantities()
+
+                #Real time viz
+                if self.visualise is True:
+                    visualise.update(self)
+                    
+                yield(self.time)  
+
+                # Reinitialise
+                self.yieldtime = 0.0                
+                self.min_timestep = max_timestep
+                self.max_timestep = min_timestep
+                self.number_of_steps = 0
+                self.number_of_first_order_steps = 0                        
+
+            
+    def evolve_to_end(self, finaltime = 1.0):
+        """Iterate evolve generator all the way to the end
+        """
+
+        for _ in self.evolve(yieldstep=None, finaltime=finaltime):            
+            pass
+        
+
+    
+    def update_boundary(self):
+        """Go through list of boundary objects and update boundary values
+        for all conserved quantities on boundary.
+        """
+
+        #FIXME: Update only those that change (if that can be worked out)
+        for i, B in enumerate(self.boundary_objects):
+            vol_id, edge_id = self.boundary_segments[i]
+            q = B.evaluate(vol_id, edge_id)
+
+            for j, name in enumerate(self.conserved_quantities):
+                Q = self.quantities[name]
+                Q.boundary_values[i] = q[j]
+
+
+    def compute_fluxes(self):
+        msg = 'Method compute_fluxes must be overridden by Domain subclass'
+        raise msg
+
+
+    def update_timestep(self, yieldstep, finaltime):
+
+        from config import min_timestep
+    
+        timestep = self.timestep
+    
+        #Record maximal and minimal values of timestep for reporting    
+        self.max_timestep = max(timestep, self.max_timestep)
+        self.min_timestep = min(timestep, self.min_timestep)
+    
+        #Protect against degenerate time steps
+        if timestep < min_timestep:
+
+            #Number of consecutive small steps taken b4 taking action
+            self.smallsteps += 1 
+
+            if self.smallsteps > self.max_smallsteps:
+                self.smallsteps = 0 #Reset
+               
+                if self.order == 1:
+                    msg = 'Minimal timestep %.16f reached ' %self.min_timestep
+                    msg += 'using 1 order scheme'
+
+                    raise msg
+                else:
+                    #Try to overcome situation by switching to 1 order
+                    self.order = 1
+
+        else:
+            self.smallsteps = 0
+            if self.order == 1:
+                if self.order != self.default_order:
+                    self.order = 2
+                    
+
+        #Ensure that final time is not exceeded
+        if finaltime is not None and self.time + timestep > finaltime:
+            timestep = finaltime-self.time
+
+        #Ensure that model time is aligned with yieldsteps
+        if self.yieldtime + timestep > yieldstep:
+            timestep = yieldstep-self.yieldtime
+
+        self.timestep = timestep
+
+
+
+    def compute_forcing_terms(self):          
+        """If there are any forcing functions driving the system
+        they should be defined in Domain subclass
+        """
+        pass
+
+
+    def update_conserved_quantities(self):
+        """Update vectors of conserved quantities using previously
+        computed fluxes specified forcing functions.
+        """
+
+        from Numeric import ones, sum, equal, Float
+   
+            
+        N = self.number_of_elements
+        d = len(self.conserved_quantities)
+        
+        timestep = self.timestep
+
+        self.compute_forcing_terms()    
+
+        #Update conserved_quantities from explicit updates
+        for i, name in enumerate(self.conserved_quantities):
+            Q = self.quantities[name]
+            Q.update(timestep)
+            
+
+    def distribute_to_vertices_and_edges(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume
+
+        Default implementation is straight first order,
+        i.e. constant values throughout each element and
+        no reference to non-conserved quantities.
+        """
+
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.distribute_to_vertices_and_edges()
+
+
+
+
+def compute_fluxes_c(domain, max_timestep):
+    """Compute all fluxes and the timestep suitable for all volumes
+    This version works directly with consecutive data structure
+    and calls C-extension.
+    Note that flux function is hardwired into C-extension.
+    """
+    
+    import sys
+
+    neighbours = Volume.neighbours
+    neighbour_faces = Volume.neighbour_faces
+    normals = Volume.normals
+    flux = Volume.explicit_update    
+
+    area = Volume.geometric[:,0]            
+    radius = Volume.geometric[:,1]            
+    edgelengths = Volume.geometric[:,2:]
+
+
+    flux[:] = 0.0 #Reset stored fluxes to zero        
+    from domain_ext import compute_fluxes
+    timestep = compute_fluxes(domain.flux_function,
+                              Volume.conserved_quantities_face0,
+                              Volume.conserved_quantities_face1,
+                              Volume.conserved_quantities_face2,
+                              Volume.field_values_face0,
+                              Volume.field_values_face1,
+                              Volume.field_values_face2,
+                              Boundary_value.conserved_quantities,
+                              Boundary_value.field_values,
+                              Vector.coordinates,
+                              neighbours,
+                              neighbour_faces,
+                              normals,
+                              flux, area, radius, edgelengths,
+                              max_timestep)
+
+    domain.timestep = timestep
+
+
+
+    
+                
+##############################################
+#Initialise module
+
+#C-extensions
+#import util
+#if util.can_use_C_extension('domain_ext.c'):
+#    compute_fluxes = compute_fluxes_c
+    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
+    #update_conserved_quantities = update_conserved_quantities_c
+#else:
+#    from shallow_water import compute_fluxes
+    #from python_versions import distribute_to_vertices_and_edges
+    #from python_versions import update_conserved_quantities
Index: inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 195)
@@ -0,0 +1,211 @@
+"""boundary.py - Classes for implementing boundary conditions
+"""
+
+from util import rotate
+
+class Boundary:
+    """Base class for boundary conditions.
+       Specific boundary conditions must provide values for
+       the conserved_quantities
+
+       A boundary object has one neighbour; the one it
+       serves.
+    """   
+
+    def __init__(self):
+        pass
+        
+    def evaluate(self, vol_id=None, edge_id=None):
+        msg = 'Generic class Boundary must be subclassed'
+        raise msg
+
+
+class Transmissive_boundary(Boundary):
+    """Transmissive boundary returns same conserved quantities as
+    those present in its neighbour volume.
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+    
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for transmissive boundary'
+            raise msg
+        
+        self.domain = domain
+        
+    def __repr__(self):
+        return 'Transmissive_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Transmissive boundaries return the edge values
+        of the volume they serve.
+        """
+        
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        return q        
+
+
+class Dirichlet_boundary(Boundary):
+    """Dirichlet boundary returns constant values for the
+    conserved quantities
+    """
+
+
+    def __init__(self, conserved_quantities=None):
+        Boundary.__init__(self)
+
+        if conserved_quantities is None:
+            msg = 'Must specify one value for each conserved quantity'
+            raise msg
+
+        from Numeric import array, Float
+        self.conserved_quantities=array(conserved_quantities).astype(Float)
+
+    def __repr__(self):
+        return 'Dirichlet boundary (%s)' %self.conserved_quantities
+    
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.conserved_quantities
+
+
+
+class Time_boundary(Boundary):
+    """Time dependent boundary returns values for the
+    conserved quantities as a function of time.
+    Must specify domain to get access to model time and a function f(t)
+    which must return conserved quantities as a function time
+    """
+
+
+    def __init__(self, domain = None, f=None):
+        Boundary.__init__(self)
+
+        try:
+            q = f(0.0)
+        except:
+            raise 'Function for time boundary could not be executed'
+
+        from Numeric import array, Float        
+        try:
+            q = array(q).astype(Float)
+        except:
+            msg = 'Return value from time boundary function could '
+            msg += 'not be converted into a Numeric array of floats.\n'
+            msg += 'Specified function should return either list or array.' 
+            raise msg
+
+        msg = 'ERROR: Time boundary function must return a 1d list or array '
+        assert len(q.shape) == 1, msg
+            
+        d = len(domain.conserved_quantities)
+        msg = 'Return value for f must be a list or an array of length %d' %d
+        assert len(q) == d, msg
+
+        self.f = f
+        self.domain = domain
+
+    def __repr__(self):
+        return 'Time boundary'
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        return self.f(self.domain.time)
+
+
+
+#THIS FAR (10/8/4)
+class Connective_boundary(Boundary):
+    """Connective boundary returns values for the
+    conserved quantities from a volume as defined by a connection table
+    mapping between tuples of (volume id, face id) for volumes that
+    have their boundaries connected. 
+
+    FIXME: Perhaps include possibility of mapping between
+    different domains as well
+
+    FIXME: In case of shallow water we may want to have a
+    special version that casts this in terms of height rather than stage
+    """
+
+
+    def __init__(self, table):
+        from domain import Volume
+        
+        Boundary.__init__(self)
+
+        self.connection_table = table
+        self.Volume = Volume
+
+    def __repr__(self):
+        return 'Connective boundary'        
+
+    #FIXME: IF we ever need to get field_values from connected volume,
+    #that method could be overridden here (using same idea as in
+    #get_conserved_quantities
+    #def get_field_values()
+    
+    def get_conserved_quantities(self, volume, face=0):
+
+        id = volume.id
+        if self.connection_table.has_key((id, face)):
+            other_id, other_face = self.connection_table[(id, face)] 
+            
+            other_volume = self.Volume.instances[other_id]
+            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
+            exec(cmd)        
+            return q
+        else:
+            msg = 'Volume, face tuple ($d, %d) has not been mapped'\
+                  %(id, face)
+            raise msg
+            
+        
+
+
+
+#FIXME: Add a boundary with a general function of x,y, and t 
+
+#FIXME: Add periodic boundaries e.g.:
+# Attempt at periodic conditions from advection_spik. Remember this
+#
+#first = 2*(N-1)*N
+#for i in range(1,2*N+1,2):
+#    k = first + i-1#
+#
+#    print i,k
+#    
+#    domain[i].faces[2].neighbour = domain[k].faces[1]
+#    domain[k].faces[1].neighbour = domain[i].faces[2]    
+    
+
+
+class General_boundary(Boundary):
+    """General boundary which can compute conserved quantities based on
+    their previous value, conserved quantities of its neighbour and model time.
+
+    Must specify initial conserved quantities,
+    neighbour,
+    domain to get access to model time 
+    a function f(q_old, neighbours_q, t) which must return
+    new conserved quantities q as a function time
+
+    FIXME: COMPLETE UNTESTED - JUST AN IDEA
+    """
+
+    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
+        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)
+
+        self.f = f
+        self.domain = domain
+
+
+    def get_conserved_quantities(self, volume=None, face=0):
+        return self.f(self.conserved_quantities,
+                      neighbour.conserved_quantities,
+                      self.domain.time)
+
+    
+    
+        
Index: inundation/ga/storm_surge/pyvolution/interface.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/interface.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/interface.py	(revision 195)
@@ -0,0 +1,2 @@
+"""Public 
+"""
Index: inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/mesh.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/mesh.py	(revision 195)
@@ -0,0 +1,566 @@
+"""Classes implementing general 2D geometrical mesh of triangles.
+
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia, 2004   
+"""
+
+
+class Mesh:
+    """Collection of triangular elements (purely geometric)
+
+    A triangular element is defined in terms of three vertex ids,
+    ordered counter clock-wise, 
+    each corresponding to a given coordinate set.
+    Vertices from different elements can point to the same
+    coordinate set.
+
+    Coordinate sets are implemented as an N x 2 Numeric array containing
+    x and y coordinates. 
+    
+
+    To instantiate:
+       Mesh(points, vertices)
+
+    where
+
+      points is either a list of 2-tuples or an Mx2 Numeric array of
+      floats representing all x, y coordinates in the mesh.
+
+      vertices is either a list of 3-tuples or an Nx3 Numeric array of
+      integers representing indices of all vertices in the mesh.
+      Each vertex is identified by its index i in [0, M-1].
+
+
+    Example:
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        e = [2.0, 2.0]
+        
+        points = [a, b, c, e]
+        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        mesh = Mesh(points, vertices)
+    
+        #creates two triangles: bac and bce
+
+
+    Mesh takes the optional third argument boundary which is a
+    dictionary mapping from (element_id, edge_id) to boundary tag.
+    The default value is None which will assign the defualt_boundary_tag
+    as specified in config.py to all boundary edges.
+    """
+
+    def __init__(self, coordinates, vertices, boundary = None):
+        """
+        Build triangles from x,y coordinates (sequence of 2-tuples or
+        Mx2 Numeric array of floats) and vertices (sequence of 3-tuples
+        or Nx3 Numeric array of non-negative integers).
+        """
+
+        from Numeric import array, zeros, Int, Float, maximum
+        from util import distance, length
+
+        #Only one of them
+        self.vertices = array(vertices).astype(Int)
+        self.coordinates = array(coordinates)
+
+
+        #Input checks
+        msg = 'Vertices must an Nx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.vertices.shape) == 2, msg
+
+        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.coordinates.shape) == 2, msg        
+
+        msg = 'Vertex indices reference non-existing coordinate sets'
+        assert max(max(self.vertices)) <= self.coordinates.shape[0], msg
+
+
+        #Register number of elements (N)
+        self.number_of_elements = N = self.vertices.shape[0]
+
+        #Allocate space for geometric quantities
+        self.centroids = zeros((N, 2), Float)
+        self.areas = zeros(N, Float)
+        self.radii = zeros(N, Float)
+        self.edgelengths = zeros((N, 3), Float)
+
+        self.neighbours = zeros((N, 3), Int)
+        self.neighbour_edges = zeros((N, 3), Int)
+        self.number_of_boundaries = zeros(N, Int) 
+        self.surrogate_neighbours = zeros((N, 3), Int)        
+        
+        self.normals = zeros((N, 6), Float)
+
+        #Get x,y coordinates for all vertices for all triangles
+        #and store
+        self.vertex_coordinates = V = self.compute_vertex_coordinates()
+        
+
+        ##print 'Initialise mesh'
+        #Initialise each triangle
+        for i in range(N):
+            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
+            
+            x0 = V[i, 0]; y0 = V[i, 1]
+            x1 = V[i, 2]; y1 = V[i, 3]
+            x2 = V[i, 4]; y2 = V[i, 5]            
+
+            #Compute centroid
+            centroid = [(x0 + x1 + x2)/3, (y0 + y1 + y2)/3]
+            self.centroids[i] = centroid
+
+            #Area
+            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+
+            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
+            msg += ' is degenerate:  area == %f' %self.areas[i]
+            assert self.areas[i] > 0.0, msg
+        
+            #Midpoints        
+            m0 = [(x1 + x2)/2, (y1 + y2)/2]
+            m1 = [(x0 + x2)/2, (y0 + y2)/2]
+            m2 = [(x1 + x0)/2, (y1 + y0)/2]
+
+            #The radius is the distance from the centroid of
+            #a triangle to the midpoint of the side of the triangle
+            #closest to the centroid
+            d0 = distance(centroid, m0)
+            d1 = distance(centroid, m1)        
+            d2 = distance(centroid, m2)
+        
+            self.radii[i] = min(d0, d1, d2)
+
+            #Normals
+            #The normal vectors
+            # - point outward from each edge
+            # - are orthogonal to the edge
+            # - have unit length
+            # - Are enumerated according to the opposite corner:
+            #   (First normal is associated with the edge opposite
+            #    the first vertex, etc)
+            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
+
+            n0 = array([x2 - x1, y2 - y1])
+            l0 = length(n0)
+
+            n1 = array([x0 - x2, y0 - y2])
+            l1 = length(n1)
+
+            n2 = array([x1 - x0, y1 - y0])
+            l2 = length(n2)
+
+            #Normalise
+            n0 /= l0
+            n1 /= l1
+            n2 /= l2
+
+            #Compute and store
+            self.normals[i, :] = [n0[1], -n0[0],
+                                  n1[1], -n1[0],
+                                  n2[1], -n2[0]]             
+
+            #Edgelengths
+            self.edgelengths[i, :] = [l0, l1, l2]        
+
+            #Initialise Neighbours (-1 means that it is a boundary neighbour)
+            self.neighbours[i, :] = [-1, -1, -1]
+
+            #Initialise edge ids of neighbours
+            #In case of boundaries this slot is not used
+            self.neighbour_edges[i, :] = [-1, -1, -1]
+
+
+        ##print 'NB'    
+        #Build neighbour structure    
+        self.build_neighbour_structure()
+
+        ##print 'VL'    
+        #Build vertex list
+        self.build_vertexlist()        
+
+        ##print 'SNB'    
+        #Build surrogate neighbour structure    
+        self.build_surrogate_neighbour_structure()        
+
+        #Build boundary dictionary mapping (id, edge) to symbolic tags        
+        self.build_boundary_dictionary(boundary)
+
+        #Update boundary indices 
+        self.build_boundary_structure()        
+
+        
+    def __len__(self):
+        return self.number_of_elements
+
+    def __repr__(self):
+        return 'Mesh: %d vertices, %d elements, %d boundary segments'\
+               %(self.coordinates.shape[0], len(self), len(self.boundary))
+
+            
+    def build_neighbour_structure(self):
+        """Update all registered triangles to point to their neighbours.
+
+        Postconditions:
+          neighbours and neighbour_edges is populated
+          number_of_boundaries integer array is defined.
+        """
+        
+        #Step 1:
+        #Build dictionary mapping from segments (2-tuple of points)
+        #to left hand side edge (facing neighbouring triangle)
+        #Also build list of vertices containing indices of triangles  
+
+        N = self.number_of_elements        
+        neighbourdict = {}
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(N):
+
+            #Register all segments as keys mapping to current triangle
+            #and segment id
+            a = self.vertices[i, 0]
+            b = self.vertices[i, 1]
+            c = self.vertices[i, 2]            
+            neighbourdict[a,b] = (i, 2) #(id, edge)
+            neighbourdict[b,c] = (i, 0) #(id, edge)
+            neighbourdict[c,a] = (i, 1) #(id, edge)
+
+            #Register the vertices as keys mapping to 
+            #(triangle, edge) tuples associated with them
+            #This is used for smoothing
+            for vertex_id, v in enumerate([a,b,c]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )    
+
+        #Step 2:
+        #Go through triangles again, but this time
+        #reverse direction of segments and lookup neighbours.
+        for i in range(N):
+            a = self.vertices[i, 0]
+            b = self.vertices[i, 1]
+            c = self.vertices[i, 2]
+
+            self.number_of_boundaries[i] = 3
+            if neighbourdict.has_key((b,a)):
+                self.neighbours[i, 2] = neighbourdict[b,a][0]
+                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
+                self.number_of_boundaries[i] -= 1
+                
+            if neighbourdict.has_key((c,b)):
+                self.neighbours[i, 0] = neighbourdict[c,b][0]
+                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
+                self.number_of_boundaries[i] -= 1               
+                
+            if neighbourdict.has_key((a,c)):
+                self.neighbours[i, 1] = neighbourdict[a,c][0]
+                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
+                self.number_of_boundaries[i] -= 1              
+
+
+        self.vertexlist = vertexlist
+    
+    def build_vertexlist(self):
+        """Build vertexlist index by vertex ids and for each entry (point id)
+        build a list of (triangles, vertex_id) pairs that use the point
+        as vertex.
+
+        Preconditions:
+          self.coordinates and self.vertices are defined  
+        
+        Postcondition:
+          self.vertexlist is build
+        """
+
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(self.number_of_elements):
+
+            a = self.vertices[i, 0]
+            b = self.vertices[i, 1]
+            c = self.vertices[i, 2]            
+
+            #Register the vertices as keys mapping to 
+            #(triangle, edge) tuples associated with them
+            #This is used for smoothing
+            for vertex_id, v in enumerate([a,b,c]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )    
+
+        self.vertexlist = vertexlist
+    
+
+    def build_surrogate_neighbour_structure(self):
+        """Build structure where each triangle edge points to its neighbours
+        if they exist. Otherwise point to the triangle itself.
+
+        Also, keep a tally of the number of boundaries for each triangle
+
+        The surrogate neighbour structure is useful for computing gradients
+        based on centroid values of neighbours.
+
+        Precondition: Neighbour structure is defined
+        Postcondition:
+          Surrogate neighbour structure is defined:
+          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
+          triangles.
+
+        """
+
+        N = self.number_of_elements
+        for i in range(N):
+            #Find all neighbouring volumes that are not boundaries
+            for k in range(3):
+                if self.neighbours[i, k] < 0:
+                    self.surrogate_neighbours[i, k] = i #Point this triangle
+                else:
+                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
+
+
+
+    def build_boundary_dictionary(self, boundary = None):
+        """Build or input the dictionary of boundary tags.
+         self.boundary is a dictionary of tags,
+         keyed by volume id and edge:
+         { (id, edge): tag, ... }
+         
+         Postconditions:
+            self.boundary is defined.
+        """
+
+        from config import default_boundary_tag
+        
+        if boundary is None:
+            boundary = {}
+            for vol_id in range(self.number_of_elements):
+                for edge_id in range(0, 3):
+                    if self.neighbours[vol_id, edge_id] < 0:
+                        boundary[(vol_id, edge_id)] = default_boundary_tag
+        else:
+            #Check that all keys in given boundary exist
+            for vol_id, edge_id in boundary.keys():
+                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+                a, b = self.neighbours.shape
+                assert vol_id < a and edge_id < b, msg
+                
+                msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
+                assert self.neighbours[vol_id, edge_id] < 0, msg
+
+            #Check that all boundary segments are assigned a value
+            for vol_id in range(self.number_of_elements):
+                for edge_id in range(0, 3):
+                    if self.neighbours[vol_id, edge_id] < 0:
+                        if not boundary.has_key( (vol_id, edge_id) ):
+                            msg = 'WARNING: Given boundary does not contain '
+                            msg += 'tags for edge (%d, %d). '\
+                                   %(vol_id, edge_id)
+                            msg += 'Assigning default tag (%s).'\
+                                   %default_boundary_tag
+
+                            #FIXME: Print only as per verbosity
+                            #print msg
+                            boundary[ (vol_id, edge_id) ] =\
+                                      default_boundary_tag
+            
+
+            
+        self.boundary = boundary
+            
+
+    def build_boundary_structure(self):
+        """Traverse boundary and 
+        enumerate neighbour indices from -1 and
+        counting down.
+
+        Precondition:
+            self.boundary is defined.
+        Post condition:
+            neighbour array has unique negative indices for boundary
+            boundary_segments array imposes an ordering on segments
+            (not otherwise available from the dictionary)             
+        """
+
+        if self.boundary is None:
+            msg = 'Boundary dictionary must be defined before '
+            msg += 'building boundary structure'
+            raise msg
+
+        
+        self.boundary_segments = self.boundary.keys()
+        self.boundary_segments.sort()
+
+        index = -1        
+        for id, edge in self.boundary_segments:
+            self.neighbours[id, edge] = index
+            index -= 1
+
+    
+    def get_boundary_tags(self):
+        """Return list of available boundary tags
+        """
+        
+        tags = {}
+        for v in self.boundary.values():
+            tags[v] = 1
+
+        return tags.keys()    
+            
+    
+
+    def check_integrity(self):
+        """Check that triangles are internally consistent e.g.
+        that area corresponds to edgelengths, that vertices
+        are arranged in a counter-clockwise order, etc etc
+        Neighbour structure will be checked by class Mesh
+        """
+        
+        from config import epsilon
+        from math import pi
+        from util import anglediff, inner_product, length
+
+        N = self.number_of_elements
+        
+        #Get x,y coordinates for all vertices for all triangles
+        V = self.get_vertex_coordinates()
+
+        #Check each triangle
+        for i in range(N):
+            x0 = V[i, 0]; y0 = V[i, 1]
+            x1 = V[i, 2]; y1 = V[i, 3]
+            x2 = V[i, 4]; y2 = V[i, 5]            
+
+            #Check that area hasn't been compromised
+            area = self.areas[i]
+            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
+                  %(x0,y0,x1,y1,x2,y2)
+            assert abs((area - ref)/area) < epsilon, msg
+        
+            #Check that points are arranged in counter clock-wise order
+            v0 = [x1-x0, y1-y0]
+            v1 = [x2-x1, y2-y1]
+            v2 = [x0-x2, y0-y2]
+
+            a0 = anglediff(v1, v0)
+            a1 = anglediff(v2, v1)
+            a2 = anglediff(v0, v2)        
+
+            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
+            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
+            assert a0 < pi and a1 < pi and a2 < pi, msg
+
+            #Check that normals are orthogonal to edge vectors
+            #Note that normal[k] lies opposite vertex k
+
+            normal0 = self.normals[i, 0:2]
+            normal1 = self.normals[i, 2:4]
+            normal2 = self.normals[i, 4:6]            
+
+            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
+                assert u[0]*v[0] + u[1]*v[1] < epsilon
+
+
+        #Check integrity of neighbour structure
+        for i in range(N):
+            for v in self.vertices[i, :]:
+                #Check that all vertices have been registered
+                assert self.vertexlist[v] is not None
+
+                #Check that this triangle is listed with at least one vertex
+                assert (i, 0) in self.vertexlist[v] or\
+                       (i, 1) in self.vertexlist[v] or\
+                       (i, 2) in self.vertexlist[v]
+                
+                
+
+            #Check neighbour structure
+            for k, neighbour_id in enumerate(self.neighbours[i,:]):
+
+                #Assert that my neighbour's neighbour is me
+                #Boundaries need not fulfill this
+                if neighbour_id >= 0:
+                    edge = self.neighbour_edges[i, k]
+                    assert self.neighbours[neighbour_id, edge] == i
+
+
+
+        #Check that all boundaries have
+        # unique, consecutive, negative indices                    
+        L = len(self.boundary)
+        for i in range(L):
+            id, edge = self.boundary_segments[i]
+            assert self.neighbours[id, edge] == -i-1
+
+
+        #NOTE: This assert doesn't hold true if there are internal boundaries
+        #FIXME: Look into this further.
+        #for id, edge in self.boundary:
+        #    assert self.neighbours[id,edge] < 0
+        
+
+    
+
+    def get_normals(self):
+        """Return all normal vectors.
+        Return normal vectors for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.normals
+
+
+    def get_normal(self, i, j):
+        """Return normal vector j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.normals[i, 2*j:2*j+2]     
+    
+
+            
+    def get_vertex_coordinates(self):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.vertex_coordinates
+
+
+    def get_vertex_coordinate(self, i, j):
+        """Return coordinates for vertex j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.vertex_coordinates[i, 2*j:2*j+2]     
+
+
+    def compute_vertex_coordinates(self):
+        """Return vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+        """
+
+        from Numeric import zeros, Float
+
+        N = self.number_of_elements
+        vertex_coordinates = zeros((N, 6), Float) 
+
+        for i in range(N):
+            for j in range(3):
+                k = self.vertices[i,j]  #Index of vertex 0
+                v_k = self.coordinates[k]
+                vertex_coordinates[i, 2*j+0] = v_k[0]
+                vertex_coordinates[i, 2*j+1] = v_k[1]
+
+        return vertex_coordinates
+  
+
+
+
+def rectangular_mesh(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+    from mesh_factory import rectangular
+    points, vertices, boundary = rectangular(m, n, len1, len2, origin)
+
+    return Mesh(points, vertices, boundary)
+
+
+#FIXME: Get oblique and contracting and circular meshes from Chris
+
Index: inundation/ga/storm_surge/pyvolution/mesh_factory.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 195)
@@ -0,0 +1,65 @@
+"""Library of standard meshes and facilities for reading various
+mesh file formats
+"""
+
+
+def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+    """Setup a rectangular grid of triangles
+    with m+1 by n+1 grid points
+    and side lengths len1, len2. If side lengths are omitted
+    the mesh defaults to the unit square.
+    
+    len1: x direction (left to right)
+    len2: y direction (bottom to top)
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements) 
+    """
+
+    from config import epsilon
+    
+    #E = m*n*2        #Number of triangular elements
+    #P = (m+1)*(n+1)  #Number of initial vertices
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n    
+
+    #Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        for j in range(n+1):
+            vertices[i,j] = len(points)
+            points.append([i*delta1 + origin[0], j*delta2 + origin[1]])
+
+
+
+    #Construct 2 triangles per rectangular element and assign tags to boundary
+    elements = []
+    boundary = {}
+    for i in range(m):
+        for j in range(n):
+            v1 = vertices[i,j+1]
+            v2 = vertices[i,j]            
+            v3 = vertices[i+1,j+1]            
+            v4 = vertices[i+1,j]            
+
+            #Update boundary dictionary and create elements
+            if i == m-1:
+                boundary[(len(elements), 2)] = 'right'
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'                
+            elements.append([v4,v3,v2]) #Lower element
+
+            if i == 0:
+                boundary[(len(elements), 2)] = 'left'
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top' 
+            elements.append([v1,v2,v3]) #Upper element   
+
+    return points, elements, boundary       
+
+
+
Index: inundation/ga/storm_surge/pyvolution/python_versions.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/python_versions.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/python_versions.py	(revision 195)
@@ -0,0 +1,2 @@
+    
+    
Index: inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/quantity.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/quantity.py	(revision 195)
@@ -0,0 +1,265 @@
+"""Class Quantity - Implements values at each triangular element
+
+To create:
+
+   Quantity(domain, vertex_values)
+
+   domain: Associated domain structure. Required.
+   
+   vertex_values: N x 3 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
+"""
+
+#FIXME: Make Conserved_quantity a subclass of Quantity 
+
+from mesh import Mesh
+
+class Quantity:
+
+    def __init__(self, domain, vertex_values=None):
+
+        from Numeric import array, zeros, Float
+
+        msg = 'First argument in Quantity.__init__ '
+        msg += 'must be of class Mesh (or a subclass thereof)'
+        assert isinstance(domain, Mesh), msg
+
+        if vertex_values is None:
+            N = domain.number_of_elements            
+            self.vertex_values = zeros((N, 3), Float)
+        else:    
+            self.vertex_values = array(vertex_values)
+
+            N, V = self.vertex_values.shape
+            assert V == 3,\
+                   'Three vertex values per element must be specified'
+            
+
+            msg = 'Number of vertex values (%d) must be consistent with'\
+                  %N
+            msg += 'number of elements in specified domain (%d).'\
+                   %domain.number_of_elements
+            
+            assert N == domain.number_of_elements, msg
+
+        self.domain = domain    
+
+        #Allocate space for other quantities
+        self.centroid_values = zeros(N, Float)
+        self.edge_values = zeros((N, 3), Float)        
+
+        #Allocate space for boundary values
+        L = len(domain.boundary)
+        self.boundary_values = zeros(L, Float)
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+        self.explicit_update = zeros(N, Float )
+        self.semi_implicit_update = zeros(N, Float )
+                
+        #Intialise centroid and edge_values
+        self.interpolate()
+
+
+    def interpolate(self):
+        """Compute interpolated values at edges and centroid
+        Pre-condition: vertex_values have been set
+        """
+
+        #FIXME: Write using vector operations (or in C)
+        N = self.vertex_values.shape[0]        
+        for i in range(N):
+            v0 = self.vertex_values[i, 0]
+            v1 = self.vertex_values[i, 1]
+            v2 = self.vertex_values[i, 2]
+            
+            self.centroid_values[i] = (v0 + v1 + v2)/3
+            
+            self.edge_values[i, 0] = 0.5*(v1 + v2)
+            self.edge_values[i, 1] = 0.5*(v0 + v2)            
+            self.edge_values[i, 2] = 0.5*(v0 + v1)
+
+    def update(self, timestep):
+        """Update centroid values based on values stored in
+        explicit_update and semi_implicit_update as well as given timestep
+        """
+
+        from Numeric import sum, equal, ones, Float
+        
+        N = self.centroid_values.shape[0]
+        
+        #Explicit updates
+        self.centroid_values += timestep*self.explicit_update
+            
+        #Semi implicit updates
+        denominator = ones(N, Float)-timestep*self.semi_implicit_update
+
+        if sum(equal(denominator, 0.0)) > 0.0:
+            msg = 'Zero division in semi implicit update. Call Stephen :-)'
+            raise msg
+        else:
+            #Update conserved_quantities from semi implicit updates
+            self.centroid_values /= denominator
+
+
+    def second_order_limiter(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using second order scheme.
+
+        precondition:
+           vertex values are estimated from gradient
+        postcondition:
+           vertex values are updated
+        """
+
+        from Numeric import zeros, Float
+
+        beta = self.domain.beta
+        
+        qc = self.centroid_values
+        qv = self.vertex_values
+        
+        #Deltas between vertex and centroid values
+        dq = zeros( qv.shape, Float)
+        for i in range(3):
+            dq[:,i] = qv[:,i] - qc
+        
+        #Find min and max at this and neighbour's centroids
+        qmax = zeros(qc.shape, Float)
+        qmin = zeros(qc.shape, Float)        
+
+        for k in range(self.domain.number_of_elements):
+            qmax[k] = qmin[k] = qc[k]
+            for i in range(3):
+                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
+  
+        for k in range(self.domain.number_of_elements):
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(3):
+                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(3):            
+                qv[k,i] = qc[k] + phi*dq[k,i]
+
+        
+    def first_order_limiter(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(3):
+            qv[:,i] = qc
+            
+
+    def distribute_to_vertices_and_edges(self, order=1):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume using
+        specified order.
+
+        This is an unsophisticated limiter as it does not take into
+        account dependencies among quantities.
+        """
+
+        if order == 1:
+            self.first_order_limiter()
+        elif order == 2:
+            self.second_order_limiter()
+        else:
+            msg = 'ERROR: unknown order %d' %order
+            raise msg
+
+        self.interpolate_from_vertices_to_edges()
+
+
+    def interpolate_from_vertices_to_edges(self):
+        #FIXME: Write using vector operations (or in C)
+        for k in range(self.vertex_values.shape[0]):
+            q0 = self.vertex_values[k, 0]
+            q1 = self.vertex_values[k, 1]
+            q2 = self.vertex_values[k, 2]
+            
+            self.edge_values[k, 0] = 0.5*(q1+q2)
+            self.edge_values[k, 1] = 0.5*(q0+q2)  
+            self.edge_values[k, 2] = 0.5*(q0+q1)
+
+            
+            
+    def set_values(self, values, location='vertices'):
+        """Set values for quantity
+
+        values: Compatible list or Numeric array (see below)
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If 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.
+        """
+
+        #FIXME: Should take functions as argument as well
+        #FIXME: Should take constants as well
+        
+        from Numeric import array, Float
+
+        values = array(values).astype(Float)
+
+        N = self.centroid_values.shape[0]
+        
+        msg = 'Number of values must match number of elements'
+        assert values.shape[0] == N, msg
+        
+        if location == 'vertices':
+            assert len(values.shape) == 2, 'Values array must be 2d'
+            msg = 'Array must be N x 3'            
+            assert values.shape[1] == 3, msg
+            
+            self.vertex_values = values
+            
+            #Intialise centroid and edge_values
+            self.interpolate()
+        elif location == 'edges':
+            assert len(values.shape) == 2, 'Values array must be 2d'
+            msg = 'Array must be N x 3'            
+            assert values.shape[1] == 3, msg
+            
+            self.edge_values = values
+            
+        elif location == 'centroids':
+            assert len(values.shape) == 1, 'Values array must be 1d'
+            self.centroid_values = values                        
+        else:
+            raise 'Invalid location: %s' %location
Index: inundation/ga/storm_surge/pyvolution/realtime_visualisation.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/realtime_visualisation.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/realtime_visualisation.py	(revision 195)
@@ -0,0 +1,268 @@
+from visual import *
+
+def normal(v1,v2):
+    """Safe computation of normalised cross product
+    """
+    try:
+        n = norm(cross(v1, v2))
+    except ZeroDivisionError:
+        print v1, v2
+        raise Exception
+        #n = vector(0,0,0)
+
+    return n    
+
+
+class Triangle:
+    def __init__(self,v0,v1,v2,color=(1,1,1), frame=None):
+        """Create one flat triangular panel with two sides
+        """
+
+        #Outward normal
+        
+        self.normal=normal(v1-v0, v2-v0)
+
+        self.panel = faces(frame = frame)
+        self.backpanel = faces(frame = frame)        
+        for v in (v0,v1,v2):
+            self.panel.append(pos=v, normal = self.normal, color=color)
+
+        for v in (v2,v1,v0):
+            self.backpanel.append(pos=v, normal = -self.normal, color=color)            
+
+            
+    def move(self,v):
+        """Move panel in direction given by vector v
+        """
+
+        for vertex in self.panel.pos:
+            vertex += v            
+
+    def update_height(self,d,fix_baseline=False):
+        """Change height of triangle by displacing the third vertex by
+        d in the direction perpendicular to the baseline (v1-v0) and the
+        outward normal vector
+        """
+
+        v0 = self.panel.pos[0]
+        v1 = self.panel.pos[1]
+        v2 = self.panel.pos[2]        
+        
+        n = normal(v1-v0, self.normal)
+
+        if fix_baseline:
+            self.panel.pos[2] -= d*n
+        else:
+            self.panel.pos[:2] += d*n        
+
+    def update_color(self, c):
+        """Change color to c
+        """
+
+        self.panel.color = c
+        
+
+    def set_vertexheights(self, heights, floor_heights = None):
+        from Numeric import zeros, Float
+
+        from config import minimum_allowed_height as hmin
+        if floor_heights is None:
+            floor_heights = zeros(heights.shape, Float)
+
+        all_vertices_below_threshold = True    
+        for k in range(3):
+            w = heights[k]
+            z = floor_heights[k]
+
+            if w-z >= hmin:
+                all_vertices_below_threshold = False
+            
+            vertex = self.panel.pos[k]
+            vertex[2] = w
+
+            vertex = self.backpanel.pos[2-k]
+            vertex[2] = w
+
+
+        #Update color to visualise dry areas 
+        if all_vertices_below_threshold: 
+            self.panel.color = self.bottom_color
+            self.backpanel.color = self.bottom_color                
+        else:
+            self.panel.color = self.top_color
+            self.backpanel.color = self.top_color
+
+
+        #update normal    
+        v0 = self.panel.pos[0]
+        v1 = self.panel.pos[1]
+        v2 = self.panel.pos[2]
+
+
+        n = normal(v1-v0, v2-v0)
+        self.panel.normal=n
+        self.backpanel.normal=-n
+
+        
+
+
+def create_surface(domain):
+    """Link surface of domains to their visual representations
+    """
+
+    fr = frame() #Default frame to contain all objects
+    s=0
+
+    Q = domain.quantities['level']
+    try:
+        Z = domain.quantities['elevation']
+        bed = True
+    except:
+        bed = False
+        
+    N = Q.vertex_values.shape[0]
+
+    domain.visuals = []
+    for i in range(N):
+
+        z0 = Q.vertex_values[i, 0]
+        z1 = Q.vertex_values[i, 1]
+        z2 = Q.vertex_values[i, 2]        
+
+        x0, y0 = domain.get_vertex_coordinate(i,0)
+        x1, y1 = domain.get_vertex_coordinate(i,1)
+        x2, y2 = domain.get_vertex_coordinate(i,2)        
+        
+        V0 = vector(x0, y0, z0)
+        V1 = vector(x1, y1, z1)
+        V2 = vector(x2, y2, z2)
+
+
+        #Top surface
+        c0 = 0.1
+        c1 = 0.4
+        c2 = 0.99
+        if s:
+            c2 = 0.7*c2   #To show triangles in slightly different shades
+
+        s = 1-s
+
+        col = (c0,c1,c2)
+            
+        visual_top = Triangle(V0,V1,V2,color=col,frame=fr)
+        visual_top.top_color = col
+
+    
+        #Bottom surface
+        #v0, v1, v2 = volume.vertices
+        #v0 = vector(v0.x, v0.y, v0.field_values[index])
+        #v1 = vector(v1.x, v1.y, v1.field_values[index])
+        #v2 = vector(v2.x, v2.y, v2.field_values[index])
+        if bed:
+            z0 = Z.vertex_values[i, 0]
+            z1 = Z.vertex_values[i, 1]
+            z2 = Z.vertex_values[i, 2]
+        else:
+            z0 = z1 = z2 = 0.0
+
+        V0 = vector(x0, y0, z0)
+        V1 = vector(x1, y1, z1)
+        V2 = vector(x2, y2, z2)        
+
+        c0 = 0.3
+        c1 = 0.3
+        c2 = 0.3
+        if s:
+            c2 = 0.4*c2   #To show triangles in slightly different shades
+
+        col = (c0,c1,c2)   
+        visual_bottom = Triangle(V0, V1, V2, color=col,frame=fr)
+        visual_top.bottom_color=col
+        
+        domain.visuals.append( (visual_top, visual_bottom) )
+        
+    update(domain)
+    #print 'Scale', scene.scale    
+
+
+def update(domain):    
+    """Update vertex heights.
+    The argument index refers to which conserved quantity to visualise.    
+    If domain.smooth is set True, vertex values will be averaged
+    yielding a smoother surface.
+    """
+
+    from Numeric import array
+
+
+    Q = domain.quantities['level']
+    N = Q.vertex_values.shape[0]
+    
+    #print scene.forward
+    #FIXME: Use smoother from pyvolution instead
+    if domain.smooth:
+        #Get all average point values
+        vertex_heights = {}
+
+        for k in range(N):
+            for i in range(3):
+                vertex = domain.vertices[k, i]
+                if vertex_heights.has_key(vertex):
+                    vertex_heights[vertex].append(
+                        Q.vertex_values[k, i])
+                else:
+                    vertex_heights[vertex] = []
+                    vertex_heights[vertex].append(
+                        Q.vertex_values[k, i])
+
+
+
+    for k in range(N):    
+        if domain.smooth:
+            #The averages
+            x = zeros(3, Float)
+            for i in range(3):            
+                vertex = domain.vertices[k, i]
+                A = array(vertex_heights[vertex])
+                x[i] = sum(A)/len(A)
+        else:        
+            #The true values
+            x = [Q.vertex_values[k, 0],
+                 Q.vertex_values[k, 1],
+                 Q.vertex_values[k, 2]]
+
+            
+        #Do it
+        floor_heights = array([pos[2] for pos in domain.visuals[k][1].panel.pos])
+        
+        domain.visuals[k][0].set_vertexheights(x, floor_heights)
+
+
+
+
+scene.width = 1000
+scene.height = 800
+
+#Original
+scene.center = (0.5,0.5,0)
+scene.forward = vector(-0.1, 0.5, -0.5)
+
+#Temporary (for bedslope)
+#scene.forward = vector(0.0006, 0.7, -0.03)    
+
+
+#Temporary for hackett - begin
+#scene.autoscale = 0
+#scene.scale = (0.002, 0.002, 0.01) #Scale z so that countours stand out more
+#scene.center = (300.0,500.0,-10)
+#Temporary for hackett - end
+
+    
+scene.ambient = 0.4
+scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
+                (-0.2, 0.2, 0.1)]
+
+
+
+
+
Index: inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 195)
@@ -0,0 +1,829 @@
+"""Class Domain -
+2D triangular domains for finite-volume computations of
+the shallow water wave equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the Shallow Water Wave Equation
+
+FIXME: Write equations here!
+
+
+Conserved quantities are w (water level or stage), uh (x momentum)
+and vh (y momentum).
+
+
+Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2004   
+"""
+
+from domain import *
+Generic_domain = Domain #Rename
+
+class Domain(Generic_domain):
+
+    def __init__(self, coordinates, vertices, boundary = None):
+
+        conserved_quantities = ['level', 'xmomentum', 'ymomentum']
+        other_quantities = ['elevation', 'friction'] 
+        
+        Generic_domain.__init__(self, coordinates, vertices, boundary,
+                                conserved_quantities, other_quantities)
+
+        from config import minimum_allowed_height
+        self.minimum_allowed_height = minimum_allowed_height
+
+        
+
+
+    def check_integrity(self):
+        Generic_domain.check_integrity(self)
+
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "level"'
+        assert self.conserved_quantities[0] == 'level', msg
+        msg = 'Second conserved quantity must be "xmomentum"'
+        assert self.conserved_quantities[1] == 'xmomentum', msg
+        msg = 'Third conserved quantity must be "ymomentum"'
+        assert self.conserved_quantities[2] == 'ymomentum', msg
+        
+
+        #Check that levels are >= bed elevation
+        from Numeric import alltrue, greater_equal
+        
+        level = self.quantities['level']
+        bed = self.quantities['elevation']        
+
+        msg = 'All water levels must be greater than the bed elevation'
+        assert alltrue( greater_equal(
+            level.vertex_values, bed.vertex_values )), msg
+
+        assert alltrue( greater_equal(
+            level.edge_values, bed.edge_values )), msg
+
+        assert alltrue( greater_equal(
+            level.centroid_values, bed.centroid_values )), msg        
+
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes(self)
+
+    def first_order_limiter(self):
+        #Call correct module function
+        first_order_limiter(self)        
+        
+####################################
+# Flux computation        
+def 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, vh are stored as elements 0, 1 and 2
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    #FIXME: Remove those and pass in height directly
+    #(unless we'll include a bed elevation discontinuity later)
+    
+    """
+
+    from config import g
+    from math import sqrt
+    from Numeric import array
+    from util import rotate
+
+    #Align momentums with x-axis
+    q_left  = rotate(ql, normal, direction = 1)
+    q_right = rotate(qr, normal, direction = 1)
+
+    z = (zl+zr)/2 #Take average of field values
+    #FIXME: Maybe allow discontinuity later
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    try:
+        u_left  = uh_left/h_left
+    except:
+        u_left = 0.0
+        h_left = 0.0
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    try:
+        u_right  = uh_right/h_right
+    except:
+        u_right = 0.0
+        h_right = 0.0
+
+        
+    vh_left  = q_left[2]
+    vh_right = q_right[2]        
+
+    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([uh_left,
+                        u_left**2*h_left + 0.5*g*h_left**2,
+                        vh_left*u_left])
+    flux_right = array([uh_right,
+                        u_right**2*h_right + 0.5*g*h_right**2,
+                        vh_right*u_right])    
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        flux = array([0.0, 0.0, 0.0])
+        max_speed = 0.0
+    else:    
+        flux = (s_max*flux_left - s_min*flux_right)/denom
+        flux += s_max*s_min*(q_right-q_left)/denom
+
+        flux = rotate(flux, normal, direction=-1)
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return flux, max_speed        
+
+
+def 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
+    domain.explicit_update
+
+    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
+    from config import max_timestep    
+
+    N = domain.number_of_elements
+    
+    neighbours = domain.neighbours
+    neighbour_edges = domain.neighbour_edges
+    normals = domain.normals
+
+    areas = domain.areas
+    radii = domain.radii
+    edgelengths = domain.edgelengths
+        
+    timestep = max_timestep #FIXME: Get rid of this
+
+    #Shortcuts
+    Level = domain.quantities['level']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']        
+
+    #Arrays
+    level = Level.edge_values
+    xmom =  Xmom.edge_values
+    ymom =  Ymom.edge_values
+    bed =   Bed.edge_values    
+
+    level_bdry = Level.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    ymom_bdry =  Ymom.boundary_values
+    
+    flux = zeros(3, Float) #Work array for summing up fluxes
+
+    #Loop
+    for k in range(N):
+        optimal_timestep = float(sys.maxint)
+
+        flux[:] = 0.  #Reset work array
+        for i in range(3):
+            #Quantities inside volume facing neighbour i
+            ql = [level[k, i], xmom[k, i], ymom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = neighbours[k,i] 
+            if n < 0:
+                m = -n-1 #Convert neg flag to index
+                qr = [level_bdry[m], xmom_bdry[m], ymom_bdry[m]]
+                zr = zl #Extend bed elevation to boundary
+            else:    
+                m = neighbour_edges[k,i]
+                qr = [level[n, m], xmom[n, m], ymom[n, m]]
+                zr = bed[n, m]                
+
+                
+            #Outward pointing normal vector   
+            normal = normals[k, 2*i:2*i+2]
+
+            #Flux computation using provided function
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            flux -= edgeflux * edgelengths[k,i]
+            
+            #Update optimal_timestep
+            try:
+                optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= areas[k]
+        Level.explicit_update[k] = flux[0]
+        Xmom.explicit_update[k] = flux[1]
+        Ymom.explicit_update[k] = flux[2]
+        
+        timestep = min(timestep, optimal_timestep)
+
+    domain.timestep = timestep    
+	
+
+####################################
+# Functions for gradient limiting (distribute_to_vertices_and_edges)
+def first_order_limiter(domain):
+    """First order limiter function, specific to the shallow water wave
+    equation.
+    
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes (see comment in code)
+    In addition, momemtums get distributed as constant values.
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+    """
+
+    #Water levels at centroids
+    wc = domain.quantities['level'].centroid_values
+
+    #Bed elevations at centroids
+    zc = domain.quantities['elevation'].centroid_values
+
+    #Water depths at centroids    
+    hc = wc - zc
+
+    #Momentums at centroids
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values        
+
+    #Water levels at vertices
+    wv = domain.quantities['level'].vertex_values
+
+    #Bed elevations at vertices
+    zv = domain.quantities['elevation'].vertex_values
+    
+    for k in range(domain.number_of_elements):
+        #Protect against infinitesimal heights and high velocities
+        if hc[k] < domain.minimum_allowed_height:
+            #Control level and height
+            if hc[k] < 0.0:
+                wc[k] = zc[k]; hc[k] = 0.0
+
+            #Control momentum
+            xmomc[k] = ymomc[k] = 0.0
+                
+                
+                
+        #Compute maximal variation in bed elevation
+        z_range = max(abs(zv[k,0]-zc[k]),
+                      abs(zv[k,1]-zc[k]),
+                      abs(zv[k,2]-zc[k]))
+
+
+        #Weighted balance between stage parallel to bed elevation 
+        #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc)
+        #where i=0,1,2 denotes the vertex ids
+        #
+        #It follows that 
+        #  wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) = 
+        #  (1-alpha)*zvi + alpha*zc + hc = 
+        #  zvi + hc + alpha*(zc-zvi);
+        #
+        #where alpha in [0,1] and defined as the ration between hc and 
+        #the maximal difference from zc to zv0, zv1 and zv2
+        
+        if z_range > 0.0:
+            alpha = min(hc[k]/z_range, 1.0)
+        else: 
+            alpha = 1.0
+            
+        #Update stage
+        for i in range(3):
+            wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
+
+    #Update momentum using straight first order
+    for name in ['xmomentum', 'ymomentum']:
+        Q = domain.quantities[name]
+        for i in range(3):
+            Q.vertex_values[:, i] = Q.centroid_values
+
+        
+def second_order_limiter(domain):
+    """Second order limiter function, specific to the shallow water wave
+    equation.
+    
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes (see comment in code)
+
+    A weighted average between shallow 
+    and deep cases is as in the first order case.
+    
+    In addition, all conserved quantities get distributed as per a
+    piecewise linear function.
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+    
+    """
+
+    #Water levels at centroids
+    wc = domain.quantities['level'].centroid_values
+
+    #Bed elevations at centroids
+    zc = domain.quantities['elevation'].centroid_values
+
+    #Water depths at centroids    
+    hc = wc - zc
+
+    #Momentums at centroids
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values        
+
+    #Water levels at vertices
+    #wv = domain.quantities['level'].vertex_values
+    #
+    #Bed elevations at vertices
+    #zv = domain.quantities['elevation'].vertex_values
+    #
+    #Momentums at vertices
+    #xmomv = domain.quantities['xmomentum'].vertex_values
+    #ymomv = domain.quantities['ymomentum'].vertex_values
+    
+    for k in range(domain.number_of_elements):
+        #Protect against infinitesimal heights and high velocities
+        if hc[k] < domain.minimum_allowed_height:
+            #Control level and height
+            if hc[k] < 0.0:
+                wc[k] = zc[k]; hc[k] = 0.0
+
+            #Control momentum as well!
+            xmomc[k] = ymomc[k] = 0.0
+                
+
+    #Compute second order limiter for each conserved quantity        
+    for name in domain.conserved_quantities:
+        domain.quantities[name].second_order_limiter()
+
+  
+#   //Adjusted water heights at vertices    
+#   hv0 = qv0[0]-zv0;
+#   hv1 = qv1[0]-zv1;
+#   hv2 = qv2[0]-zv2;
+  
+#   hmin = min( min(hv0, hv1), hv2 );
+
+#   // -----------------------------------
+#   //  Second run - Linear combination with
+#   //  first order scheme for shallow depths
+#   // -----------------------------------
+  
+  
+#   //Compute maximal variation in bed elevation
+#   //
+#   //This quantitiy is
+#   // dz = max_i abs(z_i - z_c)
+#   //and it is independent of dimension
+#   //In the 1d case zc = (z0+z1)/2
+#   //In the 2d case zc = (z0+z1+z2)/3
+#   //
+#   z_range = max( max(fabs(zv0-zc), fabs(zv1-zc)), fabs(zv2-zc) );
+  
+#   //Create alpha in [0,1], where alpha==0 means using shallow 
+#   //first order scheme and alpha==1 means using the limited
+#   //second order scheme for the stage w
+#   //If hmin < 0 then alpha=0 and first order method is 
+#   if (z_range>0.0) {
+#     alpha = max( min(2*hmin/z_range, 1.0), 0.0);
+#   } else {
+#     alpha = 1.0;
+#   }
+  
+#   if (alpha < 1) {
+#     //Update stage
+#     //
+#     //(1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
+#     //zvi + hc + alpha*(hvi - hc) 
+#     //
+#     qv0[0] = zv0 + hc + alpha*(hv0-hc);
+#     qv1[0] = zv1 + hc + alpha*(hv1-hc);
+#     qv2[0] = zv2 + hc + alpha*(hv2-hc);
+  
+  	
+#     //Update momentum as a linear combination of 
+#     //qc[j] (shallow) and momentum from gradient limiter (deep).
+#     for (j=1; j<number_of_conserved_quantities; j++) {
+#       qv0[j] = (1-alpha)*qc[j] + alpha*qv0[j];
+#       qv1[j] = (1-alpha)*qc[j] + alpha*qv1[j];
+#       qv2[j] = (1-alpha)*qc[j] + alpha*qv2[j];             
+#     }
+#   }
+
+
+#   //Finally, protect against infinitesimal heights and high speeds
+#   hv0 = qv0[0]-zv0; hv1 = qv1[0]-zv1; hv2 = qv2[0]-zv2;
+#   hmax = max(hv0, max(hv1, hv2));
+#   if (hmax < minimum_allowed_height) {
+  
+#     //Set negative heights to bed elevation        
+#     if (hc < 0.0)  {qc[0] = zc; hc = 0.0;}  
+#     if (hv0 < 0.0) {qv0[0] = zv0; hv0 = 0.0;}         
+#     if (hv1 < 0.0) {qv1[0] = zv1; hv1 = 0.0;}
+#     if (hv2 < 0.0) {qv2[0] = zv2; hv2 = 0.0;}
+
+    
+#     //FIXME: Do we need this desperately??
+#     //FIXME: We definitely do encounter negative heights, 
+#     //but what about momentum?
+    
+#     //Control momentum
+#     //qc[1] = 0.0; qc[2] = 0.0;
+#     //qv0[1] = 0.0; qv0[2] = 0.0;    
+#     //qv1[1] = 0.0; qv1[2] = 0.0;    
+#     //qv2[1] = 0.0; qv2[2] = 0.0;            
+#   }
+# }
+
+
+
+
+def limiter_org(k):
+    """limiter function, specific to the shallow water wave equation.
+       It will ensure that h is always non-negative.
+       This version works on the consecutive data structure
+    """
+
+    #This version should be used to reproduce the bad limiting when levels
+    #are close to bed elevation
+
+    import sys
+    from config import beta
+    from Numeric import ones, Float
+
+    #conserved quantities at centroid
+    qc = Volume.conserved_quantities_centroid[k, :]
+    N = len(qc)
+    
+    zc = Volume.field_values_centroid[k,0]         #bed elevation at centroid
+    wc = qc[0]                                     #Stage at centroid
+    hc = wc - zc                                   #Water depth at centroid
+
+    #conserved quant. at vertices
+    qv0 = Volume.conserved_quantities_vertex0[k,:]
+    qv1 = Volume.conserved_quantities_vertex1[k,:]
+    qv2 = Volume.conserved_quantities_vertex2[k,:]
+
+
+    #Bed elevation at vertices
+    zv0 = Volume.field_values_vertex0[k,0]
+    zv1 = Volume.field_values_vertex1[k,0]
+    zv2 = Volume.field_values_vertex2[k,0]
+
+    #deltas
+    dq0 = qv0-qc
+    dq1 = qv1-qc
+    dq2 = qv2-qc    
+
+            
+    #-------------------------------------
+    hmin = hc
+    qmin = qmax = qc
+    for i in range(3):
+        n = Volume.neighbours[k,i]
+        if n >= 0:
+            qn = Volume.conserved_quantities_centroid[n,:]     
+            zn = Volume.field_values_centroid[n,0]        
+            hn = qn[0]-zn
+
+            hmin = min(hmin, hn)
+            qmin = minimum(qmin, qn)
+            qmax = maximum(qmax, qn)
+
+    hmin = max(hmin,0.0)
+
+    
+    #-----------------------------------
+    # First run - limit on conserved quantities
+    # based on neighbouring heights
+    #-----------------------------------   
+
+    dqmax = qmax - qc 
+    dqmin = qmin - qc
+    phi = ones(len(qc), Float)
+    for j in range(N):
+        #First find the gradient limiter (phi)
+        
+        #Vertex0: 
+        r = 1.0                
+        if dq0[j] > 0: r = dqmax[j]/dq0[j]
+        if dq0[j] < 0: r = dqmin[j]/dq0[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+
+        #Vertex1: 
+        r = 1.0                
+        if dq1[j] > 0: r = dqmax[j]/dq1[j]
+        if dq1[j] < 0: r = dqmin[j]/dq1[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+
+        #Vertex2: 
+        r = 1.0                
+        if dq2[j] > 0: r = dqmax[j]/dq2[j]
+        if dq2[j] < 0: r = dqmin[j]/dq2[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+
+        #The update using phi limiter
+        qv0[j] = qc[j] + phi[j]*dq0[j] #Vertex0
+        qv1[j] = qc[j] + phi[j]*dq1[j] #Vertex1
+        qv2[j] = qc[j] + phi[j]*dq2[j] #Vertex2
+        
+
+        
+    #-----------------------------------
+    # Second run - lower limit on height
+    # based on neighbouring heights for 
+    # surfaces constructed below bed
+    #-----------------------------------
+
+    #Adjusted water heights at vertices    
+    hv0 = qv0[0]-zv0
+    hv1 = qv1[0]-zv1
+    hv2 = qv2[0]-zv2
+
+    #Deltas
+    dhmin = hmin - hc
+    dh0 = hv0 - hc
+    dh1 = hv1 - hc 
+    dh2 = hv2 - hc
+
+
+    #Work out if update is needed and calculate phi limiter
+    phi = 1.0
+    update_needed = False    
+    if hv0 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh0 < 0: r = dhmin/dh0
+        phi = min( min(beta*r, 1), phi )
+        
+    if hv1 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh1 < 0: r = dhmin/dh1
+        phi = min( min(beta*r, 1), phi )
+
+    if hv2 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh2 < 0: r = dhmin/dh2
+        phi = min( min(beta*r, 1), phi )
+
+
+    #Update stage at vertices    
+    if update_needed:
+        hv0 = hc + phi*dh0
+        if hv0 < 0.0: raise 'hv0 = %12f hc = %12f ' % (hv0 , hc)
+        qv0[0] = zv0 + hv0
+
+        hv1 = hc + phi*dh1
+        if hv1 < 0.0: raise 'hv1 = %12f hc = %12f ' % (hv1 , hc)
+        qv1[0] = zv1 + hv1        
+
+        hv2 = hc + phi*dh2
+        if hv2 < 0.0: raise 'hv2 = %12f hc = %12f ' % (hv2 , hc)
+        qv2[0] = zv2 + hv2
+        
+       
+    #In either case protect against infinitesimal heights
+    dry(k)
+
+def dry(k):
+    """Protect agains infinitesimal heights (and high speeds)
+    There are no states for any conserved quantity when h is
+    between minimum allowed height and 0.
+    """
+
+    #volume = Volume.instances[k]
+    
+    maxhv = 0.0
+
+    wv0 = Volume.conserved_quantities_vertex0[k,0]        
+    zv0 = Volume.field_values_vertex0[k,0]
+    hv0 = wv0-zv0
+
+    wv1 = Volume.conserved_quantities_vertex1[k,0]        
+    zv1 = Volume.field_values_vertex1[k,0]
+    hv1 = wv1-zv1
+
+    wv2 = Volume.conserved_quantities_vertex2[k,0]        
+    zv2 = Volume.field_values_vertex2[k,0]
+    hv2 = wv2-zv2    
+        
+    maxhv = max(hv0, hv1, hv2)
+
+    if maxhv < minimum_allowed_height:
+        #Set water level to bed elevation and zero momentums.
+        Volume.conserved_quantities_vertex0[k,:] =\
+                                                 array([Volume.field_values_vertex0[k,0], 0, 0])
+
+        Volume.conserved_quantities_vertex1[k,:] =\
+                                                 array([Volume.field_values_vertex1[k,0], 0, 0])
+
+        Volume.conserved_quantities_vertex2[k,:] =\
+                                                 array([Volume.field_values_vertex2[k,0], 0, 0])
+
+        #Set only momentums to zero
+        #Volume.conserved_quantities_vertex0[k,1:] = array([0., 0.])
+        #Volume.conserved_quantities_vertex1[k,1:] = array([0., 0.])
+        #Volume.conserved_quantities_vertex2[k,1:] = array([0., 0.]) 
+
+
+
+
+def distribute_to_vertices_and_edges(domain):
+    """Extrapolate conserved quantities from centroid to
+    vertices and edge-midpoints for each volume
+    """
+
+    from config import minimum_allowed_height
+
+    from mesh import Point        
+    order = domain.order
+        
+    limiter = domain.limiter
+    if limiter is None:
+        from shallow_water import limit_on_w_then_h as limiter
+
+
+    for k in range(Volume.number_of_instances):
+        
+        #conserved quantities at centroid        
+        q = Volume.conserved_quantities_centroid[k,:]
+
+        if order == 1:
+            ########################
+            #NOTE: This is now specific to the shallow water wave equation
+
+            zc = Volume.field_values_centroid[k,0] #bed elevation at centroid
+            wc = q[0]                              #Stage at centroid
+            hc = wc - zc                           #Water depth at centroid
+
+            #Protect against infinitesimal heights and high velocities
+            if hc < minimum_allowed_height:
+                #Set momentums to zero
+                q[0] = zc
+                q[1] = 0.0
+                q[2] = 0.0	
+            
+
+            #Bed elevation at vertices
+            zv0 = Volume.field_values_vertex0[k,0]
+            zv1 = Volume.field_values_vertex1[k,0]
+            zv2 = Volume.field_values_vertex2[k,0]
+            
+            #Update conserved quant. at vertices
+            Volume.conserved_quantities_vertex0[k,1:] = q[1:] #Momentum
+            Volume.conserved_quantities_vertex1[k,1:] = q[1:] #Momentum
+            Volume.conserved_quantities_vertex2[k,1:] = q[1:] #Momentum
+
+            #Stage
+            z_range = max( max(abs(zv0-zv1), abs(zv0-zv2)), abs(zv1-zv2) )
+            if z_range>0:
+                alpha = min(hc/z_range, 1.0)
+            else:
+                alpha = 1
+               
+            zz = hc+alpha*zc
+            Volume.conserved_quantities_vertex0[k,0] = zz + (1-alpha)*zv0
+            Volume.conserved_quantities_vertex1[k,0] = zz + (1-alpha)*zv1
+            Volume.conserved_quantities_vertex2[k,0] = zz + (1-alpha)*zv2
+
+        elif order == 2:
+            a, b = compute_gradient(k)            
+
+            #Compute extrapolated values at vertices
+            x,y = Point.coordinates[Volume.points[k,3]] #Centroid
+
+            #Get vertex coordinates
+            i = Volume.points[k,0] 
+            x0 = Point.coordinates[i,0]
+            y0 = Point.coordinates[i,1]
+
+            i = Volume.points[k,1] 
+            x1 = Point.coordinates[i,0]
+            y1 = Point.coordinates[i,1]
+
+            i = Volume.points[k,2]             
+            x2 = Point.coordinates[i,0]
+            y2 = Point.coordinates[i,1]                        
+
+
+            #Extrapolate
+            Volume.conserved_quantities_vertex0[k, :] = q + a*(x0-x) + b*(y0-y)
+            Volume.conserved_quantities_vertex1[k, :] = q + a*(x1-x) + b*(y1-y)
+            Volume.conserved_quantities_vertex2[k, :] = q + a*(x2-x) + b*(y2-y)
+            
+
+            if limiter is not None:
+                limiter(k)
+        else:
+            raise 'Unknown order'
+
+
+        #Compute conserved quantities as interpolated
+        #quantities at each edge midpoint
+        q0 = Volume.conserved_quantities_vertex0[k,:]
+        q1 = Volume.conserved_quantities_vertex1[k,:]
+        q2 = Volume.conserved_quantities_vertex2[k,:]
+
+        Volume.conserved_quantities_face0[k,:] = 0.5*(q1 + q2)
+        Volume.conserved_quantities_face1[k,:] = 0.5*(q0 + q2)
+        Volume.conserved_quantities_face2[k,:] = 0.5*(q0 + q1)        
+
+
+###############################################
+#Boundary - specific to the shallow water wave equation
+class Reflective_boundary(Boundary):
+    """Reflective boundary returns same conserved quantities as
+    those present in its neighbour volume but reflected.
+
+    This class is specific to the shallow water equation as it
+    works with the momentum quantities assumed to be the second
+    and third conserved quantities.
+    """
+    
+    def __init__(self, domain = None):        
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for reflective boundary'
+            raise msg
+        
+        self.domain = domain
+        
+        
+    def __repr__(self):
+        return 'Reflective_boundary'
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Reflective boundaries reverses the outward momentum
+        of the volume they serve.
+        """
+        
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        normal = self.domain.get_normal(vol_id, edge_id)
+        
+        r = rotate(q, normal, direction = 1)
+        r[1] = -r[1]
+        q = rotate(r, normal, direction = -1)
+
+        return q
+
+
+
+##############################################
+#Initialise module
+
Index: inundation/ga/storm_surge/pyvolution/test_advection.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_advection.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_advection.py	(revision 195)
@@ -0,0 +1,177 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+from config import g, epsilon
+from Numeric import allclose, array, zeros, ones, Float
+from advection import *
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+        
+    def tearDown(self):
+        pass
+
+    def test_init(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
+        
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+        
+        assert domain.quantities.has_key('level')
+
+        assert domain.get_conserved_quantities(0, edge=1) == 0.
+
+
+    def test_flux_1_triangle0(self):
+        a = [0.0, 0.5]
+        b = [0.0, 0.0]                
+        c = [0.5, 0.5]
+
+        points = [a, b, c]
+        vertices = [ [0,1,2] ]
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+
+ 	assert allclose(domain.neighbours, [[-1,-2,-3]])	
+		
+        #Populate boundary array with dirichlet conditions.
+        domain.quantities['level'].boundary_values[:] = 1.0
+
+	domain.order = 1
+
+	domain.distribute_to_vertices_and_edges() #Use first order default
+        
+        domain.check_integrity()
+
+	domain.compute_fluxes()
+        U = -domain.quantities['level'].explicit_update
+        R = -0.5/domain.areas[0]
+
+        assert U==R, '%s %s' %(U, R) 
+
+
+    def test_flux_1_triangle1(self):
+ 
+        a = [0.0, 0.5]
+        b = [0.0, 0.0]                
+        c = [0.5, 0.5]
+
+        points = [a, b, c]
+        vertices = [ [0,1,2] ]
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+
+        domain.set_quantity('level', [1.0], 'centroids')
+
+	domain.distribute_to_vertices_and_edges()
+        domain.check_integrity()
+        
+        
+	domain.compute_fluxes()
+        U = -domain.quantities['level'].explicit_update
+        R = 0.5/domain.areas[0]
+
+        assert U==R, '%s %s' %(U, R) 
+
+
+
+    def test_flux_1_triangle2(self):
+        
+        a = [0.0, 0.5]
+        b = [0.0, 0.0]                
+        c = [0.5, 0.5]
+
+        points = [a, b, c]
+        vertices = [ [0,1,2] ]
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+
+ 	assert allclose(domain.neighbours, [[-1,-2,-3]])	
+		
+        #Populate boundary array with dirichlet conditions.
+        domain.quantities['level'].boundary_values[0] = 1.0
+        
+	domain.distribute_to_vertices_and_edges() #Use first order default
+        
+        domain.check_integrity()
+
+	domain.compute_fluxes()
+        U = domain.quantities['level'].explicit_update
+        assert allclose(U, 0)
+
+        
+
+
+    def test_flux_2_triangles(self):
+        """Flow between two triangles
+        Check that fluxes have opposite signs
+        """
+
+        a = [0.0, 0.5]
+        b = [0.0, 0.0]                
+        c = [0.5, 0.5]
+        d = [0.5, 0.0]        
+
+        points = [a, b, c, d]
+        vertices = [ [0,1,2], [3,2,1] ]
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+
+ 	assert allclose(domain.neighbours, [[1,-1,-2], [0,-3,-4]])	
+
+        #Populate boundary array with dirichlet conditions.
+
+        domain.set_quantity('level', [1.0, 0.0], 'centroids')
+	domain.distribute_to_vertices_and_edges()
+
+	domain.compute_fluxes()	        
+
+        X = domain.quantities['level'].explicit_update
+        assert X[0] == -X[1]
+			
+
+    def test_advection_example(self):
+        #Test that system can evolve
+        
+        from mesh_factory import rectangular
+
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create advection domain with direction (1,-1)        
+        domain = Domain(points, vertices, boundary, velocity=[1.0, -1.0])
+
+        # Initial condition is zero by default
+        
+        #Boundaries
+        T = Transmissive_boundary(domain)
+        D = Dirichlet_boundary(array([3.1415]))
+        
+        domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
+        domain.check_integrity()
+
+        #Check that the boundary value gets propagated to all elements
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
+            if allclose(domain.quantities['level'].centroid_values, 3.1415):
+                break
+            
+        assert allclose(domain.quantities['level'].centroid_values, 3.1415)
+
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: inundation/ga/storm_surge/pyvolution/test_all.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_all.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_all.py	(revision 195)
@@ -0,0 +1,33 @@
+"""Regression testing framework
+This module will search for scripts in the same directory named
+test_*.py.  Each such script should be a test suite that tests a
+module through PyUnit. This script will aggregate all
+found test suites into one big test suite and run them all at once.
+"""
+
+# Author: Mark Pilgrim
+# Modified by Ole Nielsen
+
+import unittest
+
+def regressionTest():
+    import sys, os, re, unittest    
+    path = os.path.split(sys.argv[0])[0] or os.getcwd()
+    files = os.listdir(path)
+    test = re.compile("^test_[\w]*.py$", re.IGNORECASE)
+    files = filter(test.search, files)
+
+    try:
+        files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
+    except:
+        files.remove('test_all.py') 
+
+    filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    moduleNames = map(filenameToModuleName, files)
+    modules = map(__import__, moduleNames)
+    load = unittest.defaultTestLoader.loadTestsFromModule
+    return unittest.TestSuite(map(load, modules))
+
+if __name__ == "__main__":
+    unittest.main(defaultTest="regressionTest")
+    
Index: inundation/ga/storm_surge/pyvolution/test_domain.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 195)
@@ -0,0 +1,264 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt
+
+from domain import *
+from config import epsilon
+from Numeric import allclose, array
+from util import distance
+
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+
+        
+    def tearDown(self):
+        pass
+
+            
+    def test_simple(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
+
+        conserved_quantities = ['level', 'xmomentum', 'ymomentum']
+        other_quantities = ['elevation', 'friction']
+        
+        domain = Domain(points, vertices, None,
+                        conserved_quantities, other_quantities)        
+        domain.check_integrity()
+
+        for name in conserved_quantities + other_quantities:
+            assert domain.quantities.has_key(name)
+
+
+        assert domain.get_conserved_quantities(0, edge=1) == 0.    
+
+
+    def test_conserved_quantities(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
+
+        domain = Domain(points, vertices, boundary=None,
+                        conserved_quantities =\
+                        ['level', 'xmomentum', 'ymomentum'])
+
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3],
+                                      [0,0,0], [0,0,0]])        
+
+        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
+                                          [0,0,9], [-6, 3, 3],
+                                          [0,0,0], [0,0,0]])        
+
+        domain.check_integrity()
+
+        #Centroids
+        q = domain.get_conserved_quantities(0)
+        assert allclose(q, [2., 2., 0.])
+
+        q = domain.get_conserved_quantities(1)
+        assert allclose(q, [5., 5., 0.])         
+
+        q = domain.get_conserved_quantities(2)
+        assert allclose(q, [3., 3., 0.])
+
+        q = domain.get_conserved_quantities(3)
+        assert allclose(q, [0., 0., 0.])                         
+        
+
+        #Edges
+        q = domain.get_conserved_quantities(0, edge=0)
+        assert allclose(q, [2.5, 2.5, 0.])
+        q = domain.get_conserved_quantities(0, edge=1)
+        assert allclose(q, [2., 2., 0.])
+        q = domain.get_conserved_quantities(0, edge=2)
+        assert allclose(q, [1.5, 1.5, 0.])
+
+        for i in range(3):
+            q = domain.get_conserved_quantities(1, edge=i)
+            assert allclose(q, [5, 5, 0.])
+
+
+        q = domain.get_conserved_quantities(2, edge=0)
+        assert allclose(q, [4.5, 4.5, 0.])
+        q = domain.get_conserved_quantities(2, edge=1)
+        assert allclose(q, [4.5, 4.5, 0.])
+        q = domain.get_conserved_quantities(2, edge=2)
+        assert allclose(q, [0., 0., 0.])
+
+
+        q = domain.get_conserved_quantities(3, edge=0)
+        assert allclose(q, [3., 3., 0.])
+        q = domain.get_conserved_quantities(3, edge=1)
+        assert allclose(q, [-1.5, -1.5, 0.])
+        q = domain.get_conserved_quantities(3, edge=2)
+        assert allclose(q, [-1.5, -1.5, 0.])                    
+
+
+
+    def test_boundary_conditions(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'First',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Second'}                                         
+                     
+
+        domain = Domain(points, vertices, boundary,
+                        conserved_quantities =\
+                        ['level', 'xmomentum', 'ymomentum'])        
+        domain.check_integrity()
+
+
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3]])
+
+
+        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]),
+                              'Second': Transmissive_boundary(domain)} )
+
+        domain.update_boundary()
+
+        assert domain.quantities['level'].boundary_values[0] == 5. #Dirichlet
+        assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
+        assert domain.quantities['level'].boundary_values[2] ==\
+               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
+        assert domain.quantities['level'].boundary_values[3] ==\
+               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
+        assert domain.quantities['level'].boundary_values[4] ==\
+               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
+        assert domain.quantities['level'].boundary_values[5] ==\
+               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)         
+
+
+
+    def test_distribute(self):
+        """Domain implements a default first order gradient limiter
+        """
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'Third',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Third'}                                         
+                     
+
+        domain = Domain(points, vertices, boundary,
+                        conserved_quantities =\
+                        ['level', 'xmomentum', 'ymomentum'])                
+        domain.check_integrity()
+
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3]])
+
+        assert allclose( domain.quantities['level'].centroid_values,
+                         [2,5,3,0] )
+                         
+        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
+                                          [3,3,3], [4, 4, 4]])
+
+        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
+                                          [30,30,30], [40, 40, 40]])        
+
+
+        #print domain.quantities['level'].centroid_values
+        domain.distribute_to_vertices_and_edges()
+
+        #First order extrapolation
+        assert allclose( domain.quantities['level'].vertex_values,
+                         [[ 2.,  2.,  2.],
+                          [ 5.,  5.,  5.],
+                          [ 3.,  3.,  3.],
+                          [ 0.,  0.,  0.]])
+
+
+
+
+    def test_update_conserved_quantities(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'Third',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Third'}                                         
+                     
+
+        domain = Domain(points, vertices, boundary,
+                        conserved_quantities =\
+                        ['level', 'xmomentum', 'ymomentum'])                
+        domain.check_integrity()
+
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3]])
+
+        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
+                                          [3,3,3], [4, 4, 4]])
+
+        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
+                                          [30,30,30], [40, 40, 40]])        
+
+
+        #Assign some values to update vectors
+        #domain.explicit_updates[]
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py	(revision 195)
@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+from generic_boundary_conditions import *
+from config import epsilon
+from Numeric import allclose, array
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+        #print "  Setting up"
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_generic(self):
+        b = Boundary()
+
+        try:
+            b.evaluate()
+        except:
+            pass
+        else:
+            raise 'Should have raised exception'
+        
+
+    def test_dirichlet_empty(self):
+
+        try:
+            Bd = Dirichlet_boundary()
+        except:
+            pass
+        else:
+            raise 'Should have raised exception'
+
+    def test_dirichlet(self):
+        x = [3.14,0,0.1]
+        Bd = Dirichlet_boundary(x)
+
+        q = Bd.evaluate()
+        assert allclose(q, x)
+
+
+    def test_transmissive(self):
+        from domain import Domain
+        from quantity import Quantity
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        
+        domain = Domain(points, elements)
+        domain.check_integrity()
+
+        
+        domain.quantities['level'] =\
+                                   Quantity(domain, [[1,2,3], [5,5,5],
+                                                     [0,0,9], [-6, 3, 3]])
+
+        domain.quantities['ymomentum'] =\
+                                   Quantity(domain, [[1,2,3], [5,5,5],
+                                                     [0,0,9], [-6, 3, 3]])
+
+
+        domain.check_integrity()
+
+        #Test transmissve bdry
+        try:
+            T = Transmissive_boundary()
+        except:
+            pass
+        else:
+            raise 'Should have raised exception'
+
+        T = Transmissive_boundary(domain)
+
+        from config import default_boundary_tag
+        domain.set_boundary( {default_boundary_tag: T} )
+
+        assert len(domain.boundary) == len(domain.boundary_segments)
+        assert len(domain.boundary_objects) == len(domain.boundary_segments)        
+
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
+
+        
Index: inundation/ga/storm_surge/pyvolution/test_mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 195)
@@ -0,0 +1,602 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt
+
+from mesh import *
+from config import epsilon
+from Numeric import allclose, array
+from util import distance
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+        #print "  Setting up"
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_triangle_inputs(self):
+        points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
+        vertices = [0,1,2] #Wrong
+
+        try:
+            mesh = Mesh(points, vertices)
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+        
+        
+    def test_basic_triangle(self):
+
+        a = [0.0, 0.0]
+        b = [4.0, 0.0]
+        c = [0.0, 3.0]
+        
+        points = [a, b, c]
+        vertices = [[0,1,2]]
+        mesh = Mesh(points, vertices) 
+
+        #Centroid
+        centroid = mesh.centroids[0]
+        assert centroid[0] == 4.0/3
+        assert centroid[1] == 1.0
+
+        #Area
+        assert mesh.areas[0] == 6.0,\
+               'Area was %f, should have been 6.0' %triangle.get_area()        
+
+        #Normals
+        normals = mesh.get_normals()
+        assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
+        assert allclose(normals[0, 2:4], [-1.0, 0.0])
+        assert allclose(normals[0, 4:6], [0.0, -1.0])        
+
+        assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
+        assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
+        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])                
+        
+        #Edge lengths
+        assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
+
+
+        #Vertex coordinates
+        V = mesh.get_vertex_coordinates()
+        assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
+
+        V0 = mesh.get_vertex_coordinate(0, 0)
+        assert allclose(V0, [0.0, 0.0])
+
+        V1 = mesh.get_vertex_coordinate(0, 1)
+        assert allclose(V1, [4.0, 0.0])
+
+        V2 = mesh.get_vertex_coordinate(0, 2)
+        assert allclose(V2, [0.0, 3.0])                
+
+
+        #General tests:
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()
+
+        
+        #Test that the centroid is located 2/3 of the way
+        #from each vertex to the midpoint of the opposite side
+
+        V = mesh.get_vertex_coordinates()
+        
+        x0 = V[0,0]
+        y0 = V[0,1]
+        x1 = V[0,2]
+        y1 = V[0,3]
+        x2 = V[0,4]
+        y2 = V[0,5]
+        
+        m0 = [(x1 + x2)/2, (y1 + y2)/2]
+        m1 = [(x0 + x2)/2, (y0 + y2)/2]
+        m2 = [(x1 + x0)/2, (y1 + y0)/2]
+
+        d0 = distance(centroid, [x0, y0])
+        d1 = distance(m0, [x0, y0])
+        assert d0 == 2*d1/3
+        #
+        d0 = distance(centroid, [x1, y1])
+        d1 = distance(m1, [x1, y1])        
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
+
+        d0 = distance(centroid, [x2, y2])
+        d1 = distance(m2, [x2, y2])
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)        
+        
+        #Radius
+        d0 = distance(centroid, m0)
+        assert d0 == 5.0/6
+        
+        d1 = distance(centroid, m1)
+        assert d1 == sqrt(73.0/36)        
+
+        d2 = distance(centroid, m2)        
+        assert d2 == sqrt(13.0/9)        
+        
+        assert mesh.radii[0] == min(d0, d1, d2)
+        assert mesh.radii[0] == 5.0/6
+
+
+
+        #Let x be the centroid of triangle abc.
+        #Test that areas of the three triangles axc, cxb, and bxa are equal.
+
+        points = [a, b, c, centroid]
+        vertices = [[0,3,2], [2,3,1], [1,3,0]]
+        new_mesh = Mesh(points, vertices)
+        
+        assert new_mesh.areas[0] == new_mesh.areas[1]
+        assert new_mesh.areas[1] == new_mesh.areas[2]        
+        assert new_mesh.areas[1] == new_mesh.areas[2]        
+
+        assert new_mesh.areas[1] == mesh.areas[0]/3        
+
+
+
+    def test_general_triangle(self):
+        a = [2.0, 1.0]
+        b = [6.0, 2.0]
+        c = [1.0, 3.0]
+
+        points = [a, b, c]
+        vertices = [[0,1,2]]
+
+        mesh = Mesh(points, vertices) 
+        centroid = mesh.centroids[0]
+
+        
+        #Test that the centroid is located 2/3 of the way
+        #from each vertex to the midpoint of the opposite side
+
+        V = mesh.get_vertex_coordinates()
+        
+        x0 = V[0,0]
+        y0 = V[0,1]
+        x1 = V[0,2]
+        y1 = V[0,3]
+        x2 = V[0,4]
+        y2 = V[0,5]
+        
+        m0 = [(x1 + x2)/2, (y1 + y2)/2]
+        m1 = [(x0 + x2)/2, (y0 + y2)/2]
+        m2 = [(x1 + x0)/2, (y1 + y0)/2]
+
+        d0 = distance(centroid, [x0, y0])
+        d1 = distance(m0, [x0, y0])
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
+        #
+        d0 = distance(centroid, [x1, y1])
+        d1 = distance(m1, [x1, y1])        
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
+
+        d0 = distance(centroid, [x2, y2])
+        d1 = distance(m2, [x2, y2])
+        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)        
+        
+        #Radius
+        d0 = distance(centroid, m0)
+        d1 = distance(centroid, m1)
+        d2 = distance(centroid, m2)        
+        assert mesh.radii[0] == min(d0, d1, d2)
+
+
+
+        #Let x be the centroid of triangle abc.
+        #Test that areas of the three triangles axc, cxb, and bxa are equal.
+
+        points = [a, b, c, centroid]
+        vertices = [[0,3,2], [2,3,1], [1,3,0]]
+        new_mesh = Mesh(points, vertices)
+        
+        assert new_mesh.areas[0] == new_mesh.areas[1]
+        assert new_mesh.areas[1] == new_mesh.areas[2]        
+        assert new_mesh.areas[1] == new_mesh.areas[2]        
+
+        assert new_mesh.areas[1] == mesh.areas[0]/3        
+        
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()
+
+
+    def test_boundary_indices(self):
+        a = [0.0, 0.5]
+        b = [0.0, 0.0]                
+        c = [0.5, 0.5]
+
+        points = [a, b, c]
+        vertices = [ [0,1,2] ]
+        mesh = Mesh(points, vertices)        
+        mesh.check_integrity()
+
+ 	assert allclose(mesh.neighbours, [[-1,-2,-3]])	
+		
+
+
+    def test_two_triangles(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        e = [2.0, 2.0]
+        points = [a, b, c, e]
+        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        mesh = Mesh(points, vertices)
+        
+        assert mesh.areas[0] == 2.0
+
+        assert allclose(mesh.centroids[0], [2.0/3, 2.0/3])        
+
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()        
+
+
+            
+    def test_more_triangles(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
+        mesh = Mesh(points, vertices)        
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()        
+
+        assert mesh.areas[0] == 2.0
+        assert mesh.areas[1] == 2.0
+        assert mesh.areas[2] == 2.0
+        assert mesh.areas[3] == 2.0
+        assert mesh.areas[4] == 8.0
+        assert mesh.areas[5] == 4.0                                
+        
+        assert mesh.edgelengths[1,0] == 2.0
+        assert mesh.edgelengths[1,1] == 2.0
+        assert mesh.edgelengths[1,2] == sqrt(8.0)
+
+        assert allclose(mesh.centroids[0], [2.0/3, 2.0/3])
+        assert allclose(mesh.centroids[1], [4.0/3, 4.0/3]) 
+        assert allclose(mesh.centroids[2], [8.0/3, 2.0/3])
+        assert allclose(mesh.centroids[3], [2.0/3, 8.0/3])         
+
+
+    def test_mesh_and_neighbours(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        mesh = Mesh(points, vertices)
+        
+        mesh.check_integrity()
+
+
+        T = mesh
+        tid = 0
+        assert T.number_of_boundaries[tid] == 2    
+        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
+        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
+        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)        
+
+        tid = 1
+        assert T.number_of_boundaries[tid] == 0            
+        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
+        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
+        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
+
+        tid = 2
+        assert T.number_of_boundaries[tid] == 2            
+        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
+        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
+        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)        
+
+        tid = 3
+        assert T.number_of_boundaries[tid] == 2            
+        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
+        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
+        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
+
+        #Neighbouring edges
+        tid = 0
+        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
+        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
+        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)        
+
+        tid = 1        
+        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
+        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
+        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
+
+        tid = 2        
+        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
+        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
+        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)        
+
+        tid = 3        
+        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
+        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
+        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
+
+        
+    def test_rectangular_mesh_basic(self):
+        M=1
+        N=1
+        mesh = rectangular_mesh(M, N)
+
+        #Test that points are arranged in a counter clock wise order        
+        mesh.check_integrity()
+
+        M=2
+        N=2
+        mesh = rectangular_mesh(M, N)
+        
+        #Test that points are arranged in a counter clock wise order        
+        mesh.check_integrity()
+
+        #assert mesh.boundary[(7,1)] == 2 # top
+        assert mesh.boundary[(7,1)] == 'top' # top
+        assert mesh.boundary[(3,1)] == 'top' # top                
+
+        
+    def test_boundary_tags(self):
+        
+
+        mesh = rectangular_mesh(4, 4)
+
+        #Test that points are arranged in a counter clock wise order        
+        mesh.check_integrity()
+
+        #print mesh.get_boundary_tags()
+        #print mesh.boundary
+
+        for k in [1, 3, 5, 7]:
+            assert mesh.boundary[(k,2)] == 'left'
+
+        for k in [24, 26, 28, 30]:
+            assert mesh.boundary[(k,2)] == 'right'            
+            
+        for k in [7, 15, 23, 31]:
+            assert mesh.boundary[(k,1)] == 'top'
+        for k in [0, 8, 16, 24]:
+            assert mesh.boundary[(k,1)] == 'bottom'
+
+
+        
+    def test_rectangular_mesh(self):
+        M=4
+        N=16
+        len1 = 100.0
+        len2 = 17.0
+        mesh = rectangular_mesh(M, N, len1, len2)
+
+        assert len(mesh) == 2*M*N        
+
+        for i in range(len(mesh)):
+            assert mesh.areas[i] == len1*len2/(2*M*N)
+
+            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
+            assert mesh.edgelengths[i, 0] == hypo
+            assert mesh.edgelengths[i, 1] == len1/M #x direction
+            assert mesh.edgelengths[i, 2] == len2/N #y direction
+            
+        #Test that points are arranged in a counter clock wise order      
+        mesh.check_integrity()    
+
+
+    def test_rectangular_mesh2(self):
+        #Check that integers don't cause trouble
+        N = 16
+        mesh = rectangular_mesh(2*N, N, len1=10, len2=10)
+
+
+    def test_surrogate_neighbours(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        mesh = Mesh(points, vertices)        
+        mesh.check_integrity()
+
+
+        T = mesh
+        tid = 0
+        assert T.number_of_boundaries[tid] == 2    
+        assert T.surrogate_neighbours[tid, 0] == tid 
+        assert T.surrogate_neighbours[tid, 1] == 1   
+        assert T.surrogate_neighbours[tid, 2] == tid 
+
+        tid = 1
+        assert T.number_of_boundaries[tid] == 0            
+        assert T.surrogate_neighbours[tid, 0] == 2
+        assert T.surrogate_neighbours[tid, 1] == 3
+        assert T.surrogate_neighbours[tid, 2] == 0
+
+        tid = 2
+        assert T.number_of_boundaries[tid] == 2            
+        assert T.surrogate_neighbours[tid, 0] == tid 
+        assert T.surrogate_neighbours[tid, 1] == tid 
+        assert T.surrogate_neighbours[tid, 2] == 1
+
+        tid = 3
+        assert T.number_of_boundaries[tid] == 2            
+        assert T.surrogate_neighbours[tid, 0] == 1 
+        assert T.surrogate_neighbours[tid, 1] == tid
+        assert T.surrogate_neighbours[tid, 2] == tid
+
+
+    def test_boundary_inputs(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        boundary = { (0, 0): 'First',
+                     (0, 2): 'Second',
+                     (2, 0): 'Third',
+                     (2, 1): 'Fourth',
+                     (3, 1): 'Fifth',
+                     (3, 2): 'Sixth'}                                         
+                     
+        
+        mesh = Mesh(points, vertices, boundary)        
+        mesh.check_integrity()
+
+
+        #Check enumeration
+        for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):        
+            b = -k-1
+            assert mesh.neighbours[vol_id, edge_id] == b
+
+
+
+    def test_boundary_inputs_using_defaults(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        boundary = { (0, 0): 'First',
+                     (0, 2): 'Second',
+                     (2, 0): 'Third',
+                     (2, 1): 'Fourth',
+                     #(3, 1): 'Fifth',  #Skip this
+                     (3, 2): 'Sixth'}                                         
+                     
+        
+        mesh = Mesh(points, vertices, boundary)        
+        mesh.check_integrity()
+
+        from config import default_boundary_tag
+        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
+
+
+        #Check enumeration
+        for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):        
+            b = -k-1
+            assert mesh.neighbours[vol_id, edge_id] == b
+
+
+
+
+
+
+    def test_inputs(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        #Too few points
+        try:
+            mesh = Mesh([points[0]], vertices)        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
+        #Too few points - 1 element
+        try:
+            mesh = Mesh([points[0]], [vertices[0]])        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'        
+        
+        #Wrong dimension of vertices
+        try:
+            mesh = Mesh(points, vertices[0])        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
+        #Unsubscriptable coordinates object raises exception
+        try:
+            mesh = Mesh(points[0], [vertices[0]])        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
+        #Not specifying all boundary tags
+        try:
+            mesh = Mesh(points, vertices, {(3,0): 'x'})        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
+        #Specifying wrong non existing segment
+        try:
+            mesh = Mesh(points, vertices, {(5,0): 'x'})        
+        except AssertionError:
+            pass
+        else:
+            raise 'Should have raised an exception'
+
+
+
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
+    
+    
+
+
+
Index: inundation/ga/storm_surge/pyvolution/test_quantity.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 195)
@@ -0,0 +1,260 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from quantity import *
+from config import epsilon
+from Numeric import allclose, array
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        from domain import Domain
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        self.mesh1 = Domain(points, [elements[0]])        
+        self.mesh1.check_integrity()
+        
+        self.mesh4 = Domain(points, elements)        
+        self.mesh4.check_integrity()
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        
+        quantity = Quantity(self.mesh1, [[1,2,3]])
+        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
+
+        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 mesh object" error'        
+        
+
+    def test_creation_zeros(self):
+       
+        quantity = Quantity(self.mesh1)
+        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+
+        quantity = Quantity(self.mesh4)
+        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
+                                                 [0.,0.,0.], [0.,0.,0.]])        
+
+        
+    def test_interpolation(self):
+        quantity = Quantity(self.mesh1, [[1,2,3]])
+        assert allclose(quantity.centroid_values, [2.0]) #Centroid
+
+        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 
+
+
+    def test_interpolation2(self):
+        quantity = Quantity(self.mesh4,
+                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+
+        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+                                               [5., 5., 5.],
+                                               [4.5, 4.5, 0.],
+                                               [3.0, -1.5, -1.5]])
+
+    def test_boundary_allocation(self):
+        quantity = Quantity(self.mesh4,
+                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+
+        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
+
+
+    def test_set_values(self):
+        quantity = Quantity(self.mesh4)
+
+
+        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
+                            location = 'vertices')
+        assert allclose(quantity.vertex_values,
+                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+                                               [5., 5., 5.],
+                                               [4.5, 4.5, 0.],
+                                               [3.0, -1.5, -1.5]])
+
+
+        #Test default
+        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert allclose(quantity.vertex_values,
+                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])        
+        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+                                               [5., 5., 5.],
+                                               [4.5, 4.5, 0.],
+                                               [3.0, -1.5, -1.5]])
+
+        #Test centroids
+        quantity.set_values([1,2,3,4], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
+
+        #Test edges
+        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
+                            location = 'edges')
+        assert allclose(quantity.edge_values,
+                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])        
+
+        #Test exceptions
+        try:
+            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
+                                location = 'bas kamel tuba')
+        except:
+            pass
+
+
+        try:
+            quantity.set_values([[1,2,3], [0,0,9]])
+        except AssertionError:
+            pass
+        except:
+            raise 'should have raised Assertionerror'
+
+
+
+    def test_first_order_limiter(self):
+        quantity = Quantity(self.mesh4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+        #Extrapolate
+        quantity.first_order_limiter()
+
+        #Check vertices but not edge values
+        assert allclose(quantity.vertex_values,
+                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+
+
+    def test_second_order_limiter(self):
+
+        #FIXME:!!!!!!!
+        
+        quantity = Quantity(self.mesh4)
+
+        quantity.set_values([2., 4., 5., 5.], location = 'centroids')
+        #Assume gradient
+
+        #Extrapolate
+        quantity.second_order_limiter()
+        #print quantity.vertex_values
+
+
+        
+        #Check vertices but not edge values
+        #assert allclose(quantity.vertex_values,
+        #                [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+
+
+    def test_distribute_first_order(self):
+        quantity = Quantity(self.mesh4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+
+        #Extrapolate
+        quantity.distribute_to_vertices_and_edges(order=1)
+
+        assert allclose(quantity.vertex_values,
+                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
+                                               [3,3,3], [4, 4, 4]])
+
+
+
+    def test_update_explicit(self):
+        quantity = Quantity(self.mesh4)
+
+        #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.mesh4)
+
+        #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
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_both_updates(self):
+        quantity = Quantity(self.mesh4)
+
+        #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( [4.,3.,2.,1.] )
+        
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
+        x /= array( [.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)        
+
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
+
Index: inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 195)
@@ -0,0 +1,624 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+from config import g, epsilon
+from Numeric import allclose, array, zeros, ones, Float
+from shallow_water import *
+
+        
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+        #print "  Setting up"
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_flux_zero_case(self):
+        ql = zeros( 3, Float )
+        qr = zeros( 3, Float )
+        normal = zeros( 3, Float )
+        zl = zr = 0.
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [0,0,0])
+        assert max_speed == 0.
+
+    def test_flux_constants(self):
+        w = 2.0
+       
+        normal = array([1.,0])        
+        ql = array([w, 0, 0])
+        qr = array([w, 0, 0])
+        zl = zr = 0.
+        h = w - (zl+zr)/2
+        
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [0., 0.5*g*h**2, 0.])
+        assert max_speed == sqrt(g*h)
+        
+
+    def test_flux1(self):
+        #Use data from previous version of pyvolution
+        normal = array([1.,0])        
+        ql = array([-0.2, 2, 3])
+        qr = array([-0.2, 2, 3])
+        zl = zr = -0.5
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [2.,13.77433333, 20.])
+        assert allclose(max_speed, 8.38130948661)
+
+
+    def test_flux2(self):
+        #Use data from previous version of pyvolution
+        normal = array([0., -1.])        
+        ql = array([-0.075, 2, 3])
+        qr = array([-0.075, 2, 3])
+        zl = zr = -0.375
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [-3.,-20.0, -30.441])
+        assert allclose(max_speed, 11.7146428199)
+
+    def test_flux3(self):
+        #Use data from previous version of pyvolution        
+        normal = array([-sqrt(2)/2, sqrt(2)/2])        
+        ql = array([-0.075, 2, 3])
+        qr = array([-0.075, 2, 3])
+        zl = zr = -0.375
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
+        assert allclose(max_speed, 4.0716654239)
+
+    def test_flux4(self):
+        #Use data from previous version of pyvolution                
+        normal = array([-sqrt(2)/2, sqrt(2)/2])        
+        ql = array([-0.34319278, 0.10254161, 0.07273855])
+        qr = array([-0.30683287, 0.1071986, 0.05930515])
+        zl = zr = -0.375
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
+        assert allclose(max_speed, 1.31414103233)                                
+
+    def test_sw_domain_simple(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
+        
+        domain = Domain(points, vertices)        
+        domain.check_integrity()
+
+        for name in ['level', 'xmomentum', 'ymomentum',
+                     'elevation', 'friction']:
+            assert domain.quantities.has_key(name)
+
+
+        assert domain.get_conserved_quantities(0, edge=1) == 0.    
+
+
+
+    def test_boundary_conditions(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'Third',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Third'}                                         
+                     
+
+        domain = Domain(points, vertices, boundary)
+        domain.check_integrity()
+
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3]])
+
+        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
+                                          [3,3,3], [4, 4, 4]])
+
+        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
+                                          [30,30,30], [40, 40, 40]])        
+
+
+        D = Dirichlet_boundary([5,2,1])
+        T = Transmissive_boundary(domain)
+        R = Reflective_boundary(domain)
+        domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
+        
+        domain.update_boundary()
+
+        #Level
+        assert domain.quantities['level'].boundary_values[0] == 2.5
+        assert domain.quantities['level'].boundary_values[0] ==\
+               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5) 
+        assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
+        assert domain.quantities['level'].boundary_values[2] ==\
+               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
+        assert domain.quantities['level'].boundary_values[3] ==\
+               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
+        assert domain.quantities['level'].boundary_values[4] ==\
+               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
+        assert domain.quantities['level'].boundary_values[5] ==\
+               domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
+
+        #Xmomentum
+        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
+        assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
+        assert domain.quantities['xmomentum'].boundary_values[2] ==\
+               domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
+        assert domain.quantities['xmomentum'].boundary_values[3] ==\
+               domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
+        assert domain.quantities['xmomentum'].boundary_values[4] ==\
+               domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
+        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
+        
+        #Ymomentum
+        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
+        assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
+        assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
+        assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
+        assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
+        assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
+
+
+        
+
+    def test_compute_fluxes(self):
+        #Do a full triangle and check that fluxes cancel out for
+        #the constant level case
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)        
+        domain.set_quantity('level', [[2,2,2], [2,2,2],
+                                      [2,2,2], [2,2,2]])
+        domain.check_integrity()
+
+        zl=zr=0. #Assume flat bed
+        
+        #Flux across right edge of volume 1
+        normal = domain.get_normal(1,0)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+        #Check that flux seen from other triangles is inverse
+        tmp = qr; qr=ql; ql=tmp
+        normal = domain.get_normal(2,2)
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        assert allclose(flux + flux0, 0.)
+    
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1,1)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+        #Check that flux seen from other triangles is inverse
+        tmp = qr; qr=ql; ql=tmp
+        normal = domain.get_normal(3,0)
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        assert allclose(flux + flux1, 0.)
+        
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1,2)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+        #Check that flux seen from other triangles is inverse
+        tmp = qr; qr=ql; ql=tmp
+        normal = domain.get_normal(0,1)
+        flux, max_speed = flux_function(normal, ql, qr, zl, zr)        
+        assert allclose(flux + flux2, 0.)
+
+
+        #Scale by edgelengths, add up anc check that total flux is zero
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+        
+        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
+
+        #Now check that compute_flux yields zeros as well
+        domain.compute_fluxes()
+
+        for name in ['level', 'xmomentum', 'ymomentum']:
+            #print name, domain.quantities[name].explicit_update
+            assert allclose(domain.quantities[name].explicit_update[1], 0)
+
+
+
+    def test_compute_fluxes1(self):
+        #Use values from previous version
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+        
+        domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1],
+                                      [val2, val2, val2], [val3, val3, val3]])
+        domain.check_integrity()
+
+        zl=zr=0. #Assume flat bed
+
+        #Flux across right edge of volume 1
+        normal = domain.get_normal(1,0)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
+        assert allclose(max_speed, 9.21592824046)
+
+
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1,1)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+        assert allclose(flux1, [2.4098563, 0., 123.04444444])
+        assert allclose(max_speed, 7.22956891292)
+
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1,2)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
+        assert allclose(max_speed, 7.22956891292)
+        
+        #Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+        assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
+
+
+        domain.compute_fluxes()
+
+        #assert allclose(total_flux, domain.explicit_update[1,:])
+        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
+            assert allclose(total_flux[i],
+                            domain.quantities[name].explicit_update[1])
+
+        #assert allclose(domain.explicit_update, [
+        #    [0., -69.68888889, -69.68888889],
+        #    [-0.68218178, -166.6, -35.93333333],
+        #    [-111.77316251, 69.68888889, 0.],
+        #    [-35.68522449, 0., 69.68888889]])
+
+        assert allclose(domain.quantities['level'].explicit_update,
+                        [0., -0.68218178, -111.77316251, -35.68522449])
+        assert allclose(domain.quantities['xmomentum'].explicit_update,
+                        [-69.68888889, -166.6, 69.68888889, 0])
+        assert allclose(domain.quantities['ymomentum'].explicit_update,
+                        [-69.68888889, -35.93333333, 0., 69.68888889])        
+
+        
+        #assert allclose(domain.quantities[name].explicit_update
+
+        
+        
+
+
+    def test_compute_fluxes2(self):
+        #Random values, incl momentum
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=0 #Assume flat zero bed
+
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        
+        
+        domain.set_quantity('level', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1, 2, 3], [3, 4, 5],
+                             [1, -1, 0], [0, -2, 2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1, -1, 0], [0, -3, 2],
+                             [0, 1, 0], [-1, 2, 2]])                
+
+        
+        domain.check_integrity()
+
+
+
+        #Flux across right edge of volume 1
+        normal = domain.get_normal(1,0)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1,1)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1,2)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+        #Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+
+        
+        domain.compute_fluxes()
+        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
+            assert allclose(total_flux[i],
+                            domain.quantities[name].explicit_update[1])        
+        #assert allclose(total_flux, domain.explicit_update[1,:])
+        
+
+    def test_compute_fluxes3(self):
+        #Random values, incl momentum
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=-3.75 #Assume constant bed (must be less than level)
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        
+        
+        domain.set_quantity('level', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1, 2, 3], [3, 4, 5],
+                             [1, -1, 0], [0, -2, 2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1, -1, 0], [0, -3, 2],
+                             [0, 1, 0], [-1, 2, 2]])                
+
+        
+        domain.check_integrity()
+
+
+
+        #Flux across right edge of volume 1
+        normal = domain.get_normal(1,0)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1,1)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1,2)        
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+        #Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+        
+        domain.compute_fluxes()
+        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
+            assert allclose(total_flux[i],
+                            domain.quantities[name].explicit_update[1])
+
+
+
+    def test_catching_negative_heights(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=4  #Too large
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        domain.set_quantity('level', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        #Should fail
+        try:
+            domain.check_integrity()
+        except:
+            pass
+
+
+
+    def test_first_order_limiter_const_z(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=-3.75 #Assume constant bed (must be less than level)        
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        domain.set_quantity('level', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+
+
+        domain.first_order_limiter()
+
+        #Check that centroid values were distributed to vertices
+        C = domain.quantities['level'].centroid_values
+        for i in range(3): 
+            assert allclose( domain.quantities['level'].vertex_values[:,i], C)
+        
+
+    def test_first_order_limiter_variable_z(self):
+        #Check that first order limiter follows bed_slope
+        from Numeric import alltrue, greater_equal
+        from config import epsilon
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
+                                          [6,6,6], [6,6,6]])
+        domain.set_quantity('level', [[val0, val0, val0],
+                                      [val1, val1, val1],
+                                      [val2, val2, val2],
+                                      [val3, val3, val3]])
+
+        L = domain.quantities['level'].vertex_values        
+        E = domain.quantities['elevation'].vertex_values
+
+        
+        #Check that some levels are not above elevation (within eps)
+        #- so that the limiter has something to work with
+        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))                
+                            
+        domain.first_order_limiter()
+
+        #Check that all levels are above elevation (within eps) 
+        assert alltrue(alltrue(greater_equal(L,E-epsilon)))                
+        
+
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
Index: inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_util.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/test_util.py	(revision 195)
@@ -0,0 +1,186 @@
+#!/usr/bin/env python
+#
+
+import unittest
+from Numeric import zeros, array, allclose
+from math import sqrt, pi
+
+from util import *
+from config import epsilon
+
+class TestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+
+    def test_gradient(self):
+        x0 = 0.0; y0 = 0.0; z0 = 0.0
+        x1 = 1.0; y1 = 0.0; z1 = -1.0
+        x2 = 0.0; y2 = 1.0; z2 = 0.0
+        
+        zx, zy = gradient(x0, y0, x1, y1, x2, y2, array([z0]), array([z1]), array([z2]))
+
+        assert zx == -1.0
+        assert zy == 0.0
+
+
+    def test_rotate(self):
+        normal = array([0.0,-1.0])        
+
+        q = array([1.0,2.0,3.0])
+      
+        r = rotate(q, normal, direction = 1)
+        assert r[0] == 1
+        assert r[1] == -3
+        assert r[2] == 2
+
+        w = rotate(r, normal, direction = -1)        
+        assert allclose(w, q)
+
+
+
+    def test_gradient2(self):
+
+        from util import gradient
+        x0 = 2.0/3; y0 = 2.0/3
+        x1=  8.0/3; y1 = 2.0/3
+        x2 = 2.0/3; y2 = 8.0/3
+
+        q0 = array([2.0+2.0/3])
+        q1 = array([8.0+2.0/3])
+        q2 = array([2.0+8.0/3])
+        
+        #Gradient of fitted pwl surface    
+        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        
+        assert abs(a[0] - 3.0) < epsilon
+        assert abs(b[0] - 1.0) < epsilon
+    
+
+
+    def test_that_C_extension_compiles(self):
+        FN = 'util_ext.c'
+        try:
+            import util_ext
+        except:
+            from pytools import compile
+
+            try:
+                compile(FN)
+            except:
+                raise 'Could not compile %s' %FN
+            else:
+                import util_ext
+
+            
+    def test_gradient_C_extension(self):
+        from util_ext import gradient as gradient_c
+
+        x0 = 2.0/3; y0 = 2.0/3
+        x1=  8.0/3; y1 = 2.0/3
+        x2 = 2.0/3; y2 = 8.0/3
+
+        q0 = array([2.0+2.0/3])
+        q1 = array([8.0+2.0/3])
+        q2 = array([2.0+8.0/3])
+        
+        #Gradient of fitted pwl surface    
+        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        
+        assert abs(a[0] - 3.0) < epsilon
+        assert abs(b[0] - 1.0) < epsilon
+
+    def test_gradient_C_extension_scalars(self):
+        from util_ext import gradient as gradient_c
+        
+        from RandomArray import uniform, seed
+        seed(17, 53)
+
+        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
+
+        q0 = 10.0
+        q1 = 3.0
+        q2 = 7.0
+        
+        
+        #Gradient of fitted pwl surface
+	from util import gradient_python	
+        a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+
+        #print a_ref, b_ref
+        try:
+            a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        except TypeError:
+            pass
+        else:
+            raise 'gradient c-extension should throw exception when q''s are scalars'
+
+        
+    def test_gradient_C_extension3(self):
+        from util_ext import gradient as gradient_c
+        
+        from RandomArray import uniform, seed
+        seed(17, 53)
+
+        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
+
+        q0 = uniform(0.0, 10.0, 4)
+        q1 = uniform(1.0, 3.0, 4)
+        q2 = uniform(7.0, 20.0, 4)        
+        
+        
+        #Gradient of fitted pwl surface
+	from util import gradient_python	        
+        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+
+        #print a_ref, b_ref
+        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+
+        #print a, a_ref, b, b_ref
+        for i in range(4):
+            assert abs(a[i] - a_ref[i]) < epsilon
+            assert abs(b[i] - b_ref[i]) < epsilon
+        
+
+
+    #Geometric        
+    def test_distance(self):
+        from util import distance
+
+        self.failUnless( distance([4,2],[7,6]) == 5.0,
+                        'Distance is wrong!')
+        self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
+                        'distance is wrong!')
+        self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
+                        'distance is wrong!')
+
+        self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
+                        'distance is wrong!')        
+
+    def test_angle(self):
+        from util import angle
+
+        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
+
+
+    def test_anglediff(self):
+        from util import anglediff        
+
+        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
+
+
+
+      
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+    
+    
+
+
+
Index: inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/util.py	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/util.py	(revision 195)
@@ -0,0 +1,408 @@
+
+def distance(v0, v1):
+    from math import sqrt
+
+    return sqrt((v1[0]-v0[0])**2 + (v1[1]-v0[1])**2)
+
+def length(v):
+    from math import sqrt
+
+    return sqrt(v[0]**2 + v[1]**2)
+
+def inner_product(v0, v1):
+    from math import sqrt
+
+    #FIXME: Obsolete
+    try:
+        i = sqrt(v0[0]*v1[0] + v0[1]*v1[1])
+    except:
+        print v0, v1
+        import sys; sys.exit()
+
+    return i    
+
+def angle(v):
+    """Compute angle between e1 (the unit vector in the x-direction)
+    and the specified vector
+    """
+
+    from math import acos, pi
+
+    l = length(v)
+    v1 = v[0]/l
+    v2 = v[1]/l
+        
+    theta = acos(v1)
+
+    if v2 < 0:
+        #Quadrant 3 or 4
+        theta = 2*pi-theta
+
+    return theta
+
+
+def anglediff(v0, v1):
+    """Compute difference between angle of vector x0, y0 and x1, y1.
+    This is used for determining the ordering of vertices,
+    e.g. for checking if they are counter clockwise. 
+    
+    Always return a positive value
+    """
+        
+    from math import pi
+    
+    a0 = angle(v0)
+    a1 = angle(v1)
+
+    #Ensure that difference will be positive
+    if a0 < a1:
+        a0 += 2*pi
+            
+    return a0-a1
+
+
+ 
+def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
+  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
+  
+     Compile FN(s) using compiler CC (e.g. mpicc), 
+     Loader LD and shared flag SFLAG.
+     If CC is absent use default compiler dependent on platform
+     if LD is absent CC is used.
+     if SFLAG is absent platform default is used
+     FNs can be either one filename or a list of filenames
+     In the latter case, the first will be used to name so file.
+  """
+  import os, string, sys, types
+  
+  # Input check
+  #
+  assert not FNs is None, "No filename provided"
+
+  if not type(FNs) == types.ListType:
+    FNs = [FNs]
+
+
+  libext = 'so' #Default extension (Unix)
+  libs = ''
+  version = sys.version[:3]
+  
+  # Determine platform and compiler
+  #
+  if sys.platform == 'sunos5':  #Solaris
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+      
+  elif sys.platform == 'osf1V5':  #Compaq AlphaServer
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'linux2':  #Linux
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'darwin':  #Mac OS X:
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'bundle -flat_namespace -undefined suppress'
+
+  elif sys.platform == 'cygwin':  #Cygwin (compilation same as linux)
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+    libext = 'dll'
+    libs = '/lib/python%s/config/libpython%s.dll.a' %(version,version)
+      
+  elif sys.platform == 'win32':  #Windows
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+    libext = 'dll'
+
+    v = version.replace('.','')
+    dllfilename = 'python%s.dll' %(v)
+    libs = os.path.join(sys.exec_prefix,dllfilename)
+      
+      
+  else:
+    if verbose: print "Unrecognised platform %s - revert to default"\
+                %sys.platform
+    
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = 'ld'
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+
+            
+  # Find location of include files
+  #
+  if sys.platform == 'win32':  #Windows
+    include = os.path.join(sys.exec_prefix, 'include')    
+  else:  
+    include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
+                           'python'+version)
+
+  # Check existence of Python.h
+  #
+  headerfile = include + os.sep +'Python.h'
+  try:
+    open(headerfile, 'r')
+  except:
+    raise """Did not find Python header file %s.
+    Make sure files for Python C-extensions are installed. 
+    In debian linux, for example, you need to install a
+    package called something like python2.3-dev""" %headerfile
+  
+  
+  # Check filename(s)
+  #
+  object_files = ''
+  for FN in FNs:        
+    root, ext = os.path.splitext(FN)
+    if ext == '':
+      FN = FN + '.c'
+    elif ext.lower() != '.c':
+      raise Exception, "Unrecognised extension: " + FN
+    
+    try:
+      open(FN,'r')
+    except:    
+      raise Exception, "Could not open: " + FN
+
+    if not object_files: root1 = root  #Remember first filename        
+    object_files += root + '.o '  
+  
+  
+    # Compile
+    #
+    
+    s = "%s -c %s -I%s -o %s.o -Wall -O" %(compiler, FN, include, root)
+    if verbose:
+      print s
+    else:
+      s = s + ' 2> /dev/null' #Suppress errors
+  
+    try:
+      err = os.system(s)
+      if err != 0:
+          raise 'Attempting to compile %s failed - please try manually' %FN 
+    except:
+      raise 'Could not compile %s - please try manually' %FN  
+
+  
+  # Make shared library (*.so or *.dll)
+  
+  s = "%s -%s %s -o %s.%s %s -lm" %(loader, sharedflag, object_files, root1, libext, libs)
+  if verbose:
+    print s
+  else:
+    s = s + ' 2> /dev/null' #Suppress warnings
+  
+  try:  
+    err=os.system(s)
+    if err != 0:	
+	raise 'Atempting to link %s failed - please try manually' %root1     
+  except:
+    raise 'Could not link %s - please try manually' %root1
+    
+
+
+
+def can_use_C_extension(filename):
+    """Determine whether specified C-extension
+    can and should be used.
+    """
+
+    from config import use_extensions
+
+    from os.path import splitext
+
+    root, ext = splitext(filename)
+    
+    C=False
+    if use_extensions:
+        try:
+            s = 'import %s' %root
+            #print s
+            exec(s)
+        except:
+            try:
+                open(filename)
+            except:
+                msg = 'C extension %s cannot be opened' %filename
+                print msg                
+            else:    
+                print '------- Trying to compile c-extension %s' %filename
+	    
+                compile(filename)		    
+                try:
+                    compile(filename)
+                except:
+                    print 'WARNING: Could not compile C-extension %s'\
+                          %filename
+                else:
+                    try:
+                        exec('import %s' %root)
+                    except:
+                        msg = 'C extension %s seems to compile OK, '
+                        msg += 'but it can still not be imported.'
+                        raise msg
+                    else:
+                        C=True
+        else:
+            C=True
+            
+    if not C:
+        pass
+        print 'NOTICE: C-extension %s not used' %filename
+
+    return C
+
+
+
+
+
+####################################################################
+#Python versions of function that are also implemented in util_ext.c
+#
+
+def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
+    """
+    """
+    
+    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
+    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
+    a /= det
+
+    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
+    b /= det            
+
+    return a, b
+
+
+
+def rotate_python(q, normal, direction = 1):
+    """Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+
+    If direction is negative the rotation is inverted.
+    
+    Input vector is preserved
+
+    This function is specific to the shallow water wave equation
+    """
+
+    #FIXME: Needs to be tested
+
+    from Numeric import zeros, Float
+    
+    assert len(q) == 3,\
+           'Vector of conserved quantities must have length 3'\
+           'for 2D shallow water equation'
+    try:
+        l = len(normal)
+    except:
+        raise 'Normal vector must be an Numeric array'
+
+    assert l == 2, 'Normal vector must have 2 components'
+
+ 
+    n1 = normal[0]
+    n2 = normal[1]    
+    
+    r = zeros(len(q), Float) #Rotated quantities
+    r[0] = q[0]              #First quantity, height, is not rotated
+
+    if direction == -1:
+        n2 = -n2
+
+
+    r[1] =  n1*q[1] + n2*q[2]
+    r[2] = -n2*q[1] + n1*q[2]
+    
+    return r
+
+
+
+
+##############################################
+#Initialise module
+
+import util
+if util.can_use_C_extension('util_ext.c'):
+    from util_ext import gradient, rotate
+else:
+    gradient = gradient_python
+    rotate = rotate_python
+
+
+
+
+if __name__ == "__main__":
+    pass
+
Index: inundation/ga/storm_surge/pyvolution/util_ext.c
===================================================================
--- inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 195)
+++ inundation/ga/storm_surge/pyvolution/util_ext.c	(revision 195)
@@ -0,0 +1,225 @@
+// Python - C extension for finite_volumes util module.
+//
+// To compile (Python2.3):
+//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
+//  gcc -shared util_ext.o  -o util_ext.so	
+//
+// See the module util.py
+//
+//
+// Ole Nielsen, GA 2004
+	
+#include "Python.h"
+#include "Numeric/arrayobject.h"
+#include "math.h"
+
+
+int _gradient(double x0, double y0, 
+	      double x1, double y1, 
+	      double x2, double y2, 
+	      double *q0, double *q1, double *q2, 
+	      double *a, double *b, 
+	      int N) {
+
+  double det;
+  int i;
+  
+  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
+  
+  for (i=0; i<N; i++) {
+    a[i] = (y2-y0)*(q1[i]-q0[i]) - (y1-y0)*(q2[i]-q0[i]);
+    a[i] /= det;
+
+    b[i] = (x1-x0)*(q2[i]-q0[i]) - (x2-x0)*(q1[i]-q0[i]);
+    b[i] /= det;
+  }  
+
+  return 0;
+}
+
+
+
+// Computational function for rotation
+int _rotate(double *q, double *r, double *normal, int direction) {
+  /*Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+    
+    To rotate in opposite direction, call rotate with direction=1
+    Result is returned in r 
+    
+    q and r are assumed to have three components each*/
+
+
+  double n0, n1, q1, q2;
+  
+  //Determine normal and direction
+  n0 = normal[0];
+  if (direction == 1) {
+    n1 = normal[1];    
+  } else {
+    n1 = -normal[1];
+  }
+
+  //Shorthands
+  q1 = q[1];  //uh momentum
+  q2 = q[2];  //vh momentum
+
+  //Rotate
+  r[0] = q[0];
+  r[1] = n0*q1 + n1*q2;
+  r[2] = n0*q2 - n1*q1;
+
+  return 0;
+}  
+
+
+
+// Gateway to Python
+PyObject *gradient(PyObject *self, PyObject *args) {
+  //
+  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+  //
+  
+  double x0, y0, x1, y1, x2, y2;
+  PyObject *Q0, *Q1, *Q2, *result;
+  PyArrayObject *q0, *q1, *q2, *a, *b;
+  int dimensions[1];
+  int N;
+
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "ddddddOOO", &x0, &y0, &x1, &y1, &x2, &y2,
+			&Q0, &Q1, &Q2))  
+    return NULL;
+
+  //Input check
+  if PyArray_Check(Q0) {
+    q0 = (PyArrayObject *) 
+      PyArray_ContiguousFromObject(Q0, PyArray_DOUBLE, 0, 0);
+  } else {
+    PyErr_SetString(PyExc_TypeError, "q0 must be a numeric array");
+    return NULL;
+  }
+    
+  if PyArray_Check(Q1) {
+    q1 = (PyArrayObject *) 
+      PyArray_ContiguousFromObject(Q1, PyArray_DOUBLE, 0, 0);
+  } else {
+    PyErr_SetString(PyExc_TypeError, "q1 must be a numeric array");
+    return NULL;
+  }
+  
+  if PyArray_Check(Q2) {
+    q2 = (PyArrayObject *) 
+      PyArray_ContiguousFromObject(Q2, PyArray_DOUBLE, 0, 0);
+  } else {
+    PyErr_SetString(PyExc_TypeError, "q2 must be a numeric array");
+    return NULL;
+  }  
+    
+  //Allocate space for return vectors a and b  
+  
+  N = q0 -> dimensions[0];
+  dimensions[0] = N;
+  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  
+
+  
+  // Call underlying routine
+  _gradient(x0, y0, x1, y1, x2, y2, 
+	    (double *) q0 -> data, 
+	    (double *) q1 -> data, 
+	    (double *) q2 -> data,
+	    (double *) a -> data, 
+	    (double *) b-> data, N);
+
+  Py_DECREF(q0);
+  Py_DECREF(q1);
+  Py_DECREF(q2);    
+	      
+  // Return arrays a and b
+  result = Py_BuildValue("OO", a, b);
+  Py_DECREF(a);
+  Py_DECREF(b);  
+  return result;
+}
+
+
+// Gateway to Python
+PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
+  //
+  // r = rotate(q, normal, direction=0)
+  //
+  // Where q is assumed to be a Float numeric array of length 3 and
+  // normal a Float numeric array of length 2.
+
+  //FIXME: Input checks and unit tests
+  
+  PyObject *Q, *Normal;
+  PyArrayObject *q, *r, *normal;
+  
+  static char *argnames[] = {"q", "normal", "direction", NULL};
+  int dimensions[1];
+  int N, direction=0;
+
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
+				   &Q, &Normal, &direction))  
+    return NULL;  
+
+  //Input checks (convert sequenced into numeric arrays)
+  
+  q = (PyArrayObject *) 
+    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
+  normal = (PyArrayObject *) 
+    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
+  
+  //Allocate space for return vector r (don't DECREF) 
+  N = q -> dimensions[0];
+  dimensions[0] = N;
+  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+
+  //Rotate
+  _rotate((double *) q -> data,
+	  (double *) r -> data,
+	  (double *) normal -> data,	  
+	  direction);	  
+
+  //Build result and DECREF r - returning r directly causes memory leak
+  //result = Py_BuildValue("O", r);	  
+  //Py_DECREF(r);
+  
+  //return result;
+  return PyArray_Return(r);
+}    
+
+
+// Method table for python module
+//static struct PyMethodDef MethodTable[] = {
+//  {"gradient", gradient, METH_VARARGS},  
+//  {NULL, NULL}
+//};
+
+
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+  
+  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},  
+  {"gradient", gradient, METH_VARARGS, "Print out"},    
+  {NULL, NULL, 0, NULL}   /* sentinel */
+};
+
+
+
+// Module initialisation   
+void initutil_ext(void){
+  Py_InitModule("util_ext", MethodTable);
+  
+  import_array();     //Necessary for handling of NumPY structures  
+}
+
