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

dolfinx
,
6

Thanks again for your codes. Based on your code, I rewrite the codes as follows, which seems work well.

import dolfinx
from dolfinx import fem, la, geometry
from dolfinx.fem import assemble_matrix, assemble_vector
import ufl
from mpi4py import MPI


def construct_measure_matrix(V, points):
    '''
    Input:
        V: function space produced by dolfinx.fem.FunctionSpace
        points: input points, e.g.,
                points = [[x1_1, x1_2, x1_3], [x2_1, x2_2, 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

    Remark:
    The following code is based on the codes provided by Jorgen Dokken at the website:
    https://fenicsproject.discourse.group/t/how-to-construct-a-sparse-matrix-used-for-calculating-function-values-at-multiple-points/13081/4
    '''
    nx, dim = points.shape
    assert dim == 3, "The points should be shape Nx3 (N is the number of points)"
    mesh = V.mesh
    sdim = V.element.space_dimension
    # Determine what process owns a point and what cells it lies within
    # _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
    #     mesh._cpp_object, points, 1e-6)
    ## The extral parameter 1e-6 seems not allowed
    _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points)
    owning_points = np.asarray(owning_points).reshape(-1, 3)

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmaps[0]
    ref_x = np.zeros((len(cells), mesh.geometry.dim),
                     dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    # Create expression evaluating a trial function (i.e. just the basis function)
    u = ufl.TrialFunction(V)
    num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
    if len(cells) > 0:
        # NOTE: Expression lives on only this communicator rank
        expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
        # values = expr.eval(mesh, np.asarray(cells))
        ## the command np.asarry(cells) seems leading to some errors
        values = expr.eval(mesh, cells)

        # Strip out basis function values per cell
        basis_values = values[:num_dofs:num_dofs*len(cells)]
    else:
        basis_values = np.zeros(
            (0, num_dofs), dtype=dolfinx.default_scalar_type)

    rows = np.zeros(nx*sdim, dtype='int')
    cols = np.zeros(nx*sdim, dtype='int')
    vals = np.zeros(nx*sdim)
    for k in range(nx):
        jj = np.arange(sdim*k, sdim*(k+1))
        rows[jj] = k
        # Find the dofs for the cell
        cols[jj] = V.dofmap.cell_dofs(cells[k])
        vals[jj] = basis_values[0][sdim*k:sdim*(k+1)]

    ij = np.concatenate((np.array([rows]), np.array([cols])), axis=0)
    S = sp.csr_matrix((vals, ij), shape=(nx, V.tabulate_dof_coordinates().shape[0]))
    return S
1 Like