Action() command not working in the context of mixed spaces

General
1

If A is an operator defined as (\nabla \cdot, \nabla \cdot), where (\cdot, \cdot) is the usual L^2 scalar product. Let u \in V, where V is a standard conforming piecewise linear finite element space and u0 be some interpolant in the space V (i.e. a vector). I want to compute (\nabla u, \nabla v)*u0. In the primal case I am doing this using the action() command (see the following minimal code with results)

from dolphin import *

mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh, 'CG', 1)
u = Function(V)

u0 = interpolate(Expression('x[0]', degree = 2), V)
u = TrialFunction(V)
v = TestFunction(V)

a = dot(grad(u), grad(v))*dx
aV= action(a, u0)
print(assemble(aV).get_local())

output is a linear form as expected:

[-0.125 -0.25   0.    -0.25   0.     0.    -0.25   0.     0.     0.
 -0.125  0.     0.     0.     0.125  0.     0.     0.     0.25   0.
  0.     0.25   0.     0.25   0.125]

However in my problem I am working in mixed spaces and I face the following error when I try to use the action() command

Replacement expressions must have the same shape as what they replace.

when I try running following minimal code:

from dolfin import *

mesh = UnitSquareMesh(4,4)
P1 = FiniteElement('CG', triangle, 1)
P0 = FiniteElement('DG', triangle, 0)
element = MixedElement([P1, P0])
V = FunctionSpace(mesh, element)

v1, v2 = TestFunctions(V)
u, phi = TrialFunctions(V)

u0 = interpolate(Expression('x[0]', degree = 2), V.sub(0).collapse())

a = dot(grad(u), grad(v1))*dx
aV = action(a, u0)

Is there something wrong with my code or does the action() command not work for mixed spaces?
If not then is there another way to compute the residual ?

2

As the error said, you cannot use action (with in turn uses replaces) on a component of a function.

You need to replace it with a function from space (V), which can be achieved with a function assigner:

from dolfin import *

mesh = UnitSquareMesh(4, 4)
P1 = FiniteElement('CG', triangle, 1)
P0 = FiniteElement('CG', triangle, 1)
element = MixedElement([P1, P0])
V = FunctionSpace(mesh, element)

v1, v2 = TestFunctions(V)
u, phi = TrialFunctions(V)

u0 = interpolate(Expression('x[0]', degree=2), V.sub(0).collapse())
assigner = FunctionAssigner(V.sub(0), u0.function_space())
w0 = Function(V)

assigner.assign(w0.sub(0), u0)

a = dot(grad(u), grad(v1))*dx
aV = action(a, w0)
print(assemble(aV).get_local())