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(
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?
