We consider the integral equation arising as a result of heat radiation exchange in both convex and nonconvex enclosures of diffuse grey surfaces. For nonconvex geometries, the visibility function must be taken into consideration. Therefore, a geometrical algorithm has been developed to provide an efficient detection of the shadow zones. For the numerical realization of the Fredholm integral equation, a boundary element method based on Galerkin-Bubnov discretization scheme is implemented. Consequently, multigrid iteration methods, which are closely related to two-grid methods, are used to solve the system of linear equations. To demonstrate the high efficiency of these iterations, we construct some numerical experiments for different enclosure geometries.