```
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
# Geometry and material properties
L = 120.0
W = 4.0
H = 8.0
rho = 0.48
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
E = 2000000 # Modulus of elasticity
nu = 0.3 # Poisson's ratio
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
g = gamma
Nz = 121
Nx = 13
Ny = 13
# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
[121, 13, 13], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.topology.dim, )))
# Define boundary conditions
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = fem.Constant(domain, np.array([0, 0, 0], dtype=default_scalar_type))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
def free_end(x):
return np.isclose(x[0], L)
# Find facets on the boundary where x = L (free end)
fdim = domain.topology.dim - 1
# Ensure that the facets at the free end (x = L) are correctly marked
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
facet_markers = mesh.meshtags(domain, fdim, free_end_facets, np.full(len(free_end_facets), 2, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
# Define the elasticity problem
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(3) + 2 * mu * epsilon(u)
fx = 0.1
self_w = -rho * g * 0.25
w = -0.25
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0.0, -0.00011, 0.0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
'''
https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html#second-method-using-the-work-of-internal-forces
LaTeX:::
\[
M_z = \int_{x=0} (\vec{OM} \times \vec{T}) \times \vec{e}_z dS = -\int_{x=0} (y\sigma_x_x) dS = f_{j_x} - \frac{L^2}{2} = H
\]
'''
# Calculate bending moment Mz
x = ufl.SpatialCoordinate(domain)
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds()))
print(f"Bending moment Mz = {Mz:.6f}")
theoretical_moment = (w*120*120) / 2
print('theoretical_moment:::', -theoretical_moment)
'''
Bending moment Mz = -1779.556528
theoretical_moment::: -1800.0
'''
```

So far my intent is to set up a distributed load case where, small case cursive w is equal to 0.25 lbs per inch. The reason why I have modified the Y-axis force of the the force vector f to -0.00011, is that is the only way to get the Mz to come out anywhere close to the theoretical moment…

It seems so far that:

```
def free_end(x):
return np.isclose(x[0], L)
```

and:

```
def free_end(x):
return np.isclose(x[0], 0)
```

are producing the same results for Mz, however I would expect that they should be different if what I am intending to do is working then Mz should decrease the closer that the free_end facets are to the u_D boundary facets…

f = fem.Constant(domain, default_scalar_type((0.0, -0.00011, 0.0))).

Why is the force having to be scaled down so much to achieve a matching Mz and theoretical moment?

Could I get some help to set things up and get them working correctly for a distributed load case described and computed by, print(‘theoretical_moment:::’, -theoretical_moment)?