Imposing zero normal displacement through lagrange multiplier

In the following code, with an elasticity equation, i am trying to impose a zero normal displacement on a dirichlet boundary using lagrange multiplier. Can someone let me know what i should fix, probably in the mixed bilinear form (a_N in my code) to achieve this?

from dolfin import *

# dimensions and material parameters
L = 1; W = 0.2
mu = 1; lambda_ = 1.25

# mesh generation
nx = 5; ny = 2; nz = 2
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(L, W, W), nx, ny, nz)

class Dirichlet(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Neumann(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

# markers for lagrange multiplier's mesh
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Dirichlet().mark(facets, 1)

# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)

# two meshes
mesh_dmin = MeshView.create(domain, 0)
mesh_bdry = MeshView.create(facets, 1)

Velem = VectorElement("CG", mesh_dmin.ufl_cell(), 1, dim=3)
U = FunctionSpace(mesh_dmin, Velem)

Qelem = FiniteElement("CG", mesh_bdry.ufl_cell(), 1)
Q = FunctionSpace(mesh_bdry, Qelem)

# mixed space
W = MixedFunctionSpace(U, Q)

# markers for Dirichlet and Neumann BCs
facet_bdry = MeshFunction("size_t", mesh_dmin, mesh_dmin.geometry().dim()-1, 0)
Dirichlet().mark(facet_bdry, 1)
Neumann().mark(facet_bdry, 2)

# trial and test functions
(u, q) = TrialFunctions(W)
(u_t, q_t) = TestFunctions(W)

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return 2 * mu * epsilon(u) + lambda_ * tr(epsilon(u)) * Identity(len(u))

# measures
dx_  = Measure("dx", domain=W.sub_space(0).mesh())
ds_d  = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=facet_bdry, subdomain_id=1)
ds_n  = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=facet_bdry, subdomain_id=2)
# measure for LM mesh
dx_i = Measure("dx", domain=W.sub_space(1).mesh())

# first bilinear form w/o LM
a_M = inner(sigma(u), epsilon(u_t)) * dx_

g_L = Constant((1e0, 0, 0))
L_M = inner(g_L, u_t) * ds_n

w_ = Function(W)
(u_, q_) = w_.split()

## to solve elasticity problem with LM for imposing zero normal displacement on zero boundary

# dirichlet BCs
bc_Ly = DirichletBC(W.sub_space(0).sub(1), Constant(0), facet_bdry, 1)
bc_Lz = DirichletBC(W.sub_space(0).sub(2), Constant(0), facet_bdry, 1)

bcs = [bc_Ly, bc_Lz]

# mixed bilinear form --> !NOT WORKING!
a_N = q_t * inner(u, FacetNormal(mesh_bdry)) * dx_i + q * inner(u_t, FacetNormal(mesh_bdry)) * dx_i
a = a_M + a_N

problem = LinearVariationalProblem(a, L_M, w_, bcs)
nsolver = LinearVariationalSolver(problem)
nsolver.solve()