TypeError using Dirac Delta Functions on Mixed Function Spaces

[orsonhart]

I am trying to implement a dirac delta type contribution on a Mixed Function Space, following https://www.firedrakeproject.org/demos/full_waveform_inversion.py.html

When I define

 V = VectorFunctionSpace(mesh, "CG", 2)
 u = TrialFunction(V)
 v = TestFunction(V)

I can then use

source_mesh = VertexOnlyMesh(mesh, [[0.5,0.5]])
V_s = VectorFunctionSpace(source_mesh, "DG", 0)
d_s = Function(V_s)
d_s.interpolate(as_vector([4,1]))
source_cofunction = assemble(inner(d_s,TestFunction(V_s)) * dx)
q_s = Cofunction(V.dual()).interpolate(source_cofunction)

to define a linear functional $q_s$, which acts as a delta function, which I use in my variational form.

However, when I instead use a mixed function space representation, e.g.

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q
                  
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)

then I get an error with L = q_s * v, namely TypeError: Action right argument must be either Coefficient or BaseForm. I expect this is because now my test function $v$ doesn't truly live in $V$, but instead lives in some other function space, but I am not sure how to proceed. I get the same problem if I write q_s = Cofunction(W.sub(0).dual()).interpolate(source_cofunction) instead.

Perhaps it is possible to directly recover the function space that my test function $v$ lives in for the mixed formulation, and use the dual space of that when defining $q_s$, but I am not sure how (in the first case, v is of type <class 'firedrake.ufl_expr.Argument'>, but in the second is of type <class 'ufl.tensors.ListTensor'>).

Many thanks.

[orsonhart]

Figured out how to fix this, you have to instead work with a test function on the mixed space $W$ and define the Cofunction on the dual of $W$.
Example code:

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V*Q

f = Function(W)
w = TestFunction(W)

u, p = split(f)
v, q = split(w)

source_mesh = VertexOnlyMesh(mesh, [[0.5,0.5]])
V_s = VectorFunctionSpace(source_mesh, "DG", 0)
Q_s = FunctionSpace(source_mesh, "DG", 0 )
W_s = V_s * Q_s
d_s = Function(W_s)

d_s.sub(0).interpolate(1)
d_s.sub(1).interpolate(0)

source_cofunction = assemble(inner(d_s,TestFunction(W_s)) * dx)
q_s = Cofunction(W.dual()).interpolate(source_cofunction)

and include q_s * w in the functional form