How to project an element-wise constant tensor on the boundary of a 3D domain

Hi,

I have a 3D tetrahedral domain over which I solve the Navier-Stokes equations . After solving I compute a 3D stress tensor using:

sigma = 2 * eta * sym(grad(u_)) - p_ * Identity(len(u))
tensor = TensorFunctionSpace(mesh, 'CG', 1)
stress = project(sigma, tensor)

Now I want to know the average stress tensor per 2D boundary element (i.e. each boundary triangle). For this purpose I extract the mesh boundary and define a zeroth order, discontinuous Lagrange TensorFunctionSpace on this mesh. Finally I project the stress tensor onto this TensorFunctionSpace:

bmesh = BoundaryMesh(mesh, "exterior")
btensor = TensorFunctionSpace(bmesh, 'DG', 0)
stress = interpolate(stress, btensor)

Here is where I do not get the results I anticipated: I get a surface mesh with the stress but I do not get a constant value per surface element: according to paraview, the stress is defined per element node and moreover, the stress is not constant within the element. Where did I go wrong?

Thank you in advance.

P.S.

This approach did seem to work for a kind of similar 3D situation:

tensor = TensorFunctionSpace(mesh, 'DG', 0)
stress = project(sigma, tensor)

results into a constant value for the stress in each tetrahedron (i.e. stress is represented as "cell data" in paraview, i.e. a single value per element)