Index: /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py	(revision 8253)
+++ /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py	(revision 8254)
@@ -107,5 +107,5 @@
     ztplot, = pylab.plot(x, z+t)
 
-    plot1.set_ylim([-1,40])
+    plot1.set_ylim([-1,50])
     pylab.xlabel('Position')
     pylab.ylabel('Stage')
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/solver.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/solver.h	(revision 8253)
+++ 	(revision )
@@ -1,396 +1,0 @@
-/**********************************************************
-Generic params class
-**********************************************************/
-struct params;
-
-// vtable
-struct params_vtable {  
-  struct params* (*destroy) (struct params *q);
-};
-
-// structure
-struct params {
-  struct params_vtable *vtable;
-  
-  double epsilon;
-  double cfl;
-};
-
-// Public API
-
-// Destruction
-struct params* params_destroy(struct params *q);
-
-// method declarations
-
-// method implementations
-struct params* params_destroy(struct params *q) {
-  return q->vtable->destroy(q);
-}
-/**********************************************************/
-
-
-/**********************************************************
-Generic quantity class
-**********************************************************/
-struct quantity;
-
-// vtable
-struct quantity_vtable {
-  double* (*flux_formula) (struct quantity *q, double normal, struct params *p, double *quantityflux);
-  double (*sound_speed) (struct quantity *q, struct params *p);
-  double (*get_conserved) (struct quantity *q, int k, double normal);
-  struct quantity* (*destroy) (struct quantity *q);
-};
-
-// structure
-struct quantity {
-  struct quantity_vtable *vtable;
-};
-
-// Public API
-// Destruction
-struct quantity* quantity_destroy(struct quantity *q);
-
-// method declarations
-double* quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
-double quantity_sound_speed (struct quantity *q, struct params *p);
-double quantity_get_conserved (struct quantity *q, int k, double normal);
-
-// method implementations
-double* quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
-  return q->vtable->flux_formula(q, normal, p, quantityflux);
-}
-
-double quantity_sound_speed (struct quantity *q, struct params *p) {
-  return q->vtable->sound_speed(q, p);
-}
-
-double quantity_get_conserved (struct quantity *q, int k, double normal) {
-  return q->vtable->get_conserved(q, k, normal);
-}
-
-struct quantity* quantity_destroy(struct quantity *q) {
-  return q->vtable->destroy(q);
-}
-/**********************************************************/
-
-
-/**********************************************************
-Generic edge class
-**********************************************************/
-struct edge;
-
-// vtable
-struct edge_vtable{
-  double* (*flux_function) (struct edge *e, struct params *p, double *edgeflux, int num_eqns);
-  double (*extreme_sound_speeds) (struct edge *e, struct params *p);
-  struct edge* (*destroy) (struct edge *e);
-};
-
-// structure
-struct edge {
-  struct edge_vtable *vtable;
-  struct quantity *qin, *qout;
-  double normal;
-  double smin, smax;
-};
-
-// Public API
-// Destruction 
-struct edge* edge_destroy(struct edge *e);
-
-// method declarations
-double* edge_flux_function (struct edge *e, struct params *p, double *edgeflux, int num_eqns);
-double edge_extreme_sound_speeds (struct edge *e, struct params *p);
-
-// method implementations
-double* edge_flux_function (struct edge *e, struct params *p, double *edgeflux, int num_eqns) {
-  return e->vtable->flux_function(e, p, edgeflux, num_eqns);
-}
-
-double edge_extreme_sound_speeds (struct edge *e, struct params *p) {
-  return e->vtable->extreme_sound_speeds(e, p);
-}
-
-struct edge* edge_destroy(struct edge *e) {
-  return e->vtable->destroy(e);
-}
-/**********************************************************/
-
-
-/**********************************************************
-Generic cell class
-**********************************************************/
-struct cell;
-
-// vtable
-struct cell_vtable{
-  double* (*flux_function) (struct cell *c, struct params *p, double *cellflux, int num_eqns);
-  double (*extreme_sound_speeds) (struct cell *c, struct params *p);
-  double* (*forcing_terms) (struct cell *c, struct params *p, double *cellforce);
-  struct cell* (*destroy) (struct cell *c);
-};
-
-// structure
-struct cell {
-  struct cell_vtable *vtable;
-  struct edge **edges;
-  int num_edges;
-};
-
-// Public API
-
-// Destruction
-struct cell* cell_destroy(struct cell *c);
-
-// method declarations
-double* cell_flux_function (struct cell *c, struct params *p, double *cellflux, int num_eqns);
-double cell_extreme_sound_speeds (struct cell *c, struct params *p);
-double* cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
-
-// method implementations
-double* cell_flux_function (struct cell *c, struct params *p, double *cellflux, int num_eqns) {
-  return c->vtable->flux_function(c, p, cellflux, num_eqns);
-}
-
-double cell_extreme_sound_speeds (struct cell *c, struct params *p) {
-  return c->vtable->extreme_sound_speeds(c, p);
-}
-
-double* cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
-  return c->vtable->forcing_terms(c, p, cellforce);
-}
-
-struct cell* cell_destroy(struct cell *c) {
-  return c->vtable->destroy(c);
-}
-/**********************************************************/
-
-
-/**********************************************************
-Generic quantities class
-**********************************************************/
-struct quantities;
-
-// vtable
-struct quantities_vtable{
-  struct quantity* (*get_quantity) (struct quantities *qs, int i, struct quantity* q);
-  void (*update) (struct quantities *qs, double *flux, int k);
-  struct quantities* (*destroy) (struct quantities *qs);
-};
-
-// structure
-struct quantities {
-  struct quantities_vtable *vtable;
-};
-
-// Public API
-
-// Destruction 
-struct quantities* quantities_destroy(struct quantities *qs);
-
-// method declarations
-struct quantity* quantities_get_quantity(struct quantities *qs, int i, struct quantity *q); 
-void quantities_update (struct quantities *qs, double *flux, int k);
-
-// method implementations
-struct quantity* quantities_get_quantity(struct quantities *qs, int i, struct quantity *q) {
-  return qs->vtable->get_quantity(qs, i, q);
-}
-
-void quantities_update (struct quantities *qs, double *flux, int k) {
-  return qs->vtable->update(qs, flux, k);
-}
-
-struct quantities* quantities_destroy(struct quantities *qs) {
-  return qs->vtable->destroy(qs);
-}
-/**********************************************************/
-
-/**********************************************************
-Generic domain class
-**********************************************************/
-struct domain;
-
-// vtable
-struct domain_vtable {
-  struct edge* (*new_edge) (struct domain *D);
-  struct cell* (*new_cell) (struct domain *D);
-  double (*compute_fluxes) (struct domain *D);
-  struct cell* (*get_cell) (struct domain *d, int i, struct cell *c); 
-  struct domain* (*destroy) (struct domain *D);
-};
-
-// structure
-struct domain {
-  struct domain_vtable *vtable;
-
-  struct quantities *vertex_values;
-  struct quantities *boundary_values;
-  struct quantities *explicit_update;
-  long* neighbours;
-  long* neighbour_vertices;
-  double* normals;
-  double* areas;
-  int number_of_elements;
-  int number_of_equations;
-  double *max_speed_array;
-  double timestep;
-  struct params *params;
-};
-
-// Public API
-// Destruction 
-struct domain* domain_destroy(struct domain* D);
-
-// method declarations
-struct edge* domain_new_edge(struct domain *qs);
-struct cell* domain_new_cell(struct domain *qs);
-double domain_compute_fluxes (struct domain *D);
-struct cell* domain_get_cell(struct domain *d, int i, struct cell *c); 
-
-// method implementations
-struct edge* domain_new_edge(struct domain *qs) {
-  return qs->vtable->new_edge(qs);
-}
-
-struct cell* domain_new_cell(struct domain *qs) {
-  return qs->vtable->new_cell(qs);
-}
-
-double domain_compute_fluxes (struct domain *D) {
-  return D->vtable->compute_fluxes(D);
-}
-
-struct cell* domain_get_cell(struct domain *D, int i, struct cell *c) {
-  return D->vtable->get_cell(D, i, c);
-}
-
-struct domain* domain_destroy(struct domain* D) {
-  return D->vtable->destroy(D);
-}
-/**********************************************************/
-
-
-/**********************************************************
-Generic methods
-**********************************************************/
-// Declarations
-double domain_compute_fluxes_generic(struct domain *D);
-struct cell* sqpipe_domain_get_cell(struct domain *D, int i, struct cell *c);
-double* edge_flux_function_godunov (struct edge *c, struct params *p, double *edgeflux, int num_eqns);
-double* cell_flux_function_generic (struct cell *c, struct params *p, double *cellflux, int num_eqns);
-
-// Implementations
-
-// Computational function for flux computation
-// This iterates over cells.
-double domain_compute_fluxes_generic (struct domain *D) {
-  double flux[D->number_of_equations];
-  double max_speed;
-  int k, i;
-  struct cell *c;
-   
-  // pre-create a cell object to use in the loop
-  // so we need only allocate memory once
-  c = domain_new_cell(D);
-
-  for (k=0; k<D->number_of_elements; k++) {
-    c = domain_get_cell(D, k, c);
-   
-    // Compute sound speeds (this updates c) and update timestep
-    max_speed = cell_extreme_sound_speeds(c, D->params);
-    if (max_speed > D->params->epsilon) {				    
-      D->timestep = min(D->timestep, 0.5*D->params->cfl * D->areas[k]/max_speed); 
-    }
-
-    // Compute flux
-    cell_flux_function(c, D->params, flux, D->number_of_equations);    
-    for (i=0; i<D->number_of_equations; i++) {
-      flux[i] /= D->areas[k];
-    }
-    quantities_update(D->explicit_update, flux, k);
-
-    //Keep track of maximal speeds
-    D->max_speed_array[k]=max_speed;
-  }
-  
-  cell_destroy(c);
-
-  return D->timestep;	
-}
-
-struct cell* domain_get_cell_generic(struct domain *D, int k, struct cell *c) {
-  int i, ki, n, m, nm;
-
-  for (i=0; i<c->num_edges; i++) {
-    ki = k*(c->num_edges) + i;
-    c->edges[i]->normal = D->normals[ki];      
-    n = D->neighbours[ki];
-
-    c->edges[i]->qin = quantities_get_quantity(D->vertex_values, ki, c->edges[i]->qin);
-
-    if (n<0) {
-      m = -n-1;
-      c->edges[i]->qout = quantities_get_quantity(D->boundary_values, m, c->edges[i]->qout);
-    } else {
-      m = D->neighbour_vertices[ki];
-      nm = n*2+m;
-      c->edges[i]->qout = quantities_get_quantity(D->vertex_values, nm, c->edges[i]->qout);
-    }
-  } 
-
-  return c;
-}
-
-
-// Generic Godunov flux computation method
-double* edge_flux_function_godunov (struct edge *e, struct params *p, double *edgeflux, int num_eqns) {
-  int i;
-  double flux_left[num_eqns], flux_right[num_eqns];
-  double denom;
-
-  // Flux contribution from left side of  edge
-  quantity_flux_formula(e->qin, e->normal, p, flux_left);
-
-  // Flux contribution from right side of edge
-  quantity_flux_formula(e->qout, e->normal, p, flux_right);
-
-  // Flux computation
-  denom = e->smax - e->smin;
-  if (denom < p->epsilon) {
-    for (i=0; i<num_eqns; i++) edgeflux[i] = 0.0;
-  } else {
-    for (i=0; i<num_eqns; i++) {
-      edgeflux[i] = e->smax*flux_left[i] - e->smin*flux_right[i];
-      edgeflux[i] += e->smax*e->smin * (quantity_get_conserved(e->qout, i, e->normal) - quantity_get_conserved(e->qin, i, e->normal));
-      edgeflux[i] /= denom;
-    }
-  }
-
-  return edgeflux;
-}
-
-double* cell_flux_function_generic (struct cell *c, struct params *p, double *cellflux, int num_eqns) {
-  int i, j;
-  double cellforce[num_eqns];
-  double flux[num_eqns];
-
-  // Compute forcing terms first
-  cell_forcing_terms(c, p, cellforce);
-  for (i=0; i<num_eqns; i++) {
-    cellflux[i] = -cellforce[i];
-  }
-
-  // Compute flux through each edge and append to flux
-  for (j=0; j<c->num_edges; j++) {
-    edge_flux_function(c->edges[j], p, flux, num_eqns);
-    for (i=0; i<num_eqns; i++) {
-      cellflux[i] -= flux[i];
-    }
-  }
-
-  return cellflux;
-}
-/**********************************************************/
Index: /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_comp_flux.c	(revision 8254)
+++ /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_comp_flux.c	(revision 8254)
@@ -0,0 +1,61 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+
+const double pi = 3.14159265358979;
+
+// Shared code snippets
+// Generic routines
+#include "util_ext.h"
+// Solver
+#include "solver.h"
+// Python-solver interface
+#include "python_solver.h"
+// square pipe solver implementation
+#include "sqpipe_default.h"
+// sqpipe-python glue
+#include "sqpipe_python.h"
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
+  struct domain *D = (struct domain *) sqpipe_default_domain_new();
+  double timestep;
+  PyArrayObject *state;
+  PyObject *domain;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dO", &timestep, &domain)) {
+    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
+    return NULL;
+  }
+
+  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
+
+  // We need the state 
+  state = get_consecutive_array(domain, "state");
+  ((struct sqpipe_default_domain *)D)->state = (long *) state->data;
+
+  timestep = domain_compute_fluxes(D);
+  domain_destroy(D);
+  
+  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 initsqpipe_comp_flux(void){
+  Py_InitModule("sqpipe_comp_flux", MethodTable);
+  import_array();
+}
Index: /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_default.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_default.h	(revision 8254)
+++ /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_default.h	(revision 8254)
@@ -0,0 +1,305 @@
+#include "sqpipe.h"
+
+/**********************************************************
+Square pipe default quantity class
+**********************************************************/
+
+struct sqpipe_default_quantity {
+  struct sqpipe_quantity super;
+  long state;
+};
+
+// Construction
+void sqpipe_default_quantity_init (struct sqpipe_default_quantity *q);
+struct sqpipe_default_quantity *sqpipe_default_quantity_new();
+
+// Destruction
+
+// Methods
+double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
+double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p);
+double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal);
+
+// Implementation
+double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
+  if (((struct sqpipe_default_quantity *)q)->state == 0) {
+    return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
+  } else {
+    return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
+  }
+}
+
+double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p) {
+  if (((struct sqpipe_default_quantity *)q)->state == 0) {
+    return sqpipe_quantity_free_surface_sound_speed(q, p);
+  } else {
+    return sqpipe_quantity_pressurised_sound_speed(q, p);
+  }
+}
+
+double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal) {
+  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
+  double c;
+  
+  switch (k) {
+  case 0:
+    c = p->a;
+    break;
+  case 1:
+    // This should be normal^2 p->d and normal is +/- 1
+    c = p->d;
+    break;
+  default:
+    c = 0;
+  }
+  return c;
+}
+
+void sqpipe_default_quantity_init (struct sqpipe_default_quantity *q) {
+  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
+
+  sqpipe_quantity_init(p);
+
+  static struct quantity_vtable vtable = {
+    &sqpipe_default_quantity_flux_formula,
+    &sqpipe_default_quantity_sound_speed,
+    &sqpipe_default_quantity_get_conserved,
+    &sqpipe_quantity_destroy
+  };
+  p->super.vtable = &vtable;
+
+  ((struct sqpipe_default_quantity *)p)->state = 0;
+}
+
+struct sqpipe_default_quantity *sqpipe_default_quantity_new() {
+  struct sqpipe_default_quantity *p = malloc(sizeof(struct sqpipe_default_quantity));
+
+  sqpipe_default_quantity_init(p);
+
+  return p;
+}
+/**********************************************************/
+
+
+/**********************************************************
+Square pipe default edge class
+**********************************************************/
+struct sqpipe_default_edge {
+  struct sqpipe_edge super;
+};
+
+// Construction
+void sqpipe_default_edge_init (struct sqpipe_default_edge *e);
+struct sqpipe_default_edge *sqpipe_default_edge_new();
+
+// Destruction
+
+// Methods
+
+// Implementation
+
+void sqpipe_default_edge_init (struct sqpipe_default_edge *q) {
+  sqpipe_edge_init((struct sqpipe_edge *)q);
+}
+
+struct sqpipe_default_edge *sqpipe_default_edge_new() {
+  struct sqpipe_default_edge *e = (struct sqpipe_default_edge *) sqpipe_edge_new();
+
+  sqpipe_default_edge_init(e);
+
+  // Create inner and outer quantities
+  ((struct edge *)e)->qin = (struct quantity *)sqpipe_default_quantity_new();
+  ((struct edge *)e)->qout = (struct quantity *)sqpipe_default_quantity_new();
+  
+  return e;
+}
+/**********************************************************/
+
+
+/**********************************************************
+Square pipe default cell class
+**********************************************************/
+struct sqpipe_default_cell {
+  struct sqpipe_cell super;
+
+  long state;
+};
+
+// Construction
+void sqpipe_default_cell_init (struct sqpipe_default_cell *e);
+struct sqpipe_default_cell *sqpipe_default_cell_new();
+
+// Destruction
+
+// Methods
+double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
+
+// Implementation
+double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
+  if (((struct sqpipe_default_cell *)c)->state == 0) {
+    return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);    
+  } else {
+    return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
+  }
+}
+
+void sqpipe_default_cell_init (struct sqpipe_default_cell *q) {
+  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
+
+  sqpipe_cell_init(p);
+
+  static struct cell_vtable vtable = {
+    &cell_flux_function_generic,
+    &sqpipe_cell_extreme_sound_speeds,
+    &sqpipe_default_cell_forcing_terms,
+    &sqpipe_cell_destroy
+  };
+  p->super.vtable = &vtable;
+
+  ((struct sqpipe_default_cell *)p)->state = 0;
+}
+
+struct sqpipe_default_cell *sqpipe_default_cell_new() {
+  int i;
+  struct sqpipe_default_cell *c = malloc(sizeof(struct sqpipe_default_cell));
+
+  // Allocate memory for a length 2 array of pointers to edges
+  ((struct cell *)c)->edges = (struct edge  **) malloc(2 * sizeof(struct sqpipe_edge *));
+
+  sqpipe_default_cell_init(c);
+
+  // Create edges
+  for (i=0; i<((struct cell *)c)->num_edges; i++) {
+    ((struct cell *)c)->edges[i] = (struct edge *) sqpipe_default_edge_new();
+  }
+
+  return c;
+}
+/**********************************************************/
+
+
+/**********************************************************
+Square pipe default quantities class
+**********************************************************/
+struct sqpipe_default_quantities {
+  struct sqpipe_quantities super;
+};
+
+// Construction
+void sqpipe_default_quantities_init (struct sqpipe_default_quantities *q);
+struct sqpipe_default_quantities *sqpipe_default_quantities_new();
+
+// Destruction
+
+// Methods
+void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k);
+
+
+// Implementation
+void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k) {
+  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
+  ps->a[k] = flux[0];
+  ps->d[k] = flux[1];
+}
+
+void sqpipe_default_quantities_init (struct sqpipe_default_quantities *q) {
+  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
+
+  sqpipe_quantities_init(p);
+
+  static struct quantities_vtable vtable = {
+    &sqpipe_quantities_get_quantity,
+    &sqpipe_default_quantities_update,
+    &sqpipe_quantities_destroy
+  };
+  p->super.vtable = &vtable;
+}
+
+struct sqpipe_default_quantities *sqpipe_default_quantities_new() {
+  struct sqpipe_default_quantities *p = (struct sqpipe_default_quantities *) sqpipe_quantities_new();
+
+  sqpipe_default_quantities_init(p);
+
+  return p;
+}
+/**********************************************************/
+
+
+/**********************************************************
+Square pipe default domain class
+**********************************************************/
+
+struct sqpipe_default_domain {
+  struct sqpipe_domain super;
+  long *state;
+};
+
+// Construction
+void sqpipe_default_domain_init (struct sqpipe_default_domain *d);
+struct sqpipe_default_domain *sqpipe_default_domain_new();
+
+// Destruction
+
+// Methods
+struct edge* sqpipe_default_domain_new_edge(struct domain *d);
+struct cell* sqpipe_default_domain_new_cell(struct domain *d);
+struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c);
+
+// Implementation
+struct edge* sqpipe_default_domain_new_edge(struct domain *d) {
+  struct edge *e = (struct edge*) sqpipe_default_edge_new();
+
+  return e;
+}
+
+struct cell* sqpipe_default_domain_new_cell(struct domain *d) {
+  struct cell *c = (struct cell*) sqpipe_default_cell_new();
+  return c;
+}
+
+struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c) {
+  struct sqpipe_default_cell *sc = (struct sqpipe_default_cell *) c;
+  struct sqpipe_default_domain *sD = (struct sqpipe_default_domain *) D;
+  int i;
+  
+  c = domain_get_cell_generic(D, k, c);
+
+  // Set state
+  sc->state = sD->state[k];
+  for (i = 0; i<c->num_edges; i++) {
+    ((struct sqpipe_default_quantity *)(c->edges[i]->qin))->state = sc->state;
+    ((struct sqpipe_default_quantity *)(c->edges[i]->qout))->state = sc->state;
+  }
+  
+  return c;
+}
+
+void sqpipe_default_domain_init (struct sqpipe_default_domain *d) {
+  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
+
+  sqpipe_domain_init(p);
+
+  static struct domain_vtable vtable = {
+    &sqpipe_default_domain_new_edge,
+    &sqpipe_default_domain_new_cell,
+    &domain_compute_fluxes_generic,
+    &sqpipe_default_domain_get_cell,
+    &sqpipe_domain_destroy
+  };
+  p->super.vtable = &vtable;
+
+  ((struct sqpipe_default_domain *)p)->state = NULL;
+}
+
+struct sqpipe_default_domain *sqpipe_default_domain_new() {
+  struct sqpipe_default_domain *d = (struct sqpipe_default_domain *) sqpipe_domain_new();
+
+  sqpipe_default_domain_init(d);
+
+  ((struct domain *) d)->vertex_values = (struct quantities *) sqpipe_default_quantities_new();
+  ((struct domain *) d)->boundary_values = (struct quantities *) sqpipe_default_quantities_new();
+  ((struct domain *) d)->explicit_update = (struct quantities *) sqpipe_default_quantities_new();
+
+  return d;
+}
+/**********************************************************/
Index: /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py	(revision 8253)
+++ /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py	(revision 8254)
@@ -46,10 +46,12 @@
 
 from anuga_1d.base.generic_domain import Generic_domain
+from anuga_1d.sqpipe.sqpipe_forcing_terms import *
+from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
 
 #Shallow water domain
-class Sqpipe_domain(Generic_domain):
-
-    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
-
+class Domain(Generic_domain):
+
+    def __init__(self, coordinates, boundary = None, forcing_terms = [], state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
+        conserved_quantities = ['area', 'discharge']
         evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
         other_quantities = ['friction']
@@ -149,5 +151,5 @@
 
         # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_comp_flux import compute_fluxes_ext
+        from anuga_1d.sqpipe.sqpipe_comp_flux import compute_fluxes_ext
  	self.flux_timestep = compute_fluxes_ext(timestep,self)        
 
@@ -171,6 +173,83 @@
 
     def distribute_to_vertices_and_edges(self):
-        # Should be implemented by the derived class
-        raise Exception("Not implemented")
+        # Shortcuts
+        h0 = self.h0
+        epsilon = self.epsilon
+
+        area      = self.quantities['area']
+        discharge = self.quantities['discharge']
+        bed       = self.quantities['elevation']
+        height    = self.quantities['height']
+        velocity  = self.quantities['velocity']
+        width     = self.quantities['width']
+        top       = self.quantities['top']
+        stage     = self.quantities['stage']
+
+        #Arrays   
+        a_C   = area.centroid_values
+        d_C   = discharge.centroid_values
+        z_C   = bed.centroid_values
+        h_C   = height.centroid_values
+        u_C   = velocity.centroid_values
+        b_C   = width.centroid_values
+        t_C   = top.centroid_values
+        w_C   = stage.centroid_values
+
+        #Calculate height (and fix negatives)better be non-negative!
+        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
+        w_C[:] = z_C + h_C
+        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
+
+        for name in ['velocity', 'stage' ]:
+            Q = self.quantities[name]
+            if self.order == 1:
+                Q.extrapolate_first_order()
+            elif self.order == 2:
+                Q.extrapolate_second_order()
+            else:
+                raise 'Unknown order'
+
+        # These have been extrapolated
+        w_V = self.quantities['stage'].vertex_values
+        u_V = self.quantities['velocity'].vertex_values
+
+        # These are the given geometry and remain fixed
+        z_V  = bed.vertex_values
+        b_V  = width.vertex_values
+        t_V  = top.vertex_values
+
+        # These need to be updated
+        a_V  = area.vertex_values
+        d_V  = discharge.vertex_values
+        h_V  = height.vertex_values
+
+        h_V[:,:]    = w_V - z_V
+
+        # Fix up small heights
+        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
+        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
+
+        h_V[:,0] = h_0
+        h_V[:,1] = h_1
+
+        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
+        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
+
+        h_V[:,0] = h_0
+        h_V[:,1] = h_1
+
+        # It may still be possible that h is small 
+        # If we set h to zero, we should also set u to 0
+        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
+        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
+
+        # Reconstruct conserved quantities and make everything consistent
+        # with new h values
+        w_V[:,:] = z_V + h_V
+        a_V[:,:] = h_V * b_V
+        d_V[:,:] = u_V * a_V
+
+
+        return
         
     def evolve(self, yieldstep = None, finaltime = None, duration = None,
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d.h	(revision 8253)
+++ 	(revision )
@@ -1,305 +1,0 @@
-#include "sqpipe.h"
-
-/**********************************************************
-Square pipe dual fixed width quantity class
-**********************************************************/
-
-struct sqpipe_dual_fixed_width_a_d_quantity {
-  struct sqpipe_quantity super;
-  long state;
-};
-
-// Construction
-void sqpipe_dual_fixed_width_a_d_quantity_init (struct sqpipe_dual_fixed_width_a_d_quantity *q);
-struct sqpipe_dual_fixed_width_a_d_quantity *sqpipe_dual_fixed_width_a_d_quantity_new();
-
-// Destruction
-
-// Methods
-double* sqpipe_dual_fixed_width_a_d_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
-double sqpipe_dual_fixed_width_a_d_quantity_sound_speed (struct quantity *q, struct params *p);
-double sqpipe_dual_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal);
-
-// Implementation
-double* sqpipe_dual_fixed_width_a_d_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
-  if (((struct sqpipe_dual_fixed_width_a_d_quantity *)q)->state == 0) {
-    return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
-  } else {
-    return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
-  }
-}
-
-double sqpipe_dual_fixed_width_a_d_quantity_sound_speed (struct quantity *q, struct params *p) {
-  if (((struct sqpipe_dual_fixed_width_a_d_quantity *)q)->state == 0) {
-    return sqpipe_quantity_free_surface_sound_speed(q, p);
-  } else {
-    return sqpipe_quantity_pressurised_sound_speed(q, p);
-  }
-}
-
-double sqpipe_dual_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = p->a;
-    break;
-  case 1:
-    // This should be normal^2 p->d and normal is +/- 1
-    c = p->d;
-    break;
-  default:
-    c = 0;
-  }
-  return c;
-}
-
-void sqpipe_dual_fixed_width_a_d_quantity_init (struct sqpipe_dual_fixed_width_a_d_quantity *q) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
-
-  sqpipe_quantity_init(p);
-
-  static struct quantity_vtable vtable = {
-    &sqpipe_dual_fixed_width_a_d_quantity_flux_formula,
-    &sqpipe_dual_fixed_width_a_d_quantity_sound_speed,
-    &sqpipe_dual_fixed_width_a_d_quantity_get_conserved,
-    &sqpipe_quantity_destroy
-  };
-  p->super.vtable = &vtable;
-
-  ((struct sqpipe_dual_fixed_width_a_d_quantity *)p)->state = 0;
-}
-
-struct sqpipe_dual_fixed_width_a_d_quantity *sqpipe_dual_fixed_width_a_d_quantity_new() {
-  struct sqpipe_dual_fixed_width_a_d_quantity *p = malloc(sizeof(struct sqpipe_dual_fixed_width_a_d_quantity));
-
-  sqpipe_dual_fixed_width_a_d_quantity_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe dual fixed width edge class
-**********************************************************/
-struct sqpipe_dual_fixed_width_a_d_edge {
-  struct sqpipe_edge super;
-};
-
-// Construction
-void sqpipe_dual_fixed_width_a_d_edge_init (struct sqpipe_dual_fixed_width_a_d_edge *e);
-struct sqpipe_dual_fixed_width_a_d_edge *sqpipe_dual_fixed_width_a_d_edge_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-
-void sqpipe_dual_fixed_width_a_d_edge_init (struct sqpipe_dual_fixed_width_a_d_edge *q) {
-  sqpipe_edge_init((struct sqpipe_edge *)q);
-}
-
-struct sqpipe_dual_fixed_width_a_d_edge *sqpipe_dual_fixed_width_a_d_edge_new() {
-  struct sqpipe_dual_fixed_width_a_d_edge *e = (struct sqpipe_dual_fixed_width_a_d_edge *) sqpipe_edge_new();
-
-  sqpipe_dual_fixed_width_a_d_edge_init(e);
-
-  // Create inner and outer quantities
-  ((struct edge *)e)->qin = (struct quantity *)sqpipe_dual_fixed_width_a_d_quantity_new();
-  ((struct edge *)e)->qout = (struct quantity *)sqpipe_dual_fixed_width_a_d_quantity_new();
-  
-  return e;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe dual fixed width cell class
-**********************************************************/
-struct sqpipe_dual_fixed_width_a_d_cell {
-  struct sqpipe_cell super;
-
-  long state;
-};
-
-// Construction
-void sqpipe_dual_fixed_width_a_d_cell_init (struct sqpipe_dual_fixed_width_a_d_cell *e);
-struct sqpipe_dual_fixed_width_a_d_cell *sqpipe_dual_fixed_width_a_d_cell_new();
-
-// Destruction
-
-// Methods
-double* sqpipe_dual_fixed_width_a_d_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
-
-// Implementation
-double* sqpipe_dual_fixed_width_a_d_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
-  if (((struct sqpipe_dual_fixed_width_a_d_cell *)c)->state == 0) {
-    return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);    
-  } else {
-    return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
-  }
-}
-
-void sqpipe_dual_fixed_width_a_d_cell_init (struct sqpipe_dual_fixed_width_a_d_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
-
-  sqpipe_cell_init(p);
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_generic,
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_dual_fixed_width_a_d_cell_forcing_terms,
-    &sqpipe_cell_destroy
-  };
-  p->super.vtable = &vtable;
-
-  ((struct sqpipe_dual_fixed_width_a_d_cell *)p)->state = 0;
-}
-
-struct sqpipe_dual_fixed_width_a_d_cell *sqpipe_dual_fixed_width_a_d_cell_new() {
-  int i;
-  struct sqpipe_dual_fixed_width_a_d_cell *c = malloc(sizeof(struct sqpipe_dual_fixed_width_a_d_cell));
-
-  // Allocate memory for a length 2 array of pointers to edges
-  ((struct cell *)c)->edges = (struct edge  **) malloc(2 * sizeof(struct sqpipe_edge *));
-
-  sqpipe_dual_fixed_width_a_d_cell_init(c);
-
-  // Create edges
-  for (i=0; i<((struct cell *)c)->num_edges; i++) {
-    ((struct cell *)c)->edges[i] = (struct edge *) sqpipe_dual_fixed_width_a_d_edge_new();
-  }
-
-  return c;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe dual fixed width quantities class
-**********************************************************/
-struct sqpipe_dual_fixed_width_a_d_quantities {
-  struct sqpipe_quantities super;
-};
-
-// Construction
-void sqpipe_dual_fixed_width_a_d_quantities_init (struct sqpipe_dual_fixed_width_a_d_quantities *q);
-struct sqpipe_dual_fixed_width_a_d_quantities *sqpipe_dual_fixed_width_a_d_quantities_new();
-
-// Destruction
-
-// Methods
-void sqpipe_dual_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k);
-
-
-// Implementation
-void sqpipe_dual_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k) {
-  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
-  ps->a[k] = flux[0];
-  ps->d[k] = flux[1];
-}
-
-void sqpipe_dual_fixed_width_a_d_quantities_init (struct sqpipe_dual_fixed_width_a_d_quantities *q) {
-  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
-
-  sqpipe_quantities_init(p);
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_get_quantity,
-    &sqpipe_dual_fixed_width_a_d_quantities_update,
-    &sqpipe_quantities_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_dual_fixed_width_a_d_quantities *sqpipe_dual_fixed_width_a_d_quantities_new() {
-  struct sqpipe_dual_fixed_width_a_d_quantities *p = (struct sqpipe_dual_fixed_width_a_d_quantities *) sqpipe_quantities_new();
-
-  sqpipe_dual_fixed_width_a_d_quantities_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width domain class
-**********************************************************/
-
-struct sqpipe_dual_fixed_width_a_d_domain {
-  struct sqpipe_domain super;
-  long *state;
-};
-
-// Construction
-void sqpipe_dual_fixed_width_a_d_domain_init (struct sqpipe_dual_fixed_width_a_d_domain *d);
-struct sqpipe_dual_fixed_width_a_d_domain *sqpipe_dual_fixed_width_a_d_domain_new();
-
-// Destruction
-
-// Methods
-struct edge* sqpipe_dual_fixed_width_a_d_domain_new_edge(struct domain *d);
-struct cell* sqpipe_dual_fixed_width_a_d_domain_new_cell(struct domain *d);
-struct cell* sqpipe_dual_fixed_width_a_d_domain_get_cell(struct domain *D, int k, struct cell *c);
-
-// Implementation
-struct edge* sqpipe_dual_fixed_width_a_d_domain_new_edge(struct domain *d) {
-  struct edge *e = (struct edge*) sqpipe_dual_fixed_width_a_d_edge_new();
-
-  return e;
-}
-
-struct cell* sqpipe_dual_fixed_width_a_d_domain_new_cell(struct domain *d) {
-  struct cell *c = (struct cell*) sqpipe_dual_fixed_width_a_d_cell_new();
-  return c;
-}
-
-struct cell* sqpipe_dual_fixed_width_a_d_domain_get_cell(struct domain *D, int k, struct cell *c) {
-  struct sqpipe_dual_fixed_width_a_d_cell *sc = (struct sqpipe_dual_fixed_width_a_d_cell *) c;
-  struct sqpipe_dual_fixed_width_a_d_domain *sD = (struct sqpipe_dual_fixed_width_a_d_domain *) D;
-  int i;
-  
-  c = domain_get_cell_generic(D, k, c);
-
-  // Set state
-  sc->state = sD->state[k];
-  for (i = 0; i<c->num_edges; i++) {
-    ((struct sqpipe_dual_fixed_width_a_d_quantity *)(c->edges[i]->qin))->state = sc->state;
-    ((struct sqpipe_dual_fixed_width_a_d_quantity *)(c->edges[i]->qout))->state = sc->state;
-  }
-  
-  return c;
-}
-
-void sqpipe_dual_fixed_width_a_d_domain_init (struct sqpipe_dual_fixed_width_a_d_domain *d) {
-  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
-
-  sqpipe_domain_init(p);
-
-  static struct domain_vtable vtable = {
-    &sqpipe_dual_fixed_width_a_d_domain_new_edge,
-    &sqpipe_dual_fixed_width_a_d_domain_new_cell,
-    &domain_compute_fluxes_generic,
-    &sqpipe_dual_fixed_width_a_d_domain_get_cell,
-    &sqpipe_domain_destroy
-  };
-  p->super.vtable = &vtable;
-
-  ((struct sqpipe_dual_fixed_width_a_d_domain *)p)->state = NULL;
-}
-
-struct sqpipe_dual_fixed_width_a_d_domain *sqpipe_dual_fixed_width_a_d_domain_new() {
-  struct sqpipe_dual_fixed_width_a_d_domain *d = (struct sqpipe_dual_fixed_width_a_d_domain *) sqpipe_domain_new();
-
-  sqpipe_dual_fixed_width_a_d_domain_init(d);
-
-  ((struct domain *) d)->vertex_values = (struct quantities *) sqpipe_dual_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->boundary_values = (struct quantities *) sqpipe_dual_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->explicit_update = (struct quantities *) sqpipe_dual_fixed_width_a_d_quantities_new();
-
-  return d;
-}
-/**********************************************************/
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d_comp_flux.c	(revision 8253)
+++ 	(revision )
@@ -1,56 +1,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-
-const double pi = 3.14159265358979;
-
-// Shared code snippets
-#include "util_ext.h"
-#include "solver.h"
-#include "python_solver.h"
-#include "sqpipe_dual_fixed_width_a_d.h"
-#include "sqpipe_python.h"
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  struct domain *D = (struct domain *) sqpipe_dual_fixed_width_a_d_domain_new();
-  double timestep;
-  PyArrayObject *state;
-  PyObject *domain;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dO", &timestep, &domain)) {
-    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
-    return NULL;
-  }
-
-  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
-
-  // We need the state for dual models
-  state = get_consecutive_array(domain, "state");
-  ((struct sqpipe_dual_fixed_width_a_d_domain *)D)->state = (long *) state->data;
-
-  timestep = domain_compute_fluxes(D);
-  domain_destroy(D);
-  
-  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 initsqpipe_dual_fixed_width_a_d_comp_flux(void){
-  Py_InitModule("sqpipe_dual_fixed_width_a_d_comp_flux", MethodTable);
-  import_array();
-}
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_dual_fixed_width_a_d_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,129 +1,0 @@
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None, state = None, update_state = True):
-        conserved_quantities = ['area', 'discharge']
-        forcing_terms = []
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               tagged_elements      = tagged_elements)
-
-        # Allocate space for tracking state if not provided
-        if (state is None):
-            self.state = numpy.zeros(self.number_of_elements, numpy.int)
-        else:
-            self.state = state
-            if (len(self.state) != self.number_of_elements):
-                raise Exception("state has invalid length")
-
-        self.update_state = update_state
-
-	self.__doc__ = 'sqpipe_area_discharge_domain'
-
-    def compute_fluxes(self):        
-        #Call correct module function
-        #(either from this module or C-extension)
-
-        # Set initial timestep to a large value
-        import sys
-        timestep = float(sys.maxint)
-
-        # Update state before computing flux
-        if self.update_state :
-            h = self.quantities['height'].centroid_values
-            t = self.quantities['top'].centroid_values
-            self.state = numpy.where(h>t, 1, 0)
-
-        # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep,self)        
-        
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-
-        #Calculate height (and fix negatives)better be non-negative!
-        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
-        w_C[:] = z_C + h_C
-        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
-
-        for name in ['velocity', 'stage' ]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        w_V = self.quantities['stage'].vertex_values
-        u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        a_V  = area.vertex_values
-        d_V  = discharge.vertex_values
-        h_V  = height.vertex_values
-
-        h_V[:,:]    = w_V - z_V
-
-        # Fix up small heights
-        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
-        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
-        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        # It may still be possible that h is small 
-        # If we set h to zero, we should also set u to 0
-        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
-        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
-
-        # Reconstruct conserved quantities and make everything consistent
-        # with new h values
-        w_V[:,:] = z_V + h_V
-        a_V[:,:] = h_V * b_V
-        d_V[:,:] = u_V * a_V
-
-
-        return
-
-    
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d.h	(revision 8253)
+++ 	(revision )
@@ -1,245 +1,0 @@
-#include "sqpipe.h"
-
-/**********************************************************
-Square pipe free surface fixed width quantity class
-**********************************************************/
-
-struct sqpipe_free_surface_fixed_width_a_d_quantity {
-  struct sqpipe_quantity super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_a_d_quantity_init (struct sqpipe_free_surface_fixed_width_a_d_quantity *q);
-struct sqpipe_free_surface_fixed_width_a_d_quantity *sqpipe_free_surface_fixed_width_a_d_quantity_new();
-
-// Destruction
-
-// Methods
-double sqpipe_free_surface_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal);
-
-// Implementation
-double sqpipe_free_surface_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = p->a;
-    break;
-  case 1:
-    // This should be normal^2 p->d and normal is +/- 1
-    c = p->d;
-    break;
-  default:
-    c = 0;
-  }
-  return c;
-}
-
-void sqpipe_free_surface_fixed_width_a_d_quantity_init (struct sqpipe_free_surface_fixed_width_a_d_quantity *q) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
-
-  sqpipe_quantity_init(p);
-
-  static struct quantity_vtable vtable = {
-    &sqpipe_quantity_free_surface_flux_formula,
-    &sqpipe_quantity_free_surface_sound_speed,
-    &sqpipe_free_surface_fixed_width_a_d_quantity_get_conserved,
-    &sqpipe_quantity_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_a_d_quantity *sqpipe_free_surface_fixed_width_a_d_quantity_new() {
-  struct sqpipe_free_surface_fixed_width_a_d_quantity *p = (struct sqpipe_free_surface_fixed_width_a_d_quantity *) sqpipe_quantity_new();
-
-  sqpipe_free_surface_fixed_width_a_d_quantity_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width edge class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_a_d_edge {
-  struct sqpipe_edge super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_a_d_edge_init (struct sqpipe_free_surface_fixed_width_a_d_edge *e);
-struct sqpipe_free_surface_fixed_width_a_d_edge *sqpipe_free_surface_fixed_width_a_d_edge_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-
-void sqpipe_free_surface_fixed_width_a_d_edge_init (struct sqpipe_free_surface_fixed_width_a_d_edge *q) {
-  sqpipe_edge_init((struct sqpipe_edge *)q);
-}
-
-struct sqpipe_free_surface_fixed_width_a_d_edge *sqpipe_free_surface_fixed_width_a_d_edge_new() {
-  struct sqpipe_free_surface_fixed_width_a_d_edge *e = (struct sqpipe_free_surface_fixed_width_a_d_edge *) sqpipe_edge_new();
-
-  sqpipe_free_surface_fixed_width_a_d_edge_init(e);
-
-  // Create inner and outer quantities
-  ((struct edge *)e)->qin = (struct quantity *)sqpipe_free_surface_fixed_width_a_d_quantity_new();
-  ((struct edge *)e)->qout = (struct quantity *)sqpipe_free_surface_fixed_width_a_d_quantity_new();
-
-  return e;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width cell class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_a_d_cell {
-  struct sqpipe_cell super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_a_d_cell_init (struct sqpipe_free_surface_fixed_width_a_d_cell *e);
-struct sqpipe_free_surface_fixed_width_a_d_cell *sqpipe_free_surface_fixed_width_a_d_cell_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-void sqpipe_free_surface_fixed_width_a_d_cell_init (struct sqpipe_free_surface_fixed_width_a_d_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
-
-  sqpipe_cell_init(p);
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_generic, 
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_cell_free_surface_forcing_terms,
-    &sqpipe_cell_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_a_d_cell *sqpipe_free_surface_fixed_width_a_d_cell_new() {
-  int i;  
-  struct sqpipe_free_surface_fixed_width_a_d_cell *c = (struct sqpipe_free_surface_fixed_width_a_d_cell *) sqpipe_cell_new();
-
-  sqpipe_free_surface_fixed_width_a_d_cell_init(c);
-
-  // Create edges
-  for (i=0; i<((struct cell *)c)->num_edges; i++) {
-    ((struct cell *)c)->edges[i] = (struct edge *) sqpipe_free_surface_fixed_width_a_d_edge_new();
-  }
-
-  return c;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width quantities class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_a_d_quantities {
-  struct sqpipe_quantities super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_a_d_quantities_init (struct sqpipe_free_surface_fixed_width_a_d_quantities *q);
-struct sqpipe_free_surface_fixed_width_a_d_quantities *sqpipe_free_surface_fixed_width_a_d_quantities_new();
-
-// Destruction
-
-// Methods
-void sqpipe_free_surface_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k);
-
-// Implementation
-void sqpipe_free_surface_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k) {
-  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
-  ps->a[k] = flux[0];
-  ps->d[k] = flux[1];
-}
-
-void sqpipe_free_surface_fixed_width_a_d_quantities_init (struct sqpipe_free_surface_fixed_width_a_d_quantities *q) {
-  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
-
-  sqpipe_quantities_init(p);
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_get_quantity,
-    &sqpipe_free_surface_fixed_width_a_d_quantities_update,
-    &sqpipe_quantities_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_a_d_quantities *sqpipe_free_surface_fixed_width_a_d_quantities_new() {
-  struct sqpipe_free_surface_fixed_width_a_d_quantities *p = (struct sqpipe_free_surface_fixed_width_a_d_quantities *) sqpipe_quantities_new();
-
-  sqpipe_free_surface_fixed_width_a_d_quantities_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-/**********************************************************
-Square pipe free surface fixed width domain class
-**********************************************************/
-
-struct sqpipe_free_surface_fixed_width_a_d_domain {
-  struct sqpipe_domain super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_a_d_domain_init (struct sqpipe_free_surface_fixed_width_a_d_domain *d);
-struct sqpipe_free_surface_fixed_width_a_d_domain *sqpipe_free_surface_fixed_width_a_d_domain_new();
-
-// Destruction
-
-// Methods
-struct edge* sqpipe_free_surface_fixed_width_a_d_domain_new_edge(struct domain *d);
-struct cell* sqpipe_free_surface_fixed_width_a_d_domain_new_cell(struct domain *d);
-
-// Implementation
-struct edge* sqpipe_free_surface_fixed_width_a_d_domain_new_edge(struct domain *d) {
-  struct edge *e = (struct edge*) sqpipe_free_surface_fixed_width_a_d_edge_new();
-
-  return e;
-}
-
-struct cell* sqpipe_free_surface_fixed_width_a_d_domain_new_cell(struct domain *d) {
-  struct cell *c = (struct cell*) sqpipe_free_surface_fixed_width_a_d_cell_new();
-  return c;
-}
-
-void sqpipe_free_surface_fixed_width_a_d_domain_init (struct sqpipe_free_surface_fixed_width_a_d_domain *d) {
-  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
-
-  sqpipe_domain_init(p);
-
-  static struct domain_vtable vtable = {
-    &sqpipe_free_surface_fixed_width_a_d_domain_new_edge,
-    &sqpipe_free_surface_fixed_width_a_d_domain_new_cell,
-    &domain_compute_fluxes_generic,
-    &domain_get_cell_generic,
-    &sqpipe_domain_destroy
-  };
-  p->super.vtable = &vtable;  
-}
-
-struct sqpipe_free_surface_fixed_width_a_d_domain *sqpipe_free_surface_fixed_width_a_d_domain_new() {
-  struct sqpipe_free_surface_fixed_width_a_d_domain *d = (struct sqpipe_free_surface_fixed_width_a_d_domain *) sqpipe_domain_new();
-
-  sqpipe_free_surface_fixed_width_a_d_domain_init(d);
-
-  ((struct domain *) d)->vertex_values = (struct quantities *) sqpipe_free_surface_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->boundary_values = (struct quantities *) sqpipe_free_surface_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->explicit_update = (struct quantities *) sqpipe_free_surface_fixed_width_a_d_quantities_new();
-
-  return d;
-}
-/**********************************************************/
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d_comp_flux.c	(revision 8253)
+++ 	(revision )
@@ -1,51 +1,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-
-const double pi = 3.14159265358979;
-
-// Shared code snippets
-#include "util_ext.h"
-#include "solver.h"
-#include "python_solver.h"
-#include "sqpipe_free_surface_fixed_width_a_d.h"
-#include "sqpipe_python.h"
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  struct domain *D = (struct domain *) sqpipe_free_surface_fixed_width_a_d_domain_new();
-  double timestep;
-  PyObject *domain;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dO", &timestep, &domain)) {
-    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
-    return NULL;
-  }
-
-  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
-
-  timestep = domain_compute_fluxes(D);
-  domain_destroy(D);
-  
-  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 initsqpipe_free_surface_fixed_width_a_d_comp_flux(void){
-  Py_InitModule("sqpipe_free_surface_fixed_width_a_d_comp_flux", MethodTable);
-  import_array();
-}
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_a_d_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,112 +1,0 @@
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None):
-        conserved_quantities = ['area', 'discharge']
-        forcing_terms = []
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               tagged_elements      = tagged_elements)
-
-	self.__doc__ = 'sqpipe_area_discharge_domain'
-
-    def compute_fluxes(self):        
-        #Call correct module function
-        #(either from this module or C-extension)
-
-        # Set initial timestep to a large value
-        import sys
-        timestep = float(sys.maxint)
-
-        # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_free_surface_fixed_width_a_d_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep,self)        
-        
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-
-        #Calculate height (and fix negatives)better be non-negative!
-        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
-        w_C[:] = z_C + h_C
-        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
-
-        for name in ['velocity', 'stage' ]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        w_V = self.quantities['stage'].vertex_values
-        u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        a_V  = area.vertex_values
-        d_V  = discharge.vertex_values
-        h_V  = height.vertex_values
-
-        h_V[:,:]    = w_V - z_V
-
-        # Fix up small heights
-        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
-        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
-        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        # It may still be possible that h is small 
-        # If we set h to zero, we should also set u to 0
-        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
-        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
-
-        # Reconstruct conserved quantities and make everything consistent
-        # with new h and u values
-        w_V[:,:] = z_V + h_V
-        a_V[:,:] = h_V * b_V
-        d_V[:,:] = u_V * a_V
-
-        return
-
-    
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d.h	(revision 8253)
+++ 	(revision )
@@ -1,259 +1,0 @@
-#include "sqpipe.h"
-
-/**********************************************************
-Square pipe free surface fixed width quantity class
-**********************************************************/
-
-struct sqpipe_free_surface_fixed_width_w_d_quantity {
-  struct sqpipe_quantity super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_w_d_quantity_init (struct sqpipe_free_surface_fixed_width_w_d_quantity *q);
-struct sqpipe_free_surface_fixed_width_w_d_quantity *sqpipe_free_surface_fixed_width_w_d_quantity_new();
-
-// Destruction
-
-// Methods
-double* sqpipe_free_surface_fixed_width_w_d_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
-double sqpipe_free_surface_fixed_width_w_d_quantity_get_conserved (struct quantity *q, int k, double normal);
-
-// Implementation
-// Stage is the conserved quantity, but we need to use h in the flux
-// formula, hence we can't use the sqpipe_quantity_free_surface_flux_formula
-// so implement our own method here.
-double* sqpipe_free_surface_fixed_width_w_d_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
-  struct sqpipe_quantity *pq = (struct sqpipe_quantity*) q;
-  struct sqpipe_params *sp = (struct sqpipe_params *) p;  
-
-  quantityflux[0] = (normal * pq->u) * pq->h;
-  quantityflux[1] = normal * (pq->u * pq->d + 0.5 * sp->g * pq->h *pq->h);
-
-  return quantityflux;
-}
-
-double sqpipe_free_surface_fixed_width_w_d_quantity_get_conserved (struct quantity *q, int k, double normal) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = p->w;
-    break;
-  case 1:
-    // This should be normal^2 p->d and normal is +/- 1
-    c = p->d;
-    break;
-  default:
-    c = 0;
-  }
-  return c;
-}
-
-void sqpipe_free_surface_fixed_width_w_d_quantity_init (struct sqpipe_free_surface_fixed_width_w_d_quantity *q) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
-
-  sqpipe_quantity_init(p);
-
-  static struct quantity_vtable vtable = {
-    &sqpipe_free_surface_fixed_width_w_d_quantity_flux_formula,
-    &sqpipe_quantity_free_surface_sound_speed,
-    &sqpipe_free_surface_fixed_width_w_d_quantity_get_conserved,
-    &sqpipe_quantity_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_w_d_quantity *sqpipe_free_surface_fixed_width_w_d_quantity_new() {
-  struct sqpipe_free_surface_fixed_width_w_d_quantity *p = (struct sqpipe_free_surface_fixed_width_w_d_quantity *) sqpipe_quantity_new();
-
-  sqpipe_free_surface_fixed_width_w_d_quantity_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width edge class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_w_d_edge {
-  struct sqpipe_edge super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_w_d_edge_init (struct sqpipe_free_surface_fixed_width_w_d_edge *e);
-struct sqpipe_free_surface_fixed_width_w_d_edge *sqpipe_free_surface_fixed_width_w_d_edge_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-
-void sqpipe_free_surface_fixed_width_w_d_edge_init (struct sqpipe_free_surface_fixed_width_w_d_edge *q) {
-  sqpipe_edge_init((struct sqpipe_edge *)q);
-}
-
-struct sqpipe_free_surface_fixed_width_w_d_edge *sqpipe_free_surface_fixed_width_w_d_edge_new() {
-  struct sqpipe_free_surface_fixed_width_w_d_edge *e = (struct sqpipe_free_surface_fixed_width_w_d_edge *) sqpipe_edge_new();
-
-  sqpipe_free_surface_fixed_width_w_d_edge_init(e);
-
-  // Create inner and outer quantities
-  ((struct edge *)e)->qin = (struct quantity *)sqpipe_free_surface_fixed_width_w_d_quantity_new();
-  ((struct edge *)e)->qout = (struct quantity *)sqpipe_free_surface_fixed_width_w_d_quantity_new();
-
-  return e;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width cell class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_w_d_cell {
-  struct sqpipe_cell super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_w_d_cell_init (struct sqpipe_free_surface_fixed_width_w_d_cell *e);
-struct sqpipe_free_surface_fixed_width_w_d_cell *sqpipe_free_surface_fixed_width_w_d_cell_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-void sqpipe_free_surface_fixed_width_w_d_cell_init (struct sqpipe_free_surface_fixed_width_w_d_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
-
-  sqpipe_cell_init(p);
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_generic,
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_cell_free_surface_forcing_terms,
-    &sqpipe_cell_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_w_d_cell *sqpipe_free_surface_fixed_width_w_d_cell_new() {
-  int i;
-  struct sqpipe_free_surface_fixed_width_w_d_cell *c = (struct sqpipe_free_surface_fixed_width_w_d_cell *) sqpipe_cell_new();
-
-  sqpipe_free_surface_fixed_width_w_d_cell_init(c);
-
-  // Create edges
-  for (i=0; i<((struct cell *)c)->num_edges; i++) {
-    ((struct cell *)c)->edges[i] = (struct edge *) sqpipe_free_surface_fixed_width_w_d_edge_new();
-  }
-
-  return c;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width quantities class
-**********************************************************/
-struct sqpipe_free_surface_fixed_width_w_d_quantities {
-  struct sqpipe_quantities super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_w_d_quantities_init (struct sqpipe_free_surface_fixed_width_w_d_quantities *q);
-struct sqpipe_free_surface_fixed_width_w_d_quantities *sqpipe_free_surface_fixed_width_w_d_quantities_new();
-
-// Destruction
-
-// Methods
-void sqpipe_free_surface_fixed_width_w_d_quantities_update (struct quantities *qs, double *flux, int k);
-
-// Implementation
-void sqpipe_free_surface_fixed_width_w_d_quantities_update (struct quantities *qs, double *flux, int k) {
-  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
-  ps->w[k] = flux[0];
-  ps->d[k] = flux[1];
-}
-
-void sqpipe_free_surface_fixed_width_w_d_quantities_init (struct sqpipe_free_surface_fixed_width_w_d_quantities *q) {
-  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
-
-  sqpipe_quantities_init(p);
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_get_quantity,
-    &sqpipe_free_surface_fixed_width_w_d_quantities_update,
-    &sqpipe_quantities_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_w_d_quantities *sqpipe_free_surface_fixed_width_w_d_quantities_new() {
-  struct sqpipe_free_surface_fixed_width_w_d_quantities *p = (struct sqpipe_free_surface_fixed_width_w_d_quantities *) sqpipe_quantities_new();
-
-  sqpipe_free_surface_fixed_width_w_d_quantities_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-/**********************************************************
-Square pipe free surface fixed width domain class
-**********************************************************/
-
-struct sqpipe_free_surface_fixed_width_w_d_domain {
-  struct sqpipe_domain super;
-};
-
-// Construction
-void sqpipe_free_surface_fixed_width_w_d_domain_init (struct sqpipe_free_surface_fixed_width_w_d_domain *d);
-struct sqpipe_free_surface_fixed_width_w_d_domain *sqpipe_free_surface_fixed_width_w_d_domain_new();
-
-// Destruction
-
-// Methods
-struct edge* sqpipe_free_surface_fixed_width_w_d_domain_new_edge(struct domain *d);
-struct cell* sqpipe_free_surface_fixed_width_w_d_domain_new_cell(struct domain *d);
-
-// Implementation
-struct edge* sqpipe_free_surface_fixed_width_w_d_domain_new_edge(struct domain *d) {
-  struct edge *e = (struct edge*) sqpipe_free_surface_fixed_width_w_d_edge_new();
-
-  return e;
-}
-
-struct cell* sqpipe_free_surface_fixed_width_w_d_domain_new_cell(struct domain *d) {
-  struct cell *c = (struct cell*) sqpipe_free_surface_fixed_width_w_d_cell_new();
-  return c;
-}
-
-void sqpipe_free_surface_fixed_width_w_d_domain_init (struct sqpipe_free_surface_fixed_width_w_d_domain *d) {
-  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
-
-  sqpipe_domain_init(p);
-
-  static struct domain_vtable vtable = {
-    &sqpipe_free_surface_fixed_width_w_d_domain_new_edge,
-    &sqpipe_free_surface_fixed_width_w_d_domain_new_cell,
-    &domain_compute_fluxes_generic,
-    &domain_get_cell_generic,
-    &sqpipe_domain_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_free_surface_fixed_width_w_d_domain *sqpipe_free_surface_fixed_width_w_d_domain_new() {
-  struct sqpipe_free_surface_fixed_width_w_d_domain *d = (struct sqpipe_free_surface_fixed_width_w_d_domain *) sqpipe_domain_new();
-
-  sqpipe_free_surface_fixed_width_w_d_domain_init(d);
-
-  ((struct domain *) d)->vertex_values = (struct quantities *) sqpipe_free_surface_fixed_width_w_d_quantities_new();
-  ((struct domain *) d)->boundary_values = (struct quantities *) sqpipe_free_surface_fixed_width_w_d_quantities_new();
-  ((struct domain *) d)->explicit_update = (struct quantities *) sqpipe_free_surface_fixed_width_w_d_quantities_new();
-
-  return d;
-}
-/**********************************************************/
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d_comp_flux.c	(revision 8253)
+++ 	(revision )
@@ -1,51 +1,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-
-const double pi = 3.14159265358979;
-
-// Shared code snippets
-#include "util_ext.h"
-#include "solver.h"
-#include "python_solver.h"
-#include "sqpipe_free_surface_fixed_width_w_d.h"
-#include "sqpipe_python.h"
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  struct domain *D = (struct domain *) sqpipe_free_surface_fixed_width_w_d_domain_new();
-  double timestep;
-  PyObject *domain;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dO", &timestep, &domain)) {
-    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
-    return NULL;
-  }
-
-  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
-
-  timestep = domain_compute_fluxes(D);
-  domain_destroy(D);
-  
-  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 initsqpipe_free_surface_fixed_width_w_d_comp_flux(void){
-  Py_InitModule("sqpipe_free_surface_fixed_width_w_d_comp_flux", MethodTable);
-  import_array();
-}
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_free_surface_fixed_width_w_d_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,112 +1,0 @@
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None):
-        conserved_quantities = ['stage', 'discharge']
-        forcing_terms = []
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               tagged_elements      = tagged_elements)
-
-	self.__doc__ = 'sqpipe_stage_discharge_domain'
-
-    def compute_fluxes(self):        
-        #Call correct module function
-        #(either from this module or C-extension)
-
-        # Set initial timestep to a large value
-        import sys
-        timestep = float(sys.maxint)
-
-        # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_free_surface_fixed_width_w_d_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep, self)        
-        
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-
-        #Calculate height (and fix negatives)better be non-negative!        
-        w_C[:] = numpy.maximum(w_C, z_C)
-        h_C[:] = w_C - z_C
-        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
-
-        for name in ['velocity', 'stage' ]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        w_V = self.quantities['stage'].vertex_values
-        u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        a_V  = area.vertex_values
-        d_V  = discharge.vertex_values
-        h_V  = height.vertex_values
-
-        h_V[:,:]    = w_V - z_V
-
-        # Fix up small heights
-        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
-        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
-        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        # It may still be possible that h is small 
-        # If we set h to zero, we should also set u to 0
-        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
-        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
-
-        # Reconstruct conserved quantities and make everything consistent
-        # with new h and u values
-        w_V[:,:] = z_V + h_V
-        a_V[:,:] = h_V * b_V
-        d_V[:,:] = u_V * a_V
-        
-        return
-
-    
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d.h	(revision 8253)
+++ 	(revision )
@@ -1,247 +1,0 @@
-#include "sqpipe.h"
-
-/**********************************************************
-Square pipe pressurised fixed width quantity class
-**********************************************************/
-
-struct sqpipe_pressurised_fixed_width_a_d_quantity {
-  struct sqpipe_quantity super;
-};
-
-// Construction
-void sqpipe_pressurised_fixed_width_a_d_quantity_init (struct sqpipe_pressurised_fixed_width_a_d_quantity *q);
-struct sqpipe_pressurised_fixed_width_a_d_quantity *sqpipe_pressurised_fixed_width_a_d_quantity_new();
-
-// Destruction
-
-// Methods
-double sqpipe_pressurised_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal);
-
-// Implementation
-double sqpipe_pressurised_fixed_width_a_d_quantity_get_conserved (struct quantity *q, int k, double normal) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = p->a;
-    break;
-  case 1:
-    // This should be normal^2 p->d and normal is +/- 1
-    c = p->d;
-    break;
-  default:
-    c = 0;
-  }
-  return c;
-}
-
-void sqpipe_pressurised_fixed_width_a_d_quantity_init (struct sqpipe_pressurised_fixed_width_a_d_quantity *q) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
-
-  sqpipe_quantity_init(p);
-
-  static struct quantity_vtable vtable = {
-    &sqpipe_quantity_pressurised_flux_formula,
-    &sqpipe_quantity_pressurised_sound_speed,
-    &sqpipe_pressurised_fixed_width_a_d_quantity_get_conserved,
-    &sqpipe_quantity_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_pressurised_fixed_width_a_d_quantity *sqpipe_pressurised_fixed_width_a_d_quantity_new() {
-  struct sqpipe_pressurised_fixed_width_a_d_quantity *p = (struct sqpipe_pressurised_fixed_width_a_d_quantity *) sqpipe_quantity_new();
-
-  sqpipe_pressurised_fixed_width_a_d_quantity_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe pressurised fixed width edge class
-**********************************************************/
-struct sqpipe_pressurised_fixed_width_a_d_edge {
-  struct sqpipe_edge super;
-};
-
-// Construction
-void sqpipe_pressurised_fixed_width_a_d_edge_init (struct sqpipe_pressurised_fixed_width_a_d_edge *e);
-struct sqpipe_pressurised_fixed_width_a_d_edge *sqpipe_pressurised_fixed_width_a_d_edge_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-
-void sqpipe_pressurised_fixed_width_a_d_edge_init (struct sqpipe_pressurised_fixed_width_a_d_edge *q) {
-  sqpipe_edge_init((struct sqpipe_edge *)q);
-}
-
-struct sqpipe_pressurised_fixed_width_a_d_edge *sqpipe_pressurised_fixed_width_a_d_edge_new() {
-  struct sqpipe_pressurised_fixed_width_a_d_edge *e = (struct sqpipe_pressurised_fixed_width_a_d_edge *) sqpipe_edge_new();
-
-  sqpipe_pressurised_fixed_width_a_d_edge_init(e);
-
-  // Create inner and outer quantities
-  ((struct edge *)e)->qin = (struct quantity *)sqpipe_pressurised_fixed_width_a_d_quantity_new();
-  ((struct edge *)e)->qout = (struct quantity *)sqpipe_pressurised_fixed_width_a_d_quantity_new();
-  
-  return e;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe pressurised fixed width cell class
-**********************************************************/
-struct sqpipe_pressurised_fixed_width_a_d_cell {
-  struct sqpipe_cell super;
-};
-
-// Construction
-void sqpipe_pressurised_fixed_width_a_d_cell_init (struct sqpipe_pressurised_fixed_width_a_d_cell *e);
-struct sqpipe_pressurised_fixed_width_a_d_cell *sqpipe_pressurised_fixed_width_a_d_cell_new();
-
-// Destruction
-
-// Methods
-
-// Implementation
-void sqpipe_pressurised_fixed_width_a_d_cell_init (struct sqpipe_pressurised_fixed_width_a_d_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
-
-  sqpipe_cell_init(p);
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_generic,
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_cell_pressurised_forcing_terms,
-    &sqpipe_cell_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_pressurised_fixed_width_a_d_cell *sqpipe_pressurised_fixed_width_a_d_cell_new() {
-  int i;
-  struct sqpipe_pressurised_fixed_width_a_d_cell *c = (struct sqpipe_pressurised_fixed_width_a_d_cell *) sqpipe_cell_new();
-
-  sqpipe_pressurised_fixed_width_a_d_cell_init(c);
-
-  // Create edges
-  for (i=0; i<((struct cell *)c)->num_edges; i++) {
-    ((struct cell *)c)->edges[i] = (struct edge *) sqpipe_pressurised_fixed_width_a_d_edge_new();
-  }
-
-  return c;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe pressurised fixed width quantities class
-**********************************************************/
-struct sqpipe_pressurised_fixed_width_a_d_quantities {
-  struct sqpipe_quantities super;
-};
-
-// Construction
-void sqpipe_pressurised_fixed_width_a_d_quantities_init (struct sqpipe_pressurised_fixed_width_a_d_quantities *q);
-struct sqpipe_pressurised_fixed_width_a_d_quantities *sqpipe_pressurised_fixed_width_a_d_quantities_new();
-
-// Destruction
-
-// Methods
-void sqpipe_pressurised_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k);
-
-
-// Implementation
-void sqpipe_pressurised_fixed_width_a_d_quantities_update (struct quantities *qs, double *flux, int k) {
-  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
-  ps->a[k] = flux[0];
-  ps->d[k] = flux[1];
-}
-
-void sqpipe_pressurised_fixed_width_a_d_quantities_init (struct sqpipe_pressurised_fixed_width_a_d_quantities *q) {
-  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
-
-  sqpipe_quantities_init(p);
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_get_quantity,
-    &sqpipe_pressurised_fixed_width_a_d_quantities_update,
-    &sqpipe_quantities_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_pressurised_fixed_width_a_d_quantities *sqpipe_pressurised_fixed_width_a_d_quantities_new() {
-  struct sqpipe_pressurised_fixed_width_a_d_quantities *p = (struct sqpipe_pressurised_fixed_width_a_d_quantities *) sqpipe_quantities_new();
-
-  sqpipe_pressurised_fixed_width_a_d_quantities_init(p);
-
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square pipe free surface fixed width domain class
-**********************************************************/
-
-struct sqpipe_pressurised_fixed_width_a_d_domain {
-  struct sqpipe_domain super;
-};
-
-// Construction
-void sqpipe_pressurised_fixed_width_a_d_domain_init (struct sqpipe_pressurised_fixed_width_a_d_domain *d);
-struct sqpipe_pressurised_fixed_width_a_d_domain *sqpipe_pressurised_fixed_width_a_d_domain_new();
-
-// Destruction
-
-// Methods
-struct edge* sqpipe_pressurised_fixed_width_a_d_domain_new_edge(struct domain *d);
-struct cell* sqpipe_pressurised_fixed_width_a_d_domain_new_cell(struct domain *d);
-
-// Implementation
-struct edge* sqpipe_pressurised_fixed_width_a_d_domain_new_edge(struct domain *d) {
-  struct edge *e = (struct edge*) sqpipe_pressurised_fixed_width_a_d_edge_new();
-
-  return e;
-}
-
-struct cell* sqpipe_pressurised_fixed_width_a_d_domain_new_cell(struct domain *d) {
-  struct cell *c = (struct cell*) sqpipe_pressurised_fixed_width_a_d_cell_new();
-  return c;
-}
-
-void sqpipe_pressurised_fixed_width_a_d_domain_init (struct sqpipe_pressurised_fixed_width_a_d_domain *d) {
-  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
-
-  sqpipe_domain_init(p);
-
-  static struct domain_vtable vtable = {
-    &sqpipe_pressurised_fixed_width_a_d_domain_new_edge,
-    &sqpipe_pressurised_fixed_width_a_d_domain_new_cell,
-    &domain_compute_fluxes_generic,
-    &domain_get_cell_generic,
-    &sqpipe_domain_destroy
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_pressurised_fixed_width_a_d_domain *sqpipe_pressurised_fixed_width_a_d_domain_new() {
-  struct sqpipe_pressurised_fixed_width_a_d_domain *d = (struct sqpipe_pressurised_fixed_width_a_d_domain *) sqpipe_domain_new();
-
-  sqpipe_pressurised_fixed_width_a_d_domain_init(d);
-
-  ((struct domain *) d)->vertex_values = (struct quantities *) sqpipe_pressurised_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->boundary_values = (struct quantities *) sqpipe_pressurised_fixed_width_a_d_quantities_new();
-  ((struct domain *) d)->explicit_update = (struct quantities *) sqpipe_pressurised_fixed_width_a_d_quantities_new();
-
-  return d;
-}
-/**********************************************************/
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d_comp_flux.c	(revision 8253)
+++ 	(revision )
@@ -1,51 +1,0 @@
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "math.h"
-#include <stdio.h>
-
-const double pi = 3.14159265358979;
-
-// Shared code snippets
-#include "util_ext.h"
-#include "solver.h"
-#include "python_solver.h"
-#include "sqpipe_pressurised_fixed_width_a_d.h"
-#include "sqpipe_python.h"
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  struct domain *D = (struct domain *) sqpipe_pressurised_fixed_width_a_d_domain_new();
-  double timestep;
-  PyObject *domain;
-
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dO", &timestep, &domain)) {
-    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
-    return NULL;
-  }
-
-  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
-
-  timestep = domain_compute_fluxes(D);
-  domain_destroy(D);
-  
-  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 initsqpipe_pressurised_fixed_width_a_d_comp_flux(void){
-  Py_InitModule("sqpipe_pressurised_fixed_width_a_d_comp_flux", MethodTable);
-  import_array();
-}
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_pressurised_fixed_width_a_d_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,113 +1,0 @@
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None):
-        conserved_quantities = ['area', 'discharge']
-        forcing_terms = []
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               tagged_elements      = tagged_elements)
-
-	self.__doc__ = 'sqpipe_area_discharge_domain'
-
-    def compute_fluxes(self):        
-        #Call correct module function
-        #(either from this module or C-extension)
-
-        # Set initial timestep to a large value
-        import sys
-        timestep = float(sys.maxint)
-
-        # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_pressurised_fixed_width_a_d_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep,self)        
-        
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-
-        #Calculate height (and fix negatives)better be non-negative!
-        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
-        w_C[:] = z_C + h_C
-        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
-
-        for name in ['velocity', 'stage' ]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        w_V = self.quantities['stage'].vertex_values
-        u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        a_V  = area.vertex_values
-        d_V  = discharge.vertex_values
-        h_V  = height.vertex_values
-
-        h_V[:,:]    = w_V - z_V
-
-        # Fix up small heights
-        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
-        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
-        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        # It may still be possible that h is small 
-        # If we set h to zero, we should also set u to 0
-        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
-        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
-
-        # Reconstruct conserved quantities and make everything consistent
-        # with new h values
-        w_V[:,:] = z_V + h_V
-        a_V[:,:] = h_V * b_V
-        d_V[:,:] = u_V * a_V
-
-
-        return
-
-    
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template.h
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template.h	(revision 8253)
+++ 	(revision )
@@ -1,216 +1,0 @@
-/*
-Template file for implementing a square pipe domain.
-Replace the following placeholders:
-
-<name> 
-
-The name of the class. For example, we have:
-
-free_surface_fixed_width_a_d
-
-which implements a free surface flow with a fixed width channel
-and the conserved quantities are a = area, d = discharge.
-
-<flux0>
-
-The formula for the conservation of mass flux.
-
-<flux1>
-
-The formula for the conservation of momentum flux.
-
-<speed>
-
-The speed of sound.
-
-<conserved0>
-
-The first conserved quantity.
-
-<conserved1>
-
-The second conserved quantity.
-
-<forcing_term0>
-
-Formula for first forcing term.
-
-<forcing_term1>
-
-Formula for second forcing term.
- */
-
-#include "sqpipe.h"
-
-/**********************************************************
-Square Pipe <name> quantity class
-**********************************************************/
-
-struct sqpipe_<name>_quantity {
-  struct sqpipe_quantity super;
-};
-
-// Construction
-void <name>_quantity_init (struct <name>_quantity *q);
-struct <name>_quantity *<name>_quantity_new();
-
-// Destruction
-
-// Methods
-double <name>_quantity_flux_formula0 (struct quantity *q, double g, double epsilon, double h0);
-double <name>_quantity_flux_formula1 (struct quantity *q, double g, double epsilon, double h0);
-double <name>_quantity_sound_speed (struct quantity *q, double g, double epsilon, double h0);
-double <name>_quantity_get_conserved (struct quantity *q, int k);
-
-// Implementation
-double <name>_quantity_flux_formula0 (struct quantity *q, double g, double epsilon, double h0) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-
-  return <flux0>;
-}
-double <name>_quantity_flux_formula1 (struct quantity *q, double g, double epsilon, double h0) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-
-  return <flux1>;
-}
-
-double <name>_quantity_sound_speed (struct quantity *q, double g, double epsilon, double h0) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-
-  return <speed>;
-}
-
-double <name>_quantity_get_conserved (struct quantity *q, int k) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = <conserved0>;
-    break;
-  case 1:
-    c = <conserved1>;
-    break;
-  default:
-    // Raise an error?
-    c = 0;
-  }
-  return c;
-}
-
-void <name>_quantity_init (struct <name>_quantity *q) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
-
-  static struct quantity_vtable vtable = {
-    &<name>_quantity_flux_formula0,
-    &<name>_quantity_flux_formula1,
-    &<name>_quantity_sound_speed,
-    &<name>_quantity_get_conserved
-  };
-  p->super.vtable = &vtable;
-}
-
-struct <name>_quantity *<name>_quantity_new() {
-  struct <name>_quantity *p = malloc(sizeof(struct <name>_quantity));
-  <name>_quantity_init(p);
-  return p;
-}
-//**********************************************************
-
-
-/**********************************************************
-Square pipe <name> cell class
-**********************************************************/
-struct sqpipe_<name>_cell {
-  struct sqpipe_cell super;
-};
-
-// Construction
-void sqpipe_<name>_cell_init (struct sqpipe_<name>_cell *c);
-struct sqpipe_<name>_cell *sqpipe_<name>_cell_new();
-
-// Destruction
-
-// Methods
-double sqpipe_<name>_cell_forcing_term (struct cell *c, int i, double normal, double g, double epsilon, double h0, double smin, double smax);
-
-// Implementation
-double sqpipe_<name>_cell_forcing_term (struct cell *c, int i, double normal, double g, double epsilon, double h0, double smin, double smax) {
-  int f;
-
-  switch (i) {
-  case 0:
-    f = <forcing_term0>;
-    break;
-  case 1:
-    f = <forcing_term1>;
-    break;
-  default:
-    f = 0.0;
-  }
-
-  return f;
-}
-
-void sqpipe_<name>_cell_init (struct sqpipe_<name>_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_godunov, // Generic method
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_<name>_cell_forcing_term
-  };
-  p->super.vtable = &vtable;
-}
-
-struct sqpipe_<name>_cell *sqpipe_<name>_cell_new() {
-  struct sqpipe_<name>_cell *p = malloc(sizeof(struct sqpipe_<name>_cell));
-  sqpipe_<name>_cell_init(p);
-  return p;
-}
-/**********************************************************/
-
-
-/**********************************************************
-Square Pipe <name> quantities class
-**********************************************************/
-struct <name>_quantities {
-  struct sqpipe_quantities super;
-};
-
-// Construction
-void <name>_quantities_init (struct <name>_quantities *q);
-struct <name>_quantities *<name>_quantities_new();
-
-// Destruction
-
-// Methods
-struct quantity* <name>_quantities_get_quantity(struct quantities *qs, int i, double normal, double g, double epsilon, double h0);
-
-// Implementation
-struct quantity* <name>_quantities_get_quantity(struct quantities *qs, int i, double normal, double g, double epsilon, double h0) {  
-  struct quantity *q = (struct quantity *) <name>_quantity_new();
-
-  sqpipe_quantities_set_quantity(qs, q, i, normal, g, epsilon, h0);
-  return q;
-}
-
-void <name>_quantities_init (struct <name>_quantities *q) {
-  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_python_get,
-    &<name>_quantities_get_quantity,
-    &sqpipe_quantities_new_cell,
-    &sqpipe_quantities_update
-  };
-  p->super.vtable = &vtable;
-}
-
-struct <name>_quantities *<name>_quantities_new() {
-  struct <name>_quantities *p = malloc(sizeof(struct <name>_quantities));
-  <name>_quantities_init(p);
-  return p;
-}
-/**********************************************************/
-
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template_comp_flux.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template_comp_flux.c	(revision 8253)
+++ 	(revision )
@@ -1,53 +1,0 @@
-/*
-Template file for implementing a square pipe domain.
-Replace the following placeholders:
-
-<name> 
-
-The name of the class. For example, we have:
-
-free_surface_fixed_width_a_d
-
-which implements a free surface flow with a fixed width channel
-and the conserved quantities are a = area, d = discharge.
-*/
-
-#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"
-#include "solver.h"
-#include "python_solver.h"
-#include "sqpipe_<name>.h"
-
-//=========================================================================
-// Python Glue
-//=========================================================================
-
-PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  struct quantities *vertex_values = (struct quantities *) sqpipe_<name>_quantities_new();
-  struct quantities *boundary_values = (struct quantities *) sqpipe_<name>_quantities_new();
-  struct quantities *explicit_update = (struct quantities *) sqpipe_<name>_quantities_new();
-
-  return compute_fluxes_generic(self, args, vertex_values, boundary_values, explicit_update);
-}
-
-//-------------------------------
-// Method table for python module
-//-------------------------------
-
-static struct PyMethodDef MethodTable[] = {
-  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
-  {NULL, NULL}
-};
-
-// Module initialisation
-void initsqpipe_<name>_comp_flux(void){
-  Py_InitModule("sqpipe_<name>_comp_flux", MethodTable);
-  import_array();
-}
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_template_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,130 +1,0 @@
-"""
-Template file for implementing a square pipe domain.
-Replace the following placeholders:
-
-<name> 
-
-The name of the class. For example, we have:
-
-free_surface_fixed_width_a_d
-
-which implements a free surface flow with a fixed width channel
-and the conserved quantities are a = area, d = discharge.
-
-<conserved> 
-
-List of conserved quantities.
-
-<centroids>
-
-Code to clean up centroids.
-
-<extrapolate>
-
-List of quantities to extrapolate from centroids to vertices.
-
-<extrapolated>
-
-Shortcut code for quantities that have been extrapolated.
-
-<update>
-
-Shortcut code for quantities that need to be updated.
-
-<reconstruct>
-
-Code to reconstruct update quantities from extrapolated quantities.
-May need to fix up small heights.
-"""
-
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None):
-        conserved_quantities = [<conserved>]
-        forcing_terms = [gravity]
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               tagged_elements      = tagged_elements)
-
-	self.__doc__ = 'sqpipe_<name>_domain'
-
-    def compute_fluxes(self):        
-        #Call correct module function
-        #(either from this module or C-extension)
-
-        # Set initial timestep to a large value
-        import sys
-        timestep = float(sys.maxint)
-
-        # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_<name>_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep,self)        
-        
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-        state     = self.quantities['state']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-        s_C   = state.centroid_values
-
-        <centroids>
-
-        for name in [<extrapolate>]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        <extrapolated>
-        # e.g.
-        #w_V = self.quantities['stage'].vertex_values
-        #u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        <update>
-        # e.g. 
-        #a_V  = area.vertex_values
-        #d_V  = discharge.vertex_values
-        #h_V  = height.vertex_values
-
-        <reconstruct>
-
-        return
-
-    
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_test_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_test_domain.py	(revision 8253)
+++ 	(revision )
@@ -1,106 +1,0 @@
-from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain
-from anuga_1d.sqpipe.sqpipe_forcing_terms import *
-from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
-
-
-class Domain(Sqpipe_domain):
-
-    def __init__(self, coordinates, boundary = None, tagged_elements = None, bulk_modulus = 1.0):
-        conserved_quantities = ['area', 'discharge']
-        forcing_terms = []
-        state = numpy.zeros(len(coordinates)-1, numpy.int)
-        Sqpipe_domain.__init__(self,
-                               coordinates          = coordinates,
-                               conserved_quantities = conserved_quantities,
-                               forcing_terms        = forcing_terms,
-                               boundary             = boundary,
-                               state                = state,
-                               bulk_modulus         = bulk_modulus,
-        #                       update_state_flag         = False,
-                               tagged_elements      = tagged_elements)
-
-
-        self.__doc__ = 'sqpipe_area_discharge_domain'
-            
-    def distribute_to_vertices_and_edges(self):
-        # Shortcuts
-        h0 = self.h0
-        epsilon = self.epsilon
-
-        area      = self.quantities['area']
-        discharge = self.quantities['discharge']
-        bed       = self.quantities['elevation']
-        height    = self.quantities['height']
-        velocity  = self.quantities['velocity']
-        width     = self.quantities['width']
-        top       = self.quantities['top']
-        stage     = self.quantities['stage']
-
-        #Arrays   
-        a_C   = area.centroid_values
-        d_C   = discharge.centroid_values
-        z_C   = bed.centroid_values
-        h_C   = height.centroid_values
-        u_C   = velocity.centroid_values
-        b_C   = width.centroid_values
-        t_C   = top.centroid_values
-        w_C   = stage.centroid_values
-
-        #Calculate height (and fix negatives)better be non-negative!
-        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
-        w_C[:] = z_C + h_C
-        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
-
-        for name in ['velocity', 'stage' ]:
-            Q = self.quantities[name]
-            if self.order == 1:
-                Q.extrapolate_first_order()
-            elif self.order == 2:
-                Q.extrapolate_second_order()
-            else:
-                raise 'Unknown order'
-
-        # These have been extrapolated
-        w_V = self.quantities['stage'].vertex_values
-        u_V = self.quantities['velocity'].vertex_values
-
-        # These are the given geometry and remain fixed
-        z_V  = bed.vertex_values
-        b_V  = width.vertex_values
-        t_V  = top.vertex_values
-
-        # These need to be updated
-        a_V  = area.vertex_values
-        d_V  = discharge.vertex_values
-        h_V  = height.vertex_values
-
-        h_V[:,:]    = w_V - z_V
-
-        # Fix up small heights
-        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
-        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
-        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
-
-        h_V[:,0] = h_0
-        h_V[:,1] = h_1
-
-        # It may still be possible that h is small 
-        # If we set h to zero, we should also set u to 0
-        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
-        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
-
-        # Reconstruct conserved quantities and make everything consistent
-        # with new h values
-        w_V[:,:] = z_V + h_V
-        a_V[:,:] = h_V * b_V
-        d_V[:,:] = u_V * a_V
-
-
-        return
-
-    
Index: /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py	(revision 8253)
+++ /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py	(revision 8254)
@@ -1,5 +1,5 @@
 import anuga_1d.sqpipe.parabolic_canal
 import time
-import anuga_1d.sqpipe.sqpipe_test_domain as dom
+import anuga_1d.sqpipe.sqpipe_domain as dom
 
 
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_dual_fixed_width_a_d.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_dual_fixed_width_a_d.py	(revision 8253)
+++ 	(revision )
@@ -1,10 +1,0 @@
-import anuga_1d.sqpipe.parabolic_canal
-import time
-import anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_domain as dom
-
-
-domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 100000.0
-yieldstep = 500.0
-anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
-anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_free_surface_fixed_width_a_d.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_free_surface_fixed_width_a_d.py	(revision 8253)
+++ 	(revision )
@@ -1,13 +1,0 @@
-import anuga_1d.sqpipe.parabolic_canal
-import time
-import anuga_1d.sqpipe.sqpipe_free_surface_fixed_width_a_d_domain as dom
-
-
-domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 10000.0
-yieldstep = 100.0
-
-anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
-anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
-#anuga_1d.sqpipe.parabolic_canal.evolve_domain(domain)
-#anuga_1d.sqpipe.parabolic_canal.write_domain(domain, 'outfile')
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_free_surface_fixed_width_w_d.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_free_surface_fixed_width_w_d.py	(revision 8253)
+++ 	(revision )
@@ -1,13 +1,0 @@
-import anuga_1d.sqpipe.parabolic_canal
-import time
-import anuga_1d.sqpipe.sqpipe_free_surface_fixed_width_w_d_domain as dom
-
-
-domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 10000.0
-yieldstep = 100.0
-
-anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
-anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
-#anuga_1d.sqpipe.parabolic_canal.evolve_domain(domain)
-#anuga_1d.sqpipe.parabolic_canal.write_domain(domain, 'outfile')
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_pressurised_fixed_width_a_d.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_pressurised_fixed_width_a_d.py	(revision 8253)
+++ 	(revision )
@@ -1,11 +1,0 @@
-import anuga_1d.sqpipe.parabolic_canal
-import time
-import anuga_1d.sqpipe.sqpipe_pressurised_fixed_width_a_d_domain as dom
-
-
-domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 100000.0
-yieldstep = 500.0
-
-anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
-anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
Index: unk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_template.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_template.py	(revision 8253)
+++ 	(revision )
@@ -1,170 +1,0 @@
-"""
-Template file for implementing a square pipe domain.
-Replace the following placeholders:
-
-<name> 
-
-The name of the class. For example, we have:
-
-free_surface_fixed_width_a_d
-
-which implements a free surface flow with a fixed width channel
-and the conserved quantities are a = area, d = discharge.
-
-<conserved0>
-
-The first conserved quantity e.g. height, area, stage.
-"""
-import os
-from math import sqrt, pi, sin, cos
-import numpy
-import time
-
-
-from anuga_1d.sqpipe.sqpipe_<name>_domain import *
-from anuga_1d.config import g, epsilon
-from anuga_1d.base.generic_mesh import uniform_mesh
-
-#===============================================================================
-# setup problem
-#===============================================================================
-
-
-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 = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
-
-def analytic_canal(x,t):
-    u = numpy.zeros_like(x)    ## water velocity
-    h = numpy.zeros_like(x)    ## water depth
-
-    ## Define Basin Bathymetry
-    z = numpy.zeros_like(x) ## elevation of basin
-    w = numpy.zeros_like(x)   ## elevation of water surface
-
-    z[:] = z_infty*(x**2/L_x**2)
-    u[:] = -A0*omega*sin(omega*t)
-    w[:] = numpy.maximum(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t)),z)
-    h[:] = numpy.maximum(w-z, 0.0)
-
-    T = 2.0*pi/omega
-
-    return u,h,w,z, T
-
-
-def stage(x):
-    t=0.0
-    u,h,w,z,T = analytic_canal(x,t)
-
-    return w
-
-def elevation(x):
-    t=0.0
-    u,h,w,z,T = analytic_canal(x,t)
-
-    return z
-
-
-def height(x):
-    t=0.0
-    u,h,w,z,T = analytic_canal(x,t)
-
-    return h
-
-def width(x):
-    return numpy.ones_like(x)
-
-def top(x):
-    return 4.0 * numpy.ones_like(x)
-
-def state(x):
-    return numpy.ones_like(x)
-
-def area(x):
-    return height(x)*width(x)
-
-
-#===============================================================================
-finaltime = 500.0
-yieldstep = 10.0
-
-
-
-N = 100
-print "Evaluating domain with %d cells" %N
-domain = Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x))
-
-domain.set_spatial_order(2)
-domain.set_timestepping_method('rk2')
-domain.set_CFL(1.0)
-domain.set_limiter("vanleer")
-
-domain.set_beta(1.0)
-
-domain.set_quantity('<conserved0>', <conserved0>)
-domain.set_quantity('elevation',elevation)
-domain.set_quantity('width',width)
-domain.set_quantity('top',top)
-domain.set_quantity('state',state)
-
-Br = Reflective_boundary(domain)
-
-domain.set_boundary({'left': Br, 'right' : Br})
-
-t0 = time.time()
-
-for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-    domain.write_time()
-
-print 'That took %.2f seconds' %(time.time()-t0)
-
-
-N = float(N)
-HeightC    = domain.quantities['height'].centroid_values
-StageC     = domain.quantities['stage'].centroid_values
-BedC       = domain.quantities['elevation'].centroid_values
-C          = domain.centroids
-
-HeightV    = domain.quantities['height'].vertex_values
-StageV     = domain.quantities['stage'].vertex_values
-BedV       = domain.quantities['elevation'].vertex_values
-VelocityV  = domain.quantities['velocity'].vertex_values
-X          = domain.vertices
-
-
-import pylab
-
-pylab.hold(False)
-
-plot1 = pylab.subplot(311)
-
-pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
-
-plot1.set_ylim([-1,40])
-
-pylab.xlabel('Position')
-pylab.ylabel('Stage')
-pylab.legend(('Analytical Solution', 'Numerical Solution'),
-           'upper right', shadow=True)
-
-
-plot2 = pylab.subplot(312)
-
-pylab.plot(X.flat,HeightV.flat)
-
-plot2.set_ylim([-1,10])
-
-pylab.xlabel('Position')
-pylab.ylabel('Height')
-
-plot3 = pylab.subplot(313)
-
-pylab.plot(X.flat,VelocityV.flat)
-plot3.set_ylim([-15,15])
-
-pylab.xlabel('Position')
-pylab.ylabel('Velocity')
-
-pylab.show()
-
