Integrating over a shell in $d$-dimensional space. In M. E. Tuckermann's book Statistical Mechanics: Theory and Molecular Simulation, par. 3.2, the following approximation is made:
$$\int_{c<f(\vec x)< c+\epsilon} d^D x \simeq \epsilon \int_{\mathbb R^D} d^Dx \ \delta \left(f(\vec x) -c\right)$$
Basically, the integral over a $D$ dimensional shell is replaced by the integral over a $(D-1)$ dimensional area times the thickness $\epsilon$ of the shell.
I would like to know: when is this approximation good? I suspect that the above relation becomes exact in some limit, maybe the $D \to \infty$ limit, of course if $f$ is "well-behaved" enough.
To have a concrete example, we can consider a spherical shell in $D$ dimensions:
$$\int_{R<\| \vec x \|<R+\epsilon} d^D x = V_D(R+\epsilon) - V_D(R)$$
where
$$V_D(R) = \frac{R^D \pi^{D/2}} {\Gamma(D/2+1)}$$
so that
$$\int_{R<\| \vec x \|<R+\epsilon} d^D x = \frac{ \pi^{D/2}} {\Gamma(D/2+1)} \left( (R+\epsilon)^D - R^D \right)$$
while the other integral would be
$$\int_{\mathbb R^D} d^Dx \ \delta (\| \vec x \| -R)$$
but I'm not sure about how to evaluate this.

Update 
I managed to show this in the particular case of an hypersphere. 
From 
$$V_{shell} =\int_{R<\| \vec x \|<R+\epsilon} d^D x = \frac{ \pi^{D/2}} {\Gamma(D/2+1)} \left( (R+\epsilon)^D - R^D \right) = \frac{ \pi^{D/2} R^D} {\Gamma(D/2+1)} \left[\left(1+\frac{\epsilon}{R}\right)^D - 1 \right]$$
we do the following approximation for $\epsilon \ll R$
$$\left(1+\frac \epsilon R \right)^D \simeq 1 + D \frac{\epsilon}{R}$$
from which
$$V_{shell} \simeq \frac{ \pi^{D/2} R^{D-1}} {\Gamma(D/2+1)} \ \epsilon D$$
Now, since
$$\Gamma(D/2+1) = \frac D 2 \Gamma\left(\frac D 2 \right)$$
we have 
$$V_{shell} \simeq \epsilon \frac{2 \pi^{D/2} R^{D-1}}{\Gamma(D/2)} = \epsilon S_D(R) $$
where $S_D(R)$ is the surface area of the hypersphere. 
So, as all the answers and comments up until now point out, the real condition seems to be indeed that $\epsilon$ must be small.
 A: If $g: \Bbb R^D \to \Bbb R$ is a surjective submersion ("submersion" meaning in this case that $\| \nabla g (x) \| \ne 0 \ \forall x \in \Bbb R^D$), then we have the following "coarea" formula for the "$\delta$ function"
$$\int \limits_{\Bbb R ^D} \varphi (x) \delta (g(x)) \ \Bbb d ^D x = \int \limits _{g^{-1} (0)} \frac {\varphi(x)} {\| \nabla g (x) \|} \ \Bbb d \sigma ,$$
where the right-hand side is to be understood as a hyper-surface integral.
Choosing $\varphi = 1$ and $g(x) = f(x) - c$, the integral in the right-hand side of the investigated equality is
$$\int \limits_{\Bbb R ^D} \delta (f(x) - c) \ \Bbb d ^D x = \int \limits _{f^{-1} (c)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma .$$
To study the left-hand side, we shall use a convenient variant of the smooth coarea formula (Isaac Chavel, "Riemannian Geometry - A Modern Introduction", 2nd edition, Cambridge University Press, 2006, Exercise III.12.d, page 160; - I have taken the liberty to slightly modify it in order to make it more relevant to our needs):
$$\int \limits _{\Bbb R ^D} \varphi (x) \ \Bbb d ^D x = \int \limits _{\Bbb R} \Bbb d t \int \limits _{f^{-1} (t)} \frac 1 {\| \nabla f (x) \|}\varphi \Big| _{f^{-1} (t)} (x) \ \Bbb d \sigma ,$$
where $f$ is as above and $\varphi$ is either measurable and positive, or Lebesgue-integrable.
Let us apply this taking $\varphi (x) = \begin{cases} 1, & x \in f^{-1} \big( (c,c+\varepsilon) \big) \\ 0, & \text{otherwise} \end{cases}$, i.e. the characteristc function of $f^{-1} \big( (c,c+\varepsilon) \big)$. Then
$$\int \limits _{c < f(x) < c + \varepsilon} \ \Bbb d ^D x = \int \limits _{f^{-1} \big( (c,c+\varepsilon) \big)}  \ \Bbb d ^D x = \int \limits _{\Bbb R ^D} \varphi (x) \ \Bbb d ^D x = \\
\int \limits _{\Bbb R} \Bbb d t \int \limits _{f^{-1} (t)} \frac 1 {\| \nabla f (x) \|}\varphi \Big| _{f^{-1} (t)} (x) \ \Bbb d \sigma = \int \limits _c ^{c + \varepsilon} \Bbb d t \int \limits _{f^{-1} (t)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma ,$$
where we have used that if $\varphi \Big| _{f^{-1} (t)} (x) = \begin {cases} 0, & t \notin (c, c + \varepsilon) \\ 1, & t \in (c, c + \varepsilon)\end{cases}$, which allowed us to restrict the domain of integration to just $(c,c+\varepsilon)$ in the outer integral and to replace $\varphi$ by $1$ in the inner one.
Now, apply the mean value theorem for integrals to obtain some $t_\varepsilon \in (c, c + \varepsilon)$ such that
$$\int \limits _c ^{c + \varepsilon} \Bbb d t \int \limits _{f^{-1} (t)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma = \varepsilon \int \limits _{f^{-1} (t_\varepsilon)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma .$$
Finally, use an argument based upon continuity to obtain that, if $\varepsilon$ is small, then $c+\varepsilon$ is close to $c$, and since $t_\varepsilon$ sits between them, $t_\varepsilon$ is close to $c$, too, so that
$$\int \limits _{f^{-1} (t_\varepsilon)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma \simeq \int \limits _{f^{-1} (c)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma .$$
Putting everything together, we obtain that for small $\varepsilon$
$$\varepsilon \int \limits_{\Bbb R ^D} \delta (f(x) - c) \ \Bbb d ^D x = \varepsilon \int \limits _{f^{-1} (c)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma \simeq \varepsilon \int \limits _{f^{-1} (t_\varepsilon)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma = \\
\int \limits _c ^{c + \varepsilon} \Bbb d t \int \limits _{f^{-1} (t)} \frac 1 {\| \nabla f (x) \|} \ \Bbb d \sigma = 
\int \limits _{c < f(x) < c + \varepsilon} \ \Bbb d ^D x .$$
A: Consider a normalized step function $\Pi_\eta(x)$ between 0 and $\eta$:
$$
\Pi_\eta(x) \equiv \begin{cases} 1/\eta & 0 \leq x \leq \eta \\ 0 & \text{otherwise.} \end{cases}
$$
Note that by construction, $\int \Pi(x) dx = 1$ over any region containing the interval $[0, \eta]$.
We can write the integral on the left-hand side of the equation as
$$
\frac{1}{\epsilon} \int_{c \leq f(\vec{x}) \leq c + \epsilon} d^Dx = \int_{\mathbb{R}^D} \Pi_\epsilon( f(\vec{x}) - c) \, d^Dx.
$$
Note that the factor of $1/\epsilon$ on the left-hand side is necessary to balance out the $1/\epsilon$ in the definition of $\Pi_\eta$.  
Finally, we note that in a distributional sense,
$$
\lim_{\eta \to 0} \Pi_\eta(x) = \delta(x).
$$
Thus, in the limit as $\epsilon \to 0$, we have
$$
\lim_{\epsilon \to 0} \left[ \frac{1}{\epsilon} \int_{c \leq f(\vec{x}) \leq c + \epsilon} d^Dx  \right] = \int_{\mathbb{R}^D} \delta( f(\vec{x}) - c) \, d^D x
$$
from which it follows that in this limit, 
$$
\int_{c \leq f(\vec{x}) \leq c + \epsilon} d^Dx \approx \epsilon \left[\int_{\mathbb{R}^D} \delta( f(\vec{x}) - c) \, d^D x \right] + \mathcal{O}(\epsilon^2),
$$
as desired.
This wasn't the entirety of the question you were asking, though.  The real question was essentially, "how large are the $\mathcal{O}(\epsilon^2)$ terms in the above expression?"  The best answer I can give for this rather less rigorous, I'm afraid.  The formula we're looking at basically says that
$$
( \text{volume of shell of thickness $\epsilon$ surrounding a hypersurface} ) \approx \epsilon (\text{area of hypersurface}).
$$
(If you don't see why your integrals are equivalent to this, let me know and I'll edit to elaborate.)  This is often thought of saying that the volume of "paint" needed to cover a surface to a thickness $\epsilon$ can be well-approximated by assuming that the surface is flat.  
Any correction of this formula must therefore be related to the fact that the "surface" we're covering—in this case, the hypersurface $M$ given by the level set $f(\vec{x}) = c$–is curved.  I suspect, in fact, that the second-order correction is something like
$$
\epsilon^2 \int_M (\text{curvature}) d^{D-1} x,
$$
where the integrand is related to the curvature of the surface $M$ as embedded in $\mathbb{R}^D$.  
I suspect that the curvature quantity in question is (up to a numerical factor) $H$, the mean curvature of the surface;  but I haven't been able to prove it just yet.  If this is true, then the approximation is "good" so long as $\epsilon$ is significantly smaller than the radii of curvature of $M$.
