# How can point loads be applied to the end of the cantilever beam?

**URL:** <https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253>\
**Category:** mesh\
**Created:** [May 5, 2022, 2:52pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253 "2022-05-05T14:52:57Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![pshdasom](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pshdasom/32/3691_2.png) [@pshdasom](https://fenicsproject.discourse.group/u/pshdasom)\
**Post date:** [May 5, 2022, 2:52pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/1 "2022-05-05T14:52:57Z")

</div>

```auto
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.

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [May 5, 2022, 7:53pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/2 "2022-05-05T19:53:16Z")

</div>

Hello

The below illustrates how you can apply a point load on a cantilever beam:

```auto
from dolfin import *

L = 10
W = 1
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
    def inside(self, x, on_boundary):

        return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
    def inside(self, x, on_boundary):

       return on_boundary and abs(x[0] - L) < tol

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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

class PointLoad(UserExpression):

    def __init__ (self, **kwargs):
        super(). __init__ ()

        self.point = kwargs['point']
        self.value = kwargs['value']
    def eval(self, value, x):
        if near (x[0], self.point[0]) and near(x[1], self.point[1]):
            value[0] = self.value[0]
            value[1] = self.value[1]
        else:
            value[0] = 0
            value[1] = 0
    def value_shape(self):
        return (2,)
#End Load
T = PointLoad(point=(10.,1.), value=(0,-1e9), degree=1)
#Solve
a = inner(sigma(u), eps(v))*dx
L =dot(T, v)*ds

# Compute solution
u = Function(V)

solve(a == L, u, bc)
File("disp.pvd") << u

```

---

<div class="post-metadata">

**Author:** ![pshdasom](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pshdasom/32/3691_2.png) [@pshdasom](https://fenicsproject.discourse.group/u/pshdasom)\
**Post date:** [May 6, 2022, 7:57am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/3 "2022-05-06T07:57:58Z")

</div>

```auto
from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
    def inside(self, x, on_boundary):

        return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
    def inside(self, x, on_boundary):

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(v))

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

class PointLoad(UserExpression):

    def __init__ (self, **kwargs):
        super(). __init__ ()

        self.point = kwargs['point']
        self.value = kwargs['value']
    def eval(self, value, x):
        if near (x[0], self.point[0]) and near(x[1], self.point[1]):
            value[0] = self.value[0]
            value[1] = self.value[1]
        else:
            value[0] = 0
            value[1] = 0
    def value_shape(self):
        return (2,)
#End Load
T = PointLoad(point=(2, 0.067), value=(0,-P), degree=1)
#Solve
a = inner(sigma(u), eps(v))*dx
l =dot(T, v)*ds

# Compute solution
u = Function(V, name="Displacement")
solve(a == l, u, bc)

plot(1e3*u, mode="displacement")

print("Maximal deflection:", -u(L, W/2)[1])
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))

```

I checked your answer and calculated it by modifying it according to the physical properties I wanted. However, the results do not match the theoretical values. I would appreciate it if you could check if I made any mistakes.

Solving linear variational problem.  
Object cannot be plotted directly, projecting to piecewise linears.  
Maximal deflection: 6.496007594039327e-05  
Beam theory deflection: 0.0032

L = 2  
W = 0.067  
EI = 5e6  
P(point load) = 6e3

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [May 6, 2022, 10:12pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/4 "2022-05-06T22:12:41Z")

</div>

In the above implementation, `T` is interpreted as a traction which is interpolated into `V` and integrated over the external surface of the beam. The figure below shows the y-component of the traction as it is interpreted by your code. This is very different from applying a point load.

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f06722ca46937f506b4c3dbbbee9547685d1a420.png)

To implement a point load, you should use `PointSource`, e.g. as below

```python
from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
    def inside(self, x, on_boundary):

        return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
    def inside(self, x, on_boundary):

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(v))

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

#End Load
T = PointSource(V.sub(1), Point(2, 0.067), -P)

#Solve
a = inner(sigma(u), eps(v))*dx
l = dot(Constant((0.0,0.0)), v)*ds

K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)

# Compute solution
u = Function(V, name="Displacement")
solve(K, u.vector(), b)

plot(1e3*u, mode="displacement")

print("Maximal deflection:", -u(L, W/2)[1])
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))

```

producing

```nohighlight
Maximal deflection: 0.0031860791474304178
Beam theory deflection: 0.0032

```

---

<div class="post-metadata">

**Author:** ![pshdasom](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pshdasom/32/3691_2.png) [@pshdasom](https://fenicsproject.discourse.group/u/pshdasom)\
**Post date:** [May 9, 2022, 6:45am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/5 "2022-05-09T06:45:42Z")

</div>

I really appreciate. But I have another problem. My goal is to calculate the deflection amount of the cantilever beam with a uniform distribution load and a point load at the same time. However, it is difficult to apply the uniform distribution load in the code you gave me. Is there a way to apply point load and uniform distribution load at the same time?

---

<div class="post-metadata">

**Author:** ![pshdasom](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pshdasom/32/3691_2.png) [@pshdasom](https://fenicsproject.discourse.group/u/pshdasom)\
**Post date:** [May 9, 2022, 7:13am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/6 "2022-05-09T07:13:19Z")

</div>

I solved the interpretation of the uniform distribution load using the following code.

class Top(SubDomain):  
def inside(self, x, on\_boundary):  
return near(x[1], W) and on\_boundary

w = 4e3

top = Top()

top.mark(boundaries, 1)

l = dot(Constant((0, -w)), v)\*ds(1)

However, this code does not apply to the code you gave me with different solve format. I guess the solve format I used is solve (a==l, u, bc), but the solve format in the code you gave me is a matrix. I want to find a way to solve this problem.

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [May 9, 2022, 2:15pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/7 "2022-05-09T14:15:14Z")

</div>

To apply the distributed load, replace the following line in the code I provided:

```python
l = dot(Constant((0.0,0.0)), v)*ds

```

with your variational formulation, i.e.:

```python
l = dot(Constant((0, -w)), v)*ds(1)

```

Also, don’t forget to define a `Measure` for the surface integral after the boundaries have been marked, i.e.

```python
ds = Measure('ds', mesh, subdomain_data=boundaries)

```

A complete code is below:

```python
from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
    def inside(self, x, on_boundary):

        return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
    def inside(self, x, on_boundary):

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(v))

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

# Top load
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], W) and on_boundary
w = 4e3
top = Top()
top.mark(boundaries, 1)
ds = Measure('ds', mesh, subdomain_data=boundaries)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

#End Load
T = PointSource(V.sub(1), Point(2, 0.067), -P)

# Variational formulations
a = inner(sigma(u), eps(v))*dx

# Solution #1: point load only
l = dot(Constant((0,0)), v)*ds
K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)
u_pt = Function(V, name="Displacement, point load")
solve(K, u_pt.vector(), b)
print("Maximal deflection, point load only:", -u_pt(L, W/2)[1])

# Solution #2: distributed load only
l = dot(Constant((0, -w)), v)*ds(1)
K = assemble(a)
b = assemble(l)
bc.apply(K)
bc.apply(b)
u_dist = Function(V, name="Displacement, distributed load")
solve(K, u_dist.vector(), b)
print("Maximal deflection, distributed load:", -u_dist(L, W/2)[1])

# Solution #3: distributed load plus point load
l = dot(Constant((0, -w)), v)*ds(1)
K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)
u = Function(V, name="Displacement, combined load")
solve(K, u.vector(), b)
print("Maximal deflection, combined load:", -u(L, W/2)[1])

plot(1e3*u, mode="displacement")
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))

```

---

<div class="post-metadata">

**Author:** ![JBo](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jbo/32/7006_2.png) [@JBo](https://fenicsproject.discourse.group/u/JBo)\
**Post date:** [May 24, 2024, 8:53am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/8 "2024-05-24T08:53:10Z")

</div>

Good morning  
Your code has been a great help already.  
I am struggling to adapt it for a 1D beam  
Could you please suggest a version ?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 24, 2024, 10:49am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/9 "2024-05-24T10:49:33Z")

</div>

Please make a minimal attempt to add a point load to a 1D beam. Asking someone to write a new code from scratch is a big ask.

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [May 24, 2024, 10:51am UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/10 "2024-05-24T10:51:35Z")

</div>

What do you mean 1D beam ? If you are referring to 1D beam elements based on Euler-Bernoulli or Timoshenko beam theory this is very different from a 2D/3D solid formulation as above. In this case you should look at [Elastic 3D beam structures — Numerical tours of continuum mechanics using FEniCS master documentation](https://comet-fenics.readthedocs.io/en/latest/demo/beams_3D/beams_3D.html)  
A FEniCSx version will be coming soon

---

<div class="post-metadata">

**Author:** ![JBo](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jbo/32/7006_2.png) [@JBo](https://fenicsproject.discourse.group/u/JBo)\
**Post date:** [May 26, 2024, 3:18pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/12 "2024-05-26T15:18:45Z")

</div>

Good afternoon

Actually I have been researching and working on the issue for quite some time now. I apologize if it seemed like I was asking for someone to write a new code from scratch, I thought it was just a few changes from the RectangleMesh example above.

However I’ve created a new thread [here](https://fenicsproject.discourse.group/t/how-to-apply-concentrated-loads-on-beam-elements-in-3d-space/14799) with my progress so far including the code to keep things organized and helpful for the community. I would greatly appreciate it if you could take a look and provide any suggestions or guidance.

Thank you for your time and understanding.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [May 26, 2024, 3:50pm UTC](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/13 "2024-05-26T15:50:34Z")

</div>