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

dolfinx
,
4

What you are searching for is how to create a point source.
Following is code for getting the basis functions evaluated at a point in the cell:

# Author: Jørgen S. Dokken
#
# Create a point source

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl


def compute_cell_contributions(V, points):
    # 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)
    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))

        # 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)
    return cells, basis_values


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

if mesh.comm.rank == 0:
    points = np.array([[0.1, 0.2, 0],
                      [0.91, 0.93, 0]], dtype=mesh.geometry.x.dtype)
else:
    points = np.zeros((0, 3), dtype=mesh.geometry.x.dtype)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

cells, basis_values = compute_cell_contributions(V, points)
print(cells, basis_values)

This gives the list of cells on your process, and the corresponding basis values (per row), that you in turn can insert into your matrix/vector.

3 Likes