Index: uga_work/development/pipeflow/comp_flux_ext.c
===================================================================
--- /anuga_work/development/pipeflow/comp_flux_ext.c	(revision 7787)
+++ 	(revision )
@@ -1,399 +1,0 @@
-#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"
-
-
-//Innermost flux function (using w=z+h)
-int _flux_function(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;
-		
-		//printf("h0 = %f \n",h0);
-		w_left  = q_left[0];
-		uh_left = q_left[1]*normals;
-		h_left  = q_left[2];
-		z_left  = q_left[3];
-		u_left  = q_left[4]*normals;
-	
-		w_right  = q_right[0];
-		uh_right = q_right[1]*normals;
-		h_right  = q_right[2];
-		z_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*(q_right[0]-q_left[0]);
-			edgeflux[0] /= denom;
-			edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
-			edgeflux[1] += s_max*s_min*(q_right[1]-q_left[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* 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;
-		
-		
-		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];
-				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_edge_values[m];
-					qr[3] = height_edge_values[m];
-					qr[4] = velocity_edge_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(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;	
-	}
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  
-    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, "dddddOOOOOOOOOOOOOOOOiO",
-			  &cfl,
-			  &timestep,
-			  &epsilon,
-			  &g,
-			  &h0,
-			  &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,
-			  &number_of_elements,
-			  &max_speed_array)) {
-	PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
-	return NULL;
-    }
-
-
-    // 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*) 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);
-
-
-    // Return updated flux timestep
-    return Py_BuildValue("d", timestep);
-}
-
-
-PyObject *compute_fluxes_ext_short(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_ext.c: compute_fluxes_ext_short 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];
- 
-    // 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*) 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);
-    
-
-  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);
-
-
-  // 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_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
-  {NULL, NULL}
-};
-
-// Module initialisation
-void initcomp_flux_ext(void){
-  Py_InitModule("comp_flux_ext", MethodTable);
-  import_array();
-}
Index: uga_work/development/pipeflow/comp_flux_ext_wellbalanced.c
===================================================================
--- /anuga_work/development/pipeflow/comp_flux_ext_wellbalanced.c	(revision 7787)
+++ 	(revision )
@@ -1,329 +1,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-#include "util_ext.h"
-const double pi = 3.14159265358979;
-
-//Innermost flux function (using w=z+h)
-int _flux_function_wellbalanced(double *q_left, 
-								double *q_right,
-								double normals, 
-								double g, 
-								double epsilon,
-								double *edgeflux, 
-								double *max_speed) {		
-		int i;
-        double z_left, z_right;
-		double flux_left[2], flux_right[2];
-		double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
-		double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;
-		double s_max, s_min, denom;
-		
-		w_left  = q_left[0];
-		uh_left = q_left[1];
-		uh_left = uh_left*normals;
-        z_left  = q_left[2];
-		h_left  = q_left[3];
-		u_left  = q_left[4];
-	
-		w_right  = q_right[0];
-		uh_right = q_right[1];
-		uh_right = uh_right*normals;
-        z_right  = q_right[2];
-		h_right  = q_right[3];
-		u_right  = q_right[4];
-	  
-		if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
-        
-		z_star = max(z_left, z_right);		              											//equation(7)  
-				   
-		// Compute speeds in x-direction.
-		h_left_star = max(0.0, w_left-z_star);														//equation(8)
-		
-		if (h_left_star < epsilon) {
-			u_left = 0.0; h_left_star = 0.0;
-		}
-		else {
-			u_left = uh_left/h_left_star;
-		}
-	                                                                                                                                                                                           
-		h_right_star = max(0.0, w_right-z_star);													//equation(9)
-
-		if (h_right_star < epsilon) {
-			u_right = 0.0; h_right_star = 0.0; 
-		}
-		else {
-			u_right = uh_right/h_right_star;
-		}
-  
-		//printf("u_left = %g u_right = %g \n",u_left, u_right);
-		soundspeed_left = sqrt(g*h_left_star);
-		soundspeed_right = sqrt(g*h_right_star);
-		
-		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_star;
-		flux_left[1] = u_left*u_left*h_left_star + 0.5*g*h_left_star*h_left_star;
-
-		flux_right[0] = u_right*h_right_star;
-		flux_right[1] = u_right*u_right*h_right_star + 0.5*g*h_right_star*h_right_star;
-
-		// Flux computation
-		denom = s_max-s_min;
-		if (denom < epsilon && z_right > z_left) {
-			for (i=0; i<2; i++) edgeflux[i] = 0.0;
-			// Maximal wavespeed			
-			edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
-			edgeflux[1] *= normals;
-			*max_speed = 0.0;
-		} else if (denom < epsilon){
-			for (i=0; i<2; i++) edgeflux[i] = 0.0;
-			// Maximal wavespeed
-			*max_speed = 0.0;		
-		}else {
-			edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
-			edgeflux[0] += s_max*s_min*(h_right_star-h_left_star);                   //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*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
-			edgeflux[1] /= denom;
-			edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
-			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_wellbalanced(double timestep,
-		double epsilon,
-		double g,
-		
-		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* vel_edge_values,
-		
-		double* stage_boundary_values,
-		double* xmom_boundary_values,
-		double* bed_boundary_values,
-		double* height_boundary_values,
-		double* vel_boundary_values,
-		
-		double* stage_explicit_update,
-		double* xmom_explicit_update,
-		double* bed_explicit_update,
-		double* height_explicit_update,
-		double* vel_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;
-		
-		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];
-				ql[2] = bed_edge_values[ki];
-                ql[3] = height_edge_values[ki];
-				ql[4] = vel_edge_values[ki];
-                
-				n = neighbours[ki];
-				printf("neighbours[ki] = %li \n",neighbours[ki]);
-				if (n<0) {
-					m = -n-1;
-					qr[0] = stage_boundary_values[m];
-					qr[1] = xmom_boundary_values[m];
-					qr[2] = ql[2];																			
-					qr[3] = height_edge_values[m]; 														
-					qr[4] = vel_edge_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] = vel_edge_values[nm];															
-				}
-				
-				normal = normals[ki];
-				_flux_function_wellbalanced(ql, qr, normal, g, epsilon, 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*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
-					if (n>=0) {
-						timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. 
-					}
-				}
-            } // 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_wellbalanced(PyObject *self, PyObject *args) {
-  
-    PyArrayObject 
-	*neighbours, 
-	*neighbour_vertices,
-	*normals, 
-	*areas,
-	
-	*stage_edge_values,
-	*xmom_edge_values,
-	*bed_edge_values,
-	*height_edge_values,
-	*vel_edge_values,
-	
-	*stage_boundary_values,
-	*xmom_boundary_values,
-	*bed_boundary_values,
-	*height_boundary_values,
-	*vel_boundary_values,
-	
-	*stage_explicit_update,
-	*xmom_explicit_update,
-	*bed_explicit_update,
-	*height_explicit_update,
-	*vel_explicit_update,
-	
-	*max_speed_array;
-    
-  double timestep, epsilon, g;
-  int number_of_elements;
-    
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
-			&timestep,
-			&epsilon,
-			&g,
-			
-			&neighbours,
-			&neighbour_vertices,
-			&normals,
-			&areas,
-			
-			&stage_edge_values,
-			&xmom_edge_values,
-			&bed_edge_values,
-			&height_edge_values,
-			&vel_edge_values,
-			
-			&stage_boundary_values,
-			&xmom_boundary_values,
-			&bed_boundary_values,
-			&height_boundary_values,
-			&vel_boundary_values,
-			
-			&stage_explicit_update,
-			&xmom_explicit_update,
-			&bed_explicit_update,
-			&height_explicit_update,
-			&vel_explicit_update,
-			
-			&number_of_elements,
-			&max_speed_array)) {
-    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input");
-    return NULL;
-  }
-
- 
-  // Call underlying flux computation routine and update 
-  // the explicit update arrays 
-  timestep = _compute_fluxes_ext_wellbalanced(timestep,
-				 epsilon,
-				 g,
-				 
-				 (long*) neighbours -> data,
-				 (long*) neighbour_vertices -> data,
-				 (double*) normals -> data,
-				 (double*) areas -> data,
-				 
-				 (double*) stage_edge_values -> data,
-				 (double*) xmom_edge_values -> data,
-				 (double*) bed_edge_values -> data,
-				 (double*) height_edge_values -> data,
-				 (double*) vel_edge_values -> data,
-				 
-				 (double*) stage_boundary_values -> data,
-				 (double*) xmom_boundary_values -> data,
-				 (double*) bed_boundary_values -> data,
-				 (double*) height_boundary_values -> data,
-				 (double*) vel_boundary_values -> data,
-				 
-				 (double*) stage_explicit_update -> data,
-				 (double*) xmom_explicit_update -> data,
-				 (double*) bed_explicit_update -> data,
-				 (double*) height_explicit_update -> data,
-				 (double*) vel_explicit_update -> data,
-				 
-				 number_of_elements,
-				 (double*) max_speed_array -> data);
-
-  // Return updated flux timestep
-  return Py_BuildValue("d", timestep);
-}
-
-
-
-//-------------------------------
-// Method table for python module
-//-------------------------------
-
-static struct PyMethodDef MethodTable[] = {
-  {"compute_fluxes_ext_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
-  {NULL, NULL}
-};
-
-// Module initialisation
-void initcomp_flux_ext_wellbalanced(void){
-  Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
-  import_array();
-}
Index: uga_work/development/pipeflow/domain.py
===================================================================
--- /anuga_work/development/pipeflow/domain.py	(revision 7787)
+++ 	(revision )
@@ -1,1357 +1,0 @@
-"""Class Domain - 1D domains for finite-volume computations of
-   the shallow water wave equation
-
-
-   Copyright 2004
-   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-   Geoscience Australia
-"""
-
-import numpy
-from generic_boundary_conditions import *
-
-
-class Domain:
-
-    def __init__(self,
-                 coordinates,
-                 boundary = None,
-                 conserved_quantities = None,
-                 evolved_quantities = None,
-                 other_quantities = None,
-                 tagged_elements = None):
-        """
-        Build 1D elements from x coordinates
-        """
-
-        from config import timestepping_method
-        from config import CFL
-        
-        #Store Points
-        self.coordinates = numpy.array(coordinates, numpy.float)
-
-
-        #Register number of Elements
-        self.number_of_elements = N = len(self.coordinates)-1
-
-        self.beta = 1.0
-        self.set_limiter("minmod_kurganov")
-        self.set_CFL(CFL)
-        self.set_timestepping_method(timestepping_method)
-
-        self.wet_nodes = numpy.zeros((N,2), numpy.int) # should this be here
-
-        #Allocate space for neighbour and boundary structures
-        self.neighbours = numpy.zeros((N, 2), numpy.int)
-        #self.neighbour_edges = numpy.zeros((N, 2), numpy.int)
-        self.neighbour_vertices   = numpy.zeros((N, 2), numpy.int)
-        self.number_of_boundaries = numpy.zeros(N, numpy.int)
-        self.surrogate_neighbours = numpy.zeros((N, 2), numpy.int)
-        
-        #Allocate space for geometric quantities
-        self.vertices  = numpy.zeros((N, 2), numpy.float)
-        self.centroids = numpy.zeros(N, numpy.float)
-        self.areas     = numpy.zeros(N, numpy.float)
-
-        self.max_speed_array = numpy.zeros(N, numpy.float)
-      
-
-        self.normals = numpy.zeros((N, 2), numpy.float)
-        
-        for i in range(N):
-            xl = self.coordinates[i]
-            xr = self.coordinates[i+1]
-            self.vertices[i,0] = xl
-            self.vertices[i,1] = xr
-
-            centroid = (xl+xr)/2.0
-            self.centroids[i] = centroid
-
-            msg = 'Coordinates should be ordered, smallest to largest'
-            assert xr>xl, msg
-            
-            #The normal vectors
-            # - point outward from each edge
-            # - are orthogonal to the edge
-            # - have unit length
-            # - Are enumerated by left vertex then right vertex normals
-
-            nl = -1.0
-            nr =  1.0
-            self.normals[i,:] = [nl, nr]
-
-            self.areas[i] = (xr-xl)
-            
-# #         print 'N', N
-# #         print 'Centroid', self.centroids
-# #         print 'Areas', self.areas
-# #         print 'Vertex_Coordinates', self.vertices
-            
-            #Initialise Neighbours (-1 means that it is a boundary neighbour)
-            self.neighbours[i, :] = [-1, -1]
-            #Initialise edge ids of neighbours
-            #Initialise vertex ids of neighbours
-            #In case of boundaries this slot is not used
-            #self.neighbour_edges[i, :] = [-1, -1]
-            self.neighbour_vertices[i, :] = [-1, -1]
-
-        self.build_vertexlist()
-
-        #Build neighbour structure
-        self.build_neighbour_structure()
-
-        #Build surrogate neighbour structure
-        self.build_surrogate_neighbour_structure()          
-
-        #Build boundary dictionary mapping (id, edge) to symbolic tags
-        #Build boundary dictionary mapping (id, vertex) to symbolic tags
-        self.build_boundary_dictionary(boundary)
-
-        #Build tagged element  dictionary mapping (tag) to array of elements
-        self.build_tagged_elements_dictionary(tagged_elements)
-        
-        from quantity import Quantity
-        #from quantity_domain import Quantity, Conserved_quantity
-        
-        #List of quantity names entering
-        #the conservation equations
-        #(Must be a subset of quantities)
-        if conserved_quantities is None:
-            self.conserved_quantities = []
-        else:
-            self.conserved_quantities = conserved_quantities
-
-        if evolved_quantities is None:
-            self.evolved_quantities = self.conserved_quantities
-        else:
-            self.evolved_quantities = evolved_quantities
-
-        if other_quantities is None:
-            self.other_quantities = []
-        else:
-            self.other_quantities = other_quantities
-
-
-        #Build dictionary of Quantity instances keyed by quantity names
-        self.quantities = {}
-
-        #print self.conserved_quantities
-        #print self.evolved_quantities
-        
-
-        #FIXME: remove later - maybe OK, though....
-        for name in self.evolved_quantities:
-            self.quantities[name] = Quantity(self)
-        for name in self.other_quantities:
-            self.quantities[name] = Quantity(self)
-
-        #Create an empty list for explicit forcing terms
-        self.forcing_terms = []
-
-        #Defaults
-        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
-        self.beta_w = beta_w
-        self.beta_h = beta_h
-        self.epsilon = epsilon
-
-        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
-        #Or maybe get rid of order altogether and use beta_w and beta_h
-        self.default_order = 1
-        self.order = self.default_order
-
-        self.default_time_order = 1
-        self.time_order = self.default_time_order
-        
-        self.smallsteps = 0
-        self.max_smallsteps = max_smallsteps
-        self.number_of_steps = 0
-        self.number_of_first_order_steps = 0
-
-        #Model time
-        self.time = 0.0
-        self.finaltime = None
-        self.min_timestep = self.max_timestep = 0.0
-        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
-        #Checkpointing and storage
-        from config import default_datadir
-        self.set_datadir(default_datadir)
-        self.filename = 'domain'
-        self.checkpoint = False
-
-    def __len__(self):
-        return self.number_of_elements
-
-    def build_vertexlist(self):
-        """Build vertexlist index by vertex ids and for each entry (point id)
-        build a list of (triangles, vertex_id) pairs that use the point
-        as vertex.
-
-        Preconditions:
-          self.coordinates and self.triangles are defined
-
-        Postcondition:
-          self.vertexlist is built
-        """
-
-        vertexlist = [None]*len(self.coordinates)
-        for i in range(self.number_of_elements):
-
-            #a = self.triangles[i, 0]
-            #b = self.triangles[i, 1]
-            #c = self.triangles[i, 2]
-            a = i
-            b = i + 1
-
-            #Register the vertices v as lists of
-            #(triangle_id, vertex_id) tuples associated with them
-            #This is used for smoothing
-            #for vertex_id, v in enumerate([a,b,c]):
-            for vertex_id, v in enumerate([a,b]):
-                if vertexlist[v] is None:
-                    vertexlist[v] = []
-
-                vertexlist[v].append( (i, vertex_id) )
-
-        self.vertexlist = vertexlist
-
-        
-    def build_neighbour_structure(self):
-        """Update all registered triangles to point to their neighbours.
-
-        Also, keep a tally of the number of boundaries for each triangle
-
-        Postconditions:
-          neighbours and neighbour_edges is populated
-          neighbours and neighbour_vertices is populated
-          number_of_boundaries integer array is defined.
-        """
-
-        #Step 1:
-        #Build dictionary mapping from segments (2-tuple of points)
-        #to left hand side edge (facing neighbouring triangle)
-
-        N = self.number_of_elements
-        neighbourdict = {}
-        #l_edge = 0
-        #r_edge = 1
-        l_vertex = 0
-        r_vertex = 1
-        for i in range(N):
-
-            #Register all segments as keys mapping to current triangle
-            #and segment id
-            #a = self.triangles[i, 0]
-            #b = self.triangles[i, 1]
-            #c = self.triangles[i, 2]
-            a = self.vertices[i,0]
-            b = self.vertices[i,1]
-            
-            """
-            if neighbourdict.has_key((a,b)):
-                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
-                    raise msg
-            if neighbourdict.has_key((b,c)):
-                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
-                    raise msg
-            if neighbourdict.has_key((c,a)):
-                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
-                    raise msg
-                    """
-            #neighbourdict[a,b] = (i, 2) #(id, edge)
-            #neighbourdict[b,c] = (i, 0) #(id, edge)
-            #neighbourdict[c,a] = (i, 1) #(id, edge)
-            #neighbourdict[a,b] = (i, 1) #(id, edge)
-            #neighbourdict[b,a] = (i, 0) #(id, edge)
-            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
-            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
-            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
-            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
-
-
-        #Step 2:
-        #Go through triangles again, but this time
-        #reverse direction of segments and lookup neighbours.
-        for i in range(N):
-            #a = self.triangles[i, 0]
-            #b = self.triangles[i, 1]
-            #c = self.triangles[i, 2]
-            
-            a = self.vertices[i,0]
-            b = self.vertices[i,1]
-            
-            #self.number_of_boundaries[i] = 3
-            self.number_of_boundaries[i] = 2
-            #if neighbourdict.has_key((b,l_edge)):
-            if neighbourdict.has_key((b,l_vertex)):
-                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
-                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
-                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
-                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
-                self.number_of_boundaries[i] -= 1
-
-            #if neighbourdict.has_key((a,r_edge)):
-            if neighbourdict.has_key((a,r_vertex)):
-                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
-                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
-                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
-                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
-                self.number_of_boundaries[i] -= 1
-                
-            #if neighbourdict.has_key((b,a)):
-            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
-            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
-            #    self.number_of_boundaries[i] -= 1
-
-            #if neighbourdict.has_key((c,b)):
-            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
-            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
-            #    self.number_of_boundaries[i] -= 1
-
-            #if neighbourdict.has_key((a,b)):
-            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
-            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
-            #    self.number_of_boundaries[i] -= 1
-
-    def build_surrogate_neighbour_structure(self):
-        """Build structure where each triangle edge points to its neighbours
-        if they exist. Otherwise point to the triangle itself.
-
-        The surrogate neighbour structure is useful for computing gradients
-        based on centroid values of neighbours.
-
-        Precondition: Neighbour structure is defined
-        Postcondition:
-          Surrogate neighbour structure is defined:
-          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
-          triangles.
-
-        """
-
-        N = self.number_of_elements
-        for i in range(N):
-            #Find all neighbouring volumes that are not boundaries
-            #for k in range(3):
-            for k in range(2):    
-                if self.neighbours[i, k] < 0:
-                    self.surrogate_neighbours[i, k] = i #Point this triangle
-                else:
-                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
-
-    def build_boundary_dictionary(self, boundary = None):
-        """Build or check the dictionary of boundary tags.
-         self.boundary is a dictionary of tags,
-         keyed by volume id and edge:
-         { (id, edge): tag, ... }
-
-         Postconditions:
-            self.boundary is defined.
-        """
-
-        from config import default_boundary_tag
-
-        if boundary is None:
-            boundary = {}
-            for vol_id in range(self.number_of_elements):
-                #for edge_id in range(0, 3):
-                #for edge_id in range(0, 2):
-                for vertex_id in range(0, 2):    
-                    #if self.neighbours[vol_id, edge_id] < 0:
-                    if self.neighbours[vol_id, vertex_id] < 0:   
-                        #boundary[(vol_id, edge_id)] = default_boundary_tag
-                        boundary[(vol_id, vertex_id)] = default_boundary_tag
-        else:
-            #Check that all keys in given boundary exist
-            #for vol_id, edge_id in boundary.keys():
-            for vol_id, vertex_id in boundary.keys():    
-                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
-                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
-                a, b = self.neighbours.shape
-                #assert vol_id < a and edge_id < b, msg
-                assert vol_id < a and vertex_id < b, msg
-
-                #FIXME: This assert violates internal boundaries (delete it)
-                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
-                #assert self.neighbours[vol_id, edge_id] < 0, msg
-
-            #Check that all boundary segments are assigned a tag
-            for vol_id in range(self.number_of_elements):
-                #for edge_id in range(0, 3):
-                #for edge_id in range(0, 2):
-                for vertex_id in range(0, 2):    
-                    #if self.neighbours[vol_id, edge_id] < 0:
-                    if self.neighbours[vol_id, vertex_id] < 0:    
-                        #if not boundary.has_key( (vol_id, edge_id) ):
-                        if not boundary.has_key( (vol_id, vertex_id) ):    
-                            msg = 'WARNING: Given boundary does not contain '
-                            #msg += 'tags for edge (%d, %d). '\
-                            #       %(vol_id, edge_id)
-                            msg += 'tags for vertex (%d, %d). '\
-                                   %(vol_id, vertex_id)
-                            msg += 'Assigning default tag (%s).'\
-                                   %default_boundary_tag
-
-                            #FIXME: Print only as per verbosity
-                            #print msg
-
-                            #FIXME: Make this situation an error in the future
-                            #and make another function which will
-                            #enable default boundary-tags where
-                            #tags a not specified
-                            #boundary[ (vol_id, edge_id) ] =\
-                            boundary[ (vol_id, vertex_id) ] =\
-                                      default_boundary_tag
-
-
-
-        self.boundary = boundary
-
-    def build_tagged_elements_dictionary(self, tagged_elements = None):
-        """Build the dictionary of element tags.
-         self.tagged_elements is a dictionary of element arrays,
-         keyed by tag:
-         { (tag): [e1, e2, e3..] }
-
-         Postconditions:
-            self.element_tag is defined
-        """
-
-        if tagged_elements is None:
-            tagged_elements = {}
-        else:
-            #Check that all keys in given boundary exist
-            for tag in tagged_elements.keys():
-                tagged_elements[tag] = array(tagged_elements[tag]).astype(numpy.int)
-
-                msg = 'Not all elements exist. '
-                assert max(tagged_elements[tag]) < self.number_of_elements, msg
-        #print "tagged_elements", tagged_elements
-        self.tagged_elements = tagged_elements
-
-
-    def set_quantities_to_be_stored(self, q):
-        """Specify which quantities will be stored in the sww file.
-
-        q must be either:
-          - the name of a quantity
-          - a list of quantity names
-          - None
-
-        In the two first cases, the named quantities will be stored at each
-        yieldstep
-        (This is in addition to the quantities elevation and friction)  
-        If q is None, storage will be switched off altogether.
-        """
-
-
-        if q is None:
-            self.quantities_to_be_stored = []            
-            self.store = False
-            return
-
-        if isinstance(q, basestring):
-            q = [q] # Turn argument into a list
-
-        #Check correcness    
-        for quantity_name in q:
-            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
-            assert quantity_name in self.conserved_quantities, msg 
-        
-        self.quantities_to_be_stored = q
-        
-
-
-
-
-    def get_boundary_tags(self):
-        """Return list of available boundary tags
-        """
-
-        tags = {}
-        for v in self.boundary.values():
-            tags[v] = 1
-
-        return tags.keys()
-        
-    def get_vertex_coordinates(self, obj = False):
-        """Return all vertex coordinates.
-        Return all vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-
-        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
-        FIXME, we might make that the default.
-        FIXME Maybe use keyword: continuous = False for this condition?
-
-        
-        """
-
-        if obj is True:
-        
-            #V = self.vertex_coordinates
-            V = self.vertices
-            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
-
-            N = V.shape[0]
-            #return reshape(V, (3*N, 2))
-            return numpy.reshape(V, (N, 2))
-        else:    
-            #return self.vertex_coordinates
-            return self.vertices
-    
-    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
-        """Get conserved quantities at volume vol_id
-
-        If vertex is specified use it as index for vertex values
-        If edge is specified use it as index for edge values
-        If neither are specified use centroid values
-        If both are specified an exeception is raised
-
-        Return value: Vector of length == number_of_conserved quantities
-
-        """
-
-        #if not (vertex is None):# or edge is None):
-        #    msg = 'Values for both vertex and edge was specified.'
-        #    msg += 'Only one (or none) is allowed.'
-       #     raise msg
-
-        q = numpy.zeros( len(self.conserved_quantities), numpy.float)
-
-        for i, name in enumerate(self.conserved_quantities):
-            Q = self.quantities[name]
-            if vertex is not None:
-                q[i] = Q.vertex_values[vol_id, vertex]
-            #elif edge is not None:
-            #    q[i] = Q.edge_values[vol_id, edge]
-            else:
-                q[i] = Q.centroid_values[vol_id]
-
-        return q
-
-
-    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
-        """Get evolved quantities at volume vol_id
-
-        If vertex is specified use it as index for vertex values
-        If edge is specified use it as index for edge values
-        If neither are specified use centroid values
-        If both are specified an exeception is raised
-
-        Return value: Vector of length == number_of_evolved quantities
-
-        """
-
-
-
-        #if not (vertex is None):# or edge is None):
-        #    msg = 'Values for both vertex and edge was specified.'
-        #    msg += 'Only one (or none) is allowed.'
-       #     raise msg
-
-        q = numpy.zeros( len(self.evolved_quantities), numpy.float)
-
-        for i, name in enumerate(self.evolved_quantities):
-            Q = self.quantities[name]
-            if vertex is not None:
-                q[i] = Q.vertex_values[vol_id, vertex]
-            #elif edge is not None:
-            #    q[i] = Q.edge_values[vol_id, edge]
-            else:
-                q[i] = Q.centroid_values[vol_id]
-
-        return q
-        
-
-    def get_centroids(self):
-        """Return all coordinates of centroids
-        Return x coordinate of centroid for each element as a N array
-        """
-
-        return self.centroids
-
-    def get_vertices(self):
-        """Return all coordinates of centroids
-        Return x coordinate of centroid for each element as a N array
-        """
-
-        return self.vertices
-
-    def get_coordinate(self, elem_id, vertex=None):
-        """Return coordinate of centroid,
-        or left or right vertex.
-        Left vertex (vertex=0). Right vertex (vertex=1)
-        """
-
-        if vertex is None:
-            return self.centroids[elem_id]
-        else:
-            return self.vertices[elem_id,vertex]
-
-    def get_area(self, elem_id):
-        """Return area of element id
-        """
-
-        return self.areas[elem_id]
-
-    def get_quantity(self, name, location='vertices', indices = None):
-        """Get values for named quantity
-
-        name: Name of quantity
-
-        In case of location == 'centroids' the dimension values must
-        be a list of a numpyal array of length N, N being the number
-        of elements. Otherwise it must be of dimension Nx3.
-
-        Indices is the set of element ids that the operation applies to.
-
-        The values will be stored in elements following their
-        internal ordering.
-        """
-
-        return self.quantities[name].get_values( location, indices = indices)
-
-    def get_centroid_coordinates(self):
-        """Return all centroid coordinates.
-        Return all centroid coordinates for all triangles as an Nx2 array
-        (ordered as x0, y0 for each triangle)
-        """
-        return self.centroids
-
-
-    def get_timestepping_method(self):
-        return self.timestepping_method
-
-    def set_timestepping_method(self,timestepping_method):
-        
-        if timestepping_method in ['euler', 'rk2', 'rk3']:
-            self.timestepping_method = timestepping_method
-            return
-
-        msg = '%s is an incorrect timestepping type'% timestepping_method
-        raise Exception, msg
-
-
-    def set_quantity(self, name, *args, **kwargs):
-        """Set values for named quantity
-
-
-        One keyword argument is documented here:
-        expression = None, # Arbitrary expression
-
-        expression:
-          Arbitrary expression involving quantity names
-
-        See Quantity.set_values for further documentation.
-        """
-
-        #FIXME (Ole): Allow new quantities here
-        #from quantity import Quantity, Conserved_quantity
-        #Create appropriate quantity object
-        # #if name in self.conserved_quantities:
-        # #    self.quantities[name] = Conserved_quantity(self)
-        # #else:
-        # #    self.quantities[name] = Quantity(self)
-
-
-        #Do the expression stuff
-        if kwargs.has_key('expression'):
-            expression = kwargs['expression']
-            del kwargs['expression']
-
-            Q = self.create_quantity_from_expression(expression)
-            kwargs['quantity'] = Q
-        
-        #Assign values
-        self.quantities[name].set_values(*args, **kwargs)
-
-    def set_boundary(self, boundary_map):
-        """Associate boundary objects with tagged boundary segments.
-
-        Input boundary_map is a dictionary of boundary objects keyed
-        by symbolic tags to matched against tags in the internal dictionary
-        self.boundary.
-
-        As result one pointer to a boundary object is stored for each vertex
-        in the list self.boundary_objects.
-        More entries may point to the same boundary object
-
-        Schematically the mapping is from two dictionaries to one list
-        where the index is used as pointer to the boundary_values arrays
-        within each quantity.
-
-        self.boundary:          (vol_id, edge_id): tag
-        boundary_map (input):   tag: boundary_object
-        ----------------------------------------------
-        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
-
-
-        Pre-condition:
-          self.boundary has been built.
-
-        Post-condition:
-          self.boundary_objects is built
-
-        If a tag from the domain doesn't appear in the input dictionary an
-        exception is raised.
-        However, if a tag is not used to the domain, no error is thrown.
-        FIXME: This would lead to implementation of a
-        default boundary condition
-
-        Note: If a segment is listed in the boundary dictionary and if it is
-        not None, it *will* become a boundary -
-        even if there is a neighbouring triangle.
-        This would be the case for internal boundaries
-
-        Boundary objects that are None will be skipped.
-
-        FIXME: If set_boundary is called multiple times and if Boundary
-        object is changed into None, the neighbour structure will not be
-        restored!!!
-        """
-
-        self.boundary_objects = []
-        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
-
-        #FIXME: Try to remove the sorting and fix test_mesh.py
-        x = self.boundary.keys()
-        x.sort()
-
-        #Loop through edges that lie on the boundary and associate them with
-        #callable boundary objects depending on their tags
-        #for k, (vol_id, edge_id) in enumerate(x):
-        for k, (vol_id, vertex_id) in enumerate(x):
-            #tag = self.boundary[ (vol_id, edge_id) ]
-            tag = self.boundary[ (vol_id, vertex_id) ]
-
-            if boundary_map.has_key(tag):
-                B = boundary_map[tag]  #Get callable boundary object
-
-                if B is not None:
-                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
-                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
-                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
-                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
-                else:
-                    pass
-                    #FIXME: Check and perhaps fix neighbour structure
-
-            else:
-                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
-                msg += 'bound to a boundary object.\n'
-                msg += 'All boundary tags defined in domain must appear '
-                msg += 'in the supplied dictionary.\n'
-                msg += 'The tags are: %s' %self.get_boundary_tags()
-                raise msg
-
-
-
-    def check_integrity(self):
-        #Mesh.check_integrity(self)
-
-        #print self.quantities
-        #print self.conserved_quantities
-        
-        for quantity in self.conserved_quantities:
-            msg = 'Conserved quantities must be a subset of all quantities'
-            assert quantity in self.quantities, msg
-
-        for quantity in self.evolved_quantities:
-            msg = 'Evolved quantities must be a subset of all quantities'
-            assert quantity in self.quantities, msg            
-
-        # #assert hasattr(self, 'boundary_objects')
-
-    def write_time(self):
-        print self.timestepping_statistics()
-
-    def timestepping_statistics(self):
-        """Return string with time stepping statistics for printing or logging
-        """
-
-        msg = ''
-        if self.min_timestep == self.max_timestep:
-            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-                   %(self.time, self.min_timestep, self.number_of_steps,
-                     self.number_of_first_order_steps)
-        elif self.min_timestep > self.max_timestep:
-            msg += 'Time = %.4f, steps=%d (%d)'\
-                   %(self.time, self.number_of_steps,
-                     self.number_of_first_order_steps)
-        else:
-            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-                   %(self.time, self.min_timestep,
-                     self.max_timestep, self.number_of_steps,
-                     self.number_of_first_order_steps)
-
-        return msg
-
-    def get_name(self):
-        return self.filename
-
-    def set_name(self, name):
-        self.filename = name
-
-    def get_datadir(self):
-        return self.datadir
-
-    def set_datadir(self, name):
-        self.datadir = name
-
-    def set_CFL(self, cfl):
-        if cfl > 1.0:
-            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
-        self.CFL = cfl
-
-    def get_CFL(self):
-        return self.CFL
-    
-    def set_filename(self, name):
-        self.filename = name
-
-    def get_filename(self):
-        return self.filename
-
-    def get_limiter(self):
-        return self.limiter
-
-    def set_limiter(self,limiter):
-
-        possible_limiters = \
-        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
-
-        if limiter in possible_limiters:
-            self.limiter = limiter
-            return
-
-        msg = '%s is an incorrect limiter type.\n'% limiter
-        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
-        raise Exception, msg
-
-        
-    #--------------------------
-    # Main components of evolve
-    #--------------------------    
-
-    def evolve(self, yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
-        """Evolve model through time starting from self.starttime.
-
-
-        yieldstep: numpy.interval between yields where results are stored,
-                   statistics written and domain inspected or
-                   possibly modified. If omitted the internal predefined
-                   max timestep is used.
-                   numpy.internally, smaller timesteps may be taken.
-
-        duration: Duration of simulation
-
-        finaltime: Time where simulation should end. This is currently
-        relative time.  So it's the same as duration.
-
-        If both duration and finaltime are given an exception is thrown.
-
-
-        skip_initial_step: Boolean flag that decides whether the first
-        yield step is skipped or not. This is useful for example to avoid
-        duplicate steps when multiple evolve processes are dove tailed.
-
-
-        Evolve is implemented as a generator and is to be called as such, e.g.
-
-        for t in domain.evolve(yieldstep, finaltime):
-            <Do something with domain and t>
-
-
-        All times are given in seconds
-
-        """
-
-        from config import min_timestep, max_timestep, epsilon
-
-        # FIXME: Maybe lump into a larger check prior to evolving
-        msg = 'Boundary tags must be bound to boundary objects before '
-        msg += 'evolving system, '
-        msg += 'e.g. using the method set_boundary.\n'
-        msg += 'This system has the boundary tags %s '\
-               %self.get_boundary_tags()
-        assert hasattr(self, 'boundary_objects'), msg
-
-
-        if yieldstep is None:
-            yieldstep = max_timestep
-        else:
-            yieldstep = numpy.float(yieldstep)
-
-        self._order_ = self.default_order
-
-
-        if finaltime is not None and duration is not None:
-            # print 'F', finaltime, duration
-            msg = 'Only one of finaltime and duration may be specified'
-            raise msg
-        else:
-            if finaltime is not None:
-                self.finaltime = numpy.float(finaltime)
-            if duration is not None:
-                self.finaltime = self.starttime + numpy.float(duration)
-
-
-
-        N = len(self) # Number of triangles
-        self.yieldtime = 0.0 # Track time between 'yields'
-
-        # Initialise interval of timestep sizes (for reporting only)
-        self.min_timestep = max_timestep
-        self.max_timestep = min_timestep
-        self.number_of_steps = 0
-        self.number_of_first_order_steps = 0
-
-
-        # Update ghosts
-        self.update_ghosts()
-
-        # Initial update of vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update extrema if necessary (for reporting)
-        self.update_extrema()
-        
-        # Initial update boundary values
-        self.update_boundary()
-
-        # Or maybe restore from latest checkpoint
-        if self.checkpoint is True:
-            self.goto_latest_checkpoint()
-
-        if skip_initial_step is False:
-            yield(self.time)  # Yield initial values
-
-        while True:
-
-            # Evolve One Step, using appropriate timestepping method
-            if self.get_timestepping_method() == 'euler':
-                self.evolve_one_euler_step(yieldstep,finaltime)
-                
-            elif self.get_timestepping_method() == 'rk2':
-                self.evolve_one_rk2_step(yieldstep,finaltime)
-
-            elif self.get_timestepping_method() == 'rk3':
-                self.evolve_one_rk3_step(yieldstep,finaltime)               
-
-            
-            # Update extrema if necessary (for reporting)
-            self.update_extrema()            
-
-
-            
-            self.yieldtime += self.timestep
-            self.number_of_steps += 1
-            if self._order_ == 1:
-                self.number_of_first_order_steps += 1
-
-
-            # Yield results
-            if finaltime is not None and self.time >= finaltime-epsilon:
-
-                if self.time > finaltime:
-                    # FIXME (Ole, 30 April 2006): Do we need this check?
-                    # Probably not (Ole, 18 September 2008). Now changed to
-                    # Exception
-                    msg = 'WARNING (domain.py): time overshot finaltime. '
-                    msg += 'Contact Ole.Nielsen@ga.gov.au'
-                    raise Exception, msg
-                    
-
-                # Yield final time and stop
-                self.time = finaltime
-                yield(self.time)
-                break
-
-            if self.yieldtime >= yieldstep:
-                # Yield (intermediate) time and allow inspection of domain
-
-                if self.checkpoint is True:
-                    self.store_checkpoint()
-                    self.delete_old_checkpoints()
-
-                # Pass control on to outer loop for more specific actions
-
-                yield(self.time)
-
-                # Reinitialise
-                self.yieldtime = 0.0
-                self.min_timestep = max_timestep
-                self.max_timestep = min_timestep
-                self.number_of_steps = 0
-                self.number_of_first_order_steps = 0
-                #self.max_speed_array = 0.0
-
-
-    def evolve_one_euler_step(self, yieldstep, finaltime):
-        """
-        One Euler Time Step
-        Q^{n+1} = E(h) Q^n
-        """
-
-
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update timestep to fit yieldstep and finaltime
-        self.update_timestep(yieldstep, finaltime)
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-        # Update time
-        self.time += self.timestep
-
-        
-
-
-    def evolve_one_rk2_step(self, yieldstep, finaltime):
-        """
-        One 2nd order RK timestep
-        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
-        """
-
-        # Save initial conserved quantities values
-        self.backup_conserved_quantities()            
-
-        #--------------------------------------
-        # First euler step
-        #--------------------------------------
-
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update timestep to fit yieldstep and finaltime
-        self.update_timestep(yieldstep, finaltime)
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-        # Update time
-        self.time += self.timestep
-
-        #------------------------------------
-        # Second Euler step
-        #------------------------------------
-            
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        #------------------------------------
-        # Combine initial and final values
-        # of conserved quantities and cleanup
-        #------------------------------------
-        
-        # Combine steps
-        self.saxpy_conserved_quantities(0.5, 0.5)
-
-        #-----------------------------------
-        # clean up vertex values
-        #-----------------------------------
- 
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-
-
-    def evolve_one_rk3_step(self, yieldstep, finaltime):
-        """
-        One 3rd order RK timestep
-        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
-        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
-        """
-
-        # Save initial initial conserved quantities values
-        self.backup_conserved_quantities()            
-
-        initial_time = self.time
-        
-        #--------------------------------------
-        # First euler step
-        #--------------------------------------
-
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update timestep to fit yieldstep and finaltime
-        self.update_timestep(yieldstep, finaltime)
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-        # Update time
-        self.time += self.timestep
-
-        #------------------------------------
-        # Second Euler step
-        #------------------------------------
-            
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        #------------------------------------
-        #Combine steps to obtain intermediate
-        #solution at time t^n + 0.5 h
-        #------------------------------------
-
-        # Combine steps
-        self.saxpy_conserved_quantities(0.25, 0.75)
- 
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-        # Set substep time
-        self.time = initial_time + self.timestep*0.5
-
-        #------------------------------------
-        # Third Euler step
-        #------------------------------------
-            
-        # Compute fluxes across each element edge
-        self.compute_fluxes()
-
-        # Update conserved quantities
-        self.update_conserved_quantities()
-
-        #------------------------------------
-        # Combine final and initial values
-        # and cleanup
-        #------------------------------------
-        
-        # Combine steps
-        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
- 
-        # Update ghosts
-        self.update_ghosts()
-
-        # Update vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
-        # Update boundary values
-        self.update_boundary()
-
-        # Set new time
-        self.time = initial_time + self.timestep       
-        
-
-    def backup_conserved_quantities(self):
-        N = len(self) # Number_of_triangles
-
-        # Backup conserved_quantities centroid values
-        for name in self.conserved_quantities:
-            Q = self.quantities[name]
-            Q.backup_centroid_values()        
-
-    def saxpy_conserved_quantities(self,a,b):
-        N = len(self) #number_of_triangles
-
-        # Backup conserved_quantities centroid values
-        for name in self.conserved_quantities:
-            Q = self.quantities[name]
-            Q.saxpy_centroid_values(a,b)        
-
-
-
-        
-    def distribute_to_vertices_and_edges(self):
-        """Extrapolate conserved quantities from centroid to
-        vertices and edge-midpoints for each volume
-
-        Default implementation is straight first order,
-        i.e. constant values throughout each element and
-        no reference to non-conserved quantities.
-        """
-
-        for name in self.conserved_quantities:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-                #Q.limit()
-            else:
-                raise 'Unknown order'
-            #Q.interpolate_from_vertices_to_edges()
-
-
-    def update_ghosts(self):
-        pass
-    
-    def update_boundary(self):
-        """Go through list of boundary objects and update boundary values
-        for all conserved quantities on boundary.
-        """
-
-        #FIXME: Update only those that change (if that can be worked out)
-        #FIXME: Boundary objects should not include ghost nodes.
-        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
-        #    q = B.evaluate(vol_id, edge_id)
-        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
-            q = B.evaluate(vol_id, vertex_id)
-            #print 'q ',q
-            for j, name in enumerate(self.evolved_quantities):
-                #print 'name %s j = %f \n'%(name,j)
-                Q = self.quantities[name]
-                Q.boundary_values[i] = q[j]
-                #print 'Q=',Q
-    
-    def update_timestep(self, yieldstep, finaltime):
-
-        from config import min_timestep, max_timestep
-
-        # self.timestep is calculated from speed of characteristics
-        # Apply CFL condition here
-        timestep = min(self.CFL*self.flux_timestep, max_timestep)
-
-        #Record maximal and minimal values of timestep for reporting
-        self.max_timestep = max(timestep, self.max_timestep)
-        self.min_timestep = min(timestep, self.min_timestep)
-
-        #Protect against degenerate time steps
-        if timestep < min_timestep:
-
-            #Number of consecutive small steps taken b4 taking action
-            self.smallsteps += 1
-
-            if self.smallsteps > self.max_smallsteps:
-                self.smallsteps = 0 #Reset
-
-                if self.order == 1:
-                    msg = 'WARNING: Too small timestep %.16f reached '\
-                          %timestep
-                    msg += 'even after %d steps of 1 order scheme'\
-                           %self.max_smallsteps
-                    print msg
-                    timestep = min_timestep  #Try enforcing min_step
-
-                    #raise msg
-                else:
-                    #Try to overcome situation by switching to 1 order
-                    print "changing Order 1"
-                    self.order = 1
-
-        else:
-            self.smallsteps = 0
-            if self.order == 1 and self.default_order == 2:
-                self.order = 2
-
-
-        #Ensure that final time is not exceeded
-        if finaltime is not None and self.time + timestep > finaltime:
-            timestep = finaltime-self.time
-
-        #Ensure that model time is aligned with yieldsteps
-        if self.yieldtime + timestep > yieldstep:
-            timestep = yieldstep-self.yieldtime
-
-        self.timestep = timestep
-
-    def update_extrema(self):
-        pass
-
-    def compute_forcing_terms(self):
-        """If there are any forcing functions driving the system
-        they should be defined in Domain subclass and appended to
-        the list self.forcing_terms
-        """
-        #Clears explicit_update needed for second order method
-        if self.time_order == 2:
-            for name in self.conserved_quantities:
-                Q = self.quantities[name]
-                Q.explicit_update[:] = 0.0
-
-        for f in self.forcing_terms:
-            f(self)
-            
-
-    def update_derived_quantites(self):
-        pass
-    
-    #def update_conserved_quantities(self):
-    def update_conserved_quantities(self):
-        """Update vectors of conserved quantities using previously
-        computed fluxes specified forcing functions.
-        """
-
-
-        timestep = self.timestep
-
-        #Compute forcing terms
-        self.compute_forcing_terms()
-
-        #Update conserved_quantities
-        for name in self.conserved_quantities:
-            Q = self.quantities[name]
-            Q.update(timestep)
-
-
-
-if __name__ == "__main__":
-
-    points1 = [0.0, 1.0, 2.0, 3.0]
-    D1 = Domain(points1)
-
-    print D1.get_coordinate(0)
-    print D1.get_coordinate(0,1)
-    print 'Number of Elements = ',D1.number_of_elements
-
-    try:
-        print D1.get_coordinate(3)
-    except:
-        pass
-    else:
-        msg =  'Should have raised an out of bounds exception'
-        raise msg
-
-    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
-    #D2 = Domain(points2)
Index: /anuga_work/development/pipeflow/generic_domain.py
===================================================================
--- /anuga_work/development/pipeflow/generic_domain.py	(revision 7788)
+++ /anuga_work/development/pipeflow/generic_domain.py	(revision 7788)
@@ -0,0 +1,1357 @@
+"""Class Domain - 1D domains for finite-volume computations of
+   the shallow water wave equation
+
+
+   Copyright 2004
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia
+"""
+
+import numpy
+from generic_boundary_conditions import *
+
+
+class Generic_Domain:
+
+    def __init__(self,
+                 coordinates,
+                 boundary = None,
+                 conserved_quantities = None,
+                 evolved_quantities = None,
+                 other_quantities = None,
+                 tagged_elements = None):
+        """
+        Build 1D elements from x coordinates
+        """
+
+        from config import timestepping_method
+        from config import CFL
+        
+        #Store Points
+        self.coordinates = numpy.array(coordinates, numpy.float)
+
+
+        #Register number of Elements
+        self.number_of_elements = N = len(self.coordinates)-1
+
+        self.beta = 1.0
+        self.set_limiter("minmod_kurganov")
+        self.set_CFL(CFL)
+        self.set_timestepping_method(timestepping_method)
+
+        self.wet_nodes = numpy.zeros((N,2), numpy.int) # should this be here
+
+        #Allocate space for neighbour and boundary structures
+        self.neighbours = numpy.zeros((N, 2), numpy.int)
+        #self.neighbour_edges = numpy.zeros((N, 2), numpy.int)
+        self.neighbour_vertices   = numpy.zeros((N, 2), numpy.int)
+        self.number_of_boundaries = numpy.zeros(N, numpy.int)
+        self.surrogate_neighbours = numpy.zeros((N, 2), numpy.int)
+        
+        #Allocate space for geometric quantities
+        self.vertices  = numpy.zeros((N, 2), numpy.float)
+        self.centroids = numpy.zeros(N, numpy.float)
+        self.areas     = numpy.zeros(N, numpy.float)
+
+        self.max_speed_array = numpy.zeros(N, numpy.float)
+      
+
+        self.normals = numpy.zeros((N, 2), numpy.float)
+        
+        for i in range(N):
+            xl = self.coordinates[i]
+            xr = self.coordinates[i+1]
+            self.vertices[i,0] = xl
+            self.vertices[i,1] = xr
+
+            centroid = (xl+xr)/2.0
+            self.centroids[i] = centroid
+
+            msg = 'Coordinates should be ordered, smallest to largest'
+            assert xr>xl, msg
+            
+            #The normal vectors
+            # - point outward from each edge
+            # - are orthogonal to the edge
+            # - have unit length
+            # - Are enumerated by left vertex then right vertex normals
+
+            nl = -1.0
+            nr =  1.0
+            self.normals[i,:] = [nl, nr]
+
+            self.areas[i] = (xr-xl)
+            
+# #         print 'N', N
+# #         print 'Centroid', self.centroids
+# #         print 'Areas', self.areas
+# #         print 'Vertex_Coordinates', self.vertices
+            
+            #Initialise Neighbours (-1 means that it is a boundary neighbour)
+            self.neighbours[i, :] = [-1, -1]
+            #Initialise edge ids of neighbours
+            #Initialise vertex ids of neighbours
+            #In case of boundaries this slot is not used
+            #self.neighbour_edges[i, :] = [-1, -1]
+            self.neighbour_vertices[i, :] = [-1, -1]
+
+        self.build_vertexlist()
+
+        #Build neighbour structure
+        self.build_neighbour_structure()
+
+        #Build surrogate neighbour structure
+        self.build_surrogate_neighbour_structure()          
+
+        #Build boundary dictionary mapping (id, edge) to symbolic tags
+        #Build boundary dictionary mapping (id, vertex) to symbolic tags
+        self.build_boundary_dictionary(boundary)
+
+        #Build tagged element  dictionary mapping (tag) to array of elements
+        self.build_tagged_elements_dictionary(tagged_elements)
+        
+        from quantity import Quantity
+        #from quantity_domain import Quantity, Conserved_quantity
+        
+        #List of quantity names entering
+        #the conservation equations
+        #(Must be a subset of quantities)
+        if conserved_quantities is None:
+            self.conserved_quantities = []
+        else:
+            self.conserved_quantities = conserved_quantities
+
+        if evolved_quantities is None:
+            self.evolved_quantities = self.conserved_quantities
+        else:
+            self.evolved_quantities = evolved_quantities
+
+        if other_quantities is None:
+            self.other_quantities = []
+        else:
+            self.other_quantities = other_quantities
+
+
+        #Build dictionary of Quantity instances keyed by quantity names
+        self.quantities = {}
+
+        #print self.conserved_quantities
+        #print self.evolved_quantities
+        
+
+        #FIXME: remove later - maybe OK, though....
+        for name in self.evolved_quantities:
+            self.quantities[name] = Quantity(self)
+        for name in self.other_quantities:
+            self.quantities[name] = Quantity(self)
+
+        #Create an empty list for explicit forcing terms
+        self.forcing_terms = []
+
+        #Defaults
+        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
+        self.beta_w = beta_w
+        self.beta_h = beta_h
+        self.epsilon = epsilon
+
+        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
+        #Or maybe get rid of order altogether and use beta_w and beta_h
+        self.default_order = 1
+        self.order = self.default_order
+
+        self.default_time_order = 1
+        self.time_order = self.default_time_order
+        
+        self.smallsteps = 0
+        self.max_smallsteps = max_smallsteps
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+
+        #Model time
+        self.time = 0.0
+        self.finaltime = None
+        self.min_timestep = self.max_timestep = 0.0
+        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        #Checkpointing and storage
+        from config import default_datadir
+        self.set_datadir(default_datadir)
+        self.filename = 'domain'
+        self.checkpoint = False
+
+    def __len__(self):
+        return self.number_of_elements
+
+    def build_vertexlist(self):
+        """Build vertexlist index by vertex ids and for each entry (point id)
+        build a list of (triangles, vertex_id) pairs that use the point
+        as vertex.
+
+        Preconditions:
+          self.coordinates and self.triangles are defined
+
+        Postcondition:
+          self.vertexlist is built
+        """
+
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(self.number_of_elements):
+
+            #a = self.triangles[i, 0]
+            #b = self.triangles[i, 1]
+            #c = self.triangles[i, 2]
+            a = i
+            b = i + 1
+
+            #Register the vertices v as lists of
+            #(triangle_id, vertex_id) tuples associated with them
+            #This is used for smoothing
+            #for vertex_id, v in enumerate([a,b,c]):
+            for vertex_id, v in enumerate([a,b]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )
+
+        self.vertexlist = vertexlist
+
+        
+    def build_neighbour_structure(self):
+        """Update all registered triangles to point to their neighbours.
+
+        Also, keep a tally of the number of boundaries for each triangle
+
+        Postconditions:
+          neighbours and neighbour_edges is populated
+          neighbours and neighbour_vertices is populated
+          number_of_boundaries integer array is defined.
+        """
+
+        #Step 1:
+        #Build dictionary mapping from segments (2-tuple of points)
+        #to left hand side edge (facing neighbouring triangle)
+
+        N = self.number_of_elements
+        neighbourdict = {}
+        #l_edge = 0
+        #r_edge = 1
+        l_vertex = 0
+        r_vertex = 1
+        for i in range(N):
+
+            #Register all segments as keys mapping to current triangle
+            #and segment id
+            #a = self.triangles[i, 0]
+            #b = self.triangles[i, 1]
+            #c = self.triangles[i, 2]
+            a = self.vertices[i,0]
+            b = self.vertices[i,1]
+            
+            """
+            if neighbourdict.has_key((a,b)):
+                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
+                    raise msg
+            if neighbourdict.has_key((b,c)):
+                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
+                    raise msg
+            if neighbourdict.has_key((c,a)):
+                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
+                    raise msg
+                    """
+            #neighbourdict[a,b] = (i, 2) #(id, edge)
+            #neighbourdict[b,c] = (i, 0) #(id, edge)
+            #neighbourdict[c,a] = (i, 1) #(id, edge)
+            #neighbourdict[a,b] = (i, 1) #(id, edge)
+            #neighbourdict[b,a] = (i, 0) #(id, edge)
+            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
+            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
+            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
+            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
+
+
+        #Step 2:
+        #Go through triangles again, but this time
+        #reverse direction of segments and lookup neighbours.
+        for i in range(N):
+            #a = self.triangles[i, 0]
+            #b = self.triangles[i, 1]
+            #c = self.triangles[i, 2]
+            
+            a = self.vertices[i,0]
+            b = self.vertices[i,1]
+            
+            #self.number_of_boundaries[i] = 3
+            self.number_of_boundaries[i] = 2
+            #if neighbourdict.has_key((b,l_edge)):
+            if neighbourdict.has_key((b,l_vertex)):
+                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
+                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
+                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
+                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
+                self.number_of_boundaries[i] -= 1
+
+            #if neighbourdict.has_key((a,r_edge)):
+            if neighbourdict.has_key((a,r_vertex)):
+                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
+                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
+                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
+                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
+                self.number_of_boundaries[i] -= 1
+                
+            #if neighbourdict.has_key((b,a)):
+            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
+            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
+            #    self.number_of_boundaries[i] -= 1
+
+            #if neighbourdict.has_key((c,b)):
+            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
+            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
+            #    self.number_of_boundaries[i] -= 1
+
+            #if neighbourdict.has_key((a,b)):
+            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
+            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
+            #    self.number_of_boundaries[i] -= 1
+
+    def build_surrogate_neighbour_structure(self):
+        """Build structure where each triangle edge points to its neighbours
+        if they exist. Otherwise point to the triangle itself.
+
+        The surrogate neighbour structure is useful for computing gradients
+        based on centroid values of neighbours.
+
+        Precondition: Neighbour structure is defined
+        Postcondition:
+          Surrogate neighbour structure is defined:
+          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
+          triangles.
+
+        """
+
+        N = self.number_of_elements
+        for i in range(N):
+            #Find all neighbouring volumes that are not boundaries
+            #for k in range(3):
+            for k in range(2):    
+                if self.neighbours[i, k] < 0:
+                    self.surrogate_neighbours[i, k] = i #Point this triangle
+                else:
+                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
+
+    def build_boundary_dictionary(self, boundary = None):
+        """Build or check the dictionary of boundary tags.
+         self.boundary is a dictionary of tags,
+         keyed by volume id and edge:
+         { (id, edge): tag, ... }
+
+         Postconditions:
+            self.boundary is defined.
+        """
+
+        from config import default_boundary_tag
+
+        if boundary is None:
+            boundary = {}
+            for vol_id in range(self.number_of_elements):
+                #for edge_id in range(0, 3):
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:   
+                        #boundary[(vol_id, edge_id)] = default_boundary_tag
+                        boundary[(vol_id, vertex_id)] = default_boundary_tag
+        else:
+            #Check that all keys in given boundary exist
+            #for vol_id, edge_id in boundary.keys():
+            for vol_id, vertex_id in boundary.keys():    
+                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
+                a, b = self.neighbours.shape
+                #assert vol_id < a and edge_id < b, msg
+                assert vol_id < a and vertex_id < b, msg
+
+                #FIXME: This assert violates internal boundaries (delete it)
+                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
+                #assert self.neighbours[vol_id, edge_id] < 0, msg
+
+            #Check that all boundary segments are assigned a tag
+            for vol_id in range(self.number_of_elements):
+                #for edge_id in range(0, 3):
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:    
+                        #if not boundary.has_key( (vol_id, edge_id) ):
+                        if not boundary.has_key( (vol_id, vertex_id) ):    
+                            msg = 'WARNING: Given boundary does not contain '
+                            #msg += 'tags for edge (%d, %d). '\
+                            #       %(vol_id, edge_id)
+                            msg += 'tags for vertex (%d, %d). '\
+                                   %(vol_id, vertex_id)
+                            msg += 'Assigning default tag (%s).'\
+                                   %default_boundary_tag
+
+                            #FIXME: Print only as per verbosity
+                            #print msg
+
+                            #FIXME: Make this situation an error in the future
+                            #and make another function which will
+                            #enable default boundary-tags where
+                            #tags a not specified
+                            #boundary[ (vol_id, edge_id) ] =\
+                            boundary[ (vol_id, vertex_id) ] =\
+                                      default_boundary_tag
+
+
+
+        self.boundary = boundary
+
+    def build_tagged_elements_dictionary(self, tagged_elements = None):
+        """Build the dictionary of element tags.
+         self.tagged_elements is a dictionary of element arrays,
+         keyed by tag:
+         { (tag): [e1, e2, e3..] }
+
+         Postconditions:
+            self.element_tag is defined
+        """
+
+        if tagged_elements is None:
+            tagged_elements = {}
+        else:
+            #Check that all keys in given boundary exist
+            for tag in tagged_elements.keys():
+                tagged_elements[tag] = array(tagged_elements[tag]).astype(numpy.int)
+
+                msg = 'Not all elements exist. '
+                assert max(tagged_elements[tag]) < self.number_of_elements, msg
+        #print "tagged_elements", tagged_elements
+        self.tagged_elements = tagged_elements
+
+
+    def set_quantities_to_be_stored(self, q):
+        """Specify which quantities will be stored in the sww file.
+
+        q must be either:
+          - the name of a quantity
+          - a list of quantity names
+          - None
+
+        In the two first cases, the named quantities will be stored at each
+        yieldstep
+        (This is in addition to the quantities elevation and friction)  
+        If q is None, storage will be switched off altogether.
+        """
+
+
+        if q is None:
+            self.quantities_to_be_stored = []            
+            self.store = False
+            return
+
+        if isinstance(q, basestring):
+            q = [q] # Turn argument into a list
+
+        #Check correcness    
+        for quantity_name in q:
+            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
+            assert quantity_name in self.conserved_quantities, msg 
+        
+        self.quantities_to_be_stored = q
+        
+
+
+
+
+    def get_boundary_tags(self):
+        """Return list of available boundary tags
+        """
+
+        tags = {}
+        for v in self.boundary.values():
+            tags[v] = 1
+
+        return tags.keys()
+        
+    def get_vertex_coordinates(self, obj = False):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+
+        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
+        FIXME, we might make that the default.
+        FIXME Maybe use keyword: continuous = False for this condition?
+
+        
+        """
+
+        if obj is True:
+        
+            #V = self.vertex_coordinates
+            V = self.vertices
+            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
+
+            N = V.shape[0]
+            #return reshape(V, (3*N, 2))
+            return numpy.reshape(V, (N, 2))
+        else:    
+            #return self.vertex_coordinates
+            return self.vertices
+    
+    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
+        """Get conserved quantities at volume vol_id
+
+        If vertex is specified use it as index for vertex values
+        If edge is specified use it as index for edge values
+        If neither are specified use centroid values
+        If both are specified an exeception is raised
+
+        Return value: Vector of length == number_of_conserved quantities
+
+        """
+
+        #if not (vertex is None):# or edge is None):
+        #    msg = 'Values for both vertex and edge was specified.'
+        #    msg += 'Only one (or none) is allowed.'
+       #     raise msg
+
+        q = numpy.zeros( len(self.conserved_quantities), numpy.float)
+
+        for i, name in enumerate(self.conserved_quantities):
+            Q = self.quantities[name]
+            if vertex is not None:
+                q[i] = Q.vertex_values[vol_id, vertex]
+            #elif edge is not None:
+            #    q[i] = Q.edge_values[vol_id, edge]
+            else:
+                q[i] = Q.centroid_values[vol_id]
+
+        return q
+
+
+    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
+        """Get evolved quantities at volume vol_id
+
+        If vertex is specified use it as index for vertex values
+        If edge is specified use it as index for edge values
+        If neither are specified use centroid values
+        If both are specified an exeception is raised
+
+        Return value: Vector of length == number_of_evolved quantities
+
+        """
+
+
+
+        #if not (vertex is None):# or edge is None):
+        #    msg = 'Values for both vertex and edge was specified.'
+        #    msg += 'Only one (or none) is allowed.'
+       #     raise msg
+
+        q = numpy.zeros( len(self.evolved_quantities), numpy.float)
+
+        for i, name in enumerate(self.evolved_quantities):
+            Q = self.quantities[name]
+            if vertex is not None:
+                q[i] = Q.vertex_values[vol_id, vertex]
+            #elif edge is not None:
+            #    q[i] = Q.edge_values[vol_id, edge]
+            else:
+                q[i] = Q.centroid_values[vol_id]
+
+        return q
+        
+
+    def get_centroids(self):
+        """Return all coordinates of centroids
+        Return x coordinate of centroid for each element as a N array
+        """
+
+        return self.centroids
+
+    def get_vertices(self):
+        """Return all coordinates of centroids
+        Return x coordinate of centroid for each element as a N array
+        """
+
+        return self.vertices
+
+    def get_coordinate(self, elem_id, vertex=None):
+        """Return coordinate of centroid,
+        or left or right vertex.
+        Left vertex (vertex=0). Right vertex (vertex=1)
+        """
+
+        if vertex is None:
+            return self.centroids[elem_id]
+        else:
+            return self.vertices[elem_id,vertex]
+
+    def get_area(self, elem_id):
+        """Return area of element id
+        """
+
+        return self.areas[elem_id]
+
+    def get_quantity(self, name, location='vertices', indices = None):
+        """Get values for named quantity
+
+        name: Name of quantity
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a numpyal array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3.
+
+        Indices is the set of element ids that the operation applies to.
+
+        The values will be stored in elements following their
+        internal ordering.
+        """
+
+        return self.quantities[name].get_values( location, indices = indices)
+
+    def get_centroid_coordinates(self):
+        """Return all centroid coordinates.
+        Return all centroid coordinates for all triangles as an Nx2 array
+        (ordered as x0, y0 for each triangle)
+        """
+        return self.centroids
+
+
+    def get_timestepping_method(self):
+        return self.timestepping_method
+
+    def set_timestepping_method(self,timestepping_method):
+        
+        if timestepping_method in ['euler', 'rk2', 'rk3']:
+            self.timestepping_method = timestepping_method
+            return
+
+        msg = '%s is an incorrect timestepping type'% timestepping_method
+        raise Exception, msg
+
+
+    def set_quantity(self, name, *args, **kwargs):
+        """Set values for named quantity
+
+
+        One keyword argument is documented here:
+        expression = None, # Arbitrary expression
+
+        expression:
+          Arbitrary expression involving quantity names
+
+        See Quantity.set_values for further documentation.
+        """
+
+        #FIXME (Ole): Allow new quantities here
+        #from quantity import Quantity, Conserved_quantity
+        #Create appropriate quantity object
+        # #if name in self.conserved_quantities:
+        # #    self.quantities[name] = Conserved_quantity(self)
+        # #else:
+        # #    self.quantities[name] = Quantity(self)
+
+
+        #Do the expression stuff
+        if kwargs.has_key('expression'):
+            expression = kwargs['expression']
+            del kwargs['expression']
+
+            Q = self.create_quantity_from_expression(expression)
+            kwargs['quantity'] = Q
+        
+        #Assign values
+        self.quantities[name].set_values(*args, **kwargs)
+
+    def set_boundary(self, boundary_map):
+        """Associate boundary objects with tagged boundary segments.
+
+        Input boundary_map is a dictionary of boundary objects keyed
+        by symbolic tags to matched against tags in the internal dictionary
+        self.boundary.
+
+        As result one pointer to a boundary object is stored for each vertex
+        in the list self.boundary_objects.
+        More entries may point to the same boundary object
+
+        Schematically the mapping is from two dictionaries to one list
+        where the index is used as pointer to the boundary_values arrays
+        within each quantity.
+
+        self.boundary:          (vol_id, edge_id): tag
+        boundary_map (input):   tag: boundary_object
+        ----------------------------------------------
+        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
+
+
+        Pre-condition:
+          self.boundary has been built.
+
+        Post-condition:
+          self.boundary_objects is built
+
+        If a tag from the domain doesn't appear in the input dictionary an
+        exception is raised.
+        However, if a tag is not used to the domain, no error is thrown.
+        FIXME: This would lead to implementation of a
+        default boundary condition
+
+        Note: If a segment is listed in the boundary dictionary and if it is
+        not None, it *will* become a boundary -
+        even if there is a neighbouring triangle.
+        This would be the case for internal boundaries
+
+        Boundary objects that are None will be skipped.
+
+        FIXME: If set_boundary is called multiple times and if Boundary
+        object is changed into None, the neighbour structure will not be
+        restored!!!
+        """
+
+        self.boundary_objects = []
+        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
+
+        #FIXME: Try to remove the sorting and fix test_mesh.py
+        x = self.boundary.keys()
+        x.sort()
+
+        #Loop through edges that lie on the boundary and associate them with
+        #callable boundary objects depending on their tags
+        #for k, (vol_id, edge_id) in enumerate(x):
+        for k, (vol_id, vertex_id) in enumerate(x):
+            #tag = self.boundary[ (vol_id, edge_id) ]
+            tag = self.boundary[ (vol_id, vertex_id) ]
+
+            if boundary_map.has_key(tag):
+                B = boundary_map[tag]  #Get callable boundary object
+
+                if B is not None:
+                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
+                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
+                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
+                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
+                else:
+                    pass
+                    #FIXME: Check and perhaps fix neighbour structure
+
+            else:
+                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
+                msg += 'bound to a boundary object.\n'
+                msg += 'All boundary tags defined in domain must appear '
+                msg += 'in the supplied dictionary.\n'
+                msg += 'The tags are: %s' %self.get_boundary_tags()
+                raise msg
+
+
+
+    def check_integrity(self):
+        #Mesh.check_integrity(self)
+
+        #print self.quantities
+        #print self.conserved_quantities
+        
+        for quantity in self.conserved_quantities:
+            msg = 'Conserved quantities must be a subset of all quantities'
+            assert quantity in self.quantities, msg
+
+        for quantity in self.evolved_quantities:
+            msg = 'Evolved quantities must be a subset of all quantities'
+            assert quantity in self.quantities, msg            
+
+        # #assert hasattr(self, 'boundary_objects')
+
+    def write_time(self):
+        print self.timestepping_statistics()
+
+    def timestepping_statistics(self):
+        """Return string with time stepping statistics for printing or logging
+        """
+
+        msg = ''
+        if self.min_timestep == self.max_timestep:
+            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
+                   %(self.time, self.min_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
+        elif self.min_timestep > self.max_timestep:
+            msg += 'Time = %.4f, steps=%d (%d)'\
+                   %(self.time, self.number_of_steps,
+                     self.number_of_first_order_steps)
+        else:
+            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
+                   %(self.time, self.min_timestep,
+                     self.max_timestep, self.number_of_steps,
+                     self.number_of_first_order_steps)
+
+        return msg
+
+    def get_name(self):
+        return self.filename
+
+    def set_name(self, name):
+        self.filename = name
+
+    def get_datadir(self):
+        return self.datadir
+
+    def set_datadir(self, name):
+        self.datadir = name
+
+    def set_CFL(self, cfl):
+        if cfl > 1.0:
+            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
+        self.CFL = cfl
+
+    def get_CFL(self):
+        return self.CFL
+    
+    def set_filename(self, name):
+        self.filename = name
+
+    def get_filename(self):
+        return self.filename
+
+    def get_limiter(self):
+        return self.limiter
+
+    def set_limiter(self,limiter):
+
+        possible_limiters = \
+        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
+
+        if limiter in possible_limiters:
+            self.limiter = limiter
+            return
+
+        msg = '%s is an incorrect limiter type.\n'% limiter
+        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
+        raise Exception, msg
+
+        
+    #--------------------------
+    # Main components of evolve
+    #--------------------------    
+
+    def evolve(self, yieldstep = None,
+               finaltime = None,
+               duration = None,
+               skip_initial_step = False):
+        """Evolve model through time starting from self.starttime.
+
+
+        yieldstep: numpy.interval between yields where results are stored,
+                   statistics written and domain inspected or
+                   possibly modified. If omitted the internal predefined
+                   max timestep is used.
+                   numpy.internally, smaller timesteps may be taken.
+
+        duration: Duration of simulation
+
+        finaltime: Time where simulation should end. This is currently
+        relative time.  So it's the same as duration.
+
+        If both duration and finaltime are given an exception is thrown.
+
+
+        skip_initial_step: Boolean flag that decides whether the first
+        yield step is skipped or not. This is useful for example to avoid
+        duplicate steps when multiple evolve processes are dove tailed.
+
+
+        Evolve is implemented as a generator and is to be called as such, e.g.
+
+        for t in domain.evolve(yieldstep, finaltime):
+            <Do something with domain and t>
+
+
+        All times are given in seconds
+
+        """
+
+        from config import min_timestep, max_timestep, epsilon
+
+        # FIXME: Maybe lump into a larger check prior to evolving
+        msg = 'Boundary tags must be bound to boundary objects before '
+        msg += 'evolving system, '
+        msg += 'e.g. using the method set_boundary.\n'
+        msg += 'This system has the boundary tags %s '\
+               %self.get_boundary_tags()
+        assert hasattr(self, 'boundary_objects'), msg
+
+
+        if yieldstep is None:
+            yieldstep = max_timestep
+        else:
+            yieldstep = numpy.float(yieldstep)
+
+        self._order_ = self.default_order
+
+
+        if finaltime is not None and duration is not None:
+            # print 'F', finaltime, duration
+            msg = 'Only one of finaltime and duration may be specified'
+            raise msg
+        else:
+            if finaltime is not None:
+                self.finaltime = numpy.float(finaltime)
+            if duration is not None:
+                self.finaltime = self.starttime + numpy.float(duration)
+
+
+
+        N = len(self) # Number of triangles
+        self.yieldtime = 0.0 # Track time between 'yields'
+
+        # Initialise interval of timestep sizes (for reporting only)
+        self.min_timestep = max_timestep
+        self.max_timestep = min_timestep
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Initial update of vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update extrema if necessary (for reporting)
+        self.update_extrema()
+        
+        # Initial update boundary values
+        self.update_boundary()
+
+        # Or maybe restore from latest checkpoint
+        if self.checkpoint is True:
+            self.goto_latest_checkpoint()
+
+        if skip_initial_step is False:
+            yield(self.time)  # Yield initial values
+
+        while True:
+
+            # Evolve One Step, using appropriate timestepping method
+            if self.get_timestepping_method() == 'euler':
+                self.evolve_one_euler_step(yieldstep,finaltime)
+                
+            elif self.get_timestepping_method() == 'rk2':
+                self.evolve_one_rk2_step(yieldstep,finaltime)
+
+            elif self.get_timestepping_method() == 'rk3':
+                self.evolve_one_rk3_step(yieldstep,finaltime)               
+
+            
+            # Update extrema if necessary (for reporting)
+            self.update_extrema()            
+
+
+            
+            self.yieldtime += self.timestep
+            self.number_of_steps += 1
+            if self._order_ == 1:
+                self.number_of_first_order_steps += 1
+
+
+            # Yield results
+            if finaltime is not None and self.time >= finaltime-epsilon:
+
+                if self.time > finaltime:
+                    # FIXME (Ole, 30 April 2006): Do we need this check?
+                    # Probably not (Ole, 18 September 2008). Now changed to
+                    # Exception
+                    msg = 'WARNING (domain.py): time overshot finaltime. '
+                    msg += 'Contact Ole.Nielsen@ga.gov.au'
+                    raise Exception, msg
+                    
+
+                # Yield final time and stop
+                self.time = finaltime
+                yield(self.time)
+                break
+
+            if self.yieldtime >= yieldstep:
+                # Yield (intermediate) time and allow inspection of domain
+
+                if self.checkpoint is True:
+                    self.store_checkpoint()
+                    self.delete_old_checkpoints()
+
+                # Pass control on to outer loop for more specific actions
+
+                yield(self.time)
+
+                # Reinitialise
+                self.yieldtime = 0.0
+                self.min_timestep = max_timestep
+                self.max_timestep = min_timestep
+                self.number_of_steps = 0
+                self.number_of_first_order_steps = 0
+                #self.max_speed_array = 0.0
+
+
+    def evolve_one_euler_step(self, yieldstep, finaltime):
+        """
+        One Euler Time Step
+        Q^{n+1} = E(h) Q^n
+        """
+
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        
+
+
+    def evolve_one_rk2_step(self, yieldstep, finaltime):
+        """
+        One 2nd order RK timestep
+        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
+        """
+
+        # Save initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        #--------------------------------------
+        # First euler step
+        #--------------------------------------
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        #------------------------------------
+        # Second Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        # Combine initial and final values
+        # of conserved quantities and cleanup
+        #------------------------------------
+        
+        # Combine steps
+        self.saxpy_conserved_quantities(0.5, 0.5)
+
+        #-----------------------------------
+        # clean up vertex values
+        #-----------------------------------
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+
+
+    def evolve_one_rk3_step(self, yieldstep, finaltime):
+        """
+        One 3rd order RK timestep
+        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
+        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
+        """
+
+        # Save initial initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        initial_time = self.time
+        
+        #--------------------------------------
+        # First euler step
+        #--------------------------------------
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        #------------------------------------
+        # Second Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine steps to obtain intermediate
+        #solution at time t^n + 0.5 h
+        #------------------------------------
+
+        # Combine steps
+        self.saxpy_conserved_quantities(0.25, 0.75)
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Set substep time
+        self.time = initial_time + self.timestep*0.5
+
+        #------------------------------------
+        # Third Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        # Combine final and initial values
+        # and cleanup
+        #------------------------------------
+        
+        # Combine steps
+        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Set new time
+        self.time = initial_time + self.timestep       
+        
+
+    def backup_conserved_quantities(self):
+        N = len(self) # Number_of_triangles
+
+        # Backup conserved_quantities centroid values
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.backup_centroid_values()        
+
+    def saxpy_conserved_quantities(self,a,b):
+        N = len(self) #number_of_triangles
+
+        # Backup conserved_quantities centroid values
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.saxpy_centroid_values(a,b)        
+
+
+
+        
+    def distribute_to_vertices_and_edges(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume
+
+        Default implementation is straight first order,
+        i.e. constant values throughout each element and
+        no reference to non-conserved quantities.
+        """
+
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            if self.order == 1:
+                Q.extrapolate_first_order()
+            elif self.order == 2:
+                Q.extrapolate_second_order()
+                #Q.limit()
+            else:
+                raise 'Unknown order'
+            #Q.interpolate_from_vertices_to_edges()
+
+
+    def update_ghosts(self):
+        pass
+    
+    def update_boundary(self):
+        """Go through list of boundary objects and update boundary values
+        for all conserved quantities on boundary.
+        """
+
+        #FIXME: Update only those that change (if that can be worked out)
+        #FIXME: Boundary objects should not include ghost nodes.
+        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+        #    q = B.evaluate(vol_id, edge_id)
+        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
+            q = B.evaluate(vol_id, vertex_id)
+            #print 'q ',q
+            for j, name in enumerate(self.evolved_quantities):
+                #print 'name %s j = %f \n'%(name,j)
+                Q = self.quantities[name]
+                Q.boundary_values[i] = q[j]
+                #print 'Q=',Q
+    
+    def update_timestep(self, yieldstep, finaltime):
+
+        from config import min_timestep, max_timestep
+
+        # self.timestep is calculated from speed of characteristics
+        # Apply CFL condition here
+        timestep = min(self.CFL*self.flux_timestep, max_timestep)
+
+        #Record maximal and minimal values of timestep for reporting
+        self.max_timestep = max(timestep, self.max_timestep)
+        self.min_timestep = min(timestep, self.min_timestep)
+
+        #Protect against degenerate time steps
+        if timestep < min_timestep:
+
+            #Number of consecutive small steps taken b4 taking action
+            self.smallsteps += 1
+
+            if self.smallsteps > self.max_smallsteps:
+                self.smallsteps = 0 #Reset
+
+                if self.order == 1:
+                    msg = 'WARNING: Too small timestep %.16f reached '\
+                          %timestep
+                    msg += 'even after %d steps of 1 order scheme'\
+                           %self.max_smallsteps
+                    print msg
+                    timestep = min_timestep  #Try enforcing min_step
+
+                    #raise msg
+                else:
+                    #Try to overcome situation by switching to 1 order
+                    print "changing Order 1"
+                    self.order = 1
+
+        else:
+            self.smallsteps = 0
+            if self.order == 1 and self.default_order == 2:
+                self.order = 2
+
+
+        #Ensure that final time is not exceeded
+        if finaltime is not None and self.time + timestep > finaltime:
+            timestep = finaltime-self.time
+
+        #Ensure that model time is aligned with yieldsteps
+        if self.yieldtime + timestep > yieldstep:
+            timestep = yieldstep-self.yieldtime
+
+        self.timestep = timestep
+
+    def update_extrema(self):
+        pass
+
+    def compute_forcing_terms(self):
+        """If there are any forcing functions driving the system
+        they should be defined in Domain subclass and appended to
+        the list self.forcing_terms
+        """
+        #Clears explicit_update needed for second order method
+        if self.time_order == 2:
+            for name in self.conserved_quantities:
+                Q = self.quantities[name]
+                Q.explicit_update[:] = 0.0
+
+        for f in self.forcing_terms:
+            f(self)
+            
+
+    def update_derived_quantites(self):
+        pass
+    
+    #def update_conserved_quantities(self):
+    def update_conserved_quantities(self):
+        """Update vectors of conserved quantities using previously
+        computed fluxes specified forcing functions.
+        """
+
+
+        timestep = self.timestep
+
+        #Compute forcing terms
+        self.compute_forcing_terms()
+
+        #Update conserved_quantities
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.update(timestep)
+
+
+
+if __name__ == "__main__":
+
+    points1 = [0.0, 1.0, 2.0, 3.0]
+    D1 = Domain(points1)
+
+    print D1.get_coordinate(0)
+    print D1.get_coordinate(0,1)
+    print 'Number of Elements = ',D1.number_of_elements
+
+    try:
+        print D1.get_coordinate(3)
+    except:
+        pass
+    else:
+        msg =  'Should have raised an out of bounds exception'
+        raise msg
+
+    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
+    #D2 = Domain(points2)
Index: uga_work/development/pipeflow/shallow_water_domain.py
===================================================================
--- /anuga_work/development/pipeflow/shallow_water_domain.py	(revision 7787)
+++ 	(revision )
@@ -1,1391 +1,0 @@
-"""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
-"""
-
-import numpy
-
-from domain import *
-Generic_Domain = Domain #Rename
-
-#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 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"
-
-        #Realtime visualisation
-        self.visualiser = None
-        self.visualise  = False
-        self.visualise_color_stage = False
-        self.visualise_stage_range = 1.0
-        self.visualise_timer = True
-        self.visualise_range_z = None
-        
-        #Stored output
-        self.store = True
-        self.format = 'sww'
-        self.smooth = True
-
-        #Evolve parametrs
-        self.cfl = 1.0
-        
-        #Reduction operation for get_vertex_values
-        from util import mean
-        self.reduction = mean
-        #self.reduction = min  #Looks better near steep slopes
-
-        self.quantities_to_be_stored = ['stage','xmomentum']
-
-	self.__doc__ = 'shallow_water_domain'
-
-
-    def set_quantities_to_be_stored(self, q):
-        """Specify which quantities will be stored in the sww file.
-
-        q must be either:
-          - the name of a quantity
-          - a list of quantity names
-          - None
-
-        In the two first cases, the named quantities will be stored at each
-        yieldstep
-        (This is in addition to the quantities elevation and friction)  
-        If q is None, storage will be switched off altogether.
-        """
-
-
-        if q is None:
-            self.quantities_to_be_stored = []            
-            self.store = False
-            return
-
-        if isinstance(q, basestring):
-            q = [q] # Turn argument into a list
-
-        #Check correcness    
-        for quantity_name in q:
-            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
-            assert quantity_name in self.conserved_quantities, msg 
-        
-        self.quantities_to_be_stored = q
-        
-
-    def initialise_visualiser(self,scale_z=1.0,rect=None):
-        #Realtime visualisation
-        if self.visualiser is None:
-            from realtime_visualisation_new import Visualiser
-            self.visualiser = Visualiser(self,scale_z,rect)
-            self.visualiser.setup['elevation']=True
-            self.visualiser.updating['stage']=True
-        self.visualise = True
-        if self.visualise_color_stage == True:
-            self.visualiser.coloring['stage'] = True
-            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
-        print 'initialise visualiser'
-        print self.visualiser.setup
-        print self.visualiser.updating
-
-    def check_integrity(self):
-        Generic_Domain.check_integrity(self)
-        #Check that we are solving the shallow water wave equation
-
-        msg = 'First conserved quantity must be "stage"'
-        assert self.conserved_quantities[0] == 'stage', msg
-        msg = 'Second conserved quantity must be "xmomentum"'
-        assert self.conserved_quantities[1] == 'xmomentum', msg
-
-    def extrapolate_second_order_sw(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        extrapolate_second_order_sw(self)
-
-    def compute_fluxes(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        compute_fluxes_C(self)
-
-    def compute_timestep(self):
-        #Call correct module function
-        compute_timestep(self)
-
-    def distribute_to_vertices_and_edges(self):
-        #Call correct module function
-        #(either from this module or C-extension)
-        distribute_to_vertices_and_edges(self)
-        
-    def evolve(self, yieldstep = None, finaltime = None, duration = None,
-               skip_initial_step = False):
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call basic machinery from parent class
-        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
-                                       skip_initial_step):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)
-
-    def initialise_storage(self):
-        """Create and initialise self.writer object for storing data.
-        Also, save x and bed elevation
-        """
-
-        import data_manager
-
-        #Initialise writer
-        self.writer = data_manager.get_dataobject(self, mode = 'w')
-
-        #Store vertices and connectivity
-        self.writer.store_connectivity()
-
-
-    def store_timestep(self, name):
-        """Store named quantity and time.
-
-        Precondition:
-           self.write has been initialised
-        """
-        self.writer.store_timestep(name)
-
-
-#=============== End of Shallow Water Domain ===============================
-
-#Rotation of momentum vector
-def rotate(q, normal, direction = 1):
-    """Rotate the momentum component q (q[1], q[2])
-    from x,y coordinates to coordinates based on normal vector.
-
-    If direction is negative the rotation is inverted.
-
-    Input vector is preserved
-
-    This function is specific to the shallow water wave equation
-    """
-
-
-    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 numpy array'
-
-    assert l == 2, 'Normal vector must have 2 components'
-
-
-    n1 = normal[0]
-    n2 = normal[1]
-
-    r = numpy.zeros(len(q), numpy.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
-
-
-def flux_function(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numpyal vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from config import g, epsilon, h0
-    from math import sqrt
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-	
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #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 +  h0/h_left)
-
-
-    uh_left = u_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 + h0/h_right)
-
-    uh_right = u_right*h_right
-	
-
-    #We have got   h   and    u   at vertex, then  the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    
-    #Flux computation
-    flux_left  = numpy.array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left*h_left], numpy.float)
-    flux_right = numpy.array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right*h_right], numpy.float)
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = array([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[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-# 
-def compute_fluxes_python(domain):
-    """
-	Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    
-    """
-
-    import sys
-
-
-    domain.distribute_to_vertices_and_edges()
-    domain.update_boundary()
-    
-    N = domain.number_of_elements
-    Stage = domain.quantities['stage']
-    Xmom  = domain.quantities['xmomentum']
-    Bed   = domain.quantities['elevation']
-
-    stage = Stage.vertex_values
-    xmom  = Xmom.vertex_values
-    bed   = Bed.vertex_values
-    
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-
-
-	
-    flux = numpy.zeros(2, numpy.float) #Work array for summing up fluxes
-    ql   = numpy.zeros(2, numpy.float)
-    qr   = numpy.zeros(2, numpy.float)
-
-    #Loop
-    timestep = numpy.float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work array
-        #for i in range(3):
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            #ql[0] = stage[k, i]
-            #ql[1] = xmom[k, i]
-            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                qr[0] = stage[n, m]
-                qr[1] = xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-        
-            #Flux computation using provided function
-     
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-            #print 'edgeflux', edgeflux
-
-            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
-            # flux = edgefluxleft - edgefluxright
-            flux -= edgeflux #* domain.edgelengths[k,i]
-            #Update optimal_timestep
-            try:
-                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.cfl*0.5*domain.areas[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]
-        #print "flux cell",k,flux[0]
-
-    domain.flux_timestep = timestep
-    #print domain.quantities['stage'].centroid_values
-#
-def compute_timestep(domain):
-    import sys
-
-
-    N = domain.number_of_elements
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Bed = domain.quantities['elevation']
-    
-    stage = Stage.vertex_values
-    xmom =  Xmom.vertex_values
-    bed =   Bed.vertex_values
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-
-    flux = zeros(2, numpy.float) #Work array for summing up fluxes
-    ql = zeros(2, numpy.float)
-    qr = zeros(2, numpy.float)
-
-    #Loop
-    timestep = numpy.float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work array
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                qr[0] = stage[n, m]
-                qr[1] =  xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-            #Update optimal_timestep
-            try:
-                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
-            except ZeroDivisionError:
-                pass
-
-    domain.timestep = timestep
-
-# Compute flux definition
-def compute_fluxes_C_long(domain):
-
-    import sys
-
-	
-    timestep = numpy.float(sys.maxint)
-    #print 'timestep=',timestep
-    #print 'The type of timestep is',type(timestep)
-	
-    epsilon = domain.epsilon
-    #print 'epsilon=',epsilon
-    #print 'The type of epsilon is',type(epsilon)
-	
-    g = domain.g
-    #print 'g=',g
-    #print 'The type of g is',type(g)
-	
-    neighbours = domain.neighbours
-    #print 'neighbours=',neighbours
-    #print 'The type of neighbours is',type(neighbours)
-	
-    neighbour_vertices = domain.neighbour_vertices
-    #print 'neighbour_vertices=',neighbour_vertices
-    #print 'The type of neighbour_vertices is',type(neighbour_vertices)
-	
-    normals = domain.normals
-    #print 'normals=',normals
-    #print 'The type of normals is',type(normals)
-	
-    areas = domain.areas
-    #print 'areas=',areas
-    #print 'The type of areas is',type(areas)
-	
-    stage_edge_values = domain.quantities['stage'].vertex_values
-    #print 'stage_edge_values=',stage_edge_values
-    #print 'The type of stage_edge_values is',type(stage_edge_values)
-	
-    xmom_edge_values = domain.quantities['xmomentum'].vertex_values
-    #print 'xmom_edge_values=',xmom_edge_values
-    #print 'The type of xmom_edge_values is',type(xmom_edge_values)
-	
-    bed_edge_values = domain.quantities['elevation'].vertex_values
-    #print 'bed_edge_values=',bed_edge_values
-    #print 'The type of bed_edge_values is',type(bed_edge_values)
-	
-    stage_boundary_values = domain.quantities['stage'].boundary_values
-    #print 'stage_boundary_values=',stage_boundary_values
-    #print 'The type of stage_boundary_values is',type(stage_boundary_values)
-	
-    xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
-    #print 'xmom_boundary_values=',xmom_boundary_values
-    #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
-	
-    stage_explicit_update = domain.quantities['stage'].explicit_update
-    #print 'stage_explicit_update=',stage_explicit_update
-    #print 'The type of stage_explicit_update is',type(stage_explicit_update)
-	
-    xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
-    #print 'xmom_explicit_update=',xmom_explicit_update
-    #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
-	
-    number_of_elements = len(stage_edge_values)
-    #print 'number_of_elements=',number_of_elements
-    #print 'The type of number_of_elements is',type(number_of_elements)
-	
-    max_speed_array = domain.max_speed_array
-    #print 'max_speed_array=',max_speed_array
-    #print 'The type of max_speed_array is',type(max_speed_array)
-        
-        
-    from comp_flux_ext import compute_fluxes_ext
-	
-    domain.flux_timestep = compute_fluxes_ext(timestep,
-                                  epsilon,
-                                  g,
-                                  neighbours,
-                                  neighbour_vertices,
-                                  normals,
-                                  areas, 
-                                  stage_edge_values,
-                                  xmom_edge_values,
-                                  bed_edge_values,
-                                  stage_boundary_values,
-                                  xmom_boundary_values,
-                                  stage_explicit_update,
-                                  xmom_explicit_update,
-                                  number_of_elements,
-                                  max_speed_array)
-
-
-# Compute flux definition
-def compute_fluxes_C(domain):
-    import sys
-
-	
-    timestep = numpy.float(sys.maxint)
-
-    stage    = domain.quantities['stage']
-    xmom     = domain.quantities['xmomentum']
-    bed      = domain.quantities['elevation']
-    height   = domain.quantities['height']
-    velocity = domain.quantities['velocity']
-
-
-    from comp_flux_ext import compute_fluxes_ext_short
-
-    domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed,height,velocity)
-
-
-
-# ###################################
-def compute_fluxes_C_wellbalanced(domain):
-    #from numpy import zeros, numpy.float
-    #import sys
-
-	
-    #timestep = numpy.float(sys.maxint)
-    #epsilon = domain.epsilon
-    #g = domain.g
-    #neighbours = domain.neighbours
-    #neighbour_vertices = domain.neighbour_vertices
-    #normals = domain.normals
-    #areas = domain.areas
-    #stage_edge_values = domain.quantities['stage'].vertex_values
-    #xmom_edge_values = domain.quantities['xmomentum'].vertex_values
-    #bed_edge_values = domain.quantities['elevation'].vertex_values
-    #stage_boundary_values = domain.quantities['stage'].boundary_values
-    #xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
-    #stage_explicit_update = domain.quantities['stage'].explicit_update
-    #xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
-    #number_of_elements = len(stage_edge_values)
-    #max_speed_array = domain.max_speed_array
-
-    import sys
-    
-    N = domain.number_of_elements
-    timestep           = numpy.float(sys.maxint)
-    epsilon            = domain.epsilon
-    g                  = domain.g
-    neighbours         = domain.neighbours
-    neighbour_vertices = domain.neighbour_vertices
-    normals            = domain.normals
-    areas              = domain.areas
-    
-    Stage = domain.quantities['stage']
-    Xmom  = domain.quantities['xmomentum']
-    Bed   = domain.quantities['elevation']
-
-    stage_boundary_values  = Stage.boundary_values
-    xmom_boundary_values   = Xmom.boundary_values
-    stage_explicit_update  = Stage.explicit_update
-    xmom_explicit_update   = Xmom.explicit_update
-    max_speed_array        = domain.max_speed_array
-    
-    domain.distribute_to_vertices_and_edges()
-    domain.update_boundary()
-
-
-    stage_V = Stage.vertex_values
-    xmom_V  = Xmom.vertex_values
-    bed_V   = Bed.vertex_values
-    #h_V     = Height.vertex_values
-    #u_V     = Velocity.vertex_values
-
-    number_of_elements = len(stage_V)
-
-    #flux = zeros(2, numpy.float) #Work array for summing up fluxes
-    #ql = zeros(2, numpy.float)
-    #qr = zeros(2, numpy.float)
-	        
-    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced 
-	
-    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
-                                  epsilon,
-                                  g,
-                                  neighbours,
-                                  neighbour_vertices,
-                                  normals,
-                                  areas, 
-                                  stage_V,
-                                  xmom_V,
-                                  bed_V,
-                                  stage_boundary_values,
-                                  xmom_boundary_values,
-                                  stage_explicit_update,
-                                  xmom_explicit_update,
-                                  number_of_elements,
-                                  max_speed_array)
-
-# ###################################
-
-
-
-
-
-
-# Module functions for gradient limiting (distribute_to_vertices_and_edges)
-
-def distribute_to_vertices_and_edges(domain):
-    """Distribution from centroids to vertices specific to the
-    shallow water wave
-    equation.
-
-    It will ensure that h (w-z) is always non-negative even in the
-    presence of steep bed-slopes by taking a weighted average between shallow
-    and deep cases.
-
-    In addition, all conserved quantities get distributed as per either a
-    constant (order==1) or a piecewise linear function (order==2).
-
-    FIXME: more explanation about removal of artificial variability etc
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
-
-    """
-
-    #from config import optimised_gradient_limiter
-
-    #Remove very thin layers of water
-    #protect_against_infinitesimal_and_negative_heights(domain)  
-
-    import sys
-    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)
-    #FIXME replace with numpy.where
-    for i in range(N):
-	h_C[i] = w_C[i] - z_C[i]						
-        if h_C[i] <= 0.0:
-	    #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
-	    h_C[i] = 1.0e-15
-	    w_C[i] = z_C[i]
-	    uh_C[i] = 0.0
-		
-    
-##    for i in range(len(h_C)):
-##	if h_C[i] < epsilon:
-##	    u_C[i] = 0.0  #Could have been negative
-##	    h_C[i] = 0.0
-##	else:
-
-    u_C[:]  = uh_C/(h_C + h0/h_C)
-	
-    for name in [ 'velocity', 'stage' ]:
-        Q = domain.quantities[name]
-        if domain.order == 1:
-            Q.extrapolate_first_order()
-        elif domain.order == 2:
-            #print "add extrapolate second order to shallow water"
-            #if name != 'height':
-            Q.extrapolate_second_order()
-            #Q.limit()
-        else:
-            raise 'Unknown order'
-
-    stage_V = domain.quantities['stage'].vertex_values			
-    bed_V   = domain.quantities['elevation'].vertex_values	
-    h_V     = domain.quantities['height'].vertex_values
-    u_V     = domain.quantities['velocity'].vertex_values
-    xmom_V  = domain.quantities['xmomentum'].vertex_values
-
-    h_V[:]    = stage_V - bed_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
-                stage_V[i,j] = bed_V[i,j]
-                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
-                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
-                
-    xmom_V[:] = u_V * h_V
-    
-    return
-#
-
-
-
-
-
-
-
-#
-def protect_against_infinitesimal_and_negative_heights(domain):
-    """Protect against infinitesimal heights and associated high velocities
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    hc = wc - zc  #Water depths at centroids
-
-    zv = domain.quantities['elevation'].vertex_values
-    wv = domain.quantities['stage'].vertex_values
-    hv = wv-zv
-    xmomv = domain.quantities['xmomentum'].vertex_values
-    #remove the above two lines and corresponding code below
-
-    #Update
-    #FIXME replace with numpy.where
-    for k in range(domain.number_of_elements):
-
-        if hc[k] < domain.minimum_allowed_height:
-            #Control stage
-            if hc[k] < domain.epsilon:
-                wc[k] = zc[k] # Contain 'lost mass' error
-                wv[k,0] = zv[k,0]
-                wv[k,1] = zv[k,1]
-                
-            xmomc[k] = 0.0
-            
-        #N = domain.number_of_elements
-        #if (k == 0) | (k==N-1):
-        #    wc[k] = zc[k] # Contain 'lost mass' error
-        #    wv[k,0] = zv[k,0]
-        #    wv[k,1] = zv[k,1]
-    
-def h_limiter(domain):
-    """Limit slopes for each volume to eliminate artificial variance
-    introduced by e.g. second order extrapolator
-
-    limit on h = w-z
-
-    This limiter depends on two quantities (w,z) so it resides within
-    this module rather than within quantity.py
-    """
-
-    N = domain.number_of_elements
-    beta_h = domain.beta_h
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-    hvbar = zeros(hv.shape, numpy.float) #h-limited values
-
-    #Find min and max of this and neighbour's centroid values
-    hmax = zeros(hc.shape, numpy.float)
-    hmin = zeros(hc.shape, numpy.float)
-
-    for k in range(N):
-        hmax[k] = hmin[k] = hc[k]
-        #for i in range(3):
-        for i in range(2):   
-            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, numpy.float)
-    #for i in range(3):
-    for i in range(2):
-        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):
-        for i in range(2):
-            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):
-        for i in range(2):
-            hvbar[k,i] = hc[k] + phi*dh[k,i]
-
-    return hvbar
-
-def balance_deep_and_shallow(domain):
-    """Compute linear combination between stage as computed by
-    gradient-limiters limiting using w, and stage computed by
-    gradient-limiters limiting using h (h-limiter).
-    The former takes precedence when heights are large compared to the
-    bed slope while the latter takes precedence when heights are
-    relatively small.  Anything in between is computed as a balanced
-    linear combination in order to avoid numpyal disturbances which
-    would otherwise appear as a result of hard switching between
-    modes.
-
-    The h-limiter is always applied irrespective of the order.
-    """
-
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-    #Limit h
-    hvbar = h_limiter(domain)
-
-    for k in range(domain.number_of_elements):
-        #Compute maximal variation in bed elevation
-        #  This quantitiy is
-        #    dz = max_i abs(z_i - z_c)
-        #  and it is independent of dimension
-        #  In the 1d case zc = (z0+z1)/2
-        #  In the 2d case zc = (z0+z1+z2)/3
-
-        dz = max(abs(zv[k,0]-zc[k]),
-                 abs(zv[k,1]-zc[k]))#,
-#                 abs(zv[k,2]-zc[k]))
-
-
-        hmin = min( hv[k,:] )
-
-        #Create alpha in [0,1], where alpha==0 means using the h-limited
-        #stage and alpha==1 means using the w-limited stage as
-        #computed by the gradient limiter (both 1st or 2nd order)
-
-        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
-        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-
-        if dz > 0.0:
-            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
-        else:
-            #Flat bed
-            alpha = 1.0
-
-        alpha  = 0.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):
-            for i in range(2):
-                wv[k,i] = zv[k,i] + (1.0-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.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
-#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-
-
-###############################################
-#Boundaries - specific to the shallow water wave equation
-class Reflective_boundary(Boundary):
-    """Reflective boundary returns same conserved quantities as
-    those present in its neighbour volume but reflected.
-
-    This class is specific to the shallow water equation as it
-    works with the momentum quantities assumed to be the second
-    and third conserved quantities.
-    """
-
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for reflective boundary'
-            raise msg
-
-        #Handy shorthands
-        self.normals  = domain.normals
-        self.stage    = domain.quantities['stage'].vertex_values
-        self.xmom     = domain.quantities['xmomentum'].vertex_values
-        self.bed      = domain.quantities['elevation'].vertex_values
-        self.height   = domain.quantities['height'].vertex_values
-        self.velocity = domain.quantities['velocity'].vertex_values
-
-
-        self.quantities = numpy.zeros(5, numpy.float)
-
-    def __repr__(self):
-        return 'Reflective_boundary'
-
-
-    def evaluate(self, vol_id, edge_id):
-        """Reflective boundaries reverses the outward momentum
-        of the volume they serve.
-        """
-
-        q = self.quantities
-        q[0] =  self.stage[vol_id, edge_id]
-        q[1] = -self.xmom[vol_id, edge_id]
-        q[2] =  self.bed[vol_id, edge_id]
-        q[3] =  self.height[vol_id, edge_id]
-        q[4] = -self.velocity[vol_id, edge_id]
-
-        #normal = self.normals[vol_id,edge_id]
-        
-        return q
-
-class Dirichlet_boundary(Boundary):
-    """Dirichlet boundary returns constant values for the
-    conserved quantities
-    """
-
-
-    def __init__(self, quantities=None):
-        Boundary.__init__(self)
-
-        if quantities is None:
-            msg = 'Must specify one value for each evolved quantity, w,uh,z,h,u'
-            raise msg
-
-        self.quantities=numpy.array(quantities, numpy.float)
-
-    def __repr__(self):
-        return 'Dirichlet boundary (%s)' %self.quantities
-
-    def evaluate(self, vol_id=None, edge_id=None):
-        return self.quantities
-
-
-#########################
-#Standard forcing terms:
-#
-def gravity(domain):
-    """Apply gravitational pull in the presence of bed slope
-    """
-
-    from util import gradient
-
-    xmom  = domain.quantities['xmomentum'].explicit_update
-    stage = domain.quantities['stage'].explicit_update
-#    ymom = domain.quantities['ymomentum'].explicit_update
-
-    Stage = domain.quantities['stage']
-    Elevation = domain.quantities['elevation']
-    #h = Stage.edge_values - Elevation.edge_values
-    h = Stage.vertex_values - Elevation.vertex_values
-    b = Elevation.vertex_values
-    w = Stage.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-
-    for k in range(domain.number_of_elements):
-#        avg_h = sum( h[k,:] )/3
-        avg_h = sum( h[k,:] )/2
-
-        #Compute bed slope
-        #x0, y0, x1, y1, x2, y2 = x[k,:]
-        x0, x1 = x[k,:]
-        #z0, z1, z2 = v[k,:]
-        b0, b1 = b[k,:]
-
-        w0, w1 = w[k,:]
-        wx = gradient(x0, x1, w0, w1)
-
-        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
-        bx = gradient(x0, x1, b0, b1)
-        
-        #Update momentum (explicit update is reset to source values)
-        xmom[k] += -g*bx*avg_h
-        #xmom[k] = -g*bx*avg_h
-        #stage[k] = 0.0
- 
- 
-def manning_friction(domain):
-    """Apply (Manning) friction to water momentum
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-    #vh = domain.quantities['ymomentum'].centroid_values
-    eta = domain.quantities['friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g
-
-    for k in range(N):
-        if eta[k] >= eps:
-            if h[k] >= eps:
-            	#S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
-                S = -g * eta[k]**2 * uh[k]
-            	S /= h[k]**(7.0/3)
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
-            	#ymom_update[k] += S*vh[k]
-
-def linear_friction(domain):
-    """Apply linear friction to water momentum
-
-    Assumes quantity: 'linear_friction' to be present
-    """
-
-    from math import sqrt
-
-    w = domain.quantities['stage'].centroid_values
-    z = domain.quantities['elevation'].centroid_values
-    h = w-z
-
-    uh = domain.quantities['xmomentum'].centroid_values
-#    vh = domain.quantities['ymomentum'].centroid_values
-    tau = domain.quantities['linear_friction'].centroid_values
-
-    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
-
-    N = domain.number_of_elements
-    eps = domain.minimum_allowed_height
-    g = domain.g #Not necessary? Why was this added?
-
-    for k in range(N):
-        if tau[k] >= eps:
-            if h[k] >= eps:
-            	S = -tau[k]/h[k]
-
-            	#Update momentum
-            	xmom_update[k] += S*uh[k]
- #           	ymom_update[k] += S*vh[k]
-
-
-
-def check_forcefield(f):
-    """Check that f is either
-    1: a callable object f(t,x,y), where x and y are vectors
-       and that it returns an array or a list of same length
-       as x and y
-    2: a scalar
-    """
-
-
-    if callable(f):
-        #N = 3
-        N = 2
-        #x = ones(3, numpy.float)
-        #y = ones(3, numpy.float)
-        x = ones(2, numpy.float)
-        #y = ones(2, numpy.float)
-        
-        try:
-            #q = f(1.0, x=x, y=y)
-            q = f(1.0, x=x)
-        except Exception, e:
-            msg = 'Function %s could not be executed:\n%s' %(f, e)
-	    #FIXME: Reconsider this semantics
-            raise msg
-
-        try:
-            q = numpy.array(q, numpy.float)
-        except:
-            msg = 'Return value from vector function %s could ' %f
-            msg += 'not be converted into a numpy array of numpy.floats.\n'
-            msg += 'Specified function should return either list or array.'
-            raise msg
-
-        #Is this really what we want?
-        msg = 'Return vector from function %s ' %f
-        msg += 'must have same lenght as input vectors'
-        assert len(q) == N, msg
-
-    else:
-        try:
-            f = numpy.float(f)
-        except:
-            msg = 'Force field %s must be either a scalar' %f
-            msg += ' or a vector function'
-            raise msg
-    return f
-
-class Wind_stress:
-    """Apply wind stress to water momentum in terms of
-    wind speed [m/s] and wind direction [degrees]
-    """
-
-    def __init__(self, *args, **kwargs):
-        """Initialise windfield from wind speed s [m/s]
-        and wind direction phi [degrees]
-
-        Inputs v and phi can be either scalars or Python functions, e.g.
-
-        W = Wind_stress(10, 178)
-
-        #FIXME - 'normal' degrees are assumed for now, i.e. the
-        vector (1,0) has zero degrees.
-        We may need to convert from 'compass' degrees later on and also
-        map from True north to grid north.
-
-        Arguments can also be Python functions of t,x,y as in
-
-        def speed(t,x,y):
-            ...
-            return s
-
-        def angle(t,x,y):
-            ...
-            return phi
-
-        where x and y are vectors.
-
-        and then pass the functions in
-
-        W = Wind_stress(speed, angle)
-
-        The instantiated object W can be appended to the list of
-        forcing_terms as in
-
-        Alternatively, one vector valued function for (speed, angle)
-        can be applied, providing both quantities simultaneously.
-        As in
-        W = Wind_stress(F), where returns (speed, angle) for each t.
-
-        domain.forcing_terms.append(W)
-        """
-
-        from config import rho_a, rho_w, eta_w
-
-        if len(args) == 2:
-            s = args[0]
-            phi = args[1]
-        elif len(args) == 1:
-            #Assume vector function returning (s, phi)(t,x,y)
-            vector_function = args[0]
-            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
-            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
-            s = lambda t,x: vector_function(t,x=x)[0]
-            phi = lambda t,x: vector_function(t,x=x)[1]
-        else:
-           #Assume info is in 2 keyword arguments
-
-           if len(kwargs) == 2:
-               s = kwargs['s']
-               phi = kwargs['phi']
-           else:
-               raise 'Assumes two keyword arguments: s=..., phi=....'
-
-        print 'phi', phi
-        self.speed = check_forcefield(s)
-        self.phi = check_forcefield(phi)
-
-        self.const = eta_w*rho_a/rho_w
-
-
-    def __call__(self, domain):
-        """Evaluate windfield based on values found in domain
-        """
-
-        from math import pi, cos, sin, sqrt
-
-        xmom_update = domain.quantities['xmomentum'].explicit_update
-        #ymom_update = domain.quantities['ymomentum'].explicit_update
-
-        N = domain.number_of_elements
-        t = domain.time
-
-        if callable(self.speed):
-            xc = domain.get_centroid_coordinates()
-            #s_vec = self.speed(t, xc[:,0], xc[:,1])
-            s_vec = self.speed(t, xc)
-        else:
-            #Assume s is a scalar
-
-            try:
-                s_vec = self.speed * ones(N, numpy.float)
-            except:
-                msg = 'Speed must be either callable or a scalar: %s' %self.s
-                raise msg
-
-
-        if callable(self.phi):
-            xc = domain.get_centroid_coordinates()
-            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
-            phi_vec = self.phi(t, xc)
-        else:
-            #Assume phi is a scalar
-
-            try:
-                phi_vec = self.phi * ones(N, numpy.float)
-            except:
-                msg = 'Angle must be either callable or a scalar: %s' %self.phi
-                raise msg
-
-        #assign_windfield_values(xmom_update, ymom_update,
-        #                        s_vec, phi_vec, self.const)
-        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
-
-
-#def assign_windfield_values(xmom_update, ymom_update,
-#                            s_vec, phi_vec, const):
-def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
-    """Python version of assigning wind field to update vectors.
-    A c version also exists (for speed)
-    """
-    from math import pi, cos, sin, sqrt
-
-    N = len(s_vec)
-    for k in range(N):
-        s = s_vec[k]
-        phi = phi_vec[k]
-
-        #Convert to radians
-        phi = phi*pi/180
-
-        #Compute velocity vector (u, v)
-        u = s*cos(phi)
-        v = s*sin(phi)
-
-        #Compute wind stress
-        #S = const * sqrt(u**2 + v**2)
-        S = const * u
-        xmom_update[k] += S*u
-        #ymom_update[k] += S*v
Index: uga_work/development/pipeflow/test_shallow_water_domain.py
===================================================================
--- /anuga_work/development/pipeflow/test_shallow_water_domain.py	(revision 7787)
+++ 	(revision )
@@ -1,447 +1,0 @@
-#!/usr/bin/env python
-
-import unittest
-from math import sqrt, pi
-
-
-from shallow_water_domain import *
-from shallow_water_domain import flux_function as domain_flux_function
-
-
-
-import numpy
-
-
-class Test_Shallow_Water(unittest.TestCase):
-    def setUp(self):
-        self.points = [0.0, 1.0, 2.0, 3.0]
-        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
-        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
-        
-    def tearDown(self):
-        pass
-        #print "  Tearing down"
-
-
-    def test_creation(self):
-        domain = Domain(self.points)
-        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])
-
-    def test_compute_fluxes(self):
-        """
-        Compare shallow_water_domain flux calculation against a previous
-        Python implementation (defined in this file)
-        """
-        domain = Domain(self.points)
-        domain.set_quantity('stage',2.0)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        stage_ud, xmom_ud = local_compute_fluxes(domain)
-
-        domain.distribute_to_vertices_and_edges()
-        domain.compute_fluxes()
-
-        print domain.quantities['stage'].explicit_update
-        print domain.quantities['xmomentum'].explicit_update
-        print stage_ud
-        print xmom_ud
-
-        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
-        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
-
-
-    def test_local_flux_function(self):
-        normal = 1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        #This assumes h0 = 1.0e-3!!
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-        #print maxspeed
-        #print edgeflux
-        
-        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
-        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-        normal = -1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-
-
-        #print maxspeed
-        #print edgeflux        
-        
-        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
-        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-    def test_domain_flux_function(self):
-        normal = 1.0
-        ql = numpy.array([1.0, 2.0],numpy.float)
-        qr = numpy.array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-
-        #print edgeflux
-
-
-
-        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
-
-        #print domainedgeflux
-
-        assert numpy.allclose(domainedgeflux, edgeflux)
-        assert numpy.allclose(domainmaxspeed, maxspeed)
-
-    def test_gravity(self):
-        """
-        Compare shallow_water_domain gravity calculation
-        """
-
-        def slope_one(x):
-            return x
-
-        domain = Domain(self.points)
-        domain.set_quantity('stage',4.0)
-        domain.set_quantity('elevation',slope_one)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        gravity(domain)
-        
-        #print domain.quantities['stage'].vertex_values 
-        #print domain.quantities['elevation'].vertex_values 
-        #print domain.quantities['xmomentum'].explicit_update       
-
-        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
-
-
-    def xtest_evolve_first_order(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-        
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 1
-        domain.set_timestepping_method('euler')
-        
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            #pass
-        
-            print domain.quantities['stage'].vertex_values
-            print domain.quantities['elevation'].vertex_values
-            print domain.quantities['xmomentum'].vertex_values
-          
-        
-            print domain.quantities['stage'].centroid_values  
-            print domain.quantities['elevation'].centroid_values
-            print domain.quantities['xmomentum'].centroid_values    
-
-        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-    def xtest_evolve_euler_second_order_space(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 2
-        domain.set_timestepping_method('rk2')
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            pass
-        
-        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-    def xtest_evolve_second_order_space_time(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-
-        def slope_square(x):
-            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 2
-        domain.set_timestepping_method('rk3')
-
-        yieldstep=1.0
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            pass
-        
-        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-
-#==============================================================================
-
-def local_compute_fluxes(domain):
-    """Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    """
-
-    import sys
-    from Numeric import zeros, Float
-
-    N = domain.number_of_elements
-
-    tmp0 = zeros(N,Float)
-    tmp1 = zeros(N,Float)
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    #    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    #Arrays
-    #stage = Stage.edge_values
-    #xmom =  Xmom.edge_values
-    #   ymom =  Ymom.edge_values
-    #bed =   Bed.edge_values
-    
-    stage = Stage.vertex_values
-    xmom =  Xmom.vertex_values
-    bed =   Bed.vertex_values
-    
-    #print 'stage edge values', stage
-    #print 'xmom edge values', xmom
-    #print 'bed values', bed
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-    #print 'stage_bdry',stage_bdry
-    #print 'xmom_bdry', xmom_bdry
-    #    ymom_bdry =  Ymom.boundary_values
-
-    #    flux = zeros(3, Float) #Work numpy.array for summing up fluxes
-    flux = zeros(2, Float) #Work numpy.array for summing up fluxes
-    ql = zeros(2, Float)
-    qr = zeros(2, Float)
-
-    #Loop
-    timestep = float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work numpy.array
-        #for i in range(3):
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            #ql[0] = stage[k, i]
-            #ql[1] = xmom[k, i]
-            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                #print i, ' ' , m
-                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                qr[0] = stage[n, m]
-                qr[1] = xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-        
-            #Flux computation using provided function
-            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            #print 'ql',ql
-            #print 'qr',qr
-            
-
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-
-            #print 'edgeflux', edgeflux
-
-            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
-            # flux = edgefluxleft - edgefluxright
-            flux -= edgeflux #* domain.edgelengths[k,i]
-            #Update optimal_timestep
-            try:
-                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.CFL*0.5*domain.areas[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]
-        tmp0[k] = flux[0]
-        tmp1[k] = flux[1]
-
-
-    return tmp0, tmp1
-
-
-def local_flux_function(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numerical vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from config import g, epsilon, h0
-    from math import sqrt
-
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #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 +  h0/h_left)
-
-
-    uh_left = u_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 + h0/h_right)
-
-    uh_right = u_right*h_right
-    
-    #vh_left  = q_left[2]
-    #vh_right = q_right[2]
-
-    #print h_right
-    #print u_right
-    #print h_left
-    #print u_right
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    
-    #Flux computation
-
-    #flux_left  = numpy.array([u_left*h_left,
-    #                    u_left*uh_left + 0.5*g*h_left**2])
-    #flux_right = numpy.array([u_right*h_right,
-    #                    u_right*uh_right + 0.5*g*h_right**2])
-    flux_left  = numpy.array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left*h_left])
-    flux_right = numpy.array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right*h_right])
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = numpy.array([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[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
-    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
-
-
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
Index: uga_work/development/pipeflow/test_shallow_water_vel_domain.py
===================================================================
--- /anuga_work/development/pipeflow/test_shallow_water_vel_domain.py	(revision 7787)
+++ 	(revision )
@@ -1,517 +1,0 @@
-#!/usr/bin/env python
-
-import unittest
-from math import sqrt, pi
-
-
-from shallow_water_vel_domain import *
-
-
-from numpy import allclose, array, ones, numpy.float, maximum, zeros
-
-
-class Test_Shallow_Water(unittest.TestCase):
-    def setUp(self):
-        self.points = [0.0, 1.0, 2.0, 3.0]
-        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
-        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
-        
-    def tearDown(self):
-        pass
-        #print "  Tearing down"
-
-
-    def test_creation(self):
-        domain = Domain(self.points)
-        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
-
-    def test_reflective_boundary(self):
-        """
-        Test setting reflective boundary
-        """
-        
-        domain = Domain(self.points)
-        domain.set_quantity('stage',2.0)
-        domain.set_quantity('xmomentum',6.0)
-        
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-
-        domain.distribute_to_vertices_and_edges()
-        domain.update_boundary()
-
-##         print 'In test reflective'
-##         print domain.quantities['stage'].vertex_values
-##         print domain.quantities['xmomentum'].vertex_values
-##         print domain.quantities['elevation'].vertex_values
-##         print domain.quantities['height'].vertex_values
-##         print domain.quantities['velocity'].vertex_values
-        
-##         print domain.quantities['stage'].boundary_values
-##         print domain.quantities['xmomentum'].boundary_values
-##         print domain.quantities['elevation'].boundary_values
-##         print domain.quantities['height'].boundary_values
-##         print domain.quantities['velocity'].boundary_values
-
-        assert allclose( domain.quantities['stage'    ].boundary_values, [2.0, 2.0])
-        assert allclose( domain.quantities['xmomentum'].boundary_values, [-6.0, -6.0])
-        assert allclose( domain.quantities['elevation'].boundary_values, [0.0, 0.0])
-        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
-        assert allclose( domain.quantities['velocity' ].boundary_values, [-3.0, -3.0])
-
-        
-
-    def test_dirichlet_boundary(self):
-        """
-        Test setting dirichlet boundary
-        """
-        
-        domain = Domain(self.points)
-        domain.set_quantity('stage',2.0)
-        domain.set_quantity('xmomentum',6.0)
-        
-        domain.set_boundary({'exterior' : Dirichlet_boundary([3.0, 8.0, 1.0, 2.0, 4.0])})
-
-
-        domain.distribute_to_vertices_and_edges()
-        domain.update_boundary()
-
-##         print 'In test dirichlet'
-##         print domain.quantities['stage'].vertex_values
-##         print domain.quantities['xmomentum'].vertex_values
-##         print domain.quantities['elevation'].vertex_values
-##         print domain.quantities['height'].vertex_values
-##         print domain.quantities['velocity'].vertex_values
-        
-##         print domain.quantities['stage'].boundary_values
-##         print domain.quantities['xmomentum'].boundary_values
-##         print domain.quantities['elevation'].boundary_values
-##         print domain.quantities['height'].boundary_values
-##         print domain.quantities['velocity'].boundary_values
-
-        assert allclose( domain.quantities['stage'    ].boundary_values, [3.0, 3.0])
-        assert allclose( domain.quantities['xmomentum'].boundary_values, [8.0, 8.0])
-        assert allclose( domain.quantities['elevation'].boundary_values, [1.0, 1.0])
-        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
-        assert allclose( domain.quantities['velocity' ].boundary_values, [4.0, 4.0])
-
-
-    def test_compute_fluxes(self):
-        """
-        Compare shallow_water_domain flux calculation against a previous
-        Python implementation (defined in this file)
-        """
-        domain = Domain(self.points)
-        domain.set_quantity('stage',2.0)
-        domain.set_boundary({'exterior' : Dirichlet_boundary([0.0, 0.0, 0.0, 0.0, 0.0])})
-
-
-        domain.distribute_to_vertices_and_edges()
-        domain.update_boundary()
-
-        stage_ud, xmom_ud = local_compute_fluxes(domain)
-
-        domain.compute_fluxes()
-        
-        print domain.quantities['stage'].vertex_values
-        print domain.quantities['xmomentum'].vertex_values
-        print domain.quantities['elevation'].vertex_values
-        print domain.quantities['height'].vertex_values
-        print domain.quantities['velocity'].vertex_values
-
-        print domain.quantities['stage'].boundary_values
-        print domain.quantities['xmomentum'].boundary_values
-        print domain.quantities['elevation'].boundary_values
-        print domain.quantities['height'].boundary_values
-        print domain.quantities['velocity'].boundary_values
-        
-        print domain.quantities['stage'].explicit_update
-        print domain.quantities['xmomentum'].explicit_update
-        print stage_ud
-        print xmom_ud
-
-        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
-        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
-
-
-    def test_local_flux_function(self):
-        normal = 1.0
-        ql = array([1.0, 2.0],numpy.float)
-        qr = array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        #This assumes h0 = 1.0e-3!!
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-        #print maxspeed
-        #print edgeflux
-        
-        assert allclose(array([2.0, 8.9],numpy.float), edgeflux, rtol=1.0e-005)
-        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-        normal = -1.0
-        ql = array([1.0, 2.0],numpy.float)
-        qr = array([1.0, 2.0],numpy.float)
-        zl = 0.0
-        zr = 0.0
-
-        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
-
-
-        #print maxspeed
-        #print edgeflux        
-        
-        assert allclose(array([-2.0, -8.9],numpy.float), edgeflux, rtol=1.0e-005)
-        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
-
-
-    def test_gravity(self):
-        """
-        Compare shallow_water_domain gravity calculation
-        """
-
-        def slope_one(x):
-            return x
-
-        domain = Domain(self.points)
-        domain.set_quantity('stage',4.0)
-        domain.set_quantity('elevation',slope_one)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.distribute_to_vertices_and_edges()
-        domain.update_boundary()
-
-        gravity(domain)
-        
-        #print domain.quantities['stage'].vertex_values 
-        #print domain.quantities['elevation'].vertex_values 
-        #print domain.quantities['xmomentum'].explicit_update       
-
-        assert allclose( array([-34.3, -24.5, -14.7], numpy.float), domain.quantities['xmomentum'].explicit_update )
-
-
-    def test_evolve_first_order(self):
-        """
-        Compare still lake solution for various versions of shallow_water_domain
-        """
-        
-        def slope_square(x):
-            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-        domain = Domain(self.points2)
-        domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
-        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-        domain.default_order = 1
-        domain.set_timestepping_method('euler')
-        
-        yieldstep=0.25
-        finaltime=1.0
-
-        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-            pass
-        
-        ##	print domain.quantities['stage'].vertex_values 
-        ##	print domain.quantities['elevation'].vertex_values 
-        ##        print domain.quantities['xmomentum'].vertex_values
-        ##  
-        ##
-        ##        print domain.quantities['stage'].centroid_values  
-        ##	print domain.quantities['elevation'].centroid_values    
-        ##        print domain.quantities['xmomentum'].centroid_values    
-
-        #assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-        #assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-##     def test_evolve_euler_second_order_space(self):
-##         """
-##         Compare still lake solution for various versions of shallow_water_domain
-##         """
-
-##         def slope_square(x):
-##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-##         domain = Domain(self.points2)
-##         domain.set_quantity('stage',10.0)
-##         domain.set_quantity('elevation',slope_square)
-##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-##         domain.default_order = 2
-##         domain.set_timestepping_method('rk2')
-##         yieldstep=1.0
-##         finaltime=1.0
-
-##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-##             pass
-        
-##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-##     def test_evolve_second_order_space_time(self):
-##         """
-##         Compare still lake solution for various versions of shallow_water_domain
-##         """
-
-##         def slope_square(x):
-##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-
-##         domain = Domain(self.points2)
-##         domain.set_quantity('stage',10.0)
-##         domain.set_quantity('elevation',slope_square)
-##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
-
-##         domain.default_order = 2
-##         domain.set_timestepping_method('rk3')
-
-##         yieldstep=1.0
-##         finaltime=1.0
-
-##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
-##             pass
-        
-##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
-##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
-
-
-
-#==============================================================================
-
-def local_compute_fluxes(domain):
-    """Compute all fluxes and the timestep suitable for all volumes
-    in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-    """
-
-    import sys
-    from numpy import zeros, numpy.float
-
-    N = domain.number_of_elements
-
-    tmp0 = zeros(N,numpy.float)
-    tmp1 = zeros(N,numpy.float)
-
-    #Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    #    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-    #Arrays
-    #stage = Stage.edge_values
-    #xmom =  Xmom.edge_values
-    #   ymom =  Ymom.edge_values
-    #bed =   Bed.edge_values
-    
-    stage = Stage.vertex_values
-    xmom =  Xmom.vertex_values
-    bed =   Bed.vertex_values
-    
-    #print 'stage edge values', stage
-    #print 'xmom edge values', xmom
-    #print 'bed values', bed
-
-    stage_bdry = Stage.boundary_values
-    xmom_bdry =  Xmom.boundary_values
-    #print 'stage_bdry',stage_bdry
-    #print 'xmom_bdry', xmom_bdry
-    #    ymom_bdry =  Ymom.boundary_values
-
-    #    flux = zeros(3, numpy.float) #Work array for summing up fluxes
-    flux = zeros(2, numpy.float) #Work array for summing up fluxes
-    ql = zeros(2, numpy.float)
-    qr = zeros(2, numpy.float)
-
-    #Loop
-    timestep = numpy.float(sys.maxint)
-    enter = True
-    for k in range(N):
-
-        flux[:] = 0.  #Reset work array
-        #for i in range(3):
-        for i in range(2):
-            #Quantities inside volume facing neighbour i
-            #ql[0] = stage[k, i]
-            #ql[1] = xmom[k, i]
-            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
-                qr[1] = xmom_bdry[m]
-                zr = zl #Extend bed elevation to boundary
-            else:
-                #m = domain.neighbour_edges[k,i]
-                m = domain.neighbour_vertices[k,i]
-                #print i, ' ' , m
-                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                qr[0] = stage[n, m]
-                qr[1] = xmom[n, m]
-                zr = bed[n, m]
-
-
-            #Outward pointing normal vector
-            normal = domain.normals[k, i]
-        
-            #Flux computation using provided function
-            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            #print 'ql',ql
-            #print 'qr',qr
-            
-
-            edgeflux, max_speed = local_flux_function(normal, ql, qr, zl, zr)
-
-            #print 'edgeflux', edgeflux
-
-            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
-            # flux = edgefluxleft - edgefluxright
-            flux -= edgeflux #* domain.edgelengths[k,i]
-            #Update optimal_timestep
-            try:
-                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                timestep = min(timestep, domain.CFL*0.5*domain.areas[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]
-        tmp0[k] = flux[0]
-        tmp1[k] = flux[1]
-
-
-    return tmp0, tmp1
-
-
-def local_flux_function(normal, ql, qr, zl, zr):
-    """Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of w = h+z using the 'central scheme' as described in
-
-    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
-    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
-    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-
-    The implemented formula is given in equation (3.15) on page 714
-
-    Conserved quantities w, uh, are stored as elements 0 and 1
-    in the numpyal vectors ql an qr.
-
-    Bed elevations zl and zr.
-    """
-
-    from config import g, epsilon, h0
-    from math import sqrt
-    from numpy import array
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #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 +  h0/h_left)
-
-
-    uh_left = u_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 + h0/h_right)
-
-    uh_right = u_right*h_right
-    
-    #vh_left  = q_left[2]
-    #vh_right = q_right[2]
-
-    #print h_right
-    #print u_right
-    #print h_left
-    #print u_right
-
-    soundspeed_left  = sqrt(g*h_left)
-    soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    
-    #Flux computation
-
-    #flux_left  = array([u_left*h_left,
-    #                    u_left*uh_left + 0.5*g*h_left**2])
-    #flux_right = array([u_right*h_right,
-    #                    u_right*uh_right + 0.5*g*h_right**2])
-    flux_left  = array([u_left*h_left,
-                        u_left*uh_left + 0.5*g*h_left*h_left])
-    flux_right = array([u_right*h_right,
-                        u_right*uh_right + 0.5*g*h_right*h_right])
-
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = array([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[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed
-
-
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Shallow_Water, 'test')    
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
-
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
-    #print "restricted test"
-    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
Index: /anuga_work/development/pipeflow/test_sww_domain.py
===================================================================
--- /anuga_work/development/pipeflow/test_sww_domain.py	(revision 7788)
+++ /anuga_work/development/pipeflow/test_sww_domain.py	(revision 7788)
@@ -0,0 +1,428 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from sww_domain import *
+import sww_python
+
+
+import numpy
+
+
+class Test_Shallow_Water(unittest.TestCase):
+    def setUp(self):
+        self.points = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        domain = Domain(self.points)
+        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])
+
+    def test_compute_fluxes(self):
+        """
+        Compare shallow_water_domain flux calculation against a previous
+        Python implementation (defined in this file)
+        """
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        stage_ud, xmom_ud = sww_python.compute_fluxes(domain)
+
+        domain.distribute_to_vertices_and_edges()
+        domain.compute_fluxes()
+
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
+
+        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
+        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = numpy.array([1.0, 2.0],numpy.float)
+        qr = numpy.array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        #This assumes h0 = 1.0e-3!!
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        #print maxspeed
+        #print edgeflux
+        
+        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
+        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+        normal = -1.0
+        ql = numpy.array([1.0, 2.0],numpy.float)
+        qr = numpy.array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+
+        #print maxspeed
+        #print edgeflux        
+        
+        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
+        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+
+
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+
+        def slope_one(x):
+            return x
+
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
+
+
+    def xtest_evolve_first_order(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+        
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 1
+        domain.set_timestepping_method('euler')
+        
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            #pass
+        
+            print domain.quantities['stage'].vertex_values
+            print domain.quantities['elevation'].vertex_values
+            print domain.quantities['xmomentum'].vertex_values
+          
+        
+            print domain.quantities['stage'].centroid_values  
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values    
+
+        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+    def xtest_evolve_euler_second_order_space(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.set_timestepping_method('rk2')
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+    def xtest_evolve_second_order_space_time(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+
+        def slope_square(x):
+            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.set_timestepping_method('rk3')
+
+        yieldstep=1.0
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+
+#==============================================================================
+
+def local_compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    tmp0 = zeros(N,Float)
+    tmp1 = zeros(N,Float)
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    #    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
+    #   ymom =  Ymom.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    #print 'stage edge values', stage
+    #print 'xmom edge values', xmom
+    #print 'bed values', bed
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    #print 'stage_bdry',stage_bdry
+    #print 'xmom_bdry', xmom_bdry
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, Float) #Work numpy.array for summing up fluxes
+    flux = zeros(2, Float) #Work numpy.array for summing up fluxes
+    ql = zeros(2, Float)
+    qr = zeros(2, Float)
+
+    #Loop
+    timestep = float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work numpy.array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            #print 'ql',ql
+            #print 'qr',qr
+            
+
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.CFL*0.5*domain.areas[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]
+        tmp0[k] = flux[0]
+        tmp1[k] = flux[1]
+
+
+    return tmp0, tmp1
+
+
+def local_flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #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 +  h0/h_left)
+
+
+    uh_left = u_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 + h0/h_right)
+
+    uh_right = u_right*h_right
+    
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = numpy.array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = numpy.array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = numpy.array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = numpy.array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = numpy.array([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[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
+
+
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /anuga_work/development/pipeflow/test_sww_domain_vel.py
===================================================================
--- /anuga_work/development/pipeflow/test_sww_domain_vel.py	(revision 7788)
+++ /anuga_work/development/pipeflow/test_sww_domain_vel.py	(revision 7788)
@@ -0,0 +1,517 @@
+#!/usr/bin/env python
+
+import unittest
+from math import sqrt, pi
+
+
+from shallow_water_vel_domain import *
+
+
+from numpy import allclose, array, ones, numpy.float, maximum, zeros
+
+
+class Test_Shallow_Water(unittest.TestCase):
+    def setUp(self):
+        self.points = [0.0, 1.0, 2.0, 3.0]
+        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
+        
+    def tearDown(self):
+        pass
+        #print "  Tearing down"
+
+
+    def test_creation(self):
+        domain = Domain(self.points)
+        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
+
+    def test_reflective_boundary(self):
+        """
+        Test setting reflective boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test reflective'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [-6.0, -6.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [0.0, 0.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [-3.0, -3.0])
+
+        
+
+    def test_dirichlet_boundary(self):
+        """
+        Test setting dirichlet boundary
+        """
+        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_quantity('xmomentum',6.0)
+        
+        domain.set_boundary({'exterior' : Dirichlet_boundary([3.0, 8.0, 1.0, 2.0, 4.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+##         print 'In test dirichlet'
+##         print domain.quantities['stage'].vertex_values
+##         print domain.quantities['xmomentum'].vertex_values
+##         print domain.quantities['elevation'].vertex_values
+##         print domain.quantities['height'].vertex_values
+##         print domain.quantities['velocity'].vertex_values
+        
+##         print domain.quantities['stage'].boundary_values
+##         print domain.quantities['xmomentum'].boundary_values
+##         print domain.quantities['elevation'].boundary_values
+##         print domain.quantities['height'].boundary_values
+##         print domain.quantities['velocity'].boundary_values
+
+        assert allclose( domain.quantities['stage'    ].boundary_values, [3.0, 3.0])
+        assert allclose( domain.quantities['xmomentum'].boundary_values, [8.0, 8.0])
+        assert allclose( domain.quantities['elevation'].boundary_values, [1.0, 1.0])
+        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
+        assert allclose( domain.quantities['velocity' ].boundary_values, [4.0, 4.0])
+
+
+    def test_compute_fluxes(self):
+        """
+        Compare shallow_water_domain flux calculation against a previous
+        Python implementation (defined in this file)
+        """
+        domain = Domain(self.points)
+        domain.set_quantity('stage',2.0)
+        domain.set_boundary({'exterior' : Dirichlet_boundary([0.0, 0.0, 0.0, 0.0, 0.0])})
+
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        stage_ud, xmom_ud = local_compute_fluxes(domain)
+
+        domain.compute_fluxes()
+        
+        print domain.quantities['stage'].vertex_values
+        print domain.quantities['xmomentum'].vertex_values
+        print domain.quantities['elevation'].vertex_values
+        print domain.quantities['height'].vertex_values
+        print domain.quantities['velocity'].vertex_values
+
+        print domain.quantities['stage'].boundary_values
+        print domain.quantities['xmomentum'].boundary_values
+        print domain.quantities['elevation'].boundary_values
+        print domain.quantities['height'].boundary_values
+        print domain.quantities['velocity'].boundary_values
+        
+        print domain.quantities['stage'].explicit_update
+        print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
+        print xmom_ud
+
+        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
+        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],numpy.float)
+        qr = array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        #This assumes h0 = 1.0e-3!!
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+        #print maxspeed
+        #print edgeflux
+        
+        assert allclose(array([2.0, 8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],numpy.float)
+        qr = array([1.0, 2.0],numpy.float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
+
+
+        #print maxspeed
+        #print edgeflux        
+        
+        assert allclose(array([-2.0, -8.9],numpy.float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
+
+
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+
+        def slope_one(x):
+            return x
+
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.distribute_to_vertices_and_edges()
+        domain.update_boundary()
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert allclose( array([-34.3, -24.5, -14.7], numpy.float), domain.quantities['xmomentum'].explicit_update )
+
+
+    def test_evolve_first_order(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+        
+        def slope_square(x):
+            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 1
+        domain.set_timestepping_method('euler')
+        
+        yieldstep=0.25
+        finaltime=1.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        ##	print domain.quantities['stage'].vertex_values 
+        ##	print domain.quantities['elevation'].vertex_values 
+        ##        print domain.quantities['xmomentum'].vertex_values
+        ##  
+        ##
+        ##        print domain.quantities['stage'].centroid_values  
+        ##	print domain.quantities['elevation'].centroid_values    
+        ##        print domain.quantities['xmomentum'].centroid_values    
+
+        #assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        #assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+##     def test_evolve_euler_second_order_space(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk2')
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+##     def test_evolve_second_order_space_time(self):
+##         """
+##         Compare still lake solution for various versions of shallow_water_domain
+##         """
+
+##         def slope_square(x):
+##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+
+##         domain = Domain(self.points2)
+##         domain.set_quantity('stage',10.0)
+##         domain.set_quantity('elevation',slope_square)
+##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+##         domain.default_order = 2
+##         domain.set_timestepping_method('rk3')
+
+##         yieldstep=1.0
+##         finaltime=1.0
+
+##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+##             pass
+        
+##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+
+#==============================================================================
+
+def local_compute_fluxes(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+    """
+
+    import sys
+    from numpy import zeros, numpy.float
+
+    N = domain.number_of_elements
+
+    tmp0 = zeros(N,numpy.float)
+    tmp1 = zeros(N,numpy.float)
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    #    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']
+
+    #Arrays
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
+    #   ymom =  Ymom.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    #print 'stage edge values', stage
+    #print 'xmom edge values', xmom
+    #print 'bed values', bed
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    #print 'stage_bdry',stage_bdry
+    #print 'xmom_bdry', xmom_bdry
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, numpy.float) #Work array for summing up fluxes
+    flux = zeros(2, numpy.float) #Work array for summing up fluxes
+    ql = zeros(2, numpy.float)
+    qr = zeros(2, numpy.float)
+
+    #Loop
+    timestep = numpy.float(sys.maxint)
+    enter = True
+    for k in range(N):
+
+        flux[:] = 0.  #Reset work array
+        #for i in range(3):
+        for i in range(2):
+            #Quantities inside volume facing neighbour i
+            #ql[0] = stage[k, i]
+            #ql[1] = xmom[k, i]
+            ql = [stage[k, i], xmom[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[0] = stage_bdry[m]
+                qr[1] = xmom_bdry[m]
+                zr = zl #Extend bed elevation to boundary
+            else:
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
+                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                qr[0] = stage[n, m]
+                qr[1] = xmom[n, m]
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
+            normal = domain.normals[k, i]
+        
+            #Flux computation using provided function
+            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            #print 'ql',ql
+            #print 'qr',qr
+            
+
+            edgeflux, max_speed = local_flux_function(normal, ql, qr, zl, zr)
+
+            #print 'edgeflux', edgeflux
+
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
+            #Update optimal_timestep
+            try:
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, domain.CFL*0.5*domain.areas[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]
+        tmp0[k] = flux[0]
+        tmp1[k] = flux[1]
+
+
+    return tmp0, tmp1
+
+
+def local_flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numpyal vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon, h0
+    from math import sqrt
+    from numpy import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #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 +  h0/h_left)
+
+
+    uh_left = u_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 + h0/h_right)
+
+    uh_right = u_right*h_right
+    
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([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[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')    
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
+
+    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
+    #print "restricted test"
+    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
