# How to interpolate numpy array into VectorFunctionSpace?

**URL:** <https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002>\
**Category:** dolfinx\
**Created:** [August 11, 2022, 5:45pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002 "2022-08-11T17:45:34Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![mtragoza](https://avatars.discourse-cdn.com/v4/letter/m/e95f7d/32.png) [@mtragoza](https://fenicsproject.discourse.group/u/mtragoza)\
**Post date:** [August 11, 2022, 5:45pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/1 "2022-08-11T17:45:35Z")

</div>

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:

```auto
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(
    comm=MPI.COMM_WORLD,
    points=[[-1, -1], [1, 1]],
    n=[n_nodes - 1, n_nodes - 1],
    cell_type=dolfinx.mesh.CellType.triangle,
    diagonal=dolfinx.mesh.DiagonalType.right_left
)
V = dolfinx.fem.VectorFunctionSpace(mesh, ('Lagrange', 1), dim=2)

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

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')
axes[0,0].set_title('y')
axes[0,1].set_title('y_h')
axes[0,2].set_title('y_data')
axes[0,3].set_title('y_data_h')

fig.tight_layout()
fig.show()

```

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

---

<div class="post-metadata">

**Author:** ![mtragoza](https://avatars.discourse-cdn.com/v4/letter/m/e95f7d/32.png) [@mtragoza](https://fenicsproject.discourse.group/u/mtragoza)\
**Post date:** [August 11, 2022, 5:49pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/2 "2022-08-11T17:49:51Z")

</div>

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

 ![download](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fe88bfaac5934f2c06e7ba6fb8b2c7620afd1c17.png)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 11, 2022, 7:27pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/3 "2022-08-11T19:27:38Z")

</div>

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

 ![picture](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8b1269a8f81057e46a4857417df9323b1deabb07.png)

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 11, 2022, 7:32pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/4 "2022-08-11T19:32:21Z")

</div>

Is it related to [this issue](https://fenicsproject.discourse.group/t/extend-submesh-function-to-full-mesh/8006/6)?

---

<div class="post-metadata">

**Author:** ![mtragoza](https://avatars.discourse-cdn.com/v4/letter/m/e95f7d/32.png) [@mtragoza](https://fenicsproject.discourse.group/u/mtragoza)\
**Post date:** [August 11, 2022, 7:42pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/5 "2022-08-11T19:42:29Z")

</div>

@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:

```auto
# 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!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 11, 2022, 8:24pm UTC](https://fenicsproject.discourse.group/t/how-to-interpolate-numpy-array-into-vectorfunctionspace/9002/6 "2022-08-11T20:24:18Z")

</div>

> [@mtragoza](#):
>
> @dokken I installed dolfinx through conda and I am using version 0.4.1.

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.