## 2 Model Problems of Stationary Piezoelectricity
In the case of the electric field, the Neumann boundary conditions express equality of the internal charges (defined with the unknown electric displacements \(\mathbf{d}\)) and given external charges of the surface density \(c\), while the Dirichlet boundary conditions describe the given values \(\omega\) of the unknown field of electric potential \(\phi\), i.e.,  
\[\mathbf{d}^{T}\mathbf{n}=-c,\quad\mathbf{x}\in Q, \tag{6}\]  
\[\phi=\omega,\quad\mathbf{x}\in F. \tag{7}\]  
Above, \(Q\) and \(F\) represent the parts of the surface \(S\) of the piezoelectric domain \(V\), where the electric charge and electric potential are known, respectively. It should be mentioned that the general constitutive relations (3) allow any piezoelectric model based on coupling of the linear elastic and linear dielectric models. In our proposition, we apply the 3D-based hierarchical piezoelectric models proposed in Reference [15; 58] and elucidated in Reference [11; 12]. These models employ three-dimensional degrees of freedom within the electric and mechanical fields. The 3D-based elastic models may represent: three-dimensional elasticity, hierarchical and first-order shell models, as well as the solid-to-shell and shell-to-shell transition models. The 3D-based dielectric models include: three-dimensional dielectricity model and the hierarchical symmetric-thickness models. These two groups of models were introduced and thoroughly described in Reference [13; 14] and Reference [16], respectively. Note that the applied 3D-based approach allows analysis of complex piezoelectric structures within the above Formulations (1)-(3). In such structures, more than one piezoelectric model meet, including transition ones. Thanks to the reciprocity theorem [59] of linear piezoelectricity, the above strong Formulations (1)-(7) can be converted into the corresponding variational (or weak) formulation. For the model problems of linear stationary piezoelectricity considered in this paper, the weak formulation reads  
\[\left\{ \begin{array}{rcl}B(\mathbf{v},\mathbf{u})-C(\mathbf{v},\phi)&=&L(\mathbf{v})\\ -C(\mathbf{v},\mathbf{u})-b(\psi,\phi)&=&-l(\psi)\end{array} \right. \tag{8}\]  
for all admissible displacements \(\mathbf{v}\in\mathcal{V}\) and all admissible potentials \(\psi\in\Psi\), with the space definitions: \(\mathcal{V}=\{\mathbf{v}\in(H^{1}(V))^{3}:\mathbf{v}=\mathbf{0}\text{ on }W\}\) and \(\Psi=\{\psi\in H^{1}(V):\psi=0\text{ on }F\}\). The solution functions belong to the spaces: \(\mathbf{u}\in\mathbf{w}+\mathcal{V}\) and \(\phi\in\omega+\Psi\), where the given values of \(\mathbf{w}\) and \(\omega\), consistent with the Dirichlet boundary conditions (5) and (7), are called lifts (compare Reference [1]) within the displacements field and potential field, respectively. The bilinear and linear forms of the first line of the relation (8) denote the virtual strain energy of the elastic field characterized with the matrix \(\mathbf{D}\) of elasticities and the virtual work of the external body \(\mathbf{f}\) and surface \(\mathbf{p}\) loadings. The following definitions of these forms hold:  
\[B(\mathbf{v},\mathbf{u}) = \int_{V}\mathbf{\varepsilon}^{T}(\mathbf{v})\mathbf{D}\,\mathbf{\varepsilon}(\bm {u})\,dV \tag{9}\] \[L(\mathbf{v}) = \int_{V}\mathbf{v}^{T}\mathbf{f}\,dV+\int_{P}\mathbf{v}^{T}\mathbf{p}\,dS\,.\]  
In the second Equation (8), the bilinear form represents the virtual energy of the electric field characterized with the dielectricity constant matrix \(\delta\) under constant strain, while the linear form is equal to the virtual work of the surface charges of the density \(c\), namely  
\[b(\psi,\phi) = \int_{V}\mathbf{E}^{T}(\psi)\,\delta E(\phi)\,dV\] \[l(\psi) = -\int_{Q}\psi\,c\,dS\,. \tag{10}\]The mixed bilinear forms \(C\) stand for the virtual energies corresponding to the electro-mechanical (piezoelectric) coupling.