AutoDiff of interpolate_everywhere fails. Mismatch among FE space vector type and dual number types.

Hi @fverdugo,

the following MWE reproduces an error we have found when using AutoDiff (ForwardDiff) to differentiate a function that calls Gridap's `interpolate_everywhere`.

As shown in the below, the error can be workarounded using `interpolate_everywhere!` and vectors of DoFs of the appropriate `eltype` (i.e.. dual number types).

The cause of the error is a mismatch among the FE space vector type (Float64) and the result of evaluating the function to be interpolated with dual numbers.

Do you see a neat approach to solve this issue? I guess that an ad-hoc differentiation rule for interpolate_everywhere that calls interpolate_everywhere! under the hood to compute the derivative might be one possible approach. Anything else that come intos your mind? E.g. creating a FE space with dual numbers as vector type with e.g., `gradient(::FESpace)` ?

Thanks!

```julia

using Gridap
using FillArrays
using Test

function dP_dp_error(params)
function g(params)
f(x)=(params[1]+params[2]*params[2]*x[1]+(params[3]^3)*x[2])
VectorValue(get_free_dof_values(
(Gridap.FESpaces.interpolate_everywhere(f,P))))
end
gradient(g)(VectorValue(params))
end

function dP_dp_workaround(params)
function g(params)
f(x)=(params[1]+params[2]*params[2]*x[1]+(params[3]^3)*x[2])
free_values = Vector{eltype(params)}(undef,num_free_dofs(P))
dirichlet_values = Vector{eltype(params)}(undef,num_dirichlet_dofs(P))
VectorValue(get_free_dof_values(
(Gridap.FESpaces.interpolate_everywhere!(f, free_values,dirichlet_values,P))))
end
gradient(g)(VectorValue(params))
end

domain = (0,1,0,1)
ncells = (5,5)
model = CartesianDiscreteModel(domain,ncells)

# FE space to represent parameter space
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
P = TestFESpace(model,reffe; conformity=:H1)

params=[1.0,2.0,3.0]
dP_dp_params=dP_dp_workaround(params)
dP_dp_params=dP_dp_error(params)
```