Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4885)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4886)
@@ -229,4 +229,9 @@
         #(either from this module or C-extension)
         interpolate_from_vertices_to_edges(self)
+
+    def interpolate_from_edges_to_vertices(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        interpolate_from_edges_to_vertices(self)
 
 
@@ -1400,14 +1405,14 @@
 
 
-##     def limit_by_vertex(self):
-##         #Call correct module function
-##         #(either from this module or C-extension)
-##         limit_by_vertex()
-
-
-##     def limit_by_edge(self):
-##         #Call correct module function
-##         #(either from this module or C-extension)
-##         limit_by_edge()        
+    def limit_by_vertex(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        limit_by_vertex(self)
+
+
+    def limit_by_edge(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        limit_by_edge(self)        
 
     def extrapolate_second_order(self):
@@ -1439,6 +1444,9 @@
          compute_gradients,\
          limit_old,\
+         limit_by_vertex,\
+         limit_by_edge,\
          extrapolate_second_order,\
          interpolate_from_vertices_to_edges,\
+         interpolate_from_edges_to_vertices,\
          update    
 else:
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4885)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4886)
@@ -122,10 +122,10 @@
 		y2 = vertex_coordinates[k6 + 5];
 
-		// Extrapolate
+		// Extrapolate to Vertices
 		vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
 		vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
 		vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
 
-
+		// Extrapolate to Edges (midpoints)
 		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]);
@@ -137,7 +137,123 @@
 
 
-
-
-int _interpolate(int N,
+int _limit_by_vertex(int N, double beta,
+		     double* centroid_values,
+		     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];
+		qmin = qc;
+		qmax = qc;
+
+		for (i=0; i<3; i++) {
+		    n = neighbours[k3+i];
+		    if (n >= 0) {
+			qn = centroid_values[n]; //Neighbour's centroid value
+
+			qmin = min(qmin, qn);
+			qmax = max(qmax, qn);
+		    }
+		}
+
+		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
+      
+		    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 vertex and edge values using phi limiter
+		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 _limit_by_edge(int N, double beta,
+		     double* centroid_values,
+		     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];
+		qmin = qc;
+		qmax = qc;
+
+		for (i=0; i<3; i++) {
+		    n = neighbours[k3+i];
+		    if (n >= 0) {
+			qn = centroid_values[n]; //Neighbour's centroid value
+
+			qmin = min(qmin, qn);
+			qmax = max(qmax, qn);
+		    }
+		}
+
+		phi = 1.0;
+		for (i=0; i<3; i++) {    
+		    r = 1.0;
+      
+		    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+		    dqa[i] = dq;                      //Save dq for use in updating vertex values
+      
+		    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 vertex and edge 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 _interpolate_from_vertices_to_edges(int N,
 		 double* vertex_values,
 		 double* edge_values) {
@@ -157,4 +273,27 @@
 		edge_values[k3 + 1] = 0.5*(q0+q2);
 		edge_values[k3 + 2] = 0.5*(q0+q1);
+	}
+	return 0;
+}
+
+
+int _interpolate_from_edges_to_vertices(int N,
+		 double* vertex_values,
+		 double* edge_values) {
+
+	int k, k3;
+	double e0, e1, e2;
+
+
+	for (k=0; k<N; k++) {
+		k3 = 3*k;
+
+		e0 = edge_values[k3 + 0];
+		e1 = edge_values[k3 + 1];
+		e2 = edge_values[k3 + 2];
+
+		vertex_values[k3 + 0] = e1 + e2 - e0;
+		vertex_values[k3 + 1] = e2 + e0 - e1;
+		vertex_values[k3 + 2] = e0 + e1 - e2;
 	}
 	return 0;
@@ -458,5 +597,45 @@
 	N = vertex_values -> dimensions[0];
 
-	err = _interpolate(N,
+	err = _interpolate_from_vertices_to_edges(N,
+			   (double*) vertex_values -> data,
+			   (double*) edge_values -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "Interpolate could not be computed");
+	  return NULL;
+	}
+
+	// Release and return
+	Py_DECREF(vertex_values);
+	Py_DECREF(edge_values);
+
+	return Py_BuildValue("");
+}
+
+
+PyObject *interpolate_from_edges_to_vertices(PyObject *self, PyObject *args) {
+        //
+        //Compute vertex values from edge values using linear interpolation
+        // 
+	
+	PyObject *quantity;
+	PyArrayObject *vertex_values, *edge_values;
+
+	int N, err;
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "O", &quantity)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: interpolate_from_edges_to_vertices could not parse input");
+	  return NULL;
+	}
+	
+	vertex_values = get_consecutive_array(quantity, "vertex_values");
+	edge_values = get_consecutive_array(quantity, "edge_values");
+
+	N = vertex_values -> dimensions[0];
+
+	err = _interpolate_from_edges_to_vertices(N,
 			   (double*) vertex_values -> data,
 			   (double*) edge_values -> data);
@@ -795,8 +974,179 @@
 
 
+PyObject *limit_by_vertex(PyObject *self, PyObject *args) {
+  //Limit slopes for each volume to eliminate artificial variance
+  //introduced by e.g. second order extrapolator
+
+  //This is an unsophisticated limiter as it does not take into
+  //account dependencies among quantities.
+
+  //precondition:
+  //  vertex values are estimated from gradient
+  //postcondition:
+  //  vertex and edge values are updated
+  //
+
+	PyObject *quantity, *domain, *Tmp;
+	PyArrayObject
+	    *vertex_values,   //Conserved quantities at vertices
+	    *centroid_values, //Conserved quantities at centroids
+	    *edge_values,     //Conserved quantities at edges
+	    *neighbours;
+
+	double beta_w; //Safety factor
+	int N, err;
+
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "O", &quantity)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: limit_by_vertex could not parse input");
+	  return NULL;
+	}
+
+	domain = PyObject_GetAttrString(quantity, "domain");
+	if (!domain) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: limit_by_vertex 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_by_vertex could not obtain beta_w object from domain");		  	
+	  
+	  return NULL;
+	}	
+
+
+	// Get pertinent variables
+	neighbours       = get_consecutive_array(domain, "neighbours");
+	centroid_values  = get_consecutive_array(quantity, "centroid_values");
+	vertex_values    = get_consecutive_array(quantity, "vertex_values");
+	edge_values      = get_consecutive_array(quantity, "edge_values");
+	beta_w           = PyFloat_AsDouble(Tmp);
+
+
+	N = centroid_values -> dimensions[0];
+
+	err = _limit_by_vertex(N, beta_w,
+			       (double*) centroid_values -> data,
+			       (double*) vertex_values -> data,
+			       (double*) edge_values -> data,
+			       (long*)   neighbours -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError,
+			  "Internal function _limit_by_vertex failed");
+	  return NULL;
+	}	
+
+
+	// Release
+	Py_DECREF(neighbours);
+	Py_DECREF(centroid_values);
+	Py_DECREF(vertex_values);
+	Py_DECREF(edge_values);
+	Py_DECREF(Tmp);
+
+
+	return Py_BuildValue("");
+}
+
+
+
+PyObject *limit_by_edge(PyObject *self, PyObject *args) {
+  //Limit slopes for each volume to eliminate artificial variance
+  //introduced by e.g. second order extrapolator
+
+  //This is an unsophisticated limiter as it does not take into
+  //account dependencies among quantities.
+
+  //precondition:
+  //  vertex values are estimated from gradient
+  //postcondition:
+  //  vertex and edge values are updated
+  //
+
+	PyObject *quantity, *domain, *Tmp;
+	PyArrayObject
+	    *vertex_values,   //Conserved quantities at vertices
+	    *centroid_values, //Conserved quantities at centroids
+	    *edge_values,     //Conserved quantities at edges
+	    *neighbours;
+
+	double beta_w; //Safety factor
+	int N, err;
+
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "O", &quantity)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: limit_by_edge could not parse input");
+	  return NULL;
+	}
+
+	domain = PyObject_GetAttrString(quantity, "domain");
+	if (!domain) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: limit_by_edge 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_by_edge could not obtain beta_w object from domain");		  	
+	  
+	  return NULL;
+	}	
+
+
+	// Get pertinent variables
+	neighbours       = get_consecutive_array(domain, "neighbours");
+	centroid_values  = get_consecutive_array(quantity, "centroid_values");
+	vertex_values    = get_consecutive_array(quantity, "vertex_values");
+	edge_values      = get_consecutive_array(quantity, "edge_values");
+	beta_w           = PyFloat_AsDouble(Tmp);
+
+
+	N = centroid_values -> dimensions[0];
+
+	err = _limit_by_edge(N, beta_w,
+			       (double*) centroid_values -> data,
+			       (double*) vertex_values -> data,
+			       (double*) edge_values -> data,
+			       (long*)   neighbours -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError,
+			  "Internal function _limit_by_vertex failed");
+	  return NULL;
+	}	
+
+
+	// Release
+	Py_DECREF(neighbours);
+	Py_DECREF(centroid_values);
+	Py_DECREF(vertex_values);
+	Py_DECREF(edge_values);
+	Py_DECREF(Tmp);
+
+
+	return Py_BuildValue("");
+}
+
+
 
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
 	{"limit_old", limit_old, METH_VARARGS, "Print out"},
+	{"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"},
+	{"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"},
 	{"update", update, METH_VARARGS, "Print out"},
 	{"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
@@ -808,4 +1158,7 @@
 		interpolate_from_vertices_to_edges,
 		METH_VARARGS, "Print out"},
+	{"interpolate_from_edges_to_vertices",
+		interpolate_from_edges_to_vertices,
+		METH_VARARGS, "Print out"},
 	{"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},		
 	{NULL, NULL, 0, NULL}   // sentinel
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 4885)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 4886)
@@ -1445,5 +1445,5 @@
 
 
-    def test_limiter(self):
+    def test_vertex_limiter(self):
         quantity = Conserved_quantity(self.mesh4)
 
@@ -1453,5 +1453,5 @@
 
         #Limit
-        quantity.limit()
+        quantity.limit_by_vertex()
 
         #Assert that central triangle is limited by neighbours
@@ -1474,4 +1474,34 @@
 
 
+
+    def test_edge_limiter(self):
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+
+
+        #Limit
+        quantity.limit_by_edge()
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
+
+        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
+
+        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,2] <= quantity.centroid_values[0]
+
+
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/3)
+
+
     def test_limiter2(self):
         """Taken from test_shallow_water
@@ -1518,9 +1548,9 @@
 
 
-        #Extrapolate
+        #Extrapolate from centroid to vertices and edges
         quantity.extrapolate_first_order()
 
         #Interpolate
-        quantity.interpolate_from_vertices_to_edges()
+        #quantity.interpolate_from_vertices_to_edges()
 
         assert allclose(quantity.vertex_values,
@@ -1528,4 +1558,32 @@
         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
                                                [3,3,3], [4, 4, 4]])
+
+
+    def test_interpolate_from_vertices_to_edges(self):
+        quantity = Quantity(self.mesh4)
+
+        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
+
+        quantity.interpolate_from_vertices_to_edges()
+
+        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
+                                               [3., 2.5, 1.5],
+                                               [3.5, 4.5, 3.],
+                                               [2.5, 3.5, 2]])
+
+
+    def test_interpolate_from_edges_to_vertices(self):
+        quantity = Quantity(self.mesh4)
+
+        quantity.edge_values = array([[1., 1.5, 0.5],
+                                [3., 2.5, 1.5],
+                                [3.5, 4.5, 3.],
+                                [2.5, 3.5, 2]],Float)
+
+        quantity.interpolate_from_edges_to_vertices()
+
+        assert allclose(quantity.vertex_values,
+                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+
 
 
Index: /anuga_core/source/anuga/utilities/util_ext.h
===================================================================
--- /anuga_core/source/anuga/utilities/util_ext.h	(revision 4885)
+++ /anuga_core/source/anuga/utilities/util_ext.h	(revision 4886)
@@ -166,69 +166,4 @@
 
 
-void _limit_by_edge(int N, double beta, double* qc, double* qv, 
-	    double* qmin, double* qmax) { 
-
-  //N are the number of elements
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  //printf("INSIDE\n");
-  for (k=0; k<N; k++) {
-    k3 = k*3;
-    
-    //Find the gradient limiter (phi) across vertices  
-    phi = 1.0;
-    for (i=0; i<3; i++) {    
-      r = 1.0;
-      
-      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
-      dqa[i] = dq;              //Save dq for use in the next loop
-      
-      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
-      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
-  
-  
-      phi = min( min(r*beta, 1.0), phi);    
-    }
-    
-    //Then update using phi limiter
-    for (i=0; i<3; i++) {    
-      qv[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}
-
-void _limit_by_vertex(int N, double beta, double* qc, double* qv, 
-	    double* qmin, double* qmax) { 
-
-  //N are the number of elements
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  //printf("INSIDE\n");
-  for (k=0; k<N; k++) {
-    k3 = k*3;
-    
-    //Find the gradient limiter (phi) across vertices  
-    phi = 1.0;
-    for (i=0; i<3; i++) {    
-      r = 1.0;
-      
-      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
-      dqa[i] = dq;              //Save dq for use in the next loop
-      
-      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
-      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
-  
-  
-      phi = min( min(r*beta, 1.0), phi);    
-    }
-    
-    //Then update using phi limiter
-    for (i=0; i<3; i++) {    
-      qv[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}
 
 
