Computing residual directly from energy

Consider for instance a usual hyperelasticity problem (#398) which is basically finding the minimizer of a scalar potential energy (Pi):
$\arg\min_u \Pi (u)$

It would be useful to know constructs to be able to directly calculate the residual from the energy using automatic differentiation (The residual to jacobian part is already taken care of by Gridap). This is particularly useful in more complicated constitutive relations, where calculating the residual by hand would be time-consuming and prone to algebraic errors.

### MWE

```julia
using Gridap
using LineSearches: BackTracking
E = 10.
ν = 0.3
const μ = E/2/(1+ν)
const λ = 2*μ*ν/(1-2*ν)
F(∇u) = one(∇u) + ∇u
J(F) = det(F)
B(x) = VectorValue(0.0, -0.5, 0.0)
T(x) = VectorValue(-0.1, 0.0, 0.0)
ψ(F) = μ * (0.5 * (F ⊙ F - 3.) - log(J(F))) + λ/2. * (J(F) - 1.)^2
S(F) = ∇(ψ)(F)
energy(u) = ∫( ψ∘F∘∇(u) )*dΩ
res(u, v) = ∫( (S∘F∘∇(u)) ⊙ ∇(v) - B ⊙ v)*dΩ
domain = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
partition = (5, 5, 5)
model = CartesianDiscreteModel(domain, partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels, "rightFace", [2, 4, 6, 8, 14, 16, 18, 20, 26])
add_tag_from_tags!(labels, "leftFace", [1, 3, 5, 7, 13, 15, 17, 19, 25])

reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, 1)
V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags = ["leftFace", "rightFace"])

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

nls = NLSolver(
show_trace=true,
method=:newton,
linesearch=BackTracking()
)

g0(x) = VectorValue(0.0, 0.0, 0.0)
g1(x) = VectorValue(
0.0,
(0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
(0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
)

U = TrialFESpace(V, [g0, g1])

#FE problem
op = FEOperator(res, U, V)
solver = FESolver(nls)

x0 = zeros(Float64, num_free_dofs(V))
uh = FEFunction(U, x0)
uh, = solve!(uh, solver, op)
```

### Problem:
How to calculate the gradient of the energy? If I change the definition of residual to

```julia
function S(∇u)
μ * (F(∇u) - inv(F(∇u))') + λ * (J(F(∇u))) * (J(F(∇u)) - 1.) * inv(F(∇u))'
end

function res(u, v)
∫( (S∘∇(u)) ⊙ ∇(v) - B ⊙ v )*dΩ - ∫( T ⊙ v )*dΓ_t
end
```

everything works fine as it should. Also if I calculate the energy from the solution `uh`, namely
```julia
julia> sum(energy(uh))
0.11957195697992658
```

I get the correct result. It is just that the calculation of the residual from the energy using AD that is proving to be tricky.

Sorry for a bit long post, but I hope that the issue is clear :)