While in the first part of the book (until chapter 4), one gets by without differential forms, in this part they become quite useful. In fact, Chapter 5 goes full on into differential form mayhem even on the first page :), but there it is totally unavoidable. I am currently writing up some appendixes for the book and one of them will be a very short (without proofs) review of differential forms and Stokes theorem. Though that's not done yet. Best is to look through the chapter on differential forms in a book like baby Rudin.
(1) Yes, it is the differential form version of the Stokes theorem. Really in this setting it is the Green's theorem (which is Stoke's theorem in two variables, or perhaps, classical Stokes in the $xy$-plane). As to why we take out a small disc, that is because the function that we would want to integrate over $U$ has a singularity in $U$, so Stokes does not apply. So we take out a small disc around the singularity. Then the theorem does apply since everything is nicely $C^1$ on the $U \setminus \Delta_r$.
It is not a bad exercise to use the classical Green's theorem formulation and see if you can make the proof work. That is, use Green's theorem and then see if you can end up with the second integral in (**) using the classical formulation. It is slightly more tedious to write down that way, but it is not that terrible.
(2) For differential one-forms, $dx \wedge dy = - dy \wedge dx$ (let's say $x$ and $y$ are two coordinates). This means that $dx \wedge dx = -dx \wedge dx = 0$. Same thing works for the complex forms $dz$ and $d\bar{z}$ instead of $x$ and $y$ (if $z=x+iy$, then $dz = dx + i \,dy$, and $d\bar{z} = dx - i\,dy$). Then the equality just follows from the definition of the $d$ operator on differential forms:
$d(a\, dz + b\, d\bar{z}) =
\frac{\partial a}{\partial z} dz \wedge dz +
\frac{\partial a}{\partial \bar{z}} d\bar{z} \wedge dz +
\frac{\partial b}{\partial z} dz \wedge d\bar{z} +
\frac{\partial b}{\partial \bar{z}} d\bar{z} \wedge d\bar{z}
=
\frac{\partial a}{\partial \bar{z}} d\bar{z} \wedge dz +
\frac{\partial b}{\partial z} dz \wedge d\bar{z}
=
\left(\frac{\partial b}{\partial z}-\frac{\partial a}{\partial \bar{z}} \right) d\bar{z} \wedge dz$
(3) In the first guy on the left in (**) that is just $d\zeta$ (no wedge there) that's evaluated over the boundary of the disc. This is just the normal path integral from one variable complex analysis. In fact, this is precisely the same argument from the standard proof of Cauchy formula that is usually given (you reduce to a disc centered at a point). The idea is to now actually compute the integral, so parametrize by saying that $\zeta = z+re^{it}$, then $\zeta-z = re^{it}$ and $d\zeta = rie^{it}dt$. The denominator gets canceled, the $i$ also gets canceled, and you get the integral in the middle. Now this is where you use continuity, the integral is just the average of the value of $f$ over a tiny circle around $z$, so as $r \to 0$, the limit must be the value of $f$ at $z$.
That is a consequence of the “maximum modulus principle,“ and you can prove it using the Cauchy Integral Formula (compare https://proofwiki.org/wiki/Maximum_Modulus_Principle). The “trick” is to apply the CIF not to an arbitrary point, but to the point where $|F|$ attains its maximal value:
As a continuous function, $|F|$ attains it maximum on the closed disk $\bar D(0,1)$ at some point $\alpha$. Now assume that
$$
|F(\alpha)| = \max \{ |f(z)| : z \in \bar D(0,1) \} > 1 \, .
$$
It is clear that $\alpha$ lies in the open disk $D(0, 1)$. If $r > 0$ is sufficiently small then $D(\alpha,r) \subset D(0, 1)$ and we can apply the Cauchy Integral Formula to $\alpha$ and the circle $|z - \alpha| = r$.
Conclude that $|F(z)| = |F(\alpha)|$ on the circle $|z - \alpha| = r$ and consequently, $|F(z)| = |F(\alpha)|$ for $|z-\alpha| \le r$, i.e. that $|F|$ is constant in a neighborhood of $\alpha$.
Then use (for example) the Cauchy-Riemann equations to show that $F$ is constant in a neighborhood of $\alpha$.
Finally it follows from the identity principle that $F$ is constant in $D(0, 1)$ and because of the continuity, constant on $\bar D(0, 1)$. This is a contradiction because $|F(\alpha)| > 1$ but $|F(z)| \le 1$ on the boundary.
Best Answer
Sketch: Assume first $f$ is holomorphic on a neighborhood of $\overline {\mathbb D}.$ Write $f(z) = \sum_{n=0}^{\infty}a_nz^n.$ Recall that
$$\frac{1}{(1-u)^2} =\sum_{n=0}^{\infty}(n+1)u^n.$$
Hence
$$\frac{1}{(1-\bar z{\zeta})^2} =\sum_{n=0}^{\infty}(n+1)(\bar z{\zeta})^n.$$
Now compute
$$\frac{1}{\pi}\int_{\mathbb D}\left (\sum_{n=0}^{\infty}a_nz^n \right )\left (\sum_{n=0}^{\infty}(n+1)(\bar z{\zeta})^n \right )\, dA(z)$$
using polar coordinates and orthogonality of the functions $z^n$ with respect to area measure on $\mathbb D.$ This works out nicely.