Correctly updating solutions and imposing Lagrange multiplier constraints in a Mixed Formulation

Question

Hello,

In response to a comment to my previous question on operator splitting on mixed function spaces (link here) - I tried to break down the issue to further detail. I realized that I was trying to impose a distributed lagrange multiplier to impose a constraint over a domain - and doing that within the mixed formulation could be where I am getting stuck. I have attached here in a simpler code with Stokes flow in 2D around a circular patch where a specific fixed velocity V0 is to be imposed using Lagrange multipliers. I have placed here the code from after mesh generation to make the code sample shorter

# - derive a class for the side-walls for boundary conditions:
class SideWalls(SubDomain):
    def inside(self, x, on_boundary):
        checkVar = on_boundary and ( near(x[0],0.0) or near(x[0], boxW) )
        return checkVar

# - derive a class for the top-wall for boundary conditions:
class TopWall(SubDomain):
    def inside(self, x, on_boundary):
        checkVar = on_boundary and near(x[1], boxH)
        return checkVar

# - derive a class for the bottom-wall for boundary conditions:
class BottomWall(SubDomain):
    def inside(self, x, on_boundary):
        checkVar = on_boundary and near(x[1], 0.0)
        return checkVar

# - derive a class for the circle:
class CircleBall(SubDomain):
    def inside(self, x, on_boundary):
        checkVar = insideParticle(Xp_old, Yp_old, radP, x[0], x[1])
        return checkVar

Xp_old = xp0
Yp_old = yp0
Vx_old = 0.1
Vy_old = 0.1

# Define function spaces
V   = VectorFunctionSpace(mesh, "CG", 2)  # 2D velocity
Q   = FunctionSpace(mesh, "CG", 1)        # pressure
Lp  = VectorFunctionSpace(mesh, "CG", 1)  # Lagrange multiplier field

solutionSpace   = MixedFunctionSpace([V, Q, Lp])      # combined mixed function space
(U, P, L) = TrialFunctions(solutionSpace)
(W, Q, M) = TestFunctions(solutionSpace)

# the specified velocity over the domain
V0  = Function(V)                           
V0  = interpolate(Constant((0.01, 0.01)), V)

# - we will now mark the mesh boundaries:
boundaries = FacetFunction("size_t", mesh)

boundaries.set_all(0)

sidewalls = SideWalls()
sidewalls.mark(boundaries, 1)

topwall = TopWall() 
topwall.mark(boundaries, 2)

bottomwall = BottomWall() 
bottomwall.mark(boundaries, 3)

# - then we will mark the mesh domain interiors:
domains = CellFunction("size_t", mesh)
domains.set_all(0)
circleBall = CircleBall()
circleBall = circleBall.mark(domains, 1)

# - create the appropriate measures for the domains for 
# - assembling the integrals
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]

# - construct all Dirichlet boundary conditions:
noslip  = Constant((0.0,0.0))
zero    = Constant(0.0)
gravity = Constant((0.0, -gValue))
delP    = radP

# - all dirichlet boundary conditions
bc0 = DirichletBC(solutionSpace.sub(0), noslip, sidewalls)
bc1 = DirichletBC(solutionSpace.sub(0), noslip, topwall)
bc2 = DirichletBC(solutionSpace.sub(0), noslip, bottomwall)
bcD = [bc0, bc1, bc2]

# - define the variational forms for the stokes problem 
# - with a distributed lagrange multiplier inside the circle
B1 = Nu*inner(nabla_grad(U), nabla_grad(W))*dx \
        - P*nabla_div(W)*dx \
        + rhoF*dot(gravity, W)*dx \
        - dot(L, W)*dx(1) \
        - (delP**2)*inner(nabla_grad(L), nabla_grad(W))*dx(1)

B2 = Q*nabla_div(U - V0)*dx

B3 = dot(M, U - V0)*dx(1) \
        + (delP**2)*inner(nabla_grad(M), nabla_grad(U - V0))*dx(1)

variationalForm = B1 + B2 + B3

# - now solve the linear variational problem
a = rhs(variationalForm)
b = lhs(variationalForm)

solution = Function(solutionSpace)

problem = LinearVariationalProblem(a, b, solution, bcD)
solver  = LinearVariationalSolver(problem)

solver.solve()

(Usol, Psol, Lsol) = solution.split(deepcopy=True)

If i run this, I get the following error message:

*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  3b6582dfb45139c906c13b9ad57395632a2090f4
*** -------------------------------------------------------------------------

Any suggestions, tips, or advice will be very helpful.

Thanks