Obtaining surface force as function of global position

Errors
1

I am working on a problem with multiple meshes (multi-body contact) and as part of the subproblem formulation, I am:
(1) solving the displacements of the first body assuming the other body is rigid,
(2) calculating the resulting surface forces, and
(3) applying these forces to the second body and solving its displacements.

This will then iterate until the two displacement fields converge. What I am wondering is: Is there a way to calculate resulting surface tractions as a function of the global position? Specifically what I am thinking is to:
(1) calculate the surface tractions at the nodes of the first body
(2) interpolate these values given the nodes global spatial coordinate, and
(3) resample this interpolated function at the the nodes of the second body.

I have seen a number of topics that give the full reaction force on (for example) a fixed boundary, but have not been able to find quite what I need. I have a MWE that brings me up to an equivalent sticking point (ie. I have the tractions as a ComponentTensor but I’m not sure how best to evaluate this as a vector with organization akin to the displacement vector).

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, cpp, nls, la


# simple 2D Constitutive relationship (compressible NH)
C = 1e6
k = 1e10
W = lambda F: (C/2)*(ufl.tr(F.T * F) -2 - 2*ufl.ln(ufl.det(F)))+ (k/2)*((ufl.det(F)-1)**2)

# function to calculate the true stress
def sigma(r):
    I = ufl.variable(ufl.Identity(2))
    F = ufl.variable(I + ufl.grad(r))
    J = ufl.variable(ufl.det(F))
    psi = W(F)
    return (1/J)*ufl.diff(psi,F)*F.T

# defining a mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD,[[0.0,0.0],[1,1]],[10,10],mesh.CellType.triangle)
V_vec = fem.VectorFunctionSpace(domain,("CG",2)) 

# displacement function (here initialized to zero, but would eventually store a solution)
u = fem.Function(V_vec)

# normal vectors
n = ufl.FacetNormal(domain)

# stress state for the given displacement
S = sigma(u)

# tractions in the direction of the normal vectors
T = S*n

### want T as a vector in the same form as u (ie. [T_1x T_1y T_2x T_2y ...])

Any help is greatly appreciated!