0
Integrating DiracDelta in 3D in Cartesian coordinates works just fine i.e. gives vecf[{x, y, z}] vecr = {x, y, z}; vecrp = {xp, yp, zp}; Assuming[(vecr | vecrp | vecf[_]) ∈ Vectors[3, Reals], Simplify[Integrate[DiracDelta[-xp + x] DiracDelta[-yp + y] DiracDelta[-zp + z]*vecf[vecrp], ...