You can simply integrate the tractions over surface:
t = ufl.inner(P * n, n1)*ds(1)
R = dolfinx.fem.assemble_scalar(t)
where P is the stress tensor, n is the unit outward normal vector and n1 is the direction that you want to resolve the traction into (above snippet as not been tested). But, this computation is not conservative, i.e. the reactions and applied forces will not sum to exactly zero. This is a characteristic of Galerkin methods with continuous bases, and may or may not be an issue for your application.