Possibly incorrect computation of the wallShearStress field in the FO

Hi! I'm looking in at the implementation of the wallShearStress FO. The actual computation happens in wallShearStress::calcShearStress. For each patch, we do the following (comments added by me):

// Get a ref to the patch field to be computed
vectorField& ssp = shearStress.boundaryFieldRef()[patchi];

// Get normal vectors and their magnitude
const vectorField& Sfp = mesh_.Sf().boundaryField()[patchi];
const scalarField& magSfp = mesh_.magSf().boundaryField()[patchi];

// Get the stress tensor at the boundary
const symmTensorField& Reffp = Reff.boundaryField()[patchi];

// Compute the wall stress by dotting the stress tensor with the (inward) normal.
ssp = (-Sfp/magSfp) & Reffp;

I may be confused, but to me, it seems like that we got the wall stress vector, not the wall shear stress vector. To get the shear stress, we should have the additional operation of subtracting the normal component of ssp:

ssp -= (-Sfp/magSfp) * ( (-Sfp/magSfp) & ssp)

If we take a channel-like geometry with wall normals aligned with y, the resulting shear stress vector field should be exactly 0 in the second component, right? I've tested, and this is not the case.