Multiscale piezoelectric materials and structures have shown their promising potential of applications in sensor, energy harvester, actuator and so on. However, the large complexity of such materials and structures poses a significant challenge to the accurate and efficient prediction of their piezoelectric behavior which is significant for guiding and optimizing their applications. In this paper, a Direct FE2 method was proposed for concurrent multilevel modeling of the coupled electromechanical behavior of multiscale piezoelectric materials and structures. The energy equilibrium and kinematic constraints between macro- and meso-scales were satisfied by prescribing periodic boundary conditions for both displacement and electric potential, which were derived from the governing equation of FEM. Another boundary condition was prescribed to the center node of each RVE, to prevent rigid motion and shifting of electric potential of meso-scale RVEs. All these boundary conditions can be easily prescribed to the Direct FE2 using multiple points constraints (MPCs), a commonly used feature in most commercial software. The accuracy, computational efficiency and ease of numerical implementation of the proposed Direct FE2 method were validated using a series of multiscale simulations, including the quasi-static response of homogeneous piezoelectric panels, an arc honeycomb structure, a regular honeycomb structure and a piezoelectric composite panel, as well as the dynamic response of the piezoelectric composite panel. It shows that the proposed Direct FE2 method will be valuable for large-scale simulations of complicated multiscale piezoelectric materials and structures.