This is definitiely not what you want to do, as sym(grad(v)) is a matrix, div(v) is a scalar value.
Going back to sym(grad(v))
You can
print(sigma(du).ufl_shape)
which gives you
(1,4)
which I guess you want to be (2, 2), thus:
def voigt2stress(s):
return as_tensor ([[s[0],s[2]],[s[2],s[1]]])
is probably what you want