MixedFunctionSpace with boundary function spaces

Hi, I am trying to solve a non-linear three-dimensional problem coupling elasticity and solvent diffusion. The state variables, defined over a three-dimensional body, are the displacement vector u, the pressure p and the solvent concentration c. To enforce the boundary conditions I have to define additional degrees of freedom (corresponding to a Lagrange multiplier) on the boundary. So, in the end, I have a mixed function space consisting of three spaces (for u, c and p) over the three-dimensional body and one function space over the boundary. When I run the script (see below) I get this error:

Sub elements must live in the same top level domain.
Traceback (most recent call last):
  File "swelling_gel_equilibrium_clean.py", line 14, in <module>
    W = MixedFunctionSpace([V,Q,R,M])
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 502, in __init__
    element = ufl.MixedElement(*[V.ufl_element() for V in spaces])
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/ufl/finiteelement/mixedelement.py", line 64, in __init__
    error("Sub elements must live in the same top level domain.")
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/ufl/log.py", line 154, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Sub elements must live in the same top level domain.

My guess is that I cannot define a mixed function space combining function spaces over the 3D domain and function spaces over the boundary: is it right? Is there a workaround for this problem, in case this feature is not supported?

Thanks a lot! AL

==BEGIN SCRIPT==

from dolfin import *
import math
import ufl

mesh = UnitCubeMesh(1,1,1)

V = VectorFunctionSpace(mesh, 'Lagrange', 2) # displacement
Q = FunctionSpace(mesh,'Lagrange',1) # pressure
R = FunctionSpace(mesh,'Lagrange',2) # bulk concentration

boundary_mesh = BoundaryMesh(mesh, "exterior")

M = FunctionSpace(boundary_mesh,'Lagrange',2) # Lagrange multiplier
W = MixedFunctionSpace([V,Q,R,M])

# Test and trial functions
vqrm = TestFunction(W) # Test function in the mixed space
v,q,r,tcs = split(vqrm)
dupccs = TrialFunction(W) # Trial function in the mixed space (solution increment)

upccs = Function(W) # Function in the mixed space to store the solution
u,p,c,cs = split(upccs) # Function in each subspace to write the functional

# nonlinear effective diffusion coefficient
def g(c):
    return D*((2*chi-1.)*Omega*Jo*c-1.)/pow((1.+Omega*Jo*c),3)

def dmudc(c):
    return -(Rg*T)*((2.*chi-1.)*Omega*Jo-1./c)/pow((1.+Omega*Jo*c),3)

# Parameters
Omega = 6.023e-5
D = 1E-1
Gdry = 4E5
Rg = 8.314
T = 293.0
lamo = 1.001
Jo = lamo**3
cod = (lamo**3-1.0)/Omega
co = cod/Jo
chi = 0.2
po = Gdry/lamo
muo = Rg*T*(math.log(Omega*cod/(1+Omega*cod))+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
+Omega*po
mus = muo
mutop = -10000

# Kinematics
I = Identity(3)
F = I + grad(u)
J = det(F)

# Constitutive equations

S = (Gdry/lamo)*F-p*cofac(F)
mu = Rg*T*(ufl.ln(Omega*Jo*c/(1+Omega*Jo*c)) + 1./(1.+Omega*Jo*c)+ \
         chi/pow((1.+Omega*Jo*c),2))+Omega*p

h = g(c)*grad(c)+(-c*D/(Rg*T))*Omega*grad(p)

dmudp = Omega
dmu = mu-mus

# Boundary conditions

def left_boundary(x, on_boundary): # x = 0
 return on_boundary and abs(x[0]) < DOLFIN_EPS

def back_boundary(x, on_boundary): # y = 0
 return on_boundary and abs(x[1]) < DOLFIN_EPS

def bottom_boundary(x, on_boundary): # z = 0
 return on_boundary and abs(x[2]) < DOLFIN_EPS

bcl = DirichletBC((W.sub(0)).sub(0), Constant(0.), left_boundary) # u = 0 on x = 0
bcb = DirichletBC((W.sub(0)).sub(1), Constant(0.), back_boundary) # v = 0 on y = 0
bcbo = DirichletBC((W.sub(0)).sub(2), Constant(0.), bottom_boundary) # w = 0 on z = 0
bc = [bcl,bcb,bcbo]

# WEAK FORM OF BALANCE EQUATIONS

a = inner(S,grad(v))*dx + q*(J-(1./Jo+Omega*c))*dx + inner(h,grad(r))*dx + inner(dmu,tcs)*ds \
        + inner(dmudc(c)*cs,r)*ds + inner(dmudp*cs,q)*ds

Jac = derivative(a,upccs,dupccs)

problem = NonlinearVariationalProblem(a, upccs, bc, Jac)
solver = NonlinearVariationalSolver(problem)

# Initial guess for unknowns
c = interpolate(Constant(co), R)
p = interpolate(Constant(po), Q)
cs = interpolate(Constant(co),M)

# solve non-linear problem
solver.solve()

# plot solvent concentration
w = Function(V)
w = upccs.split()[2]
plot(w)
interactive()

==END SCRIPT==