Strange behavior for multiple Dirichlet boundary conditions in mixed dimensional branch?

General
1

In this example code, I am trying to solve two copies of the Poisson equation, but using the MixedFunctionSpace() function of the mixed dimensional branch by @cdaversin. I set up a MixedFunctionSpace() of two spaces, each of which will contain a solution to the Poisson equation, but with different boundary conditions.

boxlength = 1
# Define mesh and function spaces
mesh = RectangleMesh(Point(0, 0), Point(boxlength, boxlength),
                     32, 32)
V = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace(V, V)

# Create classes for boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], boxlength)
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], boxlength)

left = Left()
right = Right()
top = Top()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sLef = 1
left.mark(boundaries, sLef)
sRig = 2
right.mark(boundaries, sRig)
sTop = 3
top.mark(boundaries, sTop)
sBot = 4
bottom.mark(boundaries, sBot)

# Dirichlet boundary conditions
boundfunct = Expression('x[1]*x[1]', degree = 1)
bcsA = [DirichletBC(W.sub_space(0), 1.0, boundaries, sTop)]
bcsB = [DirichletBC(W.sub_space(1), boundfunct, boundaries, sBot)]
bcs = bcsA + bcsB

# Define model parameters
g_Lef = Expression('0', degree = 1)
g_Rig = Expression('0', degree = 1)
g_Top = Expression('0', degree = 1)
g_Bot = Expression('0', degree = 1)

# Define measures for boundary terms
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define variational form
uA, uB = TrialFunctions(W)
vA, vB = TestFunctions(W)

aA = inner(grad(uA), grad(vA))*dx
LA = - g_Lef*vA*ds(sLef) - g_Rig*vA*ds(sRig) - g_Bot*vA*ds(sBot)

aB = inner(grad(uB), grad(vB))*dx
LB = - g_Lef*vB*ds(sLef) - g_Rig*vB*ds(sRig) - g_Top*vB*ds(sTop)
a = aA + aB
L = LA + LB

sol = Function(W)   
solve(a == L, sol, bcs)

In this code, I attempt to apply a constant uA = 1 Dirichlet boundary condition to the uA solution, and the boundfunct boundary condition to the uB function. However, this code applies the constant boundary condition to both uA and uB.

How do I apply a different boundary condition to each subspace of a MixedFunctionSpace?