Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5161)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5162)
@@ -10,4 +10,5 @@
 from Numeric import allclose, argmax, zeros, Float
 from anuga.config import epsilon
+from anuga.config import beta_euler, beta_rk2
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -301,4 +302,12 @@
         self.default_order = n
         self._order_ = self.default_order
+
+        if self.default_order == 1:
+            self.set_timestepping_method('euler')
+            #self.set_all_limiters(beta_euler)
+
+        if self.default_order == 2:
+            self.set_timestepping_method('rk2')
+            #self.set_all_limiters(beta_rk2)
         
 
@@ -1550,8 +1559,8 @@
             elif self._order_ == 2:
                 Q.extrapolate_second_order()
-                Q.limit()
+                #Q.limit()
             else:
                 raise 'Unknown order'
-            Q.interpolate_from_vertices_to_edges()
+            #Q.interpolate_from_vertices_to_edges()
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5161)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5162)
@@ -62,5 +62,8 @@
         #Allocate space for Gradient
         self.x_gradient = zeros(N, Float)
-        self.y_gradient = zeros(N, Float)        
+        self.y_gradient = zeros(N, Float)
+
+        #Allocate space for Limiter Phi
+        self.phi = zeros(N, Float)        
 
         #Intialise centroid and edge_values
@@ -74,8 +77,10 @@
         #flux calculations and forcing functions
 
-        N = len(domain) # number_of_triangles
+        #Allocate space for update fields
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
-        self.centroid_backup_values = zeros(N, Float)       
+        self.centroid_backup_values = zeros(N, Float)
+
+        self.beta = 1.0
 
 
@@ -1403,9 +1408,22 @@
         compute_gradients(self)
         extrapolate_from_gradient(self)
-
+        
     def extrapolate_second_order_and_limit(self):
         #Call correct module function
         #(either from this module or C-extension)
-        extrapolate_second_order_and_limit(self)        
+        extrapolate_second_order_and_limit(self)
+
+    def bound_vertices_below_by_constant(self, bound):
+        #Call correct module function
+        #(either from this module or C-extension)
+        bound_vertices_below_by_constant(self, bound)
+
+    def bound_vertices_below_by_quantity(self, quantity):
+        #Call correct module function
+        #(either from this module or C-extension)
+
+        #check consistency
+        assert self.domain == quantity.domain
+        bound_vertices_below_by_quantity(self, quantity)                        
 
     def backup_centroid_values(self):
@@ -1450,4 +1468,7 @@
          limit_gradient_by_neighbour,\
          extrapolate_from_gradient,\
+         extrapolate_second_order_and_limit,\
+         bound_vertices_below_by_constant,\
+         bound_vertices_below_by_quantity,\
          interpolate_from_vertices_to_edges,\
          interpolate_from_edges_to_vertices,\
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5161)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5162)
@@ -145,4 +145,5 @@
 					 double* vertex_values,
 					 double* edge_values,
+					 double* phi,
 					 double* x_gradient,
 					 double* y_gradient) {
@@ -151,7 +152,7 @@
 	double x, y, x0, y0, x1, y1, x2, y2;
 	long n;
-	double qmin, qmax, qn, qc;
-	double dq, dqa[3], phi, r;
-
+	double qmin, qmax, qc;
+	double qn[3];
+	double dq, dqa[3], r;
 
 	for (k=0; k<N; k++){
@@ -161,9 +162,10 @@
 
 		// Centroid coordinates
-		x = centroids[k2]; y = centroids[k2+1];
+		x = centroids[k2+0]; 
+		y = centroids[k2+1];
 
 		// vertex coordinates
 		// x0, y0, x1, y1, x2, y2 = X[k,:]
-		x0 = vertex_coordinates[k6 + 0];
+	        x0 = vertex_coordinates[k6 + 0];
 		y0 = vertex_coordinates[k6 + 1];
 		x1 = vertex_coordinates[k6 + 2];
@@ -172,7 +174,4 @@
 		y2 = vertex_coordinates[k6 + 5];
 
-		//Centroid Value
-		qc = centroid_values[k];
-
 		// Extrapolate to Vertices
 		vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
@@ -184,7 +183,44 @@
 		edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
 		edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
-
-
-		phi = 1.0;
+	}
+
+
+
+	for (k=0; k<N; k++){
+		k6 = 6*k;
+		k3 = 3*k;
+		k2 = 2*k;
+
+
+		qc = centroid_values[k];
+		
+		qmin = qc;
+		qmax = qc;
+		
+		for (i=0; i<3; i++) {
+		    n = neighbours[k3+i];
+		    if (n < 0) {
+			qn[i] = qc;
+		    } else {
+			qn[i] = centroid_values[n];
+		    }
+
+		    qmin = min(qmin, qn[i]);
+		    qmax = max(qmax, qn[i]);
+		}
+
+		//qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
+		//qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
+
+/* 		for (i=0; i<3; i++) { */
+/* 		    n = neighbours[k3+i]; */
+/* 		    if (n < 0) { */
+/* 			qn[i] = qc; */
+/* 			qmin[i] = qtmin; */
+/* 			qmax[i] = qtmax; */
+/* 		    }  */
+/* 		} */
+			
+		phi[k] = 1.0;
 
 		for (i=0; i<3; i++) {
@@ -192,29 +228,22 @@
 		    dqa[i] = dq;                      //Save dq for use in updating vertex values
 
-		    n = neighbours[k3+i];
-		    if (n >= 0) {
-			qn = centroid_values[n];     //Neighbour's centroid value
-
-			qmin = min(qc, qn);
-			qmax = max(qc, qn);
-
-
-			r = 1.0;
+		    r = 1.0;
       
-			if (dq > 0.0) r = (qmax - qc)/dq;
-			if (dq < 0.0) r = (qmin - qc)/dq;      
+		    if (dq > 0.0) r = (qmax - qc)/dq;
+		    if (dq < 0.0) r = (qmin - qc)/dq;
 		    
-			phi = min( min(r*beta, 1.0), phi);    
-		    }
+		    phi[k] = min( min(r*beta, 1.0), phi[k]);
+		    
 		}
 
 
+
 		//Update gradient, edge and vertex values using phi limiter
-		x_gradient[k] = x_gradient[k]*phi;
-		y_gradient[k] = y_gradient[k]*phi;
-
-		edge_values[k3+0] = qc + phi*dqa[0];
-		edge_values[k3+1] = qc + phi*dqa[1];
-		edge_values[k3+2] = qc + phi*dqa[2];
+		x_gradient[k] = x_gradient[k]*phi[k];
+		y_gradient[k] = y_gradient[k]*phi[k];
+
+		edge_values[k3+0] = qc + phi[k]*dqa[0];
+		edge_values[k3+1] = qc + phi[k]*dqa[1];
+		edge_values[k3+2] = qc + phi[k]*dqa[2];
 
 
@@ -222,4 +251,5 @@
 		vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
 		vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+		
 
 	}
@@ -290,9 +320,14 @@
 }
 
+
+
+
 int _limit_edges_by_all_neighbours(int N, double beta,
-		     double* centroid_values,
-		     double* vertex_values,
-		     double* edge_values,
-		     long*   neighbours) {
+				   double* centroid_values,
+				   double* vertex_values,
+				   double* edge_values,
+				   long*   neighbours,
+				   double* x_gradient,
+				   double* y_gradient) {
 
 	int i, k, k2, k3, k6;
@@ -334,5 +369,8 @@
 		}
     
-		//Update vertex and edge values using phi limiter
+		//Update gradient, vertex and edge values using phi limiter
+		x_gradient[k] = x_gradient[k]*phi;
+		y_gradient[k] = y_gradient[k]*phi;
+
 		edge_values[k3+0] = qc + phi*dqa[0];
 		edge_values[k3+1] = qc + phi*dqa[1];
@@ -353,4 +391,61 @@
 		     double* vertex_values,
 		     double* edge_values,
+		     long*   neighbours) {
+
+	int i, k, k2, k3, k6;
+	long n;
+	double qmin, qmax, qn, qc;
+	double dq, dqa[3], phi, r;
+
+	for (k=0; k<N; k++){
+		k6 = 6*k;
+		k3 = 3*k;
+		k2 = 2*k;
+
+		qc = centroid_values[k];
+		phi = 1.0;
+
+		for (i=0; i<3; i++) {
+		    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+		    dqa[i] = dq;                      //Save dqa for use in updating vertex values
+
+		    n = neighbours[k3+i];
+		    qn = qc;
+		    if (n >= 0)  qn = centroid_values[n]; //Neighbour's centroid value
+
+		    qmin = min(qc, qn);
+		    qmax = max(qc, qn);
+
+		    r = 1.0;
+      
+		    if (dq > 0.0) r = (qmax - qc)/dq;
+		    if (dq < 0.0) r = (qmin - qc)/dq;      
+		    
+		    phi = min( min(r*beta, 1.0), phi);    
+		   
+		}
+
+
+		//Update edge and vertex values using phi limiter
+		edge_values[k3+0] = qc + phi*dqa[0];
+		edge_values[k3+1] = qc + phi*dqa[1];
+		edge_values[k3+2] = qc + phi*dqa[2];
+		
+		vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
+		vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
+		vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+
+	}
+
+	return 0;
+}
+
+
+int _limit_gradient_by_neighbour(int N, double beta,
+		     double* centroid_values,
+		     double* vertex_values,
+		     double* edge_values,
+		     double* x_gradient,
+		     double* y_gradient,
 		     long*   neighbours) {
 
@@ -403,16 +498,13 @@
 }
 
-
-int _limit_gradient_by_neighbour(int N, double beta,
+int _bound_vertices_below_by_constant(int N, double bound,
 		     double* centroid_values,
 		     double* vertex_values,
 		     double* edge_values,
 		     double* x_gradient,
-		     double* y_gradient,
-		     long*   neighbours) {
+		     double* y_gradient) {
 
 	int i, k, k2, k3, k6;
-	long n;
-	double qmin, qmax, qn, qc;
+	double qmin, qc;
 	double dq, dqa[3], phi, r;
 
@@ -423,35 +515,32 @@
 
 		qc = centroid_values[k];
+		qmin = bound;
+
+
 		phi = 1.0;
-
-		for (i=0; i<3; i++) {
-		    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+		for (i=0; i<3; i++) {    
+		    r = 1.0;
+      
+		    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
 		    dqa[i] = dq;                      //Save dq for use in updating vertex values
-
-		    n = neighbours[k3+i];
-		    if (n >= 0) {
-			qn = centroid_values[n]; //Neighbour's centroid value
-
-			qmin = min(qc, qn);
-			qmax = max(qc, qn);
-
-			r = 1.0;
       
-			if (dq > 0.0) r = (qmax - qc)/dq;
-			if (dq < 0.0) r = (qmin - qc)/dq;      
+		    if (dq < 0.0) r = (qmin - qc)/dq;      
+  
 		    
-			phi = min( min(r*beta, 1.0), phi);    
-		    }
+		    phi = min( min(r, 1.0), phi);    
 		}
-
-
-		//Update edge and vertex values using phi limiter
-		edge_values[k3+0] = qc + phi*dqa[0];
-		edge_values[k3+1] = qc + phi*dqa[1];
-		edge_values[k3+2] = qc + phi*dqa[2];
+    
+
+		//Update gradient, vertex and edge values using phi limiter
+		x_gradient[k] = x_gradient[k]*phi;
+		y_gradient[k] = y_gradient[k]*phi;
+
+		vertex_values[k3+0] = qc + phi*dqa[0];
+		vertex_values[k3+1] = qc + phi*dqa[1];
+		vertex_values[k3+2] = qc + phi*dqa[2];
 		
-		vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
-		vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
-		vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+		edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
+		edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
+		edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
 
 	}
@@ -460,5 +549,54 @@
 }
 
-
+int _bound_vertices_below_by_quantity(int N,
+				      double* bound_vertex_values,
+				      double* centroid_values,
+				      double* vertex_values,
+				      double* edge_values,
+				      double* x_gradient,
+				      double* y_gradient) {
+
+	int i, k, k2, k3, k6;
+	double qmin, qc;
+	double dq, dqa[3], phi, r;
+
+	for (k=0; k<N; k++){
+		k6 = 6*k;
+		k3 = 3*k;
+		k2 = 2*k;
+
+		qc = centroid_values[k];
+
+		phi = 1.0;
+		for (i=0; i<3; i++) {    
+		    r = 1.0;
+      
+		    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
+		    dqa[i] = dq;                      //Save dq for use in updating vertex values
+      
+		    qmin = bound_vertex_values[k3+i];
+		    if (dq < 0.0) r = (qmin - qc)/dq;      
+  
+		    
+		    phi = min( min(r, 1.0), phi);    
+		}
+    
+
+		//Update gradient, vertex and edge values using phi limiter
+		x_gradient[k] = x_gradient[k]*phi;
+		y_gradient[k] = y_gradient[k]*phi;
+
+		vertex_values[k3+0] = qc + phi*dqa[0];
+		vertex_values[k3+1] = qc + phi*dqa[1];
+		vertex_values[k3+2] = qc + phi*dqa[2];
+		
+		edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
+		edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
+		edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
+
+	}
+
+	return 0;
+}
 
 int _interpolate_from_vertices_to_edges(int N,
@@ -906,83 +1044,4 @@
 
 
-/* PyObject *compute_gradients(PyObject *self, PyObject *args) { */
-/*   //"""Compute gradients of triangle surfaces defined by centroids of */
-/*   //neighbouring volumes. */
-/*   //If one edge is on the boundary, use own centroid as neighbour centroid. */
-/*   //If two or more are on the boundary, fall back to first order scheme. */
-/*   //""" */
-
-
-/* 	PyObject *quantity, *domain, *R; */
-/* 	PyArrayObject */
-/* 		*centroids,            //Coordinates at centroids */
-/* 		*centroid_values,      //Values at centroids */
-/* 		*number_of_boundaries, //Number of boundaries for each triangle */
-/* 		*surrogate_neighbours, //True neighbours or - if one missing - self */
-/* 		*a, *b;                //Return values */
-
-/* 	int dimensions[1], N, err; */
-
-/* 	// Convert Python arguments to C */
-/* 	if (!PyArg_ParseTuple(args, "O", &quantity)) { */
-/* 	  PyErr_SetString(PyExc_RuntimeError,  */
-/* 			  "quantity_ext.c: compute_gradients could not parse input");	 */
-/* 	  return NULL; */
-/* 	} */
-
-/* 	domain = PyObject_GetAttrString(quantity, "domain"); */
-/* 	if (!domain) { */
-/* 	  PyErr_SetString(PyExc_RuntimeError,  */
-/* 			  "compute_gradients could not obtain domain object from quantity"); */
-/* 	  return NULL; */
-/* 	} */
-
-/* 	// Get pertinent variables */
-
-/* 	centroids = get_consecutive_array(domain, "centroid_coordinates"); */
-/* 	centroid_values = get_consecutive_array(quantity, "centroid_values"); */
-/* 	surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); */
-/* 	number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); */
-
-/* 	N = centroid_values -> dimensions[0]; */
-
-/* 	// Release */
-/* 	Py_DECREF(domain); */
-
-/* 	// Allocate space for return vectors a and b (don't DECREF) */
-/* 	dimensions[0] = N; */
-/* 	a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
-/* 	b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
-
-
-
-/* 	err = _compute_gradients(N, */
-/* 			(double*) centroids -> data, */
-/* 			(double*) centroid_values -> data, */
-/* 			(long*) number_of_boundaries -> data, */
-/* 			(long*) surrogate_neighbours -> data, */
-/* 			(double*) a -> data, */
-/* 			(double*) b -> data); */
-
-/* 	if (err != 0) { */
-/* 	  PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
-/* 	  return NULL; */
-/* 	} */
-
-/* 	// Release */
-/* 	Py_DECREF(centroids); */
-/* 	Py_DECREF(centroid_values); */
-/* 	Py_DECREF(number_of_boundaries); */
-/* 	Py_DECREF(surrogate_neighbours); */
-
-/* 	// Build result, release and return */
-/* 	R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); */
-/* 	Py_DECREF(a); */
-/* 	Py_DECREF(b); */
-/* 	return R; */
-/* } */
-
-
-
 PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
 
@@ -1087,4 +1146,5 @@
 	*quantity_vertex_values,        //Values at vertices
 	*quantity_edge_values,          //Values at edges
+	*quantity_phi,                  //limiter phi values
 	*quantity_x_gradient,           //x gradient
 	*quantity_y_gradient;           //y gradient
@@ -1097,10 +1157,14 @@
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOd", 
-			&domain,
-			&quantity,
-			&beta)) {
+  if (!PyArg_ParseTuple(args, "O",&quantity)) {
       PyErr_SetString(PyExc_RuntimeError, 
 		      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");	
+      return NULL;
+  }
+
+  domain = PyObject_GetAttrString(quantity, "domain");
+  if (!domain) {
+      PyErr_SetString(PyExc_RuntimeError, 
+		      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");	
       return NULL;
   }
@@ -1117,10 +1181,11 @@
   quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
   quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_phi                  = get_consecutive_array(quantity, "phi");
   quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
   quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
 
+  beta = get_python_double(quantity,"beta");
+
   ntri = quantity_centroid_values -> dimensions[0];
-
-
 
   err = _compute_gradients(ntri,
@@ -1139,17 +1204,31 @@
 
 
-  err = _extrapolate_and_limit_from_gradient(ntri, beta,
-					     (double*) domain_centroids -> data,
-					     (long*)   domain_neighbours -> data,
-					     (double*) quantity_centroid_values -> data,
-					     (double*) domain_vertex_coordinates -> data,
-					     (double*) quantity_vertex_values -> data,
-					     (double*) quantity_edge_values -> data,
-					     (double*) quantity_x_gradient -> data,
-					     (double*) quantity_y_gradient -> data);
+  err = _extrapolate_from_gradient(ntri, 
+				   (double*) domain_centroids -> data,
+				   (double*) quantity_centroid_values -> data,
+				   (double*) domain_vertex_coordinates -> data,
+				   (double*) quantity_vertex_values -> data,
+				   (double*) quantity_edge_values -> data,
+				   (double*) quantity_x_gradient -> data,
+				   (double*) quantity_y_gradient -> data);
 
   if (err != 0) {
       PyErr_SetString(PyExc_RuntimeError,
-			  "quantity_ext.c: Internal function _extrapolate_and_limit_from_gradient failed");
+			  "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
+      return NULL;
+  }
+
+
+  err = _limit_edges_by_all_neighbours(ntri, beta,
+				       (double*) quantity_centroid_values -> data,
+				       (double*) quantity_vertex_values -> data,
+				       (double*) quantity_edge_values -> data,
+				       (long*)   domain_neighbours -> data,
+				       (double*) quantity_x_gradient -> data,
+				       (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _limit_edges_by_all_neighbours failed");
       return NULL;
   }
@@ -1165,4 +1244,5 @@
   Py_DECREF(quantity_vertex_values);
   Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_phi);
   Py_DECREF(quantity_x_gradient);
   Py_DECREF(quantity_y_gradient);
@@ -1445,9 +1525,11 @@
   //
 
-	PyObject *quantity, *domain, *Tmp;
+	PyObject *quantity, *domain;
 	PyArrayObject
 	    *vertex_values,   //Conserved quantities at vertices
 	    *centroid_values, //Conserved quantities at centroids
 	    *edge_values,     //Conserved quantities at edges
+	    *x_gradient,
+	    *y_gradient,
 	    *neighbours;
 
@@ -1463,20 +1545,5 @@
 	}
 
-	domain = PyObject_GetAttrString(quantity, "domain");
-	if (!domain) {
-	  PyErr_SetString(PyExc_RuntimeError, 
-			  "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity");		  	
-	  
-	  return NULL;
-	}
-
-	// Get safety factor beta_w
-	Tmp = PyObject_GetAttrString(domain, "beta_w");
-	if (!Tmp) {
-	  PyErr_SetString(PyExc_RuntimeError, 
-			  "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain");		  	
-	  
-	  return NULL;
-	}	
+	domain = get_python_object(quantity, "domain");
 
 
@@ -1486,6 +1553,9 @@
 	vertex_values    = get_consecutive_array(quantity, "vertex_values");
 	edge_values      = get_consecutive_array(quantity, "edge_values");
-	beta_w           = PyFloat_AsDouble(Tmp);
-
+	x_gradient       = get_consecutive_array(quantity, "x_gradient");
+	y_gradient       = get_consecutive_array(quantity, "y_gradient");
+	beta_w           = get_python_double(domain,"beta_w");
+
+	Py_DECREF(domain);
 
 	N = centroid_values -> dimensions[0];
@@ -1495,9 +1565,11 @@
 					     (double*) vertex_values -> data,
 					     (double*) edge_values -> data,
-					     (long*)   neighbours -> data);
+					     (long*)   neighbours -> data,
+					     (double*) x_gradient -> data,
+					     (double*) y_gradient -> data);
 
 	if (err != 0) {
 	  PyErr_SetString(PyExc_RuntimeError,
-			  "Internal function _limit_by_vertex failed");
+			  "quantity_ect.c: limit_edges_by_all_neighbours internal function _limit_edges_by_all_neighbours failed");
 	  return NULL;
 	}	
@@ -1509,5 +1581,145 @@
 	Py_DECREF(vertex_values);
 	Py_DECREF(edge_values);
-	Py_DECREF(Tmp);
+
+
+
+	return Py_BuildValue("");
+}
+
+PyObject *bound_vertices_below_by_constant(PyObject *self, PyObject *args) {
+  //Bound a quantity below by a contant (useful for ensuring positivity
+  //precondition:
+  //  vertex values are already calulated, gradient consistent
+  //postcondition:
+  //  gradient, vertex and edge values are updated
+  //
+
+	PyObject *quantity, *domain;
+	PyArrayObject
+	    *vertex_values,   //Conserved quantities at vertices
+	    *centroid_values, //Conserved quantities at centroids
+	    *edge_values,     //Conserved quantities at edges
+	    *x_gradient,
+	    *y_gradient;
+
+	double bound; //Safety factor
+	int N, err;
+
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "Od", &quantity, &bound)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: bound_vertices_below_by_constant could not parse input");
+	  return NULL;
+	}
+
+	domain = get_python_object(quantity, "domain");
+
+	// Get pertinent variables
+	centroid_values  = get_consecutive_array(quantity, "centroid_values");
+	vertex_values    = get_consecutive_array(quantity, "vertex_values");
+	edge_values      = get_consecutive_array(quantity, "edge_values");
+	x_gradient       = get_consecutive_array(quantity, "x_gradient");
+	y_gradient       = get_consecutive_array(quantity, "y_gradient");
+
+
+
+	Py_DECREF(domain);
+
+	N = centroid_values -> dimensions[0];
+
+	err = _bound_vertices_below_by_constant(N, bound,
+					     (double*) centroid_values -> data,
+					     (double*) vertex_values -> data,
+					     (double*) edge_values -> data,
+					     (double*) x_gradient -> data,
+					     (double*) y_gradient -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ect.c: bound_vertices_below_by_constant internal function _bound_vertices_below_by_constant failed");
+	  return NULL;
+	}	
+
+
+	// Release
+	Py_DECREF(centroid_values);
+	Py_DECREF(vertex_values);
+	Py_DECREF(edge_values);
+	Py_DECREF(x_gradient);
+	Py_DECREF(y_gradient);
+
+
+
+	return Py_BuildValue("");
+}
+
+
+PyObject *bound_vertices_below_by_quantity(PyObject *self, PyObject *args) {
+  //Bound a quantity below by a contant (useful for ensuring positivity
+  //precondition:
+  //  vertex values are already calulated, gradient consistent
+  //postcondition:
+  //  gradient, vertex and edge values are updated
+  //
+
+	PyObject *quantity, *bounding_quantity, *domain;
+	PyArrayObject
+	    *vertex_values,   //Conserved quantities at vertices
+	    *centroid_values, //Conserved quantities at centroids
+	    *edge_values,     //Conserved quantities at edges
+	    *x_gradient,
+	    *y_gradient,
+	    *bound_vertex_values;
+
+	int N, err;
+
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "OO", &quantity, &bounding_quantity)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: bound_vertices_below_by_quantity could not parse input");
+	  return NULL;
+	}
+
+	domain = get_python_object(quantity, "domain");
+
+	// Get pertinent variables
+	centroid_values     = get_consecutive_array(quantity, "centroid_values");
+	vertex_values       = get_consecutive_array(quantity, "vertex_values");
+	edge_values         = get_consecutive_array(quantity, "edge_values");
+	x_gradient          = get_consecutive_array(quantity, "x_gradient");
+	y_gradient          = get_consecutive_array(quantity, "y_gradient");
+	bound_vertex_values = get_consecutive_array(bounding_quantity, "vertex_values");
+
+
+
+	Py_DECREF(domain);
+
+	N = centroid_values -> dimensions[0];
+
+	err = _bound_vertices_below_by_quantity(N, 
+						(double*) bound_vertex_values -> data,
+						(double*) centroid_values -> data,
+						(double*) vertex_values -> data,
+						(double*) edge_values -> data,
+						(double*) x_gradient -> data,
+						(double*) y_gradient -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ect.c: bound_vertices_below_by_quantity internal function _bound_vertices_below_by_quantity failed");
+	  return NULL;
+	}	
+
+
+	// Release
+	Py_DECREF(centroid_values);
+	Py_DECREF(vertex_values);
+	Py_DECREF(edge_values);
+	Py_DECREF(x_gradient);
+	Py_DECREF(y_gradient);
+	Py_DECREF(bound_vertex_values);
+
 
 
@@ -1589,4 +1801,5 @@
 
 	// Release
+	Py_DECREF(domain);
 	Py_DECREF(neighbours);
 	Py_DECREF(centroid_values);
@@ -1700,4 +1913,6 @@
 	{"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
 	{"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
+	{"bound_vertices_below_by_constant", bound_vertices_below_by_constant, METH_VARARGS, "Print out"},
+	{"bound_vertices_below_by_quantity", bound_vertices_below_by_quantity, METH_VARARGS, "Print out"},
 	{"update", update, METH_VARARGS, "Print out"},
 	{"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
Index: /anuga_core/source/anuga/advection/advection_ext.c
===================================================================
--- /anuga_core/source/anuga/advection/advection_ext.c	(revision 5161)
+++ /anuga_core/source/anuga/advection/advection_ext.c	(revision 5162)
@@ -2,7 +2,6 @@
 //
 // To compile (Python2.X):
-// use python ../utilities/compile.py
-//
-// See the module advection.py
+// use python ../utilities/compile.py to
+// compile all C files within a directory
 //
 //
@@ -19,4 +18,7 @@
 
 
+//-------------------------------------------
+// Low level routines (called from wrappers)
+//------------------------------------------
 
 double _compute_fluxes(
@@ -38,5 +40,5 @@
 
  
-        //Loop
+        //Local Variables
 
         double qr,ql;
@@ -51,4 +53,6 @@
 	int k3;
 	
+	//Loop through triangles
+
 	timestep = max_timestep;
 
@@ -59,5 +63,5 @@
             for (i=0; i<3; i++){
 	        k_i = k3+i;
-                //Quantities inside volume facing neighbour i
+                //Quantities inside triangle facing neighbour i
                 ql = quantity_edge[k_i];
 
@@ -144,4 +148,5 @@
 
   PyObject *domain, *quantity;
+ 
   PyArrayObject 
     * quantity_update,
Index: /anuga_core/source/anuga/config.py
===================================================================
--- /anuga_core/source/anuga/config.py	(revision 5161)
+++ /anuga_core/source/anuga/config.py	(revision 5162)
@@ -100,4 +100,13 @@
 #timestepping_method = 'rk2'   # 2nd Order TVD scheme
 
+# rk2 is a little more stable than euler, so rk2 timestepping
+# can deal with a larger beta when slope limiting the reconstructed
+# solution. The large beta is needed if solving problems sensitive
+# to numerical diffusion, like a small forced wave in an ocean
+beta_euler = 1.0
+beta_rk2   = 1.6
+
+
+
 # Option to search for signatures where isolated triangles are
 # responsible for a small global timestep.
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5161)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5162)
@@ -103,4 +103,5 @@
 from anuga.config import alpha_balance
 from anuga.config import optimise_dry_cells
+from anuga.config import optimised_gradient_limiter
 
 #---------------------
@@ -176,5 +177,9 @@
         self.minimum_storable_height = minimum_storable_height
         self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
-        
+
+        # Limiters
+        self.use_old_limiter = True
+
+        self.optimised_gradient_limiter = optimised_gradient_limiter
                 
 
@@ -187,8 +192,14 @@
         self.beta_w      = beta
         self.beta_w_dry  = beta
+        self.quantities['stage'].beta = beta
+        
         self.beta_uh     = beta
         self.beta_uh_dry = beta
+        self.quantities['xmomentum'].beta = beta
+        
         self.beta_vh     = beta
         self.beta_vh_dry = beta
+        self.quantities['ymomentum'].beta = beta
+        
         self.beta_h      = beta
         
@@ -384,5 +395,18 @@
         # Call correct module function
         # (either from this module or C-extension)
-        distribute_to_vertices_and_edges(self)
+        if self.use_old_limiter:
+            distribute_to_vertices_and_edges(self)
+        else:
+            for name in self.conserved_quantities:
+                Q = self.quantities[name]
+                if self._order_ == 1:
+                    Q.extrapolate_first_order()
+                elif self._order_ == 2:
+                    Q.extrapolate_second_order_and_limit()
+                    if name == 'stage':
+                        Q.bound_vertices_below_by_quantity(self.quantities['elevation'])
+                else:
+                    raise 'Unknown order'
+ 
 
 
@@ -399,8 +423,8 @@
         # self.check_integrity()
 
-        msg = 'Parameter beta_h must be in the interval [0, 1['
-        assert 0 <= self.beta_h <= 1.0, msg
-        msg = 'Parameter beta_w must be in the interval [0, 1['
-        assert 0 <= self.beta_w <= 1.0, msg
+        msg = 'Parameter beta_h must be in the interval [0, 2['
+        assert 0 <= self.beta_h <= 2.0, msg
+        msg = 'Parameter beta_w must be in the interval [0, 2['
+        assert 0 <= self.beta_w <= 2.0, msg
 
 
@@ -737,5 +761,5 @@
     """
 
-    from anuga.config import optimised_gradient_limiter
+    
 
     # Remove very thin layers of water
@@ -743,5 +767,5 @@
 
     # Extrapolate all conserved quantities
-    if optimised_gradient_limiter:
+    if domain.optimised_gradient_limiter:
         # MH090605 if second order,
         # perform the extrapolation and limiting on
@@ -764,16 +788,4 @@
                 Q.extrapolate_first_order()
             elif domain._order_ == 2:
-
-                # Experiment
-                #if name == 'stage':
-                #    #print name, 'second'
-                #    Q.extrapolate_second_order()
-                #    Q.limit()
-                #else:
-                #    #print name, 'first'                
-                #    Q.extrapolate_first_order()
-                #    #Q.extrapolate_second_order()
-                #    #Q.limit()                
-                
                 Q.extrapolate_second_order()
                 Q.limit()
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5161)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5162)
@@ -1088,5 +1088,5 @@
   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
   double dqv[3], qmin, qmax, hmin, hmax;
-  double hc, h0, h1, h2, beta_tmp;
+  double hc, h0, h1, h2, beta_tmp, hfactor;
 
 	
@@ -1185,6 +1185,11 @@
       h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
       h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-      hmin = min(h0,min(h1,h2));
-      
+      hmin = min(min(h0,min(h1,h2)),hc);
+      //hfactor = hc/(hc + 1.0);
+
+      hfactor = 0.0;
+      if (hmin > 0.1 ) {
+	  hfactor = (hmin-0.1)/(hmin+0.4);
+      }
       
       if (optimise_dry_cells) {      
@@ -1230,9 +1235,13 @@
       
       // Playing with dry wet interface
-      hmin = qmin;
-      beta_tmp = beta_w;
-      if (hmin<minimum_allowed_height)
-	beta_tmp = beta_w_dry;
+      //hmin = qmin;
+      //beta_tmp = beta_w_dry;
+      //if (hmin>minimum_allowed_height)
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
 	
+      //printf("min_alled_height = %f\n",minimum_allowed_height);
+      //printf("hmin = %f\n",hmin);
+      //printf("beta_w = %f\n",beta_w);
+      //printf("beta_tmp = %f\n",beta_tmp);
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp); 
@@ -1271,8 +1280,9 @@
       // from the centroid to the min and max
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      beta_tmp = beta_uh;
-      if (hmin<minimum_allowed_height)
-	beta_tmp = beta_uh_dry;
-	
+      //beta_tmp = beta_uh;
+      //if (hmin<minimum_allowed_height)
+      //beta_tmp = beta_uh_dry;
+      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp);
@@ -1311,9 +1321,11 @@
       // from the centroid to the min and max
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      beta_tmp = beta_vh;
-      
-      if (hmin<minimum_allowed_height)
-	beta_tmp = beta_vh_dry;
-	
+      
+      //beta_tmp = beta_vh;
+      //
+      //if (hmin<minimum_allowed_height)
+      //beta_tmp = beta_vh_dry;
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;	
+
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp);
@@ -1563,13 +1575,13 @@
   
 
-  beta_w      = get_double(domain,"beta_w");
-  beta_w_dry  = get_double(domain,"beta_w_dry");
-  beta_uh     = get_double(domain,"beta_uh");
-  beta_uh_dry = get_double(domain,"beta_uh_dry");
-  beta_vh     = get_double(domain,"beta_vh");
-  beta_vh_dry = get_double(domain,"beta_vh_dry");  
-
-  minimum_allowed_height = get_double(domain,"minimum_allowed_height");
-  epsilon = get_double(domain,"epsilon");
+  beta_w                 = get_python_double(domain,"beta_w");
+  beta_w_dry             = get_python_double(domain,"beta_w_dry");
+  beta_uh                = get_python_double(domain,"beta_uh");
+  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
+  beta_vh                = get_python_double(domain,"beta_vh");
+  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");  
+
+  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
+  epsilon                = get_python_double(domain,"epsilon");
 
   number_of_elements = stage_centroid_values -> dimensions[0];  
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5161)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5162)
@@ -248,11 +248,17 @@
         
         range = fid.variables['xmomentum_range'][:]
-        assert allclose(range,[0,0.4695096]) or\
-               allclose(range,[0,0.47790655]) # Old slope limiters
+
+
+
+        assert allclose(range,[0,0.4695096]) or \
+               allclose(range,[0,0.47790655]) or\
+               allclose(range,[0,0.46941957]) 
         
         range = fid.variables['ymomentum_range'][:]
         #assert allclose(range,[0,0.02174380])
+        
         assert allclose(range,[0,0.02174439]) or\
-               allclose(range,[0,0.02283983]) # Old slope limiters 
+               allclose(range,[0,0.02283983]) or\
+               allclose(range,[0,0.0217342]) # Old slope limiters 
         
         fid.close()
@@ -324,8 +330,8 @@
 
         extrema = fid.variables['xmomentum.extrema'][:]
-        assert allclose(extrema,[-0.06062178, 0.47886313])
+        assert allclose(extrema,[-0.06062178, 0.47873023])
         
         extrema = fid.variables['ymomentum.extrema'][:]
-        assert allclose(extrema,[0.00, 0.06241221])        
+        assert allclose(extrema,[0.00, 0.0625786])        
 
         time_interval = fid.variables['extrema.time_interval'][:]
Index: /anuga_core/source/anuga/utilities/util_ext.h
===================================================================
--- /anuga_core/source/anuga/utilities/util_ext.h	(revision 5161)
+++ /anuga_core/source/anuga/utilities/util_ext.h	(revision 5162)
@@ -261,5 +261,8 @@
     
   B = (PyArrayObject*) PyObject_GetAttrString(O, name);
-  if (!B) return NULL;     
+  if (!B) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }     
   
   //Convert to consecutive array
@@ -269,9 +272,12 @@
   Py_DECREF(B); //FIXME: Is this really needed??
   
-  if (!A) return NULL;
+  if (!A) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }
   return A;
 }
 
-double get_double(PyObject *O, char *name) {
+double get_python_double(PyObject *O, char *name) {
   PyObject *TObject;
   double tmp;
@@ -281,5 +287,5 @@
   TObject = PyObject_GetAttrString(O, name);
   if (!TObject) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_double could not obtain double from object");
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
     return 0.0;
   }  
@@ -293,3 +299,17 @@
 
 
-
+PyObject *get_python_object(PyObject *O, char *name) {
+  PyObject *Oout;
+
+  Oout = PyObject_GetAttrString(O, name);
+  if (!Oout) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
+    return NULL;
+  }
+
+  return Oout;
+
+}
+
+
+
Index: /anuga_core/source/anuga/utilities/xml_tools.py
===================================================================
--- /anuga_core/source/anuga/utilities/xml_tools.py	(revision 5161)
+++ /anuga_core/source/anuga/utilities/xml_tools.py	(revision 5162)
@@ -239,11 +239,13 @@
         fid = xml
 
-    try:
-        dom = parse(fid)
-    except Exception, e:
-        # Throw filename into dom exception
-        msg = 'XML file "%s" could not be parsed: ' %fid.name
-        msg += str(e)
-        raise Exception, msg
+    dom = parse(fid)
+
+##     try:
+##         dom = parse(fid)
+##     except Exception, e:
+##         # Throw filename into dom exception
+##         msg = 'XML file "%s" could not be parsed: ' %fid.name
+##         msg += str(e)
+##         raise Exception, msg
 
     return dom2object(dom)
Index: /anuga_validation/convergence_study/animatesww2d_alt.py
===================================================================
--- /anuga_validation/convergence_study/animatesww2d_alt.py	(revision 5161)
+++ /anuga_validation/convergence_study/animatesww2d_alt.py	(revision 5162)
@@ -1,19 +1,19 @@
-def animatesww2d_alt(sww_filename, movie_filename, range):
-    """Plot cross section of model output
-    """
+## def animatesww2d_alt(sww_filename, movie_filename, range):
+##     """Plot cross section of model output
+##     """
 
-    # Read SWW file	
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(sww_filename, 'r')
+##     # Read SWW file	
+##     from Scientific.IO.NetCDF import NetCDFFile
+##     fid = NetCDFFile(sww_filename, 'r')
     
-    x = fid.variables['x']
-    y = fid.variables['y']
-    volumes = fid.variables['volumes']
+##     x = fid.variables['x']
+##     y = fid.variables['y']
+##     volumes = fid.variables['volumes']
     
-    elevation = fid.variables['elevation']
-    time = fid.variables['time']
-    stage = fid.variables['stage']
-    xmomentum = fid.variables['xmomentum']
-    ymomentum = fid.variables['ymomentum']
+##     elevation = fid.variables['elevation']
+##     time = fid.variables['time']
+##     stage = fid.variables['stage']
+##     xmomentum = fid.variables['xmomentum']
+##     ymomentum = fid.variables['ymomentum']
     
 def animatesww2d(swwfile):
@@ -34,5 +34,5 @@
     # interpolation points are for y = 0
     # number of increments in x
-    m = 100 
+    m = 1000 
     max_x = 100000.
     x = 0
Index: /anuga_validation/convergence_study/convergence_structured.py
===================================================================
--- /anuga_validation/convergence_study/convergence_structured.py	(revision 5161)
+++ /anuga_validation/convergence_study/convergence_structured.py	(revision 5162)
@@ -17,11 +17,12 @@
 from anuga.shallow_water import Transmissive_boundary
 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from anuga.geospatial_data.geospatial_data import *
+#from anuga.geospatial_data.geospatial_data import *
 from math import cos
+from Numeric import zeros, Float
 
 #------------------------------------------------------------------------------
 # Setup computational domain
 #------------------------------------------------------------------------------
-dx = 1000.
+dx = 200.
 dy = dx
 L = 100000.
@@ -33,16 +34,14 @@
 domain = Domain(points, vertices, boundary) 
 
-domain.set_timestepping_method('euler')
+domain.set_timestepping_method('rk2')
 domain.set_default_order(2)
 domain.set_name('myexample9')                
 domain.set_datadir('.')                     # Use current directory for output
 
-domain.beta_w      = 1.0
-domain.beta_w_dry  = 0.2
-domain.beta_uh     = 1.0
-domain.beta_uh_dry = 0.2
-domain.beta_vh     = 1.0
-domain.beta_vh_dry = 0.2
-domain.beta_h      = 1.0
+domain.set_all_limiters(0.9)
+
+print domain.beta_w
+domain.use_old_limiter = False
+domain.CFL = 1.0
 
 #------------------------------------------------------------------------------
@@ -50,7 +49,22 @@
 #------------------------------------------------------------------------------
 #domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('elevation',-100)
+domain.set_quantity('elevation',0.0)
 domain.set_quantity('friction', 0.00)
-domain.set_quantity('stage', 0.0)            
+
+h0 = 10.0
+h1 = 1.0
+
+def height(x,y):
+    z = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<=50000.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+
+
+domain.set_quantity('stage', height)
+#domain.set_quantity('stage', 0.0)            
 
 #-----------------------------------------------------------------------------
@@ -65,5 +79,7 @@
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
 ## Sine wave
-                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+#                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+## Single wave
+                   f=lambda t: [h0, 0.0, 0.0])
 ## Sawtooth?
 #                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
@@ -92,5 +108,5 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 20.0, finaltime = 10*40*60.):
+for t in domain.evolve(yieldstep = 50.0, finaltime = 10*40*60.):
     domain.write_time()
     vis.update()
Index: /anuga_validation/convergence_study/dam_break.py
===================================================================
--- /anuga_validation/convergence_study/dam_break.py	(revision 5162)
+++ /anuga_validation/convergence_study/dam_break.py	(revision 5162)
@@ -0,0 +1,127 @@
+"""Simple water flow example using ANUGA
+
+Water driven up a linear slope and time varying boundary,
+similar to a beach environment
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import Transmissive_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+
+from math import cos
+from Numeric import zeros, Float
+from time import localtime, strftime, gmtime
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'dam_break_'+time
+output_file = 'dam_break'
+
+copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 250.
+dy = dx
+L = 100000.
+W = dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+domain.use_old_limiter = True
+domain.CFL = 1.0
+
+domain.beta_w      = 1.7
+domain.beta_w_dry  = 0.0
+domain.beta_uh     = 1.7
+domain.beta_uh_dry = 0.0
+domain.beta_vh     = 1.7
+domain.beta_vh_dry = 0.0
+
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+
+h0 = 10.0
+h1 = 0.0
+
+def height(x,y):
+    z = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<=50000.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+
+
+domain.set_quantity('stage', height)
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bt, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+from anuga.visualiser import RealtimeVisualiser
+vis = RealtimeVisualiser(domain)
+vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
+vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 20.0, finaltime = 60*60.):
+    domain.write_time()
+    vis.update()
+    
+vis.evolveFinished()
+
Index: /anuga_validation/convergence_study/wave.py
===================================================================
--- /anuga_validation/convergence_study/wave.py	(revision 5162)
+++ /anuga_validation/convergence_study/wave.py	(revision 5162)
@@ -0,0 +1,129 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave which showed diffusive effects of
+thefirst order and standard second order method. Problem resolved if "rk2" timestepping
+and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import Transmissive_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+
+from math import cos
+from Numeric import zeros, Float
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'wave_'+time
+output_file = 'wave'
+
+copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 500.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+#domain.set_timestepping_method('euler')
+
+#domain.set_default_order(2)
+#domain.set_timestepping_method('euler')
+#domain.set_all_limiters(2.0)
+#domain.use_old_limiter = True
+#domain.CFL = 1.0
+
+
+domain.set_default_order(2)
+#domain.set_timestepping_method('euler')
+
+domain.use_old_limiter = True
+
+domain.beta_w      = 2.0
+domain.beta_w_dry  = 0.5
+domain.beta_uh     = 2.0
+domain.beta_uh_dry = 0.5
+domain.beta_vh     = 2.0
+domain.beta_vh_dry = 0.5
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',-100.0)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)            
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+## Sine wave
+                  f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+from anuga.visualiser import RealtimeVisualiser
+vis = RealtimeVisualiser(domain)
+vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
+    domain.write_time()
+    vis.update()
+    
+vis.evolveFinished()
+
Index: /anuga_work/development/analytical_solutions/MacDonald.xya
===================================================================
--- /anuga_work/development/analytical_solutions/MacDonald.xya	(revision 5161)
+++ /anuga_work/development/analytical_solutions/MacDonald.xya	(revision 5162)
Index: /anuga_work/development/analytical_solutions/MacDonald_bed.xya
===================================================================
--- /anuga_work/development/analytical_solutions/MacDonald_bed.xya	(revision 5161)
+++ /anuga_work/development/analytical_solutions/MacDonald_bed.xya	(revision 5162)
Index: /anuga_work/development/analytical_solutions/Malpasset_new.py
===================================================================
--- /anuga_work/development/analytical_solutions/Malpasset_new.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/Malpasset_new.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/circular_hydraulic_jump.py
===================================================================
--- /anuga_work/development/analytical_solutions/circular_hydraulic_jump.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/circular_hydraulic_jump.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/contracting_channel.py
===================================================================
--- /anuga_work/development/analytical_solutions/contracting_channel.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/contracting_channel.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/dam_break_non_symmetric.py
===================================================================
--- /anuga_work/development/analytical_solutions/dam_break_non_symmetric.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/dam_break_non_symmetric.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/fit_example.py
===================================================================
--- /anuga_work/development/analytical_solutions/fit_example.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/fit_example.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/oblique_shock.py
===================================================================
--- /anuga_work/development/analytical_solutions/oblique_shock.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/oblique_shock.py	(revision 5162)
@@ -29,10 +29,10 @@
 
 
-from shallow_water import Domain, Constant_height
-from shallow_water import Transmissive_boundary, Reflective_boundary,\
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
      Dirichlet_boundary
-
+from anuga.visualiser import RealtimeVisualiser
 from math import sqrt, cos, sin, pi
-from mesh_factory import oblique
+from mesh_factory import oblique_cross
 
 
@@ -44,8 +44,8 @@
 leny = 30.
 lenx = 40.
-n = 50
-m = 60
+n = 100
+m = 120
 
-points, elements, boundary = oblique(m, n, lenx, leny)
+points, elements, boundary = oblique_cross(m, n, lenx, leny)
 domain = Domain(points, elements, boundary)
 
@@ -54,8 +54,8 @@
 
 # Store output
-domain.store=True
+#domain.store=True
 
 # Output format
-domain.format="sww" #NET.CDF binary format
+#domain.format="sww" #NET.CDF binary format
                     # "dat" for ASCII
 
@@ -87,6 +87,14 @@
 #Initial condition
 h = 0.5
-domain.set_quantity('stage', Constant_height(x_slope, h) )
+domain.set_quantity('stage', expression = 'elevation + %d'%h )
 
+#---------------------------------
+# Setup visualization
+#---------------------------------
+vis = RealtimeVisualiser(domain)
+vis.render_quantity_height("elevation", dynamic=False)
+vis.render_quantity_height("stage", zScale=10, dynamic=True)
+vis.colour_height_quantity('stage', (0.75, 0.5, 0.5))
+vis.start()
 
 ######################
@@ -96,6 +104,11 @@
 for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
     domain.write_time()
+    vis.update()
 
 print 'That took %.2f seconds' %(time.time()-t0)
+
+
+vis.evolveFinished()
+vis.join()
 
 #FIXME: Compute average water depth on either side of shock and compare
Index: /anuga_work/development/analytical_solutions/oblique_shock_01.py
===================================================================
--- /anuga_work/development/analytical_solutions/oblique_shock_01.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/oblique_shock_01.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/oblique_shock_1.py
===================================================================
--- /anuga_work/development/analytical_solutions/oblique_shock_1.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/oblique_shock_1.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/sampson_cross_mesh.py
===================================================================
--- /anuga_work/development/analytical_solutions/sampson_cross_mesh.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/sampson_cross_mesh.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/wave_runup.py
===================================================================
--- /anuga_work/development/analytical_solutions/wave_runup.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/wave_runup.py	(revision 5162)
Index: /anuga_work/development/analytical_solutions/yoon_cross_mesh.py
===================================================================
--- /anuga_work/development/analytical_solutions/yoon_cross_mesh.py	(revision 5161)
+++ /anuga_work/development/analytical_solutions/yoon_cross_mesh.py	(revision 5162)
Index: /anuga_work/development/near_shore_PMD/run_beach.py
===================================================================
--- /anuga_work/development/near_shore_PMD/run_beach.py	(revision 5161)
+++ /anuga_work/development/near_shore_PMD/run_beach.py	(revision 5162)
@@ -167,4 +167,5 @@
     domain.set_boundary( {'wall': Br, 'wave': Bwp_velocity} )
 
+
     #-------------------------------------------------------------------------
     # Evolve system through time
@@ -174,4 +175,5 @@
     for t in domain.evolve(yieldstep, finaltime):    
         domain.write_time()
+
     print 'That took %.2f seconds' %(time.time()-t0)
     print 'finished'
