```
import pyvista
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
# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
[20, 6, 6], 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], L)
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)
# Define free end facets for calculating moment
def free_end(x):
return np.isclose(x[0], L)
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
# Create a MeshTags object to mark the free end
facet_markers = mesh.meshtags(domain, fdim, free_end_facets, np.full(len(free_end_facets), 1, 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(len(u)) + 2 * mu * epsilon(u)
fx = 0.1
fy = -rho * g * 0.5
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0.0, fy)))
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()
# Calculate bending moment Mz
x = ufl.SpatialCoordinate(domain)
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds(1)))
print(f"Bending moment Mz = {Mz:.6f}")
```

This is a 3D derivation of the code found here:

What I am trying to accomplish is at some point on or near the end of the cantilever hopefully towards the bottom center to take the Mz, which in this case is an already known quantity.

In terms of:

```
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds(1)))
```

I get the feeling that ds(1)) may need to change so far now that things moved from 2D to 3D. OK… -x[1] is the co-ordinates of the Y axis of the domain, really what I want is the moment at a specific point. So far there is no errors to report however the Mz seems to be coming out to -0.00 which is an obvious misgnomer.

How can the moment about the Z be taken from a point somewhere near the bottom center of the cantilever?