Achieving high power densities through initiative designs of piezoelectric harvester in various geometries is a key point in vibration energy harvesting. State-of-the-art analytical and finite element models (FEMs) ignore structural damping, inter-laminar continuity, shear stresses, and contact layer effect between the substrate and piezoelectric layers and in addition cannot predict the performance of many recently introduced piezoelectric harvester configurations, such as non-uniform, thick piezoelectric patches, and variable thickness beams. This paper presents a comprehensive finite element formulation to calculate power generation by piezoelectric harvesters in a broader range of design cases. The presented high-order shear FEM not only is suitable for thick composite-based harvesters but also accommodates the drawbacks of the previous methods. The coupled finite element approach is verified versus experimental and analytical results. The model developed in this work is employed to analyze a non-uniform energy harvester with an E-glass fiber composite substrate layer sandwiched between piezoelectric layers with variable thickness. The numerical results show that, the advance formulation is capable of analyzing various piezoelectric harvesters including various influential parameters such as contact layer and damping dissipation. The results, furthermore, indicate that variable piezoelectric-layer thickness and an optimum fiber direction in composite substrate lamina can enhance performance of the piezoelectric harvester.