```
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
# Define the beam properties
L = 120.0 # Length in inches
H = 8.0 # Height in inches
W = 4.0 # Width in inches
I = (W * H**3) / 12 # Moment of inertia in inches^4
# Mesh properties
Nx = 160
Ny = 5
# Create the mesh
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[(0.0, -H / 2), (L, H / 2)],
(Nx, Ny),
diagonal=mesh.DiagonalType.crossed,
)
# Material properties
E = fem.Constant(domain, 2000000.0)
nu = fem.Constant(domain, 0.3)
mu = (E / (2*(1+nu)))
lmbda = ((E * nu) / ((1+nu) * (1-2*nu)))
# Define strain and stress
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return lmbda * ufl.tr(eps(v)) * ufl.Identity(2) + 2.0 * mu * eps(v)
# Applied load
fx = 0.1
fy = 50.0
f = fem.Constant(domain, (fx, -fy / (L * W)))
# Theoretical maximum moment
theoretical_moment_max = (((-fy / L) * (L**2)) / (2))
# Define function space
V = fem.functionspace(domain, ("P", 2, (2,)))
# Define variational problem
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), eps(u_)) * ufl.dx
l = ufl.inner(f, u_) * ufl.dx
# Boundary conditions
def left(x):
return np.isclose(x[0], 0.0)
fdim = domain.topology.dim - 1
marked_values = []
marked_facets = []
# Concatenate and sort the arrays based on facet indices
facets = mesh.locate_entities_boundary(domain, fdim, left)
marked_facets.append(facets)
marked_values.append(np.full_like(facets, 1))
marked_facets = np.hstack(marked_facets)
marked_values = np.hstack(marked_values)
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets]
)
# Measure for boundary integration
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
# Apply Dirichlet BCs (fixing the left end)
left_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, left_dofs)]
# Solve the problem
u = fem.Function(V, name="Displacement")
problem = LinearProblem(a, l, bcs, u=u)
problem.solve()
# Calculate reaction forces and moments
x = ufl.SpatialCoordinate(domain)
Rx = fem.assemble_scalar(fem.form(-sigma(u)[0, 0] * ds(1)))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[1, 1] * ds(1)))
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(u)[0, 0] * ds(1)))
print(f"Horizontal reaction Rx = {Rx:.6f}")
print(f" (analytic = {-L * H * fx})")
print("-" * 50)
print(f"Vertical reaction Ry = {Ry:.6f}")
print(f" (analytic = {-L * H * fy})")
print("-" * 50)
print(f"Bending moment Mz = {Mz:.6f}")
print(f" (analytic = {H * L**2 / 2 * fy})")
print("-" * 50)
print("\n")
print(f"Theoretical M max: {theoretical_moment_max:.6f} inch\lbs")
```

```
--------------------------------------------------
Bending moment Mz = -449.826670
(analytic = 2880000.0)
--------------------------------------------------
Theoretical M max: -3000.000000 inch\lbs
```

So far the Mz I am getting is nowhere near -3000.00 lbs/inch however I really am not quire sure yet…

So far I have these formulas to go by where w is pounds each inch of material and l is the length in inches of the cantilever beam.

Why this discrepancy between Mz and theoretical moment?