Hi!
Answering my own question, as shown here in 2012, we can get the normal component by solving the projection manually. The important step is to set the diagonal of the matrix of all interior nodes to 1.
# d3model.V is a 3D VectorFunctionSpace
phi = TestFunction(d3model.V)
du = TrialFunction(d3model.V)
n = d3model.N # normal vector
u = d3model.U3 # computed velocity vector
dLat = d3model.dLat_d # lateral boundary of interest
u_t = u - dot(u, n)*n # tangential component of velocity
a_n = inner(du, phi) * dLat
L_n = inner(u_t, phi) * dLat
A_n = assemble(a_n, keep_diagonal=True)
B_n = assemble(L_n)
A_n.ident_zeros() # set interior dofs from 0 to 1
u_t = Function(d3model.V)
solve(A_n, u_t.vector(), B_n)
As I really only needed the normal component to determine the tangential component, this is what I've computed and shown here, after using LagrangeInterpolator to interpolate onto a SubMesh of the 3D BoundaryMesh :