There is no way the bending moment can be computed as the curl of a displacement field…
You are using a 3D elastic code, the resulting moment over a surface S around the origin is given by:
\boldsymbol{M} = \int_S \boldsymbol{x} \times \boldsymbol{\sigma n} dS
which contains 3 values : a torsion component and two bending components. To have the bending moment in direction 1, you have to compute something like
n = ufl.FacetNormal(domain)
My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(surf_idx)))
if the surface S is indexed by surf_idx.
but this will not be strictly in equilibrium with the loading. To have consistent reaction force, I suggest you follow this procedure:
which can be easily extended to 3D.