Obtain velocity normal to boundary

2

As the FacetNormal as implied by its name, is defined on each Facet, and not for each degree of freedom, you have to approximate the normal at the degrees of freedom (for instance at corners).

For your particular case, the problem is that you are trying to interpolate an ComponentTensor, that is an expression using ufl syntax. To work around this, I would suggest using an UserExpression as follows:

n = FacetNormal(mesh)
# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)

# Approximate facet normal in a suitable space using projection
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v_ = TestFunction(V)
a = inner(u_,v_)*ds
l = inner(n, v_)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
nh = Function(V)

solve(A, nh.vector(), L)

class normal_u(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        
    def eval(self, values, x):
        n_eval = nh(x)
        u_eval = u(x)
        un = (u_eval[0]*n_eval[0] + u_eval[1]*n_eval[1])
        values[0] = un * n_eval[0]
        values[1] = un * n_eval[1]

    def value_shape(self):
        return (2,)

u_boundary.interpolate(normal_u())