We study boundary value problems of the form -u = f on and Bu = g on the boundary j , with either Dirichlet or Neumann boundary conditions, where is a smooth bounded domain in R n and the data f, g are distributions. This problem has to be first properly reformulated and, for practical applications, it is of crucial importance to obtain the continuity of the solution u in terms of f and g. For f = 0, taking advantage of the fact that u is harmonic on , we provide four formulations of this boundary value problem (one using nontangential limits of harmonic functions, one using Green functions, one using the Dirichlet-to-Neumann map, and a variational one); we show that these four formulations are equivalent. We provide a similar analysis for f = 0 and discuss the roles of f and g, which turn to be somewhat interchangeable in the low regularity case. The weak formulation is more convenient for numerical approximation, whereas the nontangential limits definition is closer to the intuition and easier to check in concrete situations. We extend the weak formulation to polygonal domains using weighted Sobolev spaces. We also point out some new phenomena for the "concentrated loads" at the vertices in the polygonal case.