How to find eigenvalues and eigenvectors of a matrix?

1

Hi all!

I went through all the posts related to finding the eigenvalues and the eigenvectors of a matrix constructed from a solution, but couldn’t implement them properly.

In the below MWE, I am finding the solution u to the Poisson equation on a disc. Next, I am creating a vector

U = as_vector((0, u))

Then I need to find the eigenvalues of F^T F, where

d = mesh.geometry().dim()
I = Identity(d)
F = I + grad(U)

Here is the MWE:

from fenics import *
from mshr import Circle, generate_mesh
import numpy as np

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary        


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

C = Circle(Point(0, 0), 3)

mesh = generate_mesh(C, 100)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)

W1 = FunctionSpace(mesh, "Lagrange", 1)

n = FacetNormal(mesh)

G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

U = as_vector((0, u))
d = mesh.geometry().dim()
I = Identity(d)
F = I + grad(U)
C = PETScMatrix()             
C = assemble(F.T*F, tensor = C);

This raises the error:

C = assemble(F.T*F, tensor = C);
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 64, in _create_dolfin_form
    raise TypeError("Invalid form type %s" % (type(form),))
TypeError: Invalid form type <class 'ufl.tensors.ComponentTensor'>

Any suggestions would be helpful.

2

The expression you have for C=F.T*F is a tensor field, but FEniCS doesn’t assemble those kinds of expressions. You can only assemble scalar quantities (for example the total energy in the domain) and variational forms.

To find C=F.T*F, you can try defining scalar Function(W1) instances to store each component of C. You will have to use the project function to project components of C into each of those functions. You could then find the eigenvalues of C at different nodal locations by reforming the appropriate matrix from the Function(W1) instances and using a library like numpy.