How to obtain the moment distribution of a beam fixed at both ends

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:

https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html#second-method-using-the-work-of-internal-forces

which can be easily extended to 3D.

1 Like