How to apply concentrated loads on beam elements in 3D space?

Hi all
I want to use FEnics to calculate sailing ship rigging.
Starting with basic problem I struggle to have a simple cantilever beam example working.
Especially having the concentrated loads case to match the beam theory.

I have been through these for inspiration:

Although I learnt from them a general understanding of the code, I am still having trouble getting my example to work for concentrated loads. I would greatly appreciate it if someone could help me figure out what I am doing wrong.

Here is the code I have so far:

# Attempt to import google.colab to check if running in Google Colab
try:
    import google.colab
except ImportError:
    # If not in Colab, import necessary FEniCS packages
    import ufl_legacy
    import dolfin
else:
    # If in Colab, ensure FEniCS is installed
    try:
        import ufl_legacy
        import dolfin
    except ImportError:
        # Download and install FEniCS if not already installed
        !wget "https://fem-on-colab.github.io/releases/fenics-install-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
        import ufl_legacy
        import dolfin

# Install meshio for mesh handling
!pip install meshio

# Importing necessary libraries from FEniCS and others
from dolfin import *
from ufl_legacy import Jacobian, diag
import numpy as np
import meshio

"""# Mesh"""

# Defining the length of the cantilever beam and the number of elements
length = 10.0
N = 80 # number of elements

# Creating an array of points for the mesh (N+1 points in 3D space)
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)

# Defining the cells (elements) of the mesh
cells = np.array([[i, i+1] for i in range(N)])

# Writing the mesh to an XDMF file using meshio
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

# Reading the mesh from the XDMF file into a FEniCS Mesh object
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

"""# Sections
Same as [CTV24_FEniCS-elastic-3D-beam-structure-tutorial_2024-05-16.ipynb](https://colab.research.google.com/drive/1PFxZcjQ5UaI9WL-cUiI8h9JzDJt5qgqf#scrollTo=hXViYLk3nLwo)
"""

# Function to calculate the tangent vector of the mesh elements
def tangent(mesh):
    t = Jacobian(mesh)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

# Calculate the tangent vector of the mesh
t = tangent(mesh)

ez = as_vector([0, 0, 1])   # unit vector in the z-direction.
a1 = cross(t, ez)           # first cross-sectional vector
a1 /= sqrt(dot(a1, a1))     # normalize
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

# Define section dimensions
thick = Constant(1)
width = thick#/3

# Define material properties
E = Constant(1e5)       # Young's modulus
nu = Constant(0.3)      # Poisson's ratio
G = E / 2 / (1 + nu)    # Shear modulus
rho = Constant(1e-3)    # Density
g = Constant(1)         # Gravity

# Define section properties
S = thick * width                 # Cross-sectional area
ES = E * S                        # Axial rigidity
EI1 = E * width * thick**3 / 12   # Bending rigidity about axis 1
EI2 = E * width**3 * thick / 12   # Bending rigidity about axis 2
GJ = G * 0.26 * thick * width**3  # Torsional rigidity
kappa = Constant(5. / 6.)         # Shear correction factor
GS1 = kappa * G * S               # Shear rigidity in direction 1
GS2 = kappa * G * S               # Shear rigidity in direction 2

"""# Functions"""

# Defining the function space for the displacement and rotation fields
Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

# Defining the test and trial functions
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

# Define generalized strains as a function of displacement and rotation
def tgrad(u):
    return dot(grad(u), t)

def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([
        dot(tgrad(w), t),                      # Axial strain
        dot(tgrad(w), a1) - dot(theta, a2),    # Shear strain in direction 1
        dot(tgrad(w), a2) + dot(theta, a1),    # Shear strain in direction 2
        dot(tgrad(theta), t),                  # Torsional strain
        dot(tgrad(theta), a1),                 # Bending strain about axis 1
        dot(tgrad(theta), a2)                  # Bending strain about axis 2
    ])

# Define generalized stresses as a function of generalized strains
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
               generalized_strains(u))

# Defining the variational form
Sig = generalized_stresses(du)  # stresses
Eps =  generalized_strains(u_)  # strain

# will be used to integrate the shear terms
dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 4})

# Define the bilinear form (stiffness matrix) of the beam
k_form = (
    sum([Sig[i] * Eps[i] * dx for i in [0, 3, 4, 5]]) # normal and bending terms
    + (Sig[1] * Eps[1] + Sig[2] * Eps[2])             # shear terms
    * dx_shear                                        # Integrate
)

"""# Boundaries"""

# Define boundary condition for clamping (fixed end at x=0))
def clamp(x, on_boundary):
    return near(x[0], 0.)

# Apply Dirichlet boundary condition
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)

"""# Concentrated load only"""

# Define concentrated load at the end of the beam
P = 1e3  # Point load in N
T = PointSource(W.sub(1), Point(length, 0, 0), -P)

# No distributed loads
l = Constant(0) * w_[1] * dx

# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
T.apply(b)      # Apply the point load to the load vector
bc.apply(K)     # Apply boundary conditions to the stiffness matrix
bc.apply(b)     # Apply boundary conditions to the load vector

# Solve for displacements and rotations
u_pt = Function(W, name="Displacement, point load")
solve(K, u_pt.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, point load only:", -u_pt(length, 0, 0)[1])

# Calculating the deflection using the beam theory formula
beam_theory_deflection = (P * length**3) / (3 * EI1)

# Print the result
print("Beam theory deflection:", float(beam_theory_deflection))

# Maximal deflection, point load only: 6.000000000035433
# Beam theory deflection: 40.0

"""# Distributed load only"""

# distributed weight
l = Constant(-rho*S*g)*w_[1]*dx

# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
bc.apply(K)     # Apply boundary conditions to the stiffness matrix
bc.apply(b)     # Apply boundary conditions to the load vector

# Solve for displacements and rotations
u_dist = Function(W, name="Displacement, distributed load")
solve(K, u_dist.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, distributed load only:", u_dist(length, 0, 0)[1])

# Calculating the deflection using the beam theory formula for distributed load
beam_theory_deflection_distributed = (-rho*S*g * length**4) / (8 * EI1)

# Print the result
print("Beam theory deflection for distributed load:", float(beam_theory_deflection_distributed))

"""# Combined loads"""

# Define concentrated load at the end of the beam
P = 1e3  # Point load in N
T = PointSource(W.sub(1), Point(length, 0, 0), -P)

# Define distributed weight
l = Constant(-rho*S*g)*w_[1]*dx

# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
T.apply(b)      # Apply the point load to the load vector
bc.apply(K)     # Apply boundary conditions to the stiffness matrix
bc.apply(b)     # Apply boundary conditions to the load vector

# Solve for displacements and rotations
u = Function(W, name="Displacement, combined load")
solve(K, u.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, combined loads:", -u(length, 0, 0)[1])

I have the same result using AutoSubDomain method instead of PointSource to apply loads.
Must be something else…

"""# Boundaries"""

# Define boundary condition for clamping (fixed end at x=0))
def clamp(x, on_boundary):
    return near(x[0], 0.)

bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)

# Define boundary condition for load ( at x=length))
def right_end(x, on_boundary):
    return near(x[0], length)

boundaries = MeshFunction("size_t", mesh, 0, 0)
AutoSubDomain(right_end).mark(boundaries, 1)        # mark 1 for load

f = 1e3  # concentrated load in N

"""# Functions"""

# Defining the function space for the displacement and rotation fields
Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

# Defining the test and trial functions
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

# Define generalized strains as a function of displacement and rotation
def tgrad(u):
    return dot(grad(u), t)

def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([
        dot(tgrad(w), t),                      # Axial strain
        dot(tgrad(w), a1) - dot(theta, a2),    # Shear strain in direction 1
        dot(tgrad(w), a2) + dot(theta, a1),    # Shear strain in direction 2
        dot(tgrad(theta), t),                  # Torsional strain
        dot(tgrad(theta), a1),                 # Bending strain about axis 1
        dot(tgrad(theta), a2)])                # Bending strain about axis 2

# Define generalized stresses as a function of generalized strains
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
               generalized_strains(u))

# Defining the variational form
Sig = generalized_stresses(du)  # stresses
Eps =  generalized_strains(u_)  # strain

# will be used to integrate the shear terms
metadata = {"quadrature_scheme": "default", "quadrature_degree": 2}
dS = Measure("dS", domain=mesh, subdomain_data=boundaries, metadata=metadata)
dx_shear = dx(scheme="default", metadata=metadata)

# Define the bilinear form (stiffness matrix) of the beam
k_form = (
    sum([Sig[i] * Eps[i] * dx for i in [0, 3, 4, 5]]) # normal and bending terms
    + (Sig[1] * Eps[1] + Sig[2] * Eps[2])             # shear terms
    * dx_shear                                        # Integrate
)

"""# Concentrated load only"""

l_form = inner(f, 2*avg(u_[2]))*dS(1)  # apply f on mark 1

# Assemble the equation
A, b = assemble_system(k_form, l_form, bc)

# Solve the equation
u = Function(W)
solve(A, u.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, point load only:", -u_pt(length, 0, 0)[1])

# Calculating the deflection using the beam theory formula
beam_theory_deflection = (f * length**3) / (3 * EI1)

# Print the result
print("Beam theory deflection:", float(beam_theory_deflection))

Finally one error was there at least …

T = PointSource(W.sub(0).sub(1), Point(length, 0, 0), -P)

Point load T now applies to the correct subspace (W.sub(0).sub(1))
with :

  • W.sub(0) : the displacement vector space Ue
  • W.sub(0).sub(1) : the y component

Now calculations match:

Maximal deflection, point load only: 40.21200000025746
Beam theory deflection: 40.0

Still interested in a neater version of this code or a fix for the AutoSubDomain method.

Did you try getting rid of the 2 in

l_form = inner(f, 2*avg(u_[2]))*dS(1)  # apply f on mark 1

and use

l_form = inner(f, avg(u_[2]))*dS(1)  # apply f on mark 1

Thank you for your previous response.
I’m still having trouble getting the AutoSubDomain method to work. When I try to apply a point load, I’m getting a displacement of 0.
I’m confident there’s a mistake in my code, and I’m still eager to resolve this. Any guidance would be greatly appreciated!

# Attempt to import google.colab to check if running in Google Colab
try:
    import google.colab
except ImportError:
    # If not in Colab, import necessary FEniCS packages
    import ufl_legacy
    import dolfin
else:
    # If in Colab, ensure FEniCS is installed
    try:
        import ufl_legacy
        import dolfin
    except ImportError:
        # Download and install FEniCS if not already installed
        !wget "https://fem-on-colab.github.io/releases/fenics-install-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
        import ufl_legacy
        import dolfin

# Install meshio for mesh handling
!pip install meshio

# !apt-get install -qq xvfb libgl1-mesa-glx
# !pip install pyvista -qq

# Importing necessary libraries from FEniCS and others
from dolfin import *
from ufl_legacy import Jacobian, diag
import numpy as np
import meshio
# import pyvista

"""# Mesh"""

# Defining the length of the cantilever beam and the number of elements
length = 10.0
N = 10 # number of elements

# Creating an array of points for the mesh (N+1 points in 3D space)
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)

# Defining the cells (elements) of the mesh
cells = np.array([[i, i+1] for i in range(N)])

# Writing the mesh to an XDMF file using meshio
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

# Reading the mesh from the XDMF file into a FEniCS Mesh object
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

"""# Sections
Same as [CTV24_FEniCS-elastic-3D-beam-structure-tutorial_2024-05-16.ipynb](https://colab.research.google.com/drive/1PFxZcjQ5UaI9WL-cUiI8h9JzDJt5qgqf#scrollTo=hXViYLk3nLwo)
"""

# Function to calculate the tangent vector of the mesh elements
def tangent(mesh):
    t = Jacobian(mesh)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

# Calculate the tangent vector of the mesh
t = tangent(mesh)

ez = as_vector([0, 0, 1])   # unit vector in the z-direction.
a1 = cross(t, ez)           # first cross-sectional vector
a1 /= sqrt(dot(a1, a1))     # normalize
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

# Define section dimensions
thick = Constant(1)
width = thick#/3

# Define material properties
E = Constant(1e5)       # Young's modulus
nu = Constant(0.3)      # Poisson's ratio
G = E / 2 / (1 + nu)    # Shear modulus
rho = Constant(1e-3)    # Density
g = Constant(1)         # Gravity

# Define section properties
S = thick * width                 # Cross-sectional area
ES = E * S                        # Axial rigidity
EI1 = E * width * thick**3 / 12   # Bending rigidity about axis 1
EI2 = E * width**3 * thick / 12   # Bending rigidity about axis 2
GJ = G * 0.26 * thick * width**3  # Torsional rigidity
kappa = Constant(5. / 6.)         # Shear correction factor
GS1 = kappa * G * S               # Shear rigidity in direction 1
GS2 = kappa * G * S               # Shear rigidity in direction 2

"""# Function Space"""

# Defining the function space for the displacement and rotation fields
Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3) # vector element of degree 2 for displacement, with three components
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3) # vector element of degree 1 for rotation, with three components
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

"""# Boundaries"""

# Define boundary condition for clamping (fixed end at x=0))
def clamp(x, on_boundary):
    return near(x[0], 0.)

# Apply Dirichlet boundary condition
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)

# Define boundary condition for load ( at x=length))
def right_end(x, on_boundary):
    return near(x[0], length)

boundaries = MeshFunction("size_t", mesh, 0, 0)
AutoSubDomain(right_end).mark(boundaries, 1)        # mark 1 for load

"""# Functions"""

# Defining the test and trial functions
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

# Define generalized strains as a function of displacement and rotation
def tgrad(u):
    return dot(grad(u), t)

def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([
        dot(tgrad(w), t),                      # Axial strain
        dot(tgrad(w), a1) - dot(theta, a2),    # Shear strain in direction 1
        dot(tgrad(w), a2) + dot(theta, a1),    # Shear strain in direction 2
        dot(tgrad(theta), t),                  # Torsional strain
        dot(tgrad(theta), a1),                 # Bending strain about axis 1
        dot(tgrad(theta), a2)])                # Bending strain about axis 2

# Define generalized stresses as a function of generalized strains
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
               generalized_strains(u))

# Defining the variational form
Sig = generalized_stresses(du)  # stresses
Eps =  generalized_strains(u_)  # strain

# will be used to integrate the shear terms
metadata = {"quadrature_scheme": "default", "quadrature_degree": 4}
dx_shear = dx(scheme="default", metadata=metadata)

# Define the bilinear form (stiffness matrix) of the beam
k_form = (
    sum([Sig[i] * Eps[i] * dx for i in [0, 3, 4, 5]]) # normal and bending terms
    + (Sig[1] * Eps[1] + Sig[2] * Eps[2])             # shear terms
    * dx_shear                                        # Integrate
)

"""# Concentrated load only"""

# concentrated load on mark 1
P = 1e3  # Point load in N
dS = Measure("dS", domain=mesh, subdomain_data=boundaries, metadata=metadata)
l_form = inner(P, avg(u_[1]))*dS(1)  # apply f along y on mark 1

# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l_form)
bc.apply(K)     # Apply boundary conditions to the stiffness matrix
bc.apply(b)     # Apply boundary conditions to the load vector

# Solve for displacements and rotations
u_pt = Function(W, name="Displacement, point load")
solve(K, u_pt.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, point load only:", -u_pt(length, 0, 0)[1])   # = -0.0

# Calculating the deflection using the beam theory formula
beam_theory_deflection = (P * length**3) / (3 * EI1)

# Print the result
print("Beam theory deflection:", float(beam_theory_deflection)) # = 40.0

Try replacing dS with ds
ds = Measure("ds", domain=mesh, subdomain_data=boundaries, metadata=metadata)
and change the following
l_form = inner(P, u_[1])*ds(1)

1 Like

Nice ! You saved my night :pray:

Maximal deflection, point load only: -40.21200000025746
Beam theory deflection: -40.0

Full code for point load with subdomain method goes like :

# Attempt to import google.colab to check if running in Google Colab
try:
    import google.colab
except ImportError:
    # If not in Colab, import necessary FEniCS packages
    import ufl_legacy
    import dolfin
else:
    # If in Colab, ensure FEniCS is installed
    try:
        import ufl_legacy
        import dolfin
    except ImportError:
        # Download and install FEniCS if not already installed
        !wget "https://fem-on-colab.github.io/releases/fenics-install-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
        import ufl_legacy
        import dolfin

# Install meshio for mesh handling
!pip install meshio

# Importing necessary libraries from FEniCS and others
from dolfin import *
from ufl_legacy import Jacobian, diag
import numpy as np
import meshio

"""# Mesh"""

# Defining the length of the cantilever beam and the number of elements
length = 10.0
N = 10 # number of elements

# Creating an array of points for the mesh (N+1 points in 3D space)
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)

# Defining the cells (elements) of the mesh
cells = np.array([[i, i+1] for i in range(N)])

# Writing the mesh to an XDMF file using meshio
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

# Reading the mesh from the XDMF file into a FEniCS Mesh object
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

"""# Sections
Same as [CTV24_FEniCS-elastic-3D-beam-structure-tutorial_2024-05-16.ipynb](https://colab.research.google.com/drive/1PFxZcjQ5UaI9WL-cUiI8h9JzDJt5qgqf#scrollTo=hXViYLk3nLwo)
"""

# Function to calculate the tangent vector of the mesh elements
def tangent(mesh):
    t = Jacobian(mesh)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))

# Calculate the tangent vector of the mesh
t = tangent(mesh)

ez = as_vector([0, 0, 1])   # unit vector in the z-direction.
a1 = cross(t, ez)           # first cross-sectional vector
a1 /= sqrt(dot(a1, a1))     # normalize
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

# Define section dimensions
thick = Constant(1)
width = thick#/3

# Define material properties
E = Constant(1e5)       # Young's modulus
nu = Constant(0.3)      # Poisson's ratio
G = E / 2 / (1 + nu)    # Shear modulus
rho = Constant(1e-3)    # Density
g = Constant(1)         # Gravity

# Define section properties
S = thick * width                 # Cross-sectional area
ES = E * S                        # Axial rigidity
EI1 = E * width * thick**3 / 12   # Bending rigidity about axis 1
EI2 = E * width**3 * thick / 12   # Bending rigidity about axis 2
GJ = G * 0.26 * thick * width**3  # Torsional rigidity
kappa = Constant(5. / 6.)         # Shear correction factor
GS1 = kappa * G * S               # Shear rigidity in direction 1
GS2 = kappa * G * S               # Shear rigidity in direction 2

"""# Function Space"""

# Defining the function space for the displacement and rotation fields
Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

"""# Boundaries"""

def clamp(x, on_boundary):
    return near(x[0], 0.)

bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)

# Define boundary condition for load (fixed end at x=length))
def right_end(x, on_boundary):
    return near(x[0], length)

boundaries = MeshFunction("size_t", mesh, 0, 0)
AutoSubDomain(right_end).mark(boundaries, 1)        # mark 1 for load

f = 1e3  # concentrated load in N

"""# Functions"""

# Defining the test and trial functions
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

# Define generalized strains as a function of displacement and rotation
def tgrad(u):
    return dot(grad(u), t)

def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([
        dot(tgrad(w), t),                      # Axial strain
        dot(tgrad(w), a1) - dot(theta, a2),    # Shear strain in direction 1
        dot(tgrad(w), a2) + dot(theta, a1),    # Shear strain in direction 2
        dot(tgrad(theta), t),                  # Torsional strain
        dot(tgrad(theta), a1),                 # Bending strain about axis 1
        dot(tgrad(theta), a2)])                # Bending strain about axis 2

# Define generalized stresses as a function of generalized strains
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
               generalized_strains(u))

# Defining the variational form
Sig = generalized_stresses(du)  # stresses
Eps =  generalized_strains(u_)  # strain

# will be used to integrate the shear terms
metadata = {"quadrature_scheme": "default", "quadrature_degree": 4}
ds = Measure("ds", domain=mesh, subdomain_data=boundaries, metadata=metadata)
dx_shear = dx(scheme="default", metadata=metadata)

# Define the bilinear form (stiffness matrix) of the beam
k_form = (
    sum([Sig[i] * Eps[i] * dx for i in [0, 3, 4, 5]]) # normal and bending terms
    + (Sig[1] * Eps[1] + Sig[2] * Eps[2])             # shear terms
    * dx_shear                                        # Integrate
)

"""# Concentrated load only"""

l_form = (inner(f,                   # apply f
                u_[1])               # along y
          *ds(1))                    # on mark 1

# Assemble the equation
K, b = assemble_system(k_form, l_form, bc)

# Solve the equation
u = Function(W)
solve(K, u.vector(), b)

# Print maximal deflection at the free end of the beam
print("Maximal deflection, point load only:", -u(length, 0, 0)[1])

# Calculating the deflection using the beam theory formula
beam_theory_deflection = (f * length**3) / (3 * EI1)

# Print the result
print("Beam theory deflection:", float(beam_theory_deflection))