Issues with the gradient of the elastic energy

[merclogyk]

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]