Keep track of the density of a region of space

I am implementing a contact model that allows particles to overlap. Is it possible to get information about the density of a region in space such that I would allow overlaps only if the density is lower than some given parameter? For example, lets say I have a system with identical particles (spherical, density = rho, radius = R). Suppose 3 particles overlap and occupy a region in space. Is it possible to get the density (3*rho, assuming they add up linearly) of the region. I am trying to build a contact model based by putting a bound on the density of a region. I feel this cannot be achieved simply by modifying/creating a new contact model and requires some fiddling around with other aspects of the code. Any advice on the functions/files/classes I need to look into?