Index: /trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py	(revision 8522)
+++ /trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py	(revision 8523)
@@ -149,5 +149,4 @@
 barrier()
 t0 = time.time()
-
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     if myid == 0:
Index: /trunk/anuga_core/validation_tests/Tests/Experimental/Isolated_Building/results.tex
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Experimental/Isolated_Building/results.tex	(revision 8522)
+++ /trunk/anuga_core/validation_tests/Tests/Experimental/Isolated_Building/results.tex	(revision 8523)
@@ -1,12 +1,27 @@
 
-\section{Isolated Building}
+\section{Isolated Obstacle}
 
-Will show a comparison of \anuga{} to the isolated building test.
+Herein \anuga{} is tested against the isolated obstacle test, which involves a dam-break flow past an isolated building. This benchmark experiment is reported in:
+
+Soares-Frazao, S. and Zech, Y. (2007) Experimental study of dam-break flow against an isolated obstacle. Journal of Hydraulic Resaerch 45(Extra Issue):27-36
+
+The data is provided in a DVD-ROM associated with the above JHR issue.
+
+Note that many models have been applied to this problem, and their performance is reviewed in Soares Frazao, S., Noel, B., Spinewine, B., and Zech, Y. (YEAR) The isolated building test case: results from the IMPACT benchmark. At the time of writing, the performance of \anuga{} on this problem is qualitatively similar to the performance of most other models -- i.e. it simulates the main aspects of the flow, but often not the details. It seems likely that these experimental flows may well partly violate the shallow water approximations, and to the extent that this is true, we would not expect perfect agreement.
 
 \subsection{Results}
+For details of the experiment which are required to run the model, consult Soares-Frazao and Zech (2007). A zoom of the model geometry around the building and reservoir mouth is shown in Figure~\ref{urbanflow_geom}, along with measurement stations. 
+\begin{figure}
+\includegraphics{XXX}
+\caption{Zoom of the flume planform geometry near the isolated building. Red areas are obstacles which are never submerged by the flow. See Soares-Frazo and Zech (2007) for further details on the flume geometry.}
+\label{urbanflow_geom}
+\end{figure}
 
-
-Will have a bunch of figures.
-
+Water was released from the reservoir, and allowed to flow past the building. Figures XX - XX2 compare the stage and velocity at each of these measurement points with the modelled result. 
+\begin{figure}
+\includegraphics{XXX}
+\caption{Comparison betweeen modelled and observed stage, x-velocity and y-velocity at sites G1-G6.}
+\label{urbanflow_geom}
+\end{figure}
 
 
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/__init__.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/__init__.py	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/__init__.py	(revision 8523)
@@ -0,0 +1,340 @@
+"""Make directory available as a Python package
+"""
+
+# Add path of package to PYTHONPATH to allow C-extensions to be loaded
+import sys
+sys.path += __path__
+
+# Make selected classes available directly
+from anuga_tsunami.swb2_domain import Domain as Domain
+from anuga_tsunami.swb2_domain import Domain
+
+#-----------------------------------------------------
+# Make selected classes available directly
+#-----------------------------------------------------
+from anuga.__metadata__ import __version__, __date__, __author__
+
+#from anuga.shallow_water.shallow_water_domain import Domain
+
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+from anuga.abstract_2d_finite_volumes.util import file_function, \
+                                        sww2timeseries, sww2csv_gauges
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross, \
+                                                    rectangular
+
+from anuga.file.csv_file import load_csv_as_building_polygons,  \
+                                load_csv_as_polygons
+
+from anuga.file.sts import create_sts_boundary
+
+from anuga.file.ungenerate import load_ungenerate
+
+from anuga.geometry.polygon import read_polygon, plot_polygons, polygon_area, \
+                                   inside_polygon, polygon_area
+from anuga.geometry.polygon_function import Polygon_function
+
+from anuga.abstract_2d_finite_volumes.pmesh2domain import \
+                                            pmesh_to_domain_instance
+
+from anuga.utilities.system_tools import file_length
+from anuga.utilities.sww_merge import sww_merge
+from anuga.utilities.file_utils import copy_code_files
+from anuga.utilities.numerical_tools import safe_acos as acos
+
+from anuga.geometry.polygon import read_polygon
+from anuga.caching import cache
+
+#-----------------------------
+# Standard Boundaries
+#-----------------------------
+from anuga.shallow_water.boundaries import File_boundary
+from anuga.shallow_water.boundaries import Reflective_boundary
+from anuga.shallow_water.boundaries import Field_boundary
+from anuga.shallow_water.boundaries import \
+                    Transmissive_stage_zero_momentum_boundary
+from anuga.shallow_water.boundaries import \
+                    Transmissive_momentum_set_stage_boundary
+from anuga.shallow_water.boundaries import \
+                    Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
+
+
+#-----------------------------
+# SWW-specific Boundaries
+#-----------------------------
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
+                            import Dirichlet_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
+                            import Time_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
+                            import Time_space_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
+                            import Transmissive_boundary
+
+
+
+#-----------------------------
+# Shalow Water Tsunamis
+#-----------------------------
+
+from anuga.shallow_water.smf import slide_tsunami, slump_tsunami
+
+
+#-----------------------------
+# Forcing
+#-----------------------------
+from anuga.shallow_water.forcing import Inflow, Rainfall, Wind_stress
+
+#-----------------------------
+# File conversion utilities
+#-----------------------------
+from anuga.file_conversion.file_conversion import sww2obj, \
+                    timefile2netcdf, tsh2sww
+from anuga.file_conversion.urs2nc import urs2nc
+from anuga.file_conversion.urs2sww import urs2sww  
+from anuga.file_conversion.urs2sts import urs2sts
+from anuga.file_conversion.dem2pts import dem2pts                    
+from anuga.file_conversion.esri2sww import esri2sww   
+from anuga.file_conversion.sww2dem import sww2dem, sww2dem_batch 
+from anuga.file_conversion.asc2dem import asc2dem     
+from anuga.file_conversion.ferret2sww import ferret2sww     
+from anuga.file_conversion.dem2dem import dem2dem
+
+
+#-----------------------------
+# Mesh API
+#-----------------------------
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+#-----------------------------
+# SWW file access
+#-----------------------------
+from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
+    
+
+
+#---------------------------
+# Operators
+#---------------------------
+from anuga.operators.base_operator import Operator
+
+#---------------------------
+# Structure Operators
+#---------------------------
+from anuga.structures.structure_operator import Structure_operator
+from anuga.structures.boyd_box_operator import Boyd_box_operator
+from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+
+#---------------------------
+# User Access Functions
+#---------------------------
+
+from anuga.utilities.system_tools import get_user_name, get_host_name, get_revision_number
+
+
+
+
+#-----------------------------
+# rectangular domains
+#-----------------------------
+def rectangular_cross_domain(*args, **kwargs):
+    """
+    Create a rectangular domain with triangulation made
+    up of m+1 by n+1 uniform rectangular cells divided
+    into 4 triangles in a cross pattern
+
+    Arguments
+    m:      number of cells in x direction
+    n:      number of cells in y direction
+    len1:   length of domain in x direction (left to right)
+            (default 1.0)
+    len2:   length of domain in y direction (bottom to top)
+            (default 1.0)
+    origin: tuple (x,y) specifying location of lower left corner
+            of domain (default (0,0))
+    """
+    points, vertices, boundary = rectangular_cross(*args, **kwargs)
+    return Domain(points, vertices, boundary)
+
+#----------------------------
+# Create domain from file
+#----------------------------
+def create_domain_from_file(file, DomainClass=Domain):
+    """
+    Create a domain from a file
+    """
+    return pmesh_to_domain_instance(file,DomainClass=DomainClass)
+
+#---------------------------
+# Create domain from regions
+#---------------------------
+
+def create_domain_from_regions(bounding_polygon,
+                               boundary_tags,
+                               maximum_triangle_area=None,
+                               mesh_filename=None,
+                               interior_regions=None,
+                               interior_holes=None,
+                               hole_tags=None,
+                               poly_geo_reference=None,
+                               mesh_geo_reference=None,
+                               minimum_triangle_angle=28.0,
+                               fail_if_polygons_outside=True,
+                               use_cache=False,
+                               verbose=True):
+    """Create domain from bounding polygons and resolutions.
+
+    bounding_polygon is a list of points in Eastings and Northings,
+    relative to the zone stated in poly_geo_reference if specified.
+    Otherwise points are just x, y coordinates with no particular 
+    association to any location.
+
+    boundary_tags is a dictionary of symbolic tags. For every tag there
+    is a list of indices referring to segments associated with that tag.
+    If a segment is omitted it will be assigned the default tag ''.
+
+    maximum_triangle_area is the maximal area per triangle
+    for the bounding polygon, excluding the  interior regions.
+
+    Interior_regions is a list of tuples consisting of (polygon,
+    resolution) for each region to be separately refined. Do not have
+    polygon lines cross or be on-top of each other.  Also do not have
+    polygon close to each other.
+    
+    NOTE: If a interior_region is outside the bounding_polygon it should 
+    throw an error
+    
+    interior_holes is a list of ploygons for each hole. These polygons do not
+    need to be closed, but their points must be specified in a counter-clockwise
+    order.
+
+    hole_tags  is a list of tag segment dictionaries.
+
+    This function does not allow segments to share points - use underlying
+    pmesh functionality for that
+
+    poly_geo_reference is the geo_reference of the bounding polygon and
+    the interior polygons.
+    If none, assume absolute.  Please pass one though, since absolute
+    references have a zone.
+    
+    mesh_geo_reference is the geo_reference of the mesh to be created.
+    If none is given one will be automatically generated.  It was use
+    the lower left hand corner of  bounding_polygon (absolute)
+    as the x and y values for the geo_ref.
+    
+    Returns the shallow water domain instance
+
+    Note, interior regions should be fully nested, as overlaps may cause
+    unintended resolutions. 
+
+    fail_if_polygons_outside: If True (the default) Exception in thrown
+    where interior polygons fall outside bounding polygon. If False, these
+    will be ignored and execution continued.
+        
+    
+    """
+
+
+    # Build arguments and keyword arguments for use with caching or apply.
+    args = (bounding_polygon,
+            boundary_tags)
+    
+    kwargs = {'maximum_triangle_area': maximum_triangle_area,
+              'mesh_filename': mesh_filename,
+              'interior_regions': interior_regions,
+              'interior_holes': interior_holes,
+              'hole_tags': hole_tags,
+              'poly_geo_reference': poly_geo_reference,
+              'mesh_geo_reference': mesh_geo_reference,
+              'minimum_triangle_angle': minimum_triangle_angle,
+              'fail_if_polygons_outside': fail_if_polygons_outside,
+              'verbose': verbose} #FIXME (Ole): See ticket:14
+
+    # Call underlying engine with or without caching
+    if use_cache is True:
+        try:
+            from anuga.caching import cache
+        except:
+            msg = 'Caching was requested, but caching module'+\
+                  'could not be imported'
+            raise (msg)
+
+
+        domain = cache(_create_domain_from_regions,
+                       args, kwargs,
+                       verbose=verbose,
+                       compression=False)
+    else:
+        domain = apply(_create_domain_from_regions,
+                       args, kwargs)
+
+    return domain
+
+        
+def _create_domain_from_regions(bounding_polygon,
+                                boundary_tags,
+                                maximum_triangle_area=None,
+                                mesh_filename=None,                           
+                                interior_regions=None,
+                                interior_holes=None,
+                                hole_tags=None,
+                                poly_geo_reference=None,
+                                mesh_geo_reference=None,
+                                minimum_triangle_angle=28.0,
+                                fail_if_polygons_outside=True,
+                                verbose=True):
+    """_create_domain_from_regions - internal function.
+
+    See create_domain_from_regions for documentation.
+    """
+
+    #from anuga.shallow_water.shallow_water_domain import Domain
+    from anuga_tsunami.swb2_domain import Domain
+    from anuga.pmesh.mesh_interface import create_mesh_from_regions
+    
+    create_mesh_from_regions(bounding_polygon,
+                             boundary_tags,
+                             maximum_triangle_area=maximum_triangle_area,
+                             interior_regions=interior_regions,
+                             filename=mesh_filename,
+                             interior_holes=interior_holes,
+                             hole_tags=hole_tags,
+                             poly_geo_reference=poly_geo_reference,
+                             mesh_geo_reference=mesh_geo_reference,
+                             minimum_triangle_angle=minimum_triangle_angle,
+                             fail_if_polygons_outside=fail_if_polygons_outside,
+                             use_cache=False,
+                             verbose=verbose)
+
+    domain = Domain(mesh_filename, use_cache=False, verbose=verbose)
+
+
+    return domain
+   
+import logging as log
+
+from anuga.config import use_psyco, g, velocity_protection
+if use_psyco:
+    # try using psyco if available
+    try:
+        import psyco
+    except:
+        import os
+        import sys
+        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
+            pass
+            # Psyco isn't supported on 64 bit systems, but it doesn't matter
+        elif sys.version[:3] == '2.7' :
+            pass
+            # Psyco isn't available for python 2.7 (16/05/2011)
+        else:
+            log.critical('WARNING: psyco (speedup) could not be imported, '
+                         'you may want to consider installing it')
+    else:
+        psyco.full() # aggressively compile everything
+        #psyco.background() # attempt to profile code - only compile most used
+        
+
+ 
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/compile_all.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/compile_all.py	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/compile_all.py	(revision 8523)
@@ -0,0 +1,23 @@
+import os
+import time
+
+buildroot = '..' + os.sep + '..' + os.sep + '..' + os.sep + '..' + os.sep + '..' + os.sep + 'anuga_core'
+
+compile_script = buildroot + os.sep + 'source' + os.sep + 'anuga' + \
+                 os.sep + 'utilities' + os.sep + 'compile.py'
+
+t0 = time.time()
+
+# Attempt to compile all ANUGA extensions
+
+execfile(compile_script)
+
+
+  
+
+
+print        
+print 'That took %.3fs' %(time.time() - t0)
+
+if sys.platform == 'win32':
+    raw_input('Press the RETURN key')
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/runup.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/runup.py	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/runup.py	(revision 8523)
@@ -0,0 +1,112 @@
+"""Runup example from the manual, slightly modified
+"""
+#---------
+#Import Modules
+#--------
+from sys import path
+import anuga
+
+import numpy
+
+import struct
+#import scipy
+
+#from Numeric import *
+
+from math import sin, pi, exp
+#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
+#from anuga.shallow_water.shallow_water_domain import Domain as Domain
+#from shallow_water_balanced_steve.swb_domain import *
+#import shallow_water_balanced_steve.swb_domain 
+#from shallow_water_balanced_steve.swb_domain import Domain as Domain
+#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/shallow_water_balanced_steve')
+#from swb_domain import *
+#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
+from balanced_basic.swb2_domain import *
+#---------
+#Setup computational domain
+#---------
+points, vertices, boundary = anuga.rectangular_cross(40,40)
+
+domain=Domain(points,vertices,boundary)    # Create Domain
+domain.set_name('runup_v2')                         # Output to file runup.sww
+domain.set_datadir('.')                          # Use current folder
+domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+domain.set_store_vertices_uniquely(True)
+#------------------
+# Define topography
+#------------------
+
+def topography(x,y):
+	return -x/2 #+0.05*numpy.sin((x+y)*50.0) #+0.1*(numpy.random.rand(len(x)) -0.5)       # Linear bed slope + small random perturbation 
+
+def stagefun(x,y):
+    #stg=-0.2*(x<0.5) -0.1*(x>=0.5)
+    stg=-0.2 # Stage
+    #topo=topography(x,y) #Bed
+    return stg #*(stg>topo) + topo*(stg<=topo)
+
+domain.set_quantity('elevation',topography)     # Use function for elevation
+domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
+
+domain.set_quantity('friction',0.00)             # Constant friction
+
+domain.set_quantity('stage', stagefun)              # Constant negative initial stage
+
+#--------------------------
+# Setup boundary conditions
+#--------------------------
+Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
+Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
+Bw=anuga.Time_boundary(domain=domain,
+	f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
+
+#----------------------------------------------
+# Associate boundary tags with boundary objects
+#----------------------------------------------
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+
+#------------------------------
+#Evolve the system through time
+#------------------------------
+
+#xwrite=open("xvel.out","wb")
+#ywrite=open("yvel.out","wb")
+## Set print options to be compatible with file writing via the 'print' statement 
+#numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
+
+for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
+    print domain.timestepping_statistics()
+
+#    momx=domain.quantities['xmomentum'].centroid_values
+#    momy=domain.quantities['ymomentum'].centroid_values
+#    #mom2=domain.quantities['xmomentum'].vertex_values
+#    dep1=(domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values+1.0e-06)
+#    #dep2=(domain.quantities['stage'].vertex_values-domain.quantities['elevation'].vertex_values+1.0e-06)
+#    velx=momx/dep1
+#    vely=momy/dep1
+#    #vel2=mom2/dep2
+#    #print vel1.max(), vel2.max()
+#    #print vel1.min(), vel2.min()
+#    #xwrite.write(velx)
+#    #print >> xwrite, str(velx)
+#    #print >> ywrite, str(vely)
+#    for i in velx:
+#        data=struct.pack('f',i)
+#        xwrite.write(data)
+#
+#    for j in vely:
+#        data2=struct.pack('f',j)
+#        ywrite.write(data2)
+#
+#    #print >> xwrite, \n
+#    #numpy.savetxt("xvel.txt",velx)
+#    #numpy.savetxt("yvel.txt",vely)
+#    #velx.tofile(xwrite," ")
+#    #vely.tofile(ywrite," ")
+#
+#xwrite.close()
+#ywrite.close()
+
+print 'Finished'
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_boundary_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_boundary_conditions.py	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_boundary_conditions.py	(revision 8523)
@@ -0,0 +1,301 @@
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Boundary, File_boundary
+
+#from swb2_domain_ext import rotate
+#####
+
+#class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        CAN BE UNSTABLE -- CAREFUL
+#        Sets boundary values of h, uh, vh as in transmissive boundary, but 
+#        also ensures:
+#        flux[0] = 0.0,
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """Transmissive boundaries return the edge values
+#        of the volume they serve.
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#        return q
+#
+######
+#
+#class  zero_mass_flux_zero_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=0. 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=vh 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=1.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=1.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+class Characteristic_boundary(Boundary):
+    """
+
+    Example:
+
+    def waveform(t):
+        return sea_level + normalized_amplitude/cosh(t-25)**2
+
+    Bf = Nudge_boundary(domain, waveform)
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain=None, function=None):
+        """ Instantiate a
+            Nudge_boundary.
+            domain is the domain containing the boundary
+            function is the function to apply
+        """
+
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise Exception, msg
+
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise Exception, msg
+
+        self.domain = domain
+        self.function = function
+
+
+    def __repr__(self):
+        """ Return a representation of this instance. """
+        msg = 'Nudge_boundary'
+        msg += '(%s)' % self.domain
+        return msg
+
+
+    def evaluate(self, vol_id, edge_id):
+        """
+        """
+
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
+        depth=max(q[0]-bed,0.0)
+        dt=self.domain.timestep
+
+        normal = self.domain.get_normal(vol_id, edge_id)
+
+
+        t = self.domain.get_time()
+
+        if hasattr(self.function, 'time'):
+            # Roll boundary over if time exceeds
+            while t > self.function.time[-1]:
+                msg = 'WARNING: domain time %.2f has exceeded' % t
+                msg += 'time provided in '
+                msg += 'Flather_boundary object.\n'
+                msg += 'I will continue, reusing the object from t==0'
+                log.critical(msg)
+                t -= self.function.time[-1]
+
+        value = self.function(t)
+        try:
+            x = float(value)
+        except:
+            x = float(value[0])
+
+        ##
+        ndotq = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
+
+        if(depth==0.):
+            q[0] = x
+            q[1] = 0. 
+            q[2] = 0. 
+
+        else:
+
+            # Asssume sub-critical flow. Set the values of the characteristics as
+            # appropriate, depending on whether we have inflow or outflow
+            tmp = (9.8/depth)**0.5
+            if(ndotq>0.):
+                # Only need to impose one characteristic
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow external velocities to be specified in the function
+                # This would be needed to use the boundary for real inflows (e.g. run_wave.py)
+                # At the moment, it seems to work okay as a 'passive' outflow boundary (e.g. run_wave.py)
+                # Here it attenuates the wave at the outflow a bit, but doesn't do too much damage compared
+                # with the other boundary conditions
+                w1 = 0. - tmp*x
+                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+                
+            else:
+                # Need to set 2 characteristics
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow this to be specified with a function
+                w1 = 0. - tmp*x
+                w2 = 0.
+                #w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+
+
+            q[0] = (w3-w1)/(2*tmp)
+            qperp= (w3+w1)/2. *depth
+            qpar=  w2*depth
+
+            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
+
+            q[1] = qperp*normal[0] + qpar*normal[1]
+            q[2] = qperp*normal[1] -qpar*normal[0]
+
+        return q
+###
+def rotate(q,normal,direction=1):
+    #Function for rotation -- better to mimic the shallow water one
+    # FIXME: Use C version
+    out=q*1.      
+
+    if(direction==1):
+        out[1] = q[1]*normal[0] + q[2]*normal[1]    
+        out[2] = -q[1]*normal[1] + q[2]*normal[0]
+    elif(direction==-1):
+        out[1] = q[1]*normal[0] - q[2]*normal[1]    
+        out[2] = +q[1]*normal[1] + q[2]*normal[0]
+    else:
+        print 'ERROR IN rotate'
+        assert 0==1
+
+    return out
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain.py	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain.py	(revision 8523)
@@ -0,0 +1,562 @@
+
+from anuga.shallow_water.shallow_water_domain import *
+from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
+
+
+##############################################################################
+# Shallow Water Balanced Domain -- alternative implementation
+# 
+# FIXME: Following the methods in CITE MODSIM PAPER
+#
+##############################################################################
+
+class Domain(Sww_domain):
+
+    def __init__(self,
+                 coordinates=None,
+                 vertices=None,
+                 boundary=None,
+                 tagged_elements=None,
+                 geo_reference=None,
+                 use_inscribed_circle=False,
+                 mesh_filename=None,
+                 use_cache=False,
+                 verbose=False,
+                 full_send_dict=None,
+                 ghost_recv_dict=None,
+                 starttime=0.0,
+                 processor=0,
+                 numproc=1,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
+
+        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
+
+        evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum']         
+        other_quantities   = [ 'elevation', 'friction', 'height', 
+                               'xvelocity', 'yvelocity', 'x', 'y' ]
+
+        
+        Sww_domain.__init__(self,
+                            coordinates = coordinates,
+                            vertices = vertices,
+                            boundary = boundary,
+                            tagged_elements = tagged_elements,
+                            geo_reference = geo_reference,
+                            use_inscribed_circle = use_inscribed_circle,
+                            mesh_filename = mesh_filename,
+                            use_cache = use_cache,
+                            verbose = verbose,
+                            conserved_quantities = conserved_quantities,
+                            evolved_quantities = evolved_quantities,
+                            other_quantities = other_quantities,
+                            full_send_dict = full_send_dict,
+                            ghost_recv_dict = ghost_recv_dict,
+                            starttime = starttime,
+                            processor = processor,
+                            numproc = numproc,
+                            number_of_full_nodes = number_of_full_nodes,
+                            number_of_full_triangles = number_of_full_triangles)
+       
+        #------------------------------------------------
+        # set some defaults
+        # Most of these override the options in config.py
+        #------------------------------------------------
+        self.set_CFL(1.0)
+        self.set_use_kinematic_viscosity(False)
+        self.timestepping_method='rk2' 
+        # The following allows storage of the negative depths associated with this method
+        self.minimum_storable_height=-99999999999.0 
+
+        self.use_edge_limiter=True
+        self.default_order=2
+        self.extrapolate_velocity_second_order=True
+
+        # Note that the extrapolation method used in quantity_ext.c (e.g.
+        # extrapolate_second_order_and_limit_by_edge) uses a constant value for
+        # all the betas.  
+        self.beta_w=1.0
+        self.beta_w_dry=1.0
+        self.beta_uh=1.0
+        self.beta_uh_dry=1.0
+        self.beta_vh=1.0
+        self.beta_vh_dry=1.0
+
+        #self.optimise_dry_cells=True
+        # "self.optimise_dry_cells=False" presently ensures that the stage is
+        # always >= minimum_bed_edge value.  Actually, we should only need to
+        # apply 'False' on the very first time step (to deal with stage values
+        # that were initialised below the bed by the user). After this, the
+        # algorithm should take care of itself, and 'True' should be okay.
+        self.optimise_dry_cells=False 
+
+        # Because gravity is treated within the flux function, 
+        # we remove it from the forcing terms.
+        #self.forcing_terms.remove(gravity)
+
+        # We need the edge_coordinates for the extrapolation
+        self.edge_coordinates=self.get_edge_midpoint_coordinates()
+
+        # We demand that vertex values are stored uniquely
+        self.set_store_vertices_smoothly(False)
+
+        self.maximum_allowed_speed=0.0
+        #self.forcing_terms.append(manning_friction_explicit)
+        #self.forcing_terms.remove(manning_friction_implicit)
+
+        print '##########################################################################'
+        print '#'
+        print '# Using anuga_tsunami solver in anuga_work/development/gareth/experimental/anuga_tsunami/'
+        print "#"
+        print "# This solver is experimental. Here are some tips on using it"
+        print "#"
+        print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
+        print "# , as the latter can be completely misleading near strong gradients in the flow. "
+        print "# There is a plot_util.py script in anuga_core/utilities/ which might help you extract"
+        print "# quantities at centroid values from sww files."
+        print "# Note that to accuractely compute centroid values from sww files, the files need to store "
+        print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
+        print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
+        print "# ANUGA viewer as well -- I expect this would be lots of work)"
+        print "#"
+        print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
+        print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
+        print "#"
+        print "# 3) This solver is not expected to perform well in problems with very"
+        print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
+        print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests"
+        print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) "
+        print "#"
+        print "# 4) This solver allows the stage_centroid_value to drop to slightly below the minimum bed_vertex_value"
+        print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
+        print "# bed_centroid_value. This means that triangles store slightly more water than they are"
+        print "# typically interpreted to store, which might have significance in some applications."
+        print "#"
+        print "# You will probably be able to tell this is causing you problems by convergence testing"
+        print "#"
+        print '# 5) Note that many options in config.py have been overridden by the solver -- I have '
+        print '# deliberately attempted to get the solver to perform well with consistent values of '
+        print '# these parameters -- so I would advise against changing them unless you at least check that '
+        print '# it really does improve things'
+        print '#'
+        print '##########################################################################'
+
+    #-------------------------------
+    # Specialisation of Evolve
+    #-------------------------------
+    # This alters the 'evolve' method from shallow_water_domain.py so that just
+    # before the file-output step, the vertex values are replaced with the
+    # centroid values. 
+    # This is useful so that we can unambigously know what the centroid values
+    # of various parameters are
+
+    def evolve(self,
+               yieldstep=None,
+               finaltime=None,
+               duration=None,
+               skip_initial_step=False):
+        """Specialisation of basic evolve method from parent class.
+        
+            Evolve the model by 1 step.
+        """
+
+        # Call check integrity here rather than from user scripts
+        # self.check_integrity()
+
+        #print 'Into Evolve'
+
+        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
+        assert 0 <= self.beta_w <= 2.0, msg
+
+        # Initial update of vertex and edge values before any STORAGE
+        # and or visualisation.
+        # This is done again in the initialisation of the Generic_Domain
+        # evolve loop but we do it here to ensure the values are ok for storage.
+        self.distribute_to_vertices_and_edges()
+
+        if self.store is True and self.get_time() == self.get_starttime():
+            self.initialise_storage()
+        #print 'Into Generic_Domain Evolve'
+        # Call basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
+                                       finaltime=finaltime, duration=duration,
+                                       skip_initial_step=skip_initial_step):
+            #print 'Out of Generic_Domain Evolve'
+            # Store model data, e.g. for subsequent visualisation
+            if self.store is True:
+                # FIXME: Extrapolation is done before writing, because I had
+                # trouble correctly computing the centroid values from the
+                # vertex values (an 'average' did not seem to work correctly in
+                # very shallow cells -- there was a discrepency between
+                # domain.quantity['blah'].centroid_values and the values
+                # computed from the sww using the vertex averge). There should
+                # be a much much more disk-efficient way to store the centroid
+                # values than this
+                self.extrapolate_second_order_edge_sw()
+
+                # Store the timestep
+                self.store_timestep()
+
+            # Pass control on to outer loop for more specific actions
+            yield(t)
+    
+    #-----------------
+    # Flux computation
+    #-----------------
+
+    ## @brief Compute fluxes and timestep suitable for all volumes in domain.
+    # @param domain The domain to calculate fluxes for.
+    def compute_fluxes(domain):
+        """Compute fluxes and timestep suitable for all volumes in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        explicit_update for each of the three conserved quantities
+        stage, xmomentum and ymomentum
+
+        The maximal allowable speed computed by the flux_function for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+          domain.explicit_update is reset to computed flux values
+          domain.timestep is set to the largest step satisfying all volumes.
+
+        This wrapper calls the underlying C version of compute fluxes
+        """
+
+        import sys
+        from swb2_domain_ext import compute_fluxes_ext_central \
+                                      as compute_fluxes_ext
+
+        #print "."
+
+        # Shortcuts
+        Stage = domain.quantities['stage']
+        Xmom = domain.quantities['xmomentum']
+        Ymom = domain.quantities['ymomentum']
+        Bed = domain.quantities['elevation']
+
+        timestep = float(sys.maxint)
+
+        flux_timestep = compute_fluxes_ext(timestep,
+                                           domain.epsilon,
+                                           domain.H0,
+                                           domain.g,
+                                           domain.neighbours,
+                                           domain.neighbour_edges,
+                                           domain.normals,
+                                           domain.edgelengths,
+                                           domain.radii,
+                                           domain.areas,
+                                           domain.tri_full_flag,
+                                           Stage.edge_values,
+                                           Xmom.edge_values,
+                                           Ymom.edge_values,
+                                           Bed.edge_values,
+                                           Stage.boundary_values,
+                                           Xmom.boundary_values,
+                                           Ymom.boundary_values,
+                                           domain.boundary_flux_type,
+                                           Stage.explicit_update,
+                                           Xmom.explicit_update,
+                                           Ymom.explicit_update,
+                                           domain.already_computed_flux,
+                                           domain.max_speed,
+                                           int(domain.optimise_dry_cells),
+                                           Stage.centroid_values,
+                                           Bed.centroid_values, 
+                                           Bed.vertex_values)
+
+        #import pdb
+        #pdb.set_trace()
+
+        domain.flux_timestep = flux_timestep
+
+
+    def protect_against_infinitesimal_and_negative_heights(domain):
+        """protect against infinitesimal heights and associated high velocities"""
+
+        from swb2_domain_ext import protect
+        #print'using swb2_protect_against ..'
+
+        # shortcuts
+        wc = domain.quantities['stage'].centroid_values
+        wv = domain.quantities['stage'].vertex_values
+        zc = domain.quantities['elevation'].centroid_values
+        zv = domain.quantities['elevation'].vertex_values
+        xmomc = domain.quantities['xmomentum'].centroid_values
+        ymomc = domain.quantities['ymomentum'].centroid_values
+        areas = domain.areas
+
+        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
+                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
+
+    def conserved_values_to_evolved_values(self, q_cons, q_evol):
+        """Mapping between conserved quantities and the evolved quantities.
+        Used where we have a boundary condition which works with conserved
+        quantities and we now want to use them for the new well balanced
+        code using the evolved quantities
+
+        Typically the old boundary condition will set the values in q_cons,
+
+        q_evol on input will have the values of the evolved quantities at the
+        edge point (useful as it provides values for evlevation).
+
+        """
+
+        wc  = q_cons[0]
+        uhc = q_cons[1]
+        vhc = q_cons[2]
+
+
+        q_evol[0]  = wc
+        q_evol[1]  = uhc
+        q_evol[2]  = vhc
+
+        return q_evol
+
+    def distribute_to_vertices_and_edges(self):
+        """ Call correct module function """
+
+        if self.use_edge_limiter:
+            #distribute_using_edge_limiter(self)
+            self.distribute_using_edge_limiter()
+        else:
+            #distribute_using_vertex_limiter(self)
+            self.distribute_using_vertex_limiter()
+
+
+    def distribute_using_vertex_limiter(domain):
+        """Distribution from centroids to vertices specific to the SWW equation.
+
+        Precondition:
+          All quantities defined at centroids and bed elevation defined at
+          vertices.
+
+        Postcondition
+          Conserved quantities defined at vertices
+        """
+        # Remove very thin layers of water
+        domain.protect_against_infinitesimal_and_negative_heights()
+
+        # Extrapolate all conserved quantities
+        if domain.optimised_gradient_limiter:
+            # MH090605 if second order,
+            # perform the extrapolation and limiting on
+            # all of the conserved quantities
+
+            if (domain._order_ == 1):
+                for name in domain.conserved_quantities:
+                    Q = domain.quantities[name]
+                    Q.extrapolate_first_order()
+            elif domain._order_ == 2:
+                domain.extrapolate_second_order_sw()
+            else:
+                raise Exception('Unknown order')
+        else:
+            # Old code:
+            for name in domain.conserved_quantities:
+                Q = domain.quantities[name]
+
+                if domain._order_ == 1:
+                    Q.extrapolate_first_order()
+                elif domain._order_ == 2:
+                    Q.extrapolate_second_order_and_limit_by_vertex()
+                else:
+                    raise Exception('Unknown order')
+
+        # Take bed elevation into account when water heights are small
+        #balance_deep_and_shallow(domain)
+
+        # Compute edge values by interpolation
+        for name in domain.conserved_quantities:
+            Q = domain.quantities[name]
+            Q.interpolate_from_vertices_to_edges()
+
+
+    def distribute_using_edge_limiter(domain):
+        """Distribution from centroids to edges specific to the SWW eqn.
+
+        Precondition:
+          All quantities defined at centroids and bed elevation defined at
+          vertices.
+
+        Postcondition
+          Conserved quantities defined at vertices
+        """
+        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
+
+        # Remove very thin layers of water
+        #print 'Before protect'
+        domain.protect_against_infinitesimal_and_negative_heights()
+        #print 'After protect'
+
+        #for name in domain.conserved_quantities:
+        #    Q = domain.quantities[name]
+        #    if domain._order_ == 1:
+        #        Q.extrapolate_first_order()
+        #    elif domain._order_ == 2:
+        #        #Q.extrapolate_second_order_and_limit_by_edge()
+        #        # FIXME: This use of 'break' is hacky
+        #        domain.extrapolate_second_order_edge_sw()
+        #        break
+        #    else:
+        #        raise Exception('Unknown order')
+        
+        #print 'Before extrapolate'
+        domain.extrapolate_second_order_edge_sw()
+        #print 'After extrapolate'
+
+        #balance_deep_and_shallow(domain)
+
+        # Compute vertex values by interpolation
+        #for name in domain.conserved_quantities:
+        #    Q = domain.quantities[name]
+        #    Q.interpolate_from_edges_to_vertices()
+        #    Q.interpolate_from_vertices_to_edges()
+
+
+    def update_other_quantities(self):
+        """ There may be a need to calculates some of the other quantities
+        based on the new values of conserved quantities
+
+        But that is not needed in this version of the solver.
+        """
+        pass
+        # The centroid values of height and x and y velocity
+        # might not have been setup
+        
+        #self.update_centroids_of_velocities_and_height()
+        #
+        #for name in ['height', 'xvelocity', 'yvelocity']:
+        #    Q = self.quantities[name]
+        #    if self._order_ == 1:
+        #        Q.extrapolate_first_order()
+        #    elif self._order_ == 2:
+        #        Q.extrapolate_second_order_and_limit_by_edge()
+        #    else:
+        #        raise Exception('Unknown order')
+
+
+    def update_centroids_of_velocities_and_height(self):
+        """Calculate the centroid values of velocities and height based
+        on the values of the quantities stage and x and y momentum
+
+        Assumes that stage and momentum are up to date
+
+        Useful for kinematic viscosity calculations
+        """
+        pass
+        ## For shallow water we need to update height xvelocity and yvelocity
+
+        ##Shortcuts
+        #W  = self.quantities['stage']
+        #UH = self.quantities['xmomentum']
+        #VH = self.quantities['ymomentum']
+        #H  = self.quantities['height']
+        #Z  = self.quantities['elevation']
+        #U  = self.quantities['xvelocity']
+        #V  = self.quantities['yvelocity']
+
+        ##print num.min(W.centroid_values)
+
+        ## Make sure boundary values of conserved quantites
+        ## are consistent with value of functions at centroids
+        ##self.distribute_to_vertices_and_edges()
+        #Z.set_boundary_values_from_edges()
+
+        ##W.set_boundary_values_from_edges()
+        ##UH.set_boundary_values_from_edges()
+        ##VH.set_boundary_values_from_edges()
+
+        ## Update height values
+        #H.set_values(W.centroid_values-Z.centroid_values, location='centroids')
+        #H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0,
+        #                                 W.boundary_values-Z.boundary_values, 0.0))
+
+        ##assert num.min(H.centroid_values) >= 0
+        ##assert num.min(H.boundary_values) >= 0
+
+        ##Aliases
+        #uh_C  = UH.centroid_values
+        #vh_C  = VH.centroid_values
+        #h_C   = H.centroid_values
+
+        #uh_B  = UH.boundary_values
+        #vh_B  = VH.boundary_values
+        #h_B   = H.boundary_values
+
+        #H0 = 1.0e-8
+        #
+        #U.set_values(uh_C/(h_C + H0/h_C), location='centroids')
+        #V.set_values(vh_C/(h_C + H0/h_C), location='centroids')
+
+        #U.set_boundary_values(uh_B/(h_B + H0/h_B))
+        #V.set_boundary_values(vh_B/(h_B + H0/h_B))
+
+
+
+    def extrapolate_second_order_edge_sw(self):
+        """Wrapper calling C version of extrapolate_second_order_sw, using
+            edges instead of vertices in the extrapolation step. The routine
+            then interpolates from edges to vertices. 
+            The momentum extrapolation / interpolation is based on either 
+            momentum or velocity, depending on the choice of 
+            extrapolate_velocity_second_order.
+
+        self  the domain to operate on
+
+        """
+
+        from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2
+
+        # Shortcuts
+        Stage = self.quantities['stage']
+        Xmom = self.quantities['xmomentum']
+        Ymom = self.quantities['ymomentum']
+        Elevation = self.quantities['elevation']
+
+        extrapol2(self,
+                  self.surrogate_neighbours,
+                  self.number_of_boundaries,
+                  self.centroid_coordinates,
+                  Stage.centroid_values,
+                  Xmom.centroid_values,
+                  Ymom.centroid_values,
+                  Elevation.centroid_values,
+                  self.edge_coordinates,
+                  Stage.edge_values,
+                  Xmom.edge_values,
+                  Ymom.edge_values,
+                  Elevation.edge_values,
+                  Stage.vertex_values,
+                  Xmom.vertex_values,
+                  Ymom.vertex_values,
+                  Elevation.vertex_values,
+                  int(self.optimise_dry_cells),
+                  int(self.extrapolate_velocity_second_order))
+
+    def set_boundary(self, boundary_map):
+        ## Specialisation of set_boundary, which also updates the 'boundary_flux_type'
+        Sww_domain.set_boundary(self,boundary_map)
+        
+        # Add a flag which can be used to distinguish flux boundaries within
+        # compute_fluxes_central
+        # Initialise to zero (which means 'not a flux_boundary')
+        self.boundary_flux_type = self.boundary_edges*0
+       
+        # HACK to set the values of domain.boundary_flux 
+        for i in range(len(self.boundary_objects)):
+            # Record the first 10 characters of the name of the boundary.
+            # FIXME: There must be a better way than this!
+            bndry_name = self.boundary_objects[i][1].__repr__()[0:10]
+
+            if(bndry_name=='zero_mass_'):
+                # Create flag 'boundary_flux_type' that can be read in
+                # compute_fluxes_central
+                self.boundary_flux_type[i]=1
+
+
Index: /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain_ext.c
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain_ext.c	(revision 8523)
+++ /trunk/anuga_work/development/gareth/experimental/anuga_tsunami/swb2_domain_ext.c	(revision 8523)
@@ -0,0 +1,1787 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.6):
+//  gcc -c swb2_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O
+//  gcc -shared swb2_domain_ext.o  -o swb2_domain_ext.so
+//
+// or use python compile.py
+//
+// See the module swb_domain.py for more documentation on 
+// how to use this module
+//
+//
+// Stephen Roberts, ANU 2009
+// Ole Nielsen, GA 2004
+// Gareth Davies, GA 2011
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+//#include "numpy_shim.h"
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+const double pi = 3.14159265358979;
+
+// Computational function for rotation
+int _rotate(double *q, double n1, double n2) {
+  /*Rotate the last  2 coordinates of q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector (n1, n2).
+
+    Result is returned in array 2x1 r
+    To rotate in opposite direction, call rotate with (q, n1, -n2)
+
+    Contents of q are changed by this function */
+
+
+  double q1, q2;
+
+  // Shorthands
+  q1 = q[1];  // x coordinate
+  q2 = q[2];  // y coordinate
+
+  // Rotate
+  q[1] =  n1*q1 + n2*q2;
+  q[2] = -n2*q1 + n1*q2;
+
+  return 0;
+}
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
+//static inline double _compute_speed(double *uh, 
+double _compute_speed(double *uh, 
+		      double *h, 
+		      double epsilon, 
+		      double h0,
+		      double limiting_threshold) {
+  
+  double u;
+
+  if (*h < limiting_threshold) {   
+    // Apply limiting of speeds according to the ANUGA manual
+    if (*h < epsilon) {
+      //*h = max(0.0,*h);  // Could have been negative
+      *h = 0.0;  // Could have been negative
+      u = 0.0;
+    } else {
+      u = *uh/(*h + h0/ *h);    
+    }
+  
+
+    // Adjust momentum to be consistent with speed
+    *uh = u * *h;
+  } else {
+    // We are in deep water - no need for limiting
+    u = *uh/ *h;
+  }
+  
+  return u;
+}
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_central(double *q_left, double *q_right,
+                           double z_left, double z_right,
+                           double n1, double n2,
+                           double epsilon, 
+                           double h0,
+                           double limiting_threshold, 
+                           double g,
+                           double *edgeflux, double *max_speed, 
+                           double *pressure_flux) 
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', 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
+  */
+
+  int i;
+
+  double w_left, h_left, uh_left, vh_left, u_left;
+  double w_right, h_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, inverse_denominator, z;
+  double uint, t1, t2, t3;
+  // Workspace (allocate once, use many)
+  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
+
+  // Copy conserved quantities to protect from modification
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
+
+  // Align x- and y-momentum with x-axis
+  _rotate(q_left_rotated, n1, n2);
+  _rotate(q_right_rotated, n1, n2);
+
+  z = 0.5*(z_left + z_right); // Average elevation values. 
+                              // Even though this will nominally allow 
+			      // for discontinuities in the elevation data, 
+			      // there is currently no numerical support for
+			      // this so results may be strange near 
+			      // jumps in the bed.
+
+  // Compute speeds in x-direction
+  w_left = q_left_rotated[0];          
+  h_left = w_left - z;
+  uh_left = q_left_rotated[1];
+  u_left = _compute_speed(&uh_left, &h_left, 
+			  epsilon, h0, limiting_threshold);
+
+  w_right = q_right_rotated[0];
+  h_right = w_right - z;
+  uh_right = q_right_rotated[1];
+  u_right = _compute_speed(&uh_right, &h_right, 
+			   epsilon, h0, limiting_threshold);
+
+  // Momentum in y-direction
+  vh_left  = q_left_rotated[2];
+  vh_right = q_right_rotated[2];
+
+  // Limit y-momentum if necessary 
+  // Leaving this out, improves speed significantly (Ole 27/5/2009)
+  // All validation tests pass, so do we really need it anymore? 
+  _compute_speed(&vh_left, &h_left, 
+		 epsilon, h0, limiting_threshold);
+  _compute_speed(&vh_right, &h_right, 
+		 epsilon, h0, limiting_threshold);
+
+  // Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);  
+  
+  // Code to use fast square root optimisation if desired.
+  // Timings on AMD 64 for the Okushiri profile gave the following timings
+  //
+  // SQRT           Total    Flux
+  //=============================
+  //
+  // Ref            405s     152s
+  // Fast (dbl)     453s     173s
+  // Fast (sng)     437s     171s
+  //
+  // Consequently, there is currently (14/5/2009) no reason to use this 
+  // approximation.
+  
+  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
+  //soundspeed_right = fast_squareroot_approximation(g*h_right);
+
+  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
+  if (s_max < 0.0) 
+  {
+    s_max = 0.0;
+  }
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+  
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
+  flux_left[2] = u_left*vh_left;
+
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
+  flux_right[2] = u_right*vh_right;
+
+  // Flux computation
+  denom = s_max - s_min;
+  if (denom < epsilon) 
+  { // FIXME (Ole): Try using h0 here
+    memset(edgeflux, 0, 3*sizeof(double));
+    //for(i=0;i<3;i++){
+    //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
+    //}
+    *max_speed = 0.0;
+  } 
+  else 
+  {
+    inverse_denominator = 1.0/denom;
+    for (i = 0; i < 3; i++) 
+    {
+      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
+      // Physics 2:141-163
+      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
+      //t1 = (q_right_rotated[i] - uint);
+      //t2 = (-q_left_rotated[i] + uint);
+      //t3 = _minmod(t1,t2);
+      t3 = 0.0;
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
+      edgeflux[i] *= inverse_denominator;
+    }
+
+    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
+
+    // Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+
+    // Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+  
+  return 0;
+}
+
+
+// Computational function for flux computation
+double _compute_fluxes_central(int number_of_elements,
+        double timestep,
+        double epsilon,
+        double H0,
+        double g,
+        long* neighbours,
+        long* neighbour_edges,
+        double* normals,
+        double* edgelengths,
+        double* radii,
+        double* areas,
+        long* tri_full_flag,
+        double* stage_edge_values,
+        double* xmom_edge_values,
+        double* ymom_edge_values,
+        double* bed_edge_values,
+        double* stage_boundary_values,
+        double* xmom_boundary_values,
+        double* ymom_boundary_values,
+        long*   boundary_flux_type,
+        double* stage_explicit_update,
+        double* xmom_explicit_update,
+        double* ymom_explicit_update,
+        long* already_computed_flux,
+        double* max_speed_array,
+        int optimise_dry_cells, 
+        double* stage_centroid_values,
+        double* bed_centroid_values,
+        double* bed_vertex_values) {
+    // Local variables
+    double max_speed, length, inv_area, zl, zr;
+    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
+
+    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
+    // threshold for performance reasons.
+    // See ANUGA manual under flux limiting
+    int k, i, m, n,j;
+    int ki, nm = 0, ki2; // Index shorthands
+
+    // Workspace (making them static actually made function slightly slower (Ole))
+    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
+    double stage_edges[3];//Work array
+    double bedslope_work;
+    int neighbours_wet[3];//Work array
+    int useint;
+    double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n, tmp;
+    static long call = 1; // Static local variable flagging already computed flux
+
+    double *max_bed_edgevalue, *min_bed_edgevalue;
+
+    max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
+    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
+    // Start computation
+    call++; // Flag 'id' of flux calculation for this timestep
+
+    // Set explicit_update to zero for all conserved_quantities.
+    // This assumes compute_fluxes called before forcing terms
+    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
+
+
+    // Compute minimum bed edge value on each triangle
+    for (k = 0; k < number_of_elements; k++){
+        max_bed_edgevalue[k] = max(bed_edge_values[3*k], 
+                                   max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+        min_bed_edgevalue[k] = min(bed_edge_values[3*k], 
+                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+    
+    }
+
+
+    // For all triangles
+    for (k = 0; k < number_of_elements; k++) {
+        // Loop through neighbours and compute edge flux for each
+        for (i = 0; i < 3; i++) {
+            ki = k * 3 + i; // Linear index to edge i of triangle k
+
+            if (already_computed_flux[ki] == call) {
+                // We've already computed the flux across this edge
+                continue;
+            }
+
+            // Get left hand side values from triangle k, edge i
+            ql[0] = stage_edge_values[ki];
+            ql[1] = xmom_edge_values[ki];
+            ql[2] = ymom_edge_values[ki];
+            zl = bed_edge_values[ki];
+            hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.0);
+
+            // Get right hand side values either from neighbouring triangle
+            // or from boundary array (Quantities at neighbour on nearest face).
+            n = neighbours[ki];
+            hc_n = hc;
+            if (n < 0) {
+                // Neighbour is a boundary condition
+                m = -n - 1; // Convert negative flag to boundary index
+
+                qr[0] = stage_boundary_values[m];
+                qr[1] = xmom_boundary_values[m];
+                qr[2] = ymom_boundary_values[m];
+                zr = zl; // Extend bed elevation to boundary
+            }
+            else {
+                // Neighbour is a real triangle
+                hc_n = max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
+                m = neighbour_edges[ki];
+                nm = n * 3 + m; // Linear index (triangle n, edge m)
+
+                qr[0] = stage_edge_values[nm];
+                qr[1] = xmom_edge_values[nm];
+                qr[2] = ymom_edge_values[nm];
+                zr = bed_edge_values[nm];
+            }
+            
+
+            if (fabs(zl-zr)>1.0e-10) {
+                report_python_error(AT,"Discontinuous Elevation");
+                return 0.0;
+            }
+            
+            // If both centroids are dry, then the cells should not exchange flux 
+            // NOTE: This also means no bedslope term -- which is arguably
+            // appropriate, since the depth is zero at the cell centroid
+            if(( hc == 0.0)&(hc_n == 0.0)&(n>=0)){
+                already_computed_flux[ki] = call; // #k Done
+                already_computed_flux[nm] = call; // #n Done
+                max_speed = 0.0;
+                continue;
+            }
+           
+            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
+            // unless the local centroid value is smaller 
+            if(n>=0){
+                if(hc==0.0){
+                    //ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                    //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl);
+                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                }
+                if(hc_n==0.0){
+                    qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
+                    //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr);
+                    //qr[0]=ql[0]; 
+                }
+            }else{
+                // Treat the boundary case
+                //if((hc==0.0)){
+                //    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                    //ql[0]=max(min(qr[0],ql[0]),zl);
+                //}
+            }
+            
+            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
+            ki2 = 2 * ki; //k*6 + i*2
+
+            // Edge flux computation (triangle k, edge i)
+            _flux_function_central(ql, qr, zl, zr,
+                    normals[ki2], normals[ki2 + 1],
+                    epsilon, h0, limiting_threshold, g,
+                    edgeflux, &max_speed, &pressure_flux);
+
+            // Prevent outflow from 'seriously' dry cells
+            // Idea: The cell will not go dry if:
+            // mass_flux <= Area_triangle*hc/(dt*edgelength)
+            //if((stage_centroid_values[k]<=max_bed_edgevalue[k])|
+            //   (ql[0]<=zl)){
+            //    if(edgeflux[0]>0.0){
+            //        tmp=min(0.5*areas[k]*(hc+bed_centroid_values[k] - min_bed_edgevalue[k])/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
+            //        tmp = min(min(edgeflux[0], tmp)/edgeflux[0], 1.0);
+            //        edgeflux[0]*=tmp;
+            //    }
+            //}
+            //if(n>=0){
+            //    if((stage_centroid_values[n]<=max_bed_edgevalue[n])|
+            //       (qr[0]<=zr)){
+            //        if(edgeflux[0]<0.0){
+            //            tmp=min(0.5*areas[n]*(hc_n+bed_centroid_values[n] - min_bed_edgevalue[n])/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
+            //            tmp = min( max(edgeflux[0], -tmp)/edgeflux[0], 1.0);
+            //            edgeflux[0]*=tmp;
+            //        }
+            //    }
+            //} 
+
+            // Multiply edgeflux by edgelength
+            length = edgelengths[ki];
+            edgeflux[0] *= length;
+            edgeflux[1] *= length;
+            edgeflux[2] *= length;
+
+            // Update triangle k with flux from edge i
+            stage_explicit_update[k] -= edgeflux[0];
+            xmom_explicit_update[k] -= edgeflux[1];
+            ymom_explicit_update[k] -= edgeflux[2];
+
+            // Compute bed slope term
+            //if(hc>-9999.0){
+                //Bedslope approx 1:
+            bedslope_work = g*length*( hc*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
+                //
+                // Bedslope approx 2
+                //stage_edge_lim = ql[0]; // Limit this to be between a constant stage and constant depth extrapolation
+                //if(stage_edge_lim > max(stage_centroid_values[k], zl +hc)){
+                //    stage_edge_lim = max(stage_centroid_values[k], zl+hc);
+                //}
+                //if(stage_edge_lim < min(stage_centroid_values[k], zl +hc)){
+                //    stage_edge_lim = min(stage_centroid_values[k], zl+hc);
+                //}
+                //bedslope_work = g*hc*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zl,0.)*(stage_edge_lim-zl)*length;
+
+                // Bedslope approx 3
+                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
+                //
+            xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
+            ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
+            //}else{
+            //   // Treat nearly dry cells 
+            //   bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
+            //   //
+            //   //bedslope_work = -pressure_flux*length;
+            //   xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
+            //   ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
+            //   //xmom_explicit_update[k] = 0.0;
+            //   //ymom_explicit_update[k] = 0.0;
+
+            ////}
+
+            already_computed_flux[ki] = call; // #k Done
+
+
+            // Update neighbour n with same flux but reversed sign
+            if (n >= 0) {
+                stage_explicit_update[n] += edgeflux[0];
+                xmom_explicit_update[n] += edgeflux[1];
+                ymom_explicit_update[n] += edgeflux[2];
+                //Add bed slope term here
+                //if(hc_n>-9999.0){
+                //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
+                    // Bedslope approx 1:
+                bedslope_work = g*length*(hc_n*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
+                    //
+                    // Bedslope approx 2:
+                    //stage_edge_lim = qr[0];
+                    //if(stage_edge_lim > max(stage_centroid_values[n], zr +hc_n)){
+                    //    stage_edge_lim = max(stage_centroid_values[n], zr+hc_n);
+                    //}
+                    //if(stage_edge_lim < min(stage_centroid_values[n], zr +hc_n)){
+                    //    stage_edge_lim = min(stage_centroid_values[n], zr+hc_n);
+                    //}
+                    //bedslope_work = g*hc_n*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zr,0.)*(stage_edge_lim-zr)*length;
+                    //
+                    // Bedslope approx 3
+                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
+                    //
+                xmom_explicit_update[n] += normals[ki2]*bedslope_work;
+                ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
+                //}else{
+                //    // Treat nearly dry cells
+                //    bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
+                //    //bedslope_work = -pressure_flux*length;
+                //    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
+                //    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
+
+                //    //xmom_explicit_update[n] = 0.0;
+                //    //ymom_explicit_update[n] = 0.0;
+                //}
+
+                already_computed_flux[nm] = call; // #n Done
+            }
+
+            // Update timestep based on edge i and possibly neighbour n
+            if (tri_full_flag[k] == 1) {
+                if (max_speed > epsilon) {
+                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                    // CFL for triangle k
+                    timestep = min(timestep, radii[k] / max_speed);
+
+                    if (n >= 0) {
+                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                        timestep = min(timestep, radii[n] / max_speed);
+                    }
+
+                    // Ted Rigby's suggested less conservative version
+                    //if(n>=0){
+                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
+                    //}else{
+                    //  timestep = min(timestep, radii[k]/max_speed);
+                    //}
+                }
+            }
+
+        } // End edge i (and neighbour n)
+
+
+        // Normalise triangle k by area and store for when all conserved
+        // quantities get updated
+        inv_area = 1.0 / areas[k];
+        stage_explicit_update[k] *= inv_area;
+        xmom_explicit_update[k] *= inv_area;
+        ymom_explicit_update[k] *= inv_area;
+
+
+        // Keep track of maximal speeds
+        max_speed_array[k] = max_speed;
+
+    } // End triangle k
+
+    free(max_bed_edgevalue);
+    free(min_bed_edgevalue);
+
+    return timestep;
+}
+
+// Protect against the water elevation falling below the triangle bed
+int _protect(int N,
+         double minimum_allowed_height,
+         double maximum_allowed_speed,
+         double epsilon,
+         double* wc,
+         double* wv,
+         double* zc,
+         double* zv,
+         double* xmomc,
+         double* ymomc,
+         double* areas) {
+
+  int k;
+  double hc, bmin, bmax;
+  double u, v, reduced_speed;
+  static double mass_error=0.; 
+  // This acts like minimum_allowed height, but scales with the vertical
+  // distance between the bed_centroid_value and the max bed_edge_value of
+  // every triangle.
+  double minimum_relative_height=0.1; 
+  int mass_added=0;
+
+  // Protect against inifintesimal and negative heights  
+  //if (maximum_allowed_speed < epsilon) {
+    for (k=0; k<N; k++) {
+      hc = wc[k] - zc[k];
+      // Definine the maximum bed edge value on triangle k.
+      bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+
+      if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
+        
+        // Set momentum to zero and ensure h is non negative
+        // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
+        //if(hc<=epsilon){
+            xmomc[k] = 0.0;
+            ymomc[k] = 0.0;
+        //}
+
+        if (hc <= 0.0){
+             // Definine the minimum bed edge value on triangle k.
+             // Setting = minimum edge value can lead to mass conservation problems
+             //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
+             // Setting = minimum vertex value seems better, but might tend to be less smooth 
+             bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
+             //bmin=zc[k]-minimum_allowed_height;
+             // Minimum allowed stage = bmin
+             // WARNING: ADDING MASS if wc[k]<bmin
+             if(wc[k]<bmin){
+                 mass_error+=(bmin-wc[k])*areas[k];
+                 mass_added=1; //Flag to warn of added mass                
+                 //printf("Adding mass to dry cell %d %f %f %f %f %f \n", k, zv[3*k], zv[3*k+1], zv[3*k+2], wc[k]- bmin, mass_error);
+             
+                 wc[k] = max(wc[k], bmin); 
+             
+
+                 // Set vertex values as well. Seems that this shouldn't be
+                 // needed. However from memory this is important at the first
+                 // time step, for 'dry' areas where the designated stage is
+                 // less than the bed centroid value
+                 wv[3*k] = max(wv[3*k], bmin);
+                 wv[3*k+1] = max(wv[3*k+1], bmin);
+                 wv[3*k+2] = max(wv[3*k+2], bmin);
+            }
+        }
+      }
+    }
+
+  if(mass_added==1){
+     printf("Cumulative mass protection: %f m^3 \n", mass_error);
+  }
+  return 0;
+}
+
+int find_qmin_and_qmax(double dq0, double dq1, double dq2, 
+               double *qmin, double *qmax){
+  // Considering the centroid of an FV triangle and the vertices of its 
+  // auxiliary triangle, find 
+  // qmin=min(q)-qc and qmax=max(q)-qc, 
+  // where min(q) and max(q) are respectively min and max over the
+  // four values (at the centroid of the FV triangle and the auxiliary 
+  // triangle vertices),
+  // and qc is the centroid
+  // dq0=q(vertex0)-q(centroid of FV triangle)
+  // dq1=q(vertex1)-q(vertex0)
+  // dq2=q(vertex2)-q(vertex0)
+
+  // This is a simple implementation 
+  *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ;
+  *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ;
+ 
+  return 0;
+}
+
+int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
+  // Given provisional jumps dqv from the FV triangle centroid to its 
+  // vertices and jumps qmin (qmax) between the centroid of the FV 
+  // triangle and the minimum (maximum) of the values at the centroid of 
+  // the FV triangle and the auxiliary triangle vertices,
+  // calculate a multiplicative factor phi by which the provisional 
+  // vertex jumps are to be limited
+  
+  int i;
+  double r=1000.0, r0=1.0, phi=1.0;
+  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
+  // FIXME: Perhaps use the epsilon used elsewhere.
+  
+  // Any provisional jump with magnitude < TINY does not contribute to 
+  // the limiting process.
+  
+  for (i=0;i<3;i++){
+    if (dqv[i]<-TINY)
+      r0=qmin/dqv[i];
+      
+    if (dqv[i]>TINY)
+      r0=qmax/dqv[i];
+      
+    r=min(r0,r);
+  }
+  
+  phi=min(r*beta_w,1.0);
+  //phi=1.;
+  dqv[0]=dqv[0]*phi;
+  dqv[1]=dqv[1]*phi;
+  dqv[2]=dqv[2]*phi;
+    
+  return 0;
+}
+
+// Computational routine
+int _extrapolate_second_order_edge_sw(int number_of_elements,
+                                 double epsilon,
+                                 double minimum_allowed_height,
+                                 double beta_w,
+                                 double beta_w_dry,
+                                 double beta_uh,
+                                 double beta_uh_dry,
+                                 double beta_vh,
+                                 double beta_vh_dry,
+                                 long* surrogate_neighbours,
+                                 long* number_of_boundaries,
+                                 double* centroid_coordinates,
+                                 double* stage_centroid_values,
+                                 double* xmom_centroid_values,
+                                 double* ymom_centroid_values,
+                                 double* elevation_centroid_values,
+                                 double* edge_coordinates,
+                                 double* stage_edge_values,
+                                 double* xmom_edge_values,
+                                 double* ymom_edge_values,
+                                 double* elevation_edge_values,
+                                 double* stage_vertex_values,
+                                 double* xmom_vertex_values,
+                                 double* ymom_vertex_values,
+                                 double* elevation_vertex_values,
+                                 int optimise_dry_cells, 
+                                 int extrapolate_velocity_second_order) {
+                  
+  // Local variables
+  double a, b; // Gradient vector used to calculate edge values from centroids
+  int k, k0, k1, k2, k3, k6, coord_index, i, ii, ktmp;
+  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
+  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
+  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
+  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp;
+  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2;
+  
+  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
+  int *count_wet_neighbours;
+
+  // Use malloc to avoid putting these variables on the stack, which can cause
+  // segfaults in large model runs
+  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+  stage_centroid_store = malloc(number_of_elements*sizeof(double));
+  min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
+  max_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
+  count_wet_neighbours = malloc(number_of_elements*sizeof(int));
+ 
+  if(extrapolate_velocity_second_order==1){
+      // Replace momentum centroid with velocity centroid to allow velocity
+      // extrapolation This will be changed back at the end of the routine
+      for (k=0; k<number_of_elements; k++){
+
+          dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
+          xmom_centroid_store[k] = xmom_centroid_values[k];
+          xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
+
+          ymom_centroid_store[k] = ymom_centroid_values[k];
+          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
+
+          min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 
+                                           min(elevation_edge_values[3*k+1],
+                                               elevation_edge_values[3*k+2]));
+          max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 
+                                           max(elevation_edge_values[3*k+1],
+                                               elevation_edge_values[3*k+2]));
+          }
+
+      }
+
+  // Count how many 'fully submerged' neighbours the cell has
+  for(k=0; k<number_of_elements;k++){ 
+      count_wet_neighbours[k]=0;
+      for (i=0; i<3; i++){
+        ktmp = surrogate_neighbours[3*k+i];              
+        if(stage_centroid_values[ktmp] > max_elevation_edgevalue[ktmp]){
+            count_wet_neighbours[k]+=1;
+        }
+      }
+  }
+
+  // Begin extrapolation routine
+  for (k = 0; k < number_of_elements; k++) 
+  {
+    k3=k*3;
+    k6=k*6;
+
+    if (number_of_boundaries[k]==3)
+    //if (0==0)
+    {
+      // No neighbours, set gradient on the triangle to zero
+      
+      stage_edge_values[k3]   = stage_centroid_values[k];
+      stage_edge_values[k3+1] = stage_centroid_values[k];
+      stage_edge_values[k3+2] = stage_centroid_values[k];
+      xmom_edge_values[k3]    = xmom_centroid_values[k];
+      xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+      xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+      ymom_edge_values[k3]    = ymom_centroid_values[k];
+      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+      
+      continue;
+    }
+    else 
+    {
+      // Triangle k has one or more neighbours. 
+      // Get centroid and edge coordinates of the triangle
+      
+      // Get the edge coordinates
+      xv0 = edge_coordinates[k6];   
+      yv0 = edge_coordinates[k6+1];
+      xv1 = edge_coordinates[k6+2]; 
+      yv1 = edge_coordinates[k6+3];
+      xv2 = edge_coordinates[k6+4]; 
+      yv2 = edge_coordinates[k6+5];
+      
+      // Get the centroid coordinates
+      coord_index = 2*k;
+      x = centroid_coordinates[coord_index];
+      y = centroid_coordinates[coord_index+1];
+      
+      // Store x- and y- differentials for the edges of 
+      // triangle k relative to the centroid
+      dxv0 = xv0 - x; 
+      dxv1 = xv1 - x; 
+      dxv2 = xv2 - x;
+      dyv0 = yv0 - y; 
+      dyv1 = yv1 - y; 
+      dyv2 = yv2 - y;
+      // Compute the minimum distance from the centroid to an edge
+      //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
+      //demin=sqrt(demin);
+    }
+
+
+            
+    if (number_of_boundaries[k]<=1)
+    {
+      //==============================================
+      // Number of boundaries <= 1
+      //==============================================    
+    
+    
+      // If no boundaries, auxiliary triangle is formed 
+      // from the centroids of the three neighbours
+      // If one boundary, auxiliary triangle is formed 
+      // from this centroid and its two neighbours
+      
+      k0 = surrogate_neighbours[k3];
+      k1 = surrogate_neighbours[k3 + 1];
+      k2 = surrogate_neighbours[k3 + 2];
+     
+      // Take note if the max neighbour bed elevation is greater than the min
+      // neighbour stage -- suggests a 'steep' bed relative to the flow
+      bedmax = max(elevation_centroid_values[k], 
+                   max(elevation_centroid_values[k0],
+                       max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
+      //bedmax = elevation_centroid_values[k];
+      stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
+                     min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
+                         min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
+                             max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
+      if(stagemin < bedmax){
+         // This will cause first order extrapolation
+         k2 = k;
+         k0 = k;
+         k1 = k;
+      } 
+     
+      // Get the auxiliary triangle's vertex coordinates 
+      // (really the centroids of neighbouring triangles)
+      coord_index = 2*k0;
+      x0 = centroid_coordinates[coord_index];
+      y0 = centroid_coordinates[coord_index+1];
+      
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index+1];
+      
+      coord_index = 2*k2;
+      x2 = centroid_coordinates[coord_index];
+      y2 = centroid_coordinates[coord_index+1];
+
+      // Store x- and y- differentials for the vertices 
+      // of the auxiliary triangle
+      dx1 = x1 - x0; 
+      dx2 = x2 - x0;
+      dy1 = y1 - y0; 
+      dy2 = y2 - y0;
+      
+      // Calculate 2*area of the auxiliary triangle
+      // The triangle is guaranteed to be counter-clockwise      
+      area2 = dy2*dx1 - dy1*dx2; 
+     
+      // Treat triangles with zero or 1 wet neighbours. 
+      if ((area2 <= 0)) //|(count_wet_neighbours[k]==0))
+      {
+          //printf("Error negative triangle area \n");
+          //report_python_error(AT, "Negative triangle area");
+          //return -1;
+          stage_edge_values[k3]   = stage_centroid_values[k];
+          stage_edge_values[k3+1] = stage_centroid_values[k];
+          stage_edge_values[k3+2] = stage_centroid_values[k];
+          // First order momentum / velocity extrapolation
+          xmom_edge_values[k3]    = xmom_centroid_values[k];
+          xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+          xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+          ymom_edge_values[k3]    = ymom_centroid_values[k];
+          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+
+          continue;
+      }  
+      
+      // Calculate heights of neighbouring cells
+      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
+      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+      hmin = min(min(h0, min(h1, h2)), hc);
+
+      hfactor = 0.0;
+      //if (hmin > 0.001) 
+      if (hmin > 0.) 
+      //if (hc>0.0)
+      {
+        hfactor = 1.0 ;//hmin/(hmin + 0.004);
+        //hfactor=hmin/(hmin + 0.004);
+      }
+      
+    
+      //-----------------------------------
+      // stage
+      //-----------------------------------      
+      
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k
+      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle (centroids of neighbouring triangles)
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
+      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+      
+      inv_area2 = 1.0/area2;
+      // Calculate the gradient of stage on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+    
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+      
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+    
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+
+       
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------            
+
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k      
+      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
+      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+      
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited      
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+      
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      
+
+      for (i=0; i < 3; i++)
+      {
+        xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+      }
+      
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------                  
+
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k
+      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
+      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+      
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+       
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+      
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
+
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      
+      for (i=0;i<3;i++)
+      {
+        ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+      }
+      
+    } // End number_of_boundaries <=1 
+    else
+    {
+
+      //==============================================
+      // Number of boundaries == 2
+      //==============================================        
+        
+      // One internal neighbour and gradient is in direction of the neighbour's centroid
+      
+      // Find the only internal neighbour (k1?)
+      for (k2 = k3; k2 < k3 + 3; k2++)
+      {
+      // Find internal neighbour of triangle k      
+      // k2 indexes the edges of triangle k   
+      
+          if (surrogate_neighbours[k2] != k)
+          {
+             break;
+          }
+      }
+      
+      if ((k2 == k3 + 3)) 
+      {
+        // If we didn't find an internal neighbour
+        report_python_error(AT, "Internal neighbour not found");
+        return -1;
+      }
+      
+      k1 = surrogate_neighbours[k2];
+      
+      // The coordinates of the triangle are already (x,y). 
+      // Get centroid of the neighbour (x1,y1)
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index + 1];
+      
+      // Compute x- and y- distances between the centroid of 
+      // triangle k and that of its neighbour
+      dx1 = x1 - x; 
+      dy1 = y1 - y;
+      
+      // Set area2 as the square of the distance
+      area2 = dx1*dx1 + dy1*dy1;
+      
+      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) 
+      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+      // respectively correspond to the x- and y- gradients 
+      // of the conserved quantities
+      dx2 = 1.0/area2;
+      dy2 = dx2*dy1;
+      dx2 *= dx1;
+      
+      
+      //-----------------------------------
+      // stage
+      //-----------------------------------            
+
+      // Compute differentials
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      //for (i=0; i < 3; i++)
+      //{
+      stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+      //}
+      
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------                        
+      
+      // Compute differentials
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1 >= 0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient      
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      //for (i=0;i<3;i++)
+      //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
+      //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
+      //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
+      
+      for (i = 0; i < 3;i++)
+      {
+          xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+      }
+      
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------                        
+
+      // Compute differentials
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      } 
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      for (i=0;i<3;i++)
+              {
+              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+              }
+    } // else [number_of_boundaries==2]
+  } // for k=0 to number_of_elements-1
+  
+  // Compute vertex values of quantities
+  for (k=0; k<number_of_elements; k++){
+      k3=3*k;
+
+      // Compute stage vertex values 
+      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ; 
+      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1]; 
+      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2]; 
+      
+     // Compute xmom vertex values 
+      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ; 
+      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1]; 
+      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2]; 
+
+      // Compute ymom vertex values 
+      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ; 
+      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1]; 
+      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2]; 
+
+      // If needed, convert from velocity to momenta
+      if(extrapolate_velocity_second_order==1){
+          //Convert velocity back to momenta at centroids
+          xmom_centroid_values[k] = xmom_centroid_store[k];
+          ymom_centroid_values[k] = ymom_centroid_store[k];
+
+          // Re-compute momenta at edges
+          for (i=0; i<3; i++){
+              de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
+              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
+              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
+          }
+
+          // Re-compute momenta at vertices
+          for (i=0; i<3; i++){
+              de[i] = max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
+              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
+              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
+          }
+
+      }
+
+   
+  } 
+
+  free(xmom_centroid_store);
+  free(ymom_centroid_store);
+  free(stage_centroid_store);
+  free(min_elevation_edgevalue);
+  free(max_elevation_edgevalue);
+  free(count_wet_neighbours);
+  return 0;
+}           
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+
+//========================================================================
+// Compute fluxes
+//========================================================================
+
+// Modified central flux function
+
+PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function_central"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
+                     domain.H0,
+                                     domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     tri_full_flag,
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
+                                     Stage.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Stage.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update,
+                                     already_computed_flux,
+                                     optimise_dry_cells,
+                                     stage.centroid_values,
+                                     bed.centroid_values)
+
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
+  */
+
+
+  PyArrayObject *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *areas,
+    *tri_full_flag,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
+    *stage_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *boundary_flux_type,
+    *stage_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *max_speed_array, //Keeps track of max speeds for each triangle
+    *stage_centroid_values,
+    *bed_centroid_values,
+    *bed_vertex_values;
+    
+  double timestep, epsilon, H0, g;
+  int optimise_dry_cells;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOO",
+            &timestep,
+            &epsilon,
+            &H0,
+            &g,
+            &neighbours,
+            &neighbour_edges,
+            &normals,
+            &edgelengths, &radii, &areas,
+            &tri_full_flag,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &bed_edge_values,
+            &stage_boundary_values,
+            &xmom_boundary_values,
+            &ymom_boundary_values,
+            &boundary_flux_type,
+            &stage_explicit_update,
+            &xmom_explicit_update,
+            &ymom_explicit_update,
+            &already_computed_flux,
+            &max_speed_array,
+            &optimise_dry_cells,
+            &stage_centroid_values,
+            &bed_centroid_values,
+            &bed_vertex_values)) {
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(neighbours);
+  CHECK_C_CONTIG(neighbour_edges);
+  CHECK_C_CONTIG(normals);
+  CHECK_C_CONTIG(edgelengths);
+  CHECK_C_CONTIG(radii);
+  CHECK_C_CONTIG(areas);
+  CHECK_C_CONTIG(tri_full_flag);
+  CHECK_C_CONTIG(stage_edge_values);
+  CHECK_C_CONTIG(xmom_edge_values);
+  CHECK_C_CONTIG(ymom_edge_values);
+  CHECK_C_CONTIG(bed_edge_values);
+  CHECK_C_CONTIG(stage_boundary_values);
+  CHECK_C_CONTIG(xmom_boundary_values);
+  CHECK_C_CONTIG(ymom_boundary_values);
+  CHECK_C_CONTIG(boundary_flux_type);
+  CHECK_C_CONTIG(stage_explicit_update);
+  CHECK_C_CONTIG(xmom_explicit_update);
+  CHECK_C_CONTIG(ymom_explicit_update);
+  CHECK_C_CONTIG(already_computed_flux);
+  CHECK_C_CONTIG(max_speed_array);
+  CHECK_C_CONTIG(stage_centroid_values);
+  CHECK_C_CONTIG(bed_centroid_values);
+  CHECK_C_CONTIG(bed_vertex_values);
+
+  int number_of_elements = stage_edge_values -> dimensions[0];
+
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_central(number_of_elements,
+                     timestep,
+                     epsilon,
+                     H0,
+                     g,
+                     (long*) neighbours -> data,
+                     (long*) neighbour_edges -> data,
+                     (double*) normals -> data,
+                     (double*) edgelengths -> data, 
+                     (double*) radii -> data, 
+                     (double*) areas -> data,
+                     (long*) tri_full_flag -> data,
+                     (double*) stage_edge_values -> data,
+                     (double*) xmom_edge_values -> data,
+                     (double*) ymom_edge_values -> data,
+                     (double*) bed_edge_values -> data,
+                     (double*) stage_boundary_values -> data,
+                     (double*) xmom_boundary_values -> data,
+                     (double*) ymom_boundary_values -> data,
+                     (long*)   boundary_flux_type -> data,
+                     (double*) stage_explicit_update -> data,
+                     (double*) xmom_explicit_update -> data,
+                     (double*) ymom_explicit_update -> data,
+                     (long*) already_computed_flux -> data,
+                     (double*) max_speed_array -> data,
+                     optimise_dry_cells, 
+                     (double*) stage_centroid_values -> data, 
+                     (double*) bed_centroid_values -> data,
+                     (double*) bed_vertex_values -> data);
+  
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+PyObject *flux_function_central(PyObject *self, PyObject *args) {
+  //
+  // Gateway to innermost flux function.
+  // This is only used by the unit tests as the c implementation is
+  // normally called by compute_fluxes in this module.
+
+
+  PyArrayObject *normal, *ql, *qr,  *edgeflux;
+  double g, epsilon, max_speed, H0, zl, zr;
+  double h0, limiting_threshold, pressure_flux;
+
+  if (!PyArg_ParseTuple(args, "OOOddOddd",
+            &normal, &ql, &qr, &zl, &zr, &edgeflux,
+            &epsilon, &g, &H0)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+
+  
+  h0 = H0*H0; // This ensures a good balance when h approaches H0.
+              // But evidence suggests that h0 can be as little as
+              // epsilon!
+	      
+  limiting_threshold = 10*H0; // Avoid applying limiter below this
+                              // threshold for performance reasons.
+                              // See ANUGA manual under flux limiting  
+ 
+  pressure_flux = 0.0; // Store the water pressure related component of the flux 
+  _flux_function_central((double*) ql -> data, 
+			 (double*) qr -> data, 
+			 zl, 
+			 zr,                         
+			 ((double*) normal -> data)[0],
+			 ((double*) normal -> data)[1],          
+			 epsilon, h0, limiting_threshold,
+			 g,
+			 (double*) edgeflux -> data, 
+			 &max_speed,
+             &pressure_flux);
+  
+  return Py_BuildValue("d", max_speed);  
+}
+
+//========================================================================
+// Gravity
+//========================================================================
+
+PyObject *gravity(PyObject *self, PyObject *args) {
+  //
+  //  gravity(g, h, v, x, xmom, ymom)
+  //
+  
+  
+  PyArrayObject *h, *v, *x, *xmom, *ymom;
+  int k, N, k3, k6;
+  double g, avg_h, zx, zy;
+  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
+  //double epsilon;
+  
+  if (!PyArg_ParseTuple(args, "dOOOOO",
+			&g, &h, &v, &x,
+			&xmom, &ymom)) {
+    //&epsilon)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
+    return NULL;
+  }
+  
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(h);
+  CHECK_C_CONTIG(v);
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+  
+  N = h -> dimensions[0];
+  for (k=0; k<N; k++) {
+    k3 = 3*k;  // base index
+    
+    // Get bathymetry
+    z0 = ((double*) v -> data)[k3 + 0];
+    z1 = ((double*) v -> data)[k3 + 1];
+    z2 = ((double*) v -> data)[k3 + 2];
+    
+    // Optimise for flat bed
+    // Note (Ole): This didn't produce measurable speed up.
+    // Revisit later
+    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
+    //  continue;
+    //} 
+    
+    // Get average depth from centroid values
+    avg_h = ((double *) h -> data)[k];
+    
+    // Compute bed slope
+    k6 = 6*k;  // base index
+    
+    x0 = ((double*) x -> data)[k6 + 0];
+    y0 = ((double*) x -> data)[k6 + 1];
+    x1 = ((double*) x -> data)[k6 + 2];
+    y1 = ((double*) x -> data)[k6 + 3];
+    x2 = ((double*) x -> data)[k6 + 4];
+    y2 = ((double*) x -> data)[k6 + 5];
+    
+    
+    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+    
+    // Update momentum
+    ((double*) xmom -> data)[k] += -g*zx*avg_h;
+    ((double*) ymom -> data)[k] += -g*zy*avg_h;
+  }
+  
+  return Py_BuildValue("");
+}
+
+
+PyObject *extrapolate_second_order_edge_sw(PyObject *self, PyObject *args) {
+  /*Compute the edge values based on a linear reconstruction 
+    on each triangle
+    
+    Python call:
+    extrapolate_second_order_sw(domain.surrogate_neighbours,
+                                domain.number_of_boundaries
+                                domain.centroid_coordinates,
+                                Stage.centroid_values
+                                Xmom.centroid_values
+                                Ymom.centroid_values
+                                domain.edge_coordinates,
+                                Stage.edge_values,
+                                Xmom.edge_values,
+                                Ymom.edge_values)
+
+    Post conditions:
+        The edges of each triangle have values from a 
+        limited linear reconstruction
+        based on centroid values
+
+  */
+  PyArrayObject *surrogate_neighbours,
+    *number_of_boundaries,
+    *centroid_coordinates,
+    *stage_centroid_values,
+    *xmom_centroid_values,
+    *ymom_centroid_values,
+    *elevation_centroid_values,
+    *edge_coordinates,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *elevation_edge_values,
+    *stage_vertex_values,
+    *xmom_vertex_values,
+    *ymom_vertex_values,
+    *elevation_vertex_values;
+  
+  PyObject *domain;
+
+  
+  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;    
+  double minimum_allowed_height, epsilon;
+  int optimise_dry_cells, number_of_elements, extrapolate_velocity_second_order, e, e2;
+  
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOii",
+            &domain,
+            &surrogate_neighbours,
+            &number_of_boundaries,
+            &centroid_coordinates,
+            &stage_centroid_values,
+            &xmom_centroid_values,
+            &ymom_centroid_values,
+            &elevation_centroid_values,
+            &edge_coordinates,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &elevation_edge_values,
+            &stage_vertex_values,
+            &xmom_vertex_values,
+            &ymom_vertex_values,
+            &elevation_vertex_values,
+            &optimise_dry_cells,
+            &extrapolate_velocity_second_order)) {         
+
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(surrogate_neighbours);
+  CHECK_C_CONTIG(number_of_boundaries);
+  CHECK_C_CONTIG(centroid_coordinates);
+  CHECK_C_CONTIG(stage_centroid_values);
+  CHECK_C_CONTIG(xmom_centroid_values);
+  CHECK_C_CONTIG(ymom_centroid_values);
+  CHECK_C_CONTIG(elevation_centroid_values);
+  CHECK_C_CONTIG(edge_coordinates);
+  CHECK_C_CONTIG(stage_edge_values);
+  CHECK_C_CONTIG(xmom_edge_values);
+  CHECK_C_CONTIG(ymom_edge_values);
+  CHECK_C_CONTIG(elevation_edge_values);
+  CHECK_C_CONTIG(stage_vertex_values);
+  CHECK_C_CONTIG(xmom_vertex_values);
+  CHECK_C_CONTIG(ymom_vertex_values);
+  CHECK_C_CONTIG(elevation_vertex_values);
+  
+  // Get the safety factor beta_w, set in the config.py file. 
+  // This is used in the limiting process
+  
+
+  beta_w                 = get_python_double(domain,"beta_w");
+  beta_w_dry             = get_python_double(domain,"beta_w_dry");
+  beta_uh                = get_python_double(domain,"beta_uh");
+  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
+  beta_vh                = get_python_double(domain,"beta_vh");
+  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");  
+
+  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
+  epsilon                = get_python_double(domain,"epsilon");
+
+  number_of_elements = stage_centroid_values -> dimensions[0];  
+
+  //printf("In C before Extrapolate");
+  //e=1;
+  // Call underlying computational routine
+  e = _extrapolate_second_order_edge_sw(number_of_elements,
+                   epsilon,
+                   minimum_allowed_height,
+                   beta_w,
+                   beta_w_dry,
+                   beta_uh,
+                   beta_uh_dry,
+                   beta_vh,
+                   beta_vh_dry,
+                   (long*) surrogate_neighbours -> data,
+                   (long*) number_of_boundaries -> data,
+                   (double*) centroid_coordinates -> data,
+                   (double*) stage_centroid_values -> data,
+                   (double*) xmom_centroid_values -> data,
+                   (double*) ymom_centroid_values -> data,
+                   (double*) elevation_centroid_values -> data,
+                   (double*) edge_coordinates -> data,
+                   (double*) stage_edge_values -> data,
+                   (double*) xmom_edge_values -> data,
+                   (double*) ymom_edge_values -> data,
+                   (double*) elevation_edge_values -> data,
+                   (double*) stage_vertex_values -> data,
+                   (double*) xmom_vertex_values -> data,
+                   (double*) ymom_vertex_values -> data,
+                   (double*) elevation_vertex_values -> data,
+                   optimise_dry_cells, 
+                   extrapolate_velocity_second_order);
+
+  if (e == -1) {
+    // Use error string set inside computational routine
+    return NULL;
+  }                   
+  
+  
+  return Py_BuildValue("");
+  
+}// extrapolate_second-order_edge_sw
+//========================================================================
+// Protect -- to prevent the water level from falling below the minimum 
+// bed_edge_value
+//========================================================================
+
+PyObject *protect(PyObject *self, PyObject *args) {
+  //
+  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
+
+
+  PyArrayObject
+  *wc,            //Stage at centroids
+  *wv,            //Stage at vertices
+  *zc,            //Elevation at centroids
+  *zv,            //Elevation at vertices
+  *xmomc,         //Momentums at centroids
+  *ymomc,
+  *areas;         // Triangle areas
+
+
+  int N;
+  double minimum_allowed_height, maximum_allowed_speed, epsilon;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddOOOOOOO",
+            &minimum_allowed_height,
+            &maximum_allowed_speed,
+            &epsilon,
+            &wc, &wv, &zc, &zv, &xmomc, &ymomc, &areas)) {
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }  
+
+  N = wc -> dimensions[0];
+
+  _protect(N,
+       minimum_allowed_height,
+       maximum_allowed_speed,
+       epsilon,
+       (double*) wc -> data,
+       (double*) wv -> data,
+       (double*) zc -> data,
+       (double*) zv -> data,
+       (double*) xmomc -> data,
+       (double*) ymomc -> data,
+       (double*) areas -> data);
+
+  return Py_BuildValue("");
+}
+
+//========================================================================
+// 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"},
+  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
+  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
+  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
+  {"extrapolate_second_order_edge_sw", extrapolate_second_order_edge_sw, METH_VARARGS, "Print out"},
+  {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {NULL, NULL, 0, NULL}
+};
+
+// Module initialisation
+void initswb2_domain_ext(void){
+  Py_InitModule("swb2_domain_ext", MethodTable);
+
+  import_array(); // Necessary for handling of NumPY structures
+}
Index: /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py	(revision 8522)
+++ /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py	(revision 8523)
@@ -54,2 +54,27 @@
         pyplot.savefig(fig_filename)
         pyplot.clf()
+
+
+### Plot of experimental set-up
+pyplot.figure(figsize=(8,3))
+pyplot.scatter(p.x, p.y, c=p.elev>0.2,edgecolors='none', s=5.0, marker='8')
+pyplot.gca().set_aspect('equal')
+pyplot.xlim((4,15))
+pyplot.title('Zoom of the flume geometry')
+#pyplot.colorbar()
+# Add gauge locations
+for i in range(Gauge_xy.shape[0]):
+    gauge_name='G' + str(i+1)
+    pyplot.scatter(Gauge_xy[i,0], Gauge_xy[i,1],marker='x',color='white')
+    pyplot.text(Gauge_xy[i,0], Gauge_xy[i,1], s=gauge_name, fontweight='black', color='white', ha='left')
+
+pyplot.text(7.0, 1.7, 'Reservoir \n Outflow',color='white', size=8)
+pyplot.quiver(7,1.5, 1.0, 0.0,color='white')
+pyplot.text(4.5, 1.0, 'Reservoir',color='white')
+
+pyplot.text(11.7, 1.5, 'Isolated \n Building',color='white',size=8)
+#pyplot.quiver(12.5,1.5, -1.0, 1.0, color='white')
+
+fig_filename=fig_dir+'Model_domain.png'
+pyplot.savefig(fig_filename, orientation='landscape')
+pyplot.clf()
