How to construct a sparse matrix used for calculating function values at multiple points

dolfinx
,
3

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