Integral constraint

I am trying to solve the PDE \Delta u= f with integral constraint \int_{\Omega}u\;dx=c_0. Following the example posted in fenics demos, I used Mixed Function space. But I am getting the following error
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ()
Can anyone point out what is wrong with my variatioanl formulation below? Thanks

from dolfin import *

mesh = UnitSquareMesh(64, 64)

V = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement("R", mesh.ufl_cell(), 0)
W_elem = MixedElement([V, R])
W = FunctionSpace(mesh, W_elem)

c0 = Constant(0.4)

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v)) * dx + c * v * dx + (u - c0) * d * dx
L = f * v * dx

# Compute solution
w = Function(W)
solve(a == L, w)
(u, c) = w.split()

plot(u, title='Final solution')
filename = 'demo.png'
plt.savefig(filename)
plt.close()

Hi, you should pass the -c0*d*dx in the linear form

I see now. You meant like this
a = inner(grad(u), grad(v)) * dx + c * v * dx + u * d * dx
L = f * v * dx - c0 * d * dx
This seems to work with no problem. But if I am passing to the linear form, I should pass +c0 * d * dx instead of -c0 * d * dx``

of course,
note that you could also have written form = inner(grad(u), grad(v)) * dx + c * v * dx + (u - c0) * d * dx - f * v * dx and use lhs(form) == rhs(form)

Thank you. It works now