## 2 Model Problems of Stationary Piezoelectricity  
First, we present the strong (or local) formulation of the model problem of stationary piezoelectricity. Due to further application of the matrix notation, convenient for finite elements, this formulation is expressed in the matrix form. The mechanical equilibrium equations and geometrical (kinematic) relations read  
\[\left. \begin{array}{c}\mathrm{div}\,\underline{\sigma}+\mathbf{f}=\mathbf{0}\\ \underline{\varepsilon}=\frac{1}{2}[u\otimes\mathrm{grad}\,+(u\otimes\mathrm{ grad}\,)^{T}]\end{array} \right\},\quad\mathbf{x}\in V, \tag{1}\]  
where the matrix representations \(\underline{\sigma}\) and \(\underline{\varepsilon}\) of the symmetric, stress, and strain tensors, \(\bar{\sigma}\) and \(\bar{\varepsilon}\), of the second rank are employed. Both relations are valid for each point \(\mathbf{x}\) of the volume domain \(V\) of the piezoelectric. The vector \(\mathbf{f}\) represents given mass loading vector, while \(\mathbf{u}\) is the unknown vector of (mechanical) displacements. Within the electric field, the electric equilibrium equation, namely the Gauss law with the assumed zero volume charges, holds together with the electric field vector \(\mathbf{E}\) definition, i.e.,  
\[\left. \begin{array}{c}\mathrm{div}\,\mathbf{d}=0\\ \mathbf{E}=-\mathrm{grad}\,\phi\end{array} \right\},\quad\mathbf{x}\in V. \tag{2}\]  
Above, the vector of electric displacements \(\mathbf{d}\), resulting from electrical charging of the piezoelectric, and the unknown scalar electric potential \(\phi\) appear. In order to express the stress tensor by the strain tensor, and the electric displacement vector by the electric potential one, the constitutive relations have to be added to the systems (1) and (2). The general form of these relations is consistent with the work of Reference [57]. The matrix form of them corresponds to the isotropic elasticity, orthotropic dielectricity, and orthotropic piezoelectricity. Because of this, the general, constitutive tensors of the elasticity, dielectricity, and piezoelectricity, \(\tilde{D}\), \(\tilde{\gamma}\), \(\tilde{C}\), of the fourth, second, and third ranks, respectively, can be represented by the matrices \(\mathbf{D}\), \(\mathbf{\gamma}\), and \(\mathbf{C}\), of the sizes: \(6\times 6\), \(3\times 3\), and \(6\times 3\), respectively, i.e.,  
\[\left. \begin{array}{c}\mathbf{\sigma}=\mathbf{D}\mathbf{\varepsilon}-\mathbf{C}\mathbf{E}\\ \mathbf{d}=\mathbf{C}^{T}\mathbf{\varepsilon}+\mathbf{\gamma}\mathbf{E}\end{array} \right\},\quad \mathbf{x}\in V, \tag{3}\]  
where the coupling of the mechanical and electric fields results from the presence of the matrix \(\mathbf{C}\), while \(\mathbf{\sigma}\) and \(\mathbf{\varepsilon}\) stand for the six-component vectorial representations of the stress and strain tensors \(\bar{\sigma}\) and \(\bar{\varepsilon}\). The above set of the partial differential, governing equations of piezoelectricity (1)-(3) has to be completed with the boundary conditions. The Neumann and Dirichlet boundary conditions, concerning, respectively, unknown functions derivatives and unknown functions themselves, become stress and displacement boundary conditions in the mechanical case, namely:  
\[\underline{\varepsilon}\mathbf{n}=\mathbf{p},\quad\mathbf{x}\in P \tag{4}\]  
\[\mathbf{u}=\mathbf{w},\quad\mathbf{x}\in W, \tag{5}\]  
where \(\mathbf{p}\) and \(\mathbf{w}\) are the given stress and displacement vectors, respectively, while \(P\) and \(W\) denote the parts of the surface \(S\equiv\partial V\) of the piezoelectric body \(V\), with the given values of the mentioned vectors. The vector \(\mathbf{n}\) represents the unit vector, outward and normal to the boundary.