Issues with the gradient of the elastic energy

I’m working on a simple linear elasticity problem using Gridap.jl. I'm trying to compute the elastic energy in two ways:
Using the global expression: 0.5u'Ku
Using the local expression: sum(∫((0.5(ε(uh_Disp) ⊙ (σfun∘(ε(uh_Disp))))))*fem_params.dΩ)
Both expressions yield the same energy values, which is expected.
However, when I compute the derivative of these energy expressions with respect to the nodal displacements:
For the global case: K * u
For the local case: EngCont_DispFn(v) = ∫(1.0 * (ε(v) ⊙ (σfun ∘ ε(uh_Disp)))) * dΩ; grad_force_dis = assemble_vector(EngCont_DispFn, U0_Disp)
I notice a mismatch between the two resulting vectors.
Could you please help me identify what might be going wrong in the way I'm assembling the vector from the local energy derivative? I’ve attached the relevant part of the code for reference. The given problem is under force loading

```
function ElasFourthOrderConstTensor(E,ν)
C1111 = (E*(1-ν*ν))/((1+ν)*(1-ν-2*ν*ν))
C1122 = (ν*E)/(1-ν-2*ν*ν)
C1112 = 0.0
C2222 = (E*(1-ν))/(1-ν-2*ν*ν)
C2212 = 0.0
C1212 =E/(2*(1+ν))
C_ten = SymFourthOrderTensorValue(C1111 ,C1112 ,C1122 ,C1112 ,C1212 ,C2212 ,C1122 ,C2212 ,C2222)
return C_ten
end

const C_mat = ElasFourthOrderConstTensor(E_mat ,ν_mat);

function σfun(ε)
σ = C_mat⊙ε
return σ
end

order = 1
f = VectorValue(0.0,-1.0)

reffe_Disp = ReferenceFE(lagrangian ,VectorValue{2,Float64},order)
V0_Disp = TestFESpace(model,reffe_Disp;conformity =:H1,
dirichlet_tags = ["LeftSupport"],
dirichlet_masks =[(true,true)])
uApp1(x) = VectorValue(0.0,0.0)
U0_Disp = TrialFESpace(V0_Disp ,[uApp1])

degree = 2*order
Ω= Triangulation(model)
dΩ= Measure(Ω,degree)

labels = get_face_labeling(model)
LoadTagId = get_tag_from_name(labels ,"LoadLine")
Γ_Load = BoundaryTriangulation(model ;tags = LoadTagId)
dΓ_Load = Measure(Γ_Load ,degree)
n_Γ_Load = get_normal_vector(Γ_Load)

fem_params = (;V0_Disp, U0_Disp, Ω, dΩ)

function stepForceLoading(fem_params)
a_Disp(u,v) = ∫(ε(v) ⊙ (σfun∘(ε(u))))fem_params.dΩ
b_Disp(v) = ∫(v ⋅ f)dΓ_Load
op_Disp = AffineFEOperator(a_Disp ,b_Disp ,fem_params.U0_Disp ,fem_params.V0_Disp)
uh_out = solve(op_Disp)
return get_free_dof_values(uh_out)
end

function KMatForce(fem_params)

return assemble_matrix(fem_params.U0_Disp, fem_params.V0_Disp) do u, v
∫((∇(u))' ⊙ (C_mat ⊙ ∇(v)))fem_params.dΩ
end
end;

u_vecDisp = stepForceLoading(fem_params)
uh_Disp = FEFunction(fem_params.U0_Disp, u_vecDisp)

KMatD = KMatForce(fem_params)
EngDisc_Disp = 0.5*u_vecDisp'*KMatD*u_vecDisp

EngCont_Disp = sum(∫((0.5*(ε(uh_Disp) ⊙ (σfun∘(ε(uh_Disp))))))*fem_params.dΩ)
EngDisc_Disp, EngCont_Disp

EngCont_DispFn(v) = ∫(((ε(v) ⊙ (σfun∘(ε(uh_Disp))))))*dΩ
grad_force_dis = assemble_vector(EngCont_DispFn,U0_Disp)
[KMatD*u_vecDisp grad_force_dis]
```

Start with the KMatForce, stepForceLoading, and EngCont_DispFn entry points, comparing their use of gradients, strains, measures, and the U0_Disp/V0_Disp spaces. Reproduce the mismatch and trace how assemble_vector handles the local derivative; done when the assembled derivative agrees with KMatD*u_vecDisp for the supplied force-loading example.

Written by the indexing model from the issue text.