I’ve used the approach that @nate suggested for large numbers of point forces corresponding to immersed boundaries before, and it’s decently efficient.
I’ll add that, if the point sources do contribute to the Jacobian, you can sometimes get away with a lumped-mass approximation of that contribution. For instance, if you have a mass-matrix-like point contribution
to your variational form, where u and v are trial and test functions, then you can do something like
ps = PointSource(V,y,C)
lmVec = Function(V).vector()
A is the LHS matrix from the rest of the variational form and
PETSc is imported from
petsc4py. This is not the exact Jacobian contribution, but can be good enough to prevent Newton iteration from diverging, and, if you converge with a consistent residual, the solution is unaffected by approximations in the Jacobian.