How to specify Lagrange Multiplier on a boundary?

My question is closely related to these topics on the old FEniCS Q&A forum (especially the first one, which has more details).

Related Post 1
Related Post 2

Essentially, I am trying to define a Lagrange multiplier on only the top surface of a cube in a mixed formulation. In the old topic linked above, it was suggested that the Lagrange multiplier space could be constrained by applying DirichletBC on redundant degrees of freedom.

I have an initial condition set (similarly to the initial condition in the Cahn-Hilliard Demo) which initializes the Lagrange multiplier to zero over the whole domain. Furthermore, as suggested in the above Q&A post, a DirichletBC sets the Lagrange multiplier to zero over the whole domain except the top surface (Code snippet included below).

So far, the solver simply diverges after 0 iterations. Specifying the method for DirichletBC seems to have no effect. I would really appreciate it if anyone has an idea on the progress of this topic since 2013-2014 or effective workarounds. I can provide more code if required. This runs on dolfin version 2017.2.0

# Class for marking everything on the domain except the top surface
class bulkWOsurf(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 1.0 – tol

# Top surface
top = CompiledSubDomain("near(x[1], side) && on_boundary", side = 1.0)

facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
facets.set_all(0)

# Top exterior facets are marked as subdomain 1, using 'top' boundary
top.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

# Set lambda to 0 in domain
lmbda_val = Constant(0.0)
t1 = DirichletBC(V.sub(2), lmbda_val, bulkWOsurf, method="pointwise")

In case anyone is interested, this can be done using the package Multiphenics.