Surface integral via differential forms I am trying to figure out why the change of variables formulas for surface integrals (for example) is what it is. I was wondering if there is a justification in terms of differential forms. 
So, let $g : V \to \mathbb{R}$ be a smooth function, andlet $S$ be the surface given by $G(V)$, where $G(x, y) = (x, y, g(x, y))$. In my current predicament, $V$ is an open subset of $\mathbb{R}^2$, but supposedly 2 may be replaced with general $n$. Let $F : S \to \mathbb{R}$ be a smooth function. I found the following formula online:
$$
\int_S F\,  dS = \int_V (F\circ G)\sqrt{g_x^2 + g_y^2 + 1}\ dA.
$$
Here $dA = dx \wedge dy$ is the standard area form. I have a few questions about this.
Could someone recommend a rigorous proof of this fact, or perhaps a proof sketch? I am familiar with the change of variables formula when we are considering a change of variables by a diffeomorphism between open sets in $\mathbb{R}^n$, but this seems to fall apart here (since $S$ is not open in $\mathbb{R}^3$). 
I was trying to justify this with the machinery of differential forms. For instance, we know that if we define $\omega = F \, dS$ to be a 2-form on $S$,
$$
\int_S \omega = \int_V G^*\omega.
$$
This is where I get really confused. As far as I know, $G^*\omega$ is defined as follows: 
\begin{align*}
G^*\omega(x, y)(v, w) &= \omega(G(x, y))(dG(x, y)(v), dG(x, y)(w)) \\
&= \omega(G(x, y))(v_1, v_2, g_xv_1 + g_yv_2, w_1, w_2, g_xw_1 + g_yw_2).
\end{align*}
But I have no idea how to proceed from here. How does this end up getting us the square root from above? I see a lot of determinants of transformations, but $dG$ is not a square matrix, so I'm not sure how to get this to appear in this case. 
Any help is much appreciated.
 A: The part you are missing is the coordinate expression for the area form $dS$. To understand this, you need to know more about Riemannian volume forms.
So, let $(M,g)$ be an oriented $n$-dimensional Riemannian manifold. The Riemannian volume form of $M$, which we denote by $\mathrm{vol}_g$, is defined geometrically as the unique $n$-form on $M$ satisfying $$\mathrm{vol}_g(e_1,\ldots,e_n)=1$$ whenever $e_1,\ldots,e_n$ is an oriented orthonormal (local) frame. The local coordinate expression of $\mathrm{vol}_g$ is $$\mathrm{vol}_g=\sqrt{\det g_{ij}}\;dx^1\wedge\ldots\wedge dx^n,$$where $g_{ij}$ is the matrix representation of $g$ in the local coordinates $x^1,\ldots,x^n$.
Now, your area form $dS$ is, in fact, the Riemannian volume form of the surface $S$, equipped with the Riemannian metric induced by the embedding in $\mathbb{R}^3$. The map $G$ is a parametrization of $S$, and the matrix representation of the Riemannian metric with respect to $G$ is given by $$g=\left(\begin{array}{cc}g_x^2+1&g_xg_y\\g_xg_y&g_y^2+1\end{array}\right).$$ Hence, $$\det g_{ij}=g_x^2+g_y^2+1,$$and your desired formula follows by putting all the pieces together.
