We study the boundary value problems for harmonic functions on open connected subsets of post-critically finite (p.c.f.) self-similar sets, on which the Laplacian is defined through a strongly recurrent self-similar local regular Dirichlet form. For a p.c.f. self-similar set $K$, we prove that for any open connected subset $Ω\subset K$ whose "geometric" boundary is a graph-directed self-similar set, there exists a finite number of matrices called $\textit{flux transfer matrices}$ whose products generate the hitting probability from a point in $Ω$ to the "resistance" boundary $\partial Ω$. The harmonic functions on $Ω$ can be expressed by integrating functions on $\partial Ω$ against the probability measures. Furthermore, we obtain a two-sided estimate of the energy of a harmonic function in terms of its values on $\partial Ω$.