Interpolator as matrix

Question

Is it possible to represent the interpolation of a solution to pre-defined points on the mesh as a matrix operator?

Answer 1

Yes, it is possible, but it is not trivial. I have a solution that I use for probing a Function at certain locations any number of times. It is much more efficient than using a regular eval, because I have precomputed some matrices and perform the new probe with a matrix vector product. An example is here:

from dolfin import *
from numpy import array
from fenicstools import Probes

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)

# Create some random Function that you want to Probe
u0 = interpolate(Expression('x[0]'), V)
x = array([[0.5, 0.5], [0.4, 0.4], [0.3, 0.3]])

# Create the Probes object
p = Probes(x.flatten(), V)

# Probe three times
p(u0)
p(u0)
p(u0)

# print solution
print p.array()

[[ 0.5 0.5 0.5]
[ 0.4 0.4 0.4]
[ 0.3 0.3 0.3]]

You would have to install fenicstools to make it work.

Answer 2

I faced the same problem. I made an interpolation matrix in the following way:

from fenics import *
import numpy as np
import scipy.sparse as sps

def make_interpolation_matrix(xs, V):
    nx, dim = xs.shape
    mesh = V.mesh()
    coords = mesh.coordinates()
    cells = mesh.cells()
    dolfin_element = V.dolfin_element()
    dofmap = V.dofmap()
    bbt = mesh.bounding_box_tree() 
    sdim = dolfin_element.space_dimension()
    v = np.zeros(sdim)
    rows = np.zeros(nx*sdim, dtype='int')
    cols = np.zeros(nx*sdim, dtype='int')  
    vals = np.zeros(nx*sdim)
    for k in range(nx):
        # Loop over all interpolation points
        x = xs[k,:]
        p = Point(x[0],x[1])
        # Find cell for the point
        cell_id = bbt.compute_first_entity_collision(p)
        # Vertex coordinates for the cell
        xvert = coords[cells[cell_id,:],:]
        # Evaluate the basis functions for the cell at x
        dolfin_element.evaluate_basis_all(v,x,xvert,cell_id)
        jj = np.arange(sdim*k,sdim*(k+1))
        rows[jj] = k
        # Find the dofs for the cell
        cols[jj] = dofmap.cell_dofs(cell_id)
        vals[jj] = v

    ij = np.concatenate((np.array([rows]), np.array([cols])),axis=0)
    M = sps.csr_matrix((vals, ij), shape=(nx,V.dim()))
    return M

mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh,'CG', 2)
xs = np.random.rand(20,mesh.geometry().dim())
M = make_interpolation_matrix(xs, V)

# Testing the matrix
f = interpolate(Expression('x[0]'),V)
fx = M*f.vector().array()
print fx-xs[:,0] # should be zero