How to find eigenvalues and eigenvectors of a tensor for post processing?

General
1

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.

2

See for instance: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/eigenvalue/demo_eigenvalue.py.html?highlight=eigenvalues

3

Hi @dokken,

Thank you for the fast reply.

Please bare with me if I sound like a beginner.

In the example you provided, the variational form is solved using the SLEPcEigenSolver. However, in my case, since I am solving a fully coupled system of equations, I already have the displacement solution, hence the stress tensor. So, would you please guide me on using SLEPcEigenSolver to solve for the eigenvalues and eigenvectors without including the variational form?

4

You can find principal components and principal directions using closed-form expression in UFL, see dolfiny/invariants.py at master · michalhabera/dolfiny · GitHub

These UFL expression have to be “evaluated” for specific FE functions as input (displacements or stress tensor is your case). In old FEniCS please use project() for this, in FEniCSX use Expressions and interpolation to get the exact point evaluations.

1 Like
5

Hi @michalhabera,

Thank you for the reply. The link you provided had a significant contribution to solving my problem. For anyone who is facing the same problem, below is the full solution:

The link provided by @michalhabera here can be used to find the eigenvalues. However, the eigenvectors solution was not included. As a result, I used the directional cosine to find the eigenvectors. The explanation of the mathematical approach can be found in the YouTube video here. Below is the code to solve for the eigenvectors:

def eigvector(A,λ):

   

    h = A-λ* ufl.Identity(3)

    a1 = as_tensor([
        [h[1,1], h[1,2]],
        [h[2,1], h[2,2]]
    ])

    b1 = as_tensor([
        [h[1,0], h[1,2]],
        [h[2,0], h[2,2]]
    ])

    c1 = as_tensor([
        [h[1,0], h[1,1]],
        [h[2,0], h[2,1]]
    ])

    A1= det(a1)
    B1= det(b1)
    C1 =det(c1)

    n1 = ufl.sqrt(A1**2+B1**2+C1**2)

    nx1 = A1/n1
    ny1 = B1/n1
    nz1 = C1/n1


    
    v= as_vector([nx1,ny1,nz1])

    # v= as_tensor([
    #     [nx1, 0.0, 0.0],
    #     [0.0, ny1, 0.0],
    #     [0.0, 0.0, nz1]
    # ])


    return v

Since I am using the legacy FEniCS, I used the function eigenstate3_legacy(A) to solve for the eigenvalues. Below is the code from here plus the eigenvectors:

def eigenstate3_legacy(A):
    """Eigenvalues and eigenprojectors of the 3x3 (real-valued) tensor A.
    Provides the spectral decomposition A = sum_{a=0}^{2} λ_a * E_a
    with eigenvalues λ_a and their associated eigenprojectors E_a = n_a^R x n_a^L
    ordered by magnitude.
    The eigenprojectors of eigenvalues with multiplicity n are returned as 1/n-fold projector.
    Note: Tensor A must not have complex eigenvalues!
    """
    if ufl.shape(A) != (3, 3):
        raise RuntimeError(f"Tensor A of shape {ufl.shape(A)} != (3, 3) is not supported!")
    #
    eps = 1.0e-10
    #
    A = ufl.variable(A)
    #
    # --- determine eigenvalues λ0, λ1, λ2
    #
    # additively decompose: A = tr(A) / 3 * I + dev(A) = q * I + B
    q = ufl.tr(A) / 3
    B = A - q * ufl.Identity(3)
    # observe: det(λI - A) = 0  with shift  λ = q + ω --> det(ωI - B) = 0 = ω**3 - j * ω - b
    j = ufl.tr(B * B) / 2  # == -I2(B) for trace-free B, j < 0 indicates A has complex eigenvalues
    b = ufl.tr(B * B * B) / 3  # == I3(B) for trace-free B
    # solve: 0 = ω**3 - j * ω - b  by substitution  ω = p * cos(phi)
    #        0 = p**3 * cos**3(phi) - j * p * cos(phi) - b  | * 4 / p**3
    #        0 = 4 * cos**3(phi) - 3 * cos(phi) - 4 * b / p**3  | --> p := sqrt(j * 4 / 3)
    #        0 = cos(3 * phi) - 4 * b / p**3
    #        0 = cos(3 * phi) - r                  with  -1 <= r <= +1
    #    phi_k = [acos(r) + (k + 1) * 2 * pi] / 3  for  k = 0, 1, 2
    p = 2 / ufl.sqrt(3) * ufl.sqrt(j + eps ** 2)  # eps: MMM
    r = 4 * b / p ** 3
    r = ufl.Max(ufl.Min(r, +1 - eps), -1 + eps)  # eps: LMM, MMH
    phi = ufl.acos(r) / 3
    # sorted eigenvalues: λ0 <= λ1 <= λ2
    λ0 = q + p * ufl.cos(phi + 2 / 3 * ufl.pi)  # low
    λ1 = q + p * ufl.cos(phi + 4 / 3 * ufl.pi)  # middle
    λ2 = q + p * ufl.cos(phi)  # high
    #

    # Determine Eigenvectors
    v0 = eigvector(A,λ0)
    v1 = eigvector(A,λ1)
    v2 = eigvector(A,λ2)

    # Eigenvalues
    u = as_tensor([
        [λ0, 0.0, 0.0],
        [0.0, λ1, 0.0],
        [0.0, 0.0, λ2]
    ])

     return u, v0, v1, v2
6

Hello @Tareq_M_Al-Shami. Thank you for providing the function to calculate the eigenvectors. Please note that you forgot to multiply the determinants A1 and C1 times -1.