We present cubature methods approximating the surface integrals arising by Galerkin discretization of boundary integral equations on surfaces in 3. This numerical integrator does not depend on the explicit form of the kernel function, the trial and test space, or the surface parametrization. Thus, it is possible to generate the system matrix for a broad class of integral equations just by replacing the sub routine for evaluating the kernel function. We will present formulae to determine the minimal order of the cubature methods for a required accuracy. Emphasis is laid on numerical experiments confirming the theoretical results.