Is there a straight forward way to compute (and plot) $\vec n \cdot \nabla u$ on a boundary?
This was asked 3 years ago here, but there didn't seem to be a simple solution. Maybe this has changed by now?
My inclination is to first build n:
n = FacetNormal(mesh)
Then, assuming you're simply trying to include this in a weak formulation:
u = TrialFunction(U) v = TestFunction(U) problem = inner(nabla_grad(u)*n,v)*dx
OR You could simply compute grad(u) directly with the Dx() function (assuming u is a scalar function, and your domain is 2D, for simplicity):
problem = Dx(u,0)*n[0]*v*dx + Dx(u,1)*n[1]*v*dx