The Boundary Element Methods (BEM), relying on Galerkin approximation of boundary integral equations, have proven efficient in the numerical simulation of electromagnetic wave scattering phenomena in homogeneous media in the frequency domain. These methods, which only necessitate a surface mesh of the domain of interest, are rapid and accurate, but can show a lack of flexibility in the design phase of meshes. Indeed, they do not allow, e.g., hanging nodes that can appear either during a local refinement of mesh in order to capture the electromagnetic details or when joining pieces of dissimilar meshes together during optimization process. For complex 3D structures, such a mesh conformity requirement can thus limit the performance of these classical methods. To overcome these weaknesses and hence make classical BEM solvers more efficient, we develop a polytopal approximation of the integral equations by using the principle of the virtual element method. In particular, we will present the construction, the theoretical and numerical analysis as well as the implementation of such an approach for the well-known EFIE (Electric Field Integral Equation) used to model the scattering by a perfectly conducting object. Some numerical experiments based on standard electromagnetism benchmarks will prove the efficiency of this new scheme to predict the Radar Cross Section (RCS) of geometrical complex objects.