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