I am trying to solve a very modified Stokes flow which contains the term [tex]\int_{V_f} -\nabla p + \nabla^2 u dV[/tex] for (scalar) unknown p and vectorial unknown u. This means I need the Laplacian of u as I cannot reduce the order of derivatives through the weak form. In practice the weak form of this term is simply [tex](\int_{V_f} -\nabla p + \nabla^2 u dV,v)[/tex].
Initially I tried having u,p trial functions from VectorH1 and H1.
My various attempt to express this term includes the very naive "div(grad(u))", which gives the message that "div" cannot be formed.

I tried constructing the second order derivatives through "hesse" as in this example:
ngsolve.org/docu/nightly/i-tutorials/uni...der/fourthorder.html
Then I need to apply it to every component of the vector u, which gives problems when I try to reassemble again. I here create u1,u2,u3 (each coming from H1), and tried then for the weak form constructions like
grad(CoefficientFunction((u1,u2,u3)))
for when I make the weak form of my problem. But apparently this is also not possible.

Am I simply using the CoefficientFunction wrong?
Is there a simpler way to include these second order derivatives of a vector function?