```
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh
import dolfinx.fem.petsc
import basix
L = 120.0 # Extending the domain to the point of interest
H = 8.0
Nx = 40 # Increase the number of elements to maintain mesh quality
Ny = 40
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[(0.0, -H / 2), (L, H / 2)],
(Nx, Ny),
diagonal=mesh.DiagonalType.crossed,
)
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)
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)
fx = 0.0
fy = -0.1
f = fem.Constant(domain, (fx, fy))
msh_cell_type = domain.basix_cell()
Ve = basix.ufl.element(
'Lagrange', msh_cell_type, degree=2, shape=(domain.geometry.dim,))
V = fem.functionspace(domain, Ve)
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), eps(u_)) * ufl.dx
l = ufl.inner(f, u_) * ufl.dx
def left(x):
return np.isclose(x[0], 0.0)
fdim = domain.topology.dim - 1
marked_values = []
marked_facets = []
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]
)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
left_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, left_dofs)]
u = fem.Function(V, name="Displacement")
problem = dolfinx.fem.petsc.LinearProblem(a, l, bcs, u=u)
problem.solve()
# Moment of inertia for the cross-section
I = (4 * 8**3) / 12
# Initialize bending moment
Mz_total = 0.0
# Loop through sections of the beam to compute the cumulative bending moment
for x_pos in np.linspace(0, L, Nx+1):
def at_section(x):
return np.isclose(x[0], x_pos)
section_facets = mesh.locate_entities_boundary(domain, fdim, at_section)
if section_facets.size > 0:
marked_section = np.hstack(section_facets)
section_tag = mesh.meshtags(
domain, fdim, marked_section, np.full_like(marked_section, 2)
)
section_ds = ufl.Measure("ds", domain=domain, subdomain_data=section_tag, subdomain_id=2)
Mz_section = fem.assemble_scalar(fem.form(-ufl.SpatialCoordinate(domain)[1] * sigma(u)[0, 0] * section_ds))
Mz_total += Mz_section
# Multiply by the moment of inertia
Mz_total *= I
print(f"Bending moment Mz at x={L} = {Mz_total:.6f}")
# Calculate the theoretical maximum bending moment for comparison
theoretical_max_moment = fy * L * L / 2
print('Theoretical max moment:', theoretical_max_moment)
```

So far I have dolfinx showing more moment for the distributed load.

I did try to make use of:

```
#x = ufl.SpatialCoordinate(domain)
#n = ufl.FacetNormal(domain)
#My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(1)))
```

however I got a long traceback. Since that didn’t work I took a guess at some method that I thought may have a chance at matching the 720.0 figure of the American Wood Council. There still is a pretty good sized discrepancy.

Have I done anything at all correct here? Could it be shown how to come close to 720.0 as the AWC says, and maybe some details why there is a difference right now of about 200 inch/lbs of moment, and maybe how to remedy that and do this correctly?