Changeset 4769

Timestamp:

Oct 30, 2007, 5:13:17 PM (17 years ago)

Author:

ole

Message:

Retired all Python code which is superseded by C implementations.
This has made the code base much smaller and there is less confusion when people try to understand the algorithms.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• abstract_2d_finite_volumes/quantity_ext.c (modified) (4 diffs)
• shallow_water/shallow_water_domain.py (modified) (39 diffs)
• shallow_water/shallow_water_ext.c (modified) (13 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (26 diffs)
• utilities/polygon.py (modified) (6 diffs)
• utilities/sparse.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -332 +332 @@
 
 
-        #General input checks
+        # General input checks
         L = [numeric, quantity, function, geospatial_data, points, filename]
         msg = 'Exactly one of the arguments '+\
@@ -360 +360 @@
 
 
-        #Determine which 'set_values_from_...' to use
+        # Determine which 'set_values_from_...' to use
 
         if numeric is not None:
@@ -1337 +1337 @@
     
 
-def update(quantity, timestep):
-    """Update centroid values based on values stored in
-    explicit_update and semi_implicit_update as well as given timestep
-
-    Function implementing forcing terms must take on argument
-    which is the domain and they must update either explicit
-    or implicit updates, e,g,:
-
-    def gravity(domain):
-        ....
-        domain.quantities['xmomentum'].explicit_update = ...
-        domain.quantities['ymomentum'].explicit_update = ...
-
-
-
-    Explicit terms must have the form
-
-        G(q, t)
-
-    and explicit scheme is
-
-       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
-
-
-    Semi implicit forcing terms are assumed to have the form
-
-       G(q, t) = H(q, t) q
-
-    and the semi implicit scheme will then be
-
-      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
-
-
-    """
-
-    from Numeric import sum, equal, ones, exp, Float
-
-    N = quantity.centroid_values.shape[0]
-
-
-    # Divide H by conserved quantity to obtain G (see docstring above)
-
-
-    for k in range(N):
-        x = quantity.centroid_values[k]
-        if x == 0.0:
-            #FIXME: Is this right
-            quantity.semi_implicit_update[k] = 0.0
-        else:
-            quantity.semi_implicit_update[k] /= x
-
-
-    # Semi implicit updates
-    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
-
-    if sum(less(denominator, 1.0)) > 0.0:
-        msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
-        raise msg
-
-    if sum(equal(denominator, 0.0)) > 0.0:
-        msg = 'Zero division in semi implicit update. Call Stephen :-)'
-        raise msg
-    else:
-        #Update conserved_quantities from semi implicit updates
-        quantity.centroid_values /= denominator
-
-#    quantity.centroid_values = exp(timestep*quantity.semi_implicit_update)*quantity.centroid_values
-
-    #Explicit updates
-    quantity.centroid_values += timestep*quantity.explicit_update
-
-def interpolate_from_vertices_to_edges(quantity):
-    """Compute edge values from vertex values using linear interpolation
-    """
-
-    for k in range(quantity.vertex_values.shape[0]):
-        q0 = quantity.vertex_values[k, 0]
-        q1 = quantity.vertex_values[k, 1]
-        q2 = quantity.vertex_values[k, 2]
-
-        quantity.edge_values[k, 0] = 0.5*(q1+q2)
-        quantity.edge_values[k, 1] = 0.5*(q0+q2)
-        quantity.edge_values[k, 2] = 0.5*(q0+q1)
-
-
-
-def backup_centroid_values(quantity):
-    """Copy centroid values to backup array"""
-
-    qc = quantity.centroid_values
-    qb = quantity.centroid_backup_values
-
-    #Check each triangle
-    for k in range(len(quantity.domain)):
-        qb[k] = qc[k]
-
-
-def saxpy_centroid_values(quantity,a,b):
-    """saxpy operation between centroid value and backup"""
-
-    qc = quantity.centroid_values
-    qb = quantity.centroid_backup_values
-
-    #Check each triangle
-    for k in range(len(quantity.domain)):
-        qc[k] = a*qc[k]+b*qb[k]       
-
-
-def compute_gradients(quantity):
-    """Compute gradients of triangle surfaces defined by centroids of
-    neighbouring volumes.
-    If one edge is on the boundary, use own centroid as neighbour centroid.
-    If two or more are on the boundary, fall back to first order scheme.
-    """
-
-    from Numeric import zeros, Float
-    from utilitites.numerical_tools import gradient
-
-    centroid_coordinates = quantity.domain.centroid_coordinates
-    surrogate_neighbours = quantity.domain.surrogate_neighbours
-    centroid_values = quantity.centroid_values
-    number_of_boundaries = quantity.domain.number_of_boundaries
-
-    N = centroid_values.shape[0]
-
-    a = zeros(N, Float)
-    b = zeros(N, Float)
-
-    for k in range(N):
-        if number_of_boundaries[k] < 2:
-            #Two or three true neighbours
-
-            #Get indices of neighbours (or self when used as surrogate)
-            k0, k1, k2 = surrogate_neighbours[k,:]
-
-            #Get data
-            q0 = centroid_values[k0]
-            q1 = centroid_values[k1]
-            q2 = centroid_values[k2]
-
-            x0, y0 = centroid_coordinates[k0] #V0 centroid
-            x1, y1 = centroid_coordinates[k1] #V1 centroid
-            x2, y2 = centroid_coordinates[k2] #V2 centroid
-
-            #Gradient
-            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-
-        elif number_of_boundaries[k] == 2:
-            #One true neighbour
-
-            #Get index of the one neighbour
-            for k0 in surrogate_neighbours[k,:]:
-                if k0 != k: break
-            assert k0 != k
-
-            k1 = k  #self
-
-            #Get data
-            q0 = centroid_values[k0]
-            q1 = centroid_values[k1]
-
-            x0, y0 = centroid_coordinates[k0] #V0 centroid
-            x1, y1 = centroid_coordinates[k1] #V1 centroid
-
-            #Gradient
-            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
-        else:
-            #No true neighbours -
-            #Fall back to first order scheme
-            pass
-
-
-    return a, b
-
-
-
-def limit(quantity):
-    """Limit slopes for each volume to eliminate artificial variance
-    introduced by e.g. second order extrapolator
-
-    This is an unsophisticated limiter as it does not take into
-    account dependencies among quantities.
-
-    precondition:
-    vertex values are estimated from gradient
-    postcondition:
-    vertex values are updated
-    """
-
-    from Numeric import zeros, Float
-
-    N = quantity.domain.number_of_nodes
-
-    beta_w = quantity.domain.beta_w
-
-    qc = quantity.centroid_values
-    qv = quantity.vertex_values
-
-    #Find min and max of this and neighbour's centroid values
-    qmax = zeros(qc.shape, Float)
-    qmin = zeros(qc.shape, Float)
-
-    for k in range(N):
-        qmax[k] = qc[k]
-        qmin[k] = qc[k]
-        for i in range(3):
-            n = quantity.domain.neighbours[k,i]
-            if n >= 0:
-                qn = qc[n] #Neighbour's centroid value
-
-                qmin[k] = min(qmin[k], qn)
-                qmax[k] = max(qmax[k], qn)
-        qmax[k] = min(qmax[k], 2.0*qc[k])
-        qmin[k] = max(qmin[k], 0.5*qc[k])
-
-
-    #Diffences between centroids and maxima/minima
-    dqmax = qmax - qc
-    dqmin = qmin - qc
-
-    #Deltas between vertex and centroid values
-    dq = zeros(qv.shape, Float)
-    for i in range(3):
-        dq[:,i] = qv[:,i] - qc
-
-    #Phi limiter
-    for k in range(N):
-
-        #Find the gradient limiter (phi) across vertices
-        phi = 1.0
-        for i in range(3):
-            r = 1.0
-            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
-            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
-
-            phi = min( min(r*beta_w, 1), phi )
-
-        #Then update using phi limiter
-        for i in range(3):
-            qv[k,i] = qc[k] + phi*dq[k,i]
-
-
 
 from anuga.utilities import compile
-if compile.can_use_C_extension('quantity_ext.c'):
-    #Replace python version with c implementations
+if compile.can_use_C_extension('quantity_ext.c'):    
+    # Underlying C implementations can be accessed 
 
     from quantity_ext import average_vertex_values, backup_centroid_values, \
@@ -1588 +1346 @@
 
     from quantity_ext import compute_gradients, limit,\
-    extrapolate_second_order, interpolate_from_vertices_to_edges, update
+    extrapolate_second_order, interpolate_from_vertices_to_edges, update    
+else:
+    msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    msg += 'Make sure compile_all.py has been run as described in '
+    msg += 'the ANUGA installation guide.'
+    raise Exception, msg
+
+
+

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -288 +288 @@
 // Gateways to Python
 PyObject *update(PyObject *self, PyObject *args) {
+  // FIXME (Ole): It would be great to turn this text into a Python DOC string
+
+  /*"""Update centroid values based on values stored in
+    explicit_update and semi_implicit_update as well as given timestep
+
+    Function implementing forcing terms must take on argument
+    which is the domain and they must update either explicit
+    or implicit updates, e,g,:
+
+    def gravity(domain):
+        ....
+        domain.quantities['xmomentum'].explicit_update = ...
+        domain.quantities['ymomentum'].explicit_update = ...
+
+
+
+    Explicit terms must have the form
+
+        G(q, t)
+
+    and explicit scheme is
+
+       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
+
+
+    Semi implicit forcing terms are assumed to have the form
+
+       G(q, t) = H(q, t) q
+
+    and the semi implicit scheme will then be
+
+      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
+
+  */
 
         PyObject *quantity;
@@ -397 +431 @@
 
 PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
-
+        //
+        //Compute edge values from vertex values using linear interpolation
+        // 
+        
         PyObject *quantity;
         PyArrayObject *vertex_values, *edge_values;
@@ -476 +513 @@
 
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
+  //"""Compute gradients of triangle surfaces defined by centroids of
+  //neighbouring volumes.
+  //If one edge is on the boundary, use own centroid as neighbour centroid.
+  //If two or more are on the boundary, fall back to first order scheme.
+  //"""
+
 
         PyObject *quantity, *domain, *R;
@@ -650 +693 @@
 
 PyObject *limit(PyObject *self, PyObject *args) {
+  //Limit slopes for each volume to eliminate artificial variance
+  //introduced by e.g. second order extrapolator
+
+  //This is an unsophisticated limiter as it does not take into
+  //account dependencies among quantities.
+
+  //precondition:
+  //  vertex values are estimated from gradient
+  //postcondition:
+  //  vertex values are updated
+  //
 
         PyObject *quantity, *domain, *Tmp;

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -74 +74 @@
 """
 
-#Subversion keywords:
+# Subversion keywords:
 #
-#$LastChangedDate$
-#$LastChangedRevision$
-#$LastChangedBy$
+# $LastChangedDate$
+# $LastChangedRevision$
+# $LastChangedBy$
 
 from Numeric import zeros, ones, Float, array, sum, size
@@ -104 +104 @@
 from anuga.config import optimise_dry_cells
 
-
-#Shallow water domain
+#---------------------
+# Shallow water domain
+#---------------------
 class Domain(Generic_Domain):
 
@@ -149 +150 @@
                                 number_of_full_triangles=number_of_full_triangles) 
 
-        #self.minimum_allowed_height = minimum_allowed_height
-        #self.H0 = minimum_allowed_height
+
         self.set_minimum_allowed_height(minimum_allowed_height)
         
@@ -167 +167 @@
         self.optimise_dry_cells = optimise_dry_cells
 
-        self.flux_function = flux_function_central
-        #self.flux_function = flux_function_kinetic
-        
-        self.forcing_terms.append(manning_friction)
+        self.forcing_terms.append(manning_friction_implicit)
         self.forcing_terms.append(gravity)
 
@@ -385 +382 @@
 
     def distribute_to_vertices_and_edges(self):
-        #Call correct module function
-        #(either from this module or C-extension)
+        # 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()
 
 
@@ -437 +396 @@
         """
 
-        #Call check integrity here rather than from user scripts
-        #self.check_integrity()
+        # Call check integrity here rather than from user scripts
+        # self.check_integrity()
 
         msg = 'Parameter beta_h must be in the interval [0, 1['
@@ -446 +405 @@
 
 
-        #Initial update of vertex and edge values before any STORAGE
-        #and or visualisation
-        #This is done again in the initialisation of the Generic_Domain
-        #evolve loop but we do it here to ensure the values are ok for storage
+        # Initial update of vertex and edge values before any STORAGE
+        # and or visualisation
+        # This is done again in the initialisation of the Generic_Domain
+        # evolve loop but we do it here to ensure the values are ok for storage
         self.distribute_to_vertices_and_edges()
 
         if self.store is True and self.time == 0.0:
             self.initialise_storage()
-            #print 'Storing results in ' + self.writer.filename
+            # 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
+            # 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
+        # Call basic machinery from parent class
         for t in Generic_Domain.evolve(self,
                                        yieldstep=yieldstep,
@@ -469 +428 @@
                                        skip_initial_step=skip_initial_step):
 
-            #Store model data, e.g. for subsequent visualisation
+            # 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
+            # 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):
@@ -486 +446 @@
         from anuga.shallow_water.data_manager import get_dataobject
 
-        #Initialise writer
+        # Initialise writer
         self.writer = get_dataobject(self, mode = 'w')
 
-        #Store vertices and connectivity
+        # Store vertices and connectivity
         self.writer.store_connectivity()
 
@@ -502 +462 @@
 
 
-#=============== 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
-    """
-
-    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
-
-
-
-####################################
+#=============== End of class Shallow Water Domain ===============================
+
+
+#-----------------
 # Flux computation
-# FIXME (Ole): Delete Python versions of code that is
-# always used in C-extensions - they get out of sync anyway
-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 and qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from math import sqrt
-
-    #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 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
-
-
+#-----------------
 
 def compute_fluxes(domain):
@@ -756 +487 @@
       domain.explicit_update is reset to computed flux values
       domain.timestep is set to the largest step satisfying all volumes.
+    
+
+    This wrapper calls the underlying C version of compute fluxes
     """
 
@@ -762 +496 @@
     N = len(domain) # number_of_triangles
 
-    #Shortcuts
+    # Shortcuts
     Stage = domain.quantities['stage']
     Xmom = domain.quantities['xmomentum']
@@ -768 +502 @@
     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 = domain.flux_function(normal, ql, qr, zl, zr)
-            flux -= edgeflux * domain.edgelengths[k,i]
-
-            #Update optimal_timestep on full cells
-            if  domain.tri_full_flag[k] == 1:
-                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.max_speed[k] = max_speed
-
-
-    domain.flux_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):
+    from shallow_water_ext import\
+         compute_fluxes_ext_central as compute_fluxes_ext
+
+
+    flux_timestep = compute_fluxes_ext(timestep,
+                                       domain.epsilon,
+                                       domain.H0,
+                                       domain.g,
+                                       domain.neighbours,
+                                       domain.neighbour_edges,
+                                       domain.normals,
+                                       domain.edgelengths,
+                                       domain.radii,
+                                       domain.areas,
+                                       domain.tri_full_flag,
+                                       Stage.edge_values,
+                                       Xmom.edge_values,
+                                       Ymom.edge_values,
+                                       Bed.edge_values,
+                                       Stage.boundary_values,
+                                       Xmom.boundary_values,
+                                       Ymom.boundary_values,
+                                       Stage.explicit_update,
+                                       Xmom.explicit_update,
+                                       Ymom.explicit_update,
+                                       domain.already_computed_flux,
+                                       domain.max_speed,
+                                       int(domain.optimise_dry_cells))
+
+    domain.flux_timestep = flux_timestep
+
+
+
+#---------------------------------------
+# Module functions for gradient limiting
+#---------------------------------------
+
+
+# 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(domain):
     """Wrapper calling C version of extrapolate_second_order_sw
     """
@@ -843 +556 @@
     Elevation = domain.quantities['elevation']
 
-    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,
-                                Elevation.centroid_values,
-                                domain.vertex_coordinates,
-                                Stage.vertex_values,
-                                Xmom.vertex_values,
-                                Ymom.vertex_values,
-                                Elevation.vertex_values,
-                                int(domain.optimise_dry_cells))
-
-def compute_fluxes_c(domain):
-    """Wrapper calling C version of compute fluxes
-    """
-
-    import sys
-
-    N = len(domain) # number_of_triangles
-
-    #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_ext_central as compute_fluxes_ext
-
-    domain.flux_timestep = compute_fluxes_ext(timestep, domain.epsilon,
-                                         domain.H0,
-                                         domain.g,
-                                         domain.neighbours,
-                                         domain.neighbour_edges,
-                                         domain.normals,
-                                         domain.edgelengths,
-                                         domain.radii,
-                                         domain.areas,
-                                         domain.tri_full_flag,
-                                         Stage.edge_values,
-                                         Xmom.edge_values,
-                                         Ymom.edge_values,
-                                         Bed.edge_values,
-                                         Stage.boundary_values,
-                                         Xmom.boundary_values,
-                                         Ymom.boundary_values,
-                                         Stage.explicit_update,
-                                         Xmom.explicit_update,
-                                         Ymom.explicit_update,
-                                         domain.already_computed_flux,
-                                         domain.max_speed,
-                                         int(domain.optimise_dry_cells))
-
-
-####################################
-# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
+    extrapol2(domain,
+              domain.surrogate_neighbours,
+              domain.number_of_boundaries,
+              domain.centroid_coordinates,
+              Stage.centroid_values,
+              Xmom.centroid_values,
+              Ymom.centroid_values,
+              Elevation.centroid_values,
+              domain.vertex_coordinates,
+              Stage.vertex_values,
+              Xmom.vertex_values,
+              Ymom.vertex_values,
+              Elevation.vertex_values,
+              int(domain.optimise_dry_cells))
+
 
 def distribute_to_vertices_and_edges(domain):
@@ -930 +598 @@
     from anuga.config import optimised_gradient_limiter
 
-    #Remove very thin layers of water
+    # Remove very thin layers of water
     protect_against_infinitesimal_and_negative_heights(domain)
 
-    #Extrapolate all conserved quantities
+    # Extrapolate all conserved quantities
     if optimised_gradient_limiter:
-        #MH090605 if second order,
-        #perform the extrapolation and limiting on
-        #all of the conserved quantities
+        # MH090605 if second order,
+        # perform the extrapolation and limiting on
+        # all of the conserved quantities
 
         if (domain._order_ == 1):
@@ -948 +616 @@
             raise 'Unknown order'
     else:
-        #old code:
+        # Old code:
         for name in domain.conserved_quantities:
             Q = domain.quantities[name]
@@ -986 +654 @@
     """
 
-    #Shortcuts
+    # 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
-    #FIXME: Modify according to c-version - or discard altogether.
-    for k in range(len(domain)):
-
-        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
@@ -1020 +664 @@
     protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
             domain.epsilon, wc, zc, xmomc, ymomc)
-
 
 
@@ -1031 +674 @@
     This limiter depends on two quantities (w,z) so it resides within
     this module rather than within quantity.py
-    """
-
-    N = len(domain)
+
+    Wrapper for c-extension
+    """
+
+    N = len(domain) # number_of_triangles
     beta_h = domain.beta_h
 
-    #Shortcuts
+    # Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
@@ -1043 +688 @@
     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
-    """
-
-    N = len(domain) # number_of_triangles
-    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)
+    from shallow_water_ext import h_limiter_sw
+    hvbar = h_limiter_sw(domain, hc, hv)
 
     return hvbar
@@ -1133 +708 @@
     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(len(domain)):
-        #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
+    Wrapper for C implementation
     """
 
@@ -1229 +717 @@
     # Omit updating xmomv 
     #
-    from shallow_water_ext import balance_deep_and_shallow
+    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
     
     # 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
@@ -1248 +734 @@
     ymomv = domain.quantities['ymomentum'].vertex_values
 
-    #Limit h
+    # Limit h
     if domain.beta_h > 0:
         hvbar = h_limiter(domain)
         
-        balance_deep_and_shallow(domain, domain.beta_h,
-                                 wc, zc, wv, zv, hvbar,
-                                 xmomc, ymomc, xmomv, ymomv)        
+        balance_deep_and_shallow_c(domain, domain.beta_h,
+                                   wc, zc, wv, zv, hvbar,
+                                   xmomc, ymomc, xmomv, ymomv)        
     else:
         # print 'Using first order h-limiter'
         # FIXME: Pass wc in for now - it will be ignored.
         
-        #This is how one would make a first order h_limited value
-        #as in the old balancer (pre 17 Feb 2005):
-        # If we wish to hard wire this, one should modify the C-code
-        #from Numeric import zeros, Float
-        #hvbar = zeros( (len(wc), 3), Float)
-        #for i in range(3):
-        #    hvbar[:,i] = wc[:] - zc[:]
-
-        balance_deep_and_shallow(domain, domain.beta_h,
-                                 wc, zc, wv, zv, wc, 
-                                 xmomc, ymomc, xmomv, ymomv)
-
-
-
-
-
-###############################################
-#Boundaries - specific to the shallow water wave equation
+        # This is how one would make a first order h_limited value
+        # as in the old balancer (pre 17 Feb 2005):
+        #  If we wish to hard wire this, one should modify the C-code
+        # from Numeric import zeros, Float
+        # hvbar = zeros( (len(wc), 3), Float)
+        # for i in range(3):
+        #     hvbar[:,i] = wc[:] - zc[:]
+
+        balance_deep_and_shallow_c(domain, domain.beta_h,
+                                   wc, zc, wv, zv, wc, 
+                                   xmomc, ymomc, xmomv, ymomv)
+
+
+
+
+#------------------------------------------------------------------
+# Boundary conditions - specific to the shallow water wave equation
+#------------------------------------------------------------------
 class Reflective_boundary(Boundary):
     """Reflective boundary returns same conserved quantities as
@@ -1293 +779 @@
             raise msg
 
-        #Handy shorthands
+        # Handy shorthands
         self.stage   = domain.quantities['stage'].edge_values
         self.xmom    = domain.quantities['xmomentum'].edge_values
@@ -1377 +863 @@
 
 
-        #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)
+        # 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)
 
 
@@ -1424 +910 @@
 
 
-        #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)
+        # 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)
 
 
@@ -1511 +997 @@
     
 
-#########################
-#Standard forcing terms:
-#
+#-----------------------
+# Standard forcing terms
+#-----------------------
+
 def gravity(domain):
     """Apply gravitational pull in the presence of bed slope
+    Wrapper calls underlying C implementation
     """
 
@@ -1521 +1009 @@
     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
+    stage = domain.quantities['stage']
+    elevation = domain.quantities['elevation']
+
+    h = stage.centroid_values - elevation.centroid_values
+    z = elevation.vertex_values
 
     x = domain.get_vertex_coordinates()
     g = domain.g
-
-    for k in range(len(domain)):
-        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.centroid_values - elevation.centroid_values
-    z = elevation.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-    
-
-    from shallow_water_ext import gravity
-    gravity(g, h, z, x, xmom, ymom) #, 1.0e-6)
-
-
-
-def manning_friction(domain):
-    """Apply (Manning) friction to water momentum
-    (Python version)
-    """
-
-    from math import sqrt
+    
+
+    from shallow_water_ext import gravity as gravity_c
+    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
+
+
+
+def manning_friction_implicit(domain):
+    """Apply (Manning) friction to water momentum    
+    Wrapper for c version
+    """
+
+
+    #print 'Implicit friction'
+
+    xmom = domain.quantities['xmomentum']
+    ymom = domain.quantities['ymomentum']
 
     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
+
+    uh = xmom.centroid_values
+    vh = ymom.centroid_values
     eta = domain.quantities['friction'].centroid_values
 
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+    xmom_update = xmom.semi_implicit_update
+    ymom_update = ymom.semi_implicit_update
 
     N = len(domain)
@@ -1587 +1049 @@
     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_implicit_c(domain):
-    """Wrapper for c version
-    """
-
-
-    #print 'Implicit friction'
+    from shallow_water_ext import manning_friction as manning_friction_c
+    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
+
+
+def manning_friction_explicit(domain):
+    """Apply (Manning) friction to water momentum    
+    Wrapper for c version
+    """
+
+    # print 'Explicit friction'
 
     xmom = domain.quantities['xmomentum']
@@ -1615 +1070 @@
     eta = domain.quantities['friction'].centroid_values
 
-    xmom_update = xmom.semi_implicit_update
-    ymom_update = ymom.semi_implicit_update
+    xmom_update = xmom.explicit_update
+    ymom_update = ymom.explicit_update
 
     N = len(domain)
@@ -1622 +1077 @@
     g = domain.g
 
-    from shallow_water_ext import manning_friction
-    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
-
-
-def manning_friction_explicit_c(domain):
-    """Wrapper for c version
-    """
-
-    #print 'Explicit friction'
-
-    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.explicit_update
-    ymom_update = ymom.explicit_update
-
-    N = len(domain)
-    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)
-
-
+    from shallow_water_ext import manning_friction as manning_friction_c
+    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
+
+
+# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
+# Is it still needed (30 Oct 2007)?
 def linear_friction(domain):
     """Apply linear friction to water momentum
@@ -1687 +1117 @@
 
 
+#---------------------------------
+# Experimental auxiliary functions
+#---------------------------------
 def check_forcefield(f):
     """Check that f is either
@@ -1703 +1136 @@
         except Exception, e:
             msg = 'Function %s could not be executed:\n%s' %(f, e)
-            #FIXME: Reconsider this semantics
+            # FIXME: Reconsider this semantics
             raise msg
 
@@ -1714 +1147 @@
             raise msg
 
-        #Is this really what we want?
+        # Is this really what we want?
         msg = 'Return vector from function %s ' %f
         msg += 'must have same lenght as input vectors'
@@ -1781 +1214 @@
             phi = args[1]
         elif len(args) == 1:
-            #Assume vector function returning (s, phi)(t,x,y)
+            # 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
+           # Assume info is in 2 keyword arguments
 
            if len(kwargs) == 2:
@@ -1817 +1250 @@
             s_vec = self.speed(t, xc[:,0], xc[:,1])
         else:
-            #Assume s is a scalar
+            # Assume s is a scalar
 
             try:
@@ -1830 +1263 @@
             phi_vec = self.phi(t, xc[:,0], xc[:,1])
         else:
-            #Assume phi is a scalar
+            # Assume phi is a scalar
 
             try:
@@ -1854 +1287 @@
         phi = phi_vec[k]
 
-        #Convert to radians
+        # Convert to radians
         phi = phi*pi/180
 
-        #Compute velocity vector (u, v)
+        # Compute velocity vector (u, v)
         u = s*cos(phi)
         v = s*sin(phi)
 
-        #Compute wind stress
+        # Compute wind stress
         S = const * sqrt(u**2 + v**2)
         xmom_update[k] += S*u
@@ -2041 +1474 @@
 #------------------
 
+
 from anuga.utilities import compile
 if compile.can_use_C_extension('shallow_water_ext.c'):
-    # Replace python version with c implementations
+    # Underlying C implementations can be accessed 
 
     from shallow_water_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_implicit_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)
-
+else:
+    msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    msg += 'Make sure compile_all.py has been run as described in '
+    msg += 'the ANUGA installation guide.'
+    raise Exception, msg
 
 

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -194 +194 @@
 
 // Innermost flux function (using stage w=z+h)
-int flux_function_central(double *q_left, double *q_right,
-                          double z_left, double z_right,
-                          double n1, double n2,
-                          double epsilon, double H0, double g,
-                          double *edgeflux, double *max_speed) {
+int _flux_function_central(double *q_left, double *q_right,
+                           double z_left, double z_right,
+                           double n1, double n2,
+                           double epsilon, double H0, double g,
+                           double *edgeflux, double *max_speed) {
 
   /*Compute fluxes between volumes for the shallow water wave equation
@@ -737 +737 @@
 ///////////////////////////////////////////////////////////////////
 // Gateways to Python
+
+
+
+PyObject *flux_function_central(PyObject *self, PyObject *args) {
+  //
+  // Gateway to innermost flux function.
+  // This is only used by the unit tests as the c implementation is
+  // normally called by compute_fluxes in this module.
+
+
+  PyArrayObject *normal, *ql, *qr,  *edgeflux;
+  double g, epsilon, max_speed, H0, zl, zr;
+
+  if (!PyArg_ParseTuple(args, "OOOddOddd",
+                        &normal, &ql, &qr, &zl, &zr, &edgeflux,
+                        &epsilon, &g, &H0)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
+    return NULL;
+  }
+
+  
+  _flux_function_central((double*) ql -> data, 
+                         (double*) qr -> data, 
+                         zl, 
+                         zr,                                             
+                         ((double*) normal -> data)[0],                                                                          
+                         ((double*) normal -> data)[1],                  
+                         epsilon, H0, g,
+                         (double*) edgeflux -> data, 
+                         &max_speed);
+  
+  return Py_BuildValue("d", max_speed);  
+}
+
+
 
 PyObject *gravity(PyObject *self, PyObject *args) {
@@ -1494 +1529 @@
   }  
 
-  //Input checks (convert sequences into numeric arrays)
+  // Input checks (convert sequences into numeric arrays)
   q = (PyArrayObject *)
     PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
@@ -1506 +1541 @@
   }
 
-  //Allocate space for return vector r (don't DECREF)
+  // Allocate space for return vector r (don't DECREF)
   dimensions[0] = 3;
   r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
 
-  //Copy
+  // Copy
   for (i=0; i<3; i++) {
     ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
   }
 
-  //Get normal and direction
+  // Get normal and direction
   n1 = ((double *) normal -> data)[0];
   n2 = ((double *) normal -> data)[1];
   if (direction == -1) n2 = -n2;
 
-  //Rotate
+  // Rotate
   _rotate((double *) r -> data, n1, n2);
 
-  //Release numeric arrays
+  // Release numeric arrays
   Py_DECREF(q);
   Py_DECREF(normal);
 
-  //return result using PyArray to avoid memory leak
+  // Return result using PyArray to avoid memory leak
   return PyArray_Return(r);
 }
@@ -1642 +1677 @@
 
       // Edge flux computation (triangle k, edge i)
-      flux_function_central(ql, qr, zl, zr,
-                            normals[ki2], normals[ki2+1],
-                            epsilon, H0, g,
-                            edgeflux, &max_speed);
+      _flux_function_central(ql, qr, zl, zr,
+                             normals[ki2], normals[ki2+1],
+                             epsilon, H0, g,
+                             edgeflux, &max_speed);
       
       
@@ -1898 +1933 @@
     *xmom_explicit_update,
     *ymom_explicit_update,
-    *already_computed_flux;//tracks whether the flux across an edge has already been computed
-
-
-  //Local variables
+    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
+
+
+  // Local variables
   double timestep, max_speed, epsilon, g, H0;
   double normal[2], ql[3], qr[3], zl, zr;
@@ -2060 +2095 @@
 
 PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
+  // 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.
   //
   //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
@@ -2160 +2204 @@
   int k, i, n, N, k3;
   int dimensions[2];
-  double beta_h; //Safety factor (see config.py)
+  double beta_h; // Safety factor (see config.py)
   double *hmin, *hmax, hn;
 
@@ -2201 +2245 @@
       n = ((long*) neighbours -> data)[k3+i];
 
-      //Initialise hvbar with values from hv
+      // 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
+        hn = ((double*) hc -> data)[n]; // Neighbour's centroid value
 
         hmin[k] = min(hmin[k], hn);
@@ -2216 +2260 @@
   _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
 
-  // // //Py_DECREF(domain); //FIXME: NEcessary?
+  // // //Py_DECREF(domain); //FIXME: Necessary?
   free(hmin);
   free(hmax);
 
-  //return result using PyArray to avoid memory leak
+  // Return result using PyArray to avoid memory leak
   return PyArray_Return(hvbar);
-  //return Py_BuildValue("");
 }
 
@@ -2337 +2380 @@
 
 
-//////////////////////////////////////////
+//-------------------------------
 // Method table for python module
+//-------------------------------
 static struct PyMethodDef MethodTable[] = {
   /* The cast of the function is necessary since PyCFunction values
@@ -2351 +2395 @@
   {"gravity", gravity, METH_VARARGS, "Print out"},
   {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
+  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
   {"balance_deep_and_shallow", balance_deep_and_shallow,
    METH_VARARGS, "Print out"},
@@ -2360 +2405 @@
   {"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}
 };

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -10 +10 @@
 
 from shallow_water_domain import *
-from shallow_water_domain import flux_function_central as flux_function
+
+# Get gateway to C implementation of flux function for direct testing
+from shallow_water_ext import flux_function_central as flux_function
 
 class Weir:
@@ -252 +254 @@
 
 
-    #FIXME (Ole): Individual flux tests do NOT test C implementation directly.    
+    # Individual flux tests
     def test_flux_zero_case(self):
         ql = zeros( 3, Float )
         qr = zeros( 3, Float )
         normal = zeros( 2, Float )
+        edgeflux = zeros( 3, Float )
         zl = zr = 0.
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [0,0,0])
+        H0 = 0.0
+        
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
+
+        assert allclose(edgeflux, [0,0,0])
         assert max_speed == 0.
 
@@ -269 +274 @@
         ql = array([w, 0, 0])
         qr = array([w, 0, 0])
+        edgeflux = zeros(3, Float)        
         zl = zr = 0.
         h = w - (zl+zr)/2
-
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [0., 0.5*g*h**2, 0.])
+        H0 = 0.0
+
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
+        assert allclose(edgeflux, [0., 0.5*g*h**2, 0.])
         assert max_speed == sqrt(g*h)
 
@@ -299 +305 @@
         qr = array([-0.2, 2, 3])
         zl = zr = -0.5
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [2.,13.77433333, 20.])
+        edgeflux = zeros(3, Float)                
+
+        H0 = 0.0
+
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
+
+        assert allclose(edgeflux, [2.,13.77433333, 20.])
         assert allclose(max_speed, 8.38130948661)
 
@@ -311 +321 @@
         qr = array([-0.075, 2, 3])
         zl = zr = -0.375
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [-3.,-20.0, -30.441])
+
+        edgeflux = zeros(3, Float)                
+        H0 = 0.0
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
+
+        assert allclose(edgeflux, [-3.,-20.0, -30.441])
         assert allclose(max_speed, 11.7146428199)
 
@@ -322 +335 @@
         qr = array([-0.075, 2, 3])
         zl = zr = -0.375
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
+
+        edgeflux = zeros(3, Float)                
+        H0 = 0.0
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
+
+        assert allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
         assert allclose(max_speed, 4.0716654239)
 
@@ -333 +349 @@
         qr = array([-0.30683287, 0.1071986, 0.05930515])
         zl = zr = -0.375
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
+
+        edgeflux = zeros(3, Float)                
+        H0 = 0.0
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                
+
+        assert allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
         assert allclose(max_speed, 1.31414103233)
 
@@ -532 +551 @@
 
     def test_compute_fluxes0(self):
-        #Do a full triangle and check that fluxes cancel out for
-        #the constant stage case
+        # Do a full triangle and check that fluxes cancel out for
+        # the constant stage case
 
         a = [0.0, 0.0]
@@ -554 +573 @@
         assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
 
-        zl=zr=0. #Assume flat bed
-
-        #Flux across right edge of volume 1
+        zl=zr=0. # Assume flat bed
+
+        edgeflux = zeros(3, Float)        
+        edgeflux0 = zeros(3, Float)
+        edgeflux1 = zeros(3, Float)
+        edgeflux2 = zeros(3, Float)                                
+        H0 = 0.0        
+
+        # Flux across right edge of volume 1
         normal = domain.get_normal(1,0)
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
-        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        #Check that flux seen from other triangles is inverse
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)                        
+
+        # Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(2,2)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-        assert allclose(flux + flux0, 0.)
-
-        #Flux across upper edge of volume 1
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                
+
+        assert allclose(edgeflux0 + edgeflux, 0.)
+
+        # Flux across upper edge of volume 1
         normal = domain.get_normal(1,1)
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
-        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        #Check that flux seen from other triangles is inverse
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                                        
+
+        # Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(3,0)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-        assert allclose(flux + flux1, 0.)
-
-        #Flux across lower left hypotenuse of volume 1
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                               
+
+        assert allclose(edgeflux1 + edgeflux, 0.)        
+        
+
+        # Flux across lower left hypotenuse of volume 1
         normal = domain.get_normal(1,2)
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
-        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        #Check that flux seen from other triangles is inverse
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)                                                               
+
+        # Check that flux seen from other triangles is inverse
         tmp = qr; qr=ql; ql=tmp
         normal = domain.get_normal(0,1)
-        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
-        assert allclose(flux + flux2, 0.)
-
-
-        #Scale by edgelengths, add up anc check that total flux is zero
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                                       
+        assert allclose(edgeflux2 + edgeflux, 0.)
+
+
+        # Scale by edgelengths, add up anc check that total flux is zero
         e0 = domain.edgelengths[1, 0]
         e1 = domain.edgelengths[1, 1]
         e2 = domain.edgelengths[1, 2]
 
-        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
-
-        #Now check that compute_flux yields zeros as well
+        assert allclose(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2, 0.)
+
+        # Now check that compute_flux yields zeros as well
         domain.compute_fluxes()
 
@@ -635 +663 @@
         zl=zr=0. #Assume flat bed
 
-        #Flux across right edge of volume 1
+        edgeflux = zeros(3, Float)        
+        edgeflux0 = zeros(3, Float)
+        edgeflux1 = zeros(3, Float)
+        edgeflux2 = zeros(3, Float)                                
+        H0 = 0.0        
+        
+
+        # Flux across right edge of volume 1
         normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
         assert allclose(normal, [1, 0])
@@ -644 +679 @@
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
         assert allclose(qr, [val2, 0, 0])
-        
-        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        #Flux across edge in the east direction (as per normal vector)
-        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
+
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)                                                       
+
+        # Flux across edge in the east direction (as per normal vector)
+        assert allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
         assert allclose(max_speed, 9.21592824046)
 
@@ -655 +690 @@
         normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
         assert allclose(normal_opposite, [-1, 0])
-        flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
-        assert allclose(flux_opposite, [-15.3598804, -253.71111111, 0.])
+
+        max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                             
+        #flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
+        assert allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
         
 
@@ -663 +700 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
-        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
-        assert allclose(flux1, [2.4098563, 0., 123.04444444])
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                                                               
+
+        assert allclose(edgeflux1, [2.4098563, 0., 123.04444444])
         assert allclose(max_speed, 7.22956891292)
 
@@ -671 +709 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
-        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)        
+
+        assert allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
         assert allclose(max_speed, 7.22956891292)
 
@@ -681 +719 @@
         e2 = domain.edgelengths[1, 2]
 
-        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+        total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
         assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
 
@@ -733 +771 @@
 
         zl=zr=0 #Assume flat zero bed
+        edgeflux = zeros(3, Float)        
+        edgeflux0 = zeros(3, Float)
+        edgeflux1 = zeros(3, Float)
+        edgeflux2 = zeros(3, Float)                                
+        H0 = 0.0        
+        
 
         domain.set_quantity('elevation', zl*ones( (4,3) ))
@@ -759 +803 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
-        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)                
 
         #Flux across upper edge of volume 1
@@ -765 +809 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
-        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                        
 
         #Flux across lower left hypotenuse of volume 1
@@ -771 +815 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
-        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)                
 
         #Scale, add up and check that compute_fluxes is correct for vol 1
@@ -778 +822 @@
         e2 = domain.edgelengths[1, 2]
 
-        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+        total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
 
 
@@ -813 +857 @@
 
 
+        edgeflux = zeros(3, Float)        
+        edgeflux0 = zeros(3, Float)
+        edgeflux1 = zeros(3, Float)
+        edgeflux2 = zeros(3, Float)                                
+        H0 = 0.0        
+        
+
+
         domain.set_quantity('stage', [[val0, val0-1, val0-2],
                                       [val1, val1+1, val1],
@@ -835 +887 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
-        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
 
         #Flux across upper edge of volume 1
@@ -841 +893 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
-        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)        
 
         #Flux across lower left hypotenuse of volume 1
@@ -847 +899 @@
         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
-        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)        
 
         #Scale, add up and check that compute_fluxes is correct for vol 1
@@ -854 +906 @@
         e2 = domain.edgelengths[1, 2]
 
-        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
+        total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
 
         domain.compute_fluxes()
@@ -5136 +5188 @@
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Shallow_Water,'test')    
+    suite = unittest.makeSuite(Test_Shallow_Water,'test')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')

anuga_core/source/anuga/utilities/polygon.py

@@ -254 +254 @@
     http://www.alienryderflex.com/polygon/
 
-    """
+    Uses underlying C-implementation in polygon_ext.c
+    """
+
+
+    if verbose: print 'Checking input to separate_points_by_polygon'
+
 
     #Input checks
+
+    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
+    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
 
 
@@ -267 +275 @@
         raise msg
 
+    #if verbose: print 'Checking input to separate_points_by_polygon 2'
     try:
         polygon = ensure_numeric(polygon, Float)
@@ -275 +284 @@
         raise msg
 
+    #if verbose: print 'check'
+
     assert len(polygon.shape) == 2,\
        'Polygon array must be a 2d array of vertices'
@@ -290 +301 @@
     M = points.shape[0]  #Number of points
 
-    px = polygon[:,0]
-    py = polygon[:,1]
-
-    #Used for an optimisation when points are far away from polygon
-    minpx = min(px); maxpx = max(px)
-    minpy = min(py); maxpy = max(py)
-
-
-    #Begin main loop
-    indices = zeros(M, Int)
-
-    inside_index = 0    #Keep track of points inside
-    outside_index = M-1 #Keep track of points outside (starting from end)
-
-    if verbose: print 'Separating %d points' %M        
-    for k in range(M):
-
-        if verbose:
-            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
-
-        #for each point
-        x = points[k, 0]
-        y = points[k, 1]
-
-        inside = False
-
-        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
-            #Check polygon
-            for i in range(N):
-                j = (i+1)%N
-
-                #Check for case where point is contained in line segment
-                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
-                    if closed:
-                        inside = True
-                    else:
-                        inside = False
-                    break
-                else:
-                    #Check if truly inside polygon
-                    if py[i] < y and py[j] >= y or\
-                       py[j] < y and py[i] >= y:
-                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
-                            inside = not inside
-
-        if inside:
-            indices[inside_index] = k
-            inside_index += 1
-        else:
-            indices[outside_index] = k
-            outside_index -= 1
-
-    return indices, inside_index
-
-
-def separate_points_by_polygon_c(points, polygon,
-                                 closed = True, verbose = False):
-    """Determine whether points are inside or outside a polygon
-
-    C-wrapper
-    """
-
-
-    if verbose: print 'Checking input to separate_points_by_polygon'
-
-
-    #Input checks
-
-    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
-    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
-
-
-    try:
-        points = ensure_numeric(points, Float)
-    except NameError, e:
-        raise NameError, e
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    #if verbose: print 'Checking input to separate_points_by_polygon 2'
-    try:
-        polygon = ensure_numeric(polygon, Float)
-    except NameError, e:
-        raise NameError, e
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-    #if verbose: print 'check'
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices'
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns'
-
-    assert len(points.shape) == 2,\
-       'Points array must be a 2d array'
-
-    assert points.shape[1] == 2,\
-       'Points array must have two columns'
-
-    N = polygon.shape[0] #Number of vertices in polygon
-    M = points.shape[0]  #Number of points
-
-    from polygon_ext import separate_points_by_polygon
+
 
     if verbose: print 'Allocating array for indices'
@@ -404 +309 @@
     #if verbose: print 'Calling C-version of inside poly'
 
-    count = separate_points_by_polygon(points, polygon, indices,
-                                       int(closed), int(verbose))
+    from polygon_ext import separate_points_by_polygon as sep_points
+    count = sep_points(points, polygon, indices,
+                       int(closed), int(verbose))
 
     if verbose: print 'Found %d points (out of %d) inside polygon'\
@@ -792 +698 @@
 from anuga.utilities.compile import can_use_C_extension
 if can_use_C_extension('polygon_ext.c'):
+    # Underlying C implementations can be accessed 
+
     from polygon_ext import point_on_line
-    separate_points_by_polygon = separate_points_by_polygon_c
+else:
+    msg = 'C implementations could not be accessed by %s.\n ' %__file__
+    msg += 'Make sure compile_all.py has been run as described in '
+    msg += 'the ANUGA installation guide.'
+    raise Exception, msg
 
 

anuga_core/source/anuga/utilities/sparse.py

@@ -320 +320 @@
 
 
-
-def csr_mv(self, B):
-    """Python version of sparse (CSR) matrix multiplication
-    """
-
-    from Numeric import zeros, Float
-
-
-    #Assume numeric types from now on
-        
-    if len(B.shape) == 0:
-        #Scalar - use __rmul__ method
-        R = B*self
-        
-    elif len(B.shape) == 1:
-        #Vector
-        assert B.shape[0] == self.N, 'Mismatching dimensions'
-        
-        R = zeros(self.M, Float) #Result
-        
-        #Multiply nonzero elements
-        for i in range(self.M):
-            for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
-                j = self.colind[ckey]
-                R[i] += self.data[ckey]*B[j]            
-                
-    elif len(B.shape) == 2:
-        
-        R = zeros((self.M, B.shape[1]), Float) #Result matrix
-        
-        #Multiply nonzero elements
-        
-        for col in range(R.shape[1]):
-            #For each column
-            for i in range(self.M):
-                for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
-                    j = self.colind[ckey]
-                    R[i, col] += self.data[ckey]*B[j,col]            
-                    
-    else:
-        raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
-    
-    return R
-
-
-
-#Setup for C extensions
+# Setup for C extensions
 from anuga.utilities import compile
 if compile.can_use_C_extension('sparse_ext.c'):
-    #Replace python version with c implementation
+    # Access underlying c implementations
     from sparse_ext import csr_mv
 
+
 if __name__ == '__main__':
-
+    # A little selftest
+    
     from Numeric import allclose, array, Float