Solution (source code)

= Solution

For an open subdomain $U\subset D$, identify $H_0^1(U)$ with the closed subspace of $H_0^1(D)$ obtained by zero extension. Its orthogonal complement is
$$
\mathcal H_U
=\{f\in H_0^1(D):(f,g)_\nabla=0
\text{ for every }g\in H_0^1(U)\}.
$$
If $f\in\mathcal H_U$, then testing against $C_0^\infty(U)$ and integrating by parts gives $\Delta f=0$ in $U$ in the <distributional derivative> sense. The <Weyl lemma> therefore gives a representative harmonic on $U$.

The <orthogonal decomposition by a closed subspace> gives
$$
H_0^1(D)=H_0^1(U)\mathbin\oplus\mathcal H_U.
$$
Project the isonormal process defining $h$ onto these two orthogonal subspaces. The projections are jointly Gaussian and uncorrelated, hence <independent random variables>[independent]. The first projection is a zero-boundary Gaussian free field $h_U$ on $U$; the second is a random distribution $h^{\mathrm{har}}$ that is harmonic on $U$. Thus
$$
h=h_U+h^{\mathrm{har}},
$$
with independent summands. This proves the <Domain Markov property of the Gaussian free field>.