Integrating using surface and volume elements As a Physics Degree undergraduate, I have been forced countless times to use a certain method to integrate over 3D surfaces and volumes, which my lecturers like to call integration through surface and volume elements. I can't stand this method, nor I could ever understand how to do it, especially because it is not mathematically rigorous.
However, at most cases, I'm not able to avoid it. For example - given the electric field of a ring of radius $r$ lying on $z=0$ at the point $(0,0,z_0)$, I am required to find the electric field of a disk of radius $R$ lying on $z=0$ at the point $(0,0,z_0)$. In order to do that, I have to use the aforesaid method - integrate the electric field I'm given with respect to the length element $dr$, from $r=0$ to $r=R$. Of course, doing that would also require to translate the charge density of a length element, to a charge density of a surface element, assuming they're both unifrom.
But - Physics is not my problem here - but the math. And that's why I came here. I tried to understand using this method, but it sometimes works - and sometimes doesn't. I would be glad to know where I'm right, and where I'm wrong.
$(\star)$ Important: The angle $\theta$ in Examples 1,2 is the angle of the polar coordinates. In Examples 3,4, it is the polar angle of the spherical coordinates (meaning it is not the azimuthal one).
$(\star)$ I will denote $\color{green}{Good}$ in green and $\color{red}{Bad}$ in red. Lowercase will be integration variables, and uppercase would be given parameters.
Example 1: Calculating the area of an empty cylinder of radius R and height H
A. With respect to $dz$
Given a perimeter of a ring $2\pi R$, the area of a ring with an infinitesimal height $dz$ would be given by $2\pi Rz$. And then:
$$S=\int\limits_{0}^{H}2\pi Rz\ dz=\color{green}{2\pi RH}$$
A correct answer, gladly.
B. With respect to $d\theta$ 
We know that if we sliced the cylinder vertically, rotating with the angle $\theta$, we would get lines of height $H$ each, multiplied by an infinitesimal width $Rd\theta$. Thus, the surface element would be given by $HRd\theta$. And then:
$$S=\int\limits_{0}^{2\pi}HR\ d\theta=\color{green}{2\pi RH}$$
Again - a good answer.
But: this is where things are going to get ugly.
Example 2: Calculating the volume of a cylinder of radius R and height H
A. With respect to $dr$
We would want to sum cylinders with infinitesimal widths $dr$, thus the volume element would be given by $2\pi H rdr$ (the perimeter of a ring of radius $r$ multiplied by the infinitesimal width $dr$ and height $H$). And then:
$$V=\int\limits_{0}^{R}2\pi Hr\ dr=\color{green}{\pi HR^2}$$
This is of course correct, but:
B. With respect to $d\theta$
We would want to sum the exact same slices we described at B. of Example 1, but now they would also have a width of $R$. Meaning: the volume element would be given by $HR^2d\theta$ (since every rectangle is of dimensions $H \times R$, and we multiply each by an infinitesimal width $Rd\theta$). And now:
$$V=\int\limits_{0}^{2\pi}HR^2\ d\theta=\color{red}{2\pi HR^2}$$
This is bad. I would show  you now 2 more examples - in the case of a sphere and a ball. It doesn't work there either.
Example 3: Calculating the area of a sphere of radius R
With respect to $d\theta$
Given a ring of radius $r$, it can be easily checked to see that, geometrically, $r$ would be given by $R\sin\theta$. The infinitesimal width of such disk, would be now $Rd\theta$, thus the surface element would be given by $2\pi R^2\sin\theta d\theta$. Therefore:
$$V=\int\limits_{0}^{\pi}2\pi R^2\sin\theta \ d\theta=\color{green}{4\pi R^2}$$
Getting optimistic, let's try to calculate the volume of the ball.
Example 4: Calculating the volume of a ball of radius R
A. With respect to $dr$
We would want to sum spheres, of radius $r$ and infinitesimal width $dr$ each. Thus, the volume element would be given by $4\pi r^2 dr$, and then:
$$V=\int\limits_{0}^{R}4\pi r^2 \ dr = \color{green}{\frac{4}{3}\pi R^3}$$
But unfortunately:
B. With respect to $d\theta$
Going again like B. of Example 3, we would want to sum the exact same rings, but now they would be disks with the infinitesimal width $Rd\theta$. The volume element would be given by $\pi (R\sin\theta)^2 Rd\theta$, which leads us to:
$$V=\int\limits_{0}^{\pi}\pi R^3 \sin^2\theta\ d\theta=\color{red}{\frac{1}{2}\pi^2 R^3}$$
I tried using the other elements too: $d\varphi$, for example. the azimuthal angle, which is much more complicated, and also tried other shapes like a cone and even paraboloid. But it just won't work right. It works sometimes - and that's not enough for me, unfortunately. I put many efforts to this post, in order to show you my way of thinking, because that's how I had been taught to do this. But maybe it is not right (it feels like it, for sure).
Thank you very much for reading all this, and I would be very glad to hear your thoughts.
P.S.: I wish I could add pictures, but I don't know any programs that I can use to draw them.
 A: In example 2B, The infinitesimal slices aren't rectangular cuboids, but prims: their heights is $H$ and their base are isosceles triangles with legs of length $R$ and angles at apex being $d\theta$. That means that the volume of each element is $\frac12HR^2d\theta$, not $HR^2 d\theta$. That will give you the correwct result.
In 4B, the infinitesimal are disks with constant width, but their width is $-dz = R\sin\theta d\theta$, not just $Rd\theta$ - that's what causes you to get an incorrect answer.
In general, you need to be careful about what shape the infinitesimal elements will have, and not assume that their volume/surface is given by the simplest possible formula. If you divide your figure in  a way that makes it difficult to calculate the volume of an infinitesimal element, than this particular method of dividing the figure may just not be very useful. Howwver, if you correctly calculate the volume/area of each infinitesimal element, the result should be the same no matter how you divide the figure. If you're getting a wrong result, that's most likely because you calculated the volume/area of infinitesimal element incorrectly.
If you want to make sure that you're calculating it right, it's better to always start with full parametrization of the figure, write the area/volume as an integral over several variables and then integrate them out one by one. The general formulas are:
If the surface is parametrized by $\vec{r} = \vec r(u,v))$ then you can use one of two equivalent formulas:
$$ dS = \left|\frac{\partial\vec r}{\partial u} \times \frac{\partial\vec r}{\partial v}\right| du dv$$
or $$ dS = \sqrt{\left|\det\begin{bmatrix}\frac{\partial\vec r}{\partial u}\cdot \frac{\partial\vec r}{\partial u} & \frac{\partial\vec r}{\partial u}\cdot \frac{\partial\vec r}{\partial v} \\ \frac{\partial\vec r}{\partial v}\cdot \frac{\partial\vec r}{\partial u} & \frac{\partial\vec r}{\partial v}\cdot \frac{\partial\vec r}{\partial v}\end{bmatrix}\right|} du dv$$
(the second formula, although more complicated works for surfaces embeded in a space of any dimension).
for the volume parametrized by $\vec{r} = \vec{r}(u,v,w) = (x(u,v,w),y(u,v,w),z(u,v,w))$ we have
$$ dV = \left|\det\begin{bmatrix}\frac{\partial x}{\partial u} & \frac{\partial x}{\partial v}&\frac{\partial x}{\partial w}\\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v}&\frac{\partial y}{\partial w}\\\frac{\partial z}{\partial u} & \frac{\partial z}{\partial v}&\frac{\partial z}{\partial w}\end{bmatrix}\right| du dv dw$$
or
$$ dV = \sqrt{\left|\det\begin{bmatrix}\frac{\partial\vec r}{\partial u}\cdot \frac{\partial\vec r}{\partial u} & \frac{\partial\vec r}{\partial u}\cdot \frac{\partial\vec r}{\partial v} & \frac{\partial\vec r}{\partial u}\cdot \frac{\partial\vec r}{\partial w}\\ \frac{\partial\vec r}{\partial v}\cdot \frac{\partial\vec r}{\partial u} & \frac{\partial\vec r}{\partial v}\cdot \frac{\partial\vec r}{\partial v} & \frac{\partial\vec r}{\partial v}\cdot \frac{\partial\vec r}{\partial w} \\\frac{\partial\vec r}{\partial w}\cdot \frac{\partial\vec r}{\partial u} & \frac{\partial\vec r}{\partial w}\cdot \frac{\partial\vec r}{\partial v} & \frac{\partial\vec r}{\partial w}\cdot \frac{\partial\vec r}{\partial w}\end{bmatrix}\right|} du dv dw$$ 
To see how it works, let us consider your examples.
Example 1: The surface of a cylinder can parametrized by coordinates $(z,\theta)$ as $\vec r(\theta,z) = (R\cos\theta, R\sin\theta, z)$. The range of the parameters are $\theta\in[0,2\pi]$,$z\in[0,H]$. The infinitezimal area, calculated from one of the two formulas above, turns out to be $dS = Rd\theta dz$
That means that the full area is given by the formula $$ S = \int_0^{2\pi} d\theta \int_0^H dz \,R$$ This is a simple case and you can choose whether to integrate over $z$ or over $\theta$. If you first integrate over $\theta$, you get to the case 1A, if you first integrate over $z$ you get to the case 1B.
Example 2: The full cylinder is parametrized by $\vec r(\rho,\theta,z) = (\rho\cos\theta, \rho\sin\theta, z)$ with $\rho\in[0,R]$,$\theta\in[0,2\pi]$,$z\in[0,H]$. The volume of infinitesimal element from the formulas above can be calculated to be $ dV = \rho d\rho d\theta dz$. The full volume is therefore $$ V=\int_0^R d\rho \int_0^{2\pi} d\theta \int_0^H dz\, \rho$$ If you integrate over $\theta$ and $z$ first, you get case 2A. If you integrate over $\rho$ and $z$ first, you get the correct frormula for case 2B: $$ V=\int_0^{2\pi} d\theta \frac12 HR^2$$
where factor $\frac12$ comes from the integral over $\rho$.
Example 3: The sphere can be parametrized by spherical coordinates $\theta$ and $\phi$, $\vec r(\theta,\varphi) = (R\sin\theta\cos\varphi, R\sin\theta\sin\varphi,R\cos\theta)$ with $\theta\in[0,\pi]$,$\varphi\in[0,2\pi]$. The surface element turns ut to be $dS = R^2\sin\theta d\theta d\varphi$. That gives the full area to be $$ S = \int_0^{\pi}d\theta \int_0^{2\pi} d\varphi R^2\sin\theta $$
If you integrate over $\varphi$ first, you get the integral that you've calculated. It is also possible to integrate over $\theta$ first obtaining formula $S = \int_0^{2\pi} d\varphi \,2R^2 $.
Example 4: The full ball has the parametrization $\vec r(r,\theta,\varphi) = (r\sin\theta\cos\varphi, r\sin\phi\theta\varphi,r\cos\theta)$ with $r\in[0,R]$,$\theta\in[0,\pi]$,$\varphi\in[0,2\pi]$. The volume element is $dV = r^2\sin\theta drd\theta d\varphi$, and the full volume is
$$ V = \int_0^R dr \int_0^{\pi}d\theta \int_0^{2\pi} d\varphi r^2\sin\theta $$
If you integrate over $\varphi$ and $\phi$ first (effectively calculating the area of the sphere from example 3), you get case 4A. You can also integrate over $r$ and $\varphi$ first, which would give you the formula $ V = \int_0^{\pi} d\theta \frac{2\pi}3 R^3 \sin\theta$.
To get case 4B, you need a different parametrization of the ball, one that allows to divide it into slices; that means that one of the coordinates must be $z$. The other two may be $\rho$ and $\varphi$, it won't matter in the end, as they will be the first ones to integrate out. The parametrization is $\vec{r} = (\rho\cos\varphi,\rho\sin\varphi,z)$ with $z \in[-R,R]$, $\phi\in[0,2\pi]$, $\rho\in[0,\sqrt{R^2-z^2}]$, $dV = \rho d\rho d\varphi dz$ and the full volume  being $$ V = \int_{-R}^R dz \int_0^{\sqrt{R^2-z^2}}d\rho \int_0^{2\pi} d\varphi \rho$$
After integrating over $\rho$ and $\varphi$ we get $$ V = \int_{-R}^R dz  \pi(R^2-z^2)$$
which after substituition $z=R\cos\theta$, $dz=-R\sin\theta d\theta$ gives the correct formula for case 4B $$ V = \int_0^\pi d\theta  \pi R^3 \sin^3\theta $$
A: Really, all of this is the subject of integration on manifolds, and for a reference, you should take a look at Spivak's Calculus on Manifolds. What you're doing is the simple, 2 and 3-dimensional case of what's really going on in higher dimensions (but the ideas you're using are really powerful, so don't dismiss them). In the general set up, one can ask: given a smooth $k$-dimensional compact manifold $M$ (or manifold with boundary) inside $\Bbb{R}^n$, what is the $k$-dimensional volume of $M$? So to give a general answer, one would have to discuss integration on manifolds.
If you can parametrize a portion of a manifold, then it is possible to give an explicit formula in terms of that parametrization. So, here's the general setup: Let $M$ be a compact $k$-dimensional smooth, oriented submanifold (with or without boundary) of $\Bbb{R}^n$ (so it's like a "nice" $k$-dimensional object sitting inside $\Bbb{R}^n$, with some regularity conditions). Let $W \subset \Bbb{R}^k$ be open, and let $\alpha:W \to M \subset \Bbb{R}^n$ be a $C^1$, orientation preserving injective map, such that for every $x \in W$, $\alpha'(x)$ is an $n \times k$ matrix with full rank. Also, let $d^kV$ be the ($k$-dimensional) volume element of $M$. Then, for any continuous function $f: M \to \Bbb{R}$, we have that (assuming all the integrals exist)
\begin{align}
\int_{\alpha(W)}f\, d^kV &= \int_W f \circ \alpha \cdot \sqrt{\det(g)} \tag{$*$} \\
& \equiv \int_W f (\alpha(x) )\cdot \sqrt{\det(g(x))} \, d^kx
\end{align}
where $\equiv$ means it's just a different notation for the same thing, $d^kx = dx_1 \cdots dx_k$, and $g$ is the $k \times k$ matrix whose $ij$ entry ($1 \leq i,j \leq k$) is given by the inner product
\begin{equation}
g_{ij} = \left\langle \dfrac{\partial \alpha}{\partial x_i}, \dfrac{\partial \alpha}{\partial x_j} \right\rangle
\end{equation}
In the integral equation above, the LHS is an integral over a subset of a manifold $M$, while the RHS is a typical multi-dimensional Riemann integral in $\Bbb{R}^k$, which we can reduce to $k$ one-dimensional integrals using Fubini's theorem.
You said you're a physics student, so a lot of this might seem foreign. But in practice, it is not that bad; it is really a mechanical process, once you know what the objects are. Right now, don't get caught up in the technicalities; understand how to use the formula. Afterwards, you can read up any book which treats integration on manifolds, to properly learn the subject.
1.)Volume of a Cylinder
Let's first consider the question of computing the volume of a cylinder of height $H$, and radius $R$. Call this cylinder $M$. In this case, we define $W = [0,R] \times [0, 2\pi] \times [0,H]$, and define $\alpha : W \to \Bbb{R}^3$ by
\begin{equation}
\alpha(r,\phi,z) =
\begin{pmatrix}
r \cos \phi \\
r \sin \phi \\
z
\end{pmatrix}
\end{equation} 
(i.e $\alpha$ is the cylindrical coordinate parametrization of the cylinder $M$ $\ddot{\smile}$). Note that $M = \alpha[W]$ (i.e the image of the set $W$ under the map $\alpha$). So, now if we want to find the volume of the cylinder $M$, we choose $f$ to be the constant function $1$, because the volume is often defined as $\int_M 1 \, dV$. Now, all that remains is to compute the matrix $g$, and the square root of its determinant. You can verify that
\begin{equation}
g(r,\phi,z) = 
\begin{pmatrix}
1 & 0 & 0 \\
0 & r^2 & 0 \\
0 & 0 & 1
\end{pmatrix}
\end{equation}
so, the square root of the determinant is $\sqrt{\det(g(r,\phi,z))} = r$. So, we have
\begin{align}
\text{vol}(M) &= \int_W (1 \circ \alpha)(r,\phi,z) \cdot r \, \, dr \, d\phi \, dz \tag{by ($*$)} \\
&= \int_W 1 \cdot r \, dr \, d \phi \, dz \\
&= \int_0^R \int_0^{2\pi} \int_0^H r \, dz \, d\phi \, dr \tag{by Fubini's theorem} \\
&= \pi R^2 H,
\end{align}
as expected. Of course, this all might seem like heavy machinery for the "simple" task of finding the volume of a cylinder, but it works in very general settings (if you have sufficient regularity so that all the integrals exist, etc... which usually in basic physics questions shouldn't be an issue).
In the above example, $M$ was the solid cylinder in $\Bbb{R}^3$, so it was a 3-dimensional object sitting inside a 3-dimensional space, and we computed its (3-dimensional) volume. i.e we had $k=n=3$
2.) Surface area of a Sphere 
Let $R>0$, and let $M = \{(x,y,z) \in \Bbb{R}^3: x^2 + y^2 + z^2 = R^2\}$ be the sphere centered at the origin, of radius $R$. (By the way, the surface area is what we might call $2$-dimensional volume, and it is customary to denote $d^2V$ by $dA$). So, in this case, we have $k=2$, while $n=3$. The power of the formula $(*)$, is that it works in every case, regardless of whether $k<n$ or $k=n$.
To compute the surface area, we follow the same procedure as in the previous example. Let $W = [0,\pi] \times [0,2\pi]$, and define $\alpha : W \to \Bbb{R}^3$ by
\begin{equation}
\alpha(\theta, \phi) = 
\begin{pmatrix}
R \sin \theta \cos \phi \\
R \sin \theta \sin \phi \\
R \cos \theta 
\end{pmatrix}
\end{equation}
Then $M = \alpha[W]$. In this example, you can compute for yourself that the matrix $g$ is given by
\begin{equation}
g(\theta, \phi) = 
\begin{pmatrix}
R^2 & 0 \\
0 & R^2 \sin^2\theta
\end{pmatrix}
\end{equation}
Hence, we have that $\sqrt{\det g(\theta,\phi)} = R^2 \sin \theta$. So, we have
\begin{align}
\text{surface area}(M) &= \int_{M} 1 \, dA \\
&= \int_{W} (1 \circ \alpha)(\theta,\phi) \cdot (R^2 \sin \theta) \, d \theta \, d\phi  \tag{by ($*$)} \\
&= \int_0^{2\pi} \int_0^{\pi} R^2 \sin \theta \, d \theta \, d\phi \tag{Fubini's theorem}\\
&= 4\pi R^2
\end{align}

These two examples were pretty easy, but in general the toughest thing will be to find the parametrizataion $\alpha: W \to M \subset \Bbb{R}^n$. After you figure that out properly, the rest is just turning the crank: calculate the partial derivatives of $\alpha$, take their inner products (dot products) to calculate the matrix $g$, take the determinant, and then use formula $(*)$. Lastly, use Fubini's theorem to reduce the $k$-dimensional integral $\int_W (\cdots)$ in $\Bbb{R}^k$ into $k$ one-dimensional integrals.
Now, the question which you might be wondering is how can one prove $(*)$? If you follow Spivak's book, then that equation can be proven by unravelling the definitions of the volume element, and how integration on manifolds was defined. That formula pretty much follows from Spivak's definitions, because all the key ideas have already been "encoded" inside the definitions, so all the reader is left to do is the extract that information.
But from a very rough heuristic point of view, it makes use of the multivariable change of variables theorem. Loosely speaking, if you have a small rectangle $R$ inside $W$ (i.e a rectangle in your "parameter space") and you fix a point $x_R \in R$, then when under the mapping $\alpha$, the volume of the rectangle $R$ gets distorted by an amount which is approximately proportional to the volume of the $k$-dimensional parallelepiped spanned by the $k$ vectors
\begin{equation}
\dfrac{\partial \alpha}{\partial x_1}(x_R) \quad, \dots,  \quad \dfrac{\partial \alpha}{\partial x_k}(x_R)
\end{equation}
this $k$-dimensional volume is equal to $\sqrt{\det (g)}$.

I'm sorry if this explanation for the formula $(*)$ wasn't very satisfying; this is really a subject which requires time to properly explain the definitions, and prove the relevant theorems etc.
