An isoparametric 3D electromechanical hexahedral finite element integrating a 3D phenomenological ferroelectric and ferroelastic constitutive law for domain switching effects is proposed. The model presents two internal variables which are the ferroelectric polarization (related to the electric field) and the ferroelastic strain (related to the mechanical stress). An implicit integration technique of the constitutive equations based on the return-mapping algorithm is developed. The mechanical strain tensor and the electric field vector are expressed in a curvilinear coordinate system in order to handle the transverse isotropy behavior of ferroelectric ceramics. The hexahedral finite element is implemented into the commercial finite element code Abaqus® via the subroutine user element. Some linear (piezoelectric) and non linear (ferroelectric and ferroelastic) benchmarks are considered as validation tests.