The p-version of the finite element method for three-dimensional curved thin walled structures

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


INTRODUCTION
The p-version of the ÿnite element method has turned out to be an e cient discretization strategy for many linear elliptic problems. To name a few consider e.g. the Poisson equation, the LamÃ e equations, the Reissner-Mindlin problem, etc. It was shown by many authors, that for this class of problems the p-version is superior to the classical h-version approach [1][2][3][4][5][6][7]. Combined with a proper mesh design, the p-version shows an exponential rate of convergence in energy norm in the preasymptotic range [8]. If a priori knowledge of the solution is used to Figure 1. Set of one-dimensional standard and hierarchic shape functions for p = 1; 2; 3. construct (geometrically) reÿned meshes towards points or lines of singularity or to resolve boundary layers, an approximation with an-in an engineering sense-acceptable error is easily obtained due to the preasymptotic range of exponential rate of convergence. In addition to high accuracy the p-version includes some further advantages. It was theoretically and numerically shown, that the p-version is free of locking e ects, if the polynomial degree is chosen to be moderately high [6; 9]. This includes, e.g. shear locking e ects appearing when Reissner-Mindlin problems are considered or Poisson locking, which plays an important role in elastoplastic problems [10][11][12]. Using the blending function method, curved boundaries can be easily considered without increasing the number of elements [13-15; 7]. Many structures can therefore be discretized using a coarse mesh, as the basis for discretizations with higher order Ansatz spaces.
In this paper an implementation of the p-version will be presented, which allows to consider three-dimensional structural problems with almost arbitrarily curved surfaces. The implementation is based on a hexahedral element formulation, being able to vary the polynomial degree for the three local directions ( ; Á; ) as well as for the Ansatz of the three components u(x) = (u x ; u y ; u z ) T of the displacement ÿeld. The use of anisotropic Ansatz spaces leads to very e cient approximations especially for thin walled structures. It will be demonstrated, that the considered implementation of the p-version allows to predict the three-dimensional structural behaviour of plates and shells with approximately the same amount of degrees of freedom as with dimensionally reduced implementations, yet signiÿcantly more accurate due to the three-dimensional model.
The outline of the paper is as follows: First we will present the implemented Ansatz spaces, based on hexahedral elements. In Section 3 a short introduction to the blending function method and the client-server concept in our software structure will be given. This approach is similar to a geometry representation used by Dey et al. [16;17] for the solution of Poisson's equation. The e ciency of the p-version for thin walled structures will then be demonstrated by several numerical examples in Section 4.

The one-dimensional hierarchic basis
Following SzabÃ o and BabuÄ ska [7] it will ÿrst be shown how hierarchic basic functions can be implemented up to any desired polynomial degree. Let us start with the standard ÿnite element basis (nodal basis) in one dimension on a standard element st = (−1; 1) (see Figure 1, left part). Obviously, every function representable by the standard basis can also be represented by the set of hierarchic basis functions (see Figure 1, right part). A principle di erence between the two bases is that in the hierarchic case all lower order shape functions are contained in the higher order basis.
Our two-and three-dimensional p-version implementation is based on the set of onedimensional hierarchic shape functions N i ( ) = i−1 ( ); i= 3; 4; : : : ; p + 1 with j ( ) = 2j − 1 2 −1 P j−1 (t) dt = 1 4j − 2 (P j ( ) − P j−2 ( )); j= 2; 3; : : : where P j ( ) are the well-known Legendre polynomials P k ( ) = 1 2 k k! d k d k ( 2 − 1) k ; k = 0; 1; : : : The linear functions N 1 ( ); N 2 ( ) are called nodal shape functions or nodal modes. Because In Reference [18] it was shown, that the condition number of the sti ness matrix for the Navier equations of elasticity is improved by an order of magnitude if a hierarchic basis of shape functions is applied. Furthermore, it is important to notice, that a hierarchic basis has an immediate consequence on the structure of the resulting sti ness matrix. If equations are ordered so that all linear modes get numbers 1 to n 1 , all quadratic modes get numbers n 1 + 1 to n 2 and so on, sti ness matrices corresponding to polynomial order 1 to p − 1 are submatrices of the sti ness matrix corresponding to polynomial order p. Shape functions for two-and three-dimensional Ansatz spaces can now be easily constructed, by simply forming the tensor product of one-dimensional hierarchic shape functions.

Hierarchic shape functions for hexahedral elements
Our implementation of the p-version in three dimensions is based on a hexahedral element formulation, using the Ansatz functions introduced by SzabÃ o and BabuÄ ska [7]. Hexahedral elements (see Figure 2) show some advantages, when being compared to tetrahedral and pentahedral element formulations: • Hexahedral element formulations lead to higher accuracy. • Hexahedral elements are especially well suited for thin walled structures. One local variable can be identiÿed to correspond with the thickness direction. Therefore it is possible to choose the polynomial degree in thickness direction di erent from those in in-plane direction. • The numerical integration of hexahedral elements can be easily performed using a Gaussian quadrature scheme.

Edge modes:
These modes are deÿned for each individual edge separately. If we consider e.g. edge 1, the corresponding edge modes read:

Face modes:
These modes are deÿned for each individual face separately. If we consider e.g. face 1, the corresponding face modes read:  ) can be varied in each local direction separately (see Figure 2). The di erence between the trunk space and the tensor product space is relevant for the face modes and the internal modes. For explanation we ÿrst consider the face modes, e.g. the modes for face 1. Indices i; j denote the polynomial degrees of the face modes in -and Á-direction, respectively.

Trunk space
Tensor product space i = 2; : : : ; p − 4 i = 2; : : : ; p j = 2; : : : ; p Á − 4 j = 2; : : : ; p Á k = 2; : : : ; p − 4 k = 2; : : : ; p i + j + k = 6; : : : ; max{p ; p Á ; p } The Ansatz space S p; p; q ( h st ) deÿnes an anisotropic set of shape functions being determined by two polynomial degrees p and q (see Figure 2). All shape functions being of higher order in and Á-direction are associated with the polynomial degree p. These shape functions correspond to the edges 1,2,3,4,9,10,11,12, to the faces 1 and 6 and to all internal modes. Shape functions for faces 1 and 6 are equal to the ones of the trunk space S p ; pÁ; p ts ( h st ) with p = p = p Á . q deÿnes the degree of all shape functions being of higher order in -direction, which are associated with the edges 5,6,7,8, with the faces 2,3,4,5 and with all internal modes. The modes corresponding to the faces 2,3,4,5 are equal to the ones of the tensor product space S p ; pÁ; p ps ( h st ) with p = p = p Á and q = p . Considering a polynomial degree p = q = p = p Á = p one observes, that the number of internal modes of S p; p; q ( h st ) is higher than the one of the trunk space S p ; pÁ; p ts ( h st ) and less than the one of the tensor product space S p ; pÁ; p ps ( h st ) (see Appendix A). Furthermore, it can be noted, that the internal modes of the space S p; p; q ( h st ) are polynomials in which the direction is preferred. Due to the built-in anisotropic behaviour of the Ansatz space S p; p; q ( h st ) it is important to consider the orientation of the local co-ordinates of a hexahedral element. In Figure 3 it is shown, how hexahedral elements should be orientated, when three-dimensional thin walled structures are to be discretized. The local co-ordinate of the hexahedral element corresponds with the thickness direction z. If the orientation of all elements is equal, then it is possible to construct discretizations where the Ansatz for the in-plane and thickness direction of thin walled structures can be treated di erently. The numerical examples in Section 4 will demonstrate that anistropic Ansatz spaces lead to e cient discretizations. Our implementation of the p-version allows not only to vary the polynomial degree for the three di erent local directions but also to choose a di erent degree for each primary variable. The following two examples illustrate how to deÿne a polynomial degree template p for a structural problem with three primary variables u = (u x ; u y ; u z ) T : will lead in case of the Ansatz space S p; p; q ( h st ) to the Ansatz: u x ∈ S 1; 1; 2 ( h st ); u y ∈ S 3; 3; 4 ( h st ) and u z ∈ S 5; 5; 6 ( h st ) This polynomial degree template p can be deÿned for each individual element or for a group of elements. Our ÿnite element code automatically takes care of interelement continuity of the Ansatz. If adjacent elements have di erent deÿnitions of p, then always the highest degree for the common edge and face modes is chosen. In Figure 4 (left part) the number of degrees of freedom of one hexahedral element with isotropic polynomial degree, i.e. for for a three-dimensional LamÃ e problem is plotted. As expected, the tensor product space supplies the highest number of degrees of freedom. In Figure 4 (right part) the ratio of internal modes to degrees of freedom is pictured. S p ; pÁ; p ps ( h st ) is the Ansatz space with the highest number of internal modes. As the internal degrees of freedom are purely local to the element, they can be eliminated by static condensation. This results in an increase of computation time on element level but drastic decrease of solution time because the condition number of the global sti ness matrix is strongly reduced. Several authors [19][20][21] have investigated these observations in detail, interpreting the internal mode condensation as a preconditioning procedure for iterative solvers.

THE BLENDING FUNCTION METHOD
One important di erence between h-and p-version ÿnite element methods lies in mapping requirements. Because in the p-version the element size is not reduced as the polynomial degree is increased, the description of the geometry has to be independent of the number of elements. This results in the necessity to construct elements with an exact representation of the boundary. The isoparametric mapping, used in standard ÿnite element formulations, can be seen as a special case of mapping using the blending function method [14-17; 7]. Following these ideas, element boundaries can be implemented as (almost) arbitrarily curved edges and faces.
Before describing an algorithm coupling a geometric model to the ÿnite element analysis and taking advantage of the blending function method, some basic concepts of element matrix is given by an integral of matrix functions M K ; M b , depending on: • Jacobian matrix, with its inverse and determinant, J; J −1 ; |J|.
Let now M be a generic matrix function corresponding to M K or M b . In general integration (8) and (9) is performed numerically in the local co-ordinate system ( ; Á; ) of a standard hexahedral element domain h st = [(−1; 1)×(−1; 1)×(−1; 1)]. Using e.g. a Gaussian quadrature the integral is replaced by a weighted sum Therefore, the terms N; C; b; Q; J −1 ; |J| have to be computed only at integration points ( l ; Á m ; n ). In the following, we will consider the blending function method and show, that Q and J can be evaluated at any interior point only from the knowledge of its local co-ordinates and information on the surface of the element. Our formulation of the blending function method follows closely the work of KirÃ alyfalvi and SzabÃ o [15]. Consider a hexahedral element, as pictured in Figure 5. X i = (X i ; Y i ; Z i ); i = 1; : : : ; 8 denote the global co-ordinates of the nodes. E i = (E ix ; E iy ; E iz ); i = 1; : : : ; 12 are functions which depend on local co-ordinates ( ; Á; ) and describe the shape of each edge. F i = (F ix ; F iy ; F iz ); i = 1; : : : ; 6 denote the functions 8 describing the shape of each face. The mapping function Q e ( ; Á; ) from local^= ( ; Á; ) T to global co-ordinates x = (x; y; z) T is obtained by The ÿrst term is the standard mapping of isoparametric eight-noded hexahedral elements. The second term is referred to as face blending (see Appendix B). Consider e.g. face 6 of the hexahedral element shown in Figure 5, where F 6 ( ; Á) describes the parametric mapping of the local ( ; Á)-plane to the surface of the element: At face 6, where = 1 the blending term (1 + )=2 equals 1 and therefore f 6 ( ; Á; = 1) describes the di erence between the curved and the bilinear face co-ordinates ( ; Á). Due to the blending term it is guaranteed, on the other hand, that this di erence (i.e. the function f 6 ( ; Á; )) decreases linearly to the opposite face 1, where = −1 such that f 6 ( ; Á; −1) = 0. The third term in Equation (11) corresponds to the edge blending. Considering e.g. edge 1 we have e i ( ; Á; ); i = 1; : : : ; 12 denote the di erence between the curved edge and the linear connection of the two end points, multiplied by a blending term (see Appendix B). The structure of the edge blending is similar to the face mapping, now with a blending term being linear in two variables. Because each edge belongs to two faces of a hexahedral element, a correction with respect to a straight line edge appears twice in the surface blending term (second sum) in Equation (11). Therefore, the corresponding edge blending term has to be subtracted. Substituting f i ( ; Á; ); i = 1; : : : ; 6 and e i ( ; Á; ); i = 1; : : : ; 12 in Equation (11) by inserting Equations (B2), (B1) of Appendix B and rearranging terms we ÿnally obtain Consider now the Jacobian matrix It contains the derivatives of the mapping function Q e with respect to the local co-ordinates ; Á and . From (14) it can be readily seen, that all coe cients of J depend linearly on the following three groups of geometric information: • the co-ordinates of the nodes X i for i = 1; : : : ; 8, • the 12 tangent vectors @E i @r where r = for i = 1; 3; 9; 11 r = Á for i = 2; 4; 10; 12 r = for i = 5; 6; 7; 8 (16) • the six tangential planes @F i @r ; @F i @s where r = ; s = Á for i = 1; 6 Therefore, the Jacobian matrix at any interior point x( ; Á; ) of the element can be computed only from nodal, edge and surface data and from the local co-ordinates ( ; Á; ) themselves. For given ( ; Á; ), edge and surface derivatives have to be sampled at points corresponding to the local co-ordinates as indicated in Figure 6, showing schematically the essential information on the edges and surfaces of a hexahedral element for computing the Jacobian matrix J.
Assuming now a tensor product integration scheme in (10) with k 3 integration points, 6k 2 corresponding sampling points on the surfaces and another 12k sampling points on the edges of a hexahedron can be collected. To compute tangential planes (17) and tangent vectors (16) in these points, one has to consider that the corresponding functions F i and E i are deÿned as  composed mappings in general. Figure 7 sketches these mappings for the upper face F 6 of the hexahedron (see also Figure 5). Deÿning u := (u; v) T = (u( ; Á; = 1); v( ; Á; = 1)) T and x := (x; y; z) T = (x(u; v); y(u; v); z(u; v)) T F 6 is given by Assuming that the mapping x is bijective, we denote by x −1 (E i ) the set of all points in the (u; v)-plane being the image of points on the edge E i under the mapping x −1 . As in general edges E i of a hexahedron are obtained by arbitrary geometric operations, their image x −1 (E i ) is not available in closed form and must be approximated by a pointwise iterative backtransformation. Once this approximation is constructed, the mapping u from ( ; Á)-to (u; v)-plane can be deÿned by (two dimensional) blending [7]. The derivatives (16) and (17) of edges and surfaces at the surface and edge sampling points with respect to the local co-ordinates can then be obtained by application of the chain rule to (18) and the corresponding expression for element edges.
Observing now this special structure of the mapping x, an implementation of the computation of element matrices in a distributed software system, strictly separating the geometric description from the deÿnition of shape functions is possible. Figure 8 sketches a client-server structure of this computation. At each integration point ( e l ; Á e m ; e n ) the mapped co-ordinates Q e ( e l ; Á e m ; e n ) and the Jacobian matrix J( e l ; Á e m ; e n ) are computed from the boundary representation data of the geometric model. We are using AutoCAD with its ACIS-kernel [31] and the ARX-interface to provide all geometric information. Surface co-ordinates and derivatives can be obtained either directly from the geometric modeller or from an intermediate interpolation surface, like the 'quasi-regional mapping' as described in Reference [15]. Details of these implementations are described in Reference [22].
For reasons of computational e ciency we do not envoke the interprocess communication shown in Figure 8 for each integration point individually, but collect all local co-ordinates, send them at once and receive computed data again by one message, only. The additional e ort due to the computation of the mapping function Q e , the Jacobian matrices J and the Summarizing, this blending function technique for mapping the geometry of p-elements o ers the possibility to completely separate all geometric computations involved in a ÿnite element analysis from the non-geometric part. This separation allows to design a distributed software system, where the geometric model of a CAD-program, although running as a different process or even on a di erent computer, is directly linked to a ÿnite element kernel. Using this software interface, the evaluation of the Jacobian matrix J at each integration point and thus the numerical integration of element matrices and load vectors is possible without any explicit knowledge on the types of surfaces or edges of the solid model. This software structure o ers the advantage of using all state-of-the-art CAD-techniques like geometric editing or parametric design in a ÿnite element analysis, immediately. The increase of e ciency for practical work may be dramatic, as such a system for computer integrated engineering relieves a user from the necessity to transfer geometric data from CAD to FEA, which is usually very time-consuming, even if only some geometric parameters of the model change.

Clamped plate under uniform load
In this ÿrst simple example of a clamped plate under uniform load we compare a Reissner-Mindlin model with three-dimensional h-and p-approximations and derive guidelines for an a priori deÿnition of the polynomial degree template p for plate-like structures.
The length of the quadratic plate is L = 12 and the thickness equals t = 0:35. The plate is loaded by a uniform pressure T z = −100 acting on the upper surface of the plate. Linear elastic material behaviour is assumed with E = 30 000 000 being Young's modulus and = 0:2 denoting Poisson's ratio. Considering the two-dimensional discretization based on the Reissner-Mindlin plate theory the warping coe cient accounting for non-uniform shear distributions is chosen to be Ä = 5=6.
Several discretizations for this problem will be investigated (see Figure 9). The two-dimensional discretization is based on a Reissner-Mindlin plate formulation with 25 quadrilateral p-elements utilizing the tensor product space S p;p ps ( q st ), see References [2; 3]. The three-dimensional discretizations for the p-version are based on the mesh consisting of 25 hexahedral elements with one element layer over the plate thickness. In order to resolve boundary layers the meshes for the p-version are reÿned towards the boundary. To draw a comparison to the h-version with trilinear hexahedral elements, discretizations with up to 22496 elements will be applied. The number of elements over the plate thickness varies in this case between one and six.
The reference value for the strain energy W = 53:09895558 of the three-dimensional problem is based on a Richardson extrapolation of values being obtained with a mesh consisting of 48 hexahedral elements in conjunction with the trunk space S p;p;p ts ( h st ) where p = 10; 11; 12. It may be interesting to note, that the two-dimensional Reissner-Mindlin converges to a strain energy of W = 53:18817430.
In order to get an impression of the di erence between the two-and three-dimensional discretizations, consider the stress distribution over plate thickness e.g. at x = y = 9:0; z of both models in Figure 10. The solutions are based on converged approximations, so that the in uence of the numerical error is negligible. The stress components xx ; yy and xy of the Reissner-Mindlin model coincide with the three-dimensional solution. Deviations of the Reissner-Mindlin model from the three-dimensional solution can be observed for zz and yz ; zx . These deviations are due to the kinematic and constitutive assumptions made in the formulation of this model. Further deviations of the Reissner-Mindlin approximation from the exact three-dimensional solution could be observed, if we would investigate the solution of both models close to the boundary. For a detailed discussion of boundary layer e ects, we refer to Reference [6].
Considering a comparison of the three-dimensional discretizations the strain energy is plotted in Figure 11. It can be seen, that the p-version with an Ansatz space S p ;pÁ;p ts ( h st ) (p = p = p Á = p is chosen in this example uniformly for all components) supplies with p = 4 and Furthermore it is obvious, that signiÿcant e ort can be saved, if anisotropic Ansatz spaces are used in the p-version. To get an indication for an economic p-distribution for platelike structures, we will ÿrst investigate the dependence of strain energy upon the polynomial degree template p using all three Ansatz spaces S p; p; q ( h st ); S p ;pÁ;p ts ( h st ); S p ;pÁ;p ps ( h st ). To test the sensitivity of the strain energy upon the discretization, p will be varied such that either the polynomial degree for the Ansatz u x and u y or for u z will be one order higher, i.e.  Figures 12-14, it is evident, that a discretization supplies less error, if the polynomial degree is raised rather for the de ection u z than for the displacement components u x ; u y . This is of course not surprising, because this example is a bending dominated problem where most of the strain energy is due to the de ection u z . Therefore, an accurate  approximation of the displacement component u z is more important than for the displacement components u x and u y .
Concerning the variation of the polynomial degree in ; Á or -direction, it has to be mentioned, that the p-version mesh is constructed such that the local coordinates ; Á and correspond with the global co-ordinates, i.e. coincides with the z-direction. Regarding the lower parts of Figures 12-14 it becomes obvious that-concerning the Ansatz spaces S p ; pÁ; p ts ( h st ) and S p ; pÁ; p ps ( h st )-for a wide range of p an increase of polynomial degree with respect to local co-ordinates and Á leads to more e cient discretizations than the ones, where the order for the local co-ordinate is raised. This e ect is even more evident, if we regard the Ansatz space S p; p; q ( h st ) (see lower part of Figure 14). This di erence between the behaviour of the three Ansatz spaces can be explained by the anisotropic nature of the set of shape functions of S p; p; q ( h st ). While the Ansatz spaces S  is less than that in thickness direction. Due to the built-in anisotropy, this is not the case for S p; p; q ( h st ), where the -direction is a priori preferred. Summarizing, it is important to handle the polynomial degree for the Ansatz of u x ; u y and u z di erently and to choose for u z a higher order approximation, e.g.: Furthermore the numerical investigation has shown, that it is also advantageous to choose the polynomial degree in ; Á-direction one order higher than in -direction: Figure 15. Comparison of the two-and three-dimensional discretizations for a clamped plate under constant load with isotropic and uniform distribution of the polynomial degree. Figure 16. Comparison of the twoand three-dimensional discretizations for a clamped plate under constant load, polynomial degrees as in Table I.
Combining these observations, we suggest a polynomial degree template for plate-like structures being of the form where p 2 ¿p 1 , p 4 ¿p 3 and p 1 ¿p 3 , p 2 ¿p 4 . Convergence towards the exact solution of the (three-dimensional) structural problem will of course only be obtained, if the lowest polynomial degree, i.e. p min = min{p 1 ; p 2 ; p 3 ; p 4 } tends to inÿnity. Considering yet plate-like structures, it is possible to gain an-in a engineering sense-acceptable error, where p min is limited to a certain value. This corresponds to accepting a modelling error inherent e.g. in all plate or shell theories, yet with the major di erence, that an increase of the polynomial degree guarantees a more and more accurate modelling of the three-dimensional structure. It should also be noted, that an optimal distribution of the polynomial degree can be obtained using anisotropic a posteriori error estimators outlined by Stein et al. [23].
In the following, the strain energies for computations with Ansatz spaces S  Figure 15 shows results of an isotropic and uniform distribution of the polynomial degree for the three-dimensional solution, whereas anisotropic Ansatz spaces were used in the computations depicted in Figure 16.
In Table I Figure 16. In in-plane direction ( -Á-plane) the polynomial degree is raised for u x and u y up to p = 7 and for u z up to p = 8. Concerning the Ansatz spaces S p ; pÁ; p ts ( h st ) and S p; p; q ( h st ) the polynomial degree in plate thickness direction ( ) is limited for u x and u y to p = 3 and for u z to p = 4. Only for the Ansatz space S p ; pÁ; p ts ( h st ) the polynomial degree  in -direction is raised for u x and u y up to p = 5 and for u z up to p = 6. A reason for this choice of p in -direction is the fact, that the Ansatz space S  ) supplies less shape functions in -direction and a higher polynomial degree is required to achieve approximately the same accuracy as with the Ansatz spaces S p ; pÁ; p ps ( h st ) and S p; p; q ( h st ). Comparing the results of the three-dimensional discretizations and the two-dimensional solution it becomes evident, that the use of anisotropic Ansatz spaces reduces the numerical e ort dramatically and that anisotropic S p ; pÁ;p ts ( h st ) as well as S p; p; q ( h st ) yields a strain energy being comparable to that of the Reissner-Mindlin model with a similar number of degrees of freedom.

Plate with columns
As a second example we consider a plate under uniform load supported by 9 columns, where signiÿcant three-dimensional e ects are to be expected near the intersections of plate and columns (see Figure 17). The thickness is set to t=L = 0:2=12 with L being the dimension of the plate. Each column has a cross-sectional area of 0:3 × 0:3 and a height of 3.0.
The structure is loaded by a uniform pressure T z = − 100 acting on the upper surface of the plate. Again, linear elastic material behaviour is assumed with E = 30 000 000 being  Young's modulus and = 0:2 denoting Poisson's ratio. In the two-dimensional case a shear correction factor Ä = 5=6 is chosen. Due to the geometric complexity of this structure and the large number of singularities in the exact solution this example is representative for engineering problems. Di erent discretizations in two as well as in three dimensions will be investigated. First we will consider the two-dimensional case, using the Reissner-Mindlin plate theory. We will apply the p-version as well as an adaptive h-version based on MITC4 elements. In both cases the columns will be modelled as elastic foundations with c w = 10 000 000 being the spring constant with respect to the de ection w (for a detailed description see Reference [4]).
In Figure 18 the discretization of the plate problem with 176 p-elements is sketched. At the boundary and reentrant corners the mesh is reÿned, allowing to resolve singularities and boundary layers of the exact solution. The adaptive h-version is based on the MITC4-element formulation [24] and an a posteriori error estimation of residual type [25; 4]. Figure 19 shows the initial mesh and a sequence of three adapted meshes. As expected, the error indicators lead to meshes being reÿned at reentrant corners and at columns modelled by elastic foundations.
In Figure 20 the relative error in energy norm for both the p-version and the adaptive h-version is plotted. For the adaptive h-version the estimated error is also shown, emphasizing the quality of the a posteriori error estimation. The reference value W = 12:16864134 for the strain energy of the two-dimensional Reissner-Mindlin model is obtained by a Richardson extrapolation based on the results of discretizations with 1963 quadrilaterals utilizing the tensor product space S p; p ps ( q st ) with p = 9; 10; 11.   From Figure 20 it becomes obvious, that the p-version is by far superior to the h-version. Applying a polynomial degree p = 10 with a resulting number of 52950 degrees of freedom a relative error in energy norm of 0.9% is obtained, while the ÿnest adapted MITC4 mesh supplies an error of 4.4 per cent with 51519 unknowns.
Next, we will consider the three-dimensional discretization of the plate with its columns. The mesh consisting of 194 hexahedral elements is constructed by sweeping the two-dimensional p-version mesh shown in Figure 21 into the third direction. The plate is discretized with one element in thickness direction, whereas each column is meshed with two elements. A series of computations using the trunk space S p ;pÁ;p ts ( h st ) is performed. The polynomial degree for the de ection u z is chosen to be one order higher than the degree for the displacements u x and u y .
The relative error in energy norm vs. the number of degrees of freedom is plotted in Figure 22. With a polynomial degree p = 7 and a corresponding number of 51393 unknowns the error in energy norm is reduced to 5.1 per cent, now with respect to the exact threedimensional solution being estimated from a Richardson extrapolation of the strain energies obtained with a mesh consisting of 202 hexahedral elements and an Ansatz space S p; p; p ts ( h st ) with polynomial degrees p = 7; 8; 9 resulting in W = 13:39080971. Note that the p-version for the three-dimensional problem yields an accuracy in the same range as that of an adaptive h-version (5.1 per cent vs 4.4 per cent) and a similar number of degrees of freedom. The signiÿcance of a three-dimensional computation is proven in Plate 1, where the deformed structure (scaling factor = 100) is plotted. The results are evaluated on a ÿne post-processing mesh being obtained from a subdivision of each p-element. It is obvious that the deformation of the structure exhibit a three-dimensional state which can of course not be captured by a two-dimensional model like the Reissner-Mindlin plate problem.
In Table II the computational cost and the error in energy norm of the two-and threedimensional discretizations are listed, all simulations being performed on a COMPAQ XP1000 machine (alpha processor ev6 21264 with 500 Mhz). To solve the overall equation system a PCG-solver either with SSOR (run 1,3) or incomplete Cholesky preconditioning (run 2) was applied. The computational time of the approximation based on the ÿnest mesh with MITC4elements amounts to 134 s with and to 128 s without the process of error estimation. The accuracy supplied by this discretization is about 4.4 per cent error in energy norm. With approximately the same number of unknowns and an accuracy of 0.9 per cent error in energy norm the p-version with (ÿ x ; ÿ y ; w) ∈ S 10; 10 ps ( q st ) results in a computational time of only 107s. From Table II it is evident, that the p-version supplies for given number of degrees of freedom not only a higher accuracy but also a lower computational time, when being compared to an adaptive h-version based on MITC4-element formulation.
A three-dimensional discretization with 194 elements of Ansatz (u x ; u y ) ∈ S 7; 7; 7 ts ( h st ); u z ∈ S 8; 8; 8 ts ( h st ) results in approximately the same number of unknowns as in the cases of the Reissner-Mindlin approximations. Although one observes that the computational time of 1260s is signiÿcantly higher than in the two-dimensional case, it is only a small fraction of the e ort which would be necessary using low-order three-dimensional elements. To e ciently integrate the (large) element sti ness matrices being the most expensive part in p-version computations we have further developed the vector-integration method ÿrst presented by Hinnant [26]. Our adaptive scheme [27; 28] results in an overall speed-up of 1.6 for this example, when being compared to the classical Gaussian integration. Further reduction of computational time could be achieved by applying more sophisticated preconditioners as being described e.g. in Reference [29].    Figure 23 shows the classical Scordelis-Lo shell having been used as a benchmark problem for shell structures by many authors [30]. The thickness of the shell is t = 0:25, the radius is R = 25 and the length equals to L = 50. The structure is loaded by a vertical shell weight = 360. At both ends the shell is supported such that u y = u z = 0. A linear elastic material law is assumed with E = 4:32 × 10 8 being Young's modulus and = 0 Poisson's ratio. The results of interest are the displacement u y at point B and the stress distribution along the cutline A-B at the middle surface of the shell. Due to symmetry only a quarter of the shell has to be discretized. Computations are performed on a mesh consisting of three hexahedral elements, being reÿned towards the free edge to resolve boundary layers (see Figure 23). In order to obtain a reference solution, a mesh with 84 hexahedral elements and a polynomial degree p = 8 in conjunction with the trunk space S p;p;p ts ( h st ) was used. This mesh contains two hexahedral elements in radial direction. The reference value of the strain energy based on this discretization is W = 1209:009850.  In Figure 24 the relative error of computations with isotropic Ansatz spaces S  ) with p = q = p = p Á = p = 1; : : : ; 8 is pictured. Again, it turns out that the trunk space is the most e cient one. Further improvement of e ciency can be obtained, if anisotropic Ansatz spaces are constructed, e.g. the polynomial degree in shell thickness direction is limited to p = q = 3 (see Figure 25). If we consider e.g. the solution obtained with the trunk space S p ;pÁ;p ts ( h st ) where p = p Á = 8; p = 3 we observe, that only 594 degrees of freedom are necessary to achieve an accuracy with approximately 1.4 per cent relative error in energy norm. Figure 26 shows the displacement u y of point B vs degrees of freedom for a discretization with

A complex shell model
Finally, we consider a more complex construction (see Figure 28), loaded by vertical weight = 100. The material behaviour is assumed to be linear elastic with E = 29 000 000 being Young's modulus and = 0:22 Poisson's ratio. It is composed of a spherical shell-like structure and a cylindrical solid. The radius of the spherical shell is R = 8 and the thickness equals t = 0:04, while the thickness of the cylindrical part is t=R = 0:18=6:9111. Displacements at the bottom of the construction are suppressed such that u x = u y = u z = 0. A mesh consisting of 45 hexahedral elements, taking advantage of symmetry is chosen to discretize the structure. The classical approach to this problem would demand for special elements in order to model the transition from shell-to solid elements (see right part of Figure 28). Due to the use of three-dimensional continuum p-elements the whole structure can be modelled with the same type of discretization and no transition elements are needed.
The relative error for a series of computations with an isotropic Ansatz space S p ;pÁ;p ts ( h st ) where p = p = p Á = p = 1; : : : ; 8 is plotted in Figure 29.
Using a polynomial degree of p = 8 and a corresponding number of 13 408 degrees of freedom an accuracy with approximately 3.1 per cent error in energy norm is achieved. This error  is estimated from an extrapolation of the strain energies obtained with Ansatz space S p;p;p ts ( h st ) and polynomial degrees p = 6; 7; 8 (see Reference [7]) resulting in W = 0:7051882141.
The deformed structure (scaling factor = 500) and the von Mises stress is plotted in Plate 2, where the results are evaluated on a ÿne post-processing mesh again, being obtained by a subdivision of each p-element.

CONCLUSIONS
In this paper an implementation of a three-dimensional p-version FEM for curved thin as well as thick walled structures is presented. It is demonstrated numerically, that the use of anisotropic Ansatz spaces in combination with the blending function method allows to e ciently compute the structural behaviour of plate and shell-like structures. Three di erent Ansatz spaces are investigated and it is shown that a exible implementation of the p-version enables to switch consistently from solid to shell-like structures, without using dimensionally reduced models. Future work will include the implementation of an anisotropic error estimation in order to automatically construct optimal anisotropic Ansatz spaces.