Calculate flux from head solution

Dear all,

I am solving a hillslope problem with standard Galerkin’s method, and I get the head solution. Then, I am trying to get the flux solution from my head solution. I know one way is use Q=-k(h)*grad(h). Is there any other inbuilt function in fenics can do the similar task?

What I do now is
grad_u=project(-K(h,z,alpha,N,M)*grad(u), VectorFunctionSpace(mesh, ‘Lagrange’, 1))
#K(h,z,alpha,N,M) is my k(h) function and u is the head solution.

If anyone has any ideas, it is very appreciated.

Thank you.