How to prescribe displacement at the mid-points of a beam element (this can't be defined variationally or by using DirichletBC)?

Linear Algebra
1

I am using the beam element formulation given in ‘Numerical tours of FEniCS’ to solve a simply supported beam with a prescribed displacement at the mid-point. However, the solution is coming out wrong; the deflected shape comes out like a triangle instead of parabolic. This is because of using “DirichletBC” to apply the displacement the mid-node.
My question is, how can we implement lagrange multipier method/penalty method or partition of stiffness matrix method in FEniCS to implement the said displacement constraint?
I am attaching the MWE below, here I need to import the geometry from gmsh since I need 3d coordinates of the line element for the code to work.
Mesh file: mesh.xdmf - Google Drive

from dolfin import *
from ufl import Jacobian, diag
import matplotlib.pyplot as plt

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
	infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

def tangent(mesh):
	t = Jacobian(mesh)
	return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

t = tangent(mesh)

ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

thick = Constant(0.14) # mm
width = 1000 #mm

Ey = 410 # KN/mm2
nu = 0.25
G = Ey / 2 / (1 + nu)
rho = Constant(0)
g = Constant(0)
Gc = 1.030e-3 #KN/mm
S = thick*width
ES = Ey * S
EI1 = Ey * width * thick ** 3 / 12
EI2 = Ey * width ** 3 * thick / 12
GJ = G*0.26*thick*width**3
kappa = Constant(5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S

delta_u =  7.5e-2
max_am_iterations = 2000
disp_current = 0.0

TFS = TensorFunctionSpace(mesh, "DG", 2)
DG = FunctionSpace(mesh, "DG", 0)
#
Ue = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, Ue*Ue)
v = TestFunction(W)
u = TrialFunction(W)
(w_, theta_) = split(v)
(dw, dtheta) = split(u)

u_new =Function(W, name="displacement")

left = CompiledSubDomain("on_boundary && near(x[0], 0, tol)", tol=1e-12)
right = CompiledSubDomain("on_boundary && near(x[0], 3, tol)", tol=1e-12)
load = CompiledSubDomain("x[1] ==0.0 && near(x[0], 1.5, tol)", tol=1e-6)
applied_disp = Expression("disp", disp=0.0, degree=2)

bc_u_left = DirichletBC(W.sub(0), Constant((0, 0, 0)), left,  )
bc_u_right_y = DirichletBC(W.sub(0).sub(1), Constant(0), right,  )
bc_u_right_z = DirichletBC(W.sub(0).sub(2), Constant(0), right,  )
bc_u_load = DirichletBC(W.sub(0).sub(1), applied_disp, load, )
bc_u = [bc_u_left,  bc_u_right_y,  bc_u_right_z, bc_u_load]
####################################################################

def tgrad(u):
	return dot(grad(u), t)

def generalized_strains(u):
	(w, theta) = split(u)
	return as_vector([dot(tgrad(w), t),
					  dot(tgrad(w), a1)-dot(theta, a2),
					  dot(tgrad(w), a2)+dot(theta, a1),
					  dot(tgrad(theta), t),
					  dot(tgrad(theta), a1),
					  dot(tgrad(theta), a2)])

def generalized_stresses(u):
	return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))

Sig = generalized_stresses(u)
Eps =  generalized_strains(v)

dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 1})
k_form = (Sig[0]*Eps[0] + Sig[3] * Eps[3] + Sig[4] * Eps[4] + Sig[5] * Eps[5] )*dx + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
l_form = Constant(-rho*S*g)*w_[1]*dx


p_disp = LinearVariationalProblem(k_form, l_form,  u_new, bc_u)
solver_disp = LinearVariationalSolver(p_disp)

disp_current += delta_u
applied_disp.disp =  -disp_current
K = assemble(k_form)
solver_disp.solve()
plt.plot(u_new.vector()[1:-1:6])