Index: inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/mesh.py	(revision 261)
+++ inundation/ga/storm_surge/pyvolution/mesh.py	(revision 262)
@@ -81,4 +81,6 @@
         #Allocate space for geometric quantities
         self.centroids = zeros((N, 2), Float)
+	#FIXME: Should be renamed to centroid_coordinates
+	
         self.areas = zeros(N, Float)
         self.radii = zeros(N, Float)
Index: inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/quantity.py	(revision 261)
+++ inundation/ga/storm_surge/pyvolution/quantity.py	(revision 262)
@@ -451,3 +451,4 @@
     #Replace python version with c implementations
     	
-    from quantity_ext import limit, compute_gradients #, extrapolate_second_order        
+    from quantity_ext import limit, compute_gradients,\
+    extrapolate_second_order        
Index: inundation/ga/storm_surge/pyvolution/quantity_ext.c
===================================================================
--- inundation/ga/storm_surge/pyvolution/quantity_ext.c	(revision 261)
+++ inundation/ga/storm_surge/pyvolution/quantity_ext.c	(revision 262)
@@ -94,14 +94,62 @@
 
 
-
-
+int _extrapolate(int N,
+		 double* centroids, 
+		 double* centroid_values, 		 
+		 double* vertex_coordinates,
+		 double* vertex_values,
+		 double* a,
+		 double* b) {
+  
+  int k, k2, k3, k6;
+  double x, y, x0, y0, x1, y1, x2, y2;
+		     
+  for (k=0; k<N; k++) {
+    k6 = 6*k;
+    k3 = 3*k;    
+    k2 = 2*k;
+        
+    //Centroid coordinates            
+    x = centroids[k2]; y = centroids[k2+1]; 
+
+    //vertex coordinates
+    //x0, y0, x1, y1, x2, y2 = X[k,:]
+    x0 = vertex_coordinates[k6 + 0];
+    y0 = vertex_coordinates[k6 + 1];    
+    x1 = vertex_coordinates[k6 + 2];
+    y1 = vertex_coordinates[k6 + 3];        
+    x2 = vertex_coordinates[k6 + 4];
+    y2 = vertex_coordinates[k6 + 5];           
+
+    //Extrapolate
+    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);
+    
+    /*
+    printf("V %f, %f, %f\n", vertex_values[k3+0], 
+	   vertex_values[k3+1], 
+	   vertex_values[k3+2]); 
+    */
+	   
+    //printf("k=%d, a=%f, b=%f\n", k, a[k], b[k]); 
+    
+	   
+    //printf("C %f, %f, %f\n", centroid_values[k], 
+    //	   centroid_values[k], 
+    //	   centroid_values[k]); 	   
+	   
+    //printf("X %f, %f, %f, %f, %f, %f\n", x0, y0, x1, y1, x2, y2);
+  }
+  return 0;
+}
+
+
+  
 /////////////////////////////////////////////////
 // Gateways to Python
-
-
-
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
   
-  PyObject *quantity, *domain;
+  PyObject *quantity, *domain, *R;
   PyArrayObject 
     *centroids,            //Coordinates at centroids
@@ -162,16 +210,20 @@
   }
   		     
-  
-//Return values
-  return Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 
+  //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_second_order(PyObject *self, PyObject *args) {
-
-
-
-PyObject *compute_gradients(PyObject *self, PyObject *args) {
   
   PyObject *quantity, *domain;
@@ -179,9 +231,13 @@
     *centroids,            //Coordinates at centroids
     *centroid_values,      //Values at centroids    
+    *vertex_coordinates,   //Coordinates at vertices            
+    *vertex_values,        //Values at vertices        
     *number_of_boundaries, //Number of boundaries for each triangle
     *surrogate_neighbours, //True neighbours or - if one missing - self
-    *a, *b;                //Return values
+    *a, *b;                //Gradients
    
-  int dimensions[1], N, err;
+  //int N, err;
+  int dimensions[1], N, err;  
+  //double *a, *b;  //Gradients
   
   // Convert Python arguments to C  
@@ -208,4 +264,12 @@
     PyObject_GetAttrString(domain, "number_of_boundaries");           
   if (!number_of_boundaries) return NULL;            
+  
+  vertex_coordinates = (PyArrayObject*) 
+    PyObject_GetAttrString(domain, "vertex_coordinates");           
+  if (!vertex_coordinates) return NULL;              
+  
+  vertex_values = (PyArrayObject*) 
+    PyObject_GetAttrString(quantity, "vertex_values");           
+  if (!vertex_values) return NULL;                
   
   N = centroid_values -> dimensions[0];
@@ -218,5 +282,11 @@
   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  
-        
+
+  //FIXME: Odd that I couldn't use normal arrays	  
+  //Allocate space for return vectors a and b (don't DECREF) 
+  //a = (double*) malloc(N * sizeof(double));
+  //if (!a) return NULL;  
+  //b = (double*) malloc(N * sizeof(double));  
+  //if (!b) return NULL;
 
   
@@ -227,5 +297,5 @@
 			   (int*) surrogate_neighbours -> data,
 			   (double*) a -> data,
-			   (double*) b -> data);     
+			   (double*) b -> data);
 			   
   if (err != 0) {
@@ -233,57 +303,37 @@
     return NULL;
   }
-
-
-  PyObject *quantity, *domain, *Tmp;
-  PyArrayObject 
-    *qv, //Conserved quantities at vertices 
-    *qc, //Conserved quantities at centroids
-    *neighbours;
-    
-  int k, i, n, N, k3;
-  double beta; //Safety factor
-  double *qmin, *qmax, qn;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "O", &quantity)) 
-    return NULL;
-
-  domain = PyObject_GetAttrString(quantity, "domain");     
-  if (!domain) 
-    return NULL;  
-    
-  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
-
-  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
-  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
-  N = qc -> dimensions[0];
-*/
-  /*
-def extrapolate_second_order(self):        
-    """Extrapolate conserved quantities from centroid to
-    vertices for each volume using
-    second order scheme.
-    """
-        
-    a, b = self.compute_gradients()
-
-    V = self.domain.get_vertex_coordinates()
-    qc = self.centroid_values
-    qv = self.vertex_values
-    
-    #Check each triangle
-    for k in range(self.domain.number_of_elements):
-        #Centroid coordinates            
-        x, y = self.domain.centroids[k]
-        
-        #vertex coordinates
-        x0, y0, x1, y1, x2, y2 = V[k,:]
-        
-        #Extrapolate
-        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
-        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
-  */
-//}    
+  
+  err = _extrapolate(N, 
+		     (double*) centroids -> data, 
+		     (double*) centroid_values -> data,		     
+		     (double*) vertex_coordinates -> data,
+		     (double*) vertex_values -> data,
+		     (double*) a -> data,
+		     (double*) b -> data);     		     
+  //a, b);
+ 		     
+		     
+  if (err != 0) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "Internal function _extrapolate failed");
+    return NULL;
+  }		     
+		
+  
+
+  //Release  
+  Py_DECREF(centroids);    
+  Py_DECREF(centroid_values);    
+  Py_DECREF(number_of_boundaries); 
+  Py_DECREF(surrogate_neighbours);
+  Py_DECREF(vertex_coordinates);                   
+  Py_DECREF(vertex_values);                 
+  Py_DECREF(a);                   
+  Py_DECREF(b);                     
+  //free(a);      
+  //free(b);        
+  
+  return Py_BuildValue("");    
+}    
 
 
@@ -359,4 +409,6 @@
   {"limit", limit, METH_VARARGS, "Print out"},    
   {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},      
+  {"extrapolate_second_order", extrapolate_second_order, 
+   METH_VARARGS, "Print out"},        
   {NULL, NULL, 0, NULL}   /* sentinel */
 };
Index: inundation/ga/storm_surge/pyvolution/test_quantity.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 261)
+++ inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 262)
@@ -194,8 +194,9 @@
         quantity.set_values([2.0, 4.0, 8.0, 2.0],
                             location = 'centroids')
-        
+
+
+
         #Gradients
         a, b = quantity.compute_gradients()
-
 
         #gradient bewteen t0 and t1 is undefined as det==0
@@ -209,5 +210,6 @@
 
         quantity.extrapolate_second_order()
-        
+
+	#print quantity.vertex_values
         assert allclose(quantity.vertex_values, [[2., 2.,  2.],
                                                  [0., 6.,  6.],
@@ -215,6 +217,5 @@
                                                  [0., 0.,  6.]])
 
-
-
+		
     def test_second_order_extrapolation2(self):
         quantity = Conserved_quantity(self.mesh4)        
