Building the matrix associated with pointwise evaluation of basis functions

[douglas-brinkerhoff]

I am in need of the sparse matrix associated with the evaluation of a firedrake function at a collection of arbitrary coordinates. In principle, I can build the matrix that I'm looking for column-wise with the following:

import firedrake as df
import numpy as np

mesh = df.UnitSquareMesh(25,25)

V_cg = df.FunctionSpace(mesh,"CG",1)
u = df.Function(V_cg)

interpolation_points = np.random.rand(100,2)

basis_values = []

for i in range(V_cg.dim()):
    u.dat.data[:] = 0.0
    u.dat.data[i] = 1.0
    basis_values.append(u.at(interpolation_points))

interpolation_operator = np.column_stack(basis_values)  ## This is the matrix I'm after
g = df.interpolate(df.cos(df.SpatialCoordinate(mesh)[0]*2*np.pi),V_cg)
g_at_points = interpolation_operator @ g.dat.data[:]

Obviously, this is horribly inefficient, and doesn't respect the sparsity of the problem. I also recognize that I can get the product of this matrix with a vector via the .at() or VertexOnlyMesh functionality - however I need the matrix itself (and its transpose) in order to solve a regularized least-squares problem for the function coefficients given some pointwise observations.

Is there a clever way to achieve this?

Thanks very much in advance to everyone for their insight.

[dham]

Hi Doug,

There is no way to do this at present. Leo Collins is currently working on this but he's just spinning up as a PhD student, so I don't expect immediate results.

Could you possibly use a matrix-free solver to solve the problem you are trying to solve?

[douglas-brinkerhoff]

Thanks, David.

A matrix-free solver would be great, but I'm solving the normal equations and am stumped as to how to adapt the built-in interpolation methods to evaluate the matvec product A^T v (where A is the operator described previously) - but I might be missing some obvious functionality to compute such a product. For what it's worth, replacing the .at() operation in the loop above with interpolation to a VertexOnlyMesh makes the code much faster, and it's straightforward enough to drop the zeros when putting the resulting columns into a sparse matrix. Obviously not an ideal solution, but one that works well enough for the size of problem that I'm working on.

On a related note, any thoughts on future functionality for the pointwise evaluation of the Mardal-Tai-Winther element?