Sorry for not posting my current parts of the code. I am not an expert for programming finite element methods, just use the functionals of FEniCSx to solve some inverse problems. FEniCSx provides a quick and convenient tool to test Bayesian or optimization inverse algorithms without too much attention on solving PDEs. The following code is my temporary solution, but it must be not efficient enough.

```
import numpy as np
import scipy
import scipy.sparse as sp
from dolfinx import fem, la, geometry
from dolfinx.fem import assemble_matrix, assemble_vector
import ufl
from mpi4py import MPI
def construct_measure_matrix(points, V):
'''
Input:
V: function space produced by dolfinx.fem.FunctionSpace
points: input points, e.g.,
points = [[x1_1,x2_1], [x1_2,x2_2], [x1_3,x2_3]] stands for
there are two measured points with coordinates:
x1=(x1_1,x1_2,_x1_3), x2=(x2_1,x2_2,x2_3)
Output:
S: a scipy.sparse.csr_matrix
'''
assert points.shape[0] == 3, "The points should be shape 3xN"
domain = V.mesh
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
ut = fem.Function(V)
S = sp.csr_matrix((len(points_on_proc), len(ut.x.array)))
for i, _ in enumerate(ut.x.array):
ut.x.array[:] = 0.0
ut.x.array[i] = 1.0
S[:, i] = ut.eval(points_on_proc, cells).squeeze()
return S
```