Derivative of stress resultant - real space again?

This is related to this question.

In dolfinx, for a somewhat unusual linear elasticity problem, I need to build a matrix defined by the derivative of the stress vector resultant, computed over a surface, wrt the displacement trial functions.

In the old dolfin this was accomplished by defining a linear form, with a test function defined over a “Real” space that allowed to compute (assemble) the stress resultant vector, and then taking the derivative of this linear form. Something like

R3_ELEMENT = VectorElement("R", mesh.ufl_cell(), 0, 3)
R3 = FunctionSpace(mesh, R3_ELEMENT)
RV3F = TestFunction(R3)
UF3_ELEMENT = VectorElement("CG", mesh.ufl_cell(), 1, 3)
UF3 = FunctionSpace(mesh, UF3_ELEMENT)
U = Function(UF3)
UT = TrialFunction(UF3)


S = dot(stress_n, RV3F) * ds
L_res_f = derivative(S, U, UT)
A = assemble(L_res_f)

where stress_n is the stress vector, computed as a function of U, and defined by the dot product between the stress tensor and the outward pointing normal.

My understanding, also from this, is that I can compute the derivative of one component of the stress vector resultant. This would give me a vector. However, what I really need is the Jacobian matrix of the whole resultant.

I’m wondering whether there is a way to avoid computing the derivatives of each component and then using them to assemble, by hand, the sought Jacobian matrix.