How do I add passive tracers to the Navier Stokes equations

I see examples solving Navier Stokes or Stokes equations, and I would like to add passive tracers into the fluid. But I can’t find such examples on the passive tracers and the plots of them.

Naively, I think I can define some real[int] Px, Py, to record the position of each particle in the fluid. And I can find out which triangle this (Px, Py) belong to, by looping through all the triangles, and find out the indices for vertices, and find out velocity values at those vertices, and then do an interpolation.

But this seems very inefficient. Is there something already programmed in the package that I can use?

I just found out that one can use something like
tmp = v(P[0],P[1])
P[0] += u(P[0],P[1]);
P[1] += tmp;
to update the position of passive tracer particles. However, how can I determine if a particle is inside the region with certain label? When the particle hits the boundary, how can I determine which boundary it hits, and the unit outer normal at that location?