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