In the last years the unique optical properties of graphene have attracted the attention of the scientific community, and novel graphene-based photonic, plasmonic and optoelectronic devices have been proposed for a plethora of applications. In this framework, numerical and analytical modelling play a crucial role both for design and analysis purposes. As a matter of fact, well-established models for the complex bi-dimensional linear conductivity of graphene have been reported in the literature. Nevertheless the analysis of graphene-based optical devices remains quite challenging from the numerical point of view. Indeed, it was demonstrated that graphene layers can be accurately modeled in conventional full-wave solvers by treating them as volumetric media with known conductivity and proper atomic thickness (< 0.5 nm), but the required fine discretization of these ultra-thin layers results in a huge computation burden. In order to overcome this limitation, a different and more efficient approach is mandatory.