eigenvalue of matrix element

Question

If I have a matrix element, is there any command in fenics to compute the eigenvalues of this?

More specifically, suppose that I just solved the elasticity equation and have a 3D vector element u. I can then compute the gradient of this, which will live in a tensor function space. I'm looking for a command that will yield a 3D vector element whose components are the eigenvalues of this tensor. Does such a command exist?

Answer 1

No but it is quite easy to construct, if you have numpy.version >= '1.8.0'

import numpy
from dolfin import *

def get_eig(hes):
  mesh = hes.function_space().mesh()
  [eigL,eigR] = numpy.linalg.eig(hes.vector().array().reshape([mesh.num_vertices(),2,2]))
  eig = Function(VectorFunctionSpace(mesh,'CG',1))
  eig.vector().set_local(eigL.flatten())
  return eig

mesh = RectangleMesh(Point(0.,0.),Point(2.,1.),20,10)
u = interpolate(Expression(("cos(x[0])","sin(x[1])")),VectorFunctionSpace(mesh,'CG',2))
hes = project(sym(grad(u)),TensorFunctionSpace(mesh,'CG',1))
eig = get_eig(hes)
plot(eig, interactive=True)

Comments

Thanks a lot. That works well when I use one processor, but it doesn't work with I use MPI. Do you know how to modify what you wrote so that it works for parallel computing.

---

You can exploit the fact that explicit expressions for symmetric 2x2 matrices.

import numpy
from dolfin import *

def get_eig(hes):
  mesh = hes.function_space().mesh()
  S00_ = hes.sub(0)
  S01_ = hes.sub(1)
  S11_ = hes.sub(3)
  eig = interpolate(Expression(("0.5*(S00+S11-sqrt((S00-S11)*(S00-S11)+4.*S01*S01))",\
                                "0.5*(S00+S11+sqrt((S00-S11)*(S00-S11)+4.*S01*S01))"),\
        S00=S00_,S01=S01_,S11=S11_),VectorFunctionSpace(mesh,'CG',1))
  return eig

mesh = RectangleMesh(Point(0.,0.),Point(2.,1.),20,10)
u = interpolate(Expression(("cos(x[0])","sin(x[1])")),VectorFunctionSpace(mesh,'CG',2))
hes = project(sym(grad(u)),TensorFunctionSpace(mesh,'CG',1))
eig = get_eig(hes)
plot(eig, interactive=True)

A similar result exists for 3x3 matrices, see https://en.wikipedia.org/wiki/Eigenvalue_algorithm

Edit: replaced projections with .sub() functionality.