In this talk we consider the coupled Hamiltonian and momentum constraints in the Einstein equations, a nonlinear geometric elliptic PDE which must be solved numerically in general relativistic simulations of massive objects. The solution theory for the constraints on manifolds with boundary (the standard computational scenario) is not completely understood, so we begin by developing some results along these lines, including well-posedness and a priori estimates in various Sobolev and Besov norms. We then establish two quasi-optimal a priori error estimates for abstract Galerkin approximations, and derive a posteriori error estimates which lead to residual and duality error indicators for driving nonlinear approximation-type numerical methods. The a priori and a posteriori estimates establish a theoretical foundation for reliable and efficient use of adaptive finite element methods for the constraints and similar geometric PDE. In the latter part of the talk we describe a number of practical adaptive numerical techniques for solving general coupled nonlinear elliptic geometric PDE systems on 2- and 3-manifolds, and illustrate the techniques by considering some examples.