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)?