I want to compute the divergence of some vector functions. a=\nabla \cdot(gf \nabla \rho)
If I create a functionspace V, and \rho,f,g are the interpolated functions of V, how to compute a? I hope that a is also an interpolated functions in V.
’ ’ ’
V = fem.functionspace(domain, (“Lagrange”, 1))

rho = fem.Function(V)

rho.interpolate(lambda x: x[0]**2 + x[1]**2) # Example expression for rho

Thanks for you reply. So is it possible to compute and restore the gradient of \rho in the 1 degree space? In that case I do not need to build a high degree function space.

You can get the gradient of \rho, which will be a piecewise constant function, thus it divergence would be zero.
If you use a second order lagrange function for \rho; the gradient will be piecewise linear (discontinuous across element boundaries), and the divergence (within each element) will be a piecewise constant.