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

**URL:** <https://fenicsproject.discourse.group/t/how-to-prescribe-displacement-at-the-mid-points-of-a-beam-element-this-cant-be-defined-variationally-or-by-using-dirichletbc/5082>\
**Category:** Linear Algebra\
**Created:** [February 12, 2021, 11:12am UTC](https://fenicsproject.discourse.group/t/how-to-prescribe-displacement-at-the-mid-points-of-a-beam-element-this-cant-be-defined-variationally-or-by-using-dirichletbc/5082 "2021-02-12T11:12:48Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ananya\_Bijaya](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ananya_bijaya/32/333_2.png) [@Ananya\_Bijaya](https://fenicsproject.discourse.group/u/Ananya_Bijaya)\
**Post date:** [February 12, 2021, 11:12am UTC](https://fenicsproject.discourse.group/t/how-to-prescribe-displacement-at-the-mid-points-of-a-beam-element-this-cant-be-defined-variationally-or-by-using-dirichletbc/5082/1 "2021-02-12T11:12:48Z")

</div>

I am using the beam element formulation given in [‘Numerical tours of FEniCS’](https://comet-fenics.readthedocs.io/en/latest/demo/beams_3D/beams_3D.html) 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](https://drive.google.com/file/d/1UsqGPknX3Nc46bE0AumQ2Lx9efAfzFSh/view?usp=sharing)

```python
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])

```