Index: /anuga_core/documentation/user_manual/anuga_user_manual.tex
===================================================================
--- /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 4470)
+++ /anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 4471)
@@ -1773,8 +1773,8 @@
 
 
-\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag, reduction=None}
+\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag}
 Decide whether vertex values should be stored uniquely as
 computed in the model or whether they should be reduced to one
-value per vertex using self.reduction.
+value per vertex using averaging.
 \end{methoddesc}
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 4471)
@@ -394,4 +394,7 @@
         """
 
+        # FIXME (Ole): Refactor this based on algorithm in test and get
+        # rid of the old vertexlist
+        
         vertexlist = [None]*self.number_of_nodes
         for i in range(self.number_of_triangles):
@@ -410,5 +413,39 @@
                 vertexlist[v].append( (i, vertex_id) )
 
+
+        # Register number of triangles touching each nodes
+        number_of_triangles_per_node = zeros(self.number_of_nodes)        
+        number_of_entries = 0        
+        for i, vertices in enumerate(vertexlist):
+            try:
+                v = len(vertices)
+            except:
+                v = 0 # Lone node - e.g. not part of the mesh
+                
+            number_of_triangles_per_node[i] = v
+            number_of_entries += v
+                
+
+        # Build vertex value index array
+        vertex_value_indices = zeros(number_of_entries)
+        k = 0
+        for vertices in vertexlist:
+            if vertices is not None:
+                for i, v_id in vertices:
+                    vertex_value_indices[k] = 3*i + v_id
+
+                    k += 1
+
+        assert k == number_of_entries    
         self.vertexlist = vertexlist
+
+        self.number_of_triangles_per_node = number_of_triangles_per_node
+        self.vertex_value_indices = vertex_value_indices
+        self.vertexlist = vertexlist
+
+        #print
+        #print number_of_triangles_per_node
+        #print vertex_value_indices
+        #print vertexlist
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 4471)
@@ -701,4 +701,28 @@
         #
         #See domain.set_boundary
+
+        # Check integrity of vertex_value_indices
+        current_node = 0
+        k = 0 # Track triangles touching on node
+        for index in self.vertex_value_indices:
+
+            if self.number_of_triangles_per_node[current_node] == 0:
+                # Node is lone - i.e. not part of the mesh
+                continue
+            
+            k += 1
+            
+            volume_id = index / 3
+            vertex_id = index % 3
+            
+            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
+                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
+            assert self.triangles[volume_id, vertex_id] == current_node, msg
+                        
+            if self.number_of_triangles_per_node[current_node] == k:
+                # Move on to next node
+                k = 0
+                current_node += 1
+
 
     def get_lone_vertices(self):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4471)
@@ -943,5 +943,4 @@
                 
                 # FIXME (Ole): Should we merge this with get_vertex_values
-                # and use the concept of a reduction operator?
                 sum = 0
                 for triangle_id, vertex_id in triangles:
@@ -1033,6 +1032,5 @@
                           xy=True,
                           smooth=None,
-                          precision=None,
-                          reduction=None):
+                          precision=None):
         """Return vertex values like an OBJ format i.e. one value per node.
 
@@ -1046,6 +1044,6 @@
 
         if smooth is True, vertex values corresponding to one common
-        coordinate set will be smoothed according to the given
-        reduction operator. In this case vertex coordinates will be
+        coordinate set will be smoothed by taking the average of vertex values for each node.
+        In this case vertex coordinates will be
         de-duplicated corresponding to the original nodes as obtained from
         the method general_mesh.get_nodes()
@@ -1078,29 +1076,48 @@
 
         if smooth is True:
-            # Ensure continuous vertex values by combining
-            # values at each node using the reduction operator
+            # Ensure continuous vertex values by averaging
+            # values at each node
             
-            if reduction is None:
-                # Take default from domain                
-                reduction = self.domain.reduction
-
             V = self.domain.get_triangles()
             N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
-            A = zeros(N, precision)
+            A = zeros(N, Float)
             points = self.domain.get_nodes()            
-
-            # Reduction loop
-            for k in range(N):
-                L = self.domain.vertexlist[k]
-
-                # Go through all triangle, vertex pairs
-                # contributing to vertex k and register vertex value
-                if L is None: continue # In case there are unused points
-                contributions = []
-                for volume_id, vertex_id in L:
+            
+            if 1:
+                # Fast C version
+                average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
+                                      ensure_numeric(self.domain.number_of_triangles_per_node),
+                                      ensure_numeric(self.vertex_values),
+                                      A)
+                A = A.astype(precision)
+            else:    
+
+                # Slow Python version
+                
+                current_node = 0
+                k = 0 # Track triangles touching on node
+                total = 0.0
+                for index in self.domain.vertex_value_indices:
+                    k += 1
+                    
+                    volume_id = index / 3
+                    vertex_id = index % 3
+                 
+                    #assert V[volume_id, vertex_id] == current_node
+                
                     v = self.vertex_values[volume_id, vertex_id]
-                    contributions.append(v)
-
-                A[k] = reduction(contributions)
+                    total += v
+
+                    #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
+                    if self.domain.number_of_triangles_per_node[current_node] == k:
+                        A[current_node] = total/k
+                
+                    
+                        # Move on to next node
+                        total = 0.0
+                        k = 0
+                        current_node += 1
+
+
 
         else:
@@ -1450,3 +1467,4 @@
 
     from quantity_ext import compute_gradients, limit,\
-    extrapolate_second_order, interpolate_from_vertices_to_edges, update
+    extrapolate_second_order, interpolate_from_vertices_to_edges, update,\
+    average_vertex_values
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4471)
@@ -218,4 +218,40 @@
 
 
+int _average_vertex_values(int N,
+			   long* vertex_value_indices,
+			   long* number_of_triangles_per_node,
+			   double* vertex_values, 
+			   double* A) {
+  //Average vertex values to obtain one value per node
+
+  int i, index; 
+  int k = 0; //Track triangles touching each node
+  int current_node = 0;
+  double total = 0.0;
+
+  for (i=0; i<N; i++) {
+    index = vertex_value_indices[i];
+    k += 1;
+            
+    //volume_id = index / 3
+    //vertex_id = index % 3
+    //total += self.vertex_values[volume_id, vertex_id]
+    total += vertex_values[index];
+      
+    //printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
+    if (number_of_triangles_per_node[current_node] == k) {
+      A[current_node] = total/k;
+                
+      // Move on to next node
+      total = 0.0;
+      k = 0;
+      current_node += 1;
+    }
+  }
+			   
+  return 0;
+}
+
+
 /////////////////////////////////////////////////
 // Gateways to Python
@@ -283,6 +319,6 @@
 
 	err = _interpolate(N,
-		     (double*) vertex_values -> data,
-		     (double*) edge_values -> data);
+			   (double*) vertex_values -> data,
+			   (double*) edge_values -> data);
 
 	if (err != 0) {
@@ -298,4 +334,45 @@
 	return Py_BuildValue("");
 }
+
+
+PyObject *average_vertex_values(PyObject *self, PyObject *args) {
+
+	PyArrayObject 
+	  *vertex_value_indices, 
+	  *number_of_triangles_per_node,
+	  *vertex_values,
+	  *A;
+	
+
+	int N, err;
+
+	// Convert Python arguments to C
+	if (!PyArg_ParseTuple(args, "OOOO",
+			      &vertex_value_indices, 
+			      &number_of_triangles_per_node,
+			      &vertex_values, 
+			      &A)) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "quantity_ext.c: average_vertex_values could not parse input");
+	  return NULL;
+	}
+	
+	N = vertex_value_indices -> dimensions[0];
+	// printf("Got parameters, N=%d\n", N);
+	err = _average_vertex_values(N,
+				     (long*) vertex_value_indices -> data,
+				     (long*) number_of_triangles_per_node -> data,
+				     (double*) vertex_values -> data, 
+				     (double*) A -> data);
+
+	if (err != 0) {
+	  PyErr_SetString(PyExc_RuntimeError, 
+			  "average_vertex_values could not be computed");
+	  return NULL;
+	}
+
+	return Py_BuildValue("");
+}
+
 
 
@@ -567,4 +644,5 @@
 		interpolate_from_vertices_to_edges,
 		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_general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 4471)
@@ -61,4 +61,52 @@
                         [triangles[0], triangles[4]])
         
+
+    def test_vertex_value_indices(self):
+        """Check that structures are correct.
+        """
+        from mesh_factory import rectangular
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        nodes, triangles, _ = rectangular(3, 6)
+        domain = General_mesh(nodes, triangles)
+
+        assert sum(domain.number_of_triangles_per_node) ==\
+               len(domain.vertex_value_indices)
+
+        # Check number of triangles per node
+        count = [0]*domain.number_of_nodes
+        for triangle in triangles:
+            for i in triangle:
+                count[i] += 1
+
+        assert allclose(count, domain.number_of_triangles_per_node)
+
+        #print nodes
+        #print triangles
+        #print domain.number_of_triangles_per_node
+        #print domain.vertex_value_indices        
+
+
+        # Check indices
+        current_node = 0
+        k = 0 # Track triangles touching on node
+        for index in domain.vertex_value_indices:
+            k += 1
+            
+            triangle = index / 3
+            vertex = index % 3
+
+            assert triangles[triangle, vertex] == current_node
+
+            if domain.number_of_triangles_per_node[current_node] == k:
+                # Move on to next node
+                k = 0
+                current_node += 1
+                
+
+
+
+
     def test_areas(self):
         from mesh_factory import rectangular
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py	(revision 4470)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py	(revision 4471)
@@ -1331,5 +1331,5 @@
     suite = unittest.makeSuite(Test_Util,'test')
 #    suite = unittest.makeSuite(Test_Util,'test_store_parameters')
-    runner = unittest.TextTestRunner()
+    runner = unittest.TextTestRunner(verbosity=0)
     runner.run(suite)
 
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 4470)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 4471)
@@ -393,9 +393,11 @@
         os.remove(sww.filename)
 
-    def test_sww_variable3(self):
+    def no_test_sww_variable3(self):
         """Test that sww information can be written correctly
         multiple timesteps using a different reduction operator (min)
         """
 
+        # Different reduction in sww files has been made obsolete.
+        
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
