I’m using this kind of technique to avoid having to write an expression both in ufl (eg for generating the exact right hand side f) and numpy (eg for pointwise evaluations):

What I cannot get to work is how to do this if u_exact is a vector-valued function. I’m trying to be clever with mixing np.stack and ufl.as_vector but without success. Is there an example somewhere I can look at?

The fem.Expression mechanism would be a solution if your points are defined via the same coordinates on the reference element. I am not sure if something is similar with points defined by their global coordinates