Maxwell related question

Dear all,

I am trying to understand a bit more of FreeFEM’s type deduction strategy in the following case. If u is a P1 function, then grad(u) is an NE0 function, say in E1. If such a P1 function is known from solving a Dirichlet-Laplace problem and its gradient needs to be used as the right hand side (rhs) of a vectorial problem that is tested against NE0 functions, then which of the following ways is the way to go?
(1) first interpolate, then plug-in as rhs

E1 [jx, jy, jz] = -k*grad(u); // k is subdomain-wise constant of type func
varf RHS([Ax, Ay, Az], [Wx, Wy, Wz]) = int3d(T)( [Jx, Jy, Jz]' * [Wx, Wy, Wz] )

(2) just include the gradient as part of the rhs, that is,

varf RHS([Ax, Ay, Az], [Wx, Wy, Wz]) = int3d(T)( -k * grad(u)' * [Wx, Wy, Wz] )

Using the notation here, in the magnetostatics module the rhs is defined as in E1, but there it is clearer, since the current is a function of the space variables. Can someone please explain what is the space of k and grad(u) during the assembly in the second case (check the test space and project/interpolate there or? a link in documentation also works)? Further, in the same magnetostatics example, is there a reason for not defining Mu as a P0 function, since it naturally is a P0 member (and this will keep the module more transparent)?

As long as your qforder is high enough, both formulation will be equivalent numerically.