I try to compute the trace of a DG0 vector field on a mesh (a surface mesh in $R^3$ or a domain in $R^2$, say). Let's forget for a moment that without further assumptions there is no well-defined trace operator.
If my vector field x is an element of VectorFunctionSpace(mesh, 'DG', 0) then by its trace I mean an element in VectorFunctionSpace(bnd, 'DG', 0) which is the tangential component of x (here bnd = BoundaryMesh(mesh, 'exterior') denotes the boundary).
I'm not quite sure how to achieve this. I've tried something like this:
Apparently it's possible to do something like
...
Vb = VectorFunctionSpace(bnd, 'DG', 0)
xb = interpolate(x, Vb)
...
at least if I do x.set_allow_extrapolation(True). See for instance this link, fenicstools provides an even more general version for inter-space interpolation. This gives me
the restriction of x to the boundary. Mathematically - in the smooth world - this would be a section in $TM \mid_{\partial M}$ whereas the trace is in the subbundle $T \partial M$.
Anyway, my next step would be to obtain the facet normal n as an element in Vb and subtract the projection $\langle n, x_b \rangle n$ from xb, and this is where I fail.
The following naive approach appears straight-forward to me, but it does not work:
n = FacetNormal(mesh)
n = project(n,Vb) # Fails
Q: How can I obtain the tangential projection of x as an element in Vb?
Thanks
