Dear Andrew,

I'm afraid that I don't know the most efficient way to do this sort of many small eigensolves, but you could get your solution in `np.ndarray` using `q_soln.dat.data` (`q_soln.dat` returns `PyOP2.Dat` object: https://op2.github.io/PyOP2/user.html#pyop2.Dat).

You could also get the mesh data calling `mesh.coordinates.dat.data` .

For instance:

>>> from firedrake import *

>>> mesh = UnitSquareMesh(2,2)

>>> V = VectorFunctionSpace(mesh, "CG", 1)

>>> f = Function(V)

>>> f.dat.data

array([[0., 0.],

[0., 0.],

[0., 0.]])

>>> mesh.coordinates.dat.data

array([[0. , 0. ],

[0. , 0.5],

[0.5, 0. ],

[0.5, 0.5],

[1. , 0. ],

[0. , 1. ],

[1. , 0.5],

[0.5, 1. ],

[1. , 1. ]])

I think you can use these data in your context.

Taking a glance at the documentation, `numpy.linalg.eig` seems to take multiple matrices, so I think you can solve the node-wise eigen problems all at once.

Thank you,

Koki

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Andrew Hicks <ahick17@lsu.edu>
Sent: Friday, March 20, 2020 4:23 AM
To: firedrake <firedrake@imperial.ac.uk>
Subject: [firedrake] Returning eigenvectors of a solution to a PDE at each node

Hello,

I am currently using firedrake to solve a PDE. The solution to this PDE is vector-valued, and I have no problem finding the solution. However, I’d like to go a step further. I have named the solution “q_soln” and this is a 2-dimesional vector. I have used this vector to construct a 2x2 matrix at the node [0.5,0.5]. I then extracted an eigenvector from this matrix, and then I printed it:

import numpy as np

from numpy import linalg as la

E1 = (1 / np.sqrt(2)) * np.array([[1,0],[0,-1]])

E2 = (1 / np.sqrt(2)) * np.array([[0,1],[1,0]])

# as an example, we evaluate the solution at the point [0.5,0.5], and pick the first eigenvector

q_soln_atpoint = q_soln.at([0.5,0.5])

Q_soln = q_soln_atpoint[0] * E1 + q_soln_atpoint[1] * E2

evals, evecs = la.eig(Q_soln) # evec is a 2x2 array with the columns being the e-vectors

q_soln_evec = np.zeros((1,2))

q_soln_evec += evecs[0,:]

print(q_soln_evec)

My question is this: how can I get firedrake to do this automatically for every node in my mesh, and then store it in as a Function “q_soln_evec”? The only reason I used numpy to get this one point is because I’m unfamiliar with firedrake. I bet firedrake has some built-in way of doing what I need here.

Thanks,

Andrew Hicks