Based on the answer I got on this question, I believe I am able to correctly compute the wall shear stress. I will put the rationale/ solution below for the sake of completeness.
The problem with the method as stated in the question is that the wall shear stress should be defined something like:
$$\int_{\Gamma}\tau_w\cdot v\,d\Gamma = \int_{\Gamma}(T - (T\cdot n) * n)\cdot v\,d\Gamma$$
where $$\tau_w$$ is the wall shear stress vector, v a test function and $$(T - (T\cdot n) * n)$$ the expression for the wall shear stress.
Trying to solve this with:
assemble(Lt = (1 / FacetArea(mesh)) * inner(w, Tt) * ds
(where Tt is the expression for the wall shear stress), only works for a zeroth-order discontinuous Galerkin method, as the mass matrix of this function space is a diagonal matrix where the non-zero entries are equal to the FacetArea. Hence, assembling this vector can be viewed as multiplying the expression for the wall shear stress (inner(w, Tt) * ds) by the inverse of the mass matrix (1 / FacetArea(mesh)). Clearly this only works when the mass matrix is a diagonal matrix. Hence, for all other function spaces this method is simply wrong. Moreover, using the DG0 function space is not the ideal in my opinion, as the nodes of this space are not located on the wall but in the center of the tetrahedrons, hence, your wall shear stress approximation will be slightly off. This effect might be neglected for very small boundary elements, but why would you accept this inaccuracy if you have all information available to get a better approximation (or maybe get the same inaccuracy with much larger elements).
The way I can now compute the wall shear stress correctly is by using the following code (Note: my TrialFunction is shear_stress, my TestFunction is w. Both are defined in the function space VectorFunctionSpace(mesh, 'CG', velocity_degree - 1))
# Define the stress tensor and vector
sigma = 2 * nu * sym(nabla_grad(uv_tmp))
T = sigma * n
# Define the wall shear stress
Tt = T - inner(T, n) * n
Lt = inner(w, Tt) * ds
# Define a left hand side for projection
# CONSTANT throughout program if mesh does not change
a = inner(shear_stress, w) * ds
A_wss = assemble(a, keep_diagonal=True)
A_wss.ident_zeros()
# Assemble the left hand side
# Should assembled each time
b = assemble(Lt)
# Compute the wall shear stress
solve(A_wss, shear_stress_.vector(), b)
Note that in assembling the A matrix I specify keep_diagonal=True and A_wss.ident_zeros(). This is necessary because the expression for the wall shear stress is only defined on the boundary domain. Consequently, the problem is ill-posed: only some columns of the mass matrix have entries (boundary nodes), whereas most parts don't. Hence, for the empty columns we put a 1 on the diagonal with ident_zeros(), which makes the matrix invertible (For a better explanation see again the answer to my other question). This procedure shouldn't affect the solution of solve(A_wss, shear_stress_.vector(), b).
