How to calculate eigenvalues and eigenvectors of the strain tensor in case 3D?

Hello dear all, I am now modeling a 3D problem using the phase-field method. How to calculate eigenvalues and eigenvectors of the strain tensor in 3D? Can anyone guide me?

Best,

zgh,

You should post the problem you are working on or at least a minimum working example of what you have tried so others can help you. See here.

I can’t directly answer your question but I can provide a temporary work around.

Personally I encountered this same exact question a few months ago when I wanted to obtain the principle components of a symmetric stress tensor after essentially solving the Fenics tutorial for elasticity example.

-\nabla * \sigma(u) = f in \Omega

The closest I got was iterating through each of the vertices, getting the correct eigenvalues and eigenvectors using Numpy, but was unable to correctly project back onto the solution space. So as for doing this programmatically I am unable to help.

However, if you need this soon, the way I got around this was to export the projected solution \sigma(u) onto a TensorFunctionSpace and output as a “.xdmf”. I then opened this in Paraview and exported the tensor from the spreadsheet view as a “.csv”. Using another software (Matlab), I then could easily grab the stress tensor - calculate the components - and then generate a spreadsheet where I project the eigen values and vectors back onto the original solution space (to be view in ParaView). I found this to be optimal especially since my application involved 500k+ elements and even with MPI this became cumbersome to try in Python. This is my personal experience take it or leave it - I’m welcome to provide further detail if this is of interest.

-Sven

Dear Sven,
Thank you so much for your nice descriptions. More precisely in one part of the implementation of the model needs to eigenvalues and eigenvectors of strain tensor in case 3D. I found it only in case of 2D as follows:

import dx
import sympy
import numpy
from dolfin import *
#----------------------------------------------------------------------#
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = 'uflacs'
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Compute the symbolic expression for eigenvalues by sympy
T = sympy.Matrix(2, 2, lambda i, j: sympy.Symbol('T[%d, %d]' % (i, j), symmetric =True, real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)
#eig_expr = T.LDLdecomposition()

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv = sympy.Matrix(eigv) # Column vector
eigw = sympy.Matrix(eigw) # Row has the elements
eigw = sympy.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sympy.Matrix([[eigv[0], 0.0],[0.0, eigv[1]]])

eigv_expr = map(str, eigv)
eigw_expr = map(str, eigw)

#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Create UFL operators for the eigenvalues and eigenvectors

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return map(eval, eigw_expr)

#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Import mesh from .xml
mesh = Mesh('mesh.xml')
# Introduce manually the material parameters
Gc =  2.7
l = 0.015
lmbda = 121.1538e3
mu = 80.7692e3
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# set up solution functions
u,unew, uold = Function(W,name='displacement'), Function(W,name='displacement'), Function(W,name='displacement')
pnew, pold, Hold = Function(V,name='phi'), Function(V,name='phi'), Function(V,name='phi')
# The way the eigenvalues are computed we cannot allow a constant value of u at start

unew_array = numpy.array(unew.vector())
unew_array = numpy.random.rand(len(unew_array))
unew.vector()[:] = unew_array

uold_array = numpy.array(uold.vector())
uold_array = numpy.random.rand(len(uold_array))
uold.vector()[:] = uold_array

u_array = numpy.array(u.vector())
u_array = numpy.random.rand(len(u_array))
u.vector()[:] = u_array
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Boundary conditions and control
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")

load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(0), load, top)
bc_u = [bcbot, bctop]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
#performing strain tensor spectral decomposition
def epsilon(u):
    return sym(grad(u))

def strn(u):
    t = sym(grad(u))

    v00, v01, v10, v11 = eigv(t)

    w00, w01, w10, w11 = eigw(t)

    w = ([w00,w01],[w10,w11])
    w = as_tensor(w)

    v = ([v00,v01],[v10,v11])

    v = as_tensor(v)

    return w*v*inv(w)

# Positive strain
def strn_p(u):
    t = sym(grad(u))
    v00, v01, v10, v11 = eigv(t)

    v00 = conditional(gt(v00,0.0),v00,0.0)
    v11 = conditional(gt(v11,0.0),v11,0.0)

    w00, w01, w10, w11 = eigw(t)
    wp = ([w00,w01],[w10,w11])
    wp = as_tensor(wp)
    vp = ([v00,v01],[v10,v11])
    vp = as_tensor(vp)
    return wp*vp*inv(wp)
# Negative strain
def strn_n(u):
    t = sym(grad(u))
    v00, v01, v10, v11 = eigv(t)
    v00 = conditional(lt(v00,0.0),v00,0.0)
    v11 = conditional(lt(v11,0.0),v11,0.0)
    w00, w01, w10, w11 = eigw(t)
    wn = ([w00,w01],[w10,w11])
    wn = as_tensor(wn)
    vn = ([v00,v01],[v10,v11])
    vn = as_tensor(vn)
    return wn*vn*inv(wn)
#----------------------------------------------------------------------#
# Define stress Constituive functions
def sigma(u):

    #return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
    # To compare the answers
    tr_p=0.5 * (tr(epsilon(u)) + abs(tr(epsilon(u))))
    tr_n=0.5 * (tr(epsilon(u)) - abs(tr(epsilon(u))))
    return ((1-pold)**2)*(lmbda * tr_p * Identity(2) + 2.0 * mu * strn_p(u))+lmbda * tr_n * Identity(2)+ 2.0 * mu * strn_n(u)
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# define positive elastic energy density and introduce a strain-history field
def psi(u):
    return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+\
           mu*inner(dev(epsilon(u)),dev(epsilon(u)))		
def H(uold,unew,Hold):
    return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
#Writing governing equations in weak form
E_du = inner(sigma(u),grad(v))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
            *inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)

Unfortunately, I do not well know your way and I am looking for a direct way to solve this problem.