# How to apply local distributed load on 1D cantilever beam in 3D space?

Still trying to learn FEniCS basics, I am trying to apply uniform distributed load on a portion of a 1D cantilever beam in 3D space.
This is how I locate the station :

``````def load_portion(x, on_boundary):
return on_boundary and start_x <= x[0] <= end_x
``````

… mark the facets :

``````facets = MeshFunction("size_t", mesh, 1, 0)         # MeshFunction for facets
``````

… define the differential element of integration :

``````ds_facets = Measure("ds",
domain=mesh,
subdomain_data=facets,
``````

… and finally apply load :

``````distributed_load = Constant(-rho*S*g)

*w_[1]                                             # y direction
*ds_facets(1))                                     # on mark 1 facets
``````

… but it doesn’t match beam theory so far :

``````maximal deflection =  -4.02e-05
beam theory =  -7.04e-05
``````

My main source is Static analysis with external distributed mass
My full code here bollow :

``````# running in Google Colab
try:
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:
!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

!pip install meshio  # mesh handling

# !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                           # mesh handling
# 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:

"""# Sections
"""

# 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))

# Values
thick_val = 1
width_val = thick_val
E_val = 1e5                                        # Young's modulus
rho_val = 1e-3                                     # Density
g_val = 1                                          # Gravity
EI1_val = E_val * width_val * thick_val**3 / 12
S_val = thick_val * width_val

# Define section dimensions
thick = Constant(thick_val)
width = Constant(width_val)

# Define material properties
E = Constant(E_val)
nu = Constant(0.3)      # Poisson's ratio
G = E / 2 / (1 + nu)    # Shear modulus
rho = Constant(rho_val)    # Density
g = Constant(g_val)         # 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"""

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

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

# Locate right_end
def right_end(x, on_boundary):
return near(x[0], length)

end_frac = 0.8                              # as a fraction of length
start_frac = 0.5                            # as a fraction of portion_load_end
end_x = end_frac * length
start_x = start_frac * end_x
return on_boundary and start_x <= x[0] <= end_x

# Create a MeshFunction objects to store the marks
boundaries = MeshFunction("size_t", mesh, 0, 0)     # MeshFunction for vertices
AutoSubDomain(right_end).mark(boundaries, 1)        # mark 1 for right_end

facets = MeshFunction("size_t", mesh, 1, 0)         # MeshFunction for facets

"""# 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 generalized_strains(u):
(w, theta) = split(u)
return as_vector([
dot(tgrad(w), a1) - dot(theta, a2),    # Shear strain in direction 1
dot(tgrad(w), a2) + dot(theta, a1),    # Shear strain in direction 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
ds_vertices = Measure("ds",
domain=mesh,
subdomain_data=boundaries,
ds_facets = Measure("ds",
domain=mesh,
subdomain_data=facets,

# 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

f = 1e3  # concentrated load

l_form = (inner(-f,                                        # f intensity
u_[1])                                     # along y
*ds_vertices(1))                                 # on mark 1 vertices

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

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

# maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]

# Calculate value with beam theory formula
beam_theory_deflection = -(f * length**3) / (3 * EI1_val)

# Print result
print("maximal deflection = ", max_deflection)           # = -40.21200000025746
print("beam theory = ",        beam_theory_deflection)   # = -40.0

# distributed weight

# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)

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

# Extract maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]
# Value for beam theory
beam_theory_deflection = (-rho_val*S_val*g_val * length**4) / (8 * EI1_val)

# Print result
print("maximal deflection = ", max_deflection)         # = -0.00015106000000100
print("beam theory = ",        beam_theory_deflection) # = -0.00015

l_form_concentrated = (inner(-f, u_[1]) *ds_vertices(1))

l_form = l_form_concentrated + l_form_distributed

# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)

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

# Extract maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]

# Value using the beam theory formula
beam_theory_deflection = (
-(f * length**3) / (3 * EI1_val) +
(-rho_val * S_val * g_val * length**4) / (8 * EI1_val))

# Print result
print("maximal deflection = ", max_deflection)           # = -40.21215106025746
print("beam theory = ",        beam_theory_deflection)   # = -40.00015

# Distributed load on a portion of the beam

*w_[1]                                             # y direction
*ds_facets(1))                                     # on mark 1 facets

# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)

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

# Calculating the deflection using the beam theory formula
# https://www.engineersedge.com/beam_bending/beam_bending24.htm

distributed_load_val = -rho_val * S_val * g_val