Delta function derivative

Hi,
I would like to include in the variational formulation of a linear problem (let’s assume for simplicity and definiteness that is just the Poisson equation on a square with zero Dirichlet boundary conditions on all the edges) a source term containing the gradient of the delta function to represent a dipole source.

The variational formulation is

\int dx \nabla v(x) \cdot \nabla u(x) = - 4\pi \int dx ( p \cdot \nabla \delta(x-x_0)) v(x) \equiv 4\pi \int dx p \cdot \nabla v(x) \delta (x-x_0) = 4\pi p \nabla v(x_0)

with $$p$$ being a constant vector with the right dimension and x_0 being the position of the source.

Up to now I solved the problem by defining a mollified version of the delta function or its gradient and implementing in FENICS either the second or third member of the above equation with the mollified delta.
Is it possible to implement the same thing using the class PointSource?

Thanks.
Iacopo

One hack might be to use several PointSources in a tiny finite difference stencil approximating \nabla v. It’s still an approximation, but, unlike the regularized delta, which needs to span a few elements for Gaussian quadrature to be effective, the finite difference stencil of PointSources could be shrunk down arbitrarily (until, of course, you eventually run into floating-point precision issues). There might be some problems when the stencil happens to straddle an interior facet, but that also sort of comes with the territory of applying Dirac distributions to discontinuous functions.