Hi all,

I am facing a weird issue while trying to calculate the Von Mises stress from a displacement field.

I obtain a negative result, which should not be as this Von mises stress is the defined as a square root of some kind of tensor L2 norm.

I am using the following code:

```
from fenics import *
from dolfin import *
from ufl import nabla_div
E = 62000
nu = 0.24
Sig_y = 7.0
mu = E/(2*(1+nu))
lambda_ = nu*E/((1+nu)*(1-2*nu))
Utop = -0.1
N = 16
h = 1/N
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), N, N, N)
tol=1E-8
class BoundaryBottom(SubDomain):
def inside(self, x, on_boundary):
return x[1] < h+tol
class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return x[1] > 1.0 - (h+tol)
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryTop = BoundaryTop()
boundaryBottom.mark(sub_domains, 1)
boundaryTop.mark(sub_domains, 2)
# mesh
vtk_file = File("VTK/Subdomains.pvd")
vtk_file << sub_domains
# Define Functional spaces V for variational formulation and projecting quantities
# sequential version
V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1)
Ws = TensorFunctionSpace(mesh, "CG", 1)
# functions and BC)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(1), Utop, boundaryTop)
bc=[bc1,bc2]
d=u.geometric_dimension()
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2*mu*epsilon(u) + lambda_*nabla_div(u)*Identity(d)
a = inner(sigma(u), epsilon(v))*dx
f= Constant((0,0,0))
L = dot(f,v)*dx
u=Function(V)
solve(a==L, u, bc, solver_parameters={"linear_solver": "petsc","preconditioner": "ilu"})
# Post treatment computations
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W, solver_type="cg")
vtkfile = File('VTK/u.pvd')
vtkfile << u
vtkfile = File('VTK/VM.pvd')
vtkfile << von_Mises
```

Any idea of where the problem could come from ?

I use a cg type projection as other type quickly lead to memory issue …

thanks in advance,

Quentin