Imposing bcs using Lagrange Multiplier

mathematics
1

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!

2

Update:

I figured out that I was defining the ‘real space’ on the whole domain(which could be wrong), so I should only define it on the boundary that I want to restrict.

So I did :

submesh = MeshView.create(mf, 1) #left boundary mesh
Lmp = FunctionSpace(submesh, 'R', 0) # space of real numbers on submesh

Also, I defined the ds and dx measures for different subspaces like this:

dx = Measure("dx", domain=W.sub_space(0))
dV = Measure("dx", domain=W.sub_space(1))

ds = Measure("ds", domain=W.sub_space(0))
dL = Measure("ds", domain=W.sub_space(1))

so the variational form becomes:

#Variational form
a = kappa*inner(u,v)*dx - p*div(v)*dx - q*div(u)*dx 
a+= lmp_test[0]*inner(u,n)*dL(1) + lmp_trial[0]*inner(v,n)*dL(1)

L = -P_i*inner(n,v)*ds(2) - P_o*inner(n,v)*ds(3)

usol = Function(W)
A = assemble(a)
b = assemble(L)

But now it gives the error:
ValueError: too many values to unpack (expected 1)

1 Like