You can use the `as_matrix`

function to build up arbitrary matrices in UFL, as demonstrated below:

```
from dolfin import *
# Set up mesh, function space, and BC:
N = 8
mesh = UnitCubeMesh(N,N,N)
V = FunctionSpace(mesh,"CG",1)
bc = DirichletBC(V,Constant(0),"on_boundary")
# Example problem data:
def D(u):
# Gradient of `u` along diagonal of a matrix:
return -as_matrix([[u.dx(0), 0 , 0 ],
[0 , u.dx(1), 0 ],
[0 , 0 , u.dx(2)]])
E = Constant((1,1,1))
g = Constant(1)
# Set up weak form and solve:
u = TrialFunction(V)
v = TestFunction(V)
# Note: With these definitions, `D(u)*E` is
# equivalent to `-grad(u)`, so we're just solving the
# Poisson equation with source term `g`.
a = -dot(D(u)*E,grad(v))*dx # (Integrated `div` by parts)
L = g*v*dx
uh = Function(V)
solve(a==L,uh,bc)
# Output results:
File("u.pvd") << uh
```