I'll call epsilon e, rho p, phi f, and let d/dp stand for the partial derivative with respect to p, etc.
Since D1 is entirely radial:
(grad)(dot)D1 = 1/p * d/dp[p*D1]
= e/p * d/dp[p^3] = 3pe
(grad)(dot)D2 = 1/p * d/dp[p*D2(p)] + 1/p * dD2(f)/df
where D2(p) is the p component of D2, etc.
= e/p * sin^2(f) * d/dp[p^3] + 3epz/p * d/df[cos^2(f)]
= 3ep*sin^2(f) + 3ez*sin(2f)
I'm assuming D is the displacement vector, so
(grad)(dot)D = (charge density)
Thus I would agree with you. The two charge densities are not the same.