Index: trunk/anuga_core/source/anuga/utilities/cg_ext.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/cg_ext.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/cg_ext.c	(revision 8691)
@@ -0,0 +1,581 @@
+// C extension for the cg_solve module. Implements a c routine of the 
+// conjugate gradient algorithm to solve Ax=b, using sparse matrix A in 
+// the csr format.
+//
+// See the module cg_solve.py
+//
+// Padarn Wilson 2012
+//
+// Note Padarn 26/11/12: Currently the input matrix for the CG solve must
+// be a Sparse_CSR matrix python object - defined in anuga/utilities/sparse.py
+//
+// Note Padarn 5/12/12: I have tried a few optimization modifications which
+// didn't seem to save any time:
+// -- conversion of the long arrays to int arrays (to save memory passing time)
+// -- taking advantage of the symmetric quality of the matrix A to reduce the zAx loop
+// -- specifying different 'chunk' sizes for the openmp loops
+// -- using blas instead of own openmp loops
+	
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include "stdio.h"
+#include "omp.h"
+
+
+// Dot product of two double vectors: a.b
+// @input N: int length of vectors a and b
+//        a: first vector of doubles
+//        b: second vector of double
+// @return: double result of a.b 
+double ddot( int N, double *a, double *b)
+{
+  double ret = 0;
+  int i;
+  #pragma omp parallel for private(i) reduction(+:ret)
+  for(i=0;i<N;i++)
+  {
+    ret+=a[i]*b[i];
+  }
+  return ret;
+
+}
+
+// In place multiplication of a double vector x by constant a: a*x
+// @input N: int length of vector x
+//        a: double scalar to multiply by
+//        x: double vector to scale
+void dscal(int N, double a, double *x)
+{
+  int i;
+  #pragma omp parallel for private(i)
+  for(i=0;i<N;i++)
+  {
+    x[i]=a*x[i];
+  }
+
+}
+
+// Copy of one vector to another - memory already allocated: y=x
+// @input N: int length of vectors x and y
+//        x: double vector to make copy of
+//        y: double vector to copy into
+void dcopy( int N, double *x, double *y)
+{
+  int i;
+  #pragma omp parallel for private(i)
+  for(i=0;i<N;i++)
+  {
+    y[i]=x[i];
+  }
+}
+
+// In place axpy operation: y = a*x + y
+// @input N: int length of vectors x and y
+//        a: double to multiply x by
+//        x: first double vector
+//        y: second double vector, stores result
+void daxpy(int N, double a, double *x, double *y)
+{
+  int i;
+  #pragma omp parallel for private(i)
+  for(i=0;i<N;i++)
+  {
+    y[i]=y[i]+a*x[i];
+  }
+}
+
+// Sparse CSR matrix-vector product: z = A*x
+// @input z: double vector to store the result
+//        data: double vector with non-zero entries of A
+//        colind: long vector of column indicies of non-zero entries of A
+//        row_ptr: long vector giving index of rows for non-zero entires of A
+//        x: double vector to be multiplied
+//        M: length of vector x
+void zAx(double * z, double * data, long * colind, long * row_ptr, double * x, int M){
+
+  
+  
+  long i, j, ckey;
+
+
+     
+    #pragma omp parallel for private(ckey,j,i)
+    for (i=0; i<M; i++){
+      z[i]=0;
+      for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
+        j = colind[ckey];
+        z[i] += data[ckey]*x[j];
+      }              
+    }
+  
+
+}
+
+// Diagonal matrix-vector product: z = D*x
+// @input z: double vector to store the result
+//        D: double vector of diagonal matrix
+//        x: double vector to be multiplied
+//        M: length of vector x
+void zDx(double * z, double * D, double * x, int M){
+
+  
+  long i, j, ckey;
+   
+    #pragma omp parallel for private(ckey,j,i)
+    for (i=0; i<M; i++){
+      z[i]=D[i]*x[i];              
+    }
+  
+
+}
+
+// Diagonal matrix-vector product: z = D*x
+// @input z: double vector to store the result
+//        D: double vector of diagonal matrix
+//        x: double vector to be multiplied
+//        M: length of vector x
+void zDinx(double * z, double * D, double * x, int M){
+
+  
+  long i, j, ckey;
+   
+    #pragma omp parallel for private(ckey,j,i)
+    for (i=0; i<M; i++){
+      z[i]=1.0/D[i]*x[i];              
+    }
+  
+
+}
+
+
+
+// Sparse CSR matrix-vector product and vector addition: z = a*A*x + y
+// @input z: double vector to store the result
+//        a: double to scale matrix-vector product by
+//        data: double vector with non-zero entries of A
+//        colind: long vector of column indicies of non-zero entries of A
+//        row_ptr: long vector giving index of rows for non-zero entires of A
+//        x: double vector to be multiplied
+//        y: double vector to add
+//        M: length of vector x
+void zaAxpy(double * z, double a, double * data, long * colind, long * row_ptr, double * x,
+      double * y,int M){
+  long i, j, ckey;
+  #pragma omp parallel for private(ckey,j,i)
+    for (i=0; i<M; i++ ){
+      z[i]=y[i];
+
+      for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
+        j = colind[ckey];
+        z[i] += a*data[ckey]*x[j];
+      }              
+
+  }
+
+}
+
+// Jacobi preconditioner for matrix, A, and right hand side, b. Mutiplies each row
+// by one divided by the diagonal element of the matrix. If the diagonal 
+// element is zero, does nothing (should nnot occur)
+//        colind: long vector of column indicies of non-zero entries of A
+//        row_ptr: long vector giving index of rows for non-zero entires of A
+//        b: double vector specifying right hand side of equation to solve
+//        M: length of vector b
+
+int _jacobi_precon_c(double* data, 
+                long* colind,
+                long* row_ptr,
+                double * precon,
+                int M){
+
+
+  long i, j, k, ckey;
+  double diag;
+
+
+     
+  #pragma omp parallel for private(diag,ckey,j,i)
+  for (i=0; i<M; i++){
+    diag = 0;
+    for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
+      j = colind[ckey];
+      if (i==j){
+        diag = data[ckey];
+      }
+    }
+    if (diag == 0){
+      diag =1;
+    }
+    precon[i]=diag;
+  }
+  
+  return 0;
+
+}
+
+// Conjugate gradient solve Ax = b for x, A given in Sparse CSR format
+// @input data: double vector with non-zero entries of A
+//        colind: long vector of column indicies of non-zero entries of A
+//        row_ptr: long vector giving index of rows for non-zero entires of A
+//        b: double vector specifying right hand side of equation to solve
+//        x: double vector with initial guess and to store result
+//        imax: maximum number of iterations
+//        tol: error tollerance for stopping criteria
+//        M: length of vectors x and b
+// @return: 0 on success  
+int _cg_solve_c(double* data, 
+                long* colind,
+                long* row_ptr,
+                double * b,
+                double * x,
+                int imax,
+                double tol,
+                double a_tol,
+                int M){
+
+  int i = 1;
+  double alpha,rTr,rTrOld,bt,rTr0;
+
+  double * d = malloc(sizeof(double)*M);
+  double * r = malloc(sizeof(double)*M);
+  double * q = malloc(sizeof(double)*M);
+  double * xold = malloc(sizeof(double)*M);
+
+  zaAxpy(r,-1.0,data,colind,row_ptr,x,b,M);
+  dcopy(M,r,d);
+
+  rTr=ddot(M,r,r);
+  rTr0 = rTr;
+  
+  while((i<imax) && (rTr>pow(tol,2)*rTr0) && (rTr > pow(a_tol,2))){
+
+    zAx(q,data,colind,row_ptr,d,M);
+    alpha = rTr/ddot(M,d,q);
+    dcopy(M,x,xold);
+    daxpy(M,alpha,d,x);
+
+    daxpy(M,-alpha,q,r);
+    rTrOld = rTr;
+    rTr = ddot(M,r,r);
+
+    bt= rTr/rTrOld;
+
+    dscal(M,bt,d);
+    daxpy(M,1.0,r,d);
+
+    i=i+1;
+
+  }
+  
+  free(d);
+  free(r);
+  free(q);
+  free(xold);
+
+  if (i>=imax){
+    return -1;
+  }
+  else{
+    return 0;
+  }
+  
+
+}          
+
+// Conjugate gradient solve Ax = b for x, A given in Sparse CSR format,
+// using a diagonal preconditioner M. 
+// @input data: double vector with non-zero entries of A
+//        colind: long vector of column indicies of non-zero entries of A
+//        row_ptr: long vector giving index of rows for non-zero entires of A
+//        b: double vector specifying right hand side of equation to solve
+//        x: double vector with initial guess and to store result
+//        imax: maximum number of iterations
+//        tol: error tollerance for stopping criteria
+//        M: length of vectors x and b
+//        precon: diagonal preconditioner given as vector
+// @return: 0 on success  
+int _cg_solve_c_precon(double* data, 
+                long* colind,
+                long* row_ptr,
+                double * b,
+                double * x,
+                int imax,
+                double tol,
+                double a_tol,
+                int M,
+                double * precon){
+
+  int i = 1;
+  double alpha,rTr,rTrOld,bt,rTr0;
+
+  double * d = malloc(sizeof(double)*M);
+  double * r = malloc(sizeof(double)*M);
+  double * q = malloc(sizeof(double)*M);
+  double * xold = malloc(sizeof(double)*M);
+  double * rhat = malloc(sizeof(double)*M);
+  double * temp = malloc(sizeof(double)*M);
+
+  zaAxpy(r,-1.0,data,colind,row_ptr,x,b,M);
+  zDinx(rhat,precon,r,M);
+  dcopy(M,rhat,d);
+
+  rTr=ddot(M,r,rhat);
+  rTr0 = rTr;
+  
+  while((i<imax) && (rTr>pow(tol,2)*rTr0) && (rTr > pow(a_tol,2))){
+
+    zAx(q,data,colind,row_ptr,d,M);
+    alpha = rTr/ddot(M,d,q);
+    dcopy(M,x,xold);
+    daxpy(M,alpha,d,x);
+
+    daxpy(M,-alpha,q,r);
+    zDinx(rhat,precon,r,M);
+    rTrOld = rTr;
+    rTr = ddot(M,r,rhat);
+
+    bt= rTr/rTrOld;
+
+    dscal(M,bt,d);
+    daxpy(M,1.0,rhat,d);
+
+    i=i+1;
+
+  }
+  free(temp);
+  free(rhat);
+  free(d);
+  free(r);
+  free(q);
+  free(xold);
+
+  if (i>=imax){
+    return -1;
+  }
+  else{
+    return 0;
+  }
+  
+
+}       
+
+		     
+/////////////////////////////////////////////////
+// Gateways to Python
+
+PyObject *jacobi_precon_c(PyObject *self, PyObject *args){
+
+  int M,err,bcols;
+  
+  PyObject *csr_sparse; // input sparse matrix (must be CSR format)
+ 
+  PyArrayObject 
+    *data,            //Non Zeros Data array
+    *colind,          //Column indices array
+    *row_ptr,         //Row pointers array
+    *precon;                //Right hand side
+
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OO", &csr_sparse, &precon)) {
+    PyErr_SetString(PyExc_RuntimeError, "jacobi_precon_c could not parse input");  
+    return NULL;
+  }
+
+  // Extract three subarrays making up the sparse matrix in CSR format.
+  data = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "data");     
+  if (!data) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Data array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  colind = (PyArrayObject*)
+    PyObject_GetAttrString(csr_sparse, "colind"); 
+  if (!colind) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Column index array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  row_ptr = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "row_ptr");   
+  if (!row_ptr) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Row pointer array could not be allocated in cg_solve_c"); 
+  }
+  
+  M = (row_ptr -> dimensions[0])-1;
+
+  // Precondition the matrix and right hand side
+  err = _jacobi_precon_c((double*) data->data, 
+                (long*) colind->data,
+                (long*) row_ptr->data,
+                (double *) precon->data,
+                M);
+
+  // Free extra references to sparse matrix parts
+  Py_DECREF(data);    
+  Py_DECREF(colind);    
+  Py_DECREF(row_ptr);     
+
+  return Py_BuildValue("");
+
+
+}
+
+PyObject *cg_solve_c(PyObject *self, PyObject *args) {
+  
+  
+  PyObject *csr_sparse; // input sparse matrix (must be CSR format)
+ 
+  int imax,M,err,bcols;
+  double tol,a_tol;
+
+  PyArrayObject 
+    *data,            //Non Zeros Data array
+    *colind,          //Column indices array
+    *row_ptr,         //Row pointers array
+    *x0,               //Initial guess - and sotrage of result.
+    *b;                //Right hand side
+
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOOiddi", &csr_sparse, &x0, &b, &imax, &tol, &a_tol, &bcols)) {
+    PyErr_SetString(PyExc_RuntimeError, "cg_solve_c could not parse input");  
+    return NULL;
+  }
+
+  // Extract three subarrays making up the sparse matrix in CSR format.
+  data = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "data");     
+  if (!data) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Data array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  colind = (PyArrayObject*)
+    PyObject_GetAttrString(csr_sparse, "colind"); 
+  if (!colind) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Column index array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  row_ptr = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "row_ptr");   
+  if (!row_ptr) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Row pointer array could not be allocated in cg_solve_c"); 
+  }
+  
+  M = (row_ptr -> dimensions[0])-1;
+  
+  // Solve system uisng conjugate gradient
+  err = _cg_solve_c((double*) data->data, 
+                (long*) colind->data,
+                (long*) row_ptr->data,
+                (double *) b->data,
+                (double *) x0->data,
+                imax,
+                tol,
+                a_tol,
+                M);
+  
+  // Free extra references to sparse matrix parts
+  Py_DECREF(data);    
+  Py_DECREF(colind);    
+  Py_DECREF(row_ptr);     
+
+  return Py_BuildValue("i",err);
+}
+
+PyObject *cg_solve_c_precon(PyObject *self, PyObject *args) {
+  
+  
+  PyObject *csr_sparse; // input sparse matrix (must be CSR format)
+ 
+  int imax,M,err,bcols;
+  double tol,a_tol;
+
+  PyArrayObject 
+    *data,            //Non Zeros Data array
+    *colind,          //Column indices array
+    *row_ptr,         //Row pointers array
+    *x0,               //Initial guess - and sotrage of result.
+    *b,                //Right hand side
+    *precon;           //diagonal preconditioner
+
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOOiddiO", &csr_sparse, &x0, &b, &imax, &tol, &a_tol, &bcols, &precon)) {
+    PyErr_SetString(PyExc_RuntimeError, "cg_solve_c_precon could not parse input");  
+    return NULL;
+  }
+
+  // Extract three subarrays making up the sparse matrix in CSR format.
+  data = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "data");     
+  if (!data) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Data array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  colind = (PyArrayObject*)
+    PyObject_GetAttrString(csr_sparse, "colind"); 
+  if (!colind) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Column index array could not be allocated in cg_solve_c");      
+    return NULL;
+  }  
+
+  row_ptr = (PyArrayObject*) 
+    PyObject_GetAttrString(csr_sparse, "row_ptr");   
+  if (!row_ptr) {
+    PyErr_SetString(PyExc_RuntimeError, 
+        "Row pointer array could not be allocated in cg_solve_c"); 
+  }
+  
+  M = (row_ptr -> dimensions[0])-1;
+  
+  // Solve system uisng conjugate gradient
+  err = _cg_solve_c_precon((double*) data->data, 
+                (long*) colind->data,
+                (long*) row_ptr->data,
+                (double *) b->data,
+                (double *) x0->data,
+                imax,
+                tol,
+                a_tol,
+                M,
+                (double *) precon->data);
+  
+  // Free extra references to sparse matrix parts
+  Py_DECREF(data);    
+  Py_DECREF(colind);    
+  Py_DECREF(row_ptr);     
+
+  return Py_BuildValue("i",err);
+}
+
+
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  {"cg_solve_c", cg_solve_c, METH_VARARGS, "Print out"},
+  {"cg_solve_c_precon", cg_solve_c_precon, METH_VARARGS, "Print out"},
+  {"jacobi_precon_c", jacobi_precon_c, METH_VARARGS, "Print out"},    
+  {NULL, NULL, 0, NULL}   /* sentinel */
+};
+
+// Module initialisation   
+void initcg_ext(void){
+  Py_InitModule("cg_ext", MethodTable);
+  
+  import_array();     //Necessary for handling of NumPY structures  
+}
+
Index: trunk/anuga_core/source/anuga/utilities/quad_tree.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/quad_tree.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/quad_tree.c	(revision 8691)
@@ -0,0 +1,522 @@
+#include "quad_tree.h"
+
+// **************** UTILITIES ***********************
+
+static void *emalloc(size_t amt,char * location)
+{
+    void *v = malloc(amt);  
+    if(!v){
+        fprintf(stderr, "out of mem in quad_tree: %s\n",location);
+        exit(EXIT_FAILURE);
+    }
+    return v;
+};
+
+double dot_points(double p1x,double p1y,double p2x,double p2y)
+{
+    
+    return p1x * p2x + p1y * p2y;
+    
+};
+
+// ***************************************************
+
+// ******************* TRIANGLE **********************
+
+triangle * new_triangle(int index,double x1,double y1,double x2,double y2,double x3,double y3)
+{
+
+	triangle * T = emalloc(sizeof(triangle),"new_triangle"); 
+	 
+	T->index = index;
+	T->x1 = x1; T->x2 = x2; T->x3 = x3;
+    T->y1 = y1; T->y2 = y2; T->y3 = y3;
+    T->next = NULL;
+
+    // Calculate the normals
+    // normal for vector from x1->x2
+    double nx_temp,ny_temp,dot_temp;
+    ny_temp = ( T->x2 - T->x1 );
+    nx_temp = -( T->y2 - T->y1 );
+    dot_temp = dot_points(nx_temp, ny_temp, nx_temp, ny_temp);
+    T->nx3 = nx_temp / sqrt(dot_temp);
+    T->ny3 = ny_temp / sqrt(dot_temp);
+    if( dot_points(T->nx3, T->ny3, T->x3 - T->x2, T->y3 - T->y2) > 0 ){
+    	T->nx3 = -T->nx3;
+    	T->ny3 = -T->ny3;
+    }
+    // normal for vector from x2->x3
+    ny_temp = ( T->x3 - T->x2 );
+    nx_temp = -(T->y3 - T ->y2 );
+    dot_temp = dot_points(nx_temp, ny_temp, nx_temp, ny_temp);
+    T->nx1 = nx_temp / sqrt(dot_temp);
+    T->ny1 = ny_temp / sqrt(dot_temp);
+    if( dot_points(T->nx1, T->ny1, T->x1 - T->x3, T->y1 - T->y3) > 0 ){
+    	T->nx1 = -T->nx1;
+    	T->ny1 = -T->ny1;
+    }
+    // normal for vector from x3->x1
+    ny_temp = ( T->x1 - T->x3 );
+    nx_temp = -( T->y1 - T->y3 );
+    dot_temp = dot_points(nx_temp, ny_temp, nx_temp, ny_temp);
+    T->nx2 = nx_temp / sqrt(dot_temp);
+    T->ny2 = ny_temp / sqrt(dot_temp);
+    if( dot_points(T->nx2, T->ny2, T->x2 - T->x1, T->y2 - T->y1) > 0 ){
+    	T->nx2 = -T->nx2;
+    	T->ny2 = -T->ny2;
+    }
+
+    return T;
+};
+
+void delete_triangle_list(triangle * T)
+{
+	while (T != NULL){
+     triangle * next = T->next;
+	   free(T);
+	   T = NULL;
+     T = next;
+  }
+};
+
+double * calculate_sigma(triangle * T,double x,double y)
+{
+
+	double  * ret_sigma = malloc(3 * sizeof(double));
+	ret_sigma[0] = dot_points(x - T->x2, y - T->y2, T->nx1, T->ny1)/
+					dot_points(T->x1 - T->x2, T->y1 - T->y2, T->nx1, T->ny1);
+	ret_sigma[1] = dot_points(x - T->x3, y - T->y3, T->nx2, T->ny2)/
+					dot_points(T->x2 - T->x3, T->y2 - T->y3, T->nx2, T->ny2);
+	ret_sigma[2] = dot_points(x - T->x1, y - T->y1, T->nx3, T->ny3)/
+					dot_points(T->x3 - T->x1, T->y3 - T->y1, T->nx3, T->ny3);
+	return ret_sigma;				
+};
+
+double dist(double x,
+	    double y) {
+  
+  return sqrt(x*x + y*y);
+}
+
+int __point_on_line(double x, double y,
+		    double x0, double y0,
+		    double x1, double y1,
+		    double rtol,
+		    double atol) {
+  /*Determine whether a point is on a line segment
+
+    Input: x, y, x0, x0, x1, y1: where
+        point is given by x, y
+	line is given by (x0, y0) and (x1, y1)
+
+  */
+
+  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
+  double nominator, denominator;
+  int is_parallel;
+
+  a0 = x - x0;
+  a1 = y - y0;
+
+  a_normal0 = a1;
+  a_normal1 = -a0;
+
+  b0 = x1 - x0;
+  b1 = y1 - y0;
+
+  nominator = fabs(a_normal0*b0 + a_normal1*b1);
+  denominator = b0*b0 + b1*b1;
+  
+  // Determine if line is parallel to point vector up to a tolerance
+  is_parallel = 0;
+  if (denominator == 0.0) {
+    // Use absolute tolerance
+    if (nominator <= atol) {
+      is_parallel = 1;
+    }
+  } else {
+    // Denominator is positive - use relative tolerance
+    if (nominator/denominator <= rtol) {
+      is_parallel = 1;
+    }    
+  }
+    
+  if (is_parallel) {
+    // Point is somewhere on the infinite extension of the line
+    // subject to specified absolute tolerance
+
+    len_a = dist(a0, a1); //sqrt(a0*a0 + a1*a1);
+    len_b = dist(b0, b1); //sqrt(b0*b0 + b1*b1);
+
+    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
+      return 1;
+    } else {
+      return 0;
+    }
+  } else {
+    return 0;
+  }
+};
+
+int __is_inside_triangle(double* point,
+			 double* triangle,
+			 int closed,
+			 double rtol,
+			 double a_tol) {
+			 
+  double vx, vy, v0x, v0y, v1x, v1y;
+  double a00, a10, a01, a11, b0, b1;
+  double denom, alpha, beta;
+  
+  double x, y; // Point coordinates
+  int i, j, res;
+
+  x = point[0];
+  y = point[1];
+  
+  // Quickly reject points that are clearly outside
+  if ((x < triangle[0]) && 
+      (x < triangle[2]) && 
+      (x < triangle[4])) return 0;       
+      
+  if ((x > triangle[0]) && 
+      (x > triangle[2]) && 
+      (x > triangle[4])) return 0;             
+  
+  if ((y < triangle[1]) && 
+      (y < triangle[3]) && 
+      (y < triangle[5])) return 0;       
+      
+  if ((y > triangle[1]) && 
+      (y > triangle[3]) && 
+      (y > triangle[5])) return 0;             
+  
+  
+  // v0 = C-A 
+  v0x = triangle[4]-triangle[0]; 
+  v0y = triangle[5]-triangle[1];
+  
+  // v1 = B-A   
+  v1x = triangle[2]-triangle[0]; 
+  v1y = triangle[3]-triangle[1];
+
+  // First check if point lies wholly inside triangle
+  a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0)
+  a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1)
+  a10 = a01;               // innerproduct(v1, v0)
+  a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1)
+    
+  denom = a11*a00 - a01*a10;
+
+  if (fabs(denom) > 0.0) {
+    // v = point-A  
+    vx = x - triangle[0]; 
+    vy = y - triangle[1];     
+    
+    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)        
+    b1 = v1x*vx + v1y*vy; // innerproduct(v1, v)            
+    
+    alpha = (b0*a11 - b1*a01)/denom;
+    beta = (b1*a00 - b0*a10)/denom;        
+    
+    if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1;
+  }
+
+  if (closed) {
+    // Check if point lies on one of the edges
+        
+    for (i=0; i<3; i++) {
+      j = (i+1) % 3; // Circular index into triangle array
+
+      res = __point_on_line(x, y,
+                            triangle[2*i], triangle[2*i+1], 
+                            triangle[2*j], triangle[2*j+1], 			    
+			    rtol, a_tol);
+      if (res) return 1;
+    }
+  };
+  // Default return if point is outside triangle			 
+  return 0;			 			 
+}			  			       	
+
+int triangle_contains_point(triangle * T,double pointx,double pointy)
+{
+    
+//    double v0x,v0y,v1x,v1y,v2x,v2y,dot00,dot01,dot02,dot11,dot12,invDenom,u,v;
+//    
+//    v0x = T->x3 - T->x1;
+//    v0y = T->y3 - T->y1;
+//    v1x = T->x2 - T->x1;
+//    v1y = T->y2 - T->y1;
+//    v2x = pointx - T->x1;
+//    v2y = pointy - T->y1;
+//    
+//    dot00 = dot_points(v0x, v0y, v0x, v0y);
+//    dot01 = dot_points(v0x, v0y, v1x, v1y);
+//    dot02 = dot_points(v0x, v0y, v2x, v2y);
+//    dot11 = dot_points(v1x, v1y, v1x, v1y);
+//    dot12 = dot_points(v1x, v1y, v2x, v2y);
+//    
+//    invDenom = 1/(dot00*dot11-dot01*dot01);
+//    u=(dot11 * dot02 - dot01 * dot12) * invDenom;
+//    v=(dot00 * dot12 - dot01 * dot02) * invDenom;
+//    
+//    if( u>=0 && v>=0 && v+u<1) return 1;
+//   
+//    else return 0;
+      double tri[6];
+      tri[0]=T->x1; tri[1]=T->y1; tri[2]=T->x2; tri[3]=T->y2; tri[4]=T->x3; tri[5]=T->y3;
+      double point[2];
+      point[0] = pointx; point[1] = pointy;
+      double rtol=1.0e-12;
+      double a_tol=1.0e-12;
+      int closed = 1;
+
+      return __is_inside_triangle(point,
+			 tri,
+			 closed,
+			 rtol,
+			 a_tol); 
+
+
+};
+
+
+
+// ***************************************************
+
+// ***************** quad_tree *********************
+
+quad_tree * new_quad_tree(double xmin, double xmax, double ymin, double ymax)
+{
+
+	quad_tree * ret = emalloc(sizeof(quad_tree),"new_quad_tree");
+	ret -> xmin = xmin; ret-> xmax = xmax; 
+	ret -> ymin = ymin; ret -> ymax = ymax;
+	ret -> parent = NULL;
+	ret -> q[0] = NULL; ret -> q[1] = NULL; ret -> q[2] = NULL; ret -> q[3] = NULL;
+	ret -> leaves = NULL;
+	ret -> end_leaves = NULL;
+  ret -> count = 0;
+	return ret; 
+
+};
+
+void delete_quad_tree(quad_tree * quadtree)
+{
+
+  quad_tree_ll * nodelist = new_quad_tree_ll(quadtree,0);
+  quad_tree_ll * last = nodelist;
+  quad_tree_ll * temp;
+  int i;
+
+  while(nodelist !=NULL){
+      
+      quadtree=nodelist->tree;
+      // if children have been added, add to the linked list
+
+      if (quadtree->q[0]!=NULL){
+          for (i=0;i<4;i++){
+              quad_tree_ll * child = new_quad_tree_ll(quadtree->q[i],0);
+              last->next=child;
+              last=child;
+          }
+      }
+      
+      if (quadtree->leaves!=NULL){
+          delete_triangle_list(quadtree->leaves);
+      }
+
+      free(quadtree);
+      quadtree=NULL;
+
+      temp = nodelist;
+      nodelist=nodelist->next;
+      free(temp);
+  }
+
+};
+
+void quad_tree_make_children(quad_tree *node){
+
+  double xmid = (node->xmin+node->xmax)/2;
+  double ymid = (node->ymin+node->ymax)/2;
+  double width = (node->xmax-node->xmin);
+  double height = (node->ymax-node->ymin);
+  // add quads 1-4
+  // include border expansion
+  double border=0.55;
+  node->q[0] = new_quad_tree(node->xmax-width*border,node->xmax,node->ymax-height*border,node->ymax);
+  node->q[0]->parent = node; 
+  
+  node->q[1] = new_quad_tree(node->xmin,node->xmin+width*border,node->ymax-height*border,node->ymax);
+  node->q[1]->parent = node; 
+  
+  node->q[2] = new_quad_tree(node->xmin,node->xmin+width*border,node->ymin,node->ymin+height*border);
+  node->q[2]->parent = node; 
+  
+  node->q[3] = new_quad_tree(node->xmax-width*border,node->xmax,node->ymin,node->ymin+height*border);
+  node->q[3]->parent = node; 
+
+}
+
+void quad_tree_add_triangle_to_list(quad_tree *node,triangle *T){
+
+  if (node->leaves == NULL){
+    // no current leaves
+    node->leaves = T;
+    node->end_leaves = T;
+  } else {
+    node->end_leaves->next = T;
+    node->end_leaves = T;
+  }
+
+}
+
+void quad_tree_insert_triangle(quad_tree *node,triangle *T)
+{
+	
+	// find the quadrant of the current node's extents in which the
+	// point lies (zero if intersects center of extents axes).
+
+	int quad = trivial_contain_split(node,T);
+
+  // always increase point count, as storing the total in tree below
+  node->count+=1;
+
+	if (quad != 0){
+		// if current node has no children yet, split:
+		if(node->q[0] == NULL){
+			
+			quad_tree_make_children(node); 
+			
+	    }
+	    // insert triangle into node corresponding to given quadrant
+	   	quad_tree_insert_triangle(node->q[quad-1],T);
+		return;
+		
+	}
+	// if triangle intersects the center axes of the node's extents, insert
+	// the triangle here
+	quad_tree_add_triangle_to_list(node,T);
+
+};
+
+
+int trivial_contain_split(quad_tree *node, triangle *T){
+
+	int p1 = trivial_contain_split_point(node,T->x1,T->y1);
+	int p2 = trivial_contain_split_point(node,T->x2,T->y2);
+	int p3 = trivial_contain_split_point(node,T->x3,T->y3);
+	if(p1 == p2 && p2 == p3){
+	 	return p1;
+	}
+	return 0;
+};
+
+int trivial_contain_split_point(quad_tree *node, double xp,double yp)
+{
+
+	double midx = (node->xmin+node->xmax)/2;
+	double midy = (node->ymin+node->ymax)/2;
+
+	int ret=0;
+	
+	if (midx < xp){
+		// quad 1 or 4
+		if (midy < yp){
+			ret = 1;
+		} else if (midy > yp){
+			ret = 4;
+		}
+	} else if (midx > xp){
+		//quad 2 or 3
+		if (midy < yp){
+			ret = 2;
+		} else if (midy > yp){
+			ret = 3;
+		}
+	}
+	return ret;
+};
+
+triangle * search_triangles_of_quad_tree(quad_tree * node,double xp,double yp){
+	
+    triangle * T = node->leaves;
+
+    while(T != NULL){
+        if(triangle_contains_point(T,xp,yp)){
+            return T; // Triangle contains point so return
+        }
+        T = T->next;
+    }
+	return T; // should be NULL if this is reached
+};
+
+// Searches the quad tree starting at 'cur_quad_tree' for the 
+// point, returning the triangle that contains it, or NULL
+// if no triangle is found.
+triangle * search(quad_tree * node, double xp, double yp){
+
+    triangle * return_T = NULL;
+
+    if(node->leaves!=NULL)
+    {
+    	return_T = search_triangles_of_quad_tree(node,xp,yp);
+    }
+    if(return_T != NULL) return return_T;
+	    
+    else
+    {
+        if(node->q[0]!=NULL) // look for child to search
+        {
+            //find correct quadrant to search
+            int quad = trivial_contain_split_point(node,xp,yp);
+            
+            if (quad!=0)
+            {
+            return_T = search(node->q[quad-1],xp,yp); // search child
+            }
+
+        return return_T; // return NULL pointer as no triangle
+        }
+	}
+	return return_T; // should not be reached
+};
+
+int quad_tree_node_count(quad_tree * tree)
+{
+  int node_count = 1;
+  if (tree->q[0]!=NULL){
+      int i;
+      for(i=0;i<4;i++){
+        node_count+=quad_tree_node_count(tree->q[i]);
+      }
+  }
+  return node_count;
+};
+
+// ***************************************************
+
+// ***************** quad_tree_ll *******************
+
+quad_tree_ll * new_quad_tree_ll(quad_tree * start,int index){
+    quad_tree_ll * list = malloc(sizeof(quad_tree_ll));
+    list->tree = start;
+    list->next = NULL;
+    list->index = index;
+    return list;
+}
+
+// ***************************************************
+
+// ***************** queue_ll *******************
+
+queue_ll * new_queue_ll(int node){
+    queue_ll * list = malloc(sizeof(queue_ll));
+    list->node=node;
+    list->next = NULL;
+    return list;
+}
+
+// ***************************************************
Index: trunk/anuga_core/source/anuga/utilities/quad_tree.h
===================================================================
--- trunk/anuga_core/source/anuga/utilities/quad_tree.h	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/quad_tree.h	(revision 8691)
@@ -0,0 +1,151 @@
+//------------------------------------------------------------------------------
+// quad_tree.h - Header for c implementation of the generic quad_tree. The 
+// 				   tree stores a set of triangles and then facilitates quickly
+//                 searching for a triangles in the list containing a given 
+//                 point.
+// 				   Tree has 'new' and 'delete' methods that should be called to
+//                 make sure memory is allocated and freed after use.
+// author: Padarn Wilson, date: 25/10/120 
+//-------------------------------------------------------------------------------
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+
+#ifndef quad_tree_H
+#define quad_tree_H
+
+
+// ***********************************************************************
+// triangle struct - Struct that stores data for triangles being inseted
+//                   into the tree. Triangles form a linked list for easy
+//                   storage.
+
+typedef struct triangle{
+
+	// points defining triangle.
+	double x1,y1;
+	double x2,y2;
+	double x3,y3;
+
+	// index stores the triangles unique id.
+	int index;
+
+	// outward normal vectors of triangles sides.
+	double nx1,ny1;
+	double nx2,ny2;
+	double nx3,ny3;
+
+	// next triangle turns the struct into a linked list.
+	struct triangle * next;  
+
+} triangle;
+
+// creates a new triangle and returns a pointer to the malloc'ed memory. 
+triangle * new_triangle(int index, double x1, double x2, double x3,
+						double y1, double y2, double y3);
+
+// deletes entire list of triangles
+void delete_triangle_list(triangle * T);
+
+// take a point and calculate 'sigma' with the given triangle. returns
+// a pointer to malloc'ed memory of a double array.
+double * calculate_sigma(triangle * T,double x,double y);
+
+// Tests to see if a triangle contains a given point,
+// returns a int value 0 false, 1 true.
+int triangle_contains_point(triangle * T,double pointx,double pointy);
+
+//**************************************************************************
+
+
+//**************************************************************************
+// quad_tree struct - Each quad_tree struct is a node in the tree, and 
+//                      can be treated as an independent quad_tree.
+
+typedef struct quad_tree{
+
+	// rectangular extents of current node
+	double xmin,xmax,ymin,ymax;
+	int count;
+	// parent and children of quad_tree
+	struct quad_tree *parent;
+	struct quad_tree *q[4];
+
+	// triangle stored in this node - leaves is a linked list of triangles
+	triangle * leaves;
+	// triangle end_leaves allows easy adding onto end of triangles
+	triangle * end_leaves;
+
+} quad_tree;
+
+// returns a new quad_tree pointer with malloc'ed memory. Inputs are the 
+// extents of the node.
+quad_tree* new_quad_tree(double xmin, double xmax, double ymin, double ymax);
+
+// delete a search tree - recursively deletes all children.
+void delete_quad_tree(quad_tree * tree);
+
+// add a new triangle to the quad_tree
+void quad_tree_insert_triangle(quad_tree *node,triangle *T);
+
+// returns the quadrant of the quad_tree containing the point, or 0 if intersects
+// center axes
+int trivial_contain_split_point(quad_tree *node, double xp,double yp);
+
+// returns the quadrant of the quad_tree containing the triangle, or 0 if intersects
+// center axes
+int trivial_contain_split(quad_tree *node, triangle *T);
+
+// returns the triangle in the quad_tree's leaves containing the point or NULL
+// if none of the triangles on the current quad_tree contain it.
+triangle * search_triangles_of_quad_tree(quad_tree * node,double xp,double yp);
+
+// search the tree for a triangle containing the given point. returns the triangle
+// if found, and NULL otherwise
+triangle * search(quad_tree * node ,double xp, double yp);
+
+// return number of noes in tree
+int quad_tree_node_count(quad_tree * tree);
+
+// split the node to make 4 children
+void quad_tree_make_children(quad_tree *node);
+
+// add a triangle to the nodes leaves
+void quad_tree_add_triangle_to_list(quad_tree *node,triangle *T);
+
+
+//**************************************************************************
+
+//**************************************************************************
+// quad_tree_ll struct - Used for serialising  
+//  
+
+typedef struct quad_tree_ll {
+
+    void * tree;
+    struct quad_tree_ll * next;
+    int index;
+
+} quad_tree_ll;
+
+quad_tree_ll * new_quad_tree_ll(quad_tree * start,int index);
+
+//**************************************************************************
+
+//**************************************************************************
+// queue_ll struct - Used for deserialising  
+//  
+
+typedef struct queue_ll {
+
+    int node;
+    struct queue_ll * next;
+
+
+} queue_ll;
+
+queue_ll * new_queue_ll(int node);
+
+//**************************************************************************
+
+#endif
Index: trunk/anuga_core/source/anuga/utilities/quad_tree_ext.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/quad_tree_ext.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/quad_tree_ext.c	(revision 8691)
@@ -0,0 +1,265 @@
+
+#include "Python.h"
+
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+#include "numpy/arrayobject.h"
+#include "math.h"
+
+#include "util_ext.h" /* in utilities */
+#include "quad_tree.h"
+
+
+
+
+static int _serialise(quad_tree * quadtree, PyObject * serial_quadtree)
+{
+
+    
+    quad_tree_ll * nodelist = new_quad_tree_ll(quadtree,0);
+    quad_tree_ll * last = nodelist;
+    quad_tree_ll * temp;
+    int i;
+    int nlist = 0;
+    while(nodelist !=NULL){
+        
+        PyObject * nodedata = PyList_New(2);
+        quadtree=nodelist->tree;
+        // if children have been added, add to the linked list
+        PyObject * child_index;
+        if (quadtree->q[0]!=NULL){
+            child_index = PyList_New(4);
+            for (i=0;i<4;i++){
+                nlist++;
+                quad_tree_ll * child = new_quad_tree_ll(quadtree->q[i],nlist);
+                last->next=child;
+                last=child;
+                PyList_SET_ITEM(child_index,i,PyInt_FromLong((long) nlist));
+            }
+        }
+        else{
+            child_index = PyList_New(1);
+            PyList_SET_ITEM(child_index,0,PyInt_FromLong(0));
+        }
+        // get list of triangle indicies
+        triangle * T = quadtree->leaves;
+        int n_tris = 0 ;
+        while(T!=NULL)
+        {
+            n_tris+=1;
+            T=T->next;
+        }
+        PyObject * tri_list = PyList_New(n_tris);
+        T = quadtree->leaves;
+        i=0;
+        while(T!=NULL)
+        {
+            PyList_SET_ITEM(tri_list,i,PyInt_FromLong((long)T->index));
+            T=T->next;
+            i++;
+        }
+
+        // set data
+        PyList_SET_ITEM(nodedata,0,child_index); 
+        PyList_SET_ITEM(nodedata,1,tri_list);
+        // add data to overall list
+        PyList_SET_ITEM(serial_quadtree, nodelist->index, nodedata);
+        
+        // set up next quad_tree to visit
+        temp = nodelist;
+        nodelist=nodelist->next;
+        free(temp);
+    }
+
+    return 0;
+}
+
+static void _deserialise_recursive(quad_tree * quadtree, PyListObject * serial_quadtree,
+        long * triangles,
+        double *vertex_coordinates,
+        long n)
+{
+
+    int k,k6;
+    double x0,y0,x1,y1,x2,y2;
+    PyObject * node = (PyObject*) PyList_GET_ITEM(serial_quadtree,n);
+    PyObject * children = (PyObject*) PyList_GET_ITEM(node,0);
+    PyObject * leaves = (PyObject*) PyList_GET_ITEM(node,1);
+
+    if((long)PyList_Size(children) == 4){
+        // has children
+        quad_tree_make_children(quadtree);
+        int i;
+        for(i=0;i<4;i++){
+            _deserialise_recursive(quadtree->q[i], serial_quadtree,
+                triangles,
+                vertex_coordinates,PyInt_AsLong(PyList_GET_ITEM(children,i)));
+        }
+    }
+
+    // add leaves
+    long n_leaves = (long) PyList_Size(leaves);
+    int i;
+    for(i=0;i<n_leaves;i++){
+        k = (int) PyInt_AsLong(PyList_GET_ITEM(leaves,i));
+        k6=k*6;
+        x0 = vertex_coordinates[k6];
+        y0 = vertex_coordinates[k6 + 1];
+        x1 = vertex_coordinates[k6 + 2];
+        y1 = vertex_coordinates[k6 + 3];
+        x2 = vertex_coordinates[k6 + 4];
+        y2 = vertex_coordinates[k6 + 5];
+        triangle * T = new_triangle(k,x0,y0,x1,y1,x2,y2);
+        quad_tree_add_triangle_to_list(quadtree,T);
+
+    }
+};
+
+static int _deserialise(quad_tree * quadtree, PyListObject * serial_quadtree,
+        long * triangles,
+        double *vertex_coordinates)
+
+{
+
+    int k,k6;
+    double x0,y0,x1,y1,x2,y2;
+    long n = 0;
+    // if the root node has children, build children, and then call recursive
+    // build on them
+    PyObject * node = (PyObject*) PyList_GET_ITEM(serial_quadtree,n);
+    PyObject * children = (PyObject*) PyList_GET_ITEM(node,0);
+    PyObject * leaves = (PyObject*) PyList_GET_ITEM(node,1);
+
+    if((long)PyList_Size(children) == 4){
+        // has children
+        quad_tree_make_children(quadtree);
+        int i;
+        for(i=0;i<4;i++){
+            _deserialise_recursive(quadtree->q[i], serial_quadtree,
+                triangles,
+                vertex_coordinates,PyInt_AsLong(PyList_GET_ITEM(children,i)));
+        }
+    }
+
+    // add leaves
+    long n_leaves = (long) PyList_Size(leaves);
+    int i;
+    for(i=0;i<n_leaves;i++){
+        k = (int) PyInt_AsLong(PyList_GET_ITEM(leaves,i));
+        k6=k*6;
+        x0 = vertex_coordinates[k6];
+        y0 = vertex_coordinates[k6 + 1];
+        x1 = vertex_coordinates[k6 + 2];
+        y1 = vertex_coordinates[k6 + 3];
+        x2 = vertex_coordinates[k6 + 4];
+        y2 = vertex_coordinates[k6 + 5];
+        triangle * T = new_triangle(k,x0,y0,x1,y1,x2,y2);
+        quad_tree_add_triangle_to_list(quadtree,T);
+    } 
+
+    return 0;
+
+}
+
+
+
+static void delete_quad_tree_cap(PyObject * cap){
+
+    quad_tree * kill = (quad_tree*) PyCapsule_GetPointer(cap,"quad tree");
+    if(kill!=NULL){
+        delete_quad_tree(kill);
+    }
+}
+
+//----------------------- PYTHON WRAPPER FUNCTION -----------------------------
+
+
+static PyObject *serialise(PyObject *self, PyObject *args) {
+
+    int err;
+    PyObject *tree;
+    PyObject *serial_quadtree;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "O",&tree
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "quad_tree_ext.serialise: could not parse input");
+      return NULL;
+    }
+
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    
+    int nodes = quad_tree_node_count(quadtree);
+    serial_quadtree = PyList_New(nodes);
+
+    err = _serialise(quadtree,serial_quadtree);
+
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "quad_tree_ext.serialise: error in serialising quad tree");
+      return NULL;
+    }
+
+    return serial_quadtree;
+
+}
+
+static PyObject *deserialise(PyObject *self, PyObject *args) {
+
+    int err;
+    PyArrayObject *triangles;
+    PyArrayObject *vertex_coordinates;
+    PyArrayObject *extents;
+    PyListObject *serial_quadtree;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOOO",&serial_quadtree,
+                                            &triangles,
+                                            &vertex_coordinates,
+                                            &extents
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "quad_tree_ext.deserialise: could not parse input");
+      return NULL;
+    }
+    double * extents_v = (double *) extents->data;
+    quad_tree *quadtree = new_quad_tree(extents_v[0],extents_v[1],extents_v[2],extents_v[3]);
+    err = _deserialise(quadtree, serial_quadtree,
+        (long*) triangles->data,
+        (double*) vertex_coordinates -> data);
+
+    return  PyCapsule_New((void*) quadtree,
+                      "quad tree",
+                      &delete_quad_tree_cap); 
+
+}
+
+
+//------------------------------------------------------------------------------
+
+
+// ------------------------------ PYTHON GLUE ----------------------------------
+
+//==============================================================================
+// Structures to allow calling from python
+//==============================================================================
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+    {"serialise",serialise, METH_VARARGS, "Print out"},
+    {"deserialise",deserialise, METH_VARARGS, "Print out"},
+	{NULL, NULL, 0, NULL}   // sentinel
+};
+
+
+// Module initialisation
+void initquad_tree_ext(void){
+  
+  Py_InitModule("quad_tree_ext", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+
+}
+
+// --------------------------------------------------------------------------------
Index: trunk/anuga_core/source/anuga/utilities/sparse_csr.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/sparse_csr.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/sparse_csr.c	(revision 8691)
@@ -0,0 +1,37 @@
+#include "sparse_csr.h"
+
+// **************** UTILITIES ***********************
+
+static void *emalloc(size_t amt,char * location)
+{
+    void *v = malloc(amt);  
+    if(!v){
+        fprintf(stderr, "out of mem in quad_tree: %s\n",location);
+        exit(EXIT_FAILURE);
+    }
+    return v;
+};
+
+// ***************************************************
+
+// 'Constructor'
+sparse_csr * make_csr(){
+
+    sparse_csr * ret = emalloc(sizeof(sparse_csr),"make_csr");
+    ret->data=NULL;
+	ret->colind=NULL;
+	ret->row_ptr=NULL;
+	ret->num_rows=0;
+	ret->num_entries=0;
+    return ret;
+}
+
+void delete_csr_matrix(sparse_csr * mat){
+
+	free(mat->data);
+	free(mat->colind);
+	free(mat->row_ptr);
+	free(mat);
+	mat=NULL;
+
+}
Index: trunk/anuga_core/source/anuga/utilities/sparse_csr.h
===================================================================
--- trunk/anuga_core/source/anuga/utilities/sparse_csr.h	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/sparse_csr.h	(revision 8691)
@@ -0,0 +1,42 @@
+/*
+* Sparse Matrix class implement in CSR format. Limited functionality.
+*
+* The CSR (compressed sparse row) format is implemented by using thee
+* double arrays, storing the data, column indicies and row pointers.
+* This format is used for fast matrix-matrix and matrix-vector 
+* multiplication
+*
+* This function is only currently used to store the result of a 
+* conversion from a DOK sparse matrix. Further functionality planned.
+*
+* Padarn Wilson, ANU 2012
+*/
+
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+#include "math.h"
+
+#ifndef SPARSE_CSR_H
+#define SPARSE_CSR_H
+
+// Struct sparse_csr, stores the three matricies which represent
+// the sparse matrix along with aiddtional information about the 
+// number of rows and the number of entries in the matrix.
+typedef struct {
+	double *data;
+	int *colind;
+	int *row_ptr;
+	int num_rows;
+	int num_entries;
+} sparse_csr;
+
+// 'Constructor' function. Returns a pointer to new malloced memory
+// All struct entries are intialised appropriately (mostly to NULL). 
+sparse_csr * make_csr();
+
+// delete_csr_contents - Free the memory associated with the struct
+// and set the pointer to NULL
+void delete_csr_matrix(sparse_csr * mat);
+
+#endif
Index: trunk/anuga_core/source/anuga/utilities/sparse_dok.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/sparse_dok.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/sparse_dok.c	(revision 8691)
@@ -0,0 +1,213 @@
+#include "sparse_dok.h"
+
+
+// **************** UTILITIES ***********************
+
+static void *emalloc(size_t amt,char * location)
+{
+    void *v = malloc(amt);  
+    if(!v){
+        fprintf(stderr, "out of mem in quad_tree: %s\n",location);
+        exit(EXIT_FAILURE);
+    }
+    return v;
+};
+
+// ***************************************************
+
+// 'Constructor'
+
+sparse_dok * make_dok(){
+
+    sparse_dok * ret = emalloc(sizeof(sparse_dok),"make_dok");
+    ret->edgetable=NULL;
+    ret->num_entries=0;
+    ret->num_rows=0;
+    return ret;
+}
+
+//-----------------------------------------------
+//            DOK HASH FUNCTIONS
+//-----------------------------------------------
+
+edge_t *find_dok_entry(sparse_dok * hashtable,edge_key_t key) {
+    edge_t *s;
+    HASH_FIND(hh, hashtable->edgetable, &key, sizeof(edge_key_t), s);  /* s: output pointer */
+    return s;
+}
+
+void add_dok_entry(sparse_dok * hashtable, edge_key_t key, double value) {
+
+    // Checking here now if there is an existing value
+    // not sure if this code will work.
+    
+    edge_t *s;
+    s = find_dok_entry(hashtable,key);
+    if (s) {
+        s->entry+=value;
+    } else {
+        hashtable->num_entries+=1;
+        if(hashtable->num_rows<key.i){
+            hashtable->num_rows = key.i;
+        }
+        s = (edge_t*) emalloc(sizeof(edge_t),"add_dok_entry");
+        memset(s, 0, sizeof(edge_t));
+        s->key.i = key.i;
+        s->key.j = key.j;
+        s->entry = value;
+        HASH_ADD(hh, hashtable->edgetable, key, sizeof(edge_key_t), s);  /* key: name of key field */    
+    }
+    if(s->entry <= 1e-10 && s->entry >= -1e-10){
+            HASH_DEL(hashtable->edgetable,s);
+            free(s);
+            hashtable->num_entries-=1;
+    }
+
+    
+}
+
+void delete_dok_entry(sparse_dok * hashtable,edge_t *edge) {
+    HASH_DEL( hashtable->edgetable, edge);  /* user: pointer to deletee */
+    free(edge);
+}
+
+void delete_all(sparse_dok * hashtable) {
+  edge_t *current_edge, *tmp;
+
+  HASH_ITER(hh, hashtable->edgetable, current_edge, tmp) {
+    HASH_DEL(hashtable->edgetable, current_edge);  /* delete it (hashtable advances to next) */
+    free(current_edge);            /* free it */
+  } 
+}
+
+void delete_dok_matrix(sparse_dok * mat) {
+
+    delete_all(mat);
+    free(mat->edgetable);
+    free(mat);
+    mat=NULL;
+
+}
+
+void print_dok_entries(sparse_dok * hashtable) {
+    edge_t *s;
+
+    for(s=hashtable->edgetable; s != NULL; s=(edge_t*)(s->hh.next)) {
+        printf("edge key i %d i %d entry %f\n",
+                      s->key.i, s->key.j, s->entry);
+    }
+}
+
+int key_sort(edge_t *a, edge_t *b) {
+    return (a->key.i - b->key.i);
+}
+
+int key_sort_2(edge_t *a, edge_t *b){
+    if(a->key.i - b->key.i==0){
+        return (a->key.j-b->key.j);
+    } else{
+        return (a->key.i-b->key.i);
+    }
+}
+
+void sort_by_key(sparse_dok * hashtable) {     
+    HASH_SORT(hashtable->edgetable, key_sort_2);
+}
+
+//----------------END--------------------
+
+//-----------------------------------------------
+//            DOK MATRIX FUNCTIONS
+//-----------------------------------------------
+
+
+void convert_to_csr_ptr(sparse_csr * new_csr, sparse_dok * hashtable){
+
+
+    sparse_csr * ret_csr = new_csr;
+
+    //entrires stores in edgetable -> end
+
+    //sort and get number of entries
+    sort_by_key(hashtable); 
+    int num_entries = hashtable->num_entries;
+    int num_rows = hashtable->num_rows+1;
+
+    //build storage matricies
+    ret_csr->data=emalloc(num_entries*sizeof(double),"convert_to_csr_ptr");
+    ret_csr->colind=emalloc(num_entries*sizeof(int),"convert_to_csr_ptr");
+    ret_csr->row_ptr=emalloc((num_rows+1)*sizeof(int),"convert_to_csr_ptr");
+
+    edge_t * edge = hashtable->edgetable;
+
+    //now convert
+    int current_row = -1;
+    int k;
+    for(k=0;k<num_entries;k++){
+        int i = edge->key.i;
+        int j = edge->key.j;
+        double value = edge->entry;
+
+        if (i!=current_row){
+            current_row=i;
+            ret_csr->row_ptr[i]=k;
+        }
+
+        ret_csr->data[k] = value;
+        ret_csr->colind[k] = j;
+        edge = edge->hh.next;
+    }
+
+    for(k=current_row+1;k<num_rows+1;k++){
+        ret_csr->row_ptr[k]=num_entries;
+    }
+
+    ret_csr -> num_rows = num_rows+1;
+    ret_csr -> num_entries = num_entries;
+
+}
+
+void add_sparse_dok(sparse_dok * dok1,double mult1,sparse_dok * dok2,double mult2){
+
+    // add both into dok1 - then leave both alone (free outside)
+    int num_entries = dok1->num_entries;
+    edge_t * edge = dok1->edgetable;
+    edge_t * edge2;
+
+    int k;
+    for(k=0;k<num_entries;k++){
+        int i = edge->key.i;
+        int j = edge->key.j;
+        double value = edge->entry;
+        edge->entry=value*mult1;
+        edge2=find_dok_entry(dok2,edge->key);
+        if(edge2!=NULL){
+            edge->entry+=edge2->entry*mult2;
+            edge2->entry=0;
+            //delete_dok_entry(dok2,edge2);
+        }
+        edge = edge->hh.next;
+    }
+
+    num_entries = dok2->num_entries;
+    edge = dok2->edgetable;
+    for(k=0;k<num_entries;k++){
+        add_dok_entry(dok1,edge->key,edge->entry*mult2);
+        edge = edge->hh.next;
+    }
+
+}
+
+int get_dok_rows(sparse_dok * dok){
+
+    int rows = 0;
+
+    edge_t *current_edge, *tmp;
+
+    HASH_ITER(hh, dok->edgetable, current_edge, tmp) {
+        if (current_edge->key.i>rows) rows = current_edge->key.i;
+    } 
+
+
+    return rows;
+}
Index: trunk/anuga_core/source/anuga/utilities/sparse_dok.h
===================================================================
--- trunk/anuga_core/source/anuga/utilities/sparse_dok.h	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/sparse_dok.h	(revision 8691)
@@ -0,0 +1,110 @@
+/*
+* Sparse Matrix class implemented in DOK format. 
+*
+* The DOK (dictionary of keys) format is implemented using 
+* uthash, which is a doubly linked list and also a hash table.
+* The hash table is populated by 'edge_t' objects which contain
+* a key specifying the (i,j) index of the matrix, and the value
+* at that index. If no entry exists, the value is assumed 0.
+* 
+* The functions in this class which create new objects return
+* pointers to new malloced memory, which much be cleaned up.
+* 
+* Padarn Wilson, ANU 2012
+*/
+
+
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+#include "math.h"
+#include "uthash.h"     /* in utilities */
+#include "sparse_csr.h"
+
+#ifndef SPARSE_DOK_H
+#define SPARSE_DOK_H
+
+// Struct edge_key_t to store the i,j position in the matrix within
+// a key of the hashtable
+typedef struct  {
+    int i;
+    int j;
+} edge_key_t;
+
+// Struct edge_t is a basic element of the hash table. By including 
+// the UT_hash_handle variable, it can be used in the methods defined
+// for uthash hashtables
+typedef struct {
+    edge_key_t key;              /* key of form i , j */
+    double entry;                /* the value stored at this entry */
+    UT_hash_handle hh;         /* makes this structure hashable */
+} edge_t;
+
+// Struct sparse_dok defining a sparse matrix. Keeps track of the 
+// number of entries and rows in the matrix, and stores the hashtable.
+// PADARN NOTE: For efficiency, it might be sensible to pre-allocated
+// the max number of rows/cols in the hash table, so that the hash table
+// can be made an appropriate size.
+typedef struct {
+	edge_t *edgetable;
+	int num_entries;
+	int num_rows;
+} sparse_dok;
+
+// 'Constructor' function. Returns pointer to new malloced memory, with 
+// appropriate initilisation.
+sparse_dok * make_dok();
+
+// --------------- Hashtable Functions -----------------
+
+// find_dok_entry - Find pointer to hash table member with given key, 
+// return NULL if member does not exist. 
+edge_t *find_dok_entry(sparse_dok * edgetable,edge_key_t key);
+
+// add_dok_entry - Add new entry to the hash table with given key, holding
+// the specified value. If entry already exists value is added to current
+// value. If entry already exists and this function causes the value to become
+// zero, the entry is removed.
+void add_dok_entry(sparse_dok * edgetable,edge_key_t key, double value);
+
+// delete_dok_entry - Remove an edge from the hash table. Implicitly sets
+// the corresponding entry in the matrix to zero.
+void delete_dok_entry(sparse_dok * edgetable,edge_t *edge);
+
+// delete_all - Remove all edges from the hash table. Used to do clean up
+void delete_all(sparse_dok * edgetable);
+
+// delete_dok_matrix - Free all the memory associated with struct and
+// set pointer to Null.
+void delete_dok_matrix(sparse_dok * mat);
+
+// print_dok_entries - Print out all of the stored entries in the hash 
+// table sorted by their key, along with their value.
+void print_dok_entries(sparse_dok * edgetable);
+
+// key_sort - Compare the relative size of two keys, used for sorting
+// PADARN NOTE: Does not need to be in header.
+int key_sort(edge_t *a, edge_t *b);
+
+// sort_by_key - Sort the linked list of the hash table by their key
+// values and the key_sort function.
+void sort_by_key(sparse_dok * hashtable);
+
+// --------------- Matrix Functions ---------------------
+
+// convert_to_csr_ptr - Convert the DOK format matrix into CSR format. The
+// new matrix is stored in the (already allocated) input pointer. The old 
+// pointer is not freed, and should be cleaned up manually.
+void convert_to_csr_ptr(sparse_csr * new_csr,sparse_dok * hashtable);
+
+// add_sparse_dok - Perform a linear addition on two dok_matricies A=a*A+b*B.
+// The result is stored in the first input sparse_dok. The second sparse_dok is 
+// not freed and should be cleaned up manually. 
+// Note: No size conditions are imposed on the matricies, so potentially different
+// sized matricies can be combined.
+void add_sparse_dok(sparse_dok * dok1,double mult1,sparse_dok * dok2,double mult2);
+
+// get_dok_rows -- Return the number of rows currently stored in the matrix
+int get_dok_rows(sparse_dok * dok);
+
+#endif
Index: trunk/anuga_core/source/anuga/utilities/sparse_matrix_ext.c
===================================================================
--- trunk/anuga_core/source/anuga/utilities/sparse_matrix_ext.c	(revision 8691)
+++ trunk/anuga_core/source/anuga/utilities/sparse_matrix_ext.c	(revision 8691)
@@ -0,0 +1,189 @@
+
+#include "Python.h"
+
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+#include "numpy/arrayobject.h"
+#include "math.h"
+
+#include "util_ext.h" /* in utilities */
+#include "sparse_dok.h"
+
+
+
+
+static int _serialise(sparse_dok * dok, PyObject * serial_dok)
+{
+
+    sort_by_key(dok);
+    int num_entries = dok->num_entries;
+    int k,i,j;
+    double val;
+
+    edge_t * edge = dok->edgetable;
+
+    for(k=0;k<num_entries;k++){
+
+        i = edge->key.i;
+        j = edge->key.j;
+        val = edge->entry;
+
+        // make new python objects to add to dict
+        PyObject * py_val = PyFloat_FromDouble(val);
+        PyObject * py_i = PyInt_FromLong((long) i);
+        PyObject * py_j = PyInt_FromLong((long) j);
+        PyObject * py_key = PyTuple_New(2);
+        PyTuple_SET_ITEM(py_key,0,py_i);
+        PyTuple_SET_ITEM(py_key,1,py_j);
+
+        PyDict_SetItem(serial_dok,py_key,py_val);
+
+        edge = edge->hh.next;
+
+    }
+   
+    return 0;
+};
+
+
+static int _deserialise(sparse_dok * dok, PyObject * serial_dok)
+
+{
+
+    PyObject * items = PyDict_Items(serial_dok);
+    PyObject * keys = PyDict_Keys(serial_dok);
+
+    int num_entries = (int) PyDict_Size(serial_dok);
+    printf("num entries %d\n",num_entries);
+    int k,i,j;
+    double val;
+    PyObject * py_val;
+    PyObject * py_key;
+
+    edge_key_t key;
+
+
+
+    for(k=0;k<num_entries;k++){
+
+        py_key = PyList_GET_ITEM(keys,k);
+        //py_val = PyList_GET_ITEM(items,k);
+        py_val = PyDict_GetItem(serial_dok,py_key);
+
+        val = PyFloat_AS_DOUBLE(py_val);
+        printf("val: %d\n",val);
+
+        i = (int) PyInt_AS_LONG(PyTuple_GET_ITEM(py_key,0));
+        j = (int) PyInt_AS_LONG(PyTuple_GET_ITEM(py_key,1));
+
+        key.i = i;
+        key.j = j;
+
+        add_dok_entry(dok,key,val);
+
+    }
+
+    Py_DECREF(keys);
+    Py_DECREF(items);
+
+    return 0;
+
+};
+
+
+// Delete capsule containing a sparse_dok - name of capsule must be exactly
+// "sparse dok".
+static void delete_sparse_dok_cap(PyObject *cap){
+
+    sparse_dok * kill = (sparse_dok*) PyCapsule_GetPointer(cap,"sparse dok");
+    if(kill!=NULL){
+        delete_dok_matrix(kill);
+    }
+
+}
+//----------------------- PYTHON WRAPPER FUNCTION -----------------------------
+
+
+static PyObject *serialise_dok(PyObject *self, PyObject *args) {
+
+    int err;
+    PyObject *sparse_dok_cap;
+    PyObject *serial_sparse_dok;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "O",&sparse_dok_cap
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "sparse_matrix_ext.serialise_dok: could not parse input");
+      return NULL;
+    }
+
+    sparse_dok * dok = (sparse_dok*) PyCapsule_GetPointer(sparse_dok_cap,"sparse dok");
+
+    serial_sparse_dok = PyDict_New();
+
+    err = _serialise(dok,serial_sparse_dok);
+
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "sparse_matrix.serialise_dok: error in serialising sparse_dok");
+      return NULL;
+    }
+
+    return serial_sparse_dok;
+
+}
+
+static PyObject *deserialise_dok(PyObject *self, PyObject *args) {
+
+    int err;
+    PyObject *sparse_dok_cap;
+    PyObject *serial_sparse_dok;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "O", &serial_sparse_dok
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "sparse_matrix.deserialise_dok: could not parse input");
+      return NULL;
+    }
+    
+    sparse_dok * dok;
+    dok = make_dok();
+
+    err = _deserialise(dok,serial_sparse_dok);
+
+    return  PyCapsule_New((void*) dok,
+                      "sparse dok",
+                      &delete_sparse_dok_cap); 
+
+}
+
+
+//------------------------------------------------------------------------------
+
+
+// ------------------------------ PYTHON GLUE ----------------------------------
+
+//==============================================================================
+// Structures to allow calling from python
+//==============================================================================
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+    {"serialise_dok",serialise_dok, METH_VARARGS, "Print out"},
+    {"deserialise_dok",deserialise_dok, METH_VARARGS, "Print out"},
+	{NULL, NULL, 0, NULL}   // sentinel
+};
+
+
+// Module initialisation
+void initsparse_matrix_ext(void){
+  
+  Py_InitModule("sparse_matrix_ext", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+
+}
+
+// --------------------------------------------------------------------------------
