Using FEniCS to generate sampling matrix

For solving inverse problems, we usually need to sample a function at some discrete points. Considering the following problem

with . Here represents a partial differential equation or other forward operator, is a sampling operator can be defined as

Usually, for linear operator , it can be discretized as a matrix also denoted as . In the following, we provide a strategy to generate the sampling matrix . For finite element method, a function can be expanded by a series of basis functions

Then we have

Hence, the sampling matrix can be obtained by the basis functions evaluated at sampling points. Notice that

so if we take and for , we find

Based on this analysis, we can generate a function as follow:

1. generate a function u in FEniCS;
2. take each component of u.vector() to be 1 and other components to be 0, and evaluate the function u at the required point.

In the following, we provide a program. The parameter V can be taken as

V = fe.FunctionSpace(mesh, 'P', 1)

The parameter points can be taken as

points = [(0.5,0.5), (0.5, 0.6)]

The program:

import numpy as np
import fenics as fe

def measureMatrix(V, points):
    len_points = len(points)
    u = fe.Function(V)
    u.vector()[:] = 0
    u_vec = u.vector()[:]
    len_u = len(u_vec)
    sampleMatrix = np.zeros((len_points, len_u))
    
    for i in range(len_points):
        for j in range(len_u):
            u_vec[j] = 1
            u.vector()[:] = u_vec
            sampleMatrix[(i, j)] = u(points[i])
            u_vec[:] = 0
    
    return sampleMatrix

The following program illustrates the effectiveness of the above program.

import fenics as fe
import numpy as np

mesh = fe.UnitSquareMesh(10, 10)
V = fe.FunctionSpace(mesh, 'P', 1)

u_D = fe.Expression('1+x[0]*x[0]+2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = fe.DirichletBC(V, u_D, boundary)

u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(-6.0)
a = fe.dot(fe.grad(u), fe.grad(v))*fe.dx
L = f*v*fe.dx

u = fe.Function(V)
fe.solve(a == L, u, bc)

points = []
for i in range(100):
    points.append((0.5, i/100.0))
    
S = measureMatrix(V, points)
val1 = S@u.vector()[:]
val2 = [u(i) for i in points]
print(np.linalg.norm(val1-val2))