Hello dear FEniCS users,

I am solving a mixed formulation using DG0 and BDM1 mixed function spaces and want to implement `DirichletBC`

using Lagrange Multiplier to relax constraint at boundaries, because otherwise if I implement the constraints using `DirichletBC`

function , then I am getting unsatisfactory results because of strict restriction at the boundary. So the idea is to do something like this (this is the MWE for a unit cube):

```
from dolfin import *
mesh = UnitSquareMesh(16, 16)
class Bottom(SubDomain): #Bottom boundary of unit cube
def inside(self,x,on_boundary):
return x[1] < DOLFIN_EPS
class Top(SubDomain): #Top boundary of unit cube
def inside(xself,x,on_boundary):
return x[1] > 1.0 - DOLFIN_EPS
class Left(SubDomain): #Left boundary of unit cube
def inside(xself,x,on_boundary):
return x[0] < DOLFIN_EPS
mf = MeshFunction("size_t", mesh, 1)
mf.set_all(0)
bottom = Bottom()
top = Top()
left = Left()
bottom.mark(mf, 3)
top.mark(mf, 2)
left.mark(mf, 1)
ds = Measure("ds")[mf]
n = FacetNormal(mesh)
plot(mesh)
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
mixed = MixedElement(BDM,DG)
W_UP = FunctionSpace(mesh, mixed)
Lmp = FunctionSpace(mesh, 'R', 0) # Space of 'real' numbers
W = MixedFunctionSpace(W_UP,*[Lmp]*1) # 1 Lagrange multiplier for the left boundary
#Test and trial functions
(up, *lmp_trial) = TrialFunctions(W)
(u,p) = split(up)
(vq, *lmp_test) = TestFunctions(W)
(v,q) = split(vq)
# Parameters of domain
kappa = 0.4
P_i = 5
P_o = 1
#Variational form
a = kappa*inner(u,v)*dx - p*div(v)*dx - q*div(u)*dx
a+= lmp_test[0]*inner(u,n)*ds(1) + lmp_trial[0]*inner(v,n)*ds(1)
L = -P_i*inner(n,v)*ds(2) - P_o*inner(n,v)*ds(3)
#Solution
usol = Function(W)
#solve(a == L, usol)
A = assemble(a)
```

I get the following error :

ArityMismatch: Adding expressions with non-matching form arguments (‘v_0^1’, ‘v_1^0’) vs (‘v_0^0’, ‘v_1^1’).

I am aware of doing something like:

```
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
Lmp = FiniteElement('R', mesh.ufl_cell(), 0) # Space of 'real' numbers
mixed = MixedElement([BDM, DG, Lmp])
W = FunctionSpace(mesh, mixed)
#Test and trial functions
u, p, lmp = TrialFunctions(W)
# (u,p) = split(up)
v, q, lmp_test = TestFunctions(W)
```

However, if I go for such an approach, then I am ultimately restricting every element in the boundary for imposing condition `inner(u,n)`

and this is what I am trying to avoid it. The goal is to enforce this `inner(u,n)`

for the whole boundary using one Lagrange Multiplier (for the case of unit cube example given above).

Does anyone have any leads or suggestions on how to do it and get rid of this error?

Thanks in advance and cheers!