Building operator to extract solution at points, discrepancy?

dolfinx
1

Hello all,

Have used legacy fenics for a long time and am finally getting around to porting my utilities over. The team has done a great job!

Particularly, I’m trying to build an operator (matrix) which acts on an element of the finite element space and extracts its values at a certain set of locations. Note, this is different from the evaluation at those locations. I need the operator as later I’ll need its transpose.

Following the excellent tutorial written here:

I’ve figured out how to obtain the action of this operator. Then, knowing what a tabulation is, I’ve written the following (serial) MWE to back out the operator corresponding to what I want.

from itertools import product

import dolfinx
from dolfinx import fem, mesh
from dolfinx.fem import functionspace
from mpi4py import MPI
import numpy as np
from scipy.sparse import csr_matrix

domain = mesh.create_unit_square(MPI.COMM_WORLD, 16, 16)
V = functionspace(domain, ("Lagrange", 1))


def u_py(x):
    return 1 + x[0] ** 2 + 2 * x[1] ** 2


u = fem.Function(V)
u.interpolate(u_py)

sensor_locations = list(product(np.linspace(0.1, 0.9, 4), np.linspace(0.1, 0.9, 4)))
true_vals = [u_py(x) for x in sensor_locations]

# Goal, obtain operator B such that B(u) = [u(p) for p in sensor_locations]
# Alternatively, obtain ufl form that, when assembled, produced B(u)

pts_3d = np.array([[p[0], p[1], 0] for p in sensor_locations])

bb_tree = dolfinx.geometry.bb_tree(domain, domain.topology.dim)
potential_colliding_cells = dolfinx.geometry.compute_collisions_points(bb_tree, pts_3d)
colliding_cells = dolfinx.geometry.compute_colliding_cells(
    domain, potential_colliding_cells, pts_3d
)
cells = []
for i, point in enumerate(pts_3d):
    if len(colliding_cells.links(i)) > 0:
        cells.append(colliding_cells.links(i)[0])

cells = np.array(cells, dtype=np.int32)

cell_global_dofs = [V.dofmap.cell_dofs(idx) for idx in cells]
tab = V.element.basix_element.tabulate(0, np.array(sensor_locations))

row_idx = []
col_idx = []
data = []
for i, target in enumerate(pts_3d):
    for local_dof, global_dof in enumerate(cell_global_dofs[i]):
        row_idx.append(i)
        col_idx.append(global_dof)
        data.append(tab[0][i][local_dof][0])

Q = csr_matrix(
    (data, (row_idx, col_idx)),
    shape=(pts_3d.shape[0], V.dofmap.index_map.size_global),
)
u_pts = Q @ u.x.array

# u_pts \not\approx true_vals !!! But it's close?

Basically,

  • Obtain the cells which contain the points
  • Back out the corresponding global degrees of freedom
  • Compute a tabulation
  • Build the operator corresponding to what I understand the interpolation operator should be doing

However, this does not coincide with what I expect the true value of interpolation to be! I expect a small amount of error (~1e-2) due to the discretization, but this error is much larger! Moreover, u.eval(pts_3d, cells) obtains what I would expect, hence I assume that there’s an issue with my implementation.

Perhaps an assumption I’ve made regarding the support of the Lagrangian elements? Or is there an ordering issue? Anyone who knows how eval works would care to help me out?

Anything is appreciated!

Best,
Abhi

2

My apologies, I belatedly found this question which I believe resolves my issue.

3

I would be great if you could post your solution/point out what errors exist in the code above.

4

Good point, soon I’ll clean up the MWE with the fix.

5

This is the functional code. There might be some optimizations possible here. Also, I haven’t tested the parallel functionality (I suspect that there’s still work on that end). The docstring cites the relevant sources.

from itertools import product
import dolfinx
from dolfinx import fem, mesh
from dolfinx.fem import functionspace
from mpi4py import MPI
import numpy as np
import scipy as sp
import ufl


domain = mesh.create_unit_square(MPI.COMM_WORLD, 16, 16)
V = functionspace(domain, ("Lagrange", 1))


def u_py(x):
    return 1 + x[0] ** 2 + 2 * x[1] ** 2


u = fem.Function(V)
u.interpolate(u_py)

sensor_locations = list(product(np.linspace(0.1, 0.9, 4), np.linspace(0.1, 0.9, 4)))
true_vals = [u_py(x) for x in sensor_locations]

# Goal, obtain operator B such that B(u) = [u(p) for p in sensor_locations]


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 and Junxiong Jia at:
        https://fenicsproject.discourse.group/t/how-to-construct-a-sparse-matrix-used-for-calculating-function-values-at-multiple-points/13081/4
        and further modified due to API changes as per the discussion at
        https://fenicsproject.discourse.group/t/point-sources-redux/13496/8
    """
    if isinstance(points, list):
        points = np.array(points)
    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
    point_ownership_data = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6
    )
    owning_points = np.asarray(point_ownership_data.dest_points).reshape(-1, 3)
    cells = point_ownership_data.dest_cells

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    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, 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.sparse.csr_matrix(
        (vals, ij), shape=(nx, V.tabulate_dof_coordinates().shape[0])
    )
    return S


B = construct_measure_matrix(V, [[p[0], p[1], 0] for p in sensor_locations])
computed_vals = B @ u.x.array
for y_true, y_computed in zip(computed_vals, true_vals):
    assert np.abs(y_true - y_computed) < 1e-2

I think with my original code there was a misunderstanding with how tabulate works. But I stopped thinking about it after this started functioning.

1 Like