Projecting Cauchy stress tensor gives different results

Hi, I am trying to solve a torsion of an hyperelastic cylinder. Then I want to compute the normal force and the torque, using the formulas
N=\int^{2\pi}_{0}\int_{0}^{r_0}\sigma_{zz}r dr d\theta and T=\int^{2\pi}_{0}\int_{0}^{r_0}\sigma_{\theta z}r^2 dr d\theta. Since I am working in Cartesian coordinates I trasform these integrals and then I need to compute \sigma_{\theta z} from the component I have. I wrote the following code. The problem is that I obtain always different values for the normal force every time I run the code. I think is something linked with the lagrangian multiplier p, maybe I am taking a wrong functional space. But when I try to change the functional spaces I am using, it gives me errors and the solver does not converge. Does someone know what is going on?

from dolfin import *
from mshr import *
from ufl import replace
import numpy as np
import csv
import matplotlib.pyplot as plt
from math import *

# Second order identity tensor
def SecondOrderIdentity(u):
d = u.geometric_dimension()
return variable(Identity(d))

I = SecondOrderIdentity(u)

# Determinant of the deformation gradient
def Jacobian(u):
return variable(det(F))

# Right Cauchy-Green tensor
def RightCauchyGreen(u):
return variable(F.T*F)

# Left Cauchy-Green tensor
def LeftCauchyGreen(u):
return variable(F*F.T)

# Invariants of an arbitrary tensor, A
def Invariants(A):
I1 = tr(A)
I2 = 0.5*(tr(A)**2 - tr(A*A))
I3 = det(A)
return [variable(I1), variable(I2), variable(I3)]

# Invariants of the (right/left) Cauchy-Green tensor
def CauchyGreenInvariants(u):
C = RightCauchyGreen(u)
[I1, I2, I3] = Invariants(C)
return [variable(I1), variable(I2), variable(I3)]

# Elastic Cauchy stress tensor
def elastic_Cauchy(u,p):
C = RightCauchyGreen(u)
B = LeftCauchyGreen(u)
J = Jacobian(u)
I = SecondOrderIdentity(u)
I1, I2, I3 = CauchyGreenInvariants(u)

psi = (mu_0/2-c2)*(I1 - 3) + c2*(I2 - 3)
W1 = diff(psi,I1)
W2 = diff(psi,I2)

Te = 2*W1*B-2*W2*inv(B)-p*I
return variable(Te)

# First Piola-Kirchhoff stress tensor
def FirstPiola(u,p):
J = Jacobian(u)
T = elastic_Cauchy(u,p)
return variable(J*T*inv(F).T)

def geometry_3d(r1,r2,h1,h2):
#Making a cylindrical geometry
#cylinder = Cylinder('coordinate of center of the top circle', 'coordinate of center of the bottom circle', 'radius of the circle at the top', 'radius of the circle at the bottom')

cylinder = Cylinder(Point(0, 0, h2), Point(0, 0, h1), r1, r2)
geometry = cylinder
# Making Mesh (30 corresponds to the mesh density)
mesh = generate_mesh(geometry, 30)

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)

bottom = AutoSubDomain(lambda x: near(x[2], h1))
top = AutoSubDomain(lambda x: near(x[2], h2))

bottom.mark(boundary_parts, 1)
top.mark(boundary_parts, 2)

return boundary_parts

# Create mesh and define function space ============================================
r1 = 12.5 #m
r2 = 12.5 #m
h1 =  0.0 #m
h2 = 10.0 #m
facet_function = geometry_3d(r1,r2,h1,h2)
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
ds = Measure("ds", domain=mesh, subdomain_data=facet_function)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())

mesh_file = File("mesh.pvd")
mesh_file << mesh

dx = dx(degree=4)
ds = ds(degree=4)

### Create function space

P2 = VectorElement("CG", mesh.ufl_cell(), 1) # displacement u
P1 = FiniteElement("CG", mesh.ufl_cell(), 1) # pressure p
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
TH = MixedElement([P2, P1])
V = FunctionSpace(mesh, TH)
W = FunctionSpace(mesh, DG)
TT = TensorElement("DG", mesh.ufl_cell(), 0)
TS = FunctionSpace(mesh, TT)

### Define functions for variational problems

# Incremental displacement and pressure
dup = TrialFunction(V)
(du, dp) = split(dup)

# Test functions for displacement and pressure
u_, p_ = TestFunctions(V)

# Displacement and pressure (current value)
up = Function(V)
(u, p) = split(up)

# Components of Cauchy stress tensor
sigma_theta_z = Function(W)
sigma_z_z = Function(W)

# Displacement from previous iteration
#h = Torsion() # Traction force on the boundary

# Time stepping parameters
dt = 0.1
t = 0.0
Tfin = 5.0

### Boundary conditions
twist = Expression(("x[0] == 0 ? 0 : sqrt(pow(x[0],2)+pow(x[1],2))*cos(atan2(x[1],x[0])+theta*x[2])-x[0]", "x[0] == 0 ? 0 : sqrt(pow(x[0],2)+pow(x[1],2))*sin(atan2(x[1],x[0])+theta*x[2])-x[1]", "0.0"), theta=0.0, degree=1)
bc_bottom = DirichletBC(V.sub(0), Constant((0, 0, 0)), facet_function, 1)
bc_top = DirichletBC(V.sub(0), twist, facet_function, 2)
bcs = [bc_bottom, bc_top]

### Input parameters
strain0 = 0.02
t_star =  0.3
mu_0 = Constant(1.0) #+ mu_2 + mu_3
c2 = Constant(2/3*mu_0) #Pa
r = Function(W)
r.interpolate(Expression("sqrt(pow(x[0],2)+pow(x[1],2))", degree=0))

J = Jacobian(u)
P = FirstPiola(u,p)

L1 = inner(P, Grad(u_))*dx #- dot(h, u_)*ds(2)
L2 = (J-1)*p_*dx
L = L1 + L2
j = derivative(L, up, dup)
problem = NonlinearVariationalProblem(L, up, bcs, J=j)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'

### Save in a file

displacement_file = File("u.pvd")
pressure_file = File("p.pvd")
normal_file = open("normal_num.csv", "w")
writer_normal = csv.writer(normal_file)
torque_file = open("torque_num.csv", "w")
writer_torque = csv.writer(torque_file)

tol = 1e-06

while t <= Tfin+tol:

t = round(t,2)
print('time: ', t)

# Increase torsion
if t <= t_star+tol:
alpha = (strain0/r2)*(t/t_star)
else:
alpha = 1/r2*strain0

twist.theta = alpha

# solve and save displacement
solver.solve()

u, p = up.split()

displacement_file << (u, t)
pressure_file << (p, t)

# Compute the normal and the torque
Tv = elastic_Cauchy(u,p)
sigma_v = project(Tv,TS) #Put here Cauchy stress tensor
sigma_values = sigma_v.vector().get_local()

for cell in cells(mesh):
i = cell.index()
x = cell.midpoint()[0]
y = cell.midpoint()[1]
z = cell.midpoint()[2]

sigma_xx=sigma_values[9*i]
sigma_yy=sigma_values[9*i+4]
sigma_zz=sigma_values[9*i+8]
sigma_xy=sigma_values[9*i+1]
sigma_xz=sigma_values[9*i+2]
sigma_yz=sigma_values[9*i+5]

theta_t = atan2(y,x)
sigma_theta_z.vector()[i] = -sin(theta_t)*sigma_xz+cos(theta_t)*sigma_yz
sigma_z_z.vector()[i] = sigma_zz

normal = assemble(sigma_z_z*ds(2))
torque = assemble(sigma_theta_z*r*ds(2))
writer_normal.writerow([t, normal])
writer_torque.writerow([t, torque])

# time increment
t += float(dt)

normal_file.close()
torque_file.close()