Index: uga_core/source/anuga/abstract_2d_finite_volumes/shallow_water_kinetic.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/shallow_water_kinetic.py	(revision 3697)
+++ 	(revision )
@@ -1,1815 +1,0 @@
-"""Class Domain -
-2D triangular domains for finite-volume computations of
-the shallow water wave equation.
-
-This module contains a specialisation of class Domain from module domain.py
-consisting of methods specific to the Shallow Water Wave Equation
-
-
-U_t + E_x + G_y = S
-
-where
-
-U = [w, uh, vh]
-E = [uh, u^2h + gh^2/2, uvh]
-G = [vh, uvh, v^2h + gh^2/2]
-S represents source terms forcing the system
-(e.g. gravity, friction, wind stress, ...)
-
-and _t, _x, _y denote the derivative with respect to t, x and y respectively.
-
-The quantities are
-
-symbol    variable name    explanation
-x         x                horizontal distance from origin [m]
-y         y                vertical distance from origin [m]
-z         elevation        elevation of bed on which flow is modelled [m]
-h         height           water height above z [m]
-w         stage            absolute water level, w = z+h [m]
-u                          speed in the x direction [m/s]
-v                          speed in the y direction [m/s]
-uh        xmomentum        momentum in the x direction [m^2/s]
-vh        ymomentum        momentum in the y direction [m^2/s]
-
-eta                        mannings friction coefficient [to appear]
-nu                         wind stress coefficient [to appear]
-
-The conserved quantities are w, uh, vh
-
-
-For details see e.g.
-Christopher Zoppou and Stephen Roberts,
-Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
-Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
-
-
-Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-Geoscience Australia, 2004
-"""
-
-#Subversion keywords:
-#
-#$LastChangedDate: 2006-01-13 16:43:01 +1100 (Fri, 13 Jan 2006) $
-#$LastChangedRevision: 2206 $
-#$LastChangedBy: steve $
-
-
-from domain import *
-from region import *
-
-Generic_Domain = Domain #Rename
-
-#Shallow water domain
-class Domain(Generic_Domain):
-
-    def __init__(self, coordinates, vertices, boundary = None,
-                 tagged_elements = None, geo_reference = None,
-                 use_inscribed_circle=False):
-
-        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
-        other_quantities = ['elevation', 'friction']
-        Generic_Domain.__init__(self, coordinates, vertices, boundary,
-                                conserved_quantities, other_quantities,
-                                tagged_elements, geo_reference, use_inscribed_circle)
-
-        from anuga.config import minimum_allowed_height, g
-        self.minimum_allowed_height = minimum_allowed_height
-        self.g = g
-
-        self.forcing_terms.append(gravity)
-        self.forcing_terms.append(manning_friction)
-
-        #Realtime visualisation
-        self.visualiser = None
-        self.visualise  = False
-        self.visualise_color_stage = False
-        self.visualise_stage_range = 1.0
-        self.visualise_timer = True
-        self.visualise_range_z = None
-
-        #Stored output
-        self.store = False
-        self.format = 'sww'
-        self.smooth = True
-
-        #Reduction operation for get_vertex_values
-        from anuga.pyvolution.util import mean
-        self.reduction = mean
-        #self.reduction = min  #Looks better near steep slopes
-
-        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
-
-
-
-    def initialise_visualiser(self,scale_z=1.0,rect=None):
-        #Realtime visualisation
-        if self.visualiser is None:
-            from realtime_visualisation_new import Visualiser
-            self.visualiser = Visualiser(self,scale_z,rect)
-            self.visualiser.setup['elevation']=True
-            self.visualiser.updating['stage']=True
-        self.visualise = True
-        if self.visualise_color_stage == True:
-            self.visualiser.coloring['stage'] = True
-            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
-        print 'initialise visualiser'
-        print self.visualiser.setup
-        print self.visualiser.updating
-
-    def check_integrity(self):
-        Generic_Domain.check_integrity(self)
-
-        #Check that we are solving the shallow water wave equation
-
-        msg = 'First conserved quantity must be "stage"'
-        assert self.conserved_quantities[0] == 'stage', msg
-        msg = 'Second conserved quantity must be "xmomentum"'
-        assert self.conserved_quantities[1] == 'xmomentum', msg
-        msg = 'Third conserved quantity must be "ymomentum"'
-        assert self.conserved_quantities[2] == 'ymomentum', msg
-
-    def extrapolate_second_order_sw(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        extrapolate_second_order_sw(self)
-
-    def compute_fluxes(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        compute_fluxes(self)
-
-    def distribute_to_vertices_and_edges(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        distribute_to_vertices_and_edges(self)
-
-
-    #FIXME: Under construction
-#     def set_defaults(self):
-#         """Set default values for uninitialised quantities.
-#         This is specific to the shallow water wave equation
-#         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum'
-#         are 0.0. Default for 'stage' is whatever the value of 'elevation'.
-#         """
-
-#         for name in self.other_quantities + self.conserved_quantities:
-#             print name
-#             print self.quantities.keys()
-#             if not self.quantities.has_key(name):
-#                 if name == 'stage':
-
-#                     if self.quantities.has_key('elevation'):
-#                         z = self.quantities['elevation'].vertex_values
-#                         self.set_quantity(name, z)
-#                     else:
-#                         self.set_quantity(name, 0.0)
-#                 else:
-#                     self.set_quantity(name, 0.0)
-
-
-
-#         #Lift negative heights up
-#         #z = self.quantities['elevation'].vertex_values
-#         #w = self.quantities['stage'].vertex_values
-
-#         #h = w-z
-
-#         #for k in range(h.shape[0]):
-#         #    for i in range(3):
-#         #        if h[k, i] < 0.0:
-#         #            w[k, i] = z[k, i]
-
-
-#         #self.quantities['stage'].interpolate()
-
-
-    def evolve(self, yieldstep = None, finaltime = None,
-               skip_initial_step = False):
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call check integrity here rather than from user scripts
-        #self.check_integrity()
-
-        msg = 'Parameter beta_h must be in the interval [0, 1['
-        assert 0 <= self.beta_h < 1.0, msg
-        msg = 'Parameter beta_w must be in the interval [0, 1['
-        assert 0 <= self.beta_w < 1.0, msg
-
-
-        #Initial update of vertex and edge values before any storage
-        #and or visualisation
-        self.distribute_to_vertices_and_edges()
-
-        #Initialise real time viz if requested
-        if self.visualise is True and self.time == 0.0:
-            if self.visualiser is None:
-                self.initialise_visualiser()
-
-            self.visualiser.update_timer()
-            self.visualiser.setup_all()
-
-        #Store model data, e.g. for visualisation
-        if self.store is True and self.time == 0.0:
-            self.initialise_storage()
-            #print 'Storing results in ' + self.writer.filename
-        else:
-            pass
-            #print 'Results will not be stored.'
-            #print 'To store results set domain.store = True'
-            #FIXME: Diagnostic output should be controlled by
-            # a 'verbose' flag living in domain (or in a parent class)
-
-        #Call basic machinery from parent class
-        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
-                                       skip_initial_step):
-            #Real time viz
-            if self.visualise is True:
-                self.visualiser.update_all()
-                self.visualiser.update_timer()
-
-
-            #Store model data, e.g. for subsequent visualisation
-            if self.store is True:
-                self.store_timestep(self.quantities_to_be_stored)
-
-            #FIXME: Could maybe be taken from specified list
-            #of 'store every step' quantities
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)
-
-    def initialise_storage(self):
-        """Create and initialise self.writer object for storing data.
-        Also, save x,y and bed elevation
-        """
-
-        import anuga.pyvolution.data_manager
-
-        #Initialise writer
-        self.writer = data_manager.get_dataobject(self, mode = 'w')
-
-        #Store vertices and connectivity
-        self.writer.store_connectivity()
-
-
-    def store_timestep(self, name):
-        """Store named quantity and time.
-
-        Precondition:
-           self.write has been initialised
-        """
-        self.writer.store_timestep(name)
-
-
-#=============== End of Shallow Water Domain ===============================
-
-
-
-#Rotation of momentum vector
-def rotate(q, normal, direction = 1):
-    """Rotate the momentum component q (q[1], q[2])
-    from x,y coordinates to coordinates based on normal vector.
-
-    If direction is negative the rotation is inverted.
-
-    Input vector is preserved
-
-    This function is specific to the shallow water wave equation
-    """
-
-    from Numeric import zeros, Float
-
-    assert len(q) == 3,\
-           'Vector of conserved quantities must have length 3'\
-           'for 2D shallow water equation'
-
-    try:
-        l = len(normal)
-    except:
-        raise 'Normal vector must be an Numeric array'
-
-    assert l == 2, 'Normal vector must have 2 components'
-
-
-    n1 = normal[0]
-    n2 = normal[1]
-
-    r = zeros(len(q), Float) #Rotated quantities
-    r[0] = q[0]              #First quantity, height, is not rotated
-
-    if direction == -1:
-        n2 = -n2
-
-
-    r[1] =  n1*q[1] + n2*q[2]
-    r[2] = -n2*q[1] + n1*q[2]
-
-    return r
-
-
-
-####################################
-# Flux computation
-def flux_function_central(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
-    in the numerical vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from anuga.config import g, epsilon
-    from math import sqrt
-    from Numeric import array
-
-
-    #Align momentums with x-axis
-    q_left  = rotate(ql, normal, direction = 1)
-    q_right = rotate(qr, normal, direction = 1)
-
-    z = (zl+zr)/2 #Take average of field values
-
-    w_left  = q_left[0]   #w=h+z
-    h_left  = w_left-z
-    uh_left = q_left[1]
-
-    if h_left < epsilon:
-        u_left = 0.0  #Could have been negative
-        h_left = 0.0
-    else:
-        u_left  = uh_left/h_left
-
-
-    w_right  = q_right[0]  #w=h+z
-    h_right  = w_right-z
-    uh_right = q_right[1]
-
-
-    if h_right < epsilon:
-        u_right = 0.0  #Could have been negative
-        h_right = 0.0
-    else:
-        u_right  = uh_right/h_right
-
-    vh_left  = q_left[2]
-    vh_right = q_right[2]
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-
-    #Flux computation
-
-    #FIXME(Ole): Why is it again that we don't
-    #use uh_left and uh_right directly in the first entries?
-    flux_left  = array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left**2,
-                        u_left*vh_left])
-    flux_right = array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right**2,
-                        u_right*vh_right])
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = array([0.0, 0.0, 0.0])
-        max_speed = 0.0
-    else:
-        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
-        edgeflux += s_max*s_min*(q_right-q_left)/denom
-
-        edgeflux = rotate(edgeflux, normal, direction=-1)
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
-def erfcc(x):
-
-    from math import fabs, exp
-
-    z=fabs(x)
-    t=1.0/(1.0+0.5*z)
-    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
-         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
-         t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
-    if x < 0.0:
-        result = 2.0-result
-
-    return result
-
-def flux_function_kinetic(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
-
-
-    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
-    in the numerical vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from anuga.config import g, epsilon
-    from math import sqrt
-    from Numeric import array
-
-    #Align momentums with x-axis
-    q_left  = rotate(ql, normal, direction = 1)
-    q_right = rotate(qr, normal, direction = 1)
-
-    z = (zl+zr)/2 #Take average of field values
-
-    w_left  = q_left[0]   #w=h+z
-    h_left  = w_left-z
-    uh_left = q_left[1]
-
-    if h_left < epsilon:
-        u_left = 0.0  #Could have been negative
-        h_left = 0.0
-    else:
-        u_left  = uh_left/h_left
-
-
-    w_right  = q_right[0]  #w=h+z
-    h_right  = w_right-z
-    uh_right = q_right[1]
-
-
-    if h_right < epsilon:
-        u_right = 0.0  #Could have been negative
-        h_right = 0.0
-    else:
-        u_right  = uh_right/h_right
-
-    vh_left  = q_left[2]
-    vh_right = q_right[2]
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-
-    #Flux computation
-
-    F_left  = 0.0
-    F_right = 0.0
-    from math import sqrt, pi, exp
-    if h_left > 0.0:
-        F_left = u_left/sqrt(g*h_left)
-    if h_right > 0.0:
-        F_right = u_right/sqrt(g*h_right)
-
-    edgeflux = array([0.0, 0.0, 0.0])
-
-    edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
-          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
-          h_right*u_right/2.0*erfcc(F_right) -  \
-          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
-
-    edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \
-          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
-          (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) -  \
-          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
-
-    edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
-          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
-          vh_right*u_right/2.0*erfcc(F_right) - \
-          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
-
-
-    edgeflux = rotate(edgeflux, normal, direction=-1)
-    max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
-flux_function = flux_function_kinetic
-
-def compute_fluxes(domain):
-    """Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    """
-
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    #Arrays
-    stage = Stage.edge_values
-    xmom =  Xmom.edge_values
-    ymom =  Ymom.edge_values
-    bed =   Bed.edge_values
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-    ymom_bdry =  Ymom.boundary_values
-
-    flux = zeros(3, Float) #Work array for summing up fluxes
-
-
-    #Loop
-    timestep = float(sys.maxint)
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work array
-        for i in range(3):
-            #Quantities inside volume facing neighbour i
-            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
-            zl = bed[k, i]
-
-            #Quantities at neighbour on nearest face
-            n = domain.neighbours[k,i]
-            if n < 0:
-                m = -n-1 #Convert negative flag to index
-                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                m = domain.neighbour_edges[k,i]
-                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, 2*i:2*i+2]
-
-            #Flux computation using provided function
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            flux -= edgeflux * domain.edgelengths[k,i]
-
-            #Update optimal_timestep
-            try:
-                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-            except ZeroDivisionError:
-                pass
-
-        #Normalise by area and store for when all conserved
-        #quantities get updated
-        flux /= domain.areas[k]
-        Stage.explicit_update[k] = flux[0]
-        Xmom.explicit_update[k] = flux[1]
-        Ymom.explicit_update[k] = flux[2]
-
-
-    domain.timestep = timestep
-
-#MH090605 The following method belongs to the shallow_water domain class
-#see comments in the corresponding method in shallow_water_ext.c
-def extrapolate_second_order_sw_c(domain):
-    """Wrapper calling C version of extrapolate_second_order_sw
-    """
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Ymom = domain.quantities['ymomentum']
-    from shallow_water_ext import extrapolate_second_order_sw
-    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
-                                domain.number_of_boundaries,
-                                domain.centroid_coordinates,
-                                Stage.centroid_values,
-                                Xmom.centroid_values,
-                                Ymom.centroid_values,
-                                domain.vertex_coordinates,
-                                Stage.vertex_values,
-                                Xmom.vertex_values,
-                                Ymom.vertex_values)
-
-def compute_fluxes_c(domain):
-    """Wrapper calling C version of compute fluxes
-    """
-
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    timestep = float(sys.maxint)
-    from shallow_water_ext import compute_fluxes
-    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     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,
-                                     domain.already_computed_flux)
-
-
-####################################
-# Module functions for gradient limiting (distribute_to_vertices_and_edges)
-
-def distribute_to_vertices_and_edges(domain):
-    """Distribution from centroids to vertices specific to the
-    shallow water wave
-    equation.
-
-    It will ensure that h (w-z) is always non-negative even in the
-    presence of steep bed-slopes by taking a weighted average between shallow
-    and deep cases.
-
-    In addition, all conserved quantities get distributed as per either a
-    constant (order==1) or a piecewise linear function (order==2).
-
-    FIXME: more explanation about removal of artificial variability etc
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
-
-    """
-
-    from anuga.config import optimised_gradient_limiter
-
-    #Remove very thin layers of water
-    protect_against_infinitesimal_and_negative_heights(domain)
-
-    #Extrapolate all conserved quantities
-    if optimised_gradient_limiter:
-        #MH090605 if second order,
-        #perform the extrapolation and limiting on
-        #all of the conserved quantitie
-
-        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 '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()
-                Q.limit()
-            else:
-                raise '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 protect_against_infinitesimal_and_negative_heights(domain):
-    """Protect against infinitesimal heights and associated high velocities
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values
-    hc = wc - zc  #Water depths at centroids
-
-    #Update
-    for k in range(domain.number_of_elements):
-
-        if hc[k] < domain.minimum_allowed_height:
-            #Control stage
-            if hc[k] < domain.epsilon:
-                wc[k] = zc[k] # Contain 'lost mass' error
-
-            #Control momentum
-            xmomc[k] = ymomc[k] = 0.0
-
-
-def protect_against_infinitesimal_and_negative_heights_c(domain):
-    """Protect against infinitesimal heights and associated high velocities
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values
-
-    from shallow_water_ext import protect
-
-    protect(domain.minimum_allowed_height, domain.epsilon,
-            wc, zc, xmomc, ymomc)
-
-
-
-def h_limiter(domain):
-    """Limit slopes for each volume to eliminate artificial variance
-    introduced by e.g. second order extrapolator
-
-    limit on h = w-z
-
-    This limiter depends on two quantities (w,z) so it resides within
-    this module rather than within quantity.py
-    """
-
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-    beta_h = domain.beta_h
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-    hvbar = zeros(hv.shape, Float) #h-limited values
-
-    #Find min and max of this and neighbour's centroid values
-    hmax = zeros(hc.shape, Float)
-    hmin = zeros(hc.shape, Float)
-
-    for k in range(N):
-        hmax[k] = hmin[k] = hc[k]
-        for i in range(3):
-            n = domain.neighbours[k,i]
-            if n >= 0:
-                hn = hc[n] #Neighbour's centroid value
-
-                hmin[k] = min(hmin[k], hn)
-                hmax[k] = max(hmax[k], hn)
-
-
-    #Diffences between centroids and maxima/minima
-    dhmax = hmax - hc
-    dhmin = hmin - hc
-
-    #Deltas between vertex and centroid values
-    dh = zeros(hv.shape, Float)
-    for i in range(3):
-        dh[:,i] = hv[:,i] - hc
-
-    #Phi limiter
-    for k in range(N):
-
-        #Find the gradient limiter (phi) across vertices
-        phi = 1.0
-        for i in range(3):
-            r = 1.0
-            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
-            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
-
-            phi = min( min(r*beta_h, 1), phi )
-
-        #Then update using phi limiter
-        for i in range(3):
-            hvbar[k,i] = hc[k] + phi*dh[k,i]
-
-    return hvbar
-
-
-
-def h_limiter_c(domain):
-    """Limit slopes for each volume to eliminate artificial variance
-    introduced by e.g. second order extrapolator
-
-    limit on h = w-z
-
-    This limiter depends on two quantities (w,z) so it resides within
-    this module rather than within quantity.py
-
-    Wrapper for c-extension
-    """
-
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-    beta_h = domain.beta_h
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv - zv
-
-    #Call C-extension
-    from shallow_water_ext import h_limiter_sw as h_limiter
-    hvbar = h_limiter(domain, hc, hv)
-
-    return hvbar
-
-
-def balance_deep_and_shallow(domain):
-    """Compute linear combination between stage as computed by
-    gradient-limiters limiting using w, and stage computed by
-    gradient-limiters limiting using h (h-limiter).
-    The former takes precedence when heights are large compared to the
-    bed slope while the latter takes precedence when heights are
-    relatively small.  Anything in between is computed as a balanced
-    linear combination in order to avoid numerical disturbances which
-    would otherwise appear as a result of hard switching between
-    modes.
-
-    The h-limiter is always applied irrespective of the order.
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-    #Limit h
-    hvbar = h_limiter(domain)
-
-    for k in range(domain.number_of_elements):
-        #Compute maximal variation in bed elevation
-        #  This quantitiy is
-        #    dz = max_i abs(z_i - z_c)
-        #  and it is independent of dimension
-        #  In the 1d case zc = (z0+z1)/2
-        #  In the 2d case zc = (z0+z1+z2)/3
-
-        dz = max(abs(zv[k,0]-zc[k]),
-                 abs(zv[k,1]-zc[k]),
-                 abs(zv[k,2]-zc[k]))
-
-
-        hmin = min( hv[k,:] )
-
-        #Create alpha in [0,1], where alpha==0 means using the h-limited
-        #stage and alpha==1 means using the w-limited stage as
-        #computed by the gradient limiter (both 1st or 2nd order)
-
-        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
-        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-
-        if dz > 0.0:
-            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
-        else:
-            #Flat bed
-            alpha = 1.0
-
-	#Let
-	#
-	#  wvi be the w-limited stage (wvi = zvi + hvi)
-	#  wvi- be the h-limited state (wvi- = zvi + hvi-)
-	#
-	#
-	#where i=0,1,2 denotes the vertex ids
-	#
-        #Weighted balance between w-limited and h-limited stage is
-        #
-	#  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
-	#
-        #It follows that the updated wvi is
-        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
-        #
-	# Momentum is balanced between constant and limited
-
-
-        #for i in range(3):
-        #    wv[k,i] = zv[k,i] + hvbar[k,i]
-
-	#return
-
-	if alpha < 1:
-
-            for i in range(3):
-                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
-
-            #Momentums at centroids
-            xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values
-
-            #Momentums at vertices
-            xmomv = domain.quantities['xmomentum'].vertex_values
-            ymomv = domain.quantities['ymomentum'].vertex_values
-
-            # Update momentum as a linear combination of
-            # xmomc and ymomc (shallow) and momentum
-            # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
-            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-
-
-def balance_deep_and_shallow_c(domain):
-    """Wrapper for C implementation
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv - zv
-
-    #Momentums at centroids
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values
-
-    #Momentums at vertices
-    xmomv = domain.quantities['xmomentum'].vertex_values
-    ymomv = domain.quantities['ymomentum'].vertex_values
-
-    #Limit h
-    hvbar = h_limiter(domain)
-
-    #This is how one would make a first order h_limited value
-    #as in the old balancer (pre 17 Feb 2005):
-    #from Numeric import zeros, Float
-    #hvbar = zeros( (len(hc), 3), Float)
-    #for i in range(3):
-    #    hvbar[:,i] = hc[:]
-
-    from shallow_water_ext import balance_deep_and_shallow
-    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
-                             xmomc, ymomc, xmomv, ymomv)
-
-
-
-
-###############################################
-#Boundaries - specific to the shallow water wave equation
-class Reflective_boundary(Boundary):
-    """Reflective boundary returns same conserved quantities as
-    those present in its neighbour volume but reflected.
-
-    This class is specific to the shallow water equation as it
-    works with the momentum quantities assumed to be the second
-    and third conserved quantities.
-    """
-
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for reflective boundary'
-            raise msg
-
-        #Handy shorthands
-        self.stage   = domain.quantities['stage'].edge_values
-        self.xmom    = domain.quantities['xmomentum'].edge_values
-        self.ymom    = domain.quantities['ymomentum'].edge_values
-        self.normals = domain.normals
-
-        from Numeric import zeros, Float
-        self.conserved_quantities = zeros(3, Float)
-
-    def __repr__(self):
-        return 'Reflective_boundary'
-
-
-    def evaluate(self, vol_id, edge_id):
-        """Reflective boundaries reverses the outward momentum
-        of the volume they serve.
-        """
-
-        q = self.conserved_quantities
-        q[0] = self.stage[vol_id, edge_id]
-        q[1] = self.xmom[vol_id, edge_id]
-        q[2] = self.ymom[vol_id, edge_id]
-
-        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
-
-
-        r = rotate(q, normal, direction = 1)
-        r[1] = -r[1]
-        q = rotate(r, normal, direction = -1)
-
-        return q
-
-
-
-class Transmissive_Momentum_Set_Stage_boundary(Boundary):
-    """Returns same momentum conserved quantities as
-    those present in its neighbour volume. Sets stage
-
-    Underlying domain must be specified when boundary is instantiated
-    """
-
-    def __init__(self, domain = None, function=None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for this type boundary'
-            raise msg
-
-        if function is None:
-            msg = 'Function must be specified for this type boundary'
-            raise msg
-
-        self.domain   = domain
-        self.function = function
-
-    def __repr__(self):
-        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
-
-    def evaluate(self, vol_id, edge_id):
-        """Transmissive Momentum Set Stage boundaries return the edge momentum
-        values of the volume they serve.
-        """
-
-        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
-        value = self.function(self.domain.time)
-        q[0] = value[0]
-        return q
-
-
-        #FIXME: Consider this (taken from File_boundary) to allow
-        #spatial variation
-        #if vol_id is not None and edge_id is not None:
-        #    i = self.boundary_indices[ vol_id, edge_id ]
-        #    return self.F(t, point_id = i)
-        #else:
-        #    return self.F(t)
-
-
-
-class Dirichlet_Discharge_boundary(Boundary):
-    """
-    Sets stage (stage0)
-    Sets momentum (wh0) in the inward normal direction.
-
-    Underlying domain must be specified when boundary is instantiated
-    """
-
-    def __init__(self, domain = None, stage0=None, wh0=None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for this type boundary'
-            raise msg
-
-        if stage0 is None:
-            raise 'set stage'
-
-        if wh0 is None:
-            wh0 = 0.0
-
-        self.domain   = domain
-        self.stage0  = stage0
-        self.wh0 = wh0
-
-    def __repr__(self):
-        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
-
-    def evaluate(self, vol_id, edge_id):
-        """Set discharge in the (inward) normal direction
-        """
-
-        normal = self.domain.get_normal(vol_id,edge_id)
-        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
-        return q
-
-
-        #FIXME: Consider this (taken from File_boundary) to allow
-        #spatial variation
-        #if vol_id is not None and edge_id is not None:
-        #    i = self.boundary_indices[ vol_id, edge_id ]
-        #    return self.F(t, point_id = i)
-        #else:
-        #    return self.F(t)
-
-
-
-#class Spatio_temporal_boundary(Boundary):
-#    """The spatio-temporal boundary, reads values for the conserved
-#    quantities from an sww NetCDF file, and returns interpolated values
-#    at the midpoints of each associated boundaty segment.
-#    Time dependency is interpolated linearly as in util.File_function.#
-#
-#    Example:
-#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
-#
-#    """
-Spatio_temporal_boundary = File_boundary
-
-
-
-
-#########################
-#Standard forcing terms:
-#
-def gravity(domain):
-    """Apply gravitational pull in the presence of bed slope
-    """
-
-    from anuga.pyvolution.util import gradient
-    from Numeric import zeros, Float, array, sum
-
-    xmom = domain.quantities['xmomentum'].explicit_update
-    ymom = domain.quantities['ymomentum'].explicit_update
-
-    Stage = domain.quantities['stage']
-    Elevation = domain.quantities['elevation']
-    h = Stage.edge_values - Elevation.edge_values
-    v = Elevation.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-
-    for k in range(domain.number_of_elements):
-        avg_h = sum( h[k,:] )/3
-
-        #Compute bed slope
-        x0, y0, x1, y1, x2, y2 = x[k,:]
-        z0, z1, z2 = v[k,:]
-
-        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
-
-        #Update momentum
-        xmom[k] += -g*zx*avg_h
-        ymom[k] += -g*zy*avg_h
-
-
-def gravity_c(domain):
-    """Wrapper calling C version
-    """
-
-    xmom = domain.quantities['xmomentum'].explicit_update
-    ymom = domain.quantities['ymomentum'].explicit_update
-
-    Stage = domain.quantities['stage']
-    Elevation = domain.quantities['elevation']
-    h = Stage.edge_values - Elevation.edge_values
-    v = Elevation.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-
-
-    from shallow_water_ext import gravity
-    gravity(g, h, v, x, xmom, ymom)
-
-
-def manning_friction(domain):
-    """Apply (Manning) friction to water momentum
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-    vh = domain.quantities['ymomentum'].centroid_values
-    eta = domain.quantities['friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g
-
-    for k in range(N):
-        if eta[k] >= eps:
-            if h[k] >= eps:
-            	S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
-            	S /= h[k]**(7.0/3)
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
-            	ymom_update[k] += S*vh[k]
-
-
-def manning_friction_c(domain):
-    """Wrapper for c version
-    """
-
-
-    xmom = domain.quantities['xmomentum']
-    ymom = domain.quantities['ymomentum']
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-
-    uh = xmom.centroid_values
-    vh = ymom.centroid_values
-    eta = domain.quantities['friction'].centroid_values
-
-    xmom_update = xmom.semi_implicit_update
-    ymom_update = ymom.semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g
-
-    from shallow_water_ext import manning_friction
-    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
-
-
-def linear_friction(domain):
-    """Apply linear friction to water momentum
-
-    Assumes quantity: 'linear_friction' to be present
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-    vh = domain.quantities['ymomentum'].centroid_values
-    tau = domain.quantities['linear_friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g #Not necessary? Why was this added?
-
-    for k in range(N):
-        if tau[k] >= eps:
-            if h[k] >= eps:
-            	S = -tau[k]/h[k]
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
-            	ymom_update[k] += S*vh[k]
-
-
-
-def check_forcefield(f):
-    """Check that f is either
-    1: a callable object f(t,x,y), where x and y are vectors
-       and that it returns an array or a list of same length
-       as x and y
-    2: a scalar
-    """
-
-    from Numeric import ones, Float, array
-
-
-    if callable(f):
-        N = 3
-        x = ones(3, Float)
-        y = ones(3, Float)
-        try:
-            q = f(1.0, x=x, y=y)
-        except Exception, e:
-            msg = 'Function %s could not be executed:\n%s' %(f, e)
-	    #FIXME: Reconsider this semantics
-            raise msg
-
-        try:
-            q = array(q).astype(Float)
-        except:
-            msg = 'Return value from vector function %s could ' %f
-            msg += 'not be converted into a Numeric array of floats.\n'
-            msg += 'Specified function should return either list or array.'
-            raise msg
-
-        #Is this really what we want?
-        msg = 'Return vector from function %s ' %f
-        msg += 'must have same lenght as input vectors'
-        assert len(q) == N, msg
-
-    else:
-        try:
-            f = float(f)
-        except:
-            msg = 'Force field %s must be either a scalar' %f
-            msg += ' or a vector function'
-            raise msg
-    return f
-
-
-class Wind_stress:
-    """Apply wind stress to water momentum in terms of
-    wind speed [m/s] and wind direction [degrees]
-    """
-
-    def __init__(self, *args, **kwargs):
-        """Initialise windfield from wind speed s [m/s]
-        and wind direction phi [degrees]
-
-        Inputs v and phi can be either scalars or Python functions, e.g.
-
-        W = Wind_stress(10, 178)
-
-        #FIXME - 'normal' degrees are assumed for now, i.e. the
-        vector (1,0) has zero degrees.
-        We may need to convert from 'compass' degrees later on and also
-        map from True north to grid north.
-
-        Arguments can also be Python functions of t,x,y as in
-
-        def speed(t,x,y):
-            ...
-            return s
-
-        def angle(t,x,y):
-            ...
-            return phi
-
-        where x and y are vectors.
-
-        and then pass the functions in
-
-        W = Wind_stress(speed, angle)
-
-        The instantiated object W can be appended to the list of
-        forcing_terms as in
-
-        Alternatively, one vector valued function for (speed, angle)
-        can be applied, providing both quantities simultaneously.
-        As in
-        W = Wind_stress(F), where returns (speed, angle) for each t.
-
-        domain.forcing_terms.append(W)
-        """
-
-        from anuga.config import rho_a, rho_w, eta_w
-        from Numeric import array, Float
-
-        if len(args) == 2:
-            s = args[0]
-            phi = args[1]
-        elif len(args) == 1:
-            #Assume vector function returning (s, phi)(t,x,y)
-            vector_function = args[0]
-            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
-            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
-        else:
-           #Assume info is in 2 keyword arguments
-
-           if len(kwargs) == 2:
-               s = kwargs['s']
-               phi = kwargs['phi']
-           else:
-               raise 'Assumes two keyword arguments: s=..., phi=....'
-
-        self.speed = check_forcefield(s)
-        self.phi = check_forcefield(phi)
-
-        self.const = eta_w*rho_a/rho_w
-
-
-    def __call__(self, domain):
-        """Evaluate windfield based on values found in domain
-        """
-
-        from math import pi, cos, sin, sqrt
-        from Numeric import Float, ones, ArrayType
-
-        xmom_update = domain.quantities['xmomentum'].explicit_update
-        ymom_update = domain.quantities['ymomentum'].explicit_update
-
-        N = domain.number_of_elements
-        t = domain.time
-
-        if callable(self.speed):
-            xc = domain.get_centroid_coordinates()
-            s_vec = self.speed(t, xc[:,0], xc[:,1])
-        else:
-            #Assume s is a scalar
-
-            try:
-                s_vec = self.speed * ones(N, Float)
-            except:
-                msg = 'Speed must be either callable or a scalar: %s' %self.s
-                raise msg
-
-
-        if callable(self.phi):
-            xc = domain.get_centroid_coordinates()
-            phi_vec = self.phi(t, xc[:,0], xc[:,1])
-        else:
-            #Assume phi is a scalar
-
-            try:
-                phi_vec = self.phi * ones(N, Float)
-            except:
-                msg = 'Angle must be either callable or a scalar: %s' %self.phi
-                raise msg
-
-        assign_windfield_values(xmom_update, ymom_update,
-                                s_vec, phi_vec, self.const)
-
-
-def assign_windfield_values(xmom_update, ymom_update,
-                            s_vec, phi_vec, const):
-    """Python version of assigning wind field to update vectors.
-    A c version also exists (for speed)
-    """
-    from math import pi, cos, sin, sqrt
-
-    N = len(s_vec)
-    for k in range(N):
-        s = s_vec[k]
-        phi = phi_vec[k]
-
-        #Convert to radians
-        phi = phi*pi/180
-
-        #Compute velocity vector (u, v)
-        u = s*cos(phi)
-        v = s*sin(phi)
-
-        #Compute wind stress
-        S = const * sqrt(u**2 + v**2)
-        xmom_update[k] += S*u
-        ymom_update[k] += S*v
-
-
-
-##############################
-#OBSOLETE STUFF
-
-def balance_deep_and_shallow_old(domain):
-    """Compute linear combination between stage as computed by
-    gradient-limiters and stage computed as constant height above bed.
-    The former takes precedence when heights are large compared to the
-    bed slope while the latter takes precedence when heights are
-    relatively small.  Anything in between is computed as a balanced
-    linear combination in order to avoid numerical disturbances which
-    would otherwise appear as a result of hard switching between
-    modes.
-    """
-
-    #OBSOLETE
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-
-    #Computed linear combination between constant stages and and
-    #stages parallel to the bed elevation.
-    for k in range(domain.number_of_elements):
-        #Compute maximal variation in bed elevation
-        #  This quantitiy is
-        #    dz = max_i abs(z_i - z_c)
-        #  and it is independent of dimension
-        #  In the 1d case zc = (z0+z1)/2
-        #  In the 2d case zc = (z0+z1+z2)/3
-
-        dz = max(abs(zv[k,0]-zc[k]),
-                 abs(zv[k,1]-zc[k]),
-                 abs(zv[k,2]-zc[k]))
-
-
-        hmin = min( hv[k,:] )
-
-        #Create alpha in [0,1], where alpha==0 means using shallow
-        #first order scheme and alpha==1 means using the stage w as
-        #computed by the gradient limiter (1st or 2nd order)
-        #
-        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
-        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-
-        if dz > 0.0:
-            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
-        else:
-            #Flat bed
-            alpha = 1.0
-
-        #Weighted balance between stage parallel to bed elevation
-        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
-	#order gradient limiter
-        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
-        #
-        #It follows that the updated wvi is
-        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
-        #  zvi + hc + alpha*(hvi - hc)
-        #
-        #Note that hvi = zc+hc-zvi in the first order case (constant).
-
-	if alpha < 1:
-            for i in range(3):
-                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-
-
-            #Momentums at centroids
-            xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values
-
-            #Momentums at vertices
-            xmomv = domain.quantities['xmomentum'].vertex_values
-            ymomv = domain.quantities['ymomentum'].vertex_values
-
-            # Update momentum as a linear combination of
-            # xmomc and ymomc (shallow) and momentum
-            # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
-            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-
-
-
-
-
-###########################
-###########################
-#Geometries
-
-
-#FIXME: Rethink this way of creating values.
-
-
-class Weir:
-    """Set a bathymetry for weir with a hole and a downstream gutter
-    x,y are assumed to be in the unit square
-    """
-
-    def __init__(self, stage):
-        self.inflow_stage = stage
-
-    def __call__(self, x, y):
-        from Numeric import zeros, Float
-        from math import sqrt
-
-        N = len(x)
-        assert N == len(y)
-
-        z = zeros(N, Float)
-        for i in range(N):
-            z[i] = -x[i]/2  #General slope
-
-            #Flattish bit to the left
-            if x[i] < 0.3:
-                z[i] = -x[i]/10
-
-            #Weir
-            if x[i] >= 0.3 and x[i] < 0.4:
-                z[i] = -x[i]+0.9
-
-            #Dip
-            x0 = 0.6
-            #depth = -1.3
-            depth = -1.0
-            #plateaux = -0.9
-            plateaux = -0.6
-            if y[i] < 0.7:
-                if x[i] > x0 and x[i] < 0.9:
-                    z[i] = depth
-
-                #RHS plateaux
-                if x[i] >= 0.9:
-                    z[i] = plateaux
-
-
-            elif y[i] >= 0.7 and y[i] < 1.5:
-                #Restrict and deepen
-                if x[i] >= x0 and x[i] < 0.8:
-                    z[i] = depth-(y[i]/3-0.3)
-                    #z[i] = depth-y[i]/5
-                    #z[i] = depth
-                elif x[i] >= 0.8:
-                    #RHS plateaux
-                    z[i] = plateaux
-
-            elif y[i] >= 1.5:
-                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
-                    #Widen up and stay at constant depth
-                    z[i] = depth-1.5/5
-                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
-                    #RHS plateaux
-                    z[i] = plateaux
-
-
-            #Hole in weir (slightly higher than inflow condition)
-            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i]+self.inflow_stage + 0.02
-
-            #Channel behind weir
-            x0 = 0.5
-            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i]+self.inflow_stage + 0.02
-
-            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
-                #Flatten it out towards the end
-                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
-
-            #Hole to the east
-            x0 = 1.1; y0 = 0.35
-            #if x[i] < -0.2 and y < 0.5:
-            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
-                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
-
-            #Tiny channel draining hole
-            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
-                z[i] = -0.9 #North south
-
-            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
-                z[i] = -1.0 + (x[i]-0.9)/3 #East west
-
-
-
-            #Stuff not in use
-
-            #Upward slope at inlet to the north west
-            #if x[i] < 0.0: # and y[i] > 0.5:
-            #    #z[i] = -y[i]+0.5  #-x[i]/2
-            #    z[i] = x[i]/4 - y[i]**2 + 0.5
-
-            #Hole to the west
-            #x0 = -0.4; y0 = 0.35 # center
-            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
-            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
-
-
-
-
-
-        return z/2
-
-class Weir_simple:
-    """Set a bathymetry for weir with a hole and a downstream gutter
-    x,y are assumed to be in the unit square
-    """
-
-    def __init__(self, stage):
-        self.inflow_stage = stage
-
-    def __call__(self, x, y):
-        from Numeric import zeros, Float
-
-        N = len(x)
-        assert N == len(y)
-
-        z = zeros(N, Float)
-        for i in range(N):
-            z[i] = -x[i]  #General slope
-
-            #Flat bit to the left
-            if x[i] < 0.3:
-                z[i] = -x[i]/10  #General slope
-
-            #Weir
-            if x[i] > 0.3 and x[i] < 0.4:
-                z[i] = -x[i]+0.9
-
-            #Dip
-            if x[i] > 0.6 and x[i] < 0.9:
-                z[i] = -x[i]-0.5  #-y[i]/5
-
-            #Hole in weir (slightly higher than inflow condition)
-            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i]+self.inflow_stage + 0.05
-
-
-        return z/2
-
-
-
-class Constant_stage:
-    """Set an initial condition with constant stage
-    """
-    def __init__(self, s):
-        self.s = s
-
-    def __call__(self, x, y):
-        return self.s
-
-class Constant_height:
-    """Set an initial condition with constant water height, e.g
-    stage s = z+h
-    """
-
-    def __init__(self, W, h):
-        self.W = W
-        self.h = h
-
-    def __call__(self, x, y):
-        if self.W is None:
-            from Numeric import ones, Float
-            return self.h*ones(len(x), Float)
-        else:
-            return self.W(x,y) + self.h
-
-
-
-
-
-
-
-##############################################
-#Initialise module
-
-
-from anuga.utilities import compile
-if compile.can_use_C_extension('shallow_water_kinetic_ext.c'):
-    #Replace python version with c implementations
-
-    from shallow_water_kinetic_ext import rotate, assign_windfield_values
-    compute_fluxes = compute_fluxes_c
-    extrapolate_second_order_sw=extrapolate_second_order_sw_c
-    gravity = gravity_c
-    manning_friction = manning_friction_c
-    h_limiter = h_limiter_c
-    balance_deep_and_shallow = balance_deep_and_shallow_c
-    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
-
-
-    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c #(like MH's)
-
-
-
-#Optimisation with psyco
-from anuga.config import use_psyco
-if use_psyco:
-    try:
-        import psyco
-    except:
-        import os
-        if os.name == 'posix' and os.uname()[4] == 'x86_64':
-            pass
-            #Psyco isn't supported on 64 bit systems, but it doesn't matter
-        else:
-            msg = 'WARNING: psyco (speedup) could not import'+\
-                  ', you may want to consider installing it'
-            print msg
-    else:
-        psyco.bind(Domain.distribute_to_vertices_and_edges)
-        psyco.bind(Domain.compute_fluxes)
-
-if __name__ == "__main__":
-    pass
-
-# Profiling stuff
-import profile
-profiler = profile.Profile()
Index: uga_core/source/anuga/abstract_2d_finite_volumes/shallow_water_kinetic_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/shallow_water_kinetic_ext.c	(revision 3697)
+++ 	(revision )
@@ -1,1442 +1,0 @@
-// Python - C extension module for shallow_water.py
-//
-// To compile (Python2.3):
-//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
-//  gcc -shared domain_ext.o  -o domain_ext.so
-//
-// or use python compile.py
-//
-// See the module shallow_water.py
-//
-//
-// Ole Nielsen, GA 2004
-
-
-#include "Python.h"
-#include "Numeric/arrayobject.h"
-#include "math.h"
-#include <stdio.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 momentum component q (q[1], q[2])
-    from x,y coordinates to coordinates based on normal vector (n1, n2).
-
-    Result is returned in array 3x1 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];  //uh momentum
-  q2 = q[2];  //vh momentum
-
-  //Rotate
-  q[1] =  n1*q1 + n2*q2;
-  q[2] = -n2*q1 + n1*q2;
-
-  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)
-  if (dq0>=0.0){
-    if (dq1>=dq2){
-      if (dq1>=0.0)
-	*qmax=dq0+dq1;
-      else
-	*qmax=dq0;
-      if ((*qmin=dq0+dq2)<0)
-	;//qmin is already set to correct value
-      else
-	*qmin=0.0;
-    }
-    else{//dq1<dq2
-      if (dq2>0)
-	*qmax=dq0+dq2;
-      else
-	*qmax=dq0;
-      if ((*qmin=dq0+dq1)<0)
-	;//qmin is the correct value
-      else
-	*qmin=0.0;
-    }
-  }
-  else{//dq0<0
-    if (dq1<=dq2){
-      if (dq1<0.0)
-	*qmin=dq0+dq1;
-      else
-	*qmin=dq0;
-      if ((*qmax=dq0+dq2)>0.0)
-	;//qmax is already set to the correct value
-      else
-	*qmax=0.0;
-    }
-    else{//dq1>dq2
-      if (dq2<0.0)
-	*qmin=dq0+dq2;
-      else
-	*qmin=dq0;
-      if ((*qmax=dq0+dq1)>0.0)
-	;//qmax is already set to the correct value
-      else
-	*qmax=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.
-  //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);
-  for (i=0;i<3;i++)
-    dqv[i]=dqv[i]*phi;
-  return 0;
-}
-
-// Computational function for flux computation (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 g,
-		  double *edgeflux, double *max_speed) {
-
-  /*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, z;
-  double q_left_copy[3], q_right_copy[3];
-  double flux_right[3], flux_left[3];
-
-  //Copy conserved quantities to protect from modification
-  for (i=0; i<3; i++) {
-    q_left_copy[i] = q_left[i];
-    q_right_copy[i] = q_right[i];
-  }
-
-  //Align x- and y-momentum with x-axis
-  _rotate(q_left_copy, n1, n2);
-  _rotate(q_right_copy, n1, n2);
-
-  z = (z_left+z_right)/2; //Take average of field values
-
-  //Compute speeds in x-direction
-  w_left = q_left_copy[0];              // h+z
-  h_left = w_left-z;
-  uh_left = q_left_copy[1];
-
-  if (h_left < epsilon) {
-    h_left = 0.0;  //Could have been negative
-    u_left = 0.0;
-  } else {
-    u_left = uh_left/h_left;
-  }
-
-  w_right = q_right_copy[0];
-  h_right = w_right-z;
-  uh_right = q_right_copy[1];
-
-  if (h_right < epsilon) {
-    h_right = 0.0; //Could have been negative
-    u_right = 0.0;
-  } else {
-    u_right = uh_right/h_right;
-  }
-
-  //Momentum in y-direction
-  vh_left  = q_left_copy[2];
-  vh_right = q_right_copy[2];
-
-
-  //Maximal and minimal wave speeds
-  soundspeed_left  = sqrt(g*h_left);
-  soundspeed_right = sqrt(g*h_right);
-
-  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
-  if (s_max < 0.0) s_max = 0.0;
-
-  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
-  if (s_min > 0.0) s_min = 0.0;
-
-  //Flux formulas
-  flux_left[0] = u_left*h_left;
-  flux_left[1] = u_left*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 == 0.0) {
-    for (i=0; i<3; i++) edgeflux[i] = 0.0;
-    *max_speed = 0.0;
-  } else {
-    for (i=0; i<3; i++) {
-      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
-      edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
-      edgeflux[i] /= denom;
-    }
-
-    //Maximal wavespeed
-    *max_speed = max(fabs(s_max), fabs(s_min));
-
-    //Rotate back
-    _rotate(edgeflux, n1, -n2);
-  }
-  return 0;
-}
-
-double erfcc(double x){
-    double z,t,result;
-
-    z=fabs(x);
-    t=1.0/(1.0+0.5*z);
-    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
-         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
-         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
-    if (x < 0.0) result = 2.0-result;
-
-    return result;
-    }
-
-
-
-// Computational function for flux computation (using stage w=z+h)
-int flux_function_kinetic(double *q_left, double *q_right,
-		  double z_left, double z_right,
-		  double n1, double n2,
-		  double epsilon, double g,
-		  double *edgeflux, double *max_speed) {
-
-  /*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
-
-    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
-  */
-
-  int i;
-
-  double w_left, h_left, uh_left, vh_left, u_left, F_left;
-  double w_right, h_right, uh_right, vh_right, u_right, F_right;
-  double s_min, s_max, soundspeed_left, soundspeed_right;
-  double z;
-  double q_left_copy[3], q_right_copy[3];
-
-
-  //Copy conserved quantities to protect from modification
-  for (i=0; i<3; i++) {
-    q_left_copy[i] = q_left[i];
-    q_right_copy[i] = q_right[i];
-  }
-
-  //Align x- and y-momentum with x-axis
-  _rotate(q_left_copy, n1, n2);
-  _rotate(q_right_copy, n1, n2);
-
-  z = (z_left+z_right)/2; //Take average of field values
-
-  //Compute speeds in x-direction
-  w_left = q_left_copy[0];              // h+z
-  h_left = w_left-z;
-  uh_left = q_left_copy[1];
-
-  if (h_left < epsilon) {
-    h_left = 0.0;  //Could have been negative
-    u_left = 0.0;
-  } else {
-    u_left = uh_left/h_left;
-  }
-
-  w_right = q_right_copy[0];
-  h_right = w_right-z;
-  uh_right = q_right_copy[1];
-
-  if (h_right < epsilon) {
-    h_right = 0.0; //Could have been negative
-    u_right = 0.0;
-  } else {
-    u_right = uh_right/h_right;
-  }
-
-  //Momentum in y-direction
-  vh_left  = q_left_copy[2];
-  vh_right = q_right_copy[2];
-
-
-  //Maximal and minimal wave speeds
-  soundspeed_left  = sqrt(g*h_left);
-  soundspeed_right = sqrt(g*h_right);
-
-  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
-  if (s_max < 0.0) s_max = 0.0;
-
-  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
-  if (s_min > 0.0) s_min = 0.0;
-
-
-  F_left  = 0.0;
-  F_right = 0.0;
-  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
-  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
-
-  for (i=0; i<3; i++) edgeflux[i] = 0.0;
-  *max_speed = 0.0;
-
-  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
-          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
-          h_right*u_right/2.0*erfcc(F_right) -  \
-          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
-
-  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
-          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
-          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
-          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
-
-  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
-          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
-          vh_right*u_right/2.0*erfcc(F_right) - \
-          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
-
-  //Maximal wavespeed
-  *max_speed = max(fabs(s_max), fabs(s_min));
-
-  //Rotate back
-  _rotate(edgeflux, n1, -n2);
-
-  return 0;
-}
-
-
-void _manning_friction(double g, double eps, int N,
-		       double* w, double* z,
-		       double* uh, double* vh,
-		       double* eta, double* xmom, double* ymom) {
-
-  int k;
-  double S, h;
-
-  for (k=0; k<N; k++) {
-    if (eta[k] > eps) {
-      h = w[k]-z[k];
-      if (h >= eps) {
-	S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-	//S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
-	S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
-	//S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
-
-
-	//Update momentum
-	xmom[k] += S*uh[k];
-	ymom[k] += S*vh[k];
-      }
-    }
-  }
-}
-
-
-
-int _balance_deep_and_shallow(int N,
-			      double* wc,
-			      double* zc,
-			      double* hc,
-			      double* wv,
-			      double* zv,
-			      double* hv,
-			      double* hvbar,
-			      double* xmomc,
-			      double* ymomc,
-			      double* xmomv,
-			      double* ymomv) {
-
-  int k, k3, i;
-  double dz, hmin, alpha;
-
-  //Compute linear combination between w-limited stages and
-  //h-limited stages close to the bed elevation.
-
-  for (k=0; k<N; k++) {
-    // Compute maximal variation in bed elevation
-    // This quantitiy is
-    //     dz = max_i abs(z_i - z_c)
-    // and it is independent of dimension
-    // In the 1d case zc = (z0+z1)/2
-    // In the 2d case zc = (z0+z1+z2)/3
-
-    k3 = 3*k;
-
-    //FIXME: Try with this one precomputed
-    dz = 0.0;
-    hmin = hv[k3];
-    for (i=0; i<3; i++) {
-      dz = max(dz, fabs(zv[k3+i]-zc[k]));
-      hmin = min(hmin, hv[k3+i]);
-    }
-
-
-    //Create alpha in [0,1], where alpha==0 means using the h-limited
-    //stage and alpha==1 means using the w-limited stage as
-    //computed by the gradient limiter (both 1st or 2nd order)
-    //
-    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
-    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
-
-
-    if (dz > 0.0)
-      //if (hmin<0.0)
-      //	alpha = 0.0;
-      //else
-      //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
-      alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 );
-    else
-      alpha = 1.0;  //Flat bed
-
-    //alpha = 1.0;
-
-    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
-
-    //	Let
-    //
-    //	  wvi be the w-limited stage (wvi = zvi + hvi)
-    //	  wvi- be the h-limited state (wvi- = zvi + hvi-)
-    //
-    //
-    //	where i=0,1,2 denotes the vertex ids
-    //
-    //  Weighted balance between w-limited and h-limited stage is
-    //
-    //	  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
-    //
-    //  It follows that the updated wvi is
-    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
-    //
-    //	 Momentum is balanced between constant and limited
-
-    if (alpha < 1) {
-      for (i=0; i<3; i++) {
-         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
-
-	//Update momentum as a linear combination of
-	//xmomc and ymomc (shallow) and momentum
-	//from extrapolator xmomv and ymomv (deep).
-	xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
-	ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-      }
-    }
-  }
-  return 0;
-}
-
-
-
-int _protect(int N,
-	     double minimum_allowed_height,
-	     double epsilon,
-	     double* wc,
-	     double* zc,
-	     double* xmomc,
-	     double* ymomc) {
-
-  int k;
-  double hc;
-
-  //Protect against initesimal and negative heights
-
-  for (k=0; k<N; k++) {
-    hc = wc[k] - zc[k];
-
-    if (hc < minimum_allowed_height) {
-      if (hc < epsilon) wc[k] = zc[k]; //Contain 'lost mass' error
-      xmomc[k] = 0.0;
-      ymomc[k] = 0.0;
-    }
-
-  }
-  return 0;
-}
-
-
-
-
-int _assign_wind_field_values(int N,
-			      double* xmom_update,
-			      double* ymom_update,
-			      double* s_vec,
-			      double* phi_vec,
-			      double cw) {
-
-  //Assign windfield values to momentum updates
-
-  int k;
-  double S, s, phi, u, v;
-
-  for (k=0; k<N; k++) {
-
-    s = s_vec[k];
-    phi = phi_vec[k];
-
-    //Convert to radians
-    phi = phi*pi/180;
-
-    //Compute velocity vector (u, v)
-    u = s*cos(phi);
-    v = s*sin(phi);
-
-    //Compute wind stress
-    S = cw * sqrt(u*u + v*v);
-    xmom_update[k] += S*u;
-    ymom_update[k] += S*v;
-  }
-  return 0;
-}
-
-
-
-///////////////////////////////////////////////////////////////////
-// Gateways to Python
-
-PyObject *gravity(PyObject *self, PyObject *args) {
-  //
-  //  gravity(g, h, v, x, xmom, ymom)
-  //
-
-
-  PyArrayObject *h, *v, *x, *xmom, *ymom;
-  int k, i, N, k3, k6;
-  double g, avg_h, zx, zy;
-  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
-
-  if (!PyArg_ParseTuple(args, "dOOOOO",
-			&g, &h, &v, &x,
-			&xmom, &ymom))
-    return NULL;
-
-  N = h -> dimensions[0];
-  for (k=0; k<N; k++) {
-    k3 = 3*k;  // base index
-    k6 = 6*k;  // base index
-
-    avg_h = 0.0;
-    for (i=0; i<3; i++) {
-      avg_h += ((double *) h -> data)[k3+i];
-    }
-    avg_h /= 3;
-
-
-    //Compute bed slope
-    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];
-
-
-    z0 = ((double*) v -> data)[k3 + 0];
-    z1 = ((double*) v -> data)[k3 + 1];
-    z2 = ((double*) v -> data)[k3 + 2];
-
-    _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 *manning_friction(PyObject *self, PyObject *args) {
-  //
-  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
-  //
-
-
-  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
-  int N;
-  double g, eps;
-
-  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-			&g, &eps, &w, &z, &uh, &vh, &eta,
-			&xmom, &ymom))
-    return NULL;
-
-  N = w -> dimensions[0];
-  _manning_friction(g, eps, N,
-		    (double*) w -> data,
-		    (double*) z -> data,
-		    (double*) uh -> data,
-		    (double*) vh -> data,
-		    (double*) eta -> data,
-		    (double*) xmom -> data,
-		    (double*) ymom -> data);
-
-  return Py_BuildValue("");
-}
-
-PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
-  /*Compute the vertex values based on a linear reconstruction on each triangle
-    These values are calculated as follows:
-    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
-    formed by the centroids of its three neighbours.
-    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
-    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
-    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
-    original triangle has gradient (a,b).
-    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
-    find_qmin_and_qmax and limit_gradient
-
-    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.vertex_coordinates,
-                                Stage.vertex_values,
-                                Xmom.vertex_values,
-                                Ymom.vertex_values)
-
-    Post conditions:
-            The vertices 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,
-    *vertex_coordinates,
-    *stage_vertex_values,
-    *xmom_vertex_values,
-    *ymom_vertex_values;
-  PyObject *domain, *Tmp;
-  //Local variables
-  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
-  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
-  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;
-  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
-  //by which these jumps are limited
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
-			&domain,
-			&surrogate_neighbours,
-			&number_of_boundaries,
-			&centroid_coordinates,
-			&stage_centroid_values,
-			&xmom_centroid_values,
-			&ymom_centroid_values,
-			&vertex_coordinates,
-			&stage_vertex_values,
-			&xmom_vertex_values,
-			&ymom_vertex_values)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
-  }
-
-  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
-  Tmp = PyObject_GetAttrString(domain, "beta_w");
-  if (!Tmp)
-    return NULL;
-  beta_w = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  number_of_elements = stage_centroid_values -> dimensions[0];
-  for (k=0; k<number_of_elements; k++) {
-    k3=k*3;
-    k6=k*6;
-
-    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
-      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
-      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
-      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
-      continue;
-    }
-    else {//we will need centroid coordinates and vertex coordinates of the triangle
-      //get the vertex coordinates of the FV triangle
-      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
-      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
-      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
-      //get the centroid coordinates of the FV triangle
-      coord_index=2*k;
-      x=((double *)centroid_coordinates->data)[coord_index];
-      y=((double *)centroid_coordinates->data)[coord_index+1];
-      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
-      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
-      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
-    }
-    if (((long *)number_of_boundaries->data)[k]<=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=((long *)surrogate_neighbours->data)[k3];
-      k1=((long *)surrogate_neighbours->data)[k3+1];
-      k2=((long *)surrogate_neighbours->data)[k3+2];
-      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
-      coord_index=2*k0;
-      x0=((double *)centroid_coordinates->data)[coord_index];
-      y0=((double *)centroid_coordinates->data)[coord_index+1];
-      coord_index=2*k1;
-      x1=((double *)centroid_coordinates->data)[coord_index];
-      y1=((double *)centroid_coordinates->data)[coord_index+1];
-      coord_index=2*k2;
-      x2=((double *)centroid_coordinates->data)[coord_index];
-      y2=((double *)centroid_coordinates->data)[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
-      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
-      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
-      if (area2<=0)
-	return NULL;
-
-      //### stage ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
-      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
-      //calculate the gradient of stage on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri 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);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
-
-      //### xmom ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
-      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
-      //calculate the gradient of xmom on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri 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);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
-
-      //### ymom ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
-      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
-      //calculate the gradient of xmom on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri 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);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
-    }//if (number_of_boundaries[k]<=1)
-    else{//number_of_boundaries==2
-      //one internal neighbour and gradient is in direction of the neighbour's centroid
-      //find the only internal neighbour
-      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
-	if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
-	  break;
-      }
-      if ((k2==k3+3))//if we didn't find an internal neighbour
-	return NULL;//error
-      k1=((long *)surrogate_neighbours->data)[k2];
-      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
-      coord_index=2*k1;
-      x1=((double *)centroid_coordinates->data)[coord_index];
-      y1=((double *)centroid_coordinates->data)[coord_index+1];
-      //compute x- and y- distances between the centroid of the FV triangle 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=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
-
-      //## xmom ###
-      //compute differentials
-      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
-
-      //## ymom ###
-      //compute differentials
-      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
-    }//else [number_of_boudaries==2]
-  }//for k=0 to number_of_elements-1
-  return Py_BuildValue("");
-}//extrapolate_second-order_sw
-
-PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
-  //
-  // r = rotate(q, normal, direction=1)
-  //
-  // Where q is assumed to be a Float numeric array of length 3 and
-  // normal a Float numeric array of length 2.
-
-
-  PyObject *Q, *Normal;
-  PyArrayObject *q, *r, *normal;
-
-  static char *argnames[] = {"q", "normal", "direction", NULL};
-  int dimensions[1], i, direction=1;
-  double n1, n2;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
-				   &Q, &Normal, &direction))
-    return NULL;
-
-  //Input checks (convert sequences into numeric arrays)
-  q = (PyArrayObject *)
-    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
-  normal = (PyArrayObject *)
-    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
-
-
-  if (normal -> dimensions[0] != 2) {
-    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
-    return NULL;
-  }
-
-  //Allocate space for return vector r (don't DECREF)
-  dimensions[0] = 3;
-  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-  //Copy
-  for (i=0; i<3; i++) {
-    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
-  }
-
-  //Get normal and direction
-  n1 = ((double *) normal -> data)[0];
-  n2 = ((double *) normal -> data)[1];
-  if (direction == -1) n2 = -n2;
-
-  //Rotate
-  _rotate((double *) r -> data, n1, n2);
-
-  //Release numeric arrays
-  Py_DECREF(q);
-  Py_DECREF(normal);
-
-  //return result using PyArray to avoid memory leak
-  return PyArray_Return(r);
-}
-
-PyObject *compute_fluxes(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"
-
-    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.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     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)
-
-
-    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,
-    *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;//tracks whether the flux across an edge has already been computed
-
-
-  //Local variables
-  double timestep, max_speed, epsilon, g;
-  double normal[2], ql[3], qr[3], zl, zr;
-  double edgeflux[3]; //Work arrays for summing up fluxes
-
-  int number_of_elements, k, i, m, n;
-  int ki, nm=0, ki2; //Index shorthands
-  static long call=1;
-
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
-			&timestep,
-			&epsilon,
-			&g,
-			&neighbours,
-			&neighbour_edges,
-			&normals,
-			&edgelengths, &radii, &areas,
-			&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)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
-  }
-  number_of_elements = stage_edge_values -> dimensions[0];
-  call++;//a static local variable to which already_computed_flux is compared
-  //set explicit_update to zero for all conserved_quantities.
-  //This assumes compute_fluxes called before forcing terms
-  for (k=0; k<number_of_elements; k++) {
-    ((double *) stage_explicit_update -> data)[k]=0.0;
-    ((double *) xmom_explicit_update -> data)[k]=0.0;
-    ((double *) ymom_explicit_update -> data)[k]=0.0;
-  }
-  //Loop through neighbours and compute edge flux for each
-  for (k=0; k<number_of_elements; k++) {
-    for (i=0; i<3; i++) {
-      ki = k*3+i;
-      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
-	continue;
-      ql[0] = ((double *) stage_edge_values -> data)[ki];
-      ql[1] = ((double *) xmom_edge_values -> data)[ki];
-      ql[2] = ((double *) ymom_edge_values -> data)[ki];
-      zl =    ((double *) bed_edge_values -> data)[ki];
-
-      //Quantities at neighbour on nearest face
-      n = ((long *) neighbours -> data)[ki];
-      if (n < 0) {
-	m = -n-1; //Convert negative flag to index
-	qr[0] = ((double *) stage_boundary_values -> data)[m];
-	qr[1] = ((double *) xmom_boundary_values -> data)[m];
-	qr[2] = ((double *) ymom_boundary_values -> data)[m];
-	zr = zl; //Extend bed elevation to boundary
-      } else {
-	m = ((long *) neighbour_edges -> data)[ki];
-	nm = n*3+m;
-	qr[0] = ((double *) stage_edge_values -> data)[nm];
-	qr[1] = ((double *) xmom_edge_values -> data)[nm];
-	qr[2] = ((double *) ymom_edge_values -> data)[nm];
-	zr =    ((double *) bed_edge_values -> data)[nm];
-      }
-      // Outward pointing normal vector
-      // normal = domain.normals[k, 2*i:2*i+2]
-      ki2 = 2*ki; //k*6 + i*2
-      normal[0] = ((double *) normals -> data)[ki2];
-      normal[1] = ((double *) normals -> data)[ki2+1];
-      //Edge flux computation
-      flux_function_kinetic(ql, qr, zl, zr,
-		    normal[0], normal[1],
-		    epsilon, g,
-		    edgeflux, &max_speed);
-      //update triangle k
-      ((long *) already_computed_flux->data)[ki]=call;
-      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
-      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
-      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
-      //update the neighbour n
-      if (n>=0){
-	((long *) already_computed_flux->data)[nm]=call;
-	((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
-	((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
-	((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
-      }
-      ///for (j=0; j<3; j++) {
-	///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
-	///}
-	//Update timestep
-	//timestep = min(timestep, domain.radii[k]/max_speed)
-	//FIXME: SR Add parameter for CFL condition
-	if (max_speed > epsilon) {
-	  timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-	  //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
-	  if (n>=0)
-	    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
-	}
-    } // end for i
-    //Normalise by area and store for when all conserved
-    //quantities get updated
-    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
-    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
-    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
-  } //end for k
-  return Py_BuildValue("d", timestep);
-}
-
-PyObject *protect(PyObject *self, PyObject *args) {
-  //
-  //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
-
-
-  PyArrayObject
-  *wc,            //Stage at centroids
-  *zc,            //Elevation at centroids
-  *xmomc,         //Momentums at centroids
-  *ymomc;
-
-
-  int N;
-  double minimum_allowed_height, epsilon;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddOOOO",
-			&minimum_allowed_height,
-			&epsilon,
-			&wc, &zc, &xmomc, &ymomc))
-    return NULL;
-
-  N = wc -> dimensions[0];
-
-  _protect(N,
-	   minimum_allowed_height,
-	   epsilon,
-	   (double*) wc -> data,
-	   (double*) zc -> data,
-	   (double*) xmomc -> data,
-	   (double*) ymomc -> data);
-
-  return Py_BuildValue("");
-}
-
-
-
-PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
-  //
-  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
-  //                             xmomc, ymomc, xmomv, ymomv)
-
-
-  PyArrayObject
-    *wc,            //Stage at centroids
-    *zc,            //Elevation at centroids
-    *hc,            //Height at centroids
-    *wv,            //Stage at vertices
-    *zv,            //Elevation at vertices
-    *hv,            //Depths at vertices
-    *hvbar,         //h-Limited depths at vertices
-    *xmomc,         //Momentums at centroids and vertices
-    *ymomc,
-    *xmomv,
-    *ymomv;
-
-  int N; //, err;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
-			&wc, &zc, &hc,
-			&wv, &zv, &hv, &hvbar,
-			&xmomc, &ymomc, &xmomv, &ymomv))
-    return NULL;
-
-  N = wc -> dimensions[0];
-
-  _balance_deep_and_shallow(N,
-			    (double*) wc -> data,
-			    (double*) zc -> data,
-			    (double*) hc -> data,
-			    (double*) wv -> data,
-			    (double*) zv -> data,
-			    (double*) hv -> data,
-			    (double*) hvbar -> data,
-                            (double*) xmomc -> data,
-			    (double*) ymomc -> data,
-			    (double*) xmomv -> data,
-			    (double*) ymomv -> data);
-
-
-  return Py_BuildValue("");
-}
-
-
-
-PyObject *h_limiter(PyObject *self, PyObject *args) {
-
-  PyObject *domain, *Tmp;
-  PyArrayObject
-    *hv, *hc, //Depth at vertices and centroids
-    *hvbar,   //Limited depth at vertices (return values)
-    *neighbours;
-
-  int k, i, n, N, k3;
-  int dimensions[2];
-  double beta_h; //Safety factor (see config.py)
-  double *hmin, *hmax, hn;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
-    return NULL;
-
-  neighbours = get_consecutive_array(domain, "neighbours");
-
-  //Get safety factor beta_h
-  Tmp = PyObject_GetAttrString(domain, "beta_h");
-  if (!Tmp)
-    return NULL;
-
-  beta_h = PyFloat_AsDouble(Tmp);
-
-  Py_DECREF(Tmp);
-
-  N = hc -> dimensions[0];
-
-  //Create hvbar
-  dimensions[0] = N;
-  dimensions[1] = 3;
-  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
-
-
-  //Find min and max of this and neighbour's centroid values
-  hmin = malloc(N * sizeof(double));
-  hmax = malloc(N * sizeof(double));
-  for (k=0; k<N; k++) {
-    k3=k*3;
-
-    hmin[k] = ((double*) hc -> data)[k];
-    hmax[k] = hmin[k];
-
-    for (i=0; i<3; i++) {
-      n = ((long*) neighbours -> data)[k3+i];
-
-      //Initialise hvbar with values from hv
-      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
-
-      if (n >= 0) {
-	hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
-
-	hmin[k] = min(hmin[k], hn);
-	hmax[k] = max(hmax[k], hn);
-      }
-    }
-  }
-
-  // Call underlying standard routine
-  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
-
-  // // //Py_DECREF(domain); //FIXME: NEcessary?
-  free(hmin);
-  free(hmax);
-
-  //return result using PyArray to avoid memory leak
-  return PyArray_Return(hvbar);
-  //return Py_BuildValue("");
-}
-
-PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
-  //a faster version of h_limiter above
-  PyObject *domain, *Tmp;
-  PyArrayObject
-    *hv, *hc, //Depth at vertices and centroids
-    *hvbar,   //Limited depth at vertices (return values)
-    *neighbours;
-
-  int k, i, N, k3,k0,k1,k2;
-  int dimensions[2];
-  double beta_h; //Safety factor (see config.py)
-  double hmin, hmax, dh[3];
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
-    return NULL;
-  neighbours = get_consecutive_array(domain, "neighbours");
-
-  //Get safety factor beta_h
-  Tmp = PyObject_GetAttrString(domain, "beta_h");
-  if (!Tmp)
-    return NULL;
-  beta_h = PyFloat_AsDouble(Tmp);
-
-  Py_DECREF(Tmp);
-
-  N = hc -> dimensions[0];
-
-  //Create hvbar
-  dimensions[0] = N;
-  dimensions[1] = 3;
-  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
-  for (k=0;k<N;k++){
-    k3=k*3;
-    //get the ids of the neighbours
-    k0 = ((long*) neighbours -> data)[k3];
-    k1 = ((long*) neighbours -> data)[k3+1];
-    k2 = ((long*) neighbours -> data)[k3+2];
-    //set hvbar provisionally
-    for (i=0;i<3;i++){
-      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
-      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
-    }
-    hmin=((double*) hc -> data)[k];
-    hmax=hmin;
-    if (k0>=0){
-      hmin=min(hmin,((double*) hc -> data)[k0]);
-      hmax=max(hmax,((double*) hc -> data)[k0]);
-    }
-    if (k1>=0){
-      hmin=min(hmin,((double*) hc -> data)[k1]);
-      hmax=max(hmax,((double*) hc -> data)[k1]);
-    }
-    if (k2>=0){
-      hmin=min(hmin,((double*) hc -> data)[k2]);
-      hmax=max(hmax,((double*) hc -> data)[k2]);
-    }
-    hmin-=((double*) hc -> data)[k];
-    hmax-=((double*) hc -> data)[k];
-    limit_gradient(dh,hmin,hmax,beta_h);
-    for (i=0;i<3;i++)
-      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
-  }
-  return PyArray_Return(hvbar);
-}
-
-PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
-  //
-  //      assign_windfield_values(xmom_update, ymom_update,
-  //                              s_vec, phi_vec, self.const)
-
-
-
-  PyArrayObject   //(one element per triangle)
-  *s_vec,         //Speeds
-  *phi_vec,       //Bearings
-  *xmom_update,   //Momentum updates
-  *ymom_update;
-
-
-  int N;
-  double cw;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOd",
-			&xmom_update,
-			&ymom_update,
-			&s_vec, &phi_vec,
-			&cw))
-    return NULL;
-
-  N = xmom_update -> dimensions[0];
-
-  _assign_wind_field_values(N,
-	   (double*) xmom_update -> data,
-	   (double*) ymom_update -> data,
-	   (double*) s_vec -> data,
-	   (double*) phi_vec -> data,
-	   cw);
-
-  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"},
-  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
-  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
-  {"gravity", gravity, METH_VARARGS, "Print out"},
-  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
-  {"balance_deep_and_shallow", balance_deep_and_shallow,
-   METH_VARARGS, "Print out"},
-  {"h_limiter", h_limiter,
-   METH_VARARGS, "Print out"},
-  {"h_limiter_sw", h_limiter_sw,
-   METH_VARARGS, "Print out"},
-  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"assign_windfield_values", assign_windfield_values,
-   METH_VARARGS | METH_KEYWORDS, "Print out"},
-  //{"distribute_to_vertices_and_edges",
-  // distribute_to_vertices_and_edges, METH_VARARGS},
-  //{"update_conserved_quantities",
-  // update_conserved_quantities, METH_VARARGS},
-  //{"set_initialcondition",
-  // set_initialcondition, METH_VARARGS},
-  {NULL, NULL}
-};
-
-// Module initialisation
-void initshallow_water_kinetic_ext(void){
-  Py_InitModule("shallow_water_kinetic_ext", MethodTable);
-
-  import_array();     //Necessary for handling of NumPY structures
-}
Index: /anuga_core/source/anuga/shallow_water/shallow_water_kinetic.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_kinetic.py	(revision 3698)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_kinetic.py	(revision 3698)
@@ -0,0 +1,1815 @@
+"""Class Domain -
+2D triangular domains for finite-volume computations of
+the shallow water wave equation.
+
+This module contains a specialisation of class Domain from module domain.py
+consisting of methods specific to the Shallow Water Wave Equation
+
+
+U_t + E_x + G_y = S
+
+where
+
+U = [w, uh, vh]
+E = [uh, u^2h + gh^2/2, uvh]
+G = [vh, uvh, v^2h + gh^2/2]
+S represents source terms forcing the system
+(e.g. gravity, friction, wind stress, ...)
+
+and _t, _x, _y denote the derivative with respect to t, x and y respectively.
+
+The quantities are
+
+symbol    variable name    explanation
+x         x                horizontal distance from origin [m]
+y         y                vertical distance from origin [m]
+z         elevation        elevation of bed on which flow is modelled [m]
+h         height           water height above z [m]
+w         stage            absolute water level, w = z+h [m]
+u                          speed in the x direction [m/s]
+v                          speed in the y direction [m/s]
+uh        xmomentum        momentum in the x direction [m^2/s]
+vh        ymomentum        momentum in the y direction [m^2/s]
+
+eta                        mannings friction coefficient [to appear]
+nu                         wind stress coefficient [to appear]
+
+The conserved quantities are w, uh, vh
+
+
+For details see e.g.
+Christopher Zoppou and Stephen Roberts,
+Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
+Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
+
+
+Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2004
+"""
+
+#Subversion keywords:
+#
+#$LastChangedDate: 2006-01-13 16:43:01 +1100 (Fri, 13 Jan 2006) $
+#$LastChangedRevision: 2206 $
+#$LastChangedBy: steve $
+
+
+from domain import *
+from region import *
+
+Generic_Domain = Domain #Rename
+
+#Shallow water domain
+class Domain(Generic_Domain):
+
+    def __init__(self, coordinates, vertices, boundary = None,
+                 tagged_elements = None, geo_reference = None,
+                 use_inscribed_circle=False):
+
+        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
+        other_quantities = ['elevation', 'friction']
+        Generic_Domain.__init__(self, coordinates, vertices, boundary,
+                                conserved_quantities, other_quantities,
+                                tagged_elements, geo_reference, use_inscribed_circle)
+
+        from anuga.config import minimum_allowed_height, g
+        self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
+
+        self.forcing_terms.append(gravity)
+        self.forcing_terms.append(manning_friction)
+
+        #Realtime visualisation
+        self.visualiser = None
+        self.visualise  = False
+        self.visualise_color_stage = False
+        self.visualise_stage_range = 1.0
+        self.visualise_timer = True
+        self.visualise_range_z = None
+
+        #Stored output
+        self.store = False
+        self.format = 'sww'
+        self.smooth = True
+
+        #Reduction operation for get_vertex_values
+        from anuga.pyvolution.util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
+
+
+
+    def initialise_visualiser(self,scale_z=1.0,rect=None):
+        #Realtime visualisation
+        if self.visualiser is None:
+            from realtime_visualisation_new import Visualiser
+            self.visualiser = Visualiser(self,scale_z,rect)
+            self.visualiser.setup['elevation']=True
+            self.visualiser.updating['stage']=True
+        self.visualise = True
+        if self.visualise_color_stage == True:
+            self.visualiser.coloring['stage'] = True
+            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
+        print 'initialise visualiser'
+        print self.visualiser.setup
+        print self.visualiser.updating
+
+    def check_integrity(self):
+        Generic_Domain.check_integrity(self)
+
+        #Check that we are solving the shallow water wave equation
+
+        msg = 'First conserved quantity must be "stage"'
+        assert self.conserved_quantities[0] == 'stage', msg
+        msg = 'Second conserved quantity must be "xmomentum"'
+        assert self.conserved_quantities[1] == 'xmomentum', msg
+        msg = 'Third conserved quantity must be "ymomentum"'
+        assert self.conserved_quantities[2] == 'ymomentum', msg
+
+    def extrapolate_second_order_sw(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_sw(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges(self)
+
+
+    #FIXME: Under construction
+#     def set_defaults(self):
+#         """Set default values for uninitialised quantities.
+#         This is specific to the shallow water wave equation
+#         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum'
+#         are 0.0. Default for 'stage' is whatever the value of 'elevation'.
+#         """
+
+#         for name in self.other_quantities + self.conserved_quantities:
+#             print name
+#             print self.quantities.keys()
+#             if not self.quantities.has_key(name):
+#                 if name == 'stage':
+
+#                     if self.quantities.has_key('elevation'):
+#                         z = self.quantities['elevation'].vertex_values
+#                         self.set_quantity(name, z)
+#                     else:
+#                         self.set_quantity(name, 0.0)
+#                 else:
+#                     self.set_quantity(name, 0.0)
+
+
+
+#         #Lift negative heights up
+#         #z = self.quantities['elevation'].vertex_values
+#         #w = self.quantities['stage'].vertex_values
+
+#         #h = w-z
+
+#         #for k in range(h.shape[0]):
+#         #    for i in range(3):
+#         #        if h[k, i] < 0.0:
+#         #            w[k, i] = z[k, i]
+
+
+#         #self.quantities['stage'].interpolate()
+
+
+    def evolve(self, yieldstep = None, finaltime = None,
+               skip_initial_step = False):
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call check integrity here rather than from user scripts
+        #self.check_integrity()
+
+        msg = 'Parameter beta_h must be in the interval [0, 1['
+        assert 0 <= self.beta_h < 1.0, msg
+        msg = 'Parameter beta_w must be in the interval [0, 1['
+        assert 0 <= self.beta_w < 1.0, msg
+
+
+        #Initial update of vertex and edge values before any storage
+        #and or visualisation
+        self.distribute_to_vertices_and_edges()
+
+        #Initialise real time viz if requested
+        if self.visualise is True and self.time == 0.0:
+            if self.visualiser is None:
+                self.initialise_visualiser()
+
+            self.visualiser.update_timer()
+            self.visualiser.setup_all()
+
+        #Store model data, e.g. for visualisation
+        if self.store is True and self.time == 0.0:
+            self.initialise_storage()
+            #print 'Storing results in ' + self.writer.filename
+        else:
+            pass
+            #print 'Results will not be stored.'
+            #print 'To store results set domain.store = True'
+            #FIXME: Diagnostic output should be controlled by
+            # a 'verbose' flag living in domain (or in a parent class)
+
+        #Call basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
+                                       skip_initial_step):
+            #Real time viz
+            if self.visualise is True:
+                self.visualiser.update_all()
+                self.visualiser.update_timer()
+
+
+            #Store model data, e.g. for subsequent visualisation
+            if self.store is True:
+                self.store_timestep(self.quantities_to_be_stored)
+
+            #FIXME: Could maybe be taken from specified list
+            #of 'store every step' quantities
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+    def initialise_storage(self):
+        """Create and initialise self.writer object for storing data.
+        Also, save x,y and bed elevation
+        """
+
+        import anuga.pyvolution.data_manager
+
+        #Initialise writer
+        self.writer = data_manager.get_dataobject(self, mode = 'w')
+
+        #Store vertices and connectivity
+        self.writer.store_connectivity()
+
+
+    def store_timestep(self, name):
+        """Store named quantity and time.
+
+        Precondition:
+           self.write has been initialised
+        """
+        self.writer.store_timestep(name)
+
+
+#=============== End of Shallow Water Domain ===============================
+
+
+
+#Rotation of momentum vector
+def rotate(q, normal, direction = 1):
+    """Rotate the momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector.
+
+    If direction is negative the rotation is inverted.
+
+    Input vector is preserved
+
+    This function is specific to the shallow water wave equation
+    """
+
+    from Numeric import zeros, Float
+
+    assert len(q) == 3,\
+           'Vector of conserved quantities must have length 3'\
+           'for 2D shallow water equation'
+
+    try:
+        l = len(normal)
+    except:
+        raise 'Normal vector must be an Numeric array'
+
+    assert l == 2, 'Normal vector must have 2 components'
+
+
+    n1 = normal[0]
+    n2 = normal[1]
+
+    r = zeros(len(q), Float) #Rotated quantities
+    r[0] = q[0]              #First quantity, height, is not rotated
+
+    if direction == -1:
+        n2 = -n2
+
+
+    r[1] =  n1*q[1] + n2*q[2]
+    r[2] = -n2*q[1] + n1*q[2]
+
+    return r
+
+
+
+####################################
+# Flux computation
+def flux_function_central(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from anuga.config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+
+    #Align momentums with x-axis
+    q_left  = rotate(ql, normal, direction = 1)
+    q_right = rotate(qr, normal, direction = 1)
+
+    z = (zl+zr)/2 #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    vh_left  = q_left[2]
+    vh_right = q_right[2]
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+
+    #Flux computation
+
+    #FIXME(Ole): Why is it again that we don't
+    #use uh_left and uh_right directly in the first entries?
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left**2,
+                        u_left*vh_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right**2,
+                        u_right*vh_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+
+        edgeflux = rotate(edgeflux, normal, direction=-1)
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+def erfcc(x):
+
+    from math import fabs, exp
+
+    z=fabs(x)
+    t=1.0/(1.0+0.5*z)
+    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
+         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
+         t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
+    if x < 0.0:
+        result = 2.0-result
+
+    return result
+
+def flux_function_kinetic(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
+
+
+    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from anuga.config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #Align momentums with x-axis
+    q_left  = rotate(ql, normal, direction = 1)
+    q_right = rotate(qr, normal, direction = 1)
+
+    z = (zl+zr)/2 #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    vh_left  = q_left[2]
+    vh_right = q_right[2]
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+
+    #Flux computation
+
+    F_left  = 0.0
+    F_right = 0.0
+    from math import sqrt, pi, exp
+    if h_left > 0.0:
+        F_left = u_left/sqrt(g*h_left)
+    if h_right > 0.0:
+        F_right = u_right/sqrt(g*h_right)
+
+    edgeflux = array([0.0, 0.0, 0.0])
+
+    edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
+          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
+          h_right*u_right/2.0*erfcc(F_right) -  \
+          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
+
+    edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \
+          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
+          (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) -  \
+          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
+
+    edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
+          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
+          vh_right*u_right/2.0*erfcc(F_right) - \
+          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
+
+
+    edgeflux = rotate(edgeflux, normal, direction=-1)
+    max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+flux_function = flux_function_kinetic
+
+def compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    stage = Stage.edge_values
+    xmom =  Xmom.edge_values
+    ymom =  Ymom.edge_values
+    bed =   Bed.edge_values
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    ymom_bdry =  Ymom.boundary_values
+
+    flux = zeros(3, Float) #Work array for summing up fluxes
+
+
+    #Loop
+    timestep = float(sys.maxint)
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        for i in range(3):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i]
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                m = domain.neighbour_edges[k,i]
+                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, 2*i:2*i+2]
+
+            #Flux computation using provided function
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            flux -= edgeflux * domain.edgelengths[k,i]
+
+            #Update optimal_timestep
+            try:
+                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+            except ZeroDivisionError:
+                pass
+
+        #Normalise by area and store for when all conserved
+        #quantities get updated
+        flux /= domain.areas[k]
+        Stage.explicit_update[k] = flux[0]
+        Xmom.explicit_update[k] = flux[1]
+        Ymom.explicit_update[k] = flux[2]
+
+
+    domain.timestep = timestep
+
+#MH090605 The following method belongs to the shallow_water domain class
+#see comments in the corresponding method in shallow_water_ext.c
+def extrapolate_second_order_sw_c(domain):
+    """Wrapper calling C version of extrapolate_second_order_sw
+    """
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    from shallow_water_ext import extrapolate_second_order_sw
+    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
+                                domain.number_of_boundaries,
+                                domain.centroid_coordinates,
+                                Stage.centroid_values,
+                                Xmom.centroid_values,
+                                Ymom.centroid_values,
+                                domain.vertex_coordinates,
+                                Stage.vertex_values,
+                                Xmom.vertex_values,
+                                Ymom.vertex_values)
+
+def compute_fluxes_c(domain):
+    """Wrapper calling C version of compute fluxes
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    timestep = float(sys.maxint)
+    from shallow_water_ext import compute_fluxes
+    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     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,
+                                     domain.already_computed_flux)
+
+
+####################################
+# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+
+def distribute_to_vertices_and_edges(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    from anuga.config import optimised_gradient_limiter
+
+    #Remove very thin layers of water
+    protect_against_infinitesimal_and_negative_heights(domain)
+
+    #Extrapolate all conserved quantities
+    if optimised_gradient_limiter:
+        #MH090605 if second order,
+        #perform the extrapolation and limiting on
+        #all of the conserved quantitie
+
+        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 '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()
+                Q.limit()
+            else:
+                raise '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 protect_against_infinitesimal_and_negative_heights(domain):
+    """Protect against infinitesimal heights and associated high velocities
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values
+    hc = wc - zc  #Water depths at centroids
+
+    #Update
+    for k in range(domain.number_of_elements):
+
+        if hc[k] < domain.minimum_allowed_height:
+            #Control stage
+            if hc[k] < domain.epsilon:
+                wc[k] = zc[k] # Contain 'lost mass' error
+
+            #Control momentum
+            xmomc[k] = ymomc[k] = 0.0
+
+
+def protect_against_infinitesimal_and_negative_heights_c(domain):
+    """Protect against infinitesimal heights and associated high velocities
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values
+
+    from shallow_water_ext import protect
+
+    protect(domain.minimum_allowed_height, domain.epsilon,
+            wc, zc, xmomc, ymomc)
+
+
+
+def h_limiter(domain):
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    hvbar = zeros(hv.shape, Float) #h-limited values
+
+    #Find min and max of this and neighbour's centroid values
+    hmax = zeros(hc.shape, Float)
+    hmin = zeros(hc.shape, Float)
+
+    for k in range(N):
+        hmax[k] = hmin[k] = hc[k]
+        for i in range(3):
+            n = domain.neighbours[k,i]
+            if n >= 0:
+                hn = hc[n] #Neighbour's centroid value
+
+                hmin[k] = min(hmin[k], hn)
+                hmax[k] = max(hmax[k], hn)
+
+
+    #Diffences between centroids and maxima/minima
+    dhmax = hmax - hc
+    dhmin = hmin - hc
+
+    #Deltas between vertex and centroid values
+    dh = zeros(hv.shape, Float)
+    for i in range(3):
+        dh[:,i] = hv[:,i] - hc
+
+    #Phi limiter
+    for k in range(N):
+
+        #Find the gradient limiter (phi) across vertices
+        phi = 1.0
+        for i in range(3):
+            r = 1.0
+            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
+            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
+
+            phi = min( min(r*beta_h, 1), phi )
+
+        #Then update using phi limiter
+        for i in range(3):
+            hvbar[k,i] = hc[k] + phi*dh[k,i]
+
+    return hvbar
+
+
+
+def h_limiter_c(domain):
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+
+    Wrapper for c-extension
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv - zv
+
+    #Call C-extension
+    from shallow_water_ext import h_limiter_sw as h_limiter
+    hvbar = h_limiter(domain, hc, hv)
+
+    return hvbar
+
+
+def balance_deep_and_shallow(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters limiting using w, and stage computed by
+    gradient-limiters limiting using h (h-limiter).
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+
+    The h-limiter is always applied irrespective of the order.
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Limit h
+    hvbar = h_limiter(domain)
+
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]),
+                 abs(zv[k,2]-zc[k]))
+
+
+        hmin = min( hv[k,:] )
+
+        #Create alpha in [0,1], where alpha==0 means using the h-limited
+        #stage and alpha==1 means using the w-limited stage as
+        #computed by the gradient limiter (both 1st or 2nd order)
+
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed
+            alpha = 1.0
+
+	#Let
+	#
+	#  wvi be the w-limited stage (wvi = zvi + hvi)
+	#  wvi- be the h-limited state (wvi- = zvi + hvi-)
+	#
+	#
+	#where i=0,1,2 denotes the vertex ids
+	#
+        #Weighted balance between w-limited and h-limited stage is
+        #
+	#  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+	#
+        #It follows that the updated wvi is
+        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
+        #
+	# Momentum is balanced between constant and limited
+
+
+        #for i in range(3):
+        #    wv[k,i] = zv[k,i] + hvbar[k,i]
+
+	#return
+
+	if alpha < 1:
+
+            for i in range(3):
+                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+            ymomc = domain.quantities['ymomentum'].centroid_values
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+            ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+
+
+def balance_deep_and_shallow_c(domain):
+    """Wrapper for C implementation
+    """
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv - zv
+
+    #Momentums at centroids
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values
+
+    #Momentums at vertices
+    xmomv = domain.quantities['xmomentum'].vertex_values
+    ymomv = domain.quantities['ymomentum'].vertex_values
+
+    #Limit h
+    hvbar = h_limiter(domain)
+
+    #This is how one would make a first order h_limited value
+    #as in the old balancer (pre 17 Feb 2005):
+    #from Numeric import zeros, Float
+    #hvbar = zeros( (len(hc), 3), Float)
+    #for i in range(3):
+    #    hvbar[:,i] = hc[:]
+
+    from shallow_water_ext import balance_deep_and_shallow
+    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
+                             xmomc, ymomc, xmomv, ymomv)
+
+
+
+
+###############################################
+#Boundaries - specific to the shallow water wave equation
+class Reflective_boundary(Boundary):
+    """Reflective boundary returns same conserved quantities as
+    those present in its neighbour volume but reflected.
+
+    This class is specific to the shallow water equation as it
+    works with the momentum quantities assumed to be the second
+    and third conserved quantities.
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for reflective boundary'
+            raise msg
+
+        #Handy shorthands
+        self.stage   = domain.quantities['stage'].edge_values
+        self.xmom    = domain.quantities['xmomentum'].edge_values
+        self.ymom    = domain.quantities['ymomentum'].edge_values
+        self.normals = domain.normals
+
+        from Numeric import zeros, Float
+        self.conserved_quantities = zeros(3, Float)
+
+    def __repr__(self):
+        return 'Reflective_boundary'
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Reflective boundaries reverses the outward momentum
+        of the volume they serve.
+        """
+
+        q = self.conserved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = self.xmom[vol_id, edge_id]
+        q[2] = self.ymom[vol_id, edge_id]
+
+        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+
+
+        r = rotate(q, normal, direction = 1)
+        r[1] = -r[1]
+        q = rotate(r, normal, direction = -1)
+
+        return q
+
+
+
+class Transmissive_Momentum_Set_Stage_boundary(Boundary):
+    """Returns same momentum conserved quantities as
+    those present in its neighbour volume. Sets stage
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain = None, function=None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise msg
+
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise msg
+
+        self.domain   = domain
+        self.function = function
+
+    def __repr__(self):
+        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Transmissive Momentum Set Stage boundaries return the edge momentum
+        values of the volume they serve.
+        """
+
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        value = self.function(self.domain.time)
+        q[0] = value[0]
+        return q
+
+
+        #FIXME: Consider this (taken from File_boundary) to allow
+        #spatial variation
+        #if vol_id is not None and edge_id is not None:
+        #    i = self.boundary_indices[ vol_id, edge_id ]
+        #    return self.F(t, point_id = i)
+        #else:
+        #    return self.F(t)
+
+
+
+class Dirichlet_Discharge_boundary(Boundary):
+    """
+    Sets stage (stage0)
+    Sets momentum (wh0) in the inward normal direction.
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain = None, stage0=None, wh0=None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise msg
+
+        if stage0 is None:
+            raise 'set stage'
+
+        if wh0 is None:
+            wh0 = 0.0
+
+        self.domain   = domain
+        self.stage0  = stage0
+        self.wh0 = wh0
+
+    def __repr__(self):
+        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Set discharge in the (inward) normal direction
+        """
+
+        normal = self.domain.get_normal(vol_id,edge_id)
+        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
+        return q
+
+
+        #FIXME: Consider this (taken from File_boundary) to allow
+        #spatial variation
+        #if vol_id is not None and edge_id is not None:
+        #    i = self.boundary_indices[ vol_id, edge_id ]
+        #    return self.F(t, point_id = i)
+        #else:
+        #    return self.F(t)
+
+
+
+#class Spatio_temporal_boundary(Boundary):
+#    """The spatio-temporal boundary, reads values for the conserved
+#    quantities from an sww NetCDF file, and returns interpolated values
+#    at the midpoints of each associated boundaty segment.
+#    Time dependency is interpolated linearly as in util.File_function.#
+#
+#    Example:
+#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
+#
+#    """
+Spatio_temporal_boundary = File_boundary
+
+
+
+
+#########################
+#Standard forcing terms:
+#
+def gravity(domain):
+    """Apply gravitational pull in the presence of bed slope
+    """
+
+    from anuga.pyvolution.util import gradient
+    from Numeric import zeros, Float, array, sum
+
+    xmom = domain.quantities['xmomentum'].explicit_update
+    ymom = domain.quantities['ymomentum'].explicit_update
+
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
+    h = Stage.edge_values - Elevation.edge_values
+    v = Elevation.vertex_values
+
+    x = domain.get_vertex_coordinates()
+    g = domain.g
+
+    for k in range(domain.number_of_elements):
+        avg_h = sum( h[k,:] )/3
+
+        #Compute bed slope
+        x0, y0, x1, y1, x2, y2 = x[k,:]
+        z0, z1, z2 = v[k,:]
+
+        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+
+        #Update momentum
+        xmom[k] += -g*zx*avg_h
+        ymom[k] += -g*zy*avg_h
+
+
+def gravity_c(domain):
+    """Wrapper calling C version
+    """
+
+    xmom = domain.quantities['xmomentum'].explicit_update
+    ymom = domain.quantities['ymomentum'].explicit_update
+
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
+    h = Stage.edge_values - Elevation.edge_values
+    v = Elevation.vertex_values
+
+    x = domain.get_vertex_coordinates()
+    g = domain.g
+
+
+    from shallow_water_ext import gravity
+    gravity(g, h, v, x, xmom, ymom)
+
+
+def manning_friction(domain):
+    """Apply (Manning) friction to water momentum
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+    vh = domain.quantities['ymomentum'].centroid_values
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g
+
+    for k in range(N):
+        if eta[k] >= eps:
+            if h[k] >= eps:
+            	S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+            	S /= h[k]**(7.0/3)
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+            	ymom_update[k] += S*vh[k]
+
+
+def manning_friction_c(domain):
+    """Wrapper for c version
+    """
+
+
+    xmom = domain.quantities['xmomentum']
+    ymom = domain.quantities['ymomentum']
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+
+    uh = xmom.centroid_values
+    vh = ymom.centroid_values
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = xmom.semi_implicit_update
+    ymom_update = ymom.semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g
+
+    from shallow_water_ext import manning_friction
+    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
+
+
+def linear_friction(domain):
+    """Apply linear friction to water momentum
+
+    Assumes quantity: 'linear_friction' to be present
+    """
+
+    from math import sqrt
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+
+    uh = domain.quantities['xmomentum'].centroid_values
+    vh = domain.quantities['ymomentum'].centroid_values
+    tau = domain.quantities['linear_friction'].centroid_values
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g #Not necessary? Why was this added?
+
+    for k in range(N):
+        if tau[k] >= eps:
+            if h[k] >= eps:
+            	S = -tau[k]/h[k]
+
+            	#Update momentum
+            	xmom_update[k] += S*uh[k]
+            	ymom_update[k] += S*vh[k]
+
+
+
+def check_forcefield(f):
+    """Check that f is either
+    1: a callable object f(t,x,y), where x and y are vectors
+       and that it returns an array or a list of same length
+       as x and y
+    2: a scalar
+    """
+
+    from Numeric import ones, Float, array
+
+
+    if callable(f):
+        N = 3
+        x = ones(3, Float)
+        y = ones(3, Float)
+        try:
+            q = f(1.0, x=x, y=y)
+        except Exception, e:
+            msg = 'Function %s could not be executed:\n%s' %(f, e)
+	    #FIXME: Reconsider this semantics
+            raise msg
+
+        try:
+            q = array(q).astype(Float)
+        except:
+            msg = 'Return value from vector function %s could ' %f
+            msg += 'not be converted into a Numeric array of floats.\n'
+            msg += 'Specified function should return either list or array.'
+            raise msg
+
+        #Is this really what we want?
+        msg = 'Return vector from function %s ' %f
+        msg += 'must have same lenght as input vectors'
+        assert len(q) == N, msg
+
+    else:
+        try:
+            f = float(f)
+        except:
+            msg = 'Force field %s must be either a scalar' %f
+            msg += ' or a vector function'
+            raise msg
+    return f
+
+
+class Wind_stress:
+    """Apply wind stress to water momentum in terms of
+    wind speed [m/s] and wind direction [degrees]
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Initialise windfield from wind speed s [m/s]
+        and wind direction phi [degrees]
+
+        Inputs v and phi can be either scalars or Python functions, e.g.
+
+        W = Wind_stress(10, 178)
+
+        #FIXME - 'normal' degrees are assumed for now, i.e. the
+        vector (1,0) has zero degrees.
+        We may need to convert from 'compass' degrees later on and also
+        map from True north to grid north.
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def speed(t,x,y):
+            ...
+            return s
+
+        def angle(t,x,y):
+            ...
+            return phi
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        W = Wind_stress(speed, angle)
+
+        The instantiated object W can be appended to the list of
+        forcing_terms as in
+
+        Alternatively, one vector valued function for (speed, angle)
+        can be applied, providing both quantities simultaneously.
+        As in
+        W = Wind_stress(F), where returns (speed, angle) for each t.
+
+        domain.forcing_terms.append(W)
+        """
+
+        from anuga.config import rho_a, rho_w, eta_w
+        from Numeric import array, Float
+
+        if len(args) == 2:
+            s = args[0]
+            phi = args[1]
+        elif len(args) == 1:
+            #Assume vector function returning (s, phi)(t,x,y)
+            vector_function = args[0]
+            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+        else:
+           #Assume info is in 2 keyword arguments
+
+           if len(kwargs) == 2:
+               s = kwargs['s']
+               phi = kwargs['phi']
+           else:
+               raise 'Assumes two keyword arguments: s=..., phi=....'
+
+        self.speed = check_forcefield(s)
+        self.phi = check_forcefield(phi)
+
+        self.const = eta_w*rho_a/rho_w
+
+
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import Float, ones, ArrayType
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        N = domain.number_of_elements
+        t = domain.time
+
+        if callable(self.speed):
+            xc = domain.get_centroid_coordinates()
+            s_vec = self.speed(t, xc[:,0], xc[:,1])
+        else:
+            #Assume s is a scalar
+
+            try:
+                s_vec = self.speed * ones(N, Float)
+            except:
+                msg = 'Speed must be either callable or a scalar: %s' %self.s
+                raise msg
+
+
+        if callable(self.phi):
+            xc = domain.get_centroid_coordinates()
+            phi_vec = self.phi(t, xc[:,0], xc[:,1])
+        else:
+            #Assume phi is a scalar
+
+            try:
+                phi_vec = self.phi * ones(N, Float)
+            except:
+                msg = 'Angle must be either callable or a scalar: %s' %self.phi
+                raise msg
+
+        assign_windfield_values(xmom_update, ymom_update,
+                                s_vec, phi_vec, self.const)
+
+
+def assign_windfield_values(xmom_update, ymom_update,
+                            s_vec, phi_vec, const):
+    """Python version of assigning wind field to update vectors.
+    A c version also exists (for speed)
+    """
+    from math import pi, cos, sin, sqrt
+
+    N = len(s_vec)
+    for k in range(N):
+        s = s_vec[k]
+        phi = phi_vec[k]
+
+        #Convert to radians
+        phi = phi*pi/180
+
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress
+        S = const * sqrt(u**2 + v**2)
+        xmom_update[k] += S*u
+        ymom_update[k] += S*v
+
+
+
+##############################
+#OBSOLETE STUFF
+
+def balance_deep_and_shallow_old(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters and stage computed as constant height above bed.
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+    """
+
+    #OBSOLETE
+
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+
+    #Computed linear combination between constant stages and and
+    #stages parallel to the bed elevation.
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]),
+                 abs(zv[k,2]-zc[k]))
+
+
+        hmin = min( hv[k,:] )
+
+        #Create alpha in [0,1], where alpha==0 means using shallow
+        #first order scheme and alpha==1 means using the stage w as
+        #computed by the gradient limiter (1st or 2nd order)
+        #
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed
+            alpha = 1.0
+
+        #Weighted balance between stage parallel to bed elevation
+        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
+	#order gradient limiter
+        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
+        #
+        #It follows that the updated wvi is
+        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
+        #  zvi + hc + alpha*(hvi - hc)
+        #
+        #Note that hvi = zc+hc-zvi in the first order case (constant).
+
+	if alpha < 1:
+            for i in range(3):
+                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+            ymomc = domain.quantities['ymomentum'].centroid_values
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+            ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+
+
+
+
+
+###########################
+###########################
+#Geometries
+
+
+#FIXME: Rethink this way of creating values.
+
+
+class Weir:
+    """Set a bathymetry for weir with a hole and a downstream gutter
+    x,y are assumed to be in the unit square
+    """
+
+    def __init__(self, stage):
+        self.inflow_stage = stage
+
+    def __call__(self, x, y):
+        from Numeric import zeros, Float
+        from math import sqrt
+
+        N = len(x)
+        assert N == len(y)
+
+        z = zeros(N, Float)
+        for i in range(N):
+            z[i] = -x[i]/2  #General slope
+
+            #Flattish bit to the left
+            if x[i] < 0.3:
+                z[i] = -x[i]/10
+
+            #Weir
+            if x[i] >= 0.3 and x[i] < 0.4:
+                z[i] = -x[i]+0.9
+
+            #Dip
+            x0 = 0.6
+            #depth = -1.3
+            depth = -1.0
+            #plateaux = -0.9
+            plateaux = -0.6
+            if y[i] < 0.7:
+                if x[i] > x0 and x[i] < 0.9:
+                    z[i] = depth
+
+                #RHS plateaux
+                if x[i] >= 0.9:
+                    z[i] = plateaux
+
+
+            elif y[i] >= 0.7 and y[i] < 1.5:
+                #Restrict and deepen
+                if x[i] >= x0 and x[i] < 0.8:
+                    z[i] = depth-(y[i]/3-0.3)
+                    #z[i] = depth-y[i]/5
+                    #z[i] = depth
+                elif x[i] >= 0.8:
+                    #RHS plateaux
+                    z[i] = plateaux
+
+            elif y[i] >= 1.5:
+                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
+                    #Widen up and stay at constant depth
+                    z[i] = depth-1.5/5
+                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
+                    #RHS plateaux
+                    z[i] = plateaux
+
+
+            #Hole in weir (slightly higher than inflow condition)
+            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.02
+
+            #Channel behind weir
+            x0 = 0.5
+            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.02
+
+            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
+                #Flatten it out towards the end
+                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
+
+            #Hole to the east
+            x0 = 1.1; y0 = 0.35
+            #if x[i] < -0.2 and y < 0.5:
+            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
+                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
+
+            #Tiny channel draining hole
+            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
+                z[i] = -0.9 #North south
+
+            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
+                z[i] = -1.0 + (x[i]-0.9)/3 #East west
+
+
+
+            #Stuff not in use
+
+            #Upward slope at inlet to the north west
+            #if x[i] < 0.0: # and y[i] > 0.5:
+            #    #z[i] = -y[i]+0.5  #-x[i]/2
+            #    z[i] = x[i]/4 - y[i]**2 + 0.5
+
+            #Hole to the west
+            #x0 = -0.4; y0 = 0.35 # center
+            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
+            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
+
+
+
+
+
+        return z/2
+
+class Weir_simple:
+    """Set a bathymetry for weir with a hole and a downstream gutter
+    x,y are assumed to be in the unit square
+    """
+
+    def __init__(self, stage):
+        self.inflow_stage = stage
+
+    def __call__(self, x, y):
+        from Numeric import zeros, Float
+
+        N = len(x)
+        assert N == len(y)
+
+        z = zeros(N, Float)
+        for i in range(N):
+            z[i] = -x[i]  #General slope
+
+            #Flat bit to the left
+            if x[i] < 0.3:
+                z[i] = -x[i]/10  #General slope
+
+            #Weir
+            if x[i] > 0.3 and x[i] < 0.4:
+                z[i] = -x[i]+0.9
+
+            #Dip
+            if x[i] > 0.6 and x[i] < 0.9:
+                z[i] = -x[i]-0.5  #-y[i]/5
+
+            #Hole in weir (slightly higher than inflow condition)
+            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.05
+
+
+        return z/2
+
+
+
+class Constant_stage:
+    """Set an initial condition with constant stage
+    """
+    def __init__(self, s):
+        self.s = s
+
+    def __call__(self, x, y):
+        return self.s
+
+class Constant_height:
+    """Set an initial condition with constant water height, e.g
+    stage s = z+h
+    """
+
+    def __init__(self, W, h):
+        self.W = W
+        self.h = h
+
+    def __call__(self, x, y):
+        if self.W is None:
+            from Numeric import ones, Float
+            return self.h*ones(len(x), Float)
+        else:
+            return self.W(x,y) + self.h
+
+
+
+
+
+
+
+##############################################
+#Initialise module
+
+
+from anuga.utilities import compile
+if compile.can_use_C_extension('shallow_water_kinetic_ext.c'):
+    #Replace python version with c implementations
+
+    from shallow_water_kinetic_ext import rotate, assign_windfield_values
+    compute_fluxes = compute_fluxes_c
+    extrapolate_second_order_sw=extrapolate_second_order_sw_c
+    gravity = gravity_c
+    manning_friction = manning_friction_c
+    h_limiter = h_limiter_c
+    balance_deep_and_shallow = balance_deep_and_shallow_c
+    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
+
+
+    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c #(like MH's)
+
+
+
+#Optimisation with psyco
+from anuga.config import use_psyco
+if use_psyco:
+    try:
+        import psyco
+    except:
+        import os
+        if os.name == 'posix' and os.uname()[4] == 'x86_64':
+            pass
+            #Psyco isn't supported on 64 bit systems, but it doesn't matter
+        else:
+            msg = 'WARNING: psyco (speedup) could not import'+\
+                  ', you may want to consider installing it'
+            print msg
+    else:
+        psyco.bind(Domain.distribute_to_vertices_and_edges)
+        psyco.bind(Domain.compute_fluxes)
+
+if __name__ == "__main__":
+    pass
+
+# Profiling stuff
+import profile
+profiler = profile.Profile()
Index: /anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c	(revision 3698)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c	(revision 3698)
@@ -0,0 +1,1442 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.3):
+//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
+//  gcc -shared domain_ext.o  -o domain_ext.so
+//
+// or use python compile.py
+//
+// See the module shallow_water.py
+//
+//
+// Ole Nielsen, GA 2004
+
+
+#include "Python.h"
+#include "Numeric/arrayobject.h"
+#include "math.h"
+#include <stdio.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 momentum component q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector (n1, n2).
+
+    Result is returned in array 3x1 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];  //uh momentum
+  q2 = q[2];  //vh momentum
+
+  //Rotate
+  q[1] =  n1*q1 + n2*q2;
+  q[2] = -n2*q1 + n1*q2;
+
+  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)
+  if (dq0>=0.0){
+    if (dq1>=dq2){
+      if (dq1>=0.0)
+	*qmax=dq0+dq1;
+      else
+	*qmax=dq0;
+      if ((*qmin=dq0+dq2)<0)
+	;//qmin is already set to correct value
+      else
+	*qmin=0.0;
+    }
+    else{//dq1<dq2
+      if (dq2>0)
+	*qmax=dq0+dq2;
+      else
+	*qmax=dq0;
+      if ((*qmin=dq0+dq1)<0)
+	;//qmin is the correct value
+      else
+	*qmin=0.0;
+    }
+  }
+  else{//dq0<0
+    if (dq1<=dq2){
+      if (dq1<0.0)
+	*qmin=dq0+dq1;
+      else
+	*qmin=dq0;
+      if ((*qmax=dq0+dq2)>0.0)
+	;//qmax is already set to the correct value
+      else
+	*qmax=0.0;
+    }
+    else{//dq1>dq2
+      if (dq2<0.0)
+	*qmin=dq0+dq2;
+      else
+	*qmin=dq0;
+      if ((*qmax=dq0+dq1)>0.0)
+	;//qmax is already set to the correct value
+      else
+	*qmax=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.
+  //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);
+  for (i=0;i<3;i++)
+    dqv[i]=dqv[i]*phi;
+  return 0;
+}
+
+// Computational function for flux computation (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 g,
+		  double *edgeflux, double *max_speed) {
+
+  /*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, z;
+  double q_left_copy[3], q_right_copy[3];
+  double flux_right[3], flux_left[3];
+
+  //Copy conserved quantities to protect from modification
+  for (i=0; i<3; i++) {
+    q_left_copy[i] = q_left[i];
+    q_right_copy[i] = q_right[i];
+  }
+
+  //Align x- and y-momentum with x-axis
+  _rotate(q_left_copy, n1, n2);
+  _rotate(q_right_copy, n1, n2);
+
+  z = (z_left+z_right)/2; //Take average of field values
+
+  //Compute speeds in x-direction
+  w_left = q_left_copy[0];              // h+z
+  h_left = w_left-z;
+  uh_left = q_left_copy[1];
+
+  if (h_left < epsilon) {
+    h_left = 0.0;  //Could have been negative
+    u_left = 0.0;
+  } else {
+    u_left = uh_left/h_left;
+  }
+
+  w_right = q_right_copy[0];
+  h_right = w_right-z;
+  uh_right = q_right_copy[1];
+
+  if (h_right < epsilon) {
+    h_right = 0.0; //Could have been negative
+    u_right = 0.0;
+  } else {
+    u_right = uh_right/h_right;
+  }
+
+  //Momentum in y-direction
+  vh_left  = q_left_copy[2];
+  vh_right = q_right_copy[2];
+
+
+  //Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);
+
+  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+  if (s_max < 0.0) s_max = 0.0;
+
+  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+  if (s_min > 0.0) s_min = 0.0;
+
+  //Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*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 == 0.0) {
+    for (i=0; i<3; i++) edgeflux[i] = 0.0;
+    *max_speed = 0.0;
+  } else {
+    for (i=0; i<3; i++) {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
+      edgeflux[i] /= denom;
+    }
+
+    //Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+
+    //Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+  return 0;
+}
+
+double erfcc(double x){
+    double z,t,result;
+
+    z=fabs(x);
+    t=1.0/(1.0+0.5*z);
+    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
+         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
+         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
+    if (x < 0.0) result = 2.0-result;
+
+    return result;
+    }
+
+
+
+// Computational function for flux computation (using stage w=z+h)
+int flux_function_kinetic(double *q_left, double *q_right,
+		  double z_left, double z_right,
+		  double n1, double n2,
+		  double epsilon, double g,
+		  double *edgeflux, double *max_speed) {
+
+  /*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
+
+    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
+  */
+
+  int i;
+
+  double w_left, h_left, uh_left, vh_left, u_left, F_left;
+  double w_right, h_right, uh_right, vh_right, u_right, F_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double z;
+  double q_left_copy[3], q_right_copy[3];
+
+
+  //Copy conserved quantities to protect from modification
+  for (i=0; i<3; i++) {
+    q_left_copy[i] = q_left[i];
+    q_right_copy[i] = q_right[i];
+  }
+
+  //Align x- and y-momentum with x-axis
+  _rotate(q_left_copy, n1, n2);
+  _rotate(q_right_copy, n1, n2);
+
+  z = (z_left+z_right)/2; //Take average of field values
+
+  //Compute speeds in x-direction
+  w_left = q_left_copy[0];              // h+z
+  h_left = w_left-z;
+  uh_left = q_left_copy[1];
+
+  if (h_left < epsilon) {
+    h_left = 0.0;  //Could have been negative
+    u_left = 0.0;
+  } else {
+    u_left = uh_left/h_left;
+  }
+
+  w_right = q_right_copy[0];
+  h_right = w_right-z;
+  uh_right = q_right_copy[1];
+
+  if (h_right < epsilon) {
+    h_right = 0.0; //Could have been negative
+    u_right = 0.0;
+  } else {
+    u_right = uh_right/h_right;
+  }
+
+  //Momentum in y-direction
+  vh_left  = q_left_copy[2];
+  vh_right = q_right_copy[2];
+
+
+  //Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);
+
+  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+  if (s_max < 0.0) s_max = 0.0;
+
+  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+  if (s_min > 0.0) s_min = 0.0;
+
+
+  F_left  = 0.0;
+  F_right = 0.0;
+  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
+  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
+
+  for (i=0; i<3; i++) edgeflux[i] = 0.0;
+  *max_speed = 0.0;
+
+  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
+          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
+          h_right*u_right/2.0*erfcc(F_right) -  \
+          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
+
+  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
+          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
+          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
+          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
+
+  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
+          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
+          vh_right*u_right/2.0*erfcc(F_right) - \
+          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
+
+  //Maximal wavespeed
+  *max_speed = max(fabs(s_max), fabs(s_min));
+
+  //Rotate back
+  _rotate(edgeflux, n1, -n2);
+
+  return 0;
+}
+
+
+void _manning_friction(double g, double eps, int N,
+		       double* w, double* z,
+		       double* uh, double* vh,
+		       double* eta, double* xmom, double* ymom) {
+
+  int k;
+  double S, h;
+
+  for (k=0; k<N; k++) {
+    if (eta[k] > eps) {
+      h = w[k]-z[k];
+      if (h >= eps) {
+	S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+	//S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+	S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+	//S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
+
+	//Update momentum
+	xmom[k] += S*uh[k];
+	ymom[k] += S*vh[k];
+      }
+    }
+  }
+}
+
+
+
+int _balance_deep_and_shallow(int N,
+			      double* wc,
+			      double* zc,
+			      double* hc,
+			      double* wv,
+			      double* zv,
+			      double* hv,
+			      double* hvbar,
+			      double* xmomc,
+			      double* ymomc,
+			      double* xmomv,
+			      double* ymomv) {
+
+  int k, k3, i;
+  double dz, hmin, alpha;
+
+  //Compute linear combination between w-limited stages and
+  //h-limited stages close to the bed elevation.
+
+  for (k=0; k<N; k++) {
+    // Compute maximal variation in bed elevation
+    // This quantitiy is
+    //     dz = max_i abs(z_i - z_c)
+    // and it is independent of dimension
+    // In the 1d case zc = (z0+z1)/2
+    // In the 2d case zc = (z0+z1+z2)/3
+
+    k3 = 3*k;
+
+    //FIXME: Try with this one precomputed
+    dz = 0.0;
+    hmin = hv[k3];
+    for (i=0; i<3; i++) {
+      dz = max(dz, fabs(zv[k3+i]-zc[k]));
+      hmin = min(hmin, hv[k3+i]);
+    }
+
+
+    //Create alpha in [0,1], where alpha==0 means using the h-limited
+    //stage and alpha==1 means using the w-limited stage as
+    //computed by the gradient limiter (both 1st or 2nd order)
+    //
+    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
+    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+
+    if (dz > 0.0)
+      //if (hmin<0.0)
+      //	alpha = 0.0;
+      //else
+      //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
+      alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 );
+    else
+      alpha = 1.0;  //Flat bed
+
+    //alpha = 1.0;
+
+    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
+
+    //	Let
+    //
+    //	  wvi be the w-limited stage (wvi = zvi + hvi)
+    //	  wvi- be the h-limited state (wvi- = zvi + hvi-)
+    //
+    //
+    //	where i=0,1,2 denotes the vertex ids
+    //
+    //  Weighted balance between w-limited and h-limited stage is
+    //
+    //	  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+    //
+    //  It follows that the updated wvi is
+    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
+    //
+    //	 Momentum is balanced between constant and limited
+
+    if (alpha < 1) {
+      for (i=0; i<3; i++) {
+         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
+
+	//Update momentum as a linear combination of
+	//xmomc and ymomc (shallow) and momentum
+	//from extrapolator xmomv and ymomv (deep).
+	xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
+	ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+      }
+    }
+  }
+  return 0;
+}
+
+
+
+int _protect(int N,
+	     double minimum_allowed_height,
+	     double epsilon,
+	     double* wc,
+	     double* zc,
+	     double* xmomc,
+	     double* ymomc) {
+
+  int k;
+  double hc;
+
+  //Protect against initesimal and negative heights
+
+  for (k=0; k<N; k++) {
+    hc = wc[k] - zc[k];
+
+    if (hc < minimum_allowed_height) {
+      if (hc < epsilon) wc[k] = zc[k]; //Contain 'lost mass' error
+      xmomc[k] = 0.0;
+      ymomc[k] = 0.0;
+    }
+
+  }
+  return 0;
+}
+
+
+
+
+int _assign_wind_field_values(int N,
+			      double* xmom_update,
+			      double* ymom_update,
+			      double* s_vec,
+			      double* phi_vec,
+			      double cw) {
+
+  //Assign windfield values to momentum updates
+
+  int k;
+  double S, s, phi, u, v;
+
+  for (k=0; k<N; k++) {
+
+    s = s_vec[k];
+    phi = phi_vec[k];
+
+    //Convert to radians
+    phi = phi*pi/180;
+
+    //Compute velocity vector (u, v)
+    u = s*cos(phi);
+    v = s*sin(phi);
+
+    //Compute wind stress
+    S = cw * sqrt(u*u + v*v);
+    xmom_update[k] += S*u;
+    ymom_update[k] += S*v;
+  }
+  return 0;
+}
+
+
+
+///////////////////////////////////////////////////////////////////
+// Gateways to Python
+
+PyObject *gravity(PyObject *self, PyObject *args) {
+  //
+  //  gravity(g, h, v, x, xmom, ymom)
+  //
+
+
+  PyArrayObject *h, *v, *x, *xmom, *ymom;
+  int k, i, N, k3, k6;
+  double g, avg_h, zx, zy;
+  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
+
+  if (!PyArg_ParseTuple(args, "dOOOOO",
+			&g, &h, &v, &x,
+			&xmom, &ymom))
+    return NULL;
+
+  N = h -> dimensions[0];
+  for (k=0; k<N; k++) {
+    k3 = 3*k;  // base index
+    k6 = 6*k;  // base index
+
+    avg_h = 0.0;
+    for (i=0; i<3; i++) {
+      avg_h += ((double *) h -> data)[k3+i];
+    }
+    avg_h /= 3;
+
+
+    //Compute bed slope
+    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];
+
+
+    z0 = ((double*) v -> data)[k3 + 0];
+    z1 = ((double*) v -> data)[k3 + 1];
+    z2 = ((double*) v -> data)[k3 + 2];
+
+    _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 *manning_friction(PyObject *self, PyObject *args) {
+  //
+  // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
+			&g, &eps, &w, &z, &uh, &vh, &eta,
+			&xmom, &ymom))
+    return NULL;
+
+  N = w -> dimensions[0];
+  _manning_friction(g, eps, N,
+		    (double*) w -> data,
+		    (double*) z -> data,
+		    (double*) uh -> data,
+		    (double*) vh -> data,
+		    (double*) eta -> data,
+		    (double*) xmom -> data,
+		    (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
+
+PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
+  /*Compute the vertex values based on a linear reconstruction on each triangle
+    These values are calculated as follows:
+    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
+    formed by the centroids of its three neighbours.
+    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
+    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
+    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
+    original triangle has gradient (a,b).
+    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
+    find_qmin_and_qmax and limit_gradient
+
+    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.vertex_coordinates,
+                                Stage.vertex_values,
+                                Xmom.vertex_values,
+                                Ymom.vertex_values)
+
+    Post conditions:
+            The vertices 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,
+    *vertex_coordinates,
+    *stage_vertex_values,
+    *xmom_vertex_values,
+    *ymom_vertex_values;
+  PyObject *domain, *Tmp;
+  //Local variables
+  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
+  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
+  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;
+  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
+  //by which these jumps are limited
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
+			&domain,
+			&surrogate_neighbours,
+			&number_of_boundaries,
+			&centroid_coordinates,
+			&stage_centroid_values,
+			&xmom_centroid_values,
+			&ymom_centroid_values,
+			&vertex_coordinates,
+			&stage_vertex_values,
+			&xmom_vertex_values,
+			&ymom_vertex_values)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+
+  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
+  Tmp = PyObject_GetAttrString(domain, "beta_w");
+  if (!Tmp)
+    return NULL;
+  beta_w = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
+  number_of_elements = stage_centroid_values -> dimensions[0];
+  for (k=0; k<number_of_elements; k++) {
+    k3=k*3;
+    k6=k*6;
+
+    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
+      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
+      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
+      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
+      continue;
+    }
+    else {//we will need centroid coordinates and vertex coordinates of the triangle
+      //get the vertex coordinates of the FV triangle
+      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
+      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
+      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
+      //get the centroid coordinates of the FV triangle
+      coord_index=2*k;
+      x=((double *)centroid_coordinates->data)[coord_index];
+      y=((double *)centroid_coordinates->data)[coord_index+1];
+      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
+      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
+      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
+    }
+    if (((long *)number_of_boundaries->data)[k]<=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=((long *)surrogate_neighbours->data)[k3];
+      k1=((long *)surrogate_neighbours->data)[k3+1];
+      k2=((long *)surrogate_neighbours->data)[k3+2];
+      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
+      coord_index=2*k0;
+      x0=((double *)centroid_coordinates->data)[coord_index];
+      y0=((double *)centroid_coordinates->data)[coord_index+1];
+      coord_index=2*k1;
+      x1=((double *)centroid_coordinates->data)[coord_index];
+      y1=((double *)centroid_coordinates->data)[coord_index+1];
+      coord_index=2*k2;
+      x2=((double *)centroid_coordinates->data)[coord_index];
+      y2=((double *)centroid_coordinates->data)[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
+      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
+      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
+      if (area2<=0)
+	return NULL;
+
+      //### stage ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
+      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
+      //calculate the gradient of stage on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri 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);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
+
+      //### xmom ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
+      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
+      //calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri 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);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
+
+      //### ymom ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
+      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
+      //calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri 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);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
+    }//if (number_of_boundaries[k]<=1)
+    else{//number_of_boundaries==2
+      //one internal neighbour and gradient is in direction of the neighbour's centroid
+      //find the only internal neighbour
+      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
+	if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
+	  break;
+      }
+      if ((k2==k3+3))//if we didn't find an internal neighbour
+	return NULL;//error
+      k1=((long *)surrogate_neighbours->data)[k2];
+      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
+      coord_index=2*k1;
+      x1=((double *)centroid_coordinates->data)[coord_index];
+      y1=((double *)centroid_coordinates->data)[coord_index+1];
+      //compute x- and y- distances between the centroid of the FV triangle 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=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
+
+      //## xmom ###
+      //compute differentials
+      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
+
+      //## ymom ###
+      //compute differentials
+      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex 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_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
+    }//else [number_of_boudaries==2]
+  }//for k=0 to number_of_elements-1
+  return Py_BuildValue("");
+}//extrapolate_second-order_sw
+
+PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
+  //
+  // r = rotate(q, normal, direction=1)
+  //
+  // Where q is assumed to be a Float numeric array of length 3 and
+  // normal a Float numeric array of length 2.
+
+
+  PyObject *Q, *Normal;
+  PyArrayObject *q, *r, *normal;
+
+  static char *argnames[] = {"q", "normal", "direction", NULL};
+  int dimensions[1], i, direction=1;
+  double n1, n2;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
+				   &Q, &Normal, &direction))
+    return NULL;
+
+  //Input checks (convert sequences into numeric arrays)
+  q = (PyArrayObject *)
+    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
+  normal = (PyArrayObject *)
+    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
+
+
+  if (normal -> dimensions[0] != 2) {
+    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
+    return NULL;
+  }
+
+  //Allocate space for return vector r (don't DECREF)
+  dimensions[0] = 3;
+  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
+
+  //Copy
+  for (i=0; i<3; i++) {
+    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
+  }
+
+  //Get normal and direction
+  n1 = ((double *) normal -> data)[0];
+  n2 = ((double *) normal -> data)[1];
+  if (direction == -1) n2 = -n2;
+
+  //Rotate
+  _rotate((double *) r -> data, n1, n2);
+
+  //Release numeric arrays
+  Py_DECREF(q);
+  Py_DECREF(normal);
+
+  //return result using PyArray to avoid memory leak
+  return PyArray_Return(r);
+}
+
+PyObject *compute_fluxes(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"
+
+    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.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     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)
+
+
+    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,
+    *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;//tracks whether the flux across an edge has already been computed
+
+
+  //Local variables
+  double timestep, max_speed, epsilon, g;
+  double normal[2], ql[3], qr[3], zl, zr;
+  double edgeflux[3]; //Work arrays for summing up fluxes
+
+  int number_of_elements, k, i, m, n;
+  int ki, nm=0, ki2; //Index shorthands
+  static long call=1;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
+			&timestep,
+			&epsilon,
+			&g,
+			&neighbours,
+			&neighbour_edges,
+			&normals,
+			&edgelengths, &radii, &areas,
+			&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)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+  number_of_elements = stage_edge_values -> dimensions[0];
+  call++;//a static local variable to which already_computed_flux is compared
+  //set explicit_update to zero for all conserved_quantities.
+  //This assumes compute_fluxes called before forcing terms
+  for (k=0; k<number_of_elements; k++) {
+    ((double *) stage_explicit_update -> data)[k]=0.0;
+    ((double *) xmom_explicit_update -> data)[k]=0.0;
+    ((double *) ymom_explicit_update -> data)[k]=0.0;
+  }
+  //Loop through neighbours and compute edge flux for each
+  for (k=0; k<number_of_elements; k++) {
+    for (i=0; i<3; i++) {
+      ki = k*3+i;
+      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
+	continue;
+      ql[0] = ((double *) stage_edge_values -> data)[ki];
+      ql[1] = ((double *) xmom_edge_values -> data)[ki];
+      ql[2] = ((double *) ymom_edge_values -> data)[ki];
+      zl =    ((double *) bed_edge_values -> data)[ki];
+
+      //Quantities at neighbour on nearest face
+      n = ((long *) neighbours -> data)[ki];
+      if (n < 0) {
+	m = -n-1; //Convert negative flag to index
+	qr[0] = ((double *) stage_boundary_values -> data)[m];
+	qr[1] = ((double *) xmom_boundary_values -> data)[m];
+	qr[2] = ((double *) ymom_boundary_values -> data)[m];
+	zr = zl; //Extend bed elevation to boundary
+      } else {
+	m = ((long *) neighbour_edges -> data)[ki];
+	nm = n*3+m;
+	qr[0] = ((double *) stage_edge_values -> data)[nm];
+	qr[1] = ((double *) xmom_edge_values -> data)[nm];
+	qr[2] = ((double *) ymom_edge_values -> data)[nm];
+	zr =    ((double *) bed_edge_values -> data)[nm];
+      }
+      // Outward pointing normal vector
+      // normal = domain.normals[k, 2*i:2*i+2]
+      ki2 = 2*ki; //k*6 + i*2
+      normal[0] = ((double *) normals -> data)[ki2];
+      normal[1] = ((double *) normals -> data)[ki2+1];
+      //Edge flux computation
+      flux_function_kinetic(ql, qr, zl, zr,
+		    normal[0], normal[1],
+		    epsilon, g,
+		    edgeflux, &max_speed);
+      //update triangle k
+      ((long *) already_computed_flux->data)[ki]=call;
+      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
+      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
+      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
+      //update the neighbour n
+      if (n>=0){
+	((long *) already_computed_flux->data)[nm]=call;
+	((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
+	((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
+	((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
+      }
+      ///for (j=0; j<3; j++) {
+	///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
+	///}
+	//Update timestep
+	//timestep = min(timestep, domain.radii[k]/max_speed)
+	//FIXME: SR Add parameter for CFL condition
+	if (max_speed > epsilon) {
+	  timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+	  //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
+	  if (n>=0)
+	    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
+	}
+    } // end for i
+    //Normalise by area and store for when all conserved
+    //quantities get updated
+    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+  } //end for k
+  return Py_BuildValue("d", timestep);
+}
+
+PyObject *protect(PyObject *self, PyObject *args) {
+  //
+  //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
+
+
+  PyArrayObject
+  *wc,            //Stage at centroids
+  *zc,            //Elevation at centroids
+  *xmomc,         //Momentums at centroids
+  *ymomc;
+
+
+  int N;
+  double minimum_allowed_height, epsilon;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddOOOO",
+			&minimum_allowed_height,
+			&epsilon,
+			&wc, &zc, &xmomc, &ymomc))
+    return NULL;
+
+  N = wc -> dimensions[0];
+
+  _protect(N,
+	   minimum_allowed_height,
+	   epsilon,
+	   (double*) wc -> data,
+	   (double*) zc -> data,
+	   (double*) xmomc -> data,
+	   (double*) ymomc -> data);
+
+  return Py_BuildValue("");
+}
+
+
+
+PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
+  //
+  //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
+  //                             xmomc, ymomc, xmomv, ymomv)
+
+
+  PyArrayObject
+    *wc,            //Stage at centroids
+    *zc,            //Elevation at centroids
+    *hc,            //Height at centroids
+    *wv,            //Stage at vertices
+    *zv,            //Elevation at vertices
+    *hv,            //Depths at vertices
+    *hvbar,         //h-Limited depths at vertices
+    *xmomc,         //Momentums at centroids and vertices
+    *ymomc,
+    *xmomv,
+    *ymomv;
+
+  int N; //, err;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
+			&wc, &zc, &hc,
+			&wv, &zv, &hv, &hvbar,
+			&xmomc, &ymomc, &xmomv, &ymomv))
+    return NULL;
+
+  N = wc -> dimensions[0];
+
+  _balance_deep_and_shallow(N,
+			    (double*) wc -> data,
+			    (double*) zc -> data,
+			    (double*) hc -> data,
+			    (double*) wv -> data,
+			    (double*) zv -> data,
+			    (double*) hv -> data,
+			    (double*) hvbar -> data,
+                            (double*) xmomc -> data,
+			    (double*) ymomc -> data,
+			    (double*) xmomv -> data,
+			    (double*) ymomv -> data);
+
+
+  return Py_BuildValue("");
+}
+
+
+
+PyObject *h_limiter(PyObject *self, PyObject *args) {
+
+  PyObject *domain, *Tmp;
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
+    *hvbar,   //Limited depth at vertices (return values)
+    *neighbours;
+
+  int k, i, n, N, k3;
+  int dimensions[2];
+  double beta_h; //Safety factor (see config.py)
+  double *hmin, *hmax, hn;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
+    return NULL;
+
+  neighbours = get_consecutive_array(domain, "neighbours");
+
+  //Get safety factor beta_h
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
+
+  N = hc -> dimensions[0];
+
+  //Create hvbar
+  dimensions[0] = N;
+  dimensions[1] = 3;
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+
+
+  //Find min and max of this and neighbour's centroid values
+  hmin = malloc(N * sizeof(double));
+  hmax = malloc(N * sizeof(double));
+  for (k=0; k<N; k++) {
+    k3=k*3;
+
+    hmin[k] = ((double*) hc -> data)[k];
+    hmax[k] = hmin[k];
+
+    for (i=0; i<3; i++) {
+      n = ((long*) neighbours -> data)[k3+i];
+
+      //Initialise hvbar with values from hv
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+
+      if (n >= 0) {
+	hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
+
+	hmin[k] = min(hmin[k], hn);
+	hmax[k] = max(hmax[k], hn);
+      }
+    }
+  }
+
+  // Call underlying standard routine
+  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
+
+  // // //Py_DECREF(domain); //FIXME: NEcessary?
+  free(hmin);
+  free(hmax);
+
+  //return result using PyArray to avoid memory leak
+  return PyArray_Return(hvbar);
+  //return Py_BuildValue("");
+}
+
+PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
+  //a faster version of h_limiter above
+  PyObject *domain, *Tmp;
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
+    *hvbar,   //Limited depth at vertices (return values)
+    *neighbours;
+
+  int k, i, N, k3,k0,k1,k2;
+  int dimensions[2];
+  double beta_h; //Safety factor (see config.py)
+  double hmin, hmax, dh[3];
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
+    return NULL;
+  neighbours = get_consecutive_array(domain, "neighbours");
+
+  //Get safety factor beta_h
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
+
+  N = hc -> dimensions[0];
+
+  //Create hvbar
+  dimensions[0] = N;
+  dimensions[1] = 3;
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  for (k=0;k<N;k++){
+    k3=k*3;
+    //get the ids of the neighbours
+    k0 = ((long*) neighbours -> data)[k3];
+    k1 = ((long*) neighbours -> data)[k3+1];
+    k2 = ((long*) neighbours -> data)[k3+2];
+    //set hvbar provisionally
+    for (i=0;i<3;i++){
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
+    }
+    hmin=((double*) hc -> data)[k];
+    hmax=hmin;
+    if (k0>=0){
+      hmin=min(hmin,((double*) hc -> data)[k0]);
+      hmax=max(hmax,((double*) hc -> data)[k0]);
+    }
+    if (k1>=0){
+      hmin=min(hmin,((double*) hc -> data)[k1]);
+      hmax=max(hmax,((double*) hc -> data)[k1]);
+    }
+    if (k2>=0){
+      hmin=min(hmin,((double*) hc -> data)[k2]);
+      hmax=max(hmax,((double*) hc -> data)[k2]);
+    }
+    hmin-=((double*) hc -> data)[k];
+    hmax-=((double*) hc -> data)[k];
+    limit_gradient(dh,hmin,hmax,beta_h);
+    for (i=0;i<3;i++)
+      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
+  }
+  return PyArray_Return(hvbar);
+}
+
+PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
+  //
+  //      assign_windfield_values(xmom_update, ymom_update,
+  //                              s_vec, phi_vec, self.const)
+
+
+
+  PyArrayObject   //(one element per triangle)
+  *s_vec,         //Speeds
+  *phi_vec,       //Bearings
+  *xmom_update,   //Momentum updates
+  *ymom_update;
+
+
+  int N;
+  double cw;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOd",
+			&xmom_update,
+			&ymom_update,
+			&s_vec, &phi_vec,
+			&cw))
+    return NULL;
+
+  N = xmom_update -> dimensions[0];
+
+  _assign_wind_field_values(N,
+	   (double*) xmom_update -> data,
+	   (double*) ymom_update -> data,
+	   (double*) s_vec -> data,
+	   (double*) phi_vec -> data,
+	   cw);
+
+  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"},
+  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {"gravity", gravity, METH_VARARGS, "Print out"},
+  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
+  {"balance_deep_and_shallow", balance_deep_and_shallow,
+   METH_VARARGS, "Print out"},
+  {"h_limiter", h_limiter,
+   METH_VARARGS, "Print out"},
+  {"h_limiter_sw", h_limiter_sw,
+   METH_VARARGS, "Print out"},
+  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"assign_windfield_values", assign_windfield_values,
+   METH_VARARGS | METH_KEYWORDS, "Print out"},
+  //{"distribute_to_vertices_and_edges",
+  // distribute_to_vertices_and_edges, METH_VARARGS},
+  //{"update_conserved_quantities",
+  // update_conserved_quantities, METH_VARARGS},
+  //{"set_initialcondition",
+  // set_initialcondition, METH_VARARGS},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initshallow_water_kinetic_ext(void){
+  Py_InitModule("shallow_water_kinetic_ext", MethodTable);
+
+  import_array();     //Necessary for handling of NumPY structures
+}
