Index: /anuga_work/development/pipeflow/swb_comp_flux_ext.c
===================================================================
--- /anuga_work/development/pipeflow/swb_comp_flux_ext.c	(revision 7782)
+++ /anuga_work/development/pipeflow/swb_comp_flux_ext.c	(revision 7782)
@@ -0,0 +1,329 @@
+#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(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(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(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(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, "swb_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(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", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initswb_comp_flux_ext(void){
+  Py_InitModule("swb_comp_flux_ext", MethodTable);
+  import_array();
+}
