Index: anuga_work/development/flow_1d/generic_1d/__init__.py
===================================================================
--- anuga_work/development/flow_1d/generic_1d/__init__.py	(revision 7830)
+++ anuga_work/development/flow_1d/generic_1d/__init__.py	(revision 7830)
@@ -0,0 +1,9 @@
+#! /usr/bin/python
+
+# To change this template, choose Tools | Templates
+# and open the template in the editor.
+
+__author__="Stephen Roberts"
+__date__ ="$14/06/2010 7:50:20 PM$"
+
+pass
Index: anuga_work/development/flow_1d/generic_1d/compile_all.py
===================================================================
--- anuga_work/development/flow_1d/generic_1d/compile_all.py	(revision 7827)
+++ 	(revision )
@@ -1,15 +1,0 @@
-import os
-
-buildroot = os.getcwd()
-
-#Attempt to compile all extensions
-
-anuga_dir = '..' + os.sep + '..' + os.sep + '..'
-utilities_dir = anuga_dir + os.sep + 'trunk' + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities'
-
-execfile( utilities_dir + os.sep + 'compile.py')
-
-
-os.chdir(buildroot)    
-
-    
Index: anuga_work/development/flow_1d/generic_1d/quantity.py
===================================================================
--- anuga_work/development/flow_1d/generic_1d/quantity.py	(revision 7827)
+++ anuga_work/development/flow_1d/generic_1d/quantity.py	(revision 7830)
@@ -1126,5 +1126,5 @@
 if __name__ == "__main__":
     #from domain import Domain
-    from sww_domain import Domain
+    from generic_domain import Generic_domain as Domain
     
     points1 = [0.0, 1.0, 2.0, 3.0]
Index: anuga_work/development/flow_1d/generic_1d/test_quantity.py
===================================================================
--- anuga_work/development/flow_1d/generic_1d/test_quantity.py	(revision 7827)
+++ anuga_work/development/flow_1d/generic_1d/test_quantity.py	(revision 7830)
@@ -473,8 +473,7 @@
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
             assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/2.0)
+                             numpy.sum(quantity.vertex_values[k,:])/2.0)
 
 
Index: anuga_work/development/flow_1d/generic_1d/util.py
===================================================================
--- anuga_work/development/flow_1d/generic_1d/util.py	(revision 7827)
+++ 	(revision )
@@ -1,175 +1,0 @@
-"""
-This module contains various auxiliary function used by pyvolution.
-"""
-
-
-import numpy
-
-def mean(x):
-    return numpy.sum(x)/len(x)
-
-
-def gradient(x0, x1, q0, q1):
-
-    if q1-q0 != 0:
-        a = (q1-q0)/(x1-x0)
-    else:
-        a = 0
-        
-    return a
-
-def  minmod(beta_p,beta_m):
-    if (abs(beta_p) < abs(beta_m)) & (beta_p*beta_m > 0.0):
-        phi = beta_p
-    elif (abs(beta_m) < abs(beta_p)) & (beta_p*beta_m > 0.0):
-        phi = beta_m
-    else:
-        phi = 0.0
-    return phi
-
-def  minmod_kurganov(a,b,c):
-    from numpy import sign
-    if sign(a)==sign(b)==sign(c):
-        return sign(a)*min(abs(a),abs(b),abs(c))
-    else:
-        return 0.0
-
-def  maxmod(a,b):
-    if (abs(a) > abs(b)) & (a*b > 0.0):
-        phi = a
-    elif (abs(b) > abs(a)) & (a*b > 0.0):
-        phi = b
-    else:
-        phi = 0.0
-    return phi
-
-def vanleer(a,b):
-    if abs(a)+abs(b) > 1e-12:
-        return (a*abs(b)+abs(a)*b)/(abs(a)+abs(b))
-    else:
-        return 0.0
-
-def vanalbada(a,b):
-    if a*a+b*b > 1e-12:
-        return (a*a*b+a*b*b)/(a*a+b*b)
-    else:
-        return 0.0
-
-def calculate_wetted_area(x1,x2,z1,z2,w1,w2):
-    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
-        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
-        A = 0.5*(w1-z1)*(x-x1)
-        L = x-x1
-    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
-        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
-        A = 0.5*(w2-z2)*(x2-x)
-        L = x2-x
-    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
-        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
-        A = 0.5*(w2-z2)*(x2-x)
-        L = x2-x
-    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
-        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
-        A = 0.5*(w1-z1)*(x-x1)
-        L = x-x1
-    elif (w1 <= z1) & (w2 <= z2):
-        A = 0.0
-    elif (w1 == z1) & (w2 > z2) & (z2 < z1):
-        A = 0.5*(x2-x1)*(w2-z2)
-    elif (w2 == z2) & (w1 > z1) & (z1 < z2):
-        A = 0.5*(x2-x1)*(w1-z1)
-    return A
-
-
-def calculate_new_wet_area(x1,x2,z1,z2,A):
-    from numpy import sqrt
-    min_centroid_height = 1.0e-3
-    # Assumes reconstructed stage flat in a wetted cell
-    M = (z2-z1)/(x2-x1)
-    L = (x2-x1)
-    min_area = min_centroid_height*L
-    max_area = 0.5*(x2-x1)*abs(z2-z1)
-    if A < max_area:
-        if (z1 < z2):
-            x = sqrt(2*A/M)+x1
-            wet_len = x-x1
-            wc = z1 + sqrt(M*2*A)
-        elif (z2 < z1):
-            x = -sqrt(-2*A/M)+x2
-            wet_len = x2-x 
-            wc = z2+sqrt(-M*2*A)
-        else:
-            wc = A/L+0.5*(z1+z2)
-            wet_len = x2-x1
-    else:
-        wc = 0.5*(z1+z2)+A/L
-        wet_len = x2-x1
-            
-    return wc,wet_len
-
-def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):
-    min_centroid_height = 1.0e-3
-    # Assumes reconstructed stage flat in a wetted cell
-    M = (z2-z1)/(x2-x1)
-    L = (x2-x1)
-    min_area = min_centroid_height*L
-    max_area = 0.5*(x2-x1)*abs(z2-z1)
-    w1,uh1 = analytic_cannal(x1,t)
-    w2,uh2 = analytic_cannal(x2,t)
-    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
-        print "test1"
-        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
-        wet_len = x-x1
-    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
-        print "test2"
-        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
-        wet_len = x2-x
-    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
-        print "test3"
-        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
-        wet_len = x2-x
-    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
-        print "test4"
-        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
-        wet_len = x-x1
-    elif (w1 >= z1) & (w2 >= z2):
-        print "test5"
-        wet_len = x2-x1 
-    else: #(w1 <= z1) & (w2 <= z2)
-        print "test5"
-        if (w1 > z1) | (w2 > z2):
-            print "ERROR"
-        wet_len = x2-x1        
-    return w1,w2,wet_len,uh1,uh2
-
-def analytic_cannal(C,t):
-
-    import math
-    #N = len(C)
-    #u = zeros(N,numpy.float)    ## water velocity
-    #h = zeros(N,numpy.float)    ## water depth
-    x = C
-    g = 9.81
-
-
-    ## Define Basin Bathymetry
-    #z_b = zeros(N,numpy.float) ## elevation of basin
-    #z = zeros(N,numpy.float)   ## elevation of water surface
-    z_infty = 10.0       ## max equilibrium water depth at lowest point.
-    L_x = 2500.0         ## width of channel
-
-    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
-    omega = math.sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
-
-    x1 = A0*cos(omega*t)-L_x # left shoreline
-    x2 = A0*cos(omega*t)+L_x # right shoreline
-    if (x >=x1) & (x <= x2):
-        z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
-        u = -A0*omega*sin(omega*t)
-        z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t))
-    else:
-       z_b = z_infty*(x**2/L_x**2)
-       u=0.0
-       z = z_b
-    h = z-z_b
-    return z,u*h
Index: anuga_work/development/flow_1d/generic_1d/util_ext.h
===================================================================
--- anuga_work/development/flow_1d/generic_1d/util_ext.h	(revision 7827)
+++ 	(revision )
@@ -1,367 +1,0 @@
-// Python - C extension for finite_volumes util module.
-//
-// To compile (Python2.3):
-//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
-//  gcc -shared util_ext.o  -o util_ext.so	
-//
-// See the module util.py
-//
-//
-// Ole Nielsen, GA 2004
-	
-#include "Python.h"	
-#include "numpy/arrayobject.h"
-#include "math.h"
-
-
-double max(double x, double y) {  
-  //Return maximum of two doubles
-  
-  if (x > y) return x;
-  else return y;
-}
-
-
-double min(double x, double y) {  
-  //Return minimum of two doubles
-  
-  if (x < y) return x;
-  else return y;
-}
-
-
-int _gradient(double x0, double y0, 
-	      double x1, double y1, 
-	      double x2, double y2, 
-	      double q0, double q1, double q2, 
-	      double *a, double *b) {
-	      
-  /*Compute gradient (a,b) based on three points (x0,y0), (x1,y1) and (x2,y2) 
-  with values q0, q1 and q2.
-  
-  Extrapolation formula (q0 is selected as an arbitrary origin)
-    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
-  
-  Substituting the known values for q1 and q2 into (1) yield the 
-  equations for a and b 
-  
-      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
-      q2-q0 = a*(x2-x0) + b*(y2-y0)                      (3)      
-      
-  or in matrix form
-  
-  /               \  /   \   /       \  
-  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
-  |               |  |   | = |       | 
-  |  x2-x0  y2-y0 |  | b |   | q2-q0 |
-  \               /  \   /   \       /
-   
-  which is solved using the standard determinant technique    
-      
-  */
-	      
-
-  double det;
-  
-  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
-
-  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
-  *a /= det;
-
-  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
-  *b /= det;
-
-  return 0;
-}
-
-
-int _gradient2(double x0, double y0, 
-	       double x1, double y1, 
-	       double q0, double q1, 
-	       double *a, double *b) {
-  /*Compute gradient (a,b) between two points (x0,y0) and (x1,y1) 
-  with values q0 and q1 such that the plane is constant in the direction 
-  orthogonal to (x1-x0, y1-y0).
-  
-  Extrapolation formula
-    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
-  
-  Substituting the known values for q1 into (1) yields an 
-  under determined  equation for a and b 
-      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
-      
-      
-  Now add the additional requirement that the gradient in the direction 
-  orthogonal to (x1-x0, y1-y0) should be zero. The orthogonal direction 
-  is given by the vector (y0-y1, x1-x0).
-  
-  Define the point (x2, y2) = (x0 + y0-y1, y0 + x1-x0) on the orthognal line. 
-  Then we know that the corresponding value q2 should be equal to q0 in order 
-  to obtain the zero gradient, hence applying (1) again    
-      q0 = q2 = q(x2, y2) = q0 + a*(x2-x0) + b*(y2-y0)
-                          = q0 + a*(x0 + y0-y1-x0) + b*(y0 + x1-x0 - y0)
-			  = q0 + a*(y0-y1) + b*(x1-x0)
-			  
-  leads to the orthogonality constraint
-     a*(y0-y1) + b*(x1-x0) = 0                           (3) 
-     
-  which closes the system and yields
-  
-  /               \  /   \   /       \  
-  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
-  |               |  |   | = |       | 
-  |  y0-y1  x1-x0 |  | b |   |   0   |
-  \               /  \   /   \       /
-   
-  which is solved using the standard determinant technique    
-      
-  */
-
-  double det, xx, yy, qq;
-  
-  xx = x1-x0;
-  yy = y1-y0;
-  qq = q1-q0;
-    
-  det = xx*xx + yy*yy;  //FIXME  catch det == 0
-  *a = xx*qq/det;
-  *b = yy*qq/det;
-        
-  return 0;
-}
-
-
-void _limit_old(int N, double beta, double* qc, double* qv, 
-	    double* qmin, double* qmax) { 
-
-  //N are the number of elements
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  //printf("INSIDE\n");
-  for (k=0; k<N; k++) {
-    k3 = k*3;
-    
-    //Find the gradient limiter (phi) across vertices  
-    phi = 1.0;
-    for (i=0; i<3; i++) {    
-      r = 1.0;
-      
-      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
-      dqa[i] = dq;              //Save dq for use in the next loop
-      
-      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
-      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
-  
-  
-      phi = min( min(r*beta, 1.0), phi);    
-    }
-    
-    //Then update using phi limiter
-    for (i=0; i<3; i++) {    
-      qv[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}
-
-
-void  print_double_array(char* name, double* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-	km = m*k;
-	printf("[");
-	for (i=0; i<m ; i++){
-	    printf("%g ",array[km+i]);
-	}
-	if (k==(n-1))
-	    printf("]");
-	else
-	    printf("]\n");
-    }
-    printf("]\n");
-}
-
-void  print_int_array(char* name, int* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-	km = m*k;
-	printf("[");
-	for (i=0; i<m ; i++){
-	    printf("%i ",array[km+i]);
-	}
-	if (k==(n-1))
-	    printf("]");
-	else
-	    printf("]\n");
-    }
-    printf("]\n");
-}
-
-
-void  print_long_array(char* name, long* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-	km = m*k;
-	printf("[");
-	for (i=0; i<m ; i++){
-	  printf("%i ",(int) array[km+i]);
-	}
-	if (k==(n-1))
-	    printf("]");
-	else
-	    printf("]\n");
-    }
-    printf("]\n");
-}
-
-void print_numeric_array(PyArrayObject *x) {  
-  int i, j;
-  for (i=0; i<x->dimensions[0]; i++) {  
-    for (j=0; j<x->dimensions[1]; j++) {
-      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
-    }
-    printf("\n");  
-  }
-  printf("\n");    
-}
-
-void print_numeric_vector(PyArrayObject *x) {  
-  int i;
-  for (i=0; i<x->dimensions[0]; i++) {
-    printf("%f ", *(double*) (x->data + i*x->strides[0]));  
-  }
-  printf("\n");  
-}
-
-PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
-  PyArrayObject *A, *B;
-  
-
-  //Get array object from attribute
-  
-  /*
-  //FIXME: THE TEST DOESN't WORK
-  printf("Err = %d\n", PyObject_HasAttrString(O, name));
-  if (PyObject_HasAttrString(O, name) == 1) {
-    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
-    if (!B) return NULL;
-  } else {
-    return NULL;
-    } 
-  */
-    
-  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
-
-  //printf("B = %p\n",(void*)B);
-  if (!B) {
-    printf("util_ext.h: get_consecutive_array could not obtain python object");
-    printf(" %s\n",name);
-    fflush(stdout);
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain python object");
-    return NULL;
-  }     
-  
-  //Convert to consecutive array
-  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
-						    B -> descr -> type, 0, 0);
-  
-  Py_DECREF(B); //FIXME: Is this really needed??
-  
-  if (!A) {
-    printf("util_ext.h: get_consecutive_array could not obtain array object");
-    printf(" %s \n",name);
-    fflush(stdout);
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
-    return NULL;
-  }
-
-
-  return A;
-}
-
-double get_python_double(PyObject *O, char *name) {
-  PyObject *TObject;
-  #define BUFFER_SIZE 80
-  char buf[BUFFER_SIZE];
-  double tmp;
-  int n;
-  
-
-  //Get double from attribute
-  TObject = PyObject_GetAttrString(O, name);
-  if (!TObject) {
-	n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
-	//printf("name = %s",name);
-    PyErr_SetString(PyExc_RuntimeError, buf);
-
-    return 0.0;
-  }  
-  
-  tmp = PyFloat_AsDouble(TObject);
-  
-  Py_DECREF(TObject);
-  
-  return tmp;
-}
-
-
-
-
-int get_python_integer(PyObject *O, char *name) {
-  PyObject *TObject;
-  #define BUFFER_SIZE 80
-  char buf[BUFFER_SIZE];
-  long tmp;
-  int n;
-  
-
-  //Get double from attribute
-  TObject = PyObject_GetAttrString(O, name);
-  if (!TObject) {
-  	n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
-	//printf("name = %s",name);
-    PyErr_SetString(PyExc_RuntimeError, buf);
-    return 0;
-  }  
-  
-  tmp = PyInt_AsLong(TObject);
-  
-  Py_DECREF(TObject);
-  
-  return tmp;
-}
-
-
-PyObject *get_python_object(PyObject *O, char *name) {
-  PyObject *Oout;
-
-  Oout = PyObject_GetAttrString(O, name);
-  if (!Oout) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
-    return NULL;
-  }
-
-  return Oout;
-
-}
-
-
-// check that numpy array objects are C contiguous memory
-#define CHECK_C_CONTIG(varname)	if (!PyArray_ISCONTIGUOUS(varname)) { \
-				    char msg[1024]; \
-				    sprintf(msg, \
-					    "%s(): file %s, line %d: " \
-				            "'%s' object is not C contiguous memory", \
-				             __func__, __FILE__, __LINE__, #varname); \
-				    PyErr_SetString(PyExc_RuntimeError, msg); \
-				    return NULL; \
-				}
