FEniCS中基函数的计算

在统计反问题的计算中，若我们利用有限元方法将正问题（PDE）离散，则我们在计算Pointwise variance field时，会需要计算有限元基函数在特定点处的值。

例如如下参考文献3.7小节中就推到了Pointwise variance field的计算：
T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler, A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with applications to global seismic inversion, SIAM Journal on Scientific Computing, 35(6), A2494-A2523, 2013.

具体来讲，在有限元方法中函数可以表示成如下形式：

计算函数在某点处的值实质上是计算如下两个向量的内积：

下面程序中的函数 make_interpolation_matrix(xs, V) 中两个参数分别是点和函数空间，例如：
xs = np.array([[0.5,0.5], [0.6, 0.6]]) V = FunctionSpace(mesh, 'Lagrange', 2)
函数的输出是一个矩阵M,

因此，函数在点处的值可以如下计算

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)
        v = dolfin_element.evaluate_basis_all(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(3, 3)
V = FunctionSpace(mesh, 'Lagrange', 2)
xs = np.array([[0.5, 0.5]])#np.random.rand(1, mesh.geometry().dim())
M = make_interpolation_matrix(xs, V)

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