Oscillations in solution to inverse problem based on mesh resolution

mesh
1

I am using dolfinx to solve an inverse elasticity problem where I have a data set of steady-state displacements and the task is to solve some PDE to recover the stiffness of the material.

I posted here previously because the first problem I ran into was how to interpolate the measured displacement data into the FEM function space. I ended up solving this (I thought) by using scipy.interpolate to create an interpolating function that I could then pass to dolfinx.Function.interpolate.

This seemed like a roundabout solution, as what I really want to do is just set the FEM basis coefficients directly from the data set, but I couldn’t figure out how to do that. One consequence of using the scipy interpolation approach is that I am not confined to using the same mesh resolution for the FEM problem as the resolution of the grid on which the data are defined.

The FEM method seems to be working now, but the current problem is that the solution (the elasticity field) exhibits oscillations that are extremely sensitive to the mesh resolution. For instance, in a 1D problem with an array of 80 complex-valued displacements, there are large oscillations in the solution if I use a mesh with 80 equally spaced nodes, which seemed to be that natural mesh resolution, since it is the data resolution. If I use slightly more or fewer nodes, the oscillations change, but they never completely disappear. See the plot below to understand what I mean.

How can I eliminate these oscillations in the solution?

download
download727×1461 167 KB

2

Here is a summary of the code I am using:

ndim = 1 # number of spatial dimensions
x = data.u.field.spatial_points() # x is an (80, 1) real array
u = data.u.field.values()         # u is an (80, 1) complex array

# discretize the domain on a mesh
mesh = dolfinx.mesh.create_interval(
    comm=MPI.COMM_WORLD,
    points=[x.min(axis=0), x.max(axis=0)],
    nx=n_mesh, # this is the parameter I'm playing with
)

# define the FEM basis function spaces
scalar_space = dolfinx.fem.FunctionSpace(mesh, ('CG', 1))

# setup physical problem
rho = 1000
omega = 2 * np.pi * 80

# function for interpolating wave image into FEM basis
u_interp = scipy.interpolate.interp1d(
    x[:,0], u[:,0], bounds_error=False, fill_value='extrapolate'
)
u_func = dolfinx.fem.Function(scalar_space)

def f(x):
    return np.ascontiguousarray(u_interp(x[:ndim].T).T)

u_func.interpolate(f)

# trial and test functions
mu_func = ufl.TrialFunction(scalar_space)
v_func = ufl.TestFunction(scalar_space)

# define the variational problem
Auv = mu_func * ufl.inner(ufl.grad(u_func), ufl.grad(v_func)) * ufl.dx
bv = self.rho * self.omega**2 * ufl.inner(u_func, v_func) * ufl.dx

problem = dolfinx.fem.petsc.LinearProblem(
    Auv, bv, bcs=[], petsc_options={'ksp_type': 'lsqr', 'pc_type': 'none'}
)
mu_pred_func = problem.solve()