  Density discontinuities cannot be precisely modelled in standard formulations
of smoothed particles hydrodynamics (SPH) because the density field is defined
smoothly as a kernel-weighted sum of neighbouring particle masses. This is a
problem when performing simulations of giant impacts between proto-planets, for
example, because planets typically do have density discontinuities both at
their surfaces and at any internal boundaries between different materials. The
inappropriate densities in these regions create artificial forces that
effectively suppress mixing between particles of different material and, as a
consequence, this problem introduces a key unknown systematic error into
studies that rely on SPH simulations. In this work we present a novel,
computationally cheap method that deals simultaneously with both of these types
of density discontinuity in SPH simulations. We perform standard hydrodynamical
tests and several example giant impact simulations, and compare the results
with standard SPH. In a simulated Moon-forming impact using $10^7$ particles,
the improved treatment at boundaries affects at least 30% of the particles at
some point during the simulation.
