Hi everyone,

I am working on a full coupling problem of three differential equations. One of the equations is the linear elasticity equation, in which we solve for the displacement *u*. Then, using the displacement solution, we can solve for the strain and eventually solve for the stress tensor. For postprocessing purposes, I want to get the principle stress and their respective directions. To do that, we have to find the eigenvalues and eigenvectors. I have tried to use NumPy to do that, but I seem not to get it right. I would appreciate it if I could get any guidance on how to find the eigenvalues and eigenvectors.

Below is an extended code from the 2D linear elasticity example in here. I have extended the 2D elasticity problem into 3D linear elasticity. Below *sig* is the stress tensor I want to use to solve for the Eigenvalues and Eigenvectors.

This is my first time asking a question on this platform, so I hope the message is properly written and formatted.

```
from fenics import *
### Domain
x0 ,x1 = 0., 25.
y0, y1 = 0., 1.
z0, z1 = 0., 1.
nx, ny, nz = 30, 10, 10
mesh = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),nx,ny,nz)
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
return sym(grad(v))
def sigma(v):
return lmbda*tr(eps(v))*Identity(len(v)) + 2.0*mu*eps(v)
rho_g = 1e-3
f = Constant((0., 0.,-rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange',2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0],0.)
bc = DirichletBC(V, Constant((0.,0., 0.)), left)
u = Function(V, name="Displacement")
solve(a == l, u, bc)
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
file_results = XDMFFile("elasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
```

I appreciate any help you can provide.