3D-Based Transition hpq/hp-Adaptive Finite Elements for Analysis of Piezoelectrics

This paper concerns the algorithm of transition piezoelectric elements for adaptive analysis of electro-mechanical systems. In addition, effectivity of the proposed elements in such an analysis is presented. The elements under consideration are assigned for joining basic elements which correspond to the mechanical models of either the first or higher order, while the electric model is of arbitrary order. In this work, three variants of the transition models are applied. The first one assures continuity of displacements between the basic models and continuity of electric potential between these models, as well. The second transition piezoelectric model guarantees additional continuity of the stress field between the basic models. The third transition model additionally enables continuous change of the strain state between the basic models. Based on the mentioned models, three types of the corresponding transition finite elements are introduced. The applied finite element approximations are hpq/hp-adaptive ones, which allows element-wise changes of the element size parameter h, and the element longitudinal and transverse orders of approximation, respectively, p and q, depending on the error level. Numerical effectiveness of the models and their approximations is investigated in the contexts of: ability to remove high stress gradients between the basic and transition models, and convergence of the numerical solutions for the model problems of piezoelectrics with and without the proposed transition elements.


Introduction
More and more common application of piezoelectrics in contemporary technology requires more and more efficient and accurate methods of their modeling and analysis. In this work, we develop the new transition elements which possess unique set of features enabling joining mechanical shell models of the first order and higher orders. Such models are incompatible due to the assumptions of the plane stresses and no elongation of the normals to the shell mid-surface, both present in the first-order model. In order to join such models, the transition models are necessary. This refers also to the piezoelectricity case, when the mechanical field modeling needs simultaneous application of the first-and higher-order mechanical models.
Adaptive capabilities of hpq/hp-type are the second important feature of the new elements, where h is the element size parameter, while p and q are the element longitudinal and transverse orders of approximation. As our implementation of these capabilities is based on the standard hpq-approximation rules, presented, respectively, in Reference [1,2] for 2D and 3D, we will not elucidate this aspect in this paper.

State of the Art
To the best of the authors' knowledge, the transition piezoelectric elements have not been proposed yet, both in the classical (non-adaptive) and adaptive versions. Because of

Hierarchical Models and Elements for Monolithic Piezoelectrics
The hierarchical two-dimensional conventional (mid-surface unknowns) models of piezoelectric plates were introduced and mathematically substantiated in Reference [9]. The mixed formulation is applied there with stresses and electric displacements completing displacements and electric potential. Carrera and co-workers [10] proposed hierarchic plate models up to the fourth-order. The formulation is displacement-based in the case of the higher-order models, and mixed in the case of the Reissner-Mindlin model. The models employ three-dimensional degrees of freedom (dofs) for the top and bottom surfaces, and conventional (mid-surface) higher-order dofs. In addition, a thermal field can be included into this formulation. In Reference [11], the 3D-based (three-dimensional unknowns only) piezoelectric hierarchy of solid, first-order (Reissner-Mindlin) and hierarchical shell, and solid-to-shell and shell-to-shell transition piezoelectric models are formulated. Their hpq-adaptive finite element approximations are presented in Reference [12]. The formulation is displacement-like with displacements and electric potential as unknowns of the problem. In the case of the electric field, this approach allows for the three-dimensional and 3D-based symmetric-thickness models of dielectricity. The approach is based on the 3Dbased hierarchical models [13] and approximations [14] for elasticity, and the corresponding models [15] and approximations [16] for dielectricity.

Models for Piezoelectric Composites
In the case of laminated piezoelectric structures, it is worth mentioning the following works. In Reference [17], the authors presented Reissner-Mindlin-type composite plate model with piezoelectric layers. The direct piezoelectric and pyroelectric effects are taken into account. In turn, the work of Reference [18] proposed the plain strain shell model with higher-order shear and normal deformation effects for analysis of cylindrical laminated piezoelectric shells with actuating and sensing piezoelectric layers. The work of Reference [6], mentioned in the context of the monolithic piezoelectrics, should be mentioned here as applied to tree-layered shell structures with two outer piezoelectric layers and an electrically neutral core. A similar sandwich plate structures are modeled as two Reissner-Mindlin piezoelectric layers and one neutral layer of the same character [19], with the mid-surface interdependent displacement dofs. The electric fields are of mixed character and are defined by the outer layers' mid-surface electric potentials and independent electric field mid-surface values of the outer layers. In order to remove locking, the reduced integration method is applied. The generalized analytical Reissner-Mindlin model of laminated piezoelectrics was presented in Reference [21]. It allows analyses of the problems where the electric potential is not prescribed through the thickness. The method assumes small thickness and thus allows decomposition of the three-dimensional Reissner-Mindlin description into two-dimensional in-plane and one-dimensional through-the-thickness problems. Subsequently, the already cited work of Reference [10] concerning hierarchical models of single layered piezoelectrics can also be applied to the numerical analysis of hierarchically modeled multilayered structures, with piezoelectric layers included. The recent example [22] of analysis of multilayered shells concerns the first-order model of the functionally graded layers. In turn, the piece-wise second-order finite element model for multilayered structures with delamination and debonding was proposed in Reference [23].

Transition Piezoelectric, Dielectric and Elastic Models And Elements Transition Elements for Three Physical Classes of Problems
It can be concluded from the previous subsection that the piezoelectric transition models and elements that serve joining the three-dimensional (or hierarchical shell) piezoelectric models with the first-order piezoelectric model are rare. The only works, which we have cited, concern the idea [15] and variational [11] and finite element [12] formulations of the same classical transition piezoelectric model and element. They guarantee continuity of the displacement and transverse deformation fields between the mentioned models. They also guarantee continuity of the electric potential between the piezoelectric models.
In Reference [16], it has been demonstrated that, in the case of dielectric field of electric potential used in electrostatics, no transition elements between the three-dimensional model and symmetric-thickness hierarchical models of the first and higher orders are necessary if the consistent finite element approximations of non-adaptive or adaptive character are applied. This is because the same definition of electric field (intensity) is applied within both dielectric models.
In the case of the elastic mechanical fields, two groups of the transition models can be distinguished. The first one guarantees continuity of the inconsistent generalized displacement fields of the neighboring models. In such inconsistent models, the kinematic unknowns are different or of different order, e.g., in the transverse direction. In the second group, additionally the definitions of the derivatives of the kinematic unknowns, e.g., strains and/or stresses, are different in the neighboring models and the corresponding finite elements. In this context, one may deal with the neighborhood of the threedimensional elasticity and shell or plate theories, the three-dimensional theory, and beam or truss theories, as well as the shell/plate theories and the one-dimensional theories of beams or trusses. Due to the paper scope, only the first case will be elucidated here. A general review, including all three cases, can be found in the thesis [24].

Non-Adaptive Transition Elements for Elasticity
We start with the transition elements of the classical (non-adaptive) finite element methods where the neighboring elements are joined through the sides of the same size and finite element interpolations on these sides are consistent. In addition, location of three-dimensional, shell/plate, and transition domains is fixed and strictly determined by the structure geometry. The first examples of the transition elements of this type can be attributed to Surana [25], who took into account the consideration of Ahmad [26] on thick shell elements and his own research [27] on general and axi-symmetric shell elements. The work, in Reference [25], of Surana presents solid-to-shell elements which enable joining the thick-shell element with the linear, quadratic or cubic solid elements. The elements possess the bigger shell part equipped with the mid-surface dofs, and the smaller solid part limited to one side of the element and equipped with the tree-dimensional dofs. Plane stress is assumed in the entire transition element. Applications of the elements of this type to axi-symmetric [28] and three-dimensional [29] problems are presented. Extension onto thermo-mechanical analysis was done in Reference [30]. Note that the independent research of Bathe and Bolourchi [31], based on the mentioned work of Ahmad, led to the analogous transition elements. Gmür and Schorderet [32] suggested application of the three-dimensional or plane stress states in the Surana's element depending on the geometrical shape of the domain the element is located in. Liao and others [33] extended application of the Surana's elements onto laminated composite structures. Both, the shell and solid parts of the transition element are modeled as such composites. In the work of Reference [34], in order to avoid the numerical locking, the method based on strain interpolation (with the points where the strain values are known) was applied instead of the reduced integration.
In the work of Dávila [35], the opposite idea is applied, i.e., the elements consist of the larger solid part and smaller shell part. The latter part comprises one or two lateral sides. The top and bottom three-dimensional dofs and mid-surface dofs are applied within the solid and shell parts, respectively. Three-dimensional stress state is assumed throughout single-layered shells or each layer of composite shells. Yet another idea was suggested in Reference [36], by Gong, who proposed generalization of the solid-to-shell transition elements for the elements joining two-dimensional and three-dimensional elements. In this approach, the shell or two-dimensional parts consist of one or two lateral sides of the transition element. The longitudinal order of approximation can be higher than three. In the entire element, either three-dimensional or plane stress state can be applied.

Adaptive Transition Elements for Elasticity
In the case of the adaptive methods, finite elements of different sizes, orders of approximation (non-conforming approximations), and mechanical models are present in finite element meshes. In the case of different element sizes resulting from h-adaptivity based on remeshing, continuity of the field of unknowns may be based on introduction of the distorted elements on the boundaries between the mesh parts of various density. In the case of different sizes resulting from h-adaptivity based on local element refinement, continuity of the fields of unknowns requires introduction of the distorted elements [37] or is based on transition elements for monolithic [38] or multilayered [39] structures. Such transition elements deliver piece-wise linear [40] or approximate [41] continuity along the entire common edge or face, or at the common nodes only, e.g., within triangles [42], quadrilaterals [43] and hexahedrons [44]. Note that application of the smart idea of the constrained approximation, proposed by Demkowicz for 2D [1] and 3D [2], removes the necessity of introduction of the transition elements of this type as the continuity is automatically guaranteed by the constraints on the boundary between the elements of different sizes.
As far as the different approximation orders of the elements are concerned, transition elements may be introduced with the consistent approximation on the boundary between specific 2D [20] or general 3D [45] elements, with the continuity guaranteed again on the entire common edge or face, or at the common nodes only. A smart alternative is the idea of 2D [1] and 3D [2] hierarchical shape functions, defined independently in the element vertices, on its edges and sides, and in the interior of the element, which allows for the same approximation order on the common edge or face of the neighboring elements and different orders in their interiors.
In the case of different models of the neighboring elements, the continuity of the field of unknowns may require transition elements which conform the model of the larger number of nodal dofs with some of these dofs constrained to zero on the boundary with the model of the smaller number of nodal dofs. The 3D-2D [46], 2D-2D [47], and 2D-1D [48] versions of this approach can be found. The alternative is the transition elements equipped with the degrees of freedom of both models (see comments in Reference [49]). Note that such transition elements are not necessary if the uni-dimensional, e.g., three-dimensional, approach is applied within each of the neighboring models. The idea lies in using the same, e.g., three-dimensional, dofs for each model. This usually requires internal constraints in the model simplified (e.g., 3D-based shell models) with respect to the three-dimensional one. The linear [50] or higher [51] order of approximation can be applied in the transverse direction on the common face of the three-dimensional and 3D-based shell elements. Even though this technique guarantees continuity of the field of unknowns, high gradients of the derivatives on the boundary between the models appear. If one wants to get rid of such gradients, the transition elements removing stress [52] and, additionally, strain, Reference [24] gradients are necessary between the basic models being joined.

Conclusions
In this paper, we will use piezoelectric transition elements combining the ideas present in the following works concerning the uncoupled problem of elasticity. Firstly, we will apply the 3D-based approach [13], where the solid [53], hierarchical shell [14], first-order shell [54], and transition [51] models are equipped with the three-dimensional dofs only. This approach will also be extended onto the dielectricity models [16]. We will also apply the idea of Dávila [35], where one deals with larger 3D part and smaller thin-walled part of the transition element. The stress state within the element will be either three-dimensional, as in Reference [35], or transition (from three-dimensional to plane one) [24]. The hpadaptive capabilities of the elements will be based solely on the ideas of Demkowicz and others [2]. Our 3D-based hierarchical finite element models will resemble the ideas presented, respectively, in Reference [55,56] for the cases of the three-dimensional and conventional (mid-surface) dofs.

The Scope and Novelty of The Paper
This paper concerns the algorithm of the solid-to-shell and shell-to-shell transition piezoelectric elements for adaptive analysis of electro-mechanical systems. Such systems can be composed of piezoelectric transducers and elastic structural elements, both of arbitrary shapes. In addition, effectivity of application of the proposed elements in such an analysis is of our interest in the paper. The elements are assigned for joining basic elements, namely piezoelectric shell elements of the first-order mechanical field and arbitrary-order electric field and piezoelectric elements of the hierarchical shell model (second-or higherorder) within the mechanical field and arbitrary-order within the electric field. Alternatively, the latter model can be replaced with the piezoelectric description based on hierarchical models of three-dimensional elasticity and three-dimensional dielectricity.
In this work, three variants of the transition models are applied. The first one, called classical, assures continuity of displacements between the basic models and continuity of electric potential between these models, as well. The second, modified transition piezoelectric model guarantees additional continuity of the stress field between the basic models. The third, enhanced model additionally enables continuous change of the strain state between the basic models. All three models are 3D-based ones, namely only threedimensional degrees of freedom are applied within the mechanical and electric fields. Based on the mentioned models, three types of the corresponding transition finite elements are introduced. The applied approximations are hpq/hp-adaptive ones, with h being the element size parameter, p standing for the longitudinal order of approximation, and q denoting the transverse order of approximation.
Due to basic character of the presented research, numerical effectiveness of these models and their finite element approximations is limited to two crucial aspects, namely ability to remove high stress gradients between the basic and transition models, and convergence of the numerical solutions for the model problems of piezoelectrics. These solutions are obtained with and without the proposed transition elements. Error estimation and adaptivity control matters are out of scope of this work and will be presented in a separate paper.
Novelty of the paper results from the above unique features of the proposed transition elements and the orignal results concerning their effectiveness.

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 div σ + f = 0 where the matrix representations σ and ε of the symmetric, stress, and strain tensors,σ and ε, of the second rank are employed. Both relations are valid for each point x of the volume domain V of the piezoelectric. The vector f represents given mass loading vector, while 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 E definition, i.e., div d = 0 Above, the vector of electric displacements d, resulting from electrical charging of the piezoelectric, and the unknown scalar electric potential φ 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,D,γ,C, of the fourth, second, and third ranks, respectively, can be represented by the matrices D, γ, and C, of the sizes: 6 × 6, 3 × 3, and 6 × 3, respectively, i.e., where the coupling of the mechanical and electric fields results from the presence of the matrix C, while σ and ε stand for the six-component vectorial representations of the stress and strain tensorsσ andε. 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: σn = p, x ∈ P (4) where p and w are the given stress and displacement vectors, respectively, while P and W denote the parts of the surface S ≡ ∂V of the piezoelectric body V, with the given values of the mentioned vectors. The vector n represents the unit vector, outward and normal to the boundary.
In the case of the electric field, the Neumann boundary conditions express equality of the internal charges (defined with the unknown electric displacements d) and given external charges of the surface density c, while the Dirichlet boundary conditions describe the given values ω of the unknown field of electric potential φ, i.e., 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 symmetricthickness 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 for all admissible displacements v ∈ V and all admissible potentials ψ ∈ Ψ, with the space The solution functions belong to the spaces: u ∈ w + V and φ ∈ ω + Ψ, where the given values of w and ω, 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 D of elasticities and the virtual work of the external body f and surface p loadings. The following definitions of these forms hold: In the second Equation (8), the bilinear form represents the virtual energy of the electric field characterized with the dielectricity constant matrix δ under constant strain, while the linear form is equal to the virtual work of the surface charges of the density c, namely The mixed bilinear forms C stand for the virtual energies corresponding to the electromechanical (piezoelectric) coupling. Their definitions read: where, as before, C is the piezoelectric constant matrix under constant strain. The tensorial version of the above Formulations (8)- (11) can be found in Reference [11,57]. The existence and uniqueness theorem for the above linear piezoelectricity problem can be found in Reference [60]. The proof takes advantage of the generalized Lax-Milgram lemma [61,62].

Transition Models For Piezoelectricity
As said in the previous sections, the transition elements under consideration serve joining the first-order shell mechanical field of a piezoelectric with the higher-order hierarchical shell models or the three-dimensional elasticity model of this field within the piezoelectric. In the electric field corresponding to the mechanical fields of both types, the same hierarchical symmetric-thickness or three-dimensional dielectric model is applied. In order to introduce the original transition models, the basic piezoelectric models that are to be joined together will be presented first.

Basic Models of Linear Piezoelectricity
The starting point for two basic models of piezoelectricity, based on the first-order and higher-order mechanical fields, is the model corresponding to the three-dimensional piezoelectricity. The matrix form of the constitutive relations of linear three-dimensional piezoelectricity, expressed by the experimentally obtainable constitutive constants, is presented in Reference [57]: Above, D −1 , δ, and c denote the isotropic inverse matrix of elasticities (or compliance matrix), the orthotropic matrix of dielectric constants under constant stress, and the coupling orthotropic matrix of piezoelectric constants under constant stress, respectively. The six-component column vectors σ, ε represent the symmetric stress and strain tensors, σ ij and ε ij , i, j = 1, 2, 3, while d and E are the column three-component electric displacement and electric field vectors, with components d j , E j , j = 1, 2, 3. The mentioned tensors and vectors are defined locally, namely in the directions consistent with the axes of piezoelectric orthotropy. In the presentation below, the third local direction is the polarization direction of the piezoelectric. In the case of the thick-or thin-walled structures, this direction coincides with the transverse one.
The standard definition of the isotropic elasticity matrix reads: where D = E/[(1 + ν)(1 − 2ν)]. As it can be seen above, the non-zero terms of the matrix are expressed by Young's modulus E and Poisson's ratio ν.
The dielectricity constant matrix under constant stress reads: where δ ii , i = 1, 2, 3 stand for the orthotropic permittivity constants.
In the case of the piezoelectricity constant matrix under constant stress, one has where c 13 , c 23 , c 33 , c 52 , and c 61 are non-zero piezoelectric constants within this matrix.

Hierarchical or Three-Dimensional Models
For both the cases, when the transverse mechanical and electric fields are polynomially constrained through the thickness of the thin-walled piezoelectric member or the transverse direction is treated in the same way as the longitudinal directions of the three-dimensional piezoelectric member, the same constitutive relations hold.
In order to present these relations, let us convert the Equation (12) to the form suitable for the variational and finite element formulation of the general piezoelectric problem. This needs inversion of the first Equation (12) and its substitution into the second one. The resultant form reads: where the following definitions of the orthotropic matrix C of piezoelectric constants under constant strain and the isotropic matrix γ of dielectric constants under constant strain were introduced: The first of the above definitions gives the following non-zero terms of the matrix γ: The second definition leads to the following structure of the matrix C:

First-Order Piezoelectric Model Stress State
In the case when the mechanical field of the symmetric-thickness thick-or thin-walled member is linearly constrained through the thickness, the plane stress assumption of the first-order shell model have to be introduced into the three-dimensional piezoelectric constitutive relations. This means that the third row equation of the first constitutive relation (16), namely: has to be combined with the assumption σ 33 = 0, so as to give Now, the condition (20) can be rewritten in the following alternative identity-equation form Substitution of the relation (21) into the first relation (16), arithmetic manipulations on the terms containing E 3 (with use of isotropic definitions of the non-zero terms of D ij , i, j = 1, 2, . . . , 6 from (13)) lead to: The alternative form of the above equation can be obtained after arithmetic manipulations on the terms containing ε 11 and ε 22 . This leads to the change of the constitutive constants of the matrix D, with zero contributions of its third row and third column. This alternative form requires either retaining ε 33 as the third term of the strain vector or removal of: the third row and the third column in D, the third row in ε, and the third row in C.
Taking the plane stress assumption into account in the second constitutive relation (16) needs: introduction of (21) into this relation, arithmetic manipulations on the terms containing E 3 , and summation of the contributions from the terms containing ε 11 and ε 22 . The latter operation leads to the change of the piezoelectric constants in C T , including zero contributions of the third column of this matrix. All these lead to: Note that the third column of C T and the third row of ε can be removed from the above relation. The alternative form of (24) does not need any changes in C T from (16) but requires replacement of ε 33 with −ν 1−ν (ε 11 +ε 22 ), as in the Equation (23).
The chosen and written down forms (23) and (24) of the constitutive relations are suitable for further considerations in the next sections of the paper.

Kinematics
The second group of assumptions required by the first-order piezoelectric symmetricthickness model is of kinematic nature and requires deformation of the normals to the midsurface onto normals, not necessarily perpendicular to the mid-surface, and no elongation of these normals during deformation. The first assumption reads: where V represents the symmetric-thickness domain or such a part of the complex domain. Above, the local displacements u can be transformed into the global ones u typical for the applied 3D-based approach, with use of the transformation matrix θ T such that u j = θ ij u i . It is obvious that i = 1, 2, 3 represent the global directions, while j = 1, 2, 3 denote the local directions of which the first two are tangent and the third one is normal to the mid-surface S m of the symmetric-thickness piezoelectric structure. It can be noticed that the above condition is expressed by the top and bottom local displacements, u t and u b that can be described with u t = u for x ∈ S t and u b = u for x ∈ S b , where S t and S b are the top and bottom surfaces within the symmetric-thickness structure. Note also that the first and second terms of the right hand-side of the above condition represent the local displacements of the mid-surface and the displacements due to rotation of the normal. The latter displacements change linearly along the normal as t is the structure thickness, while x 3 ∈ (− t 2 , t 2 ) measures the distance from the mid-surface S m . The second assumption says that the third local displacements of the top and bottom surfaces are the same. Because of this no elongation is possible in this direction, namely: where V is defined as in (25). One can express the above condition by the global displacements with use of u 3 = θ i3 u i and u t

Transition Piezoelectric Models
Three piezoelectric transition models (classical, modified, and enhanced) will be presented. In each of these three models, continuity of the displacement and electric potential fields, on the boundaries between the transition and basic models, is maintained due to three-dimensional approach which lies in application of the three-dimensional displacements and three-dimensional potential for any model of the piezoelectric structure. The hierarchical or first-order models are consistently treated as constrained threedimensional models.
The classical transition model is characterized with the mechanical field of the threedimensional model. Only on the boundary with the first-order model, the plane stress assumption is valid. In the case of the modified transition model, the stress field changes gradually from of three-dimensional one (on the boundary with the three-dimensional domain) to the plane stress one (on the boundary with the first-order model). The same stress field description is valid also for the enhanced transition model. The enhanced transition model is characterized with the gradually changing assumption of no elongation of the lines perpendicular to the mid-surface of the domain. This assumption changes from its total validity on the boundary with the first-order model to its lack on the boundary with the three-dimensional model. Note that, in the cases of the classical and modified models, this assumption is valid only on the boundary with the first-order model. Apart from this boundary, this assumption does not apply.

The Classical Approach Stress State
The mechanical and electric fields within the classical transition piezoelectric model are defined as in the case of the three-dimensional piezoelctricity characterized with the constitutive Equation (16), completed with the definitions (13), (18), and (19).

Kinematic Assumptions
In this classical transition model, the kinematic assumption of the first-order mechanical theory consisting in deformation of the normals to the mid-surface into normals during deformation and the lack of elongation of these normals are valid on the boundary R between the transition and first-order theories, i.e., where V = V ∪ ∂V, ∂V ≡ S. As in the case of the first-order model, the above conditions can be expressed by the global displacements with use of the same relations as before. Note that the mid-surface within the transition domain may only be defined on the boundary R between the first-order and transition models.

the Modified Transition Model The Assumed Stresses
The second approach is based on the transition character of the plane stress assumption. This assumption is valid on the boundary with the first-order piezoelectric model. It does not hold, however, on the boundary with the hierarchical higher-order or threedimensional piezoelectric model, i.e., three-dimensional stress state is present on this boundary. Between these two boundaries, the stress state is intermediate, namely: (28) Above, the function α = α(S m ) ∈ 0, 1 , with S m being the mid-surface of the thick-or thin-walled part of the transition domain, is a blending function equal to 1 at the boundary with the three-dimensional model, and 0 at the boundary with the first-order model. Consequently, the definition (28) becomes identical with (20) for α = 1, and with (22) for α = 0.
Taking the above equation into account in the third row of the first relation (16) leads to: σ = σ 1 − σ 2 , where the first component equals: while the second one is equal to: In addition, the second equation (16) can be divided into two components, i.e., d = d 1 + d 2 , where the first component is: while the second component reads Note that, when α = 1, the relations (29)-(32) transform into Equations (16)- (19), while, for α = 0, into Equations (23) and (24).

The Assumed Kinematics
As it comes to the kinematic assumptions within the modified transition element, they are the same as in the case of the classical elements, i.e., the conditions (27) are valid on the boundary R between the transition and first-order models.

The Enhanced Transition Model Stress State
In the case of the enhanced transition model of piezoelectricity, constitutive relations are the same as in the modified model. This means that (29) and (30), as well as (31) and (32), hold.

Displacement assumptions
The starting point for defining the kinematic assumption within this type of the transition model is the definition of the displacement field of the 3D-based hierarchical shell models [13,14]. This field can be expressed with use of the local or global displacements, u or u, which are related by the transformation matrix θ T , i.e., u = θ T u. The local displacement components are equal to: where V is the symmetric-thickness domain. In the first line above, u l j , l = 1, 2, ..., I, are the functions describing changes of displacements in the longitudinal directions. In the case of the symmetric-thickness geometry, these function are the same for each l. The functions f l (x 3 ), l = 1, 2, ..., I stand for the polynomial functions in the direction normal to the midsurface. Along this direction the coordinate x 3 ∈ (− t 2 , t 2 ) is employed. The polynomial functions in this direction span the respective polynomial space of order I. Hence, the quantity I represents the hierarchical model order. It is worth mentioning that with I → ∞ one reaches three-dimensional elasticity description of the above displacement field.
Following our propositions from Reference [14,16], in the second and third lines of (33), we define the first and the last functions, f 0 (x 3 ) and f I (x 3 ), as linear ones and the rest functions f l (x 3 ), l = 1, 2, ..., I − 1 as polynomials of order I. After some arithmetic manipulations, one can reach the final form (line) of (33).
In the case of the enhanced transition piezoelectric model, its transverse displacement field is assumed to change from the state of no elongation of the lines perpendicular to the mid-surface to the state of lack of such an assumption. In the same manner, the order of the transverse displacement field is assumed to change from I = 1 to an arbitrary value I. So as to implement such changes, the blending function α, now playing the role of a gradually switching function, is applied, namely: where V is the symmetric-thickness domain or the domain with the mid-surface defined on R at least, while j = 1, 2, 3 and k = 1, 2. The later index k corresponds to the first two, longitudinal directions, while the index 3 is used for the transverse direction. It can be noticed that with α ≡ 1, (34) becomes (33) which characterizes the hierarchical or three-dimensional models, while, for α ≡ 0, (34) becomes (25) and (26) which correspond to the first-order model. In this case, l = 0 and l = 1 refer to the bottom (l ≡ b) and top (l ≡ t) surfaces.

Algorithms of the Transition Piezoelectric Elements
The starting point for our finite element considerations is the variational Formulations (8)- (11). Derivation of the corresponding finite element formulation needs discretization of the domain V and into finite element subdomains V e and introduction of polynomial interpolation of the displacements u and electric potential φ (see Reference [11,12]) with use of the global vectors of the nodal, displacement q hpq and potential ϕ hπρ , degrees of freedom (dofs). The applied interpolations are suited for the hpqand hπρ-adaptive approximations of the displacements and potential fields, with h being the common element size for both fields and independent longitudinal and transverse approximation orders p, q and π, ρ within the mechanical and electric fields, respectively. All these lead to where the components of the first equation correspond to: the first equation (9), the second Equation (9) and the first Equation (11)

Normalized Geometry of the Transition Elements
For three (classical, modified, enhanced) transition piezoelectric models two principal geometries of the prismatic elements has to be introduced in order to be able to join the basic elements corresponding to the three-dimensional or hierarchical models with the first-order model of piezoelectricity.
The example of the first transition geometry I, presented in Figure 1, includes one vertical edge (vertex nodes: a 1 , a 4 ) corresponding to the 3D-based first-order shell theory [13,54] within the mechanical field. The other part of this field corresponds to the 3D-based hierarchical shell [13,14] or three-dimensional elasticity models [13,53]. Above, the following linear vertex (a 2 , a 3 , a 5 , a 6 ) nodes and generalized: horizontal mid-edge (a 7 ,a 8 , . . ., a 12 ), vertical mid-edge (a 13 , a 14 ), mid-base (a 15 , a 16 ), mid-side (a 17 ,a 18 ,a 19 ) and middle (a 20 ) nodes are applied within the geometry and displacement fields. For the element I, the electric field corresponds to the hierarchical symmetric-thickness or three-dimensional dielectricity models [11,16]  It is worth mentioning that the generalized (varying-order) nodes of both types may contain multiple degrees of freedom along the edges, within the bases and sides, and in the interior of the element. The related finite element approximations within the proposed element are all based on the general rules proposed by Demkowicz and others [1,2].
In the example of the second transition geometry II, shown in Figure 2, one deals with two vertical edges (nodes a 1 , a 4 and a 2 , a 5 ) and the side (nodes a 7 , a 10 ) between these edges, all corresponding to the first-order mechanical model. The rest of the element (nodes a 3 , a 6 , a 8 , a 9 , a 11 , a 12 , a 13 , a 14 , a 15 , a 16 , a 17 , a 18 ) conforms to the hierarchical shell model or the three-dimensional model of elasticity. These mechanical models are coupled with the symmetric-thickness hierarchical or three-dimensional models of dielectricity. The same nodes b 1 , b 2 , ..., b 21 , as for the geometry I, are employed in this case.  (24)

Replacement of the global quantities introduced in
where e is the element number, while e q hpq and e ϕ hπρ stand for the element, displacements and potential, dofs vectors.
The element stiffness matrix is defined in the standard way, namely: with e B relating strains ε with the element displacement dofs e q hpq and, hence, named the element strain-displacement matrix. The matrix D represents the matrix of elasticities for the isotropy case under consideration, while the term det(J) is the determinant of the Jacobian matrix transforming the global coordinates x into normalized ones ξ = [ξ 1 , ξ 2 , ξ 3 ] T . Note that the applied normalized prismatic geometry of the element e can be characterized In the case of the thick-or thin-walled geometry, the first two normalized directions coincide with the longitudinal shell directions and the third normalized direction is the same as the transverse shell direction.
The definition of the element characteristic matrix of dielectricity reads where e b allows expressing the electric field E with the element electric potential dofs e ϕ hπρ . It is called the element field-potential matrix. The matrix γ is the constitutive constant matrix of orthotropic dielectricity under constant strain.
The element characteristic matrix of piezoelectricity coupling the mechanical and electric fields is defined as follows: with C being the constitutive matrix of orthotropic piezoelectricity under constant strain. The element nodal forces vector due to the mass loading f is defined in the standard way, i.e., with e N standing for the standard shape functions matrix corresponding to the vectorial displacement field.
In addition the element nodal forces due to the surface loading p is defined in the standard way, for the bases and lateral sides, respectively: for the element bases and sides, respectively, with e n denoting the standard shape function matrix suitable for the finite element approximation of the scalar electric potential field.
In the next subsection, we will present the changes that have to be implemented in the above finite element formulation depending on the applied transition model.

Strain-Displacement Matrix
Three cases of the classical, modified, and enhanced transition elements corresponding to the analogous transition models are of interest in this section.

The Transition Element Based on the Classical Model
The strain-displacement matrix for the classical transition model is assumed in the form corresponding to the three-dimensional elasticity and the hierarchical shell models, for the three-dimensional and thick-or thin-walled geometries, respectively. This means that the three-dimensional strain and stress state is applied for this classical model. This state is characterized by ε = e B e q hpq (43) and reflects the relation between the local strains ε and global nodal displacement dofs e q hpq . The local strain directions are coincident with two tangent (the first two or two longitudinal) and one normal (the third or transverse) directions within the structure.
The strain-displacement matrix e B will also be denoted as The above blocks B k ≡ B 1 k of the strain-displacement matrix correspond to the element displacement dofs at the element nodes and are composed of the following components: Any term B k,mn , m = 1, 2, .., 6, n = 1, 2, 3 of each block B k , where m corresponds to each of six components of the strain vector and n to each global direction of the displacement dof vector, can be determined as equal to: B k,1n = β k,1 θ 1n ; B k,4n = β k,1 θ 2n + β k,2 θ 1n B k,2n = β k,2 θ 2n ; B k,5n = β k,2 θ 3n + β k,3 θ 2n with the terms θ pn representing components of the transformation matrix θ which converts the local Cartesian directions p = 1, 2, 3 of the displacement dofs into the global ones. The terms β k,p can be calculated from where the quantities P k,r = ∂N k /∂ξ r are the partial derivatives, with respect to the normalized directions r = 1, 2, 3, of the shape function N k from the corresponding shape function matrix. The auxiliary quantities A rp are equal to: with the terms j rn representing terms of the inverse Jacobian matrix J −1 and θ np being the components of the transposed transformation matrix θ T . The latter matrix transforms derivatives with respect to the global Cartesian directions into derivatives with respect to the local Cartesian directions. The terms of the transformation matrices can be obtained from: where the unit vectors components u p = U p /U, v p = V p /V, w p = W p /W, p = 1, 2, 3 can be obtained from the vectors U, V and W equal to: where the partial derivatives ∂x n /∂ξ r = j nr , n, r = 1, 2, 3 are the terms of the Jacobian matrix J.

The Modified and Enhanced Transition Elements
In the case of these two types of the transition element, the constitutive relations (29) and (30) hold within the coupled mechanical field. The strain and stress state changes from the three-dimensional one to the plane stress state. This can be expressed with the following relation where the components Due to the plane stress assumption (21), any block k of the second component needs the following modification: where the terms B k,mn , m = 1, 2, 4, 5, 6, n = 1, 2, 3 are defined in accordance with (46)- (50). Because of the modified definition (51) of the strain vector ε, also any block k of the strain-displacement matrix has to be modified in accordance with The blending function is defined in accordance with α = α(ξ 1 , ξ 2 , 1 2 ), i.e., as a function of the two normalized directions, ξ 1 and ξ 2 , tangent to the mid-surface defined with ξ 3 = 1 2 . Note that, for α ≡ 1 and α ≡ 0, the definition (54) transforms into (45) and (53), suitable for the basic (three-dimensional and first-order) piezoelectric elements, respectively.

Field-Potential Matrix
For any type of the transition models: (16), (23), (24), and (30) plus (32), the definition of the electric field vector E remains unchanged. Because of this also the field-potential matrix has the same form for the classical, modified, and enhanced transition elements. This form corresponds to the three-dimensional or hierarchical symmetric-thickness dielectric models. For these models, the following relation holds: Above, the locally defined electric field E is expressed with the scalar potential dofs where the blocks b l of the field-potential matrix correspond to the electric potential dofs of the element. The number of these blocks may be different to the number of the displacement dofs as the electric potential field is hπρ-interpolated, while the displacement field is hpqinterpolated. The form of these blocks is: The terms b l,p , where p = 1, 2, 3 are the local Cartesian directions, are equal to Note that the quantities p l,r = ∂n l /∂ξ r , r = 1, 2, 3 are the derivatives of the terms n l of the shape function matrix with respect to the normalized directions r = 1, 2, 3. The shape function n l corresponds to the electric potential dof l at the element node. The auxiliary terms A rp , r, p = 1, 2, 3 can be calculated from i.e., analogously as in the case of the strain-displacement matrix. As before, the terms θ np belong to the transformation matrix θ which transforms the local Cartesian directions p = 1, 2, 3 into global Cartesian ones n = 1, 2, 3. The quantities j rn are the terms of the inverse Jacobian matrix J −1 again.

Constitutive Constant Matrices
In the case of the classical, modified, and enhanced transition elements, the constitutive, elasticity constant matrix from (37) takes the same form (13).
The classical transition element requires the dielectricity and piezoelectricity constant matrices, in (38) and (39), in the form defined by (18) and (19), respectively. As far as the modified and enhanced transition elements are concerned, one has to apply, in relations (38) and (39), the definitions included in (32) and (30), respectively.

Blending Functions
As said above, blending functions describe the linear change of the stress state within transition elements from the three-dimensional to plane stress state. This change manifests itself by the presence of function α = α(ξ 1 , ξ 2 ) in the strain definition of the modified and enhanced transition elements and in the piezoelectric and dielectric constant matrices for these two transition elements. Note that ξ i , i = 1, 2, 3 are the normalized directions within an element and ξ 3 defines location of the reference surface within the element.
The blending function corresponding to the transition element of the geometry I is presented in Figure 3. The function is represented by the triangle spanned at the vertices a 1 , a 5 , a 6 and takes the value of 0 for the vertical edge 1 corresponding to the mechanical field of the first-order in the transverse direction ξ 3 and to the electric field of an arbitrary order in this direction. The blending function value of 1 corresponds to the mechanical and electric fields of arbitrary orders for the vertical edges 2 and 3. This function can be expressed with the affine coordinates λ 1 , i = 1, 2, 3 of the triangular element base or the normalized coordinates ξ 1 and ξ 2 as where λ 1 = 1 − ξ 1 − ξ 2 , λ 2 = ξ 1 and λ 3 = ξ 2 . The geometry II of the transition element and the respective blending function are presented in Figure 4. The function is represented by the triangle spanned at the vertices a 1 , a 2 , a 6 and equals 0 for two edges 1 and 2 corresponding to the mechanical field of the first-order in the transverse direction. This value is equal to 1 for the vertical edge 3 where the mechanical field is of arbitrary order in the transverse direction ξ 3 . The order of the electric field is arbitrary for all three vertical edges. The function can be expressed in the following way: It should be noted that due to introduction of the linear blending functions into the vectors and matrices from (29)-(32), the polynomial order the integrands of (37)- (39) increases. This needs the respective increase of the number of the Gauss points in the numerical integration of these integrands.

Shape Function Matrices
The shape function matrices where the blocks N k are diagonal and correspond to the vectorial displacement dofs k at the element node: In the case of the electric potential, the following standard form of the shape function where the blocks n l , corresponding to the scalar electric potential dof l at the element node, reduce to: It is worth noticing that the forms (44)- (45) and (64)-(65) corresponding to three proposed transition models are exactly the same as in the case of the three-dimensional and hierarchical models of piezoelectricity.

Kinematic Constraints within the Elements
In this section, the algorithms for imposition of the constraints of no elongation of the normals to the mid-surface, and for application of the constraints imposed on the generalized (varying-order) displacement dofs, as well, will be presented.

The Classical and Modified Elements Stiffness Matrix Modification
In this case the constraints of no elongation of the normals are assumed on the boundary R e between the first-order and transition elements. These constraints are imposed with the penalty method where the parameter P is employed. It represents the inverse of the penalty parameter and is a number sufficiently large to overwrite the stiffness matrix terms with the constraint terms and to replace the finite element equation with the constraint equation. Our approach takes advantage of the presentation of the method given in Reference [1,2]. Imposition of the constraints under consideration needs the following modification of the stiffness matrix: where e k P is called the penalty matrix. Its terms include the mentioned parameter P.
Note that no modification of the force vectors e f V + e f S due to penalty method is required as one deals with zero constraints in the paper, e.g., the constraints of no elongation of the normals to the mid-surface. The same refers to other constraints presented in this paper.

The arithmetic operator matrices
where the two blocks consist of the following component sub-blocks: The top and bottom dofs, k and l, are associated with the nodes belonging to R e = R| V e only.
In the case of the transition element geometry I (Figure 3), the top dofs k are associated with the node a 4 , the bottom dofs l with the node a 1 , the other dofs m with the nodes a 2 , a 3 , a 5 , a 6 , a 7 , ..., a 18 . In the case of the geometry II, one has the following assignments: the top dofs k are associated with the nodes a 4 , a 5 , a 10 , the bottom dofs l with the nodes a 1 , a 2 , a 7 , the other dofs m with the nodes a 3 , a 6 , a 8 , a 9 , a 11 , a 12 , a 13 , ..., a 20 .

Modifying Operators
Following the above assignments, one defines Within the block R 1 , the following overlapping sub-blocks can be distinguished Sub-blocks R kk , R kl are responsible for summation of the top and bottom displacement dofs, while R lk and R ll for subtraction of these dofs. In turn, the block R 2 possesses the following structure based on the identity sub-blocks I This means that neither the summation nor the subtraction is performed for the dofs m. One more remark concerns the sums and differences of displacements present in the constraints (27). In the case of the symmetric-thickness geometry, these sums define mid-surface displacements. The first two differences, after division by the thickness t, give rotations of the normal to the mid-surface, while the third difference determines a transverse elongation of the corresponding normal. Note that only the latter difference needs to be constrained (equal to zero) in accordance with the last Equation (27). Because of this, multiplication of the penalty stiffness by the diagonal matrix Z is necessary. This matrix possesses the following structure: The component block Z 1 is composed of the diagonal blocks Z k and Z l corresponding to the degrees of freedom k and l: where the degrees of freedom k and l correspond to the sums and differences of the top and bottom displacement dofs. The second component block has the simpler form which results from the fact that, for the degrees of freedom m, other than the sums and differences of the displacement dofs, the penalty terms are not necessary.

Penalty Stiffness
The respective penalty stiffness e k P , present in (66), is defined in the directions coincident with the constraint directions of which the first two are tangent and the third one is perpendicular to the mid-surface. The stiffness form corresponds to the displacement dofs that represent the mentioned sums and differences of the top and bottom dofs. The definition of the penalty stiffness reads: The penalty parameter has to be larger than the terms of e k M and the terms of e k C , i.e., P max{(k M ) ij ; i, j = 1, 2, ..., n M }, where (k M ) ij represent terms of the stiffness matrix and n M stands for the corresponding number of displacement dofs. Additionally, one has P max{(k C ) ij ; i = 1, 2, ..., n M ; j = 1, 2, ..., n E }, with (k C ) ij being terms of the piezoelectricity matrix and n E denoting the corresponding number of electric potential dofs.
The above form (75) of the penalty matrix is suitable for the geometry I of the classical and modified transition elements. In the above equation, R e stands for the element side being the subset of the boundary R between the first-order and transition model: R e ⊂ R. The normal vector n is defined as perpendicular to the mid-surface S m , i.e., n ⊥ S m . It also approximately satisfies the condition n ∈ R e -approximately due to the finite element approximation of the element geometry. Note that two tangential directions and the third normal direction are coincident with the local directions present in the constraint Equation (27). On the element level, they are defined by the relations (49) and (50). Note also that the presence of the vector n in the above definition results in imposition of the penalty terms in the third (normal) direction only, i.e., the constraints are applied in this direction only. The shape function matrix terms are non-zero on the boundary R e . Hence, the non-zero penalty terms appear for the element nodes placed on this boundary.
In the case of the classical and modified transition elements of the geometry II, the above definition has to take into consideration the fact that the constraints (27) where, due to one dimensional integration, the value of der(J) is applied for transformation of the global directions into the normalized direction ξ 3 . As previously, the constraints concern the differences of the top and bottom displacement dofs. The non-zero penalty terms appear for the element nodes and the corresponding dofs placed on the boundary L e , in the transverse direction only.

The Enhanced Transition Element
Two types of constraints are present in the case of the corresponding transition model. The first type are the gradually changing constraints of no elongation of the normals to the mid-surface, represented by the last but one term of (34). The considered type of the element requires also the additional constraints corresponding to the gradual change of the transverse order of approximation from q ≡ I = 1 to an arbitrary value q ≡ I, where I is the hierarchical model order, and represented by the last term of (34). As the constraints of both types are fully active on the boundary with the first-order model (α = 1) and gradually switched off (α → 0) towards the boundary with the hierarchical or three-dimensional model, the nodal (collocation) values of the complement of the switching functions to unity 1 − α will appear in the penalty stiffness definitions so as to activate the penalty terms.

The Modified Stiffness Matrix
The relations (66)-(68) are still valid for the enhanced transition elements and constraints of no elongation of the normals to the mid-surface S m . However, the degrees of freedom k and l corresponding to the top and bottom displacements, and the other dofs m, as well, are associated with different nodes now because the constraints of no elongation of the lines normal to the mid-surface are applied in the entire element volume V e . In the case of the geometry I, the top dofs k are associated with the nodes a 4 , a 5 , a 6 , a 10 , a 11 , a 12 , a 16 the bottom dofs l with the nodes a 1 , a 2 , a 3 , a 7 , a 8 , a 9 , a 15 , and the other dofs m with the nodes a 13 , a 14 , a 17 , a 18 , a 19 , a 20 . The geometry II needs the following assignments: the top dofs k correspond to the nodes a 4 , a 5 , a 6 , a 10 , a 11 , a 12 , a 15 , the bottom dofs l to the nodes  a 1 , a 2 , a 3 , a 7 , a 8 , a 9 , a 14 , and the other dofs m to the nodes a 13 , a 16 , a 17 , a 18 . These assignments hold also for the constraints related to the change of the transverse approximation order. In this case, the stiffness matrix modification is described with the relation Comparing (66) and (77), one can notice that there are no matrices e R and e R −1 present in the above relation. This means that the replacement of the top and bottom dofs with their sums and differences and inversely is not necessary for the constraints of the second type.

The Modifying Operators
In the case of the constraints of no elongation of the normals, the definitions of the operators e R and Z, (69)-(71) and (72)-(79), respectively, are still valid. However, the new assignments of the top k, bottom l, and other m dofs to the element nodes for geometries I and II have to be taken into account while calculating the corresponding operator matrices. In the case of the constraints of the second type, which concern the other dofs m only, the following definition is appropriate As it can be seen, the first component block of (78) is Z 1 ≡ 0. The second component block has the form which results in appearance of the penalty terms for all dofs m in all three global directions.

The Penalty Matrices
The constraints of no elongation of the normals require the following definition of the penalty stiffness matrix: where normal vector is n ⊥ S m , with S m being the mid-surface within the symmetricthickness geometry. The presence of this vector results in application of the constraints in the normal direction only. In turn, the presence of the switching function I−α guarantees the gradual change of the penalty terms. It employs the identity matrix I and the matrix of nodal values of the blending function: and the nodal blocks α m are diagonal and correspond to the dofs m, other than the top and bottom ones: and where α m ∈ 0, 1 . Its value depends on the longitudinal location of the node, the dof m is associated with. In other words, α m = α(ξ 1,m , ξ 2,m ), where ξ 1,m and ξ 2,m represent first two normalized coordinates for the dof m.
In the case of the constraints of the second type one has e k P = These constraints also change gradually and are applied in all three global directions for degrees of freedom m of the element. The definitions (81) and (82) of the switching function remain unchanged. It should stressed that the above definitions (80) and (83) are valid for both the transition elements acting between the piezoelectric elements based on the first-order shell theory and the elements based on either the hierarchical shell theory or the three-dimensional theory. The difference is that, in the first case, the entire surface S m (characterized with ξ 3 = 1 2 ) represents the shell mid-surface, while, in the second case, it represents the reference surface only (defined with ξ 3 = 1 2 ). The reference surface becomes the mid-surface only on the boundary with the first-order model R ∩ S m .

Numerical Tests
Two aspects of performance of the proposed transition piezoelectric elements will be addressed in this section. The first one is the ability to remove high stress gradients between the 3D-based transition and 3D-based basic elements in the chosen model problems. The mechanical models of the basic elements conform to the 3D-based first-order shell theory and the 3D-based hierarchical shell or three-dimensional theory. The electric field is characterized with the same dielectric model conforming to either three-dimensional or 3D-based hierarchical symmetric-thickness theory.
The second aspect will be comparison of the solution convergence for the same model problems as in the first aspect with three versions of the transition elements employed. The model problems will concern the bending-dominated plate and shell and membranedominated shell. These problems manifest different and numerically demanding stress and strain states. Explanation of the differences between these three cases of the mechanical fields can be found in Reference [53]. The electric field and electric displacement states will be typical for piezoelectricity.
It should be stressed that our model problems presented below serve the numerical assessment of the presented algorithms. Because of that the loads and charges are assumed so as the mechanical and electric parts of the potential electro-mechanical energy (or potential co-energy) are of the same order -the case very demanding from the numerical point of view as the sign of the electro-mechanical potential energy may change in the convergence studies where the error of the energy is displayed as a function of the number of dofs in the problem. Note that, in the two technologically important actuating and sensing modes of piezoelectric operation, one part of the co-energy dominates the other one. Because of this, such modes are less demanding numerically.

Model Structures and Their Geometry
Three piezoelectric domains corresponding to the bending-dominated plate, bendingdominated half-cylindrical shell and a half of the membrane-dominated cylindrical shell are presented in Figures 5-7. The plate dimensions are l × l, the length of both shells is equal to l, while the circumferential length of a half of the cylinder is l ≈ 2πR with R being the radius of the mid-surface of the shells. All the structures are of the same thickness t. The following values are set in our tests: l = 3.1415 × 10 −2 m, R = 1.0 × 10 −2 m, t = 0.03 × 10 −2 m.

Material, Load, Charge and Boundary Conditions
All three structures are made of the same piezoelectric material of material constants typical for piezoceramics [57]. The Young modulus is E = 0.5 × 10 11     The mechanical boundary conditions are as follows. The plate is clamped along its four edges. The bending-dominated shell has two straight edges clamped and two curved edges free. In the case of the membrane-dominated shell, there is no rotation along the curved edges. The electrical boundary conditions assume grounding along: four edges of the plate, along two straight and two curved edges of the half-cylindrical shell, and two curved edges of the cylindrical shell.

Discretization and Methodology
Due to symmetry of the geometry, loading, electric charge, and boundary conditions, only quarters of the bending-dominated domains and an octant of the membranedominated domain are taken into account in or numerical tests and presented in the figures below. The numerical settings are as follows. We apply the meshes 4 × 4 × 2 of the prismatic elements in the displayed quarter or octant parts of the structures in our tests. The longitudinal order of approximation within the mechanical and electric field is the same and equal p = π = 1, 2, . . . , 8. The transverse orders of approximation of both fields are equal to q = ρ = 2 in the hierarchical shell zone of the domains, q = ρ = 1 in the first-order shell zone, and the changing (from 2 to 1) order in the transition zone constituting a layer composed of couples of prismatic transition elements. The number of the first-order and hierarchical (three-dimensional) elements changes from the entirely hierarchical shell structure to the shell structure composed of the first-order shell elements only. The intermediate (or mixed) case includes a single layer of piezoelectric transition elements. The exemplary cases for the plate and shells are presented in Figures 8-10.  The presented meshes are composed of one layer of hierarchical elements (yellow and marked as MI in the figure), one layer of the transition elements (green and marked as RM/MI), and two layers of the first-order elements (blue and marked as RM in the figure).
Our parametric studies are limited to changes of the longitudinal approximation orders p ≡ π, the ratios of the numbers of elements of various types, and the types of the transition elements for three model problems. The mentioned ratios are equal to r = 1.0, 0.5, 0.0 and represent the number of the first-order elements divided by the sum of the transition and hierarchical elements along the characteristic longitudinal dimensions of the structures.

Ability to Remove High Stress Gradients
In order to present the ability of three (classical, modified, and enhanced) types of the transition elements to remove the high stress gradients between the transition elements and the neighboring basic (hierarchical and first-order) elements, exemplary (p ≡ π = 6) distributions of the global normal stress σ 33 (marked as szz) for three model problems and three types of the transition elements are presented. These distributions correspond to the mixed models from    Figure 13. Transverse stress-the bent plate with the enhanced transition elements.         Figure 18. Third normal stress-the membrane with the classical transition elements.   Figure 19. Third normal stress-the membrane with the enhanced transition elements.
In order to quantify the presented graphical results concerning stress gradients, the Table 1 is presented. Therein, gradients (jumps/differences) of the global stress σ 33 between the face centers of the neighboring elements are displayed for three model structures. Five models for each of three structures are included: the entirely hierarchical (MI) shell model (r = 0.0), three mixed (RM/MI) models (r = 0.5) with the classical, modified, and enhanced (respectively, TR1, TR2, TR3) transition models, and entirely first-order (RM) model (r = 1.0). In the case of the bending-dominated (bent) plate, the chosen elements 31 and 24 correspond to 3D model and TR1 or TR2 or TR3 model. Both elements 24 and 23 correspond to either TR1 or TR2 or TR3 model, while the elements 23 and 16 to either TR1 or TR2 or TR3 model and RM model. In the case of the bending-dominated (bent) shell, the corresponding couples of elements are: 25 and 28, 28 and 27, and 27 and 30. In the case of the membrane (membrane-dominated shell), the corresponding couples are: 6 and 13, 13 and 14, 14 and 21.          In the case of the plate, the third normal stress σ 33 corresponds to the transverse normal stress. In the entirely hierarchical (MI) shell model, the observed values of gradients are small and due to finite element discretization. In the entirely first-order (RM) model, the jumps are extremely small (correspond to the numerical zero), as the plane stress assumption holds in all elements. In the case of the mixed TR1 model, the jump between the transition and first-order elements 23 and 16 results from the three-dimensional stress state and plane stress state in the neighboring elements. Such jumps are diminished to numerical zero for the cases of the TR2 and TR3 models, as the plane stress assumption holds on both sides of the neighboring faces of the transition and first-order elements. In both shell cases, the presented jumps of σ 33 result from both discretization and plane stress assumption. The exact analysis for the transverse direction requires taking jumps for the stress components σ 11 and σ 13 = σ 31 into account. Such an analysis (not revealed here) leads to exactly the same observations as in the case of the plate. For the membrane case, where the contribution of σ 33 to the transverse normal stress is higher than the influence of two other stress components, the stress jumps in σ 33 between elements 14 and 21 are smaller in the case of the TR3 model than in the cases of the TR2 and TR1 models. In the bent shell case, the jumps in σ 33 between elements 27 and 30 grow from TR1 through TR2 to TR3, but this leads to zero transverse normal stress values in both neighboring elements for the TR2 and TR3 models. Demonstration of this needs inclusion of the stress components σ 11 and σ 13 = σ 31 again.

Solution Convergence of the Problems
In the convergence studies, the logarithm of the approximation error in potential energy is presented as a function of the logarithm of the number of dofs N of the problem. The corresponding curves are p-convergence ones as the number of dofs in the problems increase due to the change of the longitudinal order of approximation p within elements. In order to present the solution convergence for the three model problems with three types of the transition elements applied, three threesomes of pictures are displayed. The first three pictures (Figures 20-22) correspond to the bending-dominated plate problem, the second threesome (Figures 23-25) to the bending-dominated shell problem, while the last three drawings (Figures 26-28) to the membrane-dominated shell. In the figures, the first, second and third drawings correspond to the cases of the classical, modified, and enhanced transition elements applied.
The mentioned approximation error is defined as the difference of two energy measures of which E corresponds to the numerical solution under consideration and the reference solution E r plays a role of the exact solution, namely: The energy measures represent the sum of absolute values of the mechanical and electrical parts of the electro-mechanical potential energy, i.e., Due to stationarity of the solutions u and φ, resulting from (8), the potential energy can be expressed with the components of the electro-mechanical co-energy (strain, electric field, and coupling energies) corresponding to the bilinear forms B, b and C from (9)-(11). The solutions for displacements u and electric potential φ are approximated with use of the global nodal vectors of the displacement dofs q hpq and potential dofs ϕ hπρ obtained from solution of the problem (35 Figure 28. Convergence curves-the membrane shell with the enhanced transition elements. The proposed measure of the approximation error was applied as it guarantees monotonicity of the energy approximation error with changing value of p ≡ π, while the simple sum of both parts of the energy (equal to the electro-mechanical co-energy) may not. Note that such monotonicity is a must in the case of error-controlled adaptivity suggested in Reference [11,12].
In order to calculate the errors, the exact values E r of the energy measures are replaced with their best numerical approximations obtained from the so-called over-killed meshes, i.e., the meshes for which the longitudinal order of approximation is equal to p = π = 9.
In the case of two bending-dominated problems the convergence for the enhanced models is slightly higher than for the modified transition elements employed in the mixed models. Additionally, the latter elements produce slightly better convergence than the classical ones. In the case of the membrane-dominated problem, only the enhanced transition elements deliver convergence curves between the curves obtained for the basic models. The confirmation is presented in the Table 2, where the absolute values log(E r −E) of the error and relative error values (E r −E)/E r (both for p = 8, i.e., E ≡ E 8 ) are included. In addition, the averaged convergence rates are given and equal to −log(δ 8 /δ 7 )/log(N 8 /N 7 ), δ 8 = E r −E 8 , δ 7 = E r −E 7 , where E 8 , N 8 and E 7 , N 7 represent the energies and numbers of dofs for p = 8 and p = 7.

Conclusions
Three types of the transition piezoelectric elements joining the basic piezoelectric elements are possible. The joined basic elements are the three-dimensional piezoelectric elements (or hierarchical symmetric-thickness piezoelectric elements) and the first-order symmetric-thickness piezoelectric elements.
The first (classical) transition element guarantees continuity of the displacement and electric potential fields between the basic elements. It also guarantees that the assumption of no elongation of the normals to the mid-surface holds on the boundary between the transition elements and the first-order elements.
The second (modified) transition element, additionally, guarantees the gradual change from the plane stress to three-dimensional stress state between the basic elements of the first-order and three-dimensional (or hierarchical) character, respectively.
The third (enhanced) transition element, additionally, assures the gradual change of the assumption of no elongation of the normals to the mid-surface between the basic elements. The latter assumption holds on the boundary with the first-order element and is gradually switched off towards the boundary with the three-dimensional (or hierarchical) element. On this boundary, this assumption is completely off.
The gradual change of the stress state from the plane to three-dimensional is based on introduction of the blending functions that determine contributions of the basic models to the transition model. Implementation of this idea on the element level needs modifications of the strain vector, and the piezoelectric and dielectric constant matrices as well.
The gradual change of the assumption of no elongation of the normals needs utilization of the blending functions, which play the role of gradually switching functions. On the element level, the nodal values of the complement (to unity) of the gradually switching functions appear in the penalty stiffness matrix responsible for the gradual change of the constraints.
Comparing the stress distributions within the transition zones for three model problems, one can notice that high stress gradients appear on the boundary between the classical transition elements and the first-order elements. In the case of the enhanced transition elements, the gradual change of the stress state can be observed. In the case of the modified transition element, the stress distribution is intermediate with respect to the previous two distributions.
It can be observed that the discrete models of the three analyzed model structures, with three types of the transition elements employed, produce the convergence curves of slightly higher convergence for the case of the modified elements in comparison to the classical elements. The highest convergence, however, is obtained in the case of the enhanced transition elements.
It can be concluded that the enhanced transition element should be recommended in the analysis of mixed models of piezoelectrics, as this element gives the same or better convergence than the other transition models, and it removes the high stress gradients between the classical transition elements and the first-order piezoelectric elements.
Even though all three transition models can be applied to error-controlled adaptivity due to monotonicity of the applied error measure, the enhanced transition elements deliver convergence curves with the error level and convergence rates just between the curves for the basic models.
The presented research needs the following next steps. Firstly, the error estimation methods and adaptive procedures for the new transition elements have to be prepared and verified. Next, application of the proposed transition models to parametric and/or adaptive analysis of electro-mechanical systems should be performed and critical factors influencing quality of such an analysis should be determined.