Discretization error of the finite element method (FEM) is an oft-neglected form of model misspecification, which can severely affect model prediction. Particularly in contexts where partial differential equations (PDEs) need to be solved repeatedly (e.g., inverse problems, topology optimization, multiscale modeling), a painful trade-off arises between evaluation cost and model accuracy. This trade-off is central to various fields of research, like model order reduction, adaptive mesh refinement and surrogate modeling, all of which aim to minimize computational cost while keeping discretization error sufficienlty close to zero. Recently, the field of probabilistic numerics has emerged, which takes a different viewpoint to the handling of numerical error. Instead of minimizing the error, it can be expedient to accept some error, provided that the associated epistemic uncertainty is propagated consistently. In this work, we present the Bayesian finite element method (BFEM), which provides a probabilistic model for the discretization error of FEM. The method assumes a Gaussian process prior over the solution space, whose associated reproducing kernel inner product is the energy inner product. By conditioning exactly on the residuals that follow from the test functions, a posterior distribution is obtained which recovers many properties from standard FEM: the posterior mean of BFEM equals the standard FEM solution, the posterior covariance obeys Galerkin orthogonality, and the energy minimization principle in FEM is equivalent to a loss function minimization interpretation of BFEM. This Bayesian interpretation of FEM is not only theoretically sound, but also practically useful, as it allows for the propagation of model-form uncertainty due to discretization error. We demonstrate BFEM in an inverse setting, and find that it avoids the confidently wrong parameter estimates to which FEM is prone.