Index: /anuga_work/development/pipeflow/sww_comp_flux_ext.c
===================================================================
--- /anuga_work/development/pipeflow/sww_comp_flux_ext.c	(revision 7783)
+++ /anuga_work/development/pipeflow/sww_comp_flux_ext.c	(revision 7783)
@@ -0,0 +1,312 @@
+#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) {
+
+    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"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initsww_comp_flux_ext(void){
+  Py_InitModule("sww_comp_flux_ext", MethodTable);
+  import_array();
+}
