This is a re-interpretation of Ray Yang's answer, which also shows how the result can be generalized to other polygons. Introduce the function
$$v(x,y)=(1-\max(|x|,|y|))^+ ,\qquad (x,y)\in\mathbb R^2$$
This is a compactly supported Lipschitz function with support $Q=[-1,1]\times [-1,1]$. Its graph is a pyramid with $Q$ as the base.
If $u$ is harmonic in a neighborhood of $Q$, then integration by parts yields
$$0=\int_{\mathbb R^2} v\,\Delta u = \int_{\mathbb R^2} u\,\Delta v \tag{1}$$
By considering $u(\alpha x,\alpha y)$ with $\alpha\to 1^-$, we extend (1) to functions continuous in $Q$ and harmonic in its interior.
It remains to observe that $\Delta v$ is the distribution composed of
- the linear measure on $\partial Q$
- $-\sqrt{2}$ times the linear measure on the diagonals of $Q$
This follows from considering the discontinuities of the normal derivative of $v$ across the aforementioned lines; elsewhere $v$ is harmonic. One can also save the trouble of calculating the factor of $-\sqrt{2}$ by using the fact that $\int_{\mathbb R^2}\Delta v=0$.
It should be clear that there is nothing special about the square and its diagonals: any piecewise linear compactly supported function gives rise to a similar identity. For example, one can build a pyramid on top of any regular polygon $P$ and conclude that the average of a harmonic function along $\partial P$ is equal to its average over the union of segments connecting the vertices of $P$ to its center. No computations of slopes are necessary: it's clear that $(\Delta v)^+$ and $(\Delta v)^-$ are constant multiples of linear measure, and the identity $\int_{\mathbb R^2}\Delta v=0$ gives all the information we need.
Okay, I don't see that Ahlfors intended any particular method to be used (doesn't mean he didn't), so let's just use any method we can find, dominated convergence is the first that comes to mind.
First, for convenience, let's rotate it and consider
$$\int_{-\pi}^\pi \log\,\lvert 1 - e^{i\vartheta}\rvert \, d\vartheta,$$
so that we have the problematic point at $\vartheta = 0$.
A logarithmic singularity is integrable, since
$$ \int_\varepsilon^1 \log t\, dt = \left[t\log t -t\right]_\varepsilon^1 = -1 - \varepsilon\log\varepsilon + \varepsilon \to -1$$
for $\varepsilon \to 0$.
On the ray $r\cdot e^{i\vartheta}$ (for $\lvert \vartheta\rvert < \pi/2$), we consider
$$\delta(r) = \lvert 1- re^{i\vartheta}\rvert^2 = (1 - r\cos\vartheta)^2 + (r\sin\vartheta)^2 = 1+r^2 - 2r\cos\vartheta.$$
Since $\delta'(r) = 2(r - \cos\vartheta)$, $\delta(r)$ is minimised for $r = \cos\vartheta$, and $\delta(\cos\vartheta) = 1 - \cos^2\vartheta = \sin^2\vartheta$. That means that for any $0 < \varepsilon < \frac{\pi}{3}$ we have $\lvert \sin \vartheta\rvert \leqslant \lvert 1 - re^{i\vartheta}\rvert \leqslant 1$ for all $r\in [0,1]$ and $\lvert\vartheta\rvert \leqslant \varepsilon$, so the integrands $f_r(\vartheta) = \log\,\lvert 1 - re^{i\vartheta}\rvert$ are dominated by the integrable function $\log \frac{1}{\lvert\sin\vartheta\rvert}$ on $[-\varepsilon,\,\varepsilon]$. Outside that interval everything is bounded and convergence is uniform, so
$$\int_{-\pi}^\pi \log \,\lvert 1 - e^{i\vartheta}\rvert\,d\vartheta = \lim_{r \to 1} \int_{-\pi}^\pi \log\, \lvert 1 - re^{i\vartheta}\rvert\,d\vartheta = 0.$$
Another method that can be useful when you have a non-integrable limit, but such that the principal value of the integral exists (when you have a simple pole on the integration path, for example), is the removal of a small disk around the problematic point.
Now it would be good if I could draw a picture here, but since I can't, a description must do. Consider
$$B_\varepsilon := \mathbb{D}\setminus \overline{D_\varepsilon(1)} = \{z \in \mathbb{C} \colon \lvert z\rvert < 1, \lvert z-1\rvert > \varepsilon\}.$$
In a (simply connected) neighbourhood of $\overline{B}_\varepsilon$, the integrand is the real part of a holomorphic function ($\log (1-z)$), hence the (multiple of the) Cauchy integral
$$\int_{\partial B_\varepsilon} \frac{\log (1-z)}{z - z_0}\,dz$$
is $2\pi i \log (1-z_0)$ and vanishes for $z_0 = 0$.
The real part of that integral is then
$$\int_{\delta}^{2\pi-\delta} \log\, \lvert 1 - e^{i\vartheta}\rvert\, d\vartheta - \operatorname{Re} \int_{\partial B_\varepsilon \cap \mathbb{D}} \frac{\log (1-z)}{z}\,dz,$$
where $\delta = 2\arcsin(\varepsilon/2)$. The first integral is in the limit what we want, so we need to show that the second integral converges to $0$ for $\varepsilon \to 0$. The standard estimate does that,
$$\left\lvert \int_{\partial B_\varepsilon \cap \mathbb{D}} \frac{\log (1-z)}{z}\,dz\right\rvert \leqslant \pi\varepsilon \cdot \frac{\lvert\log \varepsilon\rvert + \pi/2}{1-\varepsilon}.$$
Best Answer
Basically if a function is harmonic, then its value at a point is equal to the integral over a sphere of any radius (or its boundary) centred at that point.
I.e. $u(x_0)=\frac{1}{n\alpha(n)r^{n-1}}\int_{\partial B(x_0,r)}u(x)dS(x)$
Now in your case, we consider the disk of radius $2$, so $r=2$, and in polar coordinates, you have $u(\theta)=3\sin(2\theta)+1$ So we have:
$u(0,0)=\frac{1}{4\pi}\int_0^{2\pi}3\sin(2\theta)+1d\theta=\frac{1}{4\pi}(\theta-\frac{3}{2}\cos(2\theta)|^{2\pi}_0)=\frac{2\pi}{4\pi}=\frac{1}{2}$