```
from fenics import *
from mshr import *
L = 2
H = 0.067
b = 1
EI = 5e6
Ny = 100
Nx = 100
beam = RectangleMesh(Point(0, 0), Point(L, H), Nx, Ny)
beam.init()
def eps(v):
return sym(grad(v))
E = Constant(200e9)
nu = Constant(0.3)
model = "plane_stress"
mu = E/(2*(1+nu))
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
f = 0
w = 4e3
P = 6e3
f_thernal = Constant((0, -f))
T_plane = Constant((0, -w))
T_point = Constant((0, -P))
V = VectorFunctionSpace(beam, "CG", 2)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], H) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
class Point(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]-L) < 0.01 and abs(x[1]-H) < 0.01
boundaries = MeshFunction("size_t", beam, 1)
left = Left()
pt = Point()
top = Top()
boundaries.set_all(0)
top.mark(boundaries, 1)
left.mark(boundaries, 2)
pt.mark(boundaries, 3)
ds = Measure("ds", domain=beam, subdomain_data=boundaries)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
left_end_displacement = Constant((0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, boundaries, 2)]
l = dot(T_point, u_)*ds(3)
u = Function(V, name="Displacement")
solve(a == l, u, bc)
plot(1e3*u, mode="displacement")
print("Maximal deflection:", -u(L, H/2)[1])
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))
```

I’m trying to apply a point load to the end of the cantilever beam. I kept searching for Google and tried to remove the pointwise of bc and redefine the point as a small area, but the result was not good.

I’m looking for a way to fix this. The code above works fine without errors, but the result is strange. I don’t think there’s any grammatical error.