I've been trying for days to solve the following identity, but feel like I really need a hint on how to start since I'm not yet much into how to work with submanifolds.
Let $U\in \mathbb R^n $ and $f: U \to \mathbb R$ continuously differentiable with $df_x \neq 0$ for all $x\in U$ so that, for each $t \in \mathbb R$, $M_t := f^{-1} (\{ t \})$ is a $(n-1)$-dimensional submanifold of $\mathbb R^n$. Then, for every non-negative, $\mathcal B(U)$-measureable map $g: U \to \mathbb R$:
$$\int_U g \,d\lambda_n = \int_\mathbb R \int_{M_t} \frac{g(x)}{\lVert \nabla f(x) \rVert} \omega^{M_t}(dx)\,dt,$$ where $\omega^{M_t}$ is the surface measure on $M_t$ defined by $$\omega^{M_t}(A) := \omega^{\phi^*\sigma}(\phi ^{-1} (A)) := \int_{\phi ^{-1}(A) }\sqrt{\det (\phi^*\sigma)_x} \, \lambda_n(dx)$$ for a parametrization $\phi$ and the standard inner product $\sigma$.
My attempt so far:
My first idea was to apply the implicit function theorem to find a unique function with $f(x, g(x)) = t$ for each $x$ and some neighbourhood $U_x$ and then apply Fubini to split the left hand side in two integrals. However, I lack technical knowledge to do that.
My other attempt was to use a characterization of submanifolds:
Fix $x \in M_t$. We can find a neighbourhood $x \in U_x \subseteq U$ and a $C^1$-diffeomorphism $F: U_x \to W_x$, $W_x \subseteq \mathbb R^n$ open neighbourhood of $0$. However, in this case, I don't know of suitable ways to split up the integral. Any help appreciated.
Using a partition of unity, we assume that $V$ is a sufficiently small neighborhood of $p\in U$ such that, up to renumbering of the coordinates, $\partial_n f(p)\neq 0$. As in the implicit function theorem, define the change of coordinates $\Phi(x)=(x',f(x))=(y,t)$, where $x'=(x_1,\ldots,x_{n-1})=y$. The Jacobian gives $dydt=|\partial_n f|\,dx$, hence
$$\int g\,dx = \int_{\mathbb{R}}\int g(\Phi^{-1}(y,t))\frac{dy}{|\partial_n f(\Phi(y,t))|}dt $$
Notice that $\Phi^{-1}(y,t)=(y,h(y,t))$ and $f(y,h(y,t))=t$, hence $M_t$ is the graph of the function $y\mapsto h(y,t)$.
Since the area element of $M_t$ is $\sqrt{1+|\nabla_yh(\cdot, t)|}dy$, to conclude our argument we should prove that $\sqrt{1+|\nabla_yh(\cdot, t)|}/|\nabla f|=1/|\partial_n f|$, or $|\nabla f|^2=|\partial_nf|^2+\sum_i|\partial_nf|^2|\partial_ih|^2$, but this follows from $f(y,h(y,t))=t$ by taking the $i$-derivative.
This is the Gelfand-Leray formula and admits generalizations to vector functions $f$ instead of simply scalar, this is called integration over the fibre.