Is there any efficient way to add many dirac deltas on a variational form?

[nate]

You can apply many point sources at once using this.

Furthermore, you can indeed use the NonlinearProblem class. Implement your own solver similar to here, but apply your point sources when your construct the residual (assuming they don’t constitute a component of the Jacobian).