Index: anuga_work/development/flow_1d/sww_flow/__init__.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/__init__.py	(revision 7830)
+++ anuga_work/development/flow_1d/sww_flow/__init__.py	(revision 7830)
@@ -0,0 +1,9 @@
+#! /usr/bin/python
+
+# To change this template, choose Tools | Templates
+# and open the template in the editor.
+
+__author__="Stephen Roberts"
+__date__ ="$14/06/2010 7:50:20 PM$"
+
+pass
Index: anuga_work/development/flow_1d/sww_flow/compile_all.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/compile_all.py	(revision 7827)
+++ 	(revision )
@@ -1,15 +1,0 @@
-import os
-
-buildroot = os.getcwd()
-
-#Attempt to compile all extensions
-
-anuga_dir = '..' + os.sep + '..' + os.sep + '..'
-utilities_dir = anuga_dir + os.sep + 'trunk' + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities'
-
-execfile( utilities_dir + os.sep + 'compile.py')
-
-
-os.chdir(buildroot)    
-
-    
Index: anuga_work/development/flow_1d/sww_flow/config.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/config.py	(revision 7827)
+++ 	(revision )
@@ -1,112 +1,0 @@
-"""Module where global model parameters are set for anuga_1d
-"""
-
-epsilon = 1.0e-12
-h0 = 1.0e-12
-
-default_boundary_tag = 'exterior'
-
-
-time_format = '%d/%m/%y %H:%M:%S'
-
-min_timestep = 1.0e-6 #Should be computed based on geometry
-max_timestep = 1.0e3
-max_smallsteps = 50  # Max number of degenerate steps allowed b4 trying first order
-#This is how:
-#Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)
-#Then work out minimal internal distance in mesh r_min and set
-#min_timestep = r_min/v_max
-#
-#Max speeds are calculated in the flux function as
-#
-#lambda = v +/- sqrt(gh)
-#
-# so with 500 m/s, h ~ 500^2/g = 2500 m well out of the domain of the
-# shallow water wave equation
-#
-#The actual soundspeed can be as high as 1530m/s
-#(see http://staff.washington.edu/aganse/public.projects/clustering/clustering.html),
-#but that would only happen with h>225000m in this equation. Why ?
-#The maximal speed we specify is really related to the max speed
-#of surface pertubation
-#
-
-
-v_max = 100 #For use in domain_ext.c
-sound_speed = 500
-
-
-max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order
-
-manning = 0.0  #Manning's friction coefficient
-g = 9.8       #Gravity
-#g(phi) = 9780313 * (1 + 0.0053024 sin(phi)**2 - 0.000 0059 sin(2*phi)**2) micro m/s**2, where phi is the latitude
-#The 'official' average is 9.80665
-
-
-
-
-eta_w = 3.0e-3 #Wind stress coefficient
-rho_a = 1.2e-3 #Atmospheric density
-rho_w = 1023   #Fluid density [kg/m^3] (rho_w = 1023 for salt water)
-
-
-#Betas [0;1] control the allowed steepness of gradient for second order
-#extrapolations. Values of 1 allow the steepes gradients while
-#lower values are more conservative. Values of 0 correspond to
-#1'st order extrapolations.
-#
-# Large values of beta_h may cause simulations to require more timesteps
-# as surface will 'hug' closer to the bed.
-# Small values of beta_h will make code faster, but one may experience
-# artificial momenta caused by discontinuities in water depths in
-# the presence of steep slopes. One example of this would be
-# stationary water 'lapping' upwards to a higher point on the coast.
-#
-#
-#
-#There are separate betas for the w-limiter and the h-limiter
-#
-#
-#
-#
-#Good values are:
-#beta_w = 0.9
-#beta_h = 0.2
-
-
-
-beta_w = 1.5
-beta_h = 0.2
-timestepping_method = 'euler'
-
-CFL = 1.0  #FIXME (ole): Is this in use yet??
-           #(Steve) yes, change domain.CFL to
-           #make changes
-
-
-pmesh_filename = '.\\pmesh'
-
-
-import os, sys
-
-if sys.platform == 'win32':
-    #default_datadir = 'C:\grohm_output'
-    default_datadir = '.'
-else:
-    #default_datadir = os.path.expanduser('~'+os.sep+'grohm_output')
-    default_datadir = '.'
-
-
-use_extensions = True    #Try to use C-extensions
-#use_extensions = False   #Do not use C-extensions
-
-use_psyco = True  #Use psyco optimisations
-#use_psyco = False  #Do not use psyco optimisations
-
-
-optimised_gradient_limiter = True #Use hardwired gradient limiter
-
-#Specific to shallow water W.E.
-h_min = minimum_allowed_height = 1.0e-6 #1.0e-6 #Water depth below which it is considered to be 0
-maximum_allowed_speed = 1000.0
Index: anuga_work/development/flow_1d/sww_flow/sww_boundary_conditions.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/sww_boundary_conditions.py	(revision 7827)
+++ anuga_work/development/flow_1d/sww_flow/sww_boundary_conditions.py	(revision 7830)
@@ -7,5 +7,5 @@
 __date__ ="$05/06/2010 5:44:05 PM$"
 
-from generic_domain import *
+from flow_1d.generic_1d.generic_domain import *
 
 class Dirichlet_boundary(Boundary):
Index: anuga_work/development/flow_1d/sww_flow/sww_domain.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/sww_domain.py	(revision 7827)
+++ anuga_work/development/flow_1d/sww_flow/sww_domain.py	(revision 7830)
@@ -45,5 +45,5 @@
 import numpy
 
-from generic_domain import *
+from flow_1d.generic_1d.generic_domain import Generic_domain
 from sww_boundary_conditions import *
 from sww_forcing_terms import *
@@ -65,5 +65,5 @@
                                 tagged_elements      = tagged_elements)
         
-        from config import minimum_allowed_height, g, h0
+        from flow_1d.config import minimum_allowed_height, g, h0
         self.minimum_allowed_height = minimum_allowed_height
         self.g = g
@@ -84,5 +84,5 @@
         
         #Reduction operation for get_vertex_values
-        from util import mean
+        from flow_1d.utilities.util import mean
         self.reduction = mean
         #self.reduction = min  #Looks better near steep slopes
Index: anuga_work/development/flow_1d/sww_flow/sww_forcing_terms.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/sww_forcing_terms.py	(revision 7827)
+++ anuga_work/development/flow_1d/sww_flow/sww_forcing_terms.py	(revision 7830)
@@ -14,5 +14,5 @@
     """
 
-    from util import gradient
+    from flow_1d.utilities.util import gradient
 
     xmom  = domain.quantities['xmomentum'].explicit_update
@@ -158,5 +158,5 @@
     else:
         try:
-            f = numpy.float(f)
+            f = float(f)
         except:
             msg = 'Force field %s must be either a scalar' %f
@@ -166,4 +166,2 @@
 
 
-if __name__ == "__main__":
-    print "Hello World";
Index: anuga_work/development/flow_1d/sww_flow/sww_vel_comp_flux_ext.c
===================================================================
--- anuga_work/development/flow_1d/sww_flow/sww_vel_comp_flux_ext.c	(revision 7830)
+++ anuga_work/development/flow_1d/sww_flow/sww_vel_comp_flux_ext.c	(revision 7830)
@@ -0,0 +1,638 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+const double pi = 3.14159265358979;
+
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+/* double max(double a, double b) { */
+/* 	double z; */
+/* 	z=(a>b)?a:b; */
+/* 	return z;} */
+	
+/* double min(double a, double b) { */
+/* 	double z; */
+/* 	z=(a<b)?a:b; */
+/* 	return z;} */
+
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+		      double *h, 
+		      double epsilon, 
+		      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+     u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);    
+  }
+  
+
+  // Adjust momentum to be consistent with speed
+  *uh = u * *h;
+  
+  return u;
+}
+
+//-------------------------------------------------------------
+// New vel based code
+//-------------------------------------------------------------
+//Innermost flux function (using w=z+h)
+int _flux_function_vel(double *q_left, double *q_right,
+		       double normals, double g, double epsilon, double h0,
+		       double *edgeflux, double *max_speed) {
+    
+    int i;
+    double flux_left[2], flux_right[2];
+    double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
+    double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
+    double z, s_max, s_min, denom;
+    
+
+    w_left  = q_left[0];
+    uh_left = q_left[1]*normals;
+    z_left  = q_left[2];
+    h_left  = q_left[3];
+    u_left  = q_left[4]*normals;
+
+/*     printf("w_left  = %f \n",w_left); */
+/*     printf("uh_left = %f \n",uh_left); */
+/*     printf("z_left  = %f \n",z_left); */
+/*     printf("h_left  = %f \n",h_left); */
+/*     printf("u_left  = %f \n",u_left); */
+    
+    w_right  = q_right[0];
+    uh_right = q_right[1]*normals;
+    z_right  = q_right[2];
+    h_right  = q_right[3];
+    u_right  = q_right[4]*normals;		
+    
+    z = (z_left+z_right)/2.0;		   
+    
+    soundspeed_left = sqrt(g*h_left);
+    soundspeed_right = sqrt(g*h_right);
+    
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
+    if (s_max < 0.0) s_max = 0.0;
+    
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
+    if (s_min > 0.0) s_min = 0.0;
+    
+    
+    // Flux formulas
+    flux_left[0] = u_left*h_left;
+    flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
+    
+    flux_right[0] = u_right*h_right;
+    flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
+    
+    // Flux computation
+    denom = s_max-s_min;
+    if (denom < epsilon) {
+	for (i=0; i<2; i++) edgeflux[i] = 0.0;
+	*max_speed = 0.0;
+    } else {
+	edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
+	edgeflux[0] += s_max*s_min*(w_right-w_left);
+	edgeflux[0] /= denom;
+	edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
+	edgeflux[1] += s_max*s_min*(uh_right-uh_left);
+	edgeflux[1] /= denom;
+	edgeflux[1] *= normals;
+	
+	// Maximal wavespeed
+	*max_speed = max(fabs(s_max), fabs(s_min));
+    }  
+    return 0;		
+}
+
+// Computational function for flux computation
+double _compute_fluxes_vel_ext(double cfl,
+			       double timestep,
+			       double epsilon,
+			       double g,
+			       double h0,
+			       long* neighbours,
+			       long* neighbour_vertices,
+			       double* normals,
+			       double* areas,
+			       double* stage_edge_values,
+			       double* xmom_edge_values,
+			       double* bed_edge_values,
+			       double* height_edge_values,
+			       double* velocity_edge_values,
+			       double* stage_boundary_values,
+			       double* xmom_boundary_values,
+			       double* bed_boundary_values,
+			       double* height_boundary_values,
+			       double* velocity_boundary_values,
+			       double* stage_explicit_update,
+			       double* xmom_explicit_update,
+			       int number_of_elements,
+			       double* max_speed_array) {
+		
+    double flux[2], ql[5], qr[5], edgeflux[2];
+    double max_speed, normal;
+    int k, i, ki, n, m, nm=0;
+
+    //printf("Inside _comp\n");    
+    
+    for (k=0; k<number_of_elements; k++) {
+	flux[0] = 0.0;
+	flux[1] = 0.0;
+	//printf("k = %d\n",k);
+
+	for (i=0; i<2; i++) {
+	    ki = k*2+i;
+	    
+	    ql[0] = stage_edge_values[ki];
+	    ql[1] = xmom_edge_values[ki];
+	    ql[2] = bed_edge_values[ki];
+	    ql[3] = height_edge_values[ki];
+	    ql[4] = velocity_edge_values[ki];
+	    
+	    n = neighbours[ki];
+	    if (n<0) {
+		m = -n-1;
+		qr[0] = stage_boundary_values[m];
+		qr[1] = xmom_boundary_values[m];
+		qr[2] = bed_boundary_values[m];
+		qr[3] = height_boundary_values[m];
+		qr[4] = velocity_boundary_values[m];
+	    } else {
+		m = neighbour_vertices[ki];
+		nm = n*2+m;
+		qr[0] = stage_edge_values[nm];
+		qr[1] = xmom_edge_values[nm];
+		qr[2] = bed_edge_values[nm];
+		qr[3] = height_edge_values[nm];
+		qr[4] = velocity_edge_values[nm];
+	    }
+	    
+	    normal = normals[ki];
+	    _flux_function_vel(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
+	    flux[0] -= edgeflux[0];
+	    flux[1] -= edgeflux[1];
+	    
+	    // Update timestep based on edge i and possibly neighbour n
+	    if (max_speed > epsilon) {
+		// Original CFL calculation
+		
+		timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
+		if (n>=0) {
+		    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
+		}
+	    }
+	} // End edge i (and neighbour n)
+	flux[0] /= areas[k];
+	stage_explicit_update[k] = flux[0];
+	flux[1] /= areas[k];
+	xmom_explicit_update[k] = flux[1];
+	
+	//Keep track of maximal speeds
+	max_speed_array[k]=max_speed;
+    }
+    return timestep;	
+}
+
+//-------------------------------------------------------------
+// Old code
+//------------------------------------------------------------
+//Innermost flux function (using w=z+h)
+int _flux_function(double *q_left, double *q_right,
+        double z_left, double z_right,
+		   double normals, double g, double epsilon, double h0,
+		double *edgeflux, double *max_speed) {
+		
+		int i;
+		double ql[2], qr[2], flux_left[2], flux_right[2];
+		double z, w_left, h_left, uh_left, soundspeed_left, u_left;
+		double w_right, h_right, uh_right, soundspeed_right, u_right;
+		double s_max, s_min, denom;
+		
+		//printf("h0 = %f \n",h0);
+		ql[0] = q_left[0];
+		ql[1] = q_left[1];
+		ql[1] = ql[1]*normals;
+	
+		qr[0] = q_right[0];
+		qr[1] = q_right[1];
+		qr[1] = qr[1]*normals;
+	  
+		z = (z_left+z_right)/2.0;		   
+				   
+		//w_left = ql[0];
+		//h_left = w_left-z;
+		//uh_left = ql[1];
+		
+
+
+		// Compute speeds in x-direction
+		w_left = ql[0];          
+		h_left = w_left-z;
+		uh_left = ql[1];
+
+		u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
+
+		w_right = qr[0];
+		h_right = w_right-z;
+		uh_right = qr[1];
+
+		u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
+  
+		soundspeed_left = sqrt(g*h_left);
+		soundspeed_right = sqrt(g*h_right);
+		
+		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_right[0] = u_right*h_right;
+		flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
+
+		// Flux computation
+		denom = s_max-s_min;
+		if (denom < epsilon) {
+			for (i=0; i<2; i++) edgeflux[i] = 0.0;
+			*max_speed = 0.0;
+		} else {
+			edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
+			edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
+			edgeflux[0] /= denom;
+			edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
+			edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
+			edgeflux[1] /= denom;
+			edgeflux[1] *= normals;
+    
+		// Maximal wavespeed
+        *max_speed = max(fabs(s_max), fabs(s_min));
+		}  
+    return 0;		
+	}
+		
+		
+	
+	
+// Computational function for flux computation
+double _compute_fluxes_ext(
+			   double cfl,
+			   double timestep,
+			   double epsilon,
+			   double g,
+			   double h0,
+			   long* neighbours,
+			   long* neighbour_vertices,
+			   double* normals,
+			   double* areas,
+			   double* stage_edge_values,
+			   double* xmom_edge_values,
+			   double* bed_edge_values,
+			   double* stage_boundary_values,
+			   double* xmom_boundary_values,
+			   double* stage_explicit_update,
+			   double* xmom_explicit_update,
+			   int number_of_elements,
+			   double* max_speed_array) {
+		
+                double flux[2], ql[2], qr[2], edgeflux[2];
+		double zl, zr, max_speed, normal;
+		int k, i, ki, n, m, nm=0;
+		
+		
+		for (k=0; k<number_of_elements; k++) {
+			flux[0] = 0.0;
+			flux[1] = 0.0;
+			
+			for (i=0; i<2; i++) {
+				ki = k*2+i;
+				
+				ql[0] = stage_edge_values[ki];
+				ql[1] = xmom_edge_values[ki];
+				zl = bed_edge_values[ki];
+				
+				n = neighbours[ki];
+				if (n<0) {
+					m = -n-1;
+					qr[0] = stage_boundary_values[m];
+					qr[1] = xmom_boundary_values[m];
+					zr = zl;
+				} else {
+					m = neighbour_vertices[ki];
+					nm = n*2+m;
+					qr[0] = stage_edge_values[nm];
+					qr[1] = xmom_edge_values[nm];
+					zr = bed_edge_values[nm];				
+				}
+				
+				normal = normals[ki];
+				_flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
+				flux[0] -= edgeflux[0];
+				flux[1] -= edgeflux[1];
+				
+				// Update timestep based on edge i and possibly neighbour n
+				if (max_speed > epsilon) {
+				    // Original CFL calculation
+				    
+				    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
+				    if (n>=0) {
+					timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
+				    }
+				}
+            } // End edge i (and neighbour n)
+			flux[0] /= areas[k];
+			stage_explicit_update[k] = flux[0];
+			flux[1] /= areas[k];
+			xmom_explicit_update[k] = flux[1];
+			
+			//Keep track of maximal speeds
+			max_speed_array[k]=max_speed;
+		}
+		return timestep;	
+	}
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
+  
+   PyObject 
+	*domain,
+	*stage, 
+	*xmom, 
+	*bed;
+
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_vertices,
+	*normals, 
+	*areas,
+	*stage_vertex_values,
+	*xmom_vertex_values,
+	*bed_vertex_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+    
+  double timestep, epsilon, g, h0, cfl;
+  int number_of_elements;
+
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOO",
+			&timestep,
+			&domain,
+			&stage,
+			&xmom,
+			&bed)) {
+      PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_ext could not parse input");
+      return NULL;
+  }
+
+
+    epsilon           = get_python_double(domain,"epsilon");
+    g                 = get_python_double(domain,"g");
+    h0                = get_python_double(domain,"h0");
+    cfl               = get_python_double(domain,"CFL");
+ 
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
+    normals           = get_consecutive_array(domain, "normals");
+    areas             = get_consecutive_array(domain, "areas");    
+    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
+    
+    stage_vertex_values = get_consecutive_array(stage, "vertex_values");    
+    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");    
+    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+
+
+
+    number_of_elements = stage_vertex_values -> dimensions[0];
+
+
+ 
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
+    timestep = _compute_fluxes_ext(
+	cfl,
+	timestep,
+	epsilon,
+	g,
+	h0,
+	(long*) neighbours -> data,
+	(long*) neighbour_vertices -> data,
+	(double*) normals -> data,
+	(double*) areas -> data,
+	(double*) stage_vertex_values -> data,
+	(double*) xmom_vertex_values -> data,
+	(double*) bed_vertex_values -> data,
+	(double*) stage_boundary_values -> data,
+	(double*) xmom_boundary_values -> data,
+	(double*) stage_explicit_update -> data,
+	(double*) xmom_explicit_update -> data,
+	number_of_elements,
+	(double*) max_speed_array -> data);
+
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_vertices);
+  Py_DECREF(normals);
+  Py_DECREF(areas);
+  Py_DECREF(stage_vertex_values);
+  Py_DECREF(xmom_vertex_values);
+  Py_DECREF(bed_vertex_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(max_speed_array);
+
+
+
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+//------------------------------------------------
+// New velocity based compute fluxes
+//------------------------------------------------
+PyObject *compute_fluxes_vel_ext(PyObject *self, PyObject *args) {
+  
+    PyObject 
+	*domain,
+	*stage, 
+	*xmom, 
+	*bed,
+	*height,
+	*velocity;
+
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_vertices,
+	*normals, 
+	*areas,
+	*stage_vertex_values,
+	*xmom_vertex_values,
+	*bed_vertex_values,
+	*height_vertex_values,
+	*velocity_vertex_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*bed_boundary_values,
+	*height_boundary_values,
+	*velocity_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*max_speed_array;
+    
+  double timestep, epsilon, g, h0, cfl;
+  int number_of_elements;
+
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOOOO",
+			&timestep,
+			&domain,
+			&stage,
+			&xmom,
+			&bed,
+			&height,
+			&velocity)) {
+      PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
+      printf("comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
+      return NULL;
+  }
+
+
+    epsilon           = get_python_double(domain,"epsilon");
+    g                 = get_python_double(domain,"g");
+    h0                = get_python_double(domain,"h0");
+    cfl               = get_python_double(domain,"CFL");
+ 
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
+    normals           = get_consecutive_array(domain, "normals");
+    areas             = get_consecutive_array(domain, "areas");    
+    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
+    
+    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");    
+    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");    
+    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");    
+    height_vertex_values     = get_consecutive_array(height,   "vertex_values");    
+    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");    
+
+    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");    
+    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");    
+    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");    
+    height_boundary_values    = get_consecutive_array(height,    "boundary_values");    
+    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");    
+
+    number_of_elements = stage_vertex_values -> dimensions[0];
+ 
+
+    //printf("comp_flux_vel 1 N = %d %15.5e %15.5e %15.5e %15.5e %15.5e \n",number_of_elements,cfl,timestep,epsilon,g,h0);
+
+
+
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
+    timestep = _compute_fluxes_vel_ext(cfl,
+				       timestep,
+				       epsilon,
+				       g,
+				       h0,
+				       (long*) neighbours -> data,
+				       (long*) neighbour_vertices -> data,
+				       (double*) normals -> data,
+				       (double*) areas -> data,
+				       (double*) stage_vertex_values -> data,
+				       (double*) xmom_vertex_values -> data,
+				       (double*) bed_vertex_values -> data,
+				       (double*) height_vertex_values -> data,
+				       (double*) velocity_vertex_values -> data,
+				       (double*) stage_boundary_values -> data,
+				       (double*) xmom_boundary_values -> data,
+				       (double*) bed_boundary_values -> data,
+				       (double*) height_boundary_values -> data,
+				       (double*) velocity_boundary_values -> data,
+				       (double*) stage_explicit_update -> data,
+				       (double*) xmom_explicit_update -> data,
+				       number_of_elements,
+				       (double*) max_speed_array -> data);
+
+    //printf("comp_flux_vel 2 \n");    
+    
+    Py_DECREF(neighbours);
+    Py_DECREF(neighbour_vertices);
+    Py_DECREF(normals);
+    Py_DECREF(areas);
+    Py_DECREF(stage_vertex_values);
+    Py_DECREF(xmom_vertex_values);
+    Py_DECREF(bed_vertex_values);
+    Py_DECREF(height_vertex_values);
+    Py_DECREF(velocity_vertex_values);
+    Py_DECREF(stage_boundary_values);
+    Py_DECREF(xmom_boundary_values);
+    Py_DECREF(bed_boundary_values);
+    Py_DECREF(height_boundary_values);
+    Py_DECREF(velocity_boundary_values);
+    Py_DECREF(stage_explicit_update);
+    Py_DECREF(xmom_explicit_update);
+    Py_DECREF(max_speed_array);
+    
+    //printf("comp_flux_vel 3 \n");
+
+    // Return updated flux timestep
+    return Py_BuildValue("d", timestep);
+}
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+
+static struct PyMethodDef MethodTable[] = {
+  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {"compute_fluxes_vel_ext", compute_fluxes_vel_ext, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initsww_vel_comp_flux_ext(void){
+  Py_InitModule("sww_vel_comp_flux_ext", MethodTable);
+  import_array();
+}
Index: anuga_work/development/flow_1d/sww_flow/sww_vel_domain.py
===================================================================
--- anuga_work/development/flow_1d/sww_flow/sww_vel_domain.py	(revision 7830)
+++ anuga_work/development/flow_1d/sww_flow/sww_vel_domain.py	(revision 7830)
@@ -0,0 +1,364 @@
+"""Class Domain -
+1D interval 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 = S
+
+where
+
+U = [w, uh]
+E = [uh, u^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 respectiely.
+
+The quantities are
+
+symbol    variable name    explanation
+x         x                horizontal 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]
+uh        xmomentum        momentum in the x direction [m^2/s]
+
+eta                        mannings friction coefficient [to appear]
+nu                         wind stress coefficient [to appear]
+
+The conserved quantities are w, uh
+
+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
+
+
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2006
+"""
+
+
+from flow_1d.generic_1d.generic_domain import *
+from sww_boundary_conditions import *
+from sww_forcing_terms import *
+
+#Shallow water domain
+class Domain(Generic_domain):
+
+    def __init__(self, coordinates, boundary = None, tagged_elements = None):
+
+        conserved_quantities = ['stage', 'xmomentum']
+        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
+        other_quantities = ['friction']
+        Generic_domain.__init__(self,
+                                coordinates = coordinates,
+                                boundary = boundary,
+                                conserved_quantities = conserved_quantities,
+                                evolved_quantities = evolved_quantities,
+                                other_quantities = other_quantities,
+                                tagged_elements = tagged_elements)
+        
+        from flow_1d.config import minimum_allowed_height, g, h0
+        self.minimum_allowed_height = minimum_allowed_height
+        self.g = g
+        self.h0 = h0
+
+        #forcing terms not included in 1d domain ?WHy?
+        self.forcing_terms.append(gravity)
+        #self.forcing_terms.append(manning_friction)
+        #print "\nI have Removed forcing terms line 64 1dsw"
+
+        
+        #Stored output
+        self.store = True
+        self.format = 'sww'
+        self.smooth = True
+
+        
+        #Reduction operation for get_vertex_values
+        from flow_1d.utilities.util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
+
+        self.set_quantities_to_be_stored(['stage','xmomentum'])
+
+	self.__doc__ = 'sww_vel_domain'
+
+        self.check_integrity()
+
+
+
+    def 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 = 'First evolved quantity must be "stage"'
+        assert self.evolved_quantities[0] == 'stage', msg
+        msg = 'Second evolved quantity must be "xmomentum"'
+        assert self.evolved_quantities[1] == 'xmomentum', msg
+        msg = 'Third evolved quantity must be "elevation"'
+        assert self.evolved_quantities[2] == 'elevation', msg
+        msg = 'Fourth evolved quantity must be "height"'
+        assert self.evolved_quantities[3] == 'height', msg
+        msg = 'Fifth evolved quantity must be "velocity"'
+        assert self.evolved_quantities[4] == 'velocity', msg
+
+        Generic_domain.check_integrity(self)
+
+    def compute_fluxes(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        compute_fluxes_vel(self)
+
+    def distribute_to_vertices_and_edges(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges_limit_w_u(self)
+        
+
+#=============== End of Shallow Water Domain ===============================
+#-----------------------------------
+# Compute fluxes interface
+#-----------------------------------
+def compute_fluxes(domain):
+    """
+    Python version of compute fluxes (local_compute_fluxes)
+    is available in test_shallow_water_vel_domain.py
+    """
+ 
+
+    from numpy import zeros 
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    stage = domain.quantities['stage']
+    xmom = domain.quantities['xmomentum']
+    bed = domain.quantities['elevation']
+
+
+    from sww_vel_comp_flux_ext import compute_fluxes_ext
+
+    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)
+
+
+#-----------------------------------
+# Compute flux definition with vel
+#-----------------------------------
+def compute_fluxes_vel(domain):
+    from Numeric import zeros, Float
+    import sys
+
+	
+    timestep = float(sys.maxint)
+
+    stage    = domain.quantities['stage']
+    xmom     = domain.quantities['xmomentum']
+    bed      = domain.quantities['elevation']
+    height   = domain.quantities['height']
+    velocity = domain.quantities['velocity']
+
+
+    from sww_vel_comp_flux_ext import compute_fluxes_vel_ext
+
+    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)
+
+
+
+#--------------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_w_u(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 config import optimised_gradient_limiter
+
+    #Remove very thin layers of water
+    #protect_against_infinitesimal_and_negative_heights(domain)  
+
+    import sys
+    from numpy import zeros
+    from config import epsilon, h0
+	
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+
+    #Arrays   
+    w_C   = Stage.centroid_values    
+    uh_C  = Xmom.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+	
+    #print id(h_C)
+##     for i in range(N):
+## 	h_C[i] = w_C[i] - z_C[i]						
+##         if h_C[i] <= 1.0e-12:
+## 	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
+## 	    h_C[i] = 0.0
+## 	    w_C[i] = z_C[i]
+## 	    #uh_C[i] = 0.0
+            
+##  #           u_C[i] = 0.0
+##  #       else:
+##  #           u_C[i] = uh_C[i]/h_C[i]
+		
+    h0 = 1.0e-12    
+    for i in range(N):
+        h_C[i] = w_C[i] - z_C[i]						
+	if h_C[i] < 1.0e-12:
+	    u_C[i] = 0.0  #Could have been negative
+	    h_C[i] = 0.0
+            w_C[i] = z_C[i]
+	else:
+            #u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
+            u_C[i]  = uh_C[i]/h_C[i]
+	
+    for name in [ 'velocity', 'stage' ]:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    w_V  = domain.quantities['stage'].vertex_values			
+    z_V  = domain.quantities['elevation'].vertex_values	
+    h_V  = domain.quantities['height'].vertex_values
+    u_V  = domain.quantities['velocity'].vertex_values		
+    uh_V = domain.quantities['xmomentum'].vertex_values	
+
+    #print w_V
+    #print z_V
+    
+    h_V[:,:]    = w_V - z_V
+    for i in range(len(h_C)):
+        for j in range(2):
+            if h_V[i,j] < 0.0 :
+                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
+                dh = h_V[i,j]
+                h_V[i,j] = 0.0
+                w_V[i,j] = z_V[i,j]
+                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
+                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
+                
+    uh_V[:,:] = u_V * h_V
+
+    
+    return
+
+#---------------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_w_uh(domain):
+    """Distribution from centroids to vertices specific to the
+    shallow water wave equation.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    import sys
+    from Numeric import zeros, Float
+    from config import epsilon, h0
+	
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Bed = domain.quantities['elevation']
+    Height = domain.quantities['height']
+    Velocity = domain.quantities['velocity']
+
+    #Arrays   
+    w_C   = Stage.centroid_values    
+    uh_C  = Xmom.centroid_values    
+    z_C   = Bed.centroid_values
+    h_C   = Height.centroid_values
+    u_C   = Velocity.centroid_values
+	
+
+    for i in range(N):
+	h_C[i] = w_C[i] - z_C[i]						
+        if h_C[i] <= 1.0e-6:
+	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
+	    h_C[i] = 0.0
+	    w_C[i] = z_C[i]
+	    uh_C[i] = 0.0
+            
+    for name in [ 'stage', 'xmomentum']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    w_V  = domain.quantities['stage'].vertex_values			
+    z_V  = domain.quantities['elevation'].vertex_values	
+    h_V  = domain.quantities['height'].vertex_values
+    u_V  = domain.quantities['velocity'].vertex_values		
+    uh_V = domain.quantities['xmomentum'].vertex_values	
+
+    h_V[:,:]    = w_V - z_V
+
+    for i in range(len(h_C)):
+        for j in range(2):
+            if h_V[i,j] < 0.0 :
+                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
+                dh = h_V[i,j]
+                h_V[i,j] = 0.0
+                w_V[i,j] = z_V[i,j]
+                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
+                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
+                u_V[i,j] = 0.0
+            if h_V[i,j] < h0:
+                u_V[i,j]
+            else:
+                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
+
+
