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).

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

+\int C\delta(x-y)u(x)v(x)\,dx

to your variational form, where u and v are trial and test functions, then you can do something like

where 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.

I also have the same problem with a large number of point sources in nonlinear problem. Could you please give a simple example about how to use PointSource class with nonlinear problem? Since as far as I remember, you need to do