This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

How to find eigenvalues and eigenvectors of a stress or strain tensor at an element level ?

+1 vote
asked Jan 30, 2017 by kaushikv123 FEniCS Novice (390 points)

Hi, do you mean as ufl operators so that you can use these in your forms?

commented Jan 31, 2017 by MiroK FEniCS Expert (80,920 points)

I mean are there ufl operators which i can use to compute eigenvalues and eigenvectors of stress and strain, else how would I do it ?

commented Jan 31, 2017 by kaushikv123 FEniCS Novice (390 points)

There are no such operators, but see here.

commented Jan 31, 2017 by MiroK FEniCS Expert (80,920 points)

Thank you. I will try this out and also explicitly write out modules for eigenvectors.

commented Jan 31, 2017 by kaushikv123 FEniCS Novice (390 points)

I have implemented the code for eigenvectors but i think i am making an error in the output format that i need to give for a vector in Fenics. I haven't been able to find how to do it. I will appreciate any help on this. Thanks in advance and here is the sample code below:

from dolfin import *
from sympy import Matrix
mesh = UnitSquareMesh(1, 1)
V = VectorFunctionSpace(mesh, 'CG', 2)
u = interpolate(Expression(('x[0]', 'x[0]*x[1]'), degree=2), V)

F = Identity(2) + grad(u)
C = F*F.T

# Eigenvalues are roots of characteristic polynomial
# e**2 - tr(A)*e + det(A) = 0
def eig_plus(A): 
    return (tr(A) + sqrt(tr(A)**2-4*det(A)))/2
def eig_minus(A): 
    return (tr(A) - sqrt(tr(A)**2-4*det(A)))/2
def eig_vecmat(A):
    lambdap = eig_plus(A)
    lambdan = eig_minus(A)
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]
    if c != 0:
       Eigvecmat = [[lambdap-d,lambdan-d],[c,c]]
    elif b != 0:
         Eigvecmat = [[b,b],[lambdap-a,lambdan-a]]
    elif b == 0 & c == 0:
         Eigvecmat = [[1,0],[0,1]]    
    return Eigvecmat
# Check
S = FunctionSpace(mesh, 'CG', 2)
T = TensorFunctionSpace(mesh, 'CG', 2)
f0 = project(eig_vecmat(C),T) 
plot(f0[0,0],interactive=True)
commented Jan 31, 2017 by kaushikv123 FEniCS Novice (390 points)

1 Answer

+1 vote
 
Best answer

Hi, the following shows the steps needed to define the operator (note that I didn't implement the full conditional part of your definition). Also the equality comparisons might fail with floats.

def eig_vecmat(A):
    lambdap = eig_plus(A)
    lambdan = eig_minus(A)
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]

    return conditional(ne(c, 0),
                       as_matrix(((lambdap-d,lambdan-d), (c,c))),
                       conditional(ne(b, 0),
                                   as_matrix(((b,b), (lambdap-a,lambdan-a))),
                                   Identity(2)))
answered Feb 1, 2017 by MiroK FEniCS Expert (80,920 points)
selected Feb 5, 2017 by kaushikv123

Thanks Mirok. I will add to this script that you gave me.

commented Feb 1, 2017 by kaushikv123 FEniCS Novice (390 points)
...