I am working on a existing Oasis solver for Navier Stokes (IPCS_ABCN) in which I want to introduce varying viscosity. In the original solver, the viscosity is constant, and is used in the following operation:

```
A.axpy(-0.5*nu, K, True)
```

Meaning:

```
A = nu*K + A
```

Here A and K are of types ‘dolfin.cpp.la.Matrix’ and nu, the viscosity, is of type ‘float’. A and K are defined as:

```
# Mass matrix
M = assemble_matrix(inner(u, v)*dx)
# Stiffness matrix (without viscosity coefficient)
K = assemble_matrix(inner(grad(u), grad(v))*dx)
# Allocate coefficient matrix (needs reassembling)
A = Matrix(M)
```

To introduce the varying viscosity, I defined nu as a function in stead of a constant. Ofcourse, the operation shown above with axpy doesnt work anymore, because nu is not a constant. Does anyone know how to do this operation correctly in my case?

Thanks in advance!