Is divergence of electric field always going to give you Dirac delta function? We all know that when we calculate the divergence of point charge at origin, it turns out that it's zero at all points except origin and infinite at origin, which is called Dirac delta function. refer here
$$ \nabla \cdot \mathbf E = 4\pi\delta^3(r).$$
Now let's consider a continuous charge distribution in space. We are interested in finding divergence of field at some point. For this purpose we follow the definition of divergence at a point:
$$ \nabla\cdot\mathbf E =\lim_{\Delta V \to 0} \frac{1}{\Delta V}\iint_S\mathbf E\cdot \hat{\mathbf n}dS .$$
Now if we make a sphere  and let the volume shrink arbitrarily to zero, it's obvious that continuous charge distribution will not be continuous any more and we should get the same result as we get in case of point charge. But it's not the case. Please explain where I'm wrong.
 A: 
Now let's consider a continuous charge distribution in space.

Contradicts with

it's obvious that continuous charge distribution will not be continuous anymore

That's the problem in your reasoning.
Why? This is a very usual confusion when you begin studying. Physics do not deal with reality, but with models. You are not putting a gaussian surface to measure a concrete atom of yoru concrete and real system. No! What you're doing is: "since this object is very complicated, because it has more than $10^{23} atoms, let's consider an easier thing". "Let's replae this terrible object by a comfortable IDEAL mathematical obejct". "Let's forget about that complex reality and let's imagine a nice perfectly smooth surface in which I can apply all my mathematical tools".
That's what we do in physics. We work with maths, maths are our tool, but we can only apply maths to mathematical objetcs. The point is to choose the objects so that their results are as similar to reality as possible.
In other words, we make models that behave approximately the same as reality.
But you should not be mixing the real world with our idealized models.
It's like "No, the ball I threw does not describe a parabola because the atoms are quantum objects that...." Wait, man, you're not solving a group of atoms of a ball, you're solving how a point mass would move in absence of friction and many other approximations. The result is a parabola. Reality is similar to that parabola? Yes, it is ... close enough! If then a strong wind comes, your model might be unsuitable for your situation, so the resolution will be wrong. I hope you get my point: we the physicist must be aware of what level of approximation are applying, and we must be very cautious on what are the validity limits of our model.
So, back to your problem, if you're assuming that you have a continuous medium, then forget about reality, your'e solving a continuous medium. Is it valid? It is valid as long as it matches with what you measure. If you zoom in until you see atoms, your model is probably unsuitable to describe reality, as the electric field varies a lot between atoms. But for sufficiently far, it is a good model.
A: The answer depends on how closely you want to look and what approximations you're willing to make.
At some point, what you say is true: if you look close enough, things are made of point charges, so within a small enough volume, the divergence of the electric field would just be a delta function. But only in that volume: over all of space, the charge density would be a sum of delta functions
$$\rho(\mathbf{r}) = \sum_i q_i \delta(\mathbf{r}-\mathbf{r}_i),$$
where $i$ indexes all the charges present. If we don't look too close, we can replace this charge density by an approximate version $\bar{\rho}$, which is the average of $\rho$ over a small volume surrounding $\mathbf{r}$:
$$\bar{\rho}(\mathbf{r}) = \frac{1}{V_\mathbf{r}} \int_{V_\mathbf{r}} \rho(\mathbf{r}')\, d^3\mathbf{r}' = \frac{1}{V_\mathbf{r}} \sum_{i \in V_\mathbf{r}} q_i = \frac{Q_{V_\mathbf{r}}}{V_\mathbf{r}}.$$
That is, the average charge density is the charge contained in a small volume around the point $\mathbf{r}$ divided by the volume. This works as long as the volume is neither too large nor too small; if it's too large then the approximation is too rough, and if it's too small it stops being continuous.
But... this is only the case if you look closely but not too closely, and only in some materials. Because this is when quantum mechanics steps in, and reminds us that even though particles are point charges, they don't have well defined positions, so in practice they work a bit like continuous distributions over the extent of their wavefunctions. How important this is depends: if you have a bunch of fixed ions next to each other, or some free electrons, treating them as point charges might work well. But in a crystal, and specifically in a metal, some electrons are kind of spread out over the whole material, and there they really act like a continuous distribution instead of a collection of points.
So like I said, it depends on the approximations you want to make, and what scales and situations you're considering.
A: The argument of delta function comes in picture only when the electric field goes by 1/r^2.

let's consider a continuous charge distribution in space

when you are taking a continuous charge distribution,it will remain continuous to arbitrary precision. In the view of a mathematical construct, as a continuous function always remain continuous, No matter, by how much you magnify it, so is continuous distribution. So there is no sense of a single discrete particle in a continuous distribution. a continuous function is always continuous to arbitrary precision.

Now if we make a sphere and let the volume shrink arbitrarily to zero

when you are saying you are going to arbitrary zero, you are certainly not practical, so lets go mathematical,
mathematically, since there is no notion of single charged particle, no point of discreteness, the field not necessarily varies over 1/r^2. so Dirac delta is not going to Pop up.
By practical I mean , practical considerations of V-> 0
Practically, we take a lot more particle. You just add one more charged particle.
Now you have 2 charged particle, and just imagine it a dipole, now you will get a field varying over 1/r^3.Now also Dirac delta will not come into picture.one can see by calculation that even if one takes one barn by area then also it will contain more than one charged particle (say proton). So the argument of Dirac delta is shallow practically.
say, you are taking a linear charge distribution , when you are taking a volume element for your experiment , symmetrically your volume element itself should be something like line (like in Gauss law). Its no point in saying that I am taking a linear charge distribution and taking flux of deep down on electron, if you are concerned about one electron, then their is no point in making linear,spherical, areal distribution. If you using  the properties of linear or planar distribution then you are not going to get delta function. you should consider at least thousand of electron so that you can atleast say

Oh! its seems planar, or linear

So, the flaw is You first said you took continuous charge (here,you are mathematical),

we calculate the divergence of point charge at origin

Then you took one charge (You are now practical),because mathematically, V-> 0 doesn't mean single charge. So It gave absurd result. Taking one charge is neither possible mathematically nor practically.**
