I have recently learned about Rodin and Sullivan's work that proved a conjecture of Thurston involving giving a construction for the map in the Riemann mapping theorem using circle packings and this prompted me to wonder about the possibility of a similar sort of story for harmonic functions and the Dirichlet problem.
Let $D$ be the unit disk in the plane. Then given any continuous real-valued function $f$ on the boundary of $D$, there exists a unique harmonic extension $F$ of $f$ to $D$.
For $G$ a connected graph there is a similar condition, given any real-valued function $g$ on some nonempty subset of the vertices vertices of $G$, $S \subset V(G)$ , there is a unique function $G : V(G) \to \mathbb{R}$ such that $G$ extends $g$ and is harmonic on all of the vertices not in $S$.
Namely, let $G_\epsilon$ be the square grid in the plane, parallel to the coordinate axes with edge length $\epsilon$ and let $D_\epsilon$ be the intersection of this infinite planar graph with $D$. Using $f$, we can assign values to the boundary vertices of $D_\epsilon$ (say project radially to the boundary and take the value there) which we can then extend to give a harmonic function on $F_\epsilon : V(D_\epsilon) \to \mathbb{R}$.
Do these functions $F_\epsilon$ "converge" to $F$ as $\epsilon \to 0$?
The graphs $D_\epsilon$ are embedded in the plane, so that gives one way of making sense of convergence. A similar story could maybe apply also in higher dimensions, I'd be happy to hear about that as well.
(Maybe this is easier for the honeycomb tiling?)