The theorem states that the surface integral of a vector field times the unit normal equals the volume integral of the divergence of the same vector field. As this is true for \emph{any} vectorfield, we can choose a vectorfield with special properties. Assume that $\mathbf{F}$ has constant divergence, that is $\nabla \cdot \mathbf{F} = c$. Then the volume integral becomes
\begin{equation}
  \int_V \nabla \cdot \mathbf{F}\ d\tau=\int_V c \ d\tau=c\int_V d\tau=cV
\end{equation}
We see that we can compute the volume (times a constant) with this integral. And equation \ref{eq:divergencetheorem} links this integral to the surface integral. After careful consideration we choose $\mathbf{F} =\bar{\mathbf{F}}=(x,y,z)$, and calculate the divergence
\begin{eqnarray}
  \nabla \cdot \bar{\mathbf{F}}&=& \\
  \nonumber & = & \nabla \cdot (x, y, z) \\
  \nonumber & = & \left     (\frac{\partial}{\partial{x}},\frac{\partial}{\partial{y}},\frac{\partial}{\partial{z}}   \right ) \cdot   (x,y,z)\\
  \nonumber & = &\frac{\partial{x}}{\partial{x}} +     \frac{\partial{y}}{\partial{y}}+\frac{\partial{z}}{\partial{z}}\\
  \nonumber & = & 3
\end{eqnarray}
We note that our related surface integral will compute $3V$ ($c=3$) and finally we develop a discrete formula by approximating the integral as a Riemann sum over each face:
\begin{eqnarray}
  \nonumber
  3V &=& \int_V \nabla \cdot \bar{\mathbf{F}}\ d\tau = \int_S \bar{\mathbf{F}} \cdot \mathbf{n}\ dA\\
  & \approx&\sum_{i \in S} \bar{\mathbf{F}}(f_i) \cdot \mathbf{n}(f_i) A(f_i)
  \label{eq:volume}
\end{eqnarray}
Here we used $\bar{\mathbf{F}}(f_i)$ to denote the vector field evaluated at the $i$:th face $(f_i)$, and similarly for the normal and area. At this point  we have not yet decided  where to evaluate the vector field, it is only specified that it should be done somewhere on the face. Remember that the face normal and area are constant over the face. The vector field is not though; but it can be shown that any point (on the face) will do in the limit. We chose the centroid of the triangle giving rise to
\begin{eqnarray}
  3V = \sum_{i \in S} \frac{(\mathbf{v}_1+\mathbf{v}_2+\mathbf{v}_3)_{f_i}}{3} \cdot \mathbf{n}(f_i) A(f_i),
  % V = \sum_{i \in S} (\mathbf{v}_1 +\mathbf{v}_2 + \mathbf{v}_3)_{f_i}  \mathbf{n}_i \cdot A_i
\end{eqnarray}
where $\mathbf{v}_1$, $\mathbf{v}_2$ and $\mathbf{v}_3$ are the vertices of the $i$:th triangle. The error of this approximation decreases as the area of the largest triangle goes to zero. Try to compose a formula that evaluates $\bar{\mathbf{F}}$ at other points
as well and see if you can improve on the results.