A novel efficient technique is presented for the evaluation of domain integrals that appear in the boundary element method (BEM). Herein, the source term is approximated with the use of radial basis functions, as in the dual reciprocity BEM. The proposed technique, called the Kriging Integration Method (KIM), comprises the use of the Simple Kriging Method in non-overlapping patches for obtaining the weights of the integration points located inside. As it is necessary to compute the integrals of the covariance function prior to obtaining these weights, this can be efficiently realized using the Cartesian Transformation Method. The domain integrals over all the generated partitions are then computed and added to obtain the value of the whole-domain integral. Using KIM, it is possible to evaluate approximately weakly singular domain integrals over simply or multiply connected domains without applying any transformation or regularization method to the singular integrand. The numerical results obtained in several 2D potential problems demonstrate that this integration scheme is as accurate as both the dual reciprocity method and RIM and less time consuming than the RIM.