How is Gauss's divergence theorem related to Gauss's law? I have been reading a bit about vector fields and how it relates to the EM field. However, I am confused about the nature of Gauss's law. As far as I understand, it is supposed to relate flux of the lines over a surface to the net charge inside the closed surface. But how is that relationship achieved? Gauss's divergence law allows one to find the divergence through a closed surface, but how is that connected to the charge of the surface? 
 A: Mathematical Explanation
Gauss's divergence law states that
$$\mathbf{\nabla} \cdot \mathbf{E} = \frac{\rho}{\epsilon_0}  $$
So, let's integrate this on a closed volume $V$ whose surface is $S$, it becomes
$$ \iiint\limits_{V(S)}\nabla \cdot \mathbf{E}\,\textrm{d}V = \frac{Q}{\epsilon_0} $$
where $Q$ is the total charge in $V$. However, Green-Ostrogradski theorem states that
$$ \iiint\limits_{V(S)} \nabla \cdot \mathbf{F}\,\textrm{d}V = \iint\limits_{S}\mathbf{F}\cdot\mathbf{\textrm{d}S}$$
for any field $\mathbf{F}$, so in particular
$$ \iiint\limits_{V(S)} \nabla\cdot\mathbf{E}\,\textrm{d}V = \iint\limits_{S}\mathbf{E}\cdot\mathbf{\textrm{d}S} = \Phi_S(\mathbf{E})$$
where $\Phi_S(\mathbf{E})$ is the flux of $\mathbf{E}$ through $S$. Finally, we got
$$ \Phi_S(\mathbf{E}) = \frac{Q}{\epsilon_0}$$
which is Gauss's theorem.
Intuitive Explanation
The previous explanation demonstrates the link between Gauss's divergence law and its theorem, yet we don't really understand why it works. However, once you've understood what the divergence of a field is, it will appear easy to understand.
Let's start from a possible definition of the divergence of a field $\mathbf{F}$:
$$ \nabla \cdot \mathbf{F} = \lim\limits_{\delta\tau\rightarrow 0} \frac{\iint\limits_{\delta S}\mathbf{F}\cdot\mathbf{\mathrm{d}S}}{\delta \tau}$$ 
where $\delta \tau$ is the volume delimited by $\delta S$. Well... this definition isn't really simple. Let's rewrite it as a physicist would do:
$$\nabla \cdot \mathbf{F}\, \delta \tau = \Phi_{\delta S}(\mathbf{F})$$
where $\delta \tau$ is small enough.... So here, we can see that the divergence of $\mathbf{F}$ is how much of $\mathbf{F}$ goes out of a small volume per unit time. So, if we integrate this over a big volume, we will get what comes out of this big volume per unit time (since when we cut this big volume into many small volumes, everything that goes out of a small volume goes in one of its neighbours, except for the extreme small volumes). But what goes out of this big volume is $\Phi_S(\mathbf{F})$, so we can say that
$$ \Phi_S(\mathbf{F}) = \sum\limits_{\delta \tau \in V} \Phi_{\delta S}(\mathbf{F}) = \sum\limits_{\delta\tau \in V} \nabla\cdot\mathbf{F}\,\delta \tau = \iiint\limits_{V}\nabla\cdot\mathbf{F}\,\mathbf{\textrm{d}V}$$
To a possible understanding of Gauss's law
Everyone starts electrostatic with Coulomb's law, that is two charges $q_A$ and $q_B$ with a distance $r$ between them undergo a force
$$ F_e = \frac{q_A q_B}{4\pi \epsilon_0 r^2}$$
Now, let's consider a charge $q$ in $O$. It creates a field $$\mathbf{E} = \frac{q}{4\pi\epsilon_0 \| \mathbf{OM} \|^3} \mathbf{OM}$$
Consider a small surface $\textrm{d}^2 S = r^2\textrm{d}\theta\textrm{d}\phi$ in $M$, where $\|OM\| = r$. The flux of $\mathbf{E}$ through $\textrm{d}^2 S$ is $\textrm{d}^2\Phi = \mathbf{E}\cdot\textrm{d}^2\mathbf{S}$ where $\textrm{d}^2\mathbf{S} = \frac{\textrm{d}^2S}{r}\mathbf{OM}$. This finally gives
$$ \textrm{d}^2\Phi = \frac{q}{4\pi\epsilon_0 r^3}\mathbf{OM}\cdot\frac{r^2\textrm{d}\theta\textrm{d}\phi}{r}\mathbf{OM} = \frac{q}{4\pi\epsilon_0}\textrm{d}\theta\textrm{d}\phi$$
Now, let's consider a closed surface $S$. If $O$ is in this surface, then
$$ \Phi = \iint\limits_S \frac{q}{4\pi\epsilon_0}\textrm{d}\theta\textrm{d}\phi = \frac{q}{4\pi\epsilon_0} \int\limits_{\theta = 0}^{2\pi}\int\limits_{\phi = 0}^{\pi} \textrm{d}\theta\textrm{d}\phi = \frac{q}{\epsilon_0}$$
Now, focus on the case were $O$ is outside $S$. Let $S_{+}$ (resp. $S_{-}$) be the part of $S$ where $\textrm{d}^2\Phi \geq 0$ (resp. $\textrm{d}^2\Phi < 0$), then
$$ \iint\limits_{S} \textrm{d}^2\Phi = \frac{q}{4\pi\epsilon_0}(\iint\limits_{S_{+}}\textrm{d}\theta\textrm{d}\phi  - \iint\limits_{S_{-}}\textrm{d}\theta\textrm{d}\phi) $$
However, since $S$ is a closed surface, a drawing can convince you that the preceding sum is null (the solid angle of $S_{-}$ is the same as the one of $S_{+}$...). Then, when $O$ is outside of $S$, $\Phi$ = 0.
Since the electrostatic field is linear, we can generalize this idea with the following equation:
$$ \Phi_S(\mathbf{E}) = \frac{Q}{\epsilon_0} $$
where $Q$ is the net charge inside $S$.
