In this paper, we present the multilevel Green’s function interpolation method (MLGFIM) for analyses of three-dimensional doubly periodic structures consisting of dielectric media and conducting objects. The volume integral equation (VIE) and surface integral equation (SIE) are adopted, respectively, for the inhomogeneous dielectric and conducting objects in a unit cell. Conformal basis functions defined on curvilinear hexahedron and quadrilateral elements are used to solve the volume/surface integral equation (VSIE). Periodic boundary conditions are introduced at the boundaries of the unit cell. Computation of the space-domain Green’s function is accelerated by means of Ewald’s transformation. A periodic octary-cube-tree scheme is developed to allow adaptation of the MLGFIM for analyses of doubly periodic structures. The proposed algorithm is first validated by comparison with published data in the open literature. More complex periodic structures, such as dielectric coated conducting shells, folded dielectric structures, photonic bandgap structures, and split ring resonators (SRRs), are then simulated to illustrate that the MLGFIM has a computational complexity of $O\left(N\right)$ when applied to periodic structures.