Index: /trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py	(revision 8134)
+++ /trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py	(revision 8135)
@@ -1,8 +1,9 @@
 
 from anuga import Domain
+from anuga import Quantity
 from anuga.utilities.sparse import Sparse, Sparse_CSR
 from anuga.utilities.cg_solve import conjugate_gradient
 import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+from anuga import Dirichlet_boundary
 import numpy as num
 import kinematic_viscosity_ext
@@ -11,8 +12,29 @@
 class Kinematic_Viscosity:
     """
-    Class for setting up structures and matrices for kinematic viscosity
+    Class for setting up structures and matrices for kinematic viscosity differential
+    operator using centroid values.
+
+    div ( diffusivity grad )
+
+    where diffusvity is scalar quantity (defaults to quantity with values = 1)
+    boundary values of f are used to setup entries associated with cells with boundaries
+
+    There are procedures to apply this operator, ie
+
+    (1) Calculate div( diffusivity grad u )
+    using boundary values stored in u
+
+    (2) Calculate ( u + dt div( diffusivity grad u )
+    using boundary values stored in u
+
+    (3) Solve div( diffusivity grad u ) = f
+    for quantity f and using boundary values stored in u
+
+    (4) Solve ( u + dt div( diffusivity grad u ) = f
+    for quantity f using boundary values stored in u
+
     """
 
-    def __init__(self, domain, triangle_areas=True, verbose=False):
+    def __init__(self, domain, diffusivity = None, triangle_areas=True, verbose=False):
         if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
         #Expose the domain attributes
@@ -21,25 +43,50 @@
         self.mesh = domain.mesh
         self.boundary = domain.boundary
+
+        # Setup diffusivity quantity
+        if diffusivity is None:
+            diffusivity = Quantity(domain)
+            diffusivity.set_values(1.0)
+        #check that diffusivity is a quantity associated with domain
+        assert diffusivity.domain == domain
+
+        self.diffusivity = diffusivity
+        self.diffusivity_bdry_data = self.diffusivity.boundary_values
+        self.diffusivity_data = self.diffusivity.centroid_values
+
+
         self.n = len(self.domain)
         self.dt = 1e-6 #Need to set to domain.timestep
         self.boundary_len = len(domain.boundary)
         self.tot_len = self.n + self.boundary_len
+
+        # FIXME SR: maybe this should already exist in mesh!
         self.boundary_enum = self.enumerate_boundary()
+        
         self.verbose = verbose
 
         #Geometric Information
         if verbose: log.critical('Kinematic Viscosity: Building geometric structure')
+
         self.geo_structure_indices = num.zeros((self.n, 3), num.int)
         self.geo_structure_values = num.zeros((self.n, 3), num.float)
+
+        # Only needs to built once, doesn't change
         kinematic_viscosity_ext.build_geo_structure(self, self.n, self.tot_len)
+
+        # Setup type of scaling
         self.apply_triangle_areas = triangle_areas
-        self.triangle_areas = Sparse(self.n, self.n)
+
+        # FIXME SR: should this really be a matrix?
+        temp  = Sparse(self.n, self.n)
         for i in range(self.n):
-            self.triangle_areas[i, i] = 1.0 / self.mesh.areas[i]
+            temp[i, i] = 1.0 / self.mesh.areas[i]
             
-        self.triangle_areas = Sparse_CSR(self.triangle_areas)
-
-        #Application Information
+        self.triangle_areas = Sparse_CSR(temp)
+
+        # FIXME SR: More to do with solving equation
         self.qty_considered = 1 #1 or 2 (uh or vh respectively)
+
+        #FIXME SR: Do we need Sparse or should we just go with Sparse_CSR
         self.operator_matrix = Sparse(self.n, self.tot_len)
 
@@ -49,8 +96,11 @@
         self.operator_colind = num.zeros((4 * self.n, ), num.int)
         #Sparse_CSR.rowptr (4 entries in every row, we know this already) = [0,4,8,...,4*n]
-        self.operator_rowptr = 4 * num.indices((self.n + 1, ))[0, :]
-        self.stage_heights = num.zeros((self.n, 1), num.float)
-        self.stage_heights_scaling = Sparse(self.n, self.n)
-        self.boundary_vector = num.zeros((self.n, ), num.float)
+        self.operator_rowptr = 4 * num.arange(self.n + 1)
+
+        # Build matrix self.elliptic_operator_matrix [A B]
+        self.build_elliptic_operator()
+
+        # Build self.boundary_term
+        self.build_operator_boundary_term()
 
         self.parabolic_solve = False #Are we doing a parabolic solve at the moment?
@@ -61,4 +111,5 @@
     def enumerate_boundary(self):
         #Enumerate by boundary index
+        # FIXME SR: Probably should be part of mesh
         enumeration = {}
         for i in range(self.n):
@@ -70,4 +121,5 @@
 
     def set_qty_considered(self, qty):
+        # FIXME SR: Probably should just be set by quantity to which operation is applied
         if qty == 1 or qty == 'u':
             self.qty_considered = 1
@@ -78,46 +130,68 @@
             assert 0 == 1, msg
 
-    def apply_stage_heights(self, h):
-        msg = "(KV_Operator.apply_stage_heights) h vector has incorrect length"
-        assert h.size == self.n, msg
-
-        self.operator_matrix = Sparse(self.n, self.tot_len)
-
-        #Evaluate the boundary stage heights - use generic_domain.update_boundary? (check enumeration)
-        boundary_heights = num.zeros((self.boundary_len, 1), num.float)
-        for key in self.boundary_enum.keys():
-            boundary_heights[self.boundary_enum[key]] = self.boundary[key].evaluate()[0]
-
-        kinematic_viscosity_ext.build_operator_matrix(self, self.n, self.tot_len, h, boundary_heights)
-        self.operator_matrix = \
-            Sparse_CSR(None, self.operator_data, self.operator_colind, self.operator_rowptr, self.n, self.tot_len)
-
-        self.stage_heights = h
-        #Set up the scaling matrix
-        data = h
-        num.putmask(data, data != 0, 1 / data) #take the reciprocal of each entry unless it is zero
-        self.stage_heights_scaling = \
-            Sparse_CSR(None, data, num.arange(self.n), num.append(num.arange(self.n), self.n), self.n, self.n)
-
-        self.build_boundary_vector()
-
-        return self.operator_matrix
-
-    def build_boundary_vector(self):
-        #Require stage heights applied
-        X = num.zeros((self.tot_len, 2), num.float)
-        for key in self.boundary_enum.keys():
-            quantities = self.boundary[key].evaluate()
-            h_i = quantities[0]
-            if h_i == 0.0:
-                X[self.n + self.boundary_enum[key], :] = num.array([0.0, 0.0])
-            else:
-                X[self.n + self.boundary_enum[key], :] = quantities[1:] / h_i
-        self.boundary_vector = self.operator_matrix * X
+    def build_elliptic_operator(self):
+        """
+        Builds matrix representing
+
+        div ( diffusivity grad )
+
+        which has the form [ A B ]
+        """
+
+        #Arrays self.operator_data, self.operator_colind, self.operator_rowptr
+        # are setup via this call
+        kinematic_viscosity_ext.build_operator_matrix(self,self.diffusivity_data, self.diffusivity_bdry_data)
+
+        self.elliptic_operator_matrix = Sparse_CSR(None, \
+                self.operator_data, self.operator_colind, self.operator_rowptr, \
+                self.n, self.tot_len)
+
+#        #Set up the scaling matrix
+#        data = h
+#        num.putmask(data, data != 0, 1 / data) #take the reciprocal of each entry unless it is zero
+#        self.stage_heights_scaling = \
+#            Sparse_CSR(None, self.diffusivity_data, num.arange(self.n), num.arange(self.n +1), self.n, self.n)
+
+    def update_elliptic_operator(self):
+        """
+        Updates the data values of matrix representing
+
+        div ( diffusivity grad )
+
+        (as when diffusivitiy has changed)
+        """
+
+        #Array self.operator_data is changed by this call, which should flow
+        # through to the Sparse_CSR matrix.
+
+        kinematic_viscosity_ext.update_operator_matrix(self, self.n, self.tot_len, \
+          self.diffusivity_data, self.diffusivity_bdry_data)
+
+
+
+
+    def build_operator_boundary_term(self):
+        """
+        Operator has form [A B] and U = [ u ; b]
+
+        This procedure calculates B b which can be calculated as
+
+        [A B] [ 0 ; b]
+        """
+
+        n = self.n
+        tot_len = self.tot_len
+
+        self.boundary_term = num.zeros((n, ), num.float)
+        
+        X = num.zeros((tot_len,), num.float)
+
+        X[n:] = self.diffusivity_bdry_data
+        self.boundary_term[:] = self.operator_matrix * X
+
         #Tidy up
         if self.apply_triangle_areas:
-            self.boundary_vector = self.triangle_areas * self.boundary_vector
-
-        return self.boundary_vector
+            self.boundary_term[:] = self.triangle_areas * self.boundary_term
+
 
     def elliptic_multiply(self, V, qty_considered=None, include_boundary=True):
Index: /trunk/anuga_core/source/anuga/operators/kinematic_viscosity_ext.c
===================================================================
--- /trunk/anuga_core/source/anuga/operators/kinematic_viscosity_ext.c	(revision 8134)
+++ /trunk/anuga_core/source/anuga/operators/kinematic_viscosity_ext.c	(revision 8135)
@@ -85,11 +85,12 @@
 }
 
-int build_operator_matrix(int n, 
-                          int tot_len, 
-                          long *geo_indices, 
-                          double *geo_values, 
-                          double *h, 
-                          double *boundary_heights, 
-                          double *data, 
+
+int build_operator_matrix(int n,
+                          int tot_len,
+                          long *geo_indices,
+                          double *geo_values,
+                          double *cell_data,
+                          double *bdry_data,
+                          double *data,
                           long *colind) {
 	int i, k, edge, j[4], *sorted_j, this_index;
@@ -102,10 +103,10 @@
 			j[edge] = geo_indices[3*i+edge];
 			if (j[edge]<n) { //interior
-				h_j = h[j[edge]];
+				h_j = cell_data[j[edge]];
 			} else { //boundary
-				h_j = boundary_heights[j[edge]-n];
-			}
-			v[edge] = -0.5*(h[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
-			v_i += 0.5*(h[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
+				h_j = bdry_data[j[edge]-n];
+			}
+			v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
+			v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
 		}
 		//Organise the set of 4 values/indices into the data and colind arrays
@@ -131,4 +132,50 @@
 }
 
+int update_operator_matrix(int n,
+                          int tot_len,
+                          long *geo_indices,
+                          double *geo_values,
+                          double *cell_data,
+                          double *bdry_data,
+                          double *data,
+                          long *colind) {
+	int i, k, edge, j[4], *sorted_j, this_index;
+	double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
+	for (i=0; i<n; i++) {
+		v_i = 0.0;
+		j[3] = i;
+		//Get the values of each interaction, and the column index at which they occur
+		for (edge=0; edge<3; edge++) {
+			j[edge] = geo_indices[3*i+edge];
+			if (j[edge]<n) { //interior
+				h_j = cell_data[j[edge]];
+			} else { //boundary
+				h_j = bdry_data[j[edge]-n];
+			}
+			v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
+			v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
+		}
+		//Organise the set of 4 values/indices into the data and colind arrays
+		sorted_j = quicksort((int *)j,4);
+		for (k=0; k<4; k++) { //loop through the nonzero indices
+			this_index = sorted_j[k];
+			if (this_index == i) {
+				data[4*i+k] = v_i;
+				colind[4*i+k] = i;
+			} else if (this_index == j[0]) {
+				data[4*i+k] = v[0];
+				colind[4*i+k] = j[0];
+			} else if (this_index == j[1]) {
+				data[4*i+k] = v[1];
+				colind[4*i+k] = j[1];
+			} else { //this_index == j[2]
+				data[4*i+k] = v[2];
+				colind[4*i+k] = j[2];
+			}
+		}
+	}
+	return 0;
+}
+
 /**
 * Python wrapper methods
@@ -147,5 +194,5 @@
     
     //Convert Python arguments to C
-    if (!PyArg_ParseTuple(args, "Oii", &kv_operator, &n, &tot_len)) {
+    if (!PyArg_ParseTuple(args, "Oii", &kv_operator)) {
         PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input");
         return NULL;
@@ -158,4 +205,7 @@
 
     //Extract parameters
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
     centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
     neighbours = get_consecutive_array(mesh,"neighbours");
@@ -195,6 +245,6 @@
 	int n, tot_len, err;
 	PyArrayObject
-		*h,
-		*boundary_heights,
+		*cell_data,
+		*bdry_data,
 		*geo_indices,
 		*geo_values,
@@ -203,9 +253,14 @@
 	
 	//Convert Python arguments to C
-    if (!PyArg_ParseTuple(args, "OiiOO", &kv_operator, &n, &tot_len, &h, &boundary_heights)) {
-        PyErr_SetString(PyExc_RuntimeError, "get_stage_height_interactions could not parse input");
-        return NULL;
-    }
-	
+    if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
+        PyErr_SetString(PyExc_RuntimeError, "build_operator_matrix could not parse input");
+        return NULL;
+    }
+
+
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
+
 	geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
 	geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
@@ -216,6 +271,6 @@
 								(long *)geo_indices -> data, 
 								(double *)geo_values -> data, 
-								(double *)h -> data, 
-								(double *)boundary_heights -> data, 
+								(double *)cell_data -> data,
+								(double *)bdry_data -> data,
 								(double *)_data -> data, 
 								(long *)colind -> data);
@@ -233,7 +288,57 @@
 }
 
+
+static PyObject *py_update_operator_matrix(PyObject *self, PyObject *args) {
+	PyObject *kv_operator;
+	int n, tot_len, err;
+	PyArrayObject
+		*cell_data,
+		*bdry_data,
+		*geo_indices,
+		*geo_values,
+		*_data,
+		*colind;
+
+	//Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
+        PyErr_SetString(PyExc_RuntimeError, "update_operator_matrix could not parse input");
+        return NULL;
+    }
+
+
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
+
+	geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
+	geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
+	_data = get_consecutive_array(kv_operator,"operator_data");
+	colind = get_consecutive_array(kv_operator,"operator_colind");
+
+	err = update_operator_matrix(n,tot_len,
+								(long *)geo_indices -> data,
+								(double *)geo_values -> data,
+								(double *)cell_data -> data,
+								(double *)bdry_data -> data,
+								(double *)_data -> data,
+								(long *)colind -> data);
+    if (err != 0) {
+        PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
+        return NULL;
+    }
+
+	Py_DECREF(geo_indices);
+	Py_DECREF(geo_values);
+	Py_DECREF(_data);
+	Py_DECREF(colind);
+
+    return Py_BuildValue("");
+
+}
+
 static struct PyMethodDef MethodTable[] = {
         {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
 		{"build_operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},
+        {"update_operator_matrix",py_update_operator_matrix,METH_VARARGS,"Print out"},
         {NULL,NULL,0,NULL} // sentinel
 };
Index: /trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py	(revision 8134)
+++ /trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py	(revision 8135)
@@ -29,8 +29,16 @@
         domain.set_boundary({'edge0': D0, 'edge1': D1, 'edge2': D2})
 
+        domain.set_quantity('stage', lambda x,y : x+2*y )
+        domain.set_quantity('elevation', lambda x,y : 3*x+5*y )
+
+
+        #print domain.quantities['stage'].vertex_values
+
+        #print domain.quantities['stage'].edge_values
+
         domain.update_boundary()
 
 
-        print domain.quantities['elevation'].boundary_values
+        #print domain.quantities['stage'].boundary_values
         
         return Kinematic_Viscosity(domain)
@@ -88,14 +96,21 @@
         assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]]))
 
-    def test_apply_heights(self):
+    def test_elliptic_matrix(self):
         operator1 = self.operator1()
         operator2 = self.operator2()
 
-        h = num.array([1.0])
-        A = operator1.apply_stage_heights(h)
+        domain1 = operator1.domain
+        domain2 = operator2.domain
+
+
+        A = operator1.elliptic_operator_matrix
+
+        print A
+        print num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+        
         assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
 
-        h = num.array([2.0])
-        A = operator1.apply_stage_heights(h)
+
+        A = operator1.elliptic_operator_matrix
         assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]))
 
