We present a thermodynamically consistent approach to describing the process of ferroelectric polarization in piezoceramics of PZT type. In our phenomenological model, quantities such as displacement, strain and stress, as well as electric potential, electric field density, dielectric displacement and polarization are used to define the material's local state. Internal strain and polarization are introduced to account for non-reversible processes. Based on energy, enthalpy and dissipation densities, the material's characteristic response to electric and mechanic stimuli is then specified. While energy or enthalpy are used to define the reversible, electromechanically coupled behavior often termed Voigt's theory of linear piezoelasticity, a dissipation function governs the evolution of internal strain and polarization. Flow rules are implied which are intrinsically electromechanically coupled; we discuss proper choices to reproduce not only electric polarization hysteresis, but also asymmetric ferroelectric straining and mechanical depolarization [1]. We provide various formulations based on different thermodynamic potentials, where stress and strain, electric field and dielectric displacement are understood either as primary, independent or dual, dependent fields. In a finite element context, these choices give rise to either classical or mixed methods. The tangential-displacment normal-normal-stress method for elasticity has been shown to be of advantage when dicretizing thin structures using anisotropic, flat elements [2]. Also, a reduced set of shape functions has been proposed for the approximation of source-free fields where the condition of vanishing divergence and continuous normal flux is met exactly [3]. These have been used for the dielectric displacement in [4]. Here, we explore the benefits of such mixed elements in the proposed thermodynamic framework. Computational results are presented to assess the capabilities of different choices.