So, after having a first look: in the NGSolve version, "source" does not have boundary conditions applied, i.e. you still have reaction forces corresponding to nodes with Dirichlet boundary conditions.
You can kick these out with
Code:
nRes = Norm(bc_projector * source.vec)
print("Norm res:", nRes)
where is defined upfront as
Code:
bc_projector = Projector(V.FreeDofs(), True)
This will give you vanishing residual in the elastic regime (i.e. Newton converges after first step).
In the inelastic regime, there seems to be an issue with "Jacobian" stress as Newton does not converge.
I did not spot the bug though.
Besides that, I found it unusual, that you have used (u...trial func, v... tst func)
Code:
InnerProduct(Strain(u), Jacobian_stress(Strain(v)))
This is consistent with your FEniCS implementation, but I would have guessed that you need
Code:
InnerProduct(Strain(v), Jacobian_stress(Strain(u)))
instead. In my test, it didn't make a difference though, maybe shadowed by "the" bug?
(didn't run FEniCS yet to see whether it makes a difference there)