interpolation inside a system matrix doesn't work #5477

The interpolation features let us put a delta function on the right-hand side of a PDE (see the full waveform demo) but we can't have an interpolation operator inside the system matrix itself. The simplest example of where you'd want to do this is sparse data interpolation. The objective functional you want to minimize looks like

J(u) = sum over point cloud of |u(x_k) - observed value at x_k|^2 +
           integral over mesh of |grad(grad(u))|^2 dx

Desired feature: I'd like to be able to form this objective and then call solve(derivative(J, u) == 0, u). See the "minimization" case in the attached script.

The problem: If I is the interpolation operator from the mesh into the point cloud, the offending term in the system matrix looks like adjoint(I) * M * I or just adjoint(I) * I, which breaks form assembly.

Workarounds: I can form and solve the linear system that's equivalent to minimizing J. I was able to figure out how to do this using petsc4py; you assemble I and then create adjoint(I) * I as a PETSc matrix object. See the "petsc4py" case in the script. Claude also found another way to do this using only Firedrake, but it involves some funny business in the matrix assembly. See the "firedrake" case in the script. Both of these are fine for now but I'd still like to be able to set this up as a minimization problem and not as a linear system.

You could do this with pyadjoint and ROL or TAO but here we don't get the second derivative of J, only its action on a vector. I don't want to worry about iterative methods if I can avoid it.

Script is below. You can run it as python test.py --method minimization to get the failing case of what I'd like to do, --method petsc4py for the low-level way, and --method firedrake for the dodgy one that Claude came up with. The two working methods give the same result.

test.py