How to interpolate numpy array into VectorFunctionSpace?

I am trying to solve a parameter identification problem using dolfinx in which I have access to measured 2D displacement data on a uniform rectangular grid (i.e. a numpy array of shape (n_x, n_y, 2)) and I want to solve a PDE to recover material parameters. However, I am stuck on how to represent the displacement data in the FEM basis.

From reading the forums, the recommended approach is to use scipy interpolation to obtain a continuous function that interpolates the array data, then use Function.interpolate to represent the data in the FEM basis. This works fine if I use a scalar FunctionSpace to interpolate single component of the displacement data, but I am getting weird results when I try it with a VectorFunctionSpace on all components at once. Please see the minimally reproducing code below to understand the problem:

from scipy.interpolate import LinearNDInterpolator

# define a vector-valued function f: R^2 \to R^2
#   NOTE that in reality, I do not have an explicit function like this
f = lambda x: np.stack([
    (x[0] + x[1])**2,
    (x[0] - x[1])**2
], axis=0)

def sample_domain(n):
    x = np.linspace(-1, 1, n)
    return np.stack(np.meshgrid(x, x), axis=-1).reshape(-1, 2)

# evaluate f on a set of domain samples
n = 100
x = sample_domain(n)
y = f(x.T).T

# create an FEM mesh and basis
n_nodes = 5 # per dimension
mesh = dolfinx.mesh.create_rectangle(
    points=[[-1, -1], [1, 1]],
    n=[n_nodes - 1, n_nodes - 1],
V = dolfinx.fem.VectorFunctionSpace(mesh, ('Lagrange', 1), dim=2)

# interpolate f into the FEM basis
f_h = dolfinx.fem.Function(V)

def eval_func(f_h, x):
    '''Evaluate a dolfinx function on a set of 2D points.'''
    x = np.concatenate([x, np.zeros((len(x), 1))], axis=1)
    tree = dolfinx.geometry.BoundingBoxTree(mesh, mesh.geometry.dim)
    cells = dolfinx.geometry.compute_collisions(tree, x)
    cells = dolfinx.geometry.compute_colliding_cells(mesh, cells, x)
    cells = [cells.links(i)[0] for i in range(x.shape[0])]
    return f_h.eval(x, cells)

# evaluate the FEM function on the domain
y_h = eval_func(f_h, x)  # this looks fine

# create data set of samples from f
#   NOTE this is what I have in reality that I want to represent in the FEM basis
n_data = 5
x_data = sample_domain(n_data)
y_data = f(x_data.T).T

# create function that interpolates data samples
f_data = LinearNDInterpolator(x_data, y_data)
y_data = f_data(x) # this looks fine

# interpolate f_data into the FEM basis
f_data_h = dolfinx.fem.Function(V)
f_data_h.interpolate(lambda x: f_data(x[:2].T).T)
y_data_h = eval_func(f_data_h, x) # this looks incorrect...

# reshape into images for plotting
y = y.reshape(n, n, 2)
y_h = y_h.reshape(n, n, 2)
y_data = y_data.reshape(n, n, 2)
y_data_h = y_data_h.reshape(n, n, 2)

# plot the components of y, y_h, y_data, and y_data_h
fig, axes = plt.subplots(2, 4, figsize=(8, 4))

axes[0,0].imshow(y[...,0].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[1,0].imshow(y[...,1].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[0,1].imshow(y_h[...,0].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[1,1].imshow(y_h[...,1].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[0,2].imshow(y_data[...,0].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[1,2].imshow(y_data[...,1].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[0,3].imshow(y_data_h[...,0].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])
axes[1,3].imshow(y_data_h[...,1].real, cmap='Greys', vmin=0, vmax=4, extent=[-1, 1, -1, 1])

axes[0,0].set_ylabel('component 0')
axes[1,0].set_ylabel('component 1')


Can someone please help me understand why y_data_h does not look like the other images in the plots above?

Here is what the plot looks like for the above code:

What version of DOLFINx are you using? I cannot reproduce this with dolfinx/dolfinx:nightly docker images.
I got the following plot:

Is it related to this issue?

@dokken I installed dolfinx through conda and I am using version 0.4.1.

@nate It does seem to be related to that. I tried the following modification:

# interpolate f_data into the FEM basis
f_data_h = dolfinx.fem.Function(V)
f_data_h.interpolate(lambda x: np.ascontiguousarray(f_data(x[:2].T).T))
y_data_h = eval_func(f_data_h, x)

And the plots look correct now. Thanks for the help!

This seems like it has been fixed since the 0.4.1 release.

For now you can use @nate’s solution, but keep in mind that the behaviour will change once you update dolfinx.