How to construct a matrix of gradients depending on nodal locations in gridap

I want to calculate the gradients of residuals(displacement) with respect to another vector. based on mathematical calculations, as residuals are a vector this gradient should be a matrix. This matrix should has value when we are calculating the gradient of residuals wrt variable at the same node and it should be zero in other nodes. I couldnt find a way to directly get the node number and compare with other nodes to make others zero. so i wrote a function to logically do this as below

function vector_to_matrixut(vecu)
    mat = zeros(sr, ur)
    for i in 1:sr
        mat[i, (i-1)*2+1] = vecu[(i-1)*2+1]
        mat[i, 2*i] = vecu[2*i]
    return mat

the input parameter is the global vector of the gradients. as below:

function dσfun(ε_in,s_in)
    σ_elas = C_mat⊙ε_in
    if tr(ε_in) >= 0
        dσ = (2 * s_in)*σ_elas
    elseif  tr(ε_in) < 0
        dσ = (2 * s_in)*P_dev ⊙ σ_elas
    return  dσ
function  drus(fem_params,pth, uh_in,sh_in, vApp)
    uApp1(x) = VectorValue(0.0,0.0)
    uApp2(x) = VectorValue(0.0,0.0)
    uApp3(x) = VectorValue(0.0,-vApp)
    U_Disp = TrialFESpace(V0_Disp ,[uApp1 ,uApp2 ,uApp3])
    dA_Disp(v,pth,uh_in,sh_in) =  ((p->Em(p))∘pth) * (ε(v) ⊙ (dσfun∘(ε(uh_in), sh_in)))
    return  assemble_vector(v -> ∫(dA_Disp(v,pth,uh_in,sh_in))fem_params.dΩ ,U_Disp)

Is there any way that gridap do this calculation and considering nodal number without an explicit function?
The issue that i am facing here is that due to boundary conditions the number of that global vector is not matching with the number of nodes so some values of that vector can not be used.
let me know if you have any questions to clarify
thank you