Solving for displacement and lagrange multiplier with 2D body

I’m trying to find the solution of displacements (u) and Lagrange multipliers (p) in 2D rectangle body with fixed left edge (u = 0) and traction specified on the right edge. All other edges are traction-free. I have never dealt with Lagrange multipliers and I’m struggling to set up the problem in Fenics.

Here’s the Euler-Lagrnage equations I’ve found:
∇∙σ=0 (1)
∇∙u=0 (2)
σ=2μϵ+pI (3)

and by definition of strain:
ϵ=1/2 (∇u+(∇u)^T ) (4)

To solve for u and p, I did the following:

Define function space

V1 = VectorFunctionSpace(mesh, ‘P’, 1)
V2 = FunctionSpace(mesh, ‘R’, 0)

Assign left/right boundaries

def left(x,on_boundary):
return near(x[0],0) and on_boundary
def right(x,on_boundary):
return near(x[0],3) and on_boundary

Mark the boundary on the right

boundary_subdomains = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
AutoSubDomain(left).mark(boundary_subdomains, 1)
AutoSubDomain(right).mark(boundary_subdomains, 2) # This is the boundary where the traction is applied
dss = ds(subdomain_data=boundary_subdomains)

BC: u_x = u_y = 0 at the left

bc = DirichletBC(V,Constant((0,0)),left)

Define strain

def epsilon(u):
return 0.5*(grad(u) + grad(u).T)

Define stress

def sigma(u,p):
return 2muepsilon(u) + p*Identity(d)

Define variational problem

u = TrialFunction(V1)
p = TrialFunction(V2)
d = u.geometric_dimension() # space dimension of u
v = TestFunction(V1)
q = TestFunction(V2)
f = Constant((0,0)) # External force = 0
T = Constant((0,0)) # Traction free BC

I’m planning to add Petrov-Galerkin term in the weak form of the equation (Why I introduced Lagrange multiplier here). Am I going the right way so far? How should I derive the weak form to solve for u and p?